Changeset 26731

Timestamp:

12/12/21 19:36:55 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with transient

Location:

issm/trunk-jpl/src/jl/core

Files:

• analyses/masstransportanalysis.jl (modified) (9 diffs)
• analyses/stressbalanceanalysis.jl (modified) (1 diff)
• analyses/transientanalysis.jl (modified) (2 diffs)
• elements.jl (modified) (1 diff)
• inputs.jl (modified) (2 diffs)
• modules.jl (modified) (3 diffs)
• solutionsequences.jl (modified) (1 diff)
• solve.jl (modified) (1 diff)
• toolkits.jl (modified) (2 diffs)

issm/trunk-jpl/src/jl/core/analyses/masstransportanalysis.jl ¶

@@ -38 +38 @@
         FetchDataToInput(md,inputs,elements,md.basalforcings.groundedice_melting_rate./md.constants.yts,BasalforcingsGroundediceMeltingRateEnum)
         FetchDataToInput(md,inputs,elements,md.basalforcings.floatingice_melting_rate./md.constants.yts,BasalforcingsFloatingiceMeltingRateEnum)
+        FetchDataToInput(md,inputs,elements,md.smb.mass_balance./md.constants.yts,SmbMassBalanceEnum)
         FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum)
         FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum)
@@ -54 +55 @@
 function Core(analysis::MasstransportAnalysis,femmodel::FemModel)# {{{
 
-        error("not supported yet")
+        println("   computing mass transport")
+        SetCurrentConfiguration!(femmodel, analysis)
+
+        InputDuplicatex(femmodel, ThicknessEnum, ThicknessOldEnum)
+        InputDuplicatex(femmodel, BaseEnum, BaseOldEnum)
+        InputDuplicatex(femmodel, SurfaceEnum, SurfaceOldEnum)
+
+        solutionsequence_linear(femmodel,analysis)
 
         #Save output
-        RequestedOutputsx(femmodel, [ThicknessEnum, SurfaceEnum, BedEnum])
+        RequestedOutputsx(femmodel, [ThicknessEnum, SurfaceEnum, BaseEnum])
 
 end #}}}
 function CreateKMatrix(analysis::MasstransportAnalysis,element::Tria)# {{{
-
-        error("not supported")
 
         #Internmediaries
@@ -70 +76 @@
         Ke = ElementMatrix(element.nodes)
         dbasis = Matrix{Float64}(undef,numnodes,2)
+        basis  = Vector{Float64}(undef,numnodes)
 
         #Retrieve all inputs and parameters
         xyz_list = GetVerticesCoordinates(element.vertices)
-        H_input  = GetInput(element, ThicknessEnum)
-
-        #Prepare material object
-        material = Matice(element)
-        
+        vx_input      = GetInput(element, VxEnum)
+        vy_input      = GetInput(element, VyEnum)
+        dt            = FindParam(element, TimesteppingTimeStepEnum)
+        stabilization = FindParam(element, MasstransportStabilizationEnum)
+
+        h = CharacteristicLength(element)
+
         #Start integrating
         gauss = GaussTria(2)
@@ -83 +92 @@
 
                 Jdet = JacobianDeterminant(xyz_list, gauss)
-                NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss)
-
-                H  = GetInputValue(H_input, gauss, ig)
-                mu = ViscositySSA(material, xyz_list, gauss, ig)
-
+                NodalFunctionsDerivatives(element, dbasis, xyz_list, gauss)
+                NodalFunctions(element, basis, gauss, ig, P1Enum)
+
+      #Transient term
                 for i in 1:numnodes
                         for j in 1:numnodes
-                                Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,1]*dbasis[i,1] + dbasis[j,2]*dbasis[i,2])
-                                Ke.values[2*i-1,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,2]*dbasis[i,1] + dbasis[j,1]*dbasis[i,2])
-                                Ke.values[2*i  ,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,1]*dbasis[i,2] + dbasis[j,2]*dbasis[i,1])
-                                Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,2]*dbasis[i,2] + dbasis[j,1]*dbasis[i,1])
+                                Ke.values[i ,j] += gauss.weights[ig]*Jdet*basis[i]*basis[j]
                         end
                 end
-        end
-
-        #Add basal friction
-        phi=GetGroundedPortion(element, xyz_list)
-
-        if(phi>0)
-                basis = Vector{Float64}(undef,numnodes)
-                friction = CoreFriction(element)
-
-                #Start integrating
-                gauss = GaussTria(2)
-                for ig in 1:gauss.numgauss
-
-                        Jdet = JacobianDeterminant(xyz_list, gauss)
-                        NodalFunctions(element, basis, gauss, ig, P1Enum)
-
-                        alpha2 = Alpha2(friction, gauss, ig)
-
-                        for i in 1:numnodes
-                                for j in 1:numnodes
-                                        Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j]
-                                        Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j]
-                                end
-                        end
+
+      #Advection term
+                vx  = GetInputValue(vx_input, gauss, ig)
+      vy  = GetInputValue(vy_input, gauss, ig)
+                dvx = GetInputDerivativeValue(vx_input, xyz_list, gauss, ig)
+      dvy = GetInputDerivativeValue(vy_input, xyz_list, gauss, ig)
+      for i in 1:numnodes
+         for j in 1:numnodes
+            #\phi_i \phi_j \nabla\cdot v
+            Ke.values[i ,j] += dt*gauss.weights[ig]*Jdet*basis[i]*basis[j]*(dvx[1] + dvy[2])
+            #\phi_i v\cdot\nabla\phi_j
+            Ke.values[i ,j] += dt*gauss.weights[ig]*Jdet*basis[i]*(vx*dbasis[j,1] + vy*dbasis[j,2])
+         end
+      end
+
+                #Stabilization
+                if(stabilization==0)
+                        #do nothing
+                elseif (stabilization==1)
+                        vx = GetInputAverageValue(vx_input)
+                        vy = GetInputAverageValue(vy_input)
+                        D  = dt*gauss.weights[ig]*Jdet*[h/2*abs(vx) 0; 0 h/2*abs(vy)]
+                        for i in 1:numnodes; for j in 1:numnodes
+                                        Ke.values[i ,j] += (dbasis[i,1]*(D[1,1]*dbasis[j,1] + D[1,2]*dbasis[j,2]) +
+                                                                                          dbasis[i,2]*(D[2,1]*dbasis[j,1] + D[2,2]*dbasis[j,2]))
+                        end end
+                else
+                        error("Stabilization ",stabilization, " not supported yet")
                 end
         end
@@ -126 +135 @@
 end #}}}
 function CreatePVector(analysis::MasstransportAnalysis,element::Tria)# {{{
-
-        error("not supported")
 
         #Internmediaries
@@ -138 +145 @@
         #Retrieve all inputs and parameters
         xyz_list = GetVerticesCoordinates(element.vertices)
-        H_input  = GetInput(element, ThicknessEnum)
-        s_input  = GetInput(element, SurfaceEnum)
-        rho_ice  = FindParam(element, MaterialsRhoIceEnum)
-        g        = FindParam(element, ConstantsGEnum)
+        H_input         = GetInput(element, ThicknessEnum)
+        gmb_input       = GetInput(element, BasalforcingsGroundediceMeltingRateEnum)
+        fmb_input       = GetInput(element, BasalforcingsFloatingiceMeltingRateEnum)
+   smb_input       = GetInput(element, SmbMassBalanceEnum)
+        olevelset_input = GetInput(element, MaskOceanLevelsetEnum)
+        dt            = FindParam(element, TimesteppingTimeStepEnum)
+        stabilization = FindParam(element, MasstransportStabilizationEnum)
+
+   #How much is actually grounded?
+        phi=GetGroundedPortion(element, xyz_list)
 
         #Start integrating
-        gauss = GaussTria(2)
+        gauss = GaussTria(3)
         for ig in 1:gauss.numgauss
 
@@ -150 +163 @@
                 NodalFunctions(element, basis, gauss, ig, P1Enum)
 
-                H  = GetInputValue(H_input, gauss, ig)
-                ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig)
+                smb = GetInputValue(smb_input, gauss, ig)
+      H   = GetInputValue(H_input, gauss, ig)
+
+      #Only apply melt on fully floating cells
+      if(phi<0.00000001)
+         mb = GetInputValue(fmb_input, gauss, ig)
+      else
+         mb = GetInputValue(gmb_input, gauss, ig)
+      end
 
                 for i in 1:numnodes
-                        pe.values[2*i-1] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[1]*basis[i]
-                        pe.values[2*i  ] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[2]*basis[i]
-                end
-        end
-
-        if(IsIcefront(element))
-
-                #Get additional parameters and inputs
-                b_input   = GetInput(element, BaseEnum)
-                rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
-
-                #Get normal and ice front coordinates
-                xyz_list_front = Matrix{Float64}(undef,2,3)
-                GetIcefrontCoordinates!(element, xyz_list_front, xyz_list, MaskIceLevelsetEnum)
-                nx, ny = NormalSection(element, xyz_list_front)
-
-                gauss = GaussTria(element, xyz_list, xyz_list_front, 3)
-                for ig in 1:gauss.numgauss
-
-                        Jdet = JacobianDeterminantSurface(xyz_list_front, gauss)
-                        NodalFunctions(element, basis, gauss, ig, P1Enum)
-
-                        H  = GetInputValue(H_input, gauss, ig)
-                        b  = GetInputValue(b_input, gauss, ig)
-                        sl = 0
-
-                        term = 0.5*g*rho_ice*H^2 + 0.5*g*rho_water*(min(0, H+b-sl)^2 - min(0, b-sl)^2)
-
-                        for i in 1:numnodes
-                                pe.values[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i]
-                                pe.values[2*i  ] += gauss.weights[ig]*Jdet*term*ny*basis[i]
-                        end
+         pe.values[i] += gauss.weights[ig]*Jdet*(H + dt*(smb - mb))*basis[i]
                 end
         end
@@ -237 +226 @@
 
                 #Enforce minimum thickness
-                if(thickness[i]<min_thickness)
-                        thickness[i] = min_thickness
+                if(thickness[i]<H_min)
+                        thickness[i] = H_min
                 end
         end
@@ -255 +244 @@
    GetInputListOnVertices!(element, phi, MaskOceanLevelsetEnum)
    sealevel = zeros(3)
+        newsurface = Vector{Float64}(undef,3)
+        newbase    = Vector{Float64}(undef,3)
 
    for i in 1:3

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

@@ -82 +82 @@
 
         #Call solution sequence to compute new speeds
+        println("   computing stress balance");
         solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol)
 

issm/trunk-jpl/src/jl/core/analyses/transientanalysis.jl ¶

@@ -6 +6 @@
 
         AddParam(parameters, md.constants.yts, ConstantsYtsEnum)
-        AddParam(parameters, md.timestepping.start_time/md.constants.yts, TimeEnum)
-        AddParam(parameters, md.timestepping.final_time/md.constants.yts, TimesteppingFinalTimeEnum)
-        AddParam(parameters, md.timestepping.time_step/md.constants.yts,  TimesteppingTimeStepEnum)
+        AddParam(parameters, md.timestepping.start_time*md.constants.yts, TimeEnum)
+        AddParam(parameters, md.timestepping.final_time*md.constants.yts, TimesteppingFinalTimeEnum)
+        AddParam(parameters, md.timestepping.time_step*md.constants.yts,  TimesteppingTimeStepEnum)
         AddParam(parameters, md.transient.isstressbalance, TransientIsstressbalanceEnum)
         AddParam(parameters, md.transient.ismasstransport, TransientIsmasstransportEnum)
@@ -26 +26 @@
    while(time < finaltime - (yts*eps(Float64))) #make sure we run up to finaltime.
 
-                step+=1
                 time+=dt
-                AddParam(femmodel.parameters, time, TimeEnum);
-                AddParam(femmodel.parameters, step, StepEnum);
-      println("iteration ", step, "/", ceil((finaltime-time)/dt)+step," time [yr]: ", time/yts, " (time step: ",  dt/yts, ")")
+                AddParam(femmodel.parameters, time, TimeEnum)
+                AddParam(femmodel.parameters, step, StepEnum)
+                println("iteration ", step, "/", Int(ceil((finaltime-time)/dt))+step," time [yr]: ", time/yts, " (time step: ",  dt/yts, " [yr])")
 
       if(isstressbalance) Core(StressbalanceAnalysis(), femmodel) end
       if(ismasstransport) Core(MasstransportAnalysis(), femmodel) end
+
+                step+=1
    end
 

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

@@ -253 +253 @@
         return nx, ny
 end#}}}
+function CharacteristicLength(element::Tria) #{{{
+
+        return sqrt(2*GetArea(element))
+end#}}}
 
 #Finite Element stuff

issm/trunk-jpl/src/jl/core/inputs.jl ¶

@@ -97 +97 @@
         input.values[indices] = values
 end#}}}
-function GetInputAverageValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
+function GetInputAverageValue(input::ElementInput) #{{{
 
         numnodes = NumberofNodesTria(input.interp)
@@ -151 +151 @@
 
 end#}}}
+function DuplicateInput(inputs::Inputs, old::IssmEnum, new::IssmEnum)#{{{
+
+        #Fetch input that needs to be copied
+        oldinput = inputs.lookup[old]
+
+        if typeof(oldinput)==ElementInput
+                inputs.lookup[new] = ElementInput(new, oldinput.interp, copy(oldinput.values), copy(oldinput.element_values))
+        end
+
+end#}}}

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

@@ -115 +115 @@
         FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum)
 end# }}}
-function OutputResultsx(femmodel::FemModel, md::model)# {{{
-
-        output = Dict()
-        for i in length(femmodel.results)
-                result = femmodel.results[i]
-                output[EnumToString(result.enum)] = result.value
-        end
-        md.results["Stressbalance"] = output
+function OutputResultsx(femmodel::FemModel, md::model, solution::String)# {{{
+
+
+        if solution=="TransientSolution"
+
+                #Compute maximum number of steps
+                maxstep = 0
+                for i in length(femmodel.results)
+                        if(femmodel.results[i].step>maxstep) maxstep = femmodel.results[i].step end
+                end
+
+                #Initialize vector now that we know the size
+                output = Vector{Dict}(undef, maxstep)
+                for i in 1:maxstep; output[i] = Dict() end
+
+                #Insert results in vector
+                for i in 1:length(femmodel.results)
+                        result = femmodel.results[i]
+                        step   = femmodel.results[i].step
+                        (output[step])[EnumToString(result.enum)] = result.value
+                end
+        else
+                output = Dict()
+                for i in length(femmodel.results)
+                        result = femmodel.results[i]
+                        output[EnumToString(result.enum)] = result.value
+                end
+        end
+
+        md.results[solution] = output
 
 end# }}}
@@ -182 +204 @@
         return ug
 
+end#}}}
+function InputDuplicatex(femmodel::FemModel, oldenum::IssmEnum, newenum::IssmEnum) #{{{
+        DuplicateInput(femmodel.inputs, oldenum, newenum)
 end#}}}
 function Reducevectorgtofx(ug::IssmVector,nodes::Vector{Node}) #{{{
@@ -281 +306 @@
 
                         #Add to femmodel.results dataset
-                        push!(femmodel.results,result)
+                        push!(femmodel.results, result)
 
                 else

issm/trunk-jpl/src/jl/core/solutionsequences.jl ¶

@@ +1 @@
+function solutionsequence_linear(femmodel::FemModel,analysis::Analysis) # {{{
+
+        #Get new matrices
+        Kff, Kfs, pf = SystemMatricesx(femmodel,analysis)
+
+        #Enforce constraints
+        ys = CreateNodalConstraintsx(femmodel.nodes)
+        Reduceloadx!(pf, Kfs, ys)
+
+        #Solve!
+        uf = Solverx(Kff, pf)
+
+        #Merge uf with ys and update inputs
+        gsize = NumberOfDofs(femmodel.nodes,GsetEnum)
+        ug = IssmVector(gsize)
+        Mergesolutionfromftogx(ug, uf, ys, femmodel.nodes)
+
+        InputUpdateFromSolutionx(analysis, ug, femmodel)
+
+end# }}}
 function solutionsequence_nonlinear(femmodel::FemModel,analysis::Analysis,maxiter::Int64,restol::Float64,reltol::Float64,abstol::Float64) # {{{
 

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

@@ -41 +41 @@
 
         #move results to md
-        OutputResultsx(femmodel, md)
+        OutputResultsx(femmodel, md, solutionstring)
 
 end# }}}

issm/trunk-jpl/src/jl/core/toolkits.jl ¶

@@ -125 +125 @@
 
 end#}}}
-function Solverx(A::IssmMatrix,b::IssmVector,xold::IssmVector) #{{{
+function Solverx(A::IssmMatrix, b::IssmVector, xold::IssmVector) #{{{
 
         #Initialize output
         #x = IssmVector(GetSize(xold))
+        
+        return Solverx(A, b)
+
+end#}}}
+function Solverx(A::IssmMatrix, b::IssmVector) #{{{
+
+        #Initialize output
         x = IssmVector(0)
 
@@ -136 +143 @@
         return x
 
-
 end#}}}