Changeset 26705

Timestamp:

12/03/21 11:18:50 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added subelement GL calculation

Location:

issm/trunk-jpl/src/jl/solve

Files:

• analyses/stressbalanceanalysis.jl (modified) (3 diffs)
• elements.jl (modified) (2 diffs)
• friction.jl (modified) (2 diffs)

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

@@ -44 +44 @@
         FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum)
         FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum)
-
-        FetchDataToInput(md,inputs,elements,md.friction.coefficient,FrictionCoefficientEnum)
+        FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum)
+        FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum)
+
+        #Deal with friction
+        if typeof(md.friction)==BuddFriction
+                FetchDataToInput(md,inputs,elements,md.friction.coefficient,FrictionCoefficientEnum)
+        else
+                error("Friction ", typeof(md.friction), " not supported yet")
+        end
 
 end#}}}
@@ -53 +60 @@
         AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum)
         AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum)
+
+        #Deal with friction
+        if typeof(md.friction)==BuddFriction
+                AddParam(parameters, 1, FrictionLawEnum)
+        else
+                error("Friction ", typeof(md.friction), " not supported yet")
+        end
 end#}}}
 
@@ -104 +118 @@
 
         #Add basal friction
-        basis = Vector{Float64}(undef,numnodes)
-        friction = Friction(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*alpha2*basis[i]*basis[j]
-                                Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*alpha2*basis[i]*basis[j]
+        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
                 end
         end
-        
 
         return Ke

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

@@ -89 +89 @@
 
 end # }}}
+function GetInputListOnVertices!(element::Tria, vector::Vector{Float64}, enum::IssmEnum) # {{{
+
+        #Get Input first 
+        input = GetInput(element, enum)
+
+        #Get value at each vertex (i.e. P1 Nodes)
+        gauss=GaussTria(P1Enum)
+        for i in 1:gauss.numgauss
+                vector[i] = GetInputValue(input, gauss, i)
+        end
+
+end # }}}
 function FindParam(element::Tria,enum::IssmEnum) # {{{
 
@@ -107 +119 @@
 
 end # }}}
+function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{
+
+        level = Vector{Float64}(undef,3)
+        GetInputListOnVertices!(element, level, MaskOceanLevelsetEnum)
+
+        #Be sure that values are not zero
+        epsilon = 1.e-15
+        for i in 1:3
+                if(level[i]==0.) level[i]=level[i]+epsilon end
+        end
+
+        if level[1]>0 && level[2]>0 && level[3]>0
+      #Completely grounded
+                phi = 1.0
+   elseif level[1]<0 && level[2]<0 && level[3]<0
+      #Completely floating
+      phi = 0.0
+   else
+                #Partially floating,
+                if(level[1]*level[2]>0) #Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+                        s1=level[3]/(level[3]-level[2]);
+                        s2=level[3]/(level[3]-level[1]);
+                elseif(level[2]*level[3]>0) #Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+                        s1=level[1]/(level[1]-level[2]);
+                        s2=level[1]/(level[1]-level[3]);
+                elseif(level[1]*level[3]>0) #Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+                        s1=level[2]/(level[2]-level[1]);
+                        s2=level[2]/(level[2]-level[3]);
+                else
+                        error("not supposed to be here...")
+                end
+
+                if(level[1]*level[2]*level[3]>0)
+                        phi = s1*s2
+                else
+                        phi = (1-s1*s2)
+                end
+        end
+
+        return phi
+end#}}}
 
 #Finite Element stuff

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

@@ -1 +1 @@
 #Friction class definition
-struct Friction#{{{
+
+abstract type CoreFriction end
+struct CoreBuddFriction#{{{
         H_input::ElementInput
         b_input::ElementInput
@@ -8 +10 @@
         g::Float64
 end# }}}
-function Friction(element::Tria) #{{{
+function CoreFriction(element::Tria) #{{{
 
-        H_input  = GetInput(element, ThicknessEnum)
-        b_input  = GetInput(element, BaseEnum)
-        c_input  = GetInput(element, FrictionCoefficientEnum)
+        frictionlaw = FindParam(element, FrictionLawEnum)
 
-        rho_ice   = FindParam(element, MaterialsRhoIceEnum)
-        rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
-        g         = FindParam(element, ConstantsGEnum)
+        if frictionlaw==1
+                H_input  = GetInput(element, ThicknessEnum)
+                b_input  = GetInput(element, BaseEnum)
+                c_input  = GetInput(element, FrictionCoefficientEnum)
 
-        return Friction(H_input, b_input, c_input, rho_ice, rho_water, g)
+                rho_ice   = FindParam(element, MaterialsRhoIceEnum)
+                rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
+                g         = FindParam(element, ConstantsGEnum)
+
+                return CoreBuddFriction(H_input, b_input, c_input, rho_ice, rho_water, g)
+        else
+                error("Friction ",typeof(md.friction)," not supported yet")
+        end
 end#}}}
 
 #vertices functions
-function Alpha2(friction::Friction, gauss::GaussTria, i::Int64) #{{{
+function Alpha2(friction::CoreBuddFriction, gauss::GaussTria, i::Int64) #{{{
 
         #Get effective pressure