Changeset 24139

Timestamp:

09/11/19 06:18:57 (6 years ago)

Author:

rueckamp

Message:

CHG: added ThicknessResidual as possible output to MasstransportAnalysis

Location:

issm/trunk-jpl/src/c

Files:

• analyses/MasstransportAnalysis.cpp (modified) (6 diffs)
• shared/Enum/EnumDefinitions.h (modified) (1 diff)
• shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• shared/Enum/StringToEnumx.cpp (modified) (5 diffs)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -157 +157 @@
         /*Initialize cumdeltalthickness input*/
         InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
+        /*Initialize ThicknessResidual input*/
+        InputUpdateFromConstantx(elements,0.,ThicknessResidualEnum);
 
         /*Get what we need for ocean-induced basal melting*/
@@ -851 +853 @@
         int numnodes = element->GetNumberOfNodes();
         IssmDouble* newthickness = xNew<IssmDouble>(numnodes);
+        IssmDouble* thicknessresidual = xNew<IssmDouble>(numnodes);
 
         /*Use the dof list to index into the solution vector: */
-   IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
+        IssmDouble minthickness = element->FindParam(MasstransportMinThicknessEnum);
         for(int i=0;i<numnodes;i++){
                 newthickness[i]=solution[doflist[i]];
+                thicknessresidual[i]=0.;        
                 /*Check solution*/
                 if(xIsNan<IssmDouble>(newthickness[i])) _error_("NaN found in solution vector");
                 if(xIsInf<IssmDouble>(newthickness[i])) _error_("Inf found in solution vector");
-      if(newthickness[i]<minthickness) newthickness[i]=minthickness;
+                if(newthickness[i]<minthickness){
+                        thicknessresidual[i]=minthickness-newthickness[i];
+                        newthickness[i]=minthickness;
+                }
         }
         element->AddBasalInput(ThicknessEnum,newthickness,element->GetElementType());
+        element->AddBasalInput(ThicknessResidualEnum,thicknessresidual,element->GetElementType());
+        
         xDelete<int>(doflist);
         xDelete<IssmDouble>(newthickness);
+        xDelete<IssmDouble>(thicknessresidual);
 
         /*Update bed and surface accordingly*/
@@ -869 +879 @@
         /*Get basal element*/
         int domaintype; element->FindParam(&domaintype,DomainTypeEnum);
-   Element* basalelement=element;
+        Element* basalelement=element;
         if(domaintype!=Domain2DhorizontalEnum) basalelement = element->SpawnBasalElement();
 
         /*Fetch number of nodes and dof for this finite element*/
-   int numvertices = basalelement->GetNumberOfVertices();
+        int numvertices = basalelement->GetNumberOfVertices();
 
         /*Now, we need to do some "processing"*/
-   newthickness  = xNew<IssmDouble>(numvertices);
+        newthickness  = xNew<IssmDouble>(numvertices);
         IssmDouble* oldthickness      = xNew<IssmDouble>(numvertices);
         IssmDouble* cumdeltathickness = xNew<IssmDouble>(numvertices);
@@ -889 +899 @@
 
         /*Get previous base, thickness, surfac and current sealevel and bed:*/
-   basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
+        basalelement->GetInputListOnVertices(&newthickness[0],ThicknessEnum);
         basalelement->GetInputListOnVertices(&oldthickness[0],ThicknessOldEnum);
         basalelement->GetInputListOnVertices(&oldbase[0],BaseEnum);
@@ -897 +907 @@
         basalelement->GetInputListOnVertices(&cumdeltathickness[0],SealevelriseCumDeltathicknessEnum);
 
-   /*Do we do grounding line migration?*/
+        /*Do we do grounding line migration?*/
         bool isgroundingline;
         element->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
@@ -909 +919 @@
 
         /*Find MasstransportHydrostaticAdjustment to figure out how to update the geometry:*/
-   int hydroadjustment;
+        int hydroadjustment;
         basalelement->FindParam(&hydroadjustment,MasstransportHydrostaticAdjustmentEnum);
         IssmDouble rho_ice   = basalelement->FindParam(MaterialsRhoIceEnum);

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

@@ -779 +779 @@
         ThicknessOldEnum,
         ThicknessPositiveEnum,
+        ThicknessResidualEnum,
         VelEnum,
         VxAverageEnum,

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

@@ -785 +785 @@
                 case ThicknessOldEnum : return "ThicknessOld";
                 case ThicknessPositiveEnum : return "ThicknessPositive";
+                case ThicknessResidualEnum : return "ThicknessResidual";
                 case VelEnum : return "Vel";
                 case VxAverageEnum : return "VxAverage";

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

@@ -803 +803 @@
               else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
               else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
+              else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
               else if (strcmp(name,"Vel")==0) return VelEnum;
               else if (strcmp(name,"VxAverage")==0) return VxAverageEnum;
@@ -874 +875 @@
               else if (strcmp(name,"Outputdefinition54")==0) return Outputdefinition54Enum;
               else if (strcmp(name,"Outputdefinition55")==0) return Outputdefinition55Enum;
-              else if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
+              if (strcmp(name,"Outputdefinition56")==0) return Outputdefinition56Enum;
+              else if (strcmp(name,"Outputdefinition57")==0) return Outputdefinition57Enum;
               else if (strcmp(name,"Outputdefinition58")==0) return Outputdefinition58Enum;
               else if (strcmp(name,"Outputdefinition59")==0) return Outputdefinition59Enum;
@@ -997 +998 @@
               else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum;
               else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum;
-              else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
+              if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum;
+              else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum;
               else if (strcmp(name,"DoubleMatParam")==0) return DoubleMatParamEnum;
               else if (strcmp(name,"DoubleParam")==0) return DoubleParamEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"Materials")==0) return MaterialsEnum;
               else if (strcmp(name,"Matestar")==0) return MatestarEnum;
-              else if (strcmp(name,"Matice")==0) return MaticeEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
+              if (strcmp(name,"Matice")==0) return MaticeEnum;
+              else if (strcmp(name,"Matlitho")==0) return MatlithoEnum;
               else if (strcmp(name,"MatrixParam")==0) return MatrixParamEnum;
               else if (strcmp(name,"MaxAbsVx")==0) return MaxAbsVxEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"StringArrayParam")==0) return StringArrayParamEnum;
               else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
-              else if (strcmp(name,"StringParam")==0) return StringParamEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
+              if (strcmp(name,"StringParam")==0) return StringParamEnum;
+              else if (strcmp(name,"SubelementFriction1")==0) return SubelementFriction1Enum;
               else if (strcmp(name,"SubelementFriction2")==0) return SubelementFriction2Enum;
               else if (strcmp(name,"SubelementMelt1")==0) return SubelementMelt1Enum;