Changeset 27282

Timestamp:

09/14/22 12:59:00 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: preparing ground for GlaDS bigger domain

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (9 diffs)
• classes/Elements/Element.cpp (modified) (1 diff)
• classes/Loads/Channel.cpp (modified) (2 diffs)
• modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp (modified) (8 diffs)
• modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h (modified) (1 diff)
• shared/Enum/Enum.vim (modified) (1 diff)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/Enumjl.vim (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (8 diffs)

issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp ¶

@@ -233 +233 @@
 }/*}}}*/
 ElementMatrix* HydrologyGlaDSAnalysis::CreateKMatrix(Element* element){/*{{{*/
+
+        /*Skip if water or ice shelf element*/
+        if(element->IsAllFloating() || !element->IsIceInElement()) return NULL;
 
         /*Intermediaries */
@@ -332 +335 @@
 
         /*Skip if water or ice shelf element*/
-        if(element->IsAllFloating()) return NULL;
+        if(element->IsAllFloating() || !element->IsIceInElement()) return NULL;
 
         /*Intermediaries */
@@ -443 +446 @@
 }/*}}}*/
 void           HydrologyGlaDSAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-        SetActiveNodesLSMx(femmodel);
+
+        /*Update active elements based on ice levelset and ocean levelset*/
+        GetMaskOfIceVerticesLSMx(femmodel,true);
+        SetActiveNodesLSMx(femmodel,true);
+
+        IssmDouble rho_ice   = femmodel->parameters->FindParam(MaterialsRhoIceEnum);
+        IssmDouble rho_water = femmodel->parameters->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g         = femmodel->parameters->FindParam(ConstantsGEnum);
+
+        /*Constrain all nodes that are grounded and unconstrain the ones that float*/
+        for(Object* & object : femmodel->elements->objects){
+                Element    *element  = xDynamicCast<Element*>(object);
+                int         numnodes  = element->GetNumberOfNodes();
+                IssmDouble *mask      = xNew<IssmDouble>(numnodes);
+                IssmDouble *bed       = xNew<IssmDouble>(numnodes);
+                IssmDouble *thickness = xNew<IssmDouble>(numnodes);
+                IssmDouble *ls_active = xNew<IssmDouble>(numnodes);
+
+                element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum);
+                element->GetInputListOnNodes(&bed[0],BaseEnum);
+                element->GetInputListOnNodes(&thickness[0],BaseEnum);
+                element->GetInputListOnNodes(&ls_active[0],HydrologyMaskNodeActivationEnum);
+
+                for(int in=0;in<numnodes;in++){
+                        Node* node=element->GetNode(in);
+                        if(mask[in]>0. && ls_active[in]==1.){
+                                node->Activate(); //Not sure if we need this!
+                        }
+                        else{
+                                IssmDouble phi =  rho_ice*g*thickness[in] + rho_water*g*bed[in]; //FIXME this is correct!
+                                node->Deactivate();// Not sure if we need this
+                                node->ApplyConstraint(0,phi);
+                        }
+                }
+                xDelete<IssmDouble>(mask);
+                xDelete<IssmDouble>(bed);
+                xDelete<IssmDouble>(thickness);
+                xDelete<IssmDouble>(ls_active);
+        }
+
+        return;
 }/*}}}*/
 
@@ -471 +514 @@
 void HydrologyGlaDSAnalysis::UpdateSheetThickness(Element* element){/*{{{*/
 
-        /*Skip if water or ice shelf element*/
-        //if(element->IsAllFloating()) return;
-
         /*Intermediaries */
         IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r,H,b;
@@ -484 +524 @@
         /*Initialize new sheet thickness*/
         IssmDouble* h_new = xNew<IssmDouble>(numvertices);
+
+        /*Set to 0 if inactive element*/
+        if(element->IsAllFloating() || !element->IsIceInElement()){
+                for(int iv=0;iv<numvertices;iv++) h_new[iv] = 0.;
+                element->AddInput(HydrologySheetThicknessEnum,h_new,P1Enum);
+                xDelete<IssmDouble>(h_new);
+                return;
+        }
 
         /*Retrieve all inputs and parameters*/
@@ -545 +593 @@
         
         /*Force floating ice to have zero sheet thickness*/
-        if(element->IsAllFloating()){
+        if(element->IsAllFloating() || !element->IsIceInElement()){
                 for(int iv=0;iv<numvertices;iv++) h_new[iv] = 0.;
         }
@@ -564 +612 @@
 }/*}}}*/
 void HydrologyGlaDSAnalysis::UpdateEffectivePressure(Element* element){/*{{{*/
-
-        /*Skip if water or ice shelf element*/
-        //if(element->IsAllFloating()) return;
 
         /*Intermediary*/
@@ -583 +628 @@
         Input* phi_input = element->GetInput(HydraulicPotentialEnum); _assert_(phi_input);
 
+        /*Set to 0 if inactive element*/
+        if(element->IsAllFloating() || !element->IsIceInElement()){
+                for(int iv=0;iv<numnodes;iv++) N[iv] = 0.;
+                element->AddInput(EffectivePressureEnum,N,P1Enum);
+                xDelete<IssmDouble>(N);
+                return;
+        }
+
         /* Start  looping on the number of gaussian points: */
         Gauss* gauss=element->NewGauss();
@@ -607 +660 @@
                 if(xIsNan<IssmDouble>(N[iv])) _error_("NaN found in solution vector");
                 if(xIsInf<IssmDouble>(N[iv])) _error_("Inf found in solution vector");
-        }
-
-   /*Force floating ice to have zero effective pressure*/
-   if(element->IsAllFloating()){
-                for(int iv=0;iv<numnodes;iv++) N[iv] = 0.;
         }
 

issm/trunk-jpl/src/c/classes/Elements/Element.cpp ¶

@@ -2366 +2366 @@
 }/*}}}*/
 bool       Element::IsGrounded(){/*{{{*/
+        /*At least ONE node is grounded (partially grounded returns true)*/
 
         Input* input=this->GetInput(MaskOceanLevelsetEnum); _assert_(input);

issm/trunk-jpl/src/c/classes/Loads/Channel.cpp ¶

@@ -603 +603 @@
         /*Initialize Element matrix and return if necessary*/
         Tria*  tria=(Tria*)element;
-        if(!tria->IsIceInElement() || tria->IsAllFloating() || this->boundary){
+        if(this->boundary){
                 this->S = 0.;
                 return;
         }
         _assert_(tria->FiniteElement()==P1Enum); 
+
+        /*Evaluate all fields on center of edge*/
         int index1=tria->GetVertexIndex(vertices[0]);
         int index2=tria->GetVertexIndex(vertices[1]);
+        GaussTria* gauss=new GaussTria();
+        gauss->GaussEdgeCenter(index1,index2);
+
+        /*Set to 0 if inactive*/
+        IssmDouble active;
+        Input* active_input = element->GetInput(HydrologyMaskNodeActivationEnum); _assert_(active_input);
+        active_input->GetInputValue(&active,gauss);
+        if(active!=1.){
+                this->S = 0.;
+                delete gauss;
+                return;
+        }
 
         /*Intermediaries */
@@ -644 +658 @@
         tx = tx/Lt;
         ty = ty/Lt;
-
-        /*Evaluate fields on center of edge*/
-        GaussTria* gauss=new GaussTria();
-        gauss->GaussEdgeCenter(index1,index2);
 
         /*Get input values at gauss points*/

issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp ¶

@@ -10 +10 @@
 #include "../modules.h"
 
-void SetActiveNodesLSMx(FemModel* femmodel){/*{{{*/
+void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
         /* activate/deactivate nodes for levelset method according to IceMaskNodeActivation */
 
-        /* find parameters */
-        int domaintype;
-        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
+        /*Determine which node activation mask to pull from*/
+        int nodeactivationmask = IceMaskNodeActivationEnum;
+        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
+
 
         for(Object* & object : femmodel->elements->objects){
@@ -28 +29 @@
                         case TaylorHoodEnum:case XTaylorHoodEnum:case LATaylorHoodEnum:
                         case CrouzeixRaviartEnum:case LACrouzeixRaviartEnum:case OneLayerP4zEnum:{
-                                Input* input=element->GetInput(IceMaskNodeActivationEnum);
-                                if(!input) _error_("Input " << EnumToStringx(IceMaskNodeActivationEnum) << " not found in element");
+                                Input* input=element->GetInput(nodeactivationmask);
+                                if(!input) _error_("Input " << EnumToStringx(nodeactivationmask) << " not found in element");
 
                                 /* Start looping on the number of vertices: */
@@ -45 +46 @@
                         }
                         default:
-                                element->GetInputListOnNodes(&mask[0],IceMaskNodeActivationEnum);
+                                element->GetInputListOnNodes(&mask[0],nodeactivationmask);
                                 break;
                 }
@@ -58 +59 @@
 }/*}}}*/
 
-void GetMaskOfIceVerticesLSMx0(FemModel* femmodel){/*{{{*/
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology){/*{{{*/
+
+        /*Determine which node activation to construct*/
+        int nodeactivationmask = IceMaskNodeActivationEnum;
+        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
 
         /*Initialize vector with number of vertices*/
@@ -68 +73 @@
 
         /*Fill vector with values: */
-        for(Object* & object : femmodel->elements->objects){
-                Element* element=xDynamicCast<Element*>(object);
-                if(element->IsIceInElement()){
-                        int nbv = element->GetNumberOfVertices();
-                        for(int iv=0;iv<nbv;iv++){
-                                vec_mask_ice->SetValue(element->vertices[iv]->Pid(),1.,INS_VAL);
+        if(ishydrology){
+                for(Object* & object : femmodel->elements->objects){
+                        Element* element=xDynamicCast<Element*>(object);
+                        if(element->IsIceInElement()){
+                                int nbv = element->GetNumberOfVertices();
+                                for(int iv=0;iv<nbv;iv++){
+                                        vec_mask_ice->SetValue(element->vertices[iv]->Pid(),1.,INS_VAL);
+                                }
+                        }
+                }
+        }
+        else{
+                for(Object* & object : femmodel->elements->objects){
+                        Element* element=xDynamicCast<Element*>(object);
+                        if(element->IsIceInElement() && element->IsGrounded()){
+                                int nbv = element->GetNumberOfVertices();
+                                for(int iv=0;iv<nbv;iv++){
+                                        vec_mask_ice->SetValue(element->vertices[iv]->Pid(),1.,INS_VAL);
+                                }
                         }
                 }
@@ -80 +98 @@
         /*Assemble vector and serialize */
         vec_mask_ice->Assemble();
-        InputUpdateFromVectorx(femmodel,vec_mask_ice,IceMaskNodeActivationEnum,VertexPIdEnum);
+        InputUpdateFromVectorx(femmodel,vec_mask_ice,nodeactivationmask,VertexPIdEnum);
         delete vec_mask_ice;
 }/*}}}*/
-void GetMaskOfIceVerticesLSMx(FemModel* femmodel){/*{{{*/
+void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology){/*{{{*/
 
-        femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
+        /*Set configuration to levelset*/
+        if(ishydrology){
+                /*We may not be running with ismovingfront so we can't assume LevelsetAnalysis is active*/
+                femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
+        }
+        else{
+                femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum);
+        }
+
+        /*Determine which node activation to construct*/
+        int nodeactivationmask = IceMaskNodeActivationEnum;
+        if(ishydrology) nodeactivationmask = HydrologyMaskNodeActivationEnum;
 
         /*Create vector on gset*/
@@ -97 +126 @@
                 Element* element=xDynamicCast<Element*>(object);
 
-                if(element->IsIceInElement()){
-                        int numnodes = element->GetNumberOfNodes();
-                        int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
-                        int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
-                        IssmDouble* ones = xNew<IssmDouble>(gsize_local);
-                        for(int n=0;n<gsize_local;n++) ones[n] = 1.;
-                        vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
-                        xDelete<IssmDouble>(ones);
-                        xDelete<int>(glist_local);
+                if(ishydrology){
+                        if(element->IsIceInElement() && element->IsGrounded()){
+                                int numnodes = element->GetNumberOfNodes();
+                                int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                IssmDouble* ones = xNew<IssmDouble>(gsize_local);
+                                for(int n=0;n<gsize_local;n++) ones[n] = 1.;
+                                vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
+                                xDelete<IssmDouble>(ones);
+                                xDelete<int>(glist_local);
+                        }
+                }
+                else{
+                        if(element->IsIceInElement()){
+                                int numnodes = element->GetNumberOfNodes();
+                                int  gsize_local=GetNumberOfDofs(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                int* glist_local=GetGlobalDofList(element->nodes,numnodes,GsetEnum,NoneEnum);
+                                IssmDouble* ones = xNew<IssmDouble>(gsize_local);
+                                for(int n=0;n<gsize_local;n++) ones[n] = 1.;
+                                vec_mask_ice->SetValues(gsize_local,glist_local,ones,INS_VAL);
+                                xDelete<IssmDouble>(ones);
+                                xDelete<int>(glist_local);
+                        }
                 }
         }
@@ -120 +163 @@
         for(Object* & object : femmodel->elements->objects){
                 Element* element=xDynamicCast<Element*>(object);
-                element->InputUpdateFromSolutionOneDof(local_ug,IceMaskNodeActivationEnum);
+                element->InputUpdateFromSolutionOneDof(local_ug,nodeactivationmask);
         }
 

issm/trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.h ¶

@@ -8 +8 @@
 #include "../../classes/classes.h"
 
-void SetActiveNodesLSMx(FemModel* femmodel);
-void GetMaskOfIceVerticesLSMx0(FemModel* femmodel);
-void GetMaskOfIceVerticesLSMx(FemModel* femmodel);
+void SetActiveNodesLSMx(FemModel* femmodel,bool ishydrology=false);
+void GetMaskOfIceVerticesLSMx0(FemModel* femmodel,bool ishydrology=false);
+void GetMaskOfIceVerticesLSMx(FemModel* femmodel,bool ishydrology=false);
 #endif  /* _SETACTIVENODESLSMX_H*/

issm/trunk-jpl/src/c/shared/Enum/Enum.vim ¶

@@ -820 +820 @@
 syn keyword cConstant HydrologyWaterVxEnum
 syn keyword cConstant HydrologyWaterVyEnum
+syn keyword cConstant HydrologyMaskNodeActivationEnum
 syn keyword cConstant IceEnum
 syn keyword cConstant IceMaskNodeActivationEnum

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h ¶

@@ -816 +816 @@
         HydrologyWaterVxEnum,
         HydrologyWaterVyEnum,
+        HydrologyMaskNodeActivationEnum,
         IceEnum,
         IceMaskNodeActivationEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp ¶

@@ -822 +822 @@
                 case HydrologyWaterVxEnum : return "HydrologyWaterVx";
                 case HydrologyWaterVyEnum : return "HydrologyWaterVy";
+                case HydrologyMaskNodeActivationEnum : return "HydrologyMaskNodeActivation";
                 case IceEnum : return "Ice";
                 case IceMaskNodeActivationEnum : return "IceMaskNodeActivation";

issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim ¶

@@ -813 +813 @@
 syn keyword juliaConstC HydrologyWaterVxEnum
 syn keyword juliaConstC HydrologyWaterVyEnum
+syn keyword juliaConstC HydrologyMaskNodeActivationEnum
 syn keyword juliaConstC IceEnum
 syn keyword juliaConstC IceMaskNodeActivationEnum

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp ¶

@@ -840 +840 @@
               else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
               else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
+              else if (strcmp(name,"HydrologyMaskNodeActivation")==0) return HydrologyMaskNodeActivationEnum;
               else if (strcmp(name,"Ice")==0) return IceEnum;
               else if (strcmp(name,"IceMaskNodeActivation")==0) return IceMaskNodeActivationEnum;
@@ -874 +875 @@
               else if (strcmp(name,"Misfit")==0) return MisfitEnum;
               else if (strcmp(name,"MovingFrontalVx")==0) return MovingFrontalVxEnum;
-              else if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
+              if (strcmp(name,"MovingFrontalVy")==0) return MovingFrontalVyEnum;
+              else if (strcmp(name,"Neumannflux")==0) return NeumannfluxEnum;
               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
               else if (strcmp(name,"Node")==0) return NodeEnum;
@@ -997 +998 @@
               else if (strcmp(name,"SmbCciceValue")==0) return SmbCciceValueEnum;
               else if (strcmp(name,"SmbCotValue")==0) return SmbCotValueEnum;
-              else if (strcmp(name,"SmbD")==0) return SmbDEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
+              if (strcmp(name,"SmbD")==0) return SmbDEnum;
+              else if (strcmp(name,"SmbDailyairdensity")==0) return SmbDailyairdensityEnum;
               else if (strcmp(name,"SmbDailyairhumidity")==0) return SmbDailyairhumidityEnum;
               else if (strcmp(name,"SmbDailydlradiation")==0) return SmbDailydlradiationEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"Surface")==0) return SurfaceEnum;
               else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum;
-              else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
+              if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum;
+              else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum;
               else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum;
               else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"Outputdefinition69")==0) return Outputdefinition69Enum;
               else if (strcmp(name,"Outputdefinition6")==0) return Outputdefinition6Enum;
-              else if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
+              if (strcmp(name,"Outputdefinition70")==0) return Outputdefinition70Enum;
+              else if (strcmp(name,"Outputdefinition71")==0) return Outputdefinition71Enum;
               else if (strcmp(name,"Outputdefinition72")==0) return Outputdefinition72Enum;
               else if (strcmp(name,"Outputdefinition73")==0) return Outputdefinition73Enum;
@@ -1366 +1367 @@
               else if (strcmp(name,"ElementHook")==0) return ElementHookEnum;
               else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
-              else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
+              if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum;
+              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
               else if (strcmp(name,"EsaSolution")==0) return EsaSolutionEnum;
               else if (strcmp(name,"EsaTransitions")==0) return EsaTransitionsEnum;
@@ -1489 +1490 @@
               else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
               else if (strcmp(name,"MasstransportAnalysis")==0) return MasstransportAnalysisEnum;
-              else if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
          else stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
+              if (strcmp(name,"MasstransportSolution")==0) return MasstransportSolutionEnum;
+              else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum;
               else if (strcmp(name,"Matenhancedice")==0) return MatenhancediceEnum;
               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
@@ -1612 +1613 @@
               else if (strcmp(name,"SoftMigration")==0) return SoftMigrationEnum;
               else if (strcmp(name,"SpatialLinearFloatingMeltRate")==0) return SpatialLinearFloatingMeltRateEnum;
-              else if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
          else stage=14;
    }
    if(stage==14){
-              if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
+              if (strcmp(name,"SpcDynamic")==0) return SpcDynamicEnum;
+              else if (strcmp(name,"SpcStatic")==0) return SpcStaticEnum;
               else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum;
               else if (strcmp(name,"Sset")==0) return SsetEnum;