Changeset 26703

Timestamp:

12/03/21 07:51:10 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added results and starting test301 with basal friction

Location:

issm/trunk-jpl/src/jl

Files:

• md/classes.jl (modified) (6 diffs)
• md/utils.jl (modified) (1 diff)
• solve/analyses/stressbalanceanalysis.jl (modified) (2 diffs)
• solve/femmodel.jl (modified) (1 diff)
• solve/modules.jl (modified) (5 diffs)
• solve/results.jl (added)
• solve/solve.jl (modified) (3 diffs)
• test101.jl (modified) (1 diff)
• test301.jl (added)

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

@@ -15 +15 @@
 end# }}}
 function Base.show(io::IO, this::Mesh2dTriangle)# {{{
-
-        println(io,typeof(this),":")
-        for name in fieldnames(typeof(this))
-                a=getfield(this,name)
-                #print(io,"   $(name) = ")
-                @printf "%19s: " name
-                if isa(a,String)
-                        println(io, a)
-                elseif length(a)>1
-                        if !isempty(a)
-                                println(io, typeof(a), " of size ", size(a))
-                        else
-                                println(io,"empty")
-                        end
-                else
-                        println(io, a)
-                end
-        end
+        IssmStructDisp(io, this)
 end# }}}
 mutable struct Mesh3dPrism{T} <: AbstractMesh
@@ -66 +49 @@
         return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
 end# }}}
+function Base.show(io::IO, this::Mask)# {{{
+        IssmStructDisp(io, this)
+end# }}}
 
 mutable struct Initialization
@@ -87 +73 @@
 end# }}}
 function Base.show(io::IO, this::Stressbalance)# {{{
-
-        println(io,typeof(this),":")
-        for name in fieldnames(typeof(this))
-                a=getfield(this,name)
-                #print(io,"   $(name) = ")
-                @printf "%19s: " name
-                if isa(a,String)
-                        println(io, a)
-                elseif length(a)>1
-                        if !isempty(a)
-                                println(io, typeof(a), " of size ", size(a))
-                        else
-                                println(io,"empty")
-                        end
-                else
-                        println(io, a)
-                end
-        end
+        IssmStructDisp(io, this)
 end# }}}
 
@@ -137 +106 @@
 end# }}}
 
+abstract type AbstractFriction end
+mutable struct BuddFriction <: AbstractFriction
+        coefficient::Vector{Float64}
+end
+function BuddFriction() #{{{
+        return BuddFriction(Vector{Float64}(undef,0))
+end# }}}
+function Base.show(io::IO, this::BuddFriction)# {{{
+        IssmStructDisp(io, this)
+end# }}}
+mutable struct WeertmanFriction <: AbstractFriction
+        C::Vector{Float64}
+        m::Vector{Float64}
+end
+
 mutable struct model
         mesh::AbstractMesh
@@ -145 +129 @@
         stressbalance::Stressbalance
         constants::Constants
+        results
+        friction::AbstractFriction
 end
 function model() #{{{
-        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
+        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), [], BuddFriction())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -161 +147 @@
         @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# }}}

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

@@ -107 +107 @@
         return md
 end #}}}
+function IssmStructDisp(io::IO, modelfield::Any) # {{{
+        println(io,typeof(modelfield),":")
+        for name in fieldnames(typeof(modelfield))
+                a=getfield(modelfield,name)
+                #print(io,"   $(name) = ")
+                @printf "%19s: " name
+                if isa(a,String)
+                        println(io, a)
+                elseif length(a)>1
+                        if !isempty(a)
+                                println(io, typeof(a), " of size ", size(a))
+                        else
+                                println(io,"empty")
+                        end
+                else
+                        println(io, a)
+                end
+        end
+end #}}}

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

@@ -176 +176 @@
         #Now split solution vector into x and y components
         numnodes = 3
-        vx= Vector{Float64}(undef,numnodes)
-        vy= Vector{Float64}(undef,numnodes)
+        vx  = Vector{Float64}(undef,numnodes)
+        vy  = Vector{Float64}(undef,numnodes)
+        vel = Vector{Float64}(undef,numnodes)
         for i in 1:numnodes 
                 vx[i]=values[2*i-1] 
@@ -183 +184 @@
                 @assert isfinite(vx[i])
                 @assert isfinite(vy[i])
+
+                vel[i] =sqrt(vx[i]^2 + vy[i]^2)
         end
 
-        AddInput(element,VxEnum,vx,P1Enum)
-        AddInput(element,VyEnum,vy,P1Enum)
+        AddInput(element, VxEnum,  vx,  P1Enum)
+        AddInput(element, VyEnum,  vy,  P1Enum)
+        AddInput(element, VelEnum, vel, P1Enum)
 end#}}}

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

@@ -8 +8 @@
         constraints::Vector{Constraint}
         #loads::Vector{Loads}
+        results::Vector{Result}
 end#}}}
 

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

@@ -12 +12 @@
         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}())
@@ -36 +37 @@
 
         #Build FemModel
-        return FemModel(elements,vertices,nodes,parameters,inputs,constraints)
+        return FemModel(elements, vertices, nodes, parameters, inputs, constraints, results)
 end# }}}
 function CreateElements(elements::Vector{Tria},md::model) #{{{
@@ -80 +81 @@
         AddParam(parameters,md.materials.rho_freshwater,MaterialsRhoFreshwaterEnum)
         AddParam(parameters,md.constants.g,ConstantsGEnum)
+        
+        #Set step and time, this will be overwritten if we run a transient
+        AddParam(parameters,1,StepEnum)
+        AddParam(parameters,0.0,TimeEnum)
 end# }}}
 function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{
@@ -86 +91 @@
         FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
         FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum)
+end# }}}
+function OutputResultsx(femmodel::FemModel, md::model)# {{{
+
+        md.results = Dict("Stressbalance"=> femmodel.results)
+
 end# }}}
 
@@ -224 +234 @@
 function RequestedOutputsx(femmodel::FemModel,outputlist::Vector{IssmEnum})# {{{
 
-        error("not implemented")
-end# }}}
+        #Get Step and Time from parameters
+        step = FindParam(femmodel.parameters,StepEnum)
+        time = FindParam(femmodel.parameters,TimeEnum)
+        
+        #Now fetch results
+        for i in 1:length(outputlist)
+
+                #See if outputlist[i] is an input
+                if outputlist[i]>InputsSTARTEnum && outputlist[i]<InputsENDEnum
+
+                        #Create Result
+                        input  = GetInput(femmodel.inputs,outputlist[i])
+                        result = Result(step, time, copy(input.values))
+
+                        #Add to femmodel.results dataset
+                        push!(femmodel.results,result)
+
+                else
+                        println("Output ",outputlist[i]," not supported yet")
+                end
+        end
+end# }}}

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

@@ -9 +9 @@
 include("./elements.jl")
 include("./constraints.jl")
+include("./results.jl")
 include("./matice.jl")
 include("./femmodel.jl")
@@ -20 +21 @@
 function IssmCore(md::model) #{{{
 
+
         #Construct FemModel
         femmodel=ModelProcessor(md)
@@ -25 +27 @@
         #Solve
         analysis = StressbalanceAnalysis()
-        Core(analysis,femmodel)
+        Core(analysis, femmodel)
 
         #Save output
-        RequestedOutputsx(femmodel,[VxEnum,VyEnum])
-        error("STOP")
+        RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum])
+
+        #move results to md
+        OutputResultsx(femmodel, md)
 
 end# }}}

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

@@ -4 +4 @@
 
 md = ISSM.model()
-#md = ISSM.triangle(md,"../../test/Exp/Square.exp",180000.)
-md = ISSM.triangle(md,"../../test/Exp/Square.exp",250000.)
+md = ISSM.triangle(md,"../../test/Exp/Square.exp",50000.)
+#md = ISSM.triangle(md,"../../test/Exp/Square.exp",250000.)
 #md = ISSM.triangle(md,"../../test/Exp/Square.exp",1000.)
 md = ISSM.setmask(md,"all","")