Changeset 27165

Timestamp:

07/28/22 08:11:55 (3 years ago)

Author:

vverjans

Message:

CHG: allowing stochastic water pressure for frictioncoulomb

Location:

issm/trunk-jpl/src

Files:

• c/classes/Loads/Friction.cpp (modified) (1 diff)
• c/modules/StochasticForcingx/StochasticForcingx.cpp (modified) (3 diffs)
• c/shared/Enum/Enum.vim (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/Enumjl.vim (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (8 diffs)
• m/classes/stochasticforcing.m (modified) (3 diffs)
• m/classes/stochasticforcing.py (modified) (2 diffs)

issm/trunk-jpl/src/c/classes/Loads/Friction.cpp

@@ -344 +344 @@
 
         /*Get effective pressure*/
-        IssmDouble Neff = EffectivePressure(gauss);
+        bool ispwStochastic;
+        IssmDouble Neff;
+        element->parameters->FindParam(&ispwStochastic,StochasticForcingIsWaterPressureEnum);
+        if(ispwStochastic){
+                /*Retrieve stochastic water pressure and compute ice pressure*/
+                IssmDouble p_ice,p_water,Neff_limit;
+                element->GetInputValue(&p_water,gauss,FrictionCoulombWaterPressureEnum);
+                element->parameters->FindParam(&Neff_limit,FrictionEffectivePressureLimitEnum);
+                p_ice = IcePressure(gauss);
+                Neff  = max(Neff_limit*p_ice, p_ice - p_water);
+        }       
+        else{
+                /*Compute effective pressure directly*/
+                Neff = EffectivePressure(gauss);
+        }
+        
+        /*Get velocity magnitude*/
         IssmDouble vmag = VelMag(gauss);
 

issm/trunk-jpl/src/c/modules/StochasticForcingx/StochasticForcingx.cpp

@@ -207 +207 @@
                                                         gauss->GaussVertex(i);
                                                         p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
-                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
+                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; 
                                                 }
                                                 element->AddInput(FrictionWaterPressureEnum,p_water,P1DGEnum);
@@ -228 +228 @@
                                                         gauss->GaussVertex(i);
                                                         p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
-                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; //make sure positive (29Nov2021)
+                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; 
                                                 }
                                                 element->AddInput(FrictionSchoofWaterPressureEnum,p_water,P1DGEnum);
@@ -234 +234 @@
                                                 delete friction;
                                         }
-                                        break;                          default:
+                                        break;
+                                case FrictionCoulombWaterPressureEnum:
+                                        /*Specify that WaterPressure is stochastic*/ 
+                                        femmodel->parameters->SetParam(true,StochasticForcingIsWaterPressureEnum);
+                                        for(Object* &object:femmodel->elements->objects){
+                  Element* element = xDynamicCast<Element*>(object);
+                  int numvertices  = element->GetNumberOfVertices();
+                  IssmDouble p_water_deterministic[numvertices];
+                  IssmDouble p_water[numvertices];
+                                                element->GetInputValue(&dimensionid,StochasticForcingDefaultIdEnum);
+                                                Gauss* gauss=element->NewGauss();
+                                                Friction* friction = new Friction(element);
+                                                for(int i=0;i<numvertices;i++){
+                                                        gauss->GaussVertex(i);
+                                                        p_water_deterministic[i] = friction->SubglacialWaterPressure(gauss);
+                                                        p_water[i]               = p_water_deterministic[i] + noisefield[dimensionid]; 
+                                                }
+                                                element->AddInput(FrictionCoulombWaterPressureEnum,p_water,P1DGEnum);
+                                                delete gauss;
+                                                delete friction;
+                                        }
+                                        break;
+                                default:
                                         _error_("Field "<<EnumToStringx(fields[j])<<" does not support stochasticity yet.");
                         }

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

@@ -747 +747 @@
 syn keyword cConstant FrictionCoefficientEnum
 syn keyword cConstant FrictionCoefficientcoulombEnum
+syn keyword cConstant FrictionCoulombWaterPressureEnum
 syn keyword cConstant FrictionEffectivePressureEnum
 syn keyword cConstant FrictionMEnum

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

@@ -743 +743 @@
         FrictionCoefficientEnum,
         FrictionCoefficientcoulombEnum,
+        FrictionCoulombWaterPressureEnum,
         FrictionEffectivePressureEnum,
         FrictionMEnum,

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

@@ -749 +749 @@
                 case FrictionCoefficientEnum : return "FrictionCoefficient";
                 case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
+                case FrictionCoulombWaterPressureEnum : return "FrictionCoulombWaterPressure";
                 case FrictionEffectivePressureEnum : return "FrictionEffectivePressure";
                 case FrictionMEnum : return "FrictionM";

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

@@ -740 +740 @@
 syn keyword juliaConstC FrictionCoefficientEnum
 syn keyword juliaConstC FrictionCoefficientcoulombEnum
+syn keyword juliaConstC FrictionCoulombWaterPressureEnum
 syn keyword juliaConstC FrictionEffectivePressureEnum
 syn keyword juliaConstC FrictionMEnum

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

@@ -767 +767 @@
               else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
               else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
+              else if (strcmp(name,"FrictionCoulombWaterPressure")==0) return FrictionCoulombWaterPressureEnum;
               else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
               else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
@@ -874 +875 @@
               else if (strcmp(name,"RadarAttenuationWolff")==0) return RadarAttenuationWolffEnum;
               else if (strcmp(name,"RadarIcePeriod")==0) return RadarIcePeriodEnum;
-              else if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
+              if (strcmp(name,"RadarPowerMacGregor")==0) return RadarPowerMacGregorEnum;
+              else if (strcmp(name,"RadarPowerWolff")==0) return RadarPowerWolffEnum;
               else if (strcmp(name,"RheologyBAbsGradient")==0) return RheologyBAbsGradientEnum;
               else if (strcmp(name,"RheologyBInitialguess")==0) return RheologyBInitialguessEnum;
@@ -997 +998 @@
               else if (strcmp(name,"SmbDswdiffrf")==0) return SmbDswdiffrfEnum;
               else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum;
-              else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
+              if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
+              else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum;
               else if (strcmp(name,"SmbDzTop")==0) return SmbDzTopEnum;
               else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum;
               else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum;
-              else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
+              if (strcmp(name,"Thickness")==0) return ThicknessEnum;
+              else if (strcmp(name,"ThicknessOld")==0) return ThicknessOldEnum;
               else if (strcmp(name,"ThicknessPositive")==0) return ThicknessPositiveEnum;
               else if (strcmp(name,"ThicknessResidual")==0) return ThicknessResidualEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum;
               else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum;
-              else if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
+              if (strcmp(name,"Outputdefinition87")==0) return Outputdefinition87Enum;
+              else if (strcmp(name,"Outputdefinition88")==0) return Outputdefinition88Enum;
               else if (strcmp(name,"Outputdefinition89")==0) return Outputdefinition89Enum;
               else if (strcmp(name,"Outputdefinition8")==0) return Outputdefinition8Enum;
@@ -1366 +1367 @@
               else if (strcmp(name,"FloatingArea")==0) return FloatingAreaEnum;
               else if (strcmp(name,"FloatingAreaScaled")==0) return FloatingAreaScaledEnum;
-              else if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"Free")==0) return FreeEnum;
+              if (strcmp(name,"FloatingMeltRate")==0) return FloatingMeltRateEnum;
+              else if (strcmp(name,"Free")==0) return FreeEnum;
               else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum;
               else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum;
@@ -1489 +1490 @@
               else if (strcmp(name,"MeltingAnalysis")==0) return MeltingAnalysisEnum;
               else if (strcmp(name,"MeshElements")==0) return MeshElementsEnum;
-              else if (strcmp(name,"MeshX")==0) return MeshXEnum;
          else stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"MeshY")==0) return MeshYEnum;
+              if (strcmp(name,"MeshX")==0) return MeshXEnum;
+              else if (strcmp(name,"MeshY")==0) return MeshYEnum;
               else if (strcmp(name,"MinVel")==0) return MinVelEnum;
               else if (strcmp(name,"MinVx")==0) return MinVxEnum;
@@ -1612 +1613 @@
               else if (strcmp(name,"SubelementMelt2")==0) return SubelementMelt2Enum;
               else if (strcmp(name,"SubelementMigration")==0) return SubelementMigrationEnum;
-              else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
          else stage=14;
    }
    if(stage==14){
-              if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
+              if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
+              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
               else if (strcmp(name,"Tetra")==0) return TetraEnum;
               else if (strcmp(name,"TetraInput")==0) return TetraInputEnum;

issm/trunk-jpl/src/m/classes/stochasticforcing.m

@@ -99 +99 @@
                   error('md.friction does not agree with stochasticforcing field %s', char(field));
                end
-                                        if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof'))
+                                        if(strcmp(class(md.friction),'friction') || strcmp(class(md.friction),'frictionschoof') || strcmp(class(md.friction),'frictioncoulomb'))
                   if(md.friction.coupling~=0 && md.friction.coupling~=1 && md.friction.coupling~=2)
                      error('stochasticforcing field %s is only implemented for cases md.friction.coupling 0 or 1 or 2', char(field));
@@ -271 +271 @@
                 'FloatingMeltRate',...
                 'FrictionWaterPressure',...
+                'FrictionCoulombWaterPressure',...
                 'FrictionSchoofWaterPressure',...
                 'FrontalForcingsRignotAutoregression',...
@@ -282 +283 @@
                 'basalforcings',...
                 'friction',...
+                'frictioncoulomb',...
                 'frictionschoof',...
                 'frontalforcingsrignotautoregression',...

issm/trunk-jpl/src/m/classes/stochasticforcing.py

@@ -107 +107 @@
                 if (type(md.friction).__name__ != mdname):
                     raise TypeError('md.friction does not agree with stochasticforcing field {}'.format(field))
-                if (type(md.friction).__name__=='friction' or type(md.friction).__name__=='frictionschoof'):
+                if (type(md.friction).__name__=='friction' or type(md.friction).__name__=='frictionschoof' or type(md.friction).__name__=='frictioncoulomb'):
                     if md.friction.coupling not in[0, 1, 2]:
                         raise TypeError('stochasticforcing field {} is only implemented for cases md.friction.coupling 0 or 1 or 2'.format(field))
@@ -234 +234 @@
                      'FloatingMeltRate': 'basalforcings',
                      'FrictionWaterPressure': 'friction',
+                     'FrictionCoulombWaterPressure': 'frictioncoulomb',
                      'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotAutoregression': 'frontalforcingsrignotautoregression',