Changeset 26730

Timestamp:

12/12/21 14:16:13 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: adding mass transport, not quite working yet

Location:

issm/trunk-jpl/src/jl

Files:

• core/analyses/masstransportanalysis.jl (added)
• core/analyses/stressbalanceanalysis.jl (modified) (3 diffs)
• core/analyses/transientanalysis.jl (added)
• core/constraints.jl (modified) (1 diff)
• core/elements.jl (modified) (3 diffs)
• core/femmodel.jl (modified) (1 diff)
• core/modules.jl (modified) (5 diffs)
• core/parameters.jl (modified) (3 diffs)
• core/solve.jl (modified) (2 diffs)
• tests/test.jl (modified) (1 diff)
• tests/test101.jl (modified) (1 diff)
• tests/test208.jl (added)
• usr/classes.jl (modified) (4 diffs)
• usr/utils.jl (modified) (1 diff)

issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl

@@ -33 +33 @@
 
         #Provide node indices to element
-        for i in 1:md.mesh.numberofvertices
+        for i in 1:md.mesh.numberofelements
                 Update(elements[i],inputs,i,md,P1Enum)
         end
@@ -72 +72 @@
 function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{
 
+        #Set current analysis to Stressnalance
+        SetCurrentConfiguration!(femmodel, analysis)
+
         #Fetch parameters relevant to solution sequence
         maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum)
@@ -78 +81 @@
         abstol  = FindParam(femmodel.parameters,StressbalanceAbstolEnum)
 
+        #Call solution sequence to compute new speeds
         solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol)
+
+        #Save output
+        RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum])
 
 end #}}}

issm/trunk-jpl/src/jl/core/constraints.jl

@@ -14 +14 @@
 
         #Apply Constraint
-        ApplyConstraint(node,constraint.dof,constraint.value)
+        ApplyConstraint(node, constraint.dof, constraint.value)
 
 end# }}}

issm/trunk-jpl/src/jl/core/elements.jl

@@ -3 +3 @@
         sid::Int64
         pid::Int64
+
         #vertexids::Int64[3]
         #vertices::Vertex[3]
         vertexids::Vector{Int64}
         vertices::Vector{Vertex}
-        nodeids::Vector{Int64}
+
         nodes::Vector{Node}
+        nodes_list::Vector{Vector{Node}}
+        nodes_ids_list::Vector{Vector{Int64}}
+
         parameters::Parameters
         inputs::Inputs
 end# }}}
+function Tria(sid::Int64, pid::Int64, vertexids::Vector{Int64}) #{{{
+
+        #This is the default constructor, at this point we don't have much information
+        tempparams   = Parameters(Dict{IssmEnum,Parameter}())
+        tempinputs   = Inputs(-1,-1,Dict{IssmEnum,Input}())
+        return Tria(sid, pid,
+                                        vertexids, Vector{Vertex}(undef,3),
+                                        Vector{Node}(undef,0), Vector{Vector{Node}}(undef,0), Vector{Vector{Int64}}(undef,0),
+                                        tempparams, tempinputs)
+
+end #}}}
 
 #Element functions
@@ -31 +46 @@
         end
 end #}}}
-function Update(element::Tria,inputs::Inputs,index::Int64,md::model,finiteelement::IssmEnum) #{{{
+function Update(element::Tria, inputs::Inputs, index::Int64, md::model, finiteelement::IssmEnum) #{{{
 
         if finiteelement==P1Enum
                 numnodes = 3
-                element.nodeids    = Vector{Int64}(undef,numnodes)
-                element.nodeids[1] = md.mesh.elements[index,1]
-                element.nodeids[2] = md.mesh.elements[index,2]
-                element.nodeids[3] = md.mesh.elements[index,3]
+                nodeids    = Vector{Int64}(undef,numnodes)
+                nodeids[1] = md.mesh.elements[index,1]
+                nodeids[2] = md.mesh.elements[index,2]
+                nodeids[3] = md.mesh.elements[index,3]
+
+                push!(element.nodes_ids_list, nodeids)
+                push!(element.nodes_list, Vector{Node}(undef, numnodes))
         else
                 error("not supported yet")
         end
 end #}}}
-function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) # {{{
+function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs,index::Int64) # {{{
 
    #Configure vertices
@@ -51 +69 @@
 
         #Configure nodes (assuming P1 finite elements)
-        for i in 1:3
-                element.nodes[i] = nodes[element.nodeids[i]]
+        nodes_list     = element.nodes_list[index]
+        nodes_ids_list = element.nodes_ids_list[index]
+        for i in 1:3
+                nodes_list[i] = nodes[nodes_ids_list[i]]
         end
 
         #Point to real datasets
+        element.nodes      = element.nodes_list[index]
         element.parameters = parameters
         element.inputs     = inputs

issm/trunk-jpl/src/jl/core/femmodel.jl

@@ -1 +1 @@
 #femmodel class definition
 mutable struct FemModel #{{{
+        analyses::Vector{Analysis}
         elements::Vector{Tria}
         vertices::Vector{Vertex}
+
         nodes::Vector{Node}
+        nodes_list::Vector{Vector{Node}}
+
         parameters::Parameters
         inputs::Inputs
+
         constraints::Vector{Constraint}
+        constraints_list::Vector{Vector{Constraint}}
+
         #loads::Vector{Loads}
+
         results::Vector{Result}
 end#}}}
 
 #femmodel functions
+function SetCurrentConfiguration!(femmodel::FemModel, analysis::Analysis) #{{{
 
+        #Find the index of this analysis
+        index = -1
+        for i in 1:length(femmodel.analyses)
+                if(typeof(femmodel.analyses[i]) == typeof(analysis)) index = i end
+        end
+        if(index<1) error("Could not find analysis ",analysis, " in femmodel") end
+
+        #Plug right nodes onto element
+        for i in 1:length(femmodel.elements)
+                femmodel.elements[i].nodes = femmodel.elements[i].nodes_list[index]
+        end
+
+        #Plug in nodes and other datasets
+        femmodel.nodes       = femmodel.nodes_list[index]
+        femmodel.constraints = femmodel.constraints_list[index]
+
+end#}}}

issm/trunk-jpl/src/jl/core/modules.jl

@@ -5 +5 @@
         end
 end#}}}
-function ModelProcessor(md::model) #{{{
+function ModelProcessor(md::model, solutionstring::String) #{{{
 
         #Initialize structures
         elements    = Vector{Tria}(undef,0)
         vertices    = Vector{Vertex}(undef,0)
-        nodes       = Vector{Node}(undef,0)
-        constraints = Vector{Constraint}(undef,0)
         results     = Vector{Result}(undef,0)
         parameters  = Parameters(Dict{IssmEnum,Parameter}())
-        inputs      = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
+        inputs      = Inputs(md.mesh.numberofelements, md.mesh.numberofvertices, Dict{IssmEnum,Input}())
 
         #Create  elements, vertices and materials (independent of the analysis)
-        CreateElements(elements,md)
-        CreateVertices(vertices,md)
-        CreateParameters(parameters,md)
-        CreateInputs(inputs,elements,md)
+        CreateElements(elements, md)
+        CreateVertices(vertices, md)
+        CreateParameters(parameters, md)
+        CreateInputs(inputs,elements, md)
+        if solutionstring=="TransientSolution"
+                UpdateParameters(TransientAnalysis(), parameters, md)
+        end
 
         #Now create analysis specific data structure
-        analysis = StressbalanceAnalysis()
-        UpdateParameters(analysis,parameters,md)
-        CreateNodes(analysis,nodes,md)
-        UpdateElements(analysis,elements,inputs,md)
-        CreateConstraints(analysis,constraints,md)
-
-        #Configure objects
-        ConfigureObjectx(elements,nodes,vertices,parameters,inputs)
-
-        #Constrain and Number nodes
-        SpcNodesx(nodes,constraints,parameters)
-        NodesDofx(nodes,parameters)
+        if solutionstring=="StressbalanceSolution"
+                analyses = [StressbalanceAnalysis()]
+        elseif solutionstring=="TransientSolution"
+                analyses = [StressbalanceAnalysis(), MasstransportAnalysis()]
+        else
+                error(solutionstring, " not supported by ModelProcessor")
+        end
+
+        #Initialize analysis specific datasets
+        numanalyses = length(analyses)
+        nodes       = Vector{Vector{Node}}(undef,numanalyses)
+        constraints = Vector{Vector{Constraint}}(undef,numanalyses)
+        for i in 1:numanalyses
+                analysis = analyses[i]
+                println("   creating datasets for analysis ", typeof(analysis))
+                nodes[i]       = Vector{Node}(undef,0)
+                constraints[i] = Vector{Constraint}(undef,0)
+
+                UpdateParameters(analysis, parameters, md)
+                UpdateElements(analysis, elements, inputs, md)
+                CreateNodes(analysis, nodes[i], md)
+                CreateConstraints(analysis, constraints[i], md)
+
+                #Configure objects
+                ConfigureObjectx(elements, nodes[i], vertices, parameters, inputs, i)
+
+                #Constrain and Number nodes
+                SpcNodesx(nodes[i], constraints[i], parameters)
+                NodesDofx(nodes[i], parameters)
+        end
 
         #Build FemModel
-        return FemModel(elements, vertices, nodes, parameters, inputs, constraints, results)
+        return FemModel(analyses, elements, vertices,
+                                                 Vector{Node}(undef,0), nodes,
+                                                 parameters, inputs,
+                                                 Vector{Constraint}(undef,0), constraints,
+                                                 results)
 end# }}}
 function CreateElements(elements::Vector{Tria},md::model) #{{{
@@ -44 +67 @@
         @assert length(elements)==0
 
-        tempparams   = Parameters(Dict{IssmEnum,Parameter}())
-        tempinputs   = Inputs(-1,-1,Dict{IssmEnum,Input}())
+
 
         count = 0
@@ -55 +77 @@
 
                 #Call constructor and add to dataset elements
-                push!(elements,Tria(i,count,vertexids,Vector{Vertex}(undef,3),nodeids,Vector{Node}(undef,3),tempparams,tempinputs))
+                push!(elements,Tria(i,count, vertexids))
         end
 
@@ -85 +107 @@
         AddParam(parameters,1,StepEnum)
         AddParam(parameters,0.0,TimeEnum)
+
 end# }}}
 function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{
@@ -104 +127 @@
 
 #Other modules
-function ConfigureObjectx(elements::Vector{Tria},nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) #{{{
+function ConfigureObjectx(elements::Vector{Tria}, nodes::Vector{Node}, vertices::Vector{Vertex}, parameters::Parameters, inputs::Inputs, analysis::Int64) #{{{
 
         for i in 1:length(elements)
-                Configure(elements[i],nodes,vertices,parameters,inputs)
+                Configure(elements[i], nodes, vertices, parameters, inputs, analysis)
         end
 

issm/trunk-jpl/src/jl/core/parameters.jl

@@ -8 +8 @@
         enum::IssmEnum
         value::Int64
+end# }}}
+struct BoolParam <: Parameter #{{{
+        enum::IssmEnum
+        value::Bool
 end# }}}
 
@@ -22 +26 @@
         return param.value
 end#}}}
+function GetParameterValue(param::BoolParam) #{{{
+        return param.value
+end#}}}
 
 #Parameters functions
@@ -29 +36 @@
 
 end#}}}
-function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
+function AddParam(parameters::Parameters,value::Int64, enum::IssmEnum) #{{{
 
         parameters.lookup[enum] = IntParam(enum,value)
 
 end#}}}
-function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
+function AddParam(parameters::Parameters,value::Bool, enum::IssmEnum) #{{{
 
-        parameters.lookup[enum] = IntParam(enum,value)
+        parameters.lookup[enum] = BoolParam(enum,value)
 
 end#}}}

issm/trunk-jpl/src/jl/core/solve.jl

@@ -12 +12 @@
 include("./matice.jl")
 include("./friction.jl")
+include("./analyses/analysis.jl")
 include("./femmodel.jl")
-include("./analyses/analysis.jl")
+
+#All analyses
 include("./analyses/stressbalanceanalysis.jl")
+include("./analyses/masstransportanalysis.jl")
+include("./analyses/transientanalysis.jl")
+
 include("./solutionsequences.jl")
 include("./modules.jl")
@@ -20 +25 @@
 include("./utils.jl")
 
-function IssmCore(md::model) #{{{
-
+function IssmCore(md::model, solutionstring::String) #{{{
 
         #Construct FemModel
-        femmodel=ModelProcessor(md)
+        femmodel=ModelProcessor(md, solutionstring)
 
-        #Solve
-        analysis = StressbalanceAnalysis()
+        #Solve FIXME: to be improved later...
+        if(solutionstring=="StressbalanceSolution")
+                analysis = StressbalanceAnalysis()
+        elseif (solutionstring=="TransientSolution")
+                analysis = TransientAnalysis()
+        else
+                error("not supported")
+        end
         Core(analysis, femmodel)
-
-        #Save output
-        RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum])
 
         #move results to md

issm/trunk-jpl/src/jl/tests/test.jl

@@ -5 +5 @@
 
 #Load model from MATLAB file
-file = matopen("./Pig.mat")
+#file = matopen("./Pig.mat")
+file = matopen("/Users/mmorligh/Desktop/issmuci/trunk-jpl/test/NightlyRun/temp.mat")
 mat  = read(file, "md")
 close(file)

issm/trunk-jpl/src/jl/tests/test101.jl

@@ -25 +25 @@
 vy    = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy")
 index = Int.(archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index"))
-md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
-md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vx=0 .*InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0)
+md.initialization.vy=0 .*InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0)
 
 md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices)

issm/trunk-jpl/src/jl/usr/classes.jl

@@ -142 +142 @@
 end
 # }}}
+#Basalforcings {{{
+mutable struct Basalforcings
+        groundedice_melting_rate::Vector{Float64}
+        floatingice_melting_rate::Vector{Float64}
+end
+function Basalforcings() #{{{
+        return Basalforcings( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
+function Base.show(io::IO, this::Basalforcings)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+# }}}
+#Surfaceforcings {{{
+mutable struct SMBforcings
+        mass_balance::Vector{Float64}
+end
+function SMBforcings() #{{{
+        return SMBforcings( Vector{Float64}(undef,0))
+end# }}}
+function Base.show(io::IO, this::SMBforcings)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+# }}}
+#Timestepping{{{
+abstract type AbstractTimestepping end
+mutable struct Timestepping <: AbstractTimestepping
+        start_time::Float64
+        final_time::Float64
+        time_step::Float64
+end
+function Timestepping() #{{{
+        return Timestepping( 0., 0., 0.)
+end# }}}
+function Base.show(io::IO, this::Timestepping)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+# }}}
+#Masstransport {{{
+mutable struct Masstransport
+        spcthickness::Vector{Float64}
+        min_thickness::Float64
+        stabilization::Int64
+end
+function Masstransport() #{{{
+        return Masstransport( Vector{Float64}(undef,0), 10.0, 1)
+end# }}}
+function Base.show(io::IO, this::Masstransport)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+# }}}
+#Transient {{{
+mutable struct Transient
+        issmb::Bool
+        ismasstransport::Bool
+        isstressbalance::Bool
+        isgroundingline::Bool
+        ismovingfront::Bool
+end
+function Transient() #{{{
+        return Transient( true, true, true, true, true)
+end# }}}
+function Base.show(io::IO, this::Transient)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+# }}}
 
 #Model structure
@@ -154 +219 @@
         results::Dict
         friction::AbstractFriction
+        basalforcings::Basalforcings
+        smb::SMBforcings
+        timestepping::Timestepping
+        masstransport::Masstransport
+        transient::Transient
 end
 function model() #{{{
-        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), Dict(), BuddFriction())
+        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(),
+                                         Initialization(),Stressbalance(), Constants(), Dict(),
+                                         BuddFriction(), Basalforcings(), SMBforcings(), Timestepping(),
+                                         Masstransport(), Transient())
 end#}}}
 function model(matmd::Dict) #{{{
@@ -185 +258 @@
                                 setfield!(mdfield, Symbol(name2), Int64(value_matlab))
 
+                        elseif typeof(value_matlab)==Float64 && typeof(value_julia)==Bool
+                                setfield!(mdfield, Symbol(name2), Bool(value_matlab))
+
                         elseif typeof(value_matlab)==Matrix{Float64} && typeof(value_julia)==Vector{Float64}
                                 if(size(value_matlab,2)!=1) error("only one column expected") end
@@ -218 +294 @@
         @printf "%19s: %-26s -- %s\n" "materials" typeof(md.materials) "material properties"
         @printf "%19s: %-26s -- %s\n" "initialization" typeof(md.initialization) "initial state"
-        @printf "%19s: %-26s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
         @printf "%19s: %-26s -- %s\n" "constants" typeof(md.constants) "physical constants"
         @printf "%19s: %-26s -- %s\n" "friction" typeof(md.friction) "basal friction"
-
-end# }}}
-
-
+        @printf "%19s: %-26s -- %s\n" "basalforcings" typeof(md.basalforcings) "basal forcings"
+        @printf "%19s: %-26s -- %s\n" "smb" typeof(md.smb) "surface mass balance"
+        @printf "%19s: %-26s -- %s\n" "timestepping" typeof(md.timestepping) "time stepping for transient simulations"
+        @printf "%19s: %-26s -- %s\n" "stressbalance" typeof(md.stressbalance) "parameters stress balance simulations"
+        @printf "%19s: %-26s -- %s\n" "masstransport" typeof(md.masstransport) "parameters mass transport simulations"
+        @printf "%19s: %-26s -- %s\n" "transient" typeof(md.transient) "parameters for transient simulations"
+
+end# }}}

issm/trunk-jpl/src/jl/usr/utils.jl

@@ -207 +207 @@
         if solution=="sb" || solution=="Stressbalance"
                 solutionstring = "StressbalanceSolution"
+        elseif solution=="tr" || solution=="Transient"
+                        solutionstring = "TransientSolution"
         else
                 error("solutionstring "*solution*" not supported!");
         end
 
-        IssmCore(md)
+        IssmCore(md, solutionstring)
 
         return md