Changeset 27508

Timestamp:

01/10/23 15:49:52 (2 years ago)

Author:

youngmc3

Message:

ADD: Regularised Coulomb friction law 2

Location:

issm/trunk-jpl/src/c

Files:

• classes/Loads/Friction.cpp (modified) (4 diffs)
• classes/Loads/Friction.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) (9 diffs)

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

@@ -400 +400 @@
                         case 14:
                                 GetAlpha2RegCoulomb(palpha2,gauss);
+                                break;
+                        case 15:
+                                GetAlpha2RegCoulomb2(palpha2,gauss);
                                 break;
                         default:
@@ -1015 +1018 @@
         *palpha2=alpha2;
 }/*}}}*/
+void Friction::GetAlpha2RegCoulomb2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient
+         *
+         *               C N |u_b|^(1/m-1)
+         * alpha2= __________________________
+         *          (|u_b| + K N^m )^(1/m)
+         *
+         * */
+
+        /*diverse: */
+        IssmDouble  C,K,m,alpha2;
+
+        /*Recover parameters: */
+        element->GetInputValue(&C,gauss,FrictionCEnum);
+        element->GetInputValue(&m,gauss,FrictionMEnum);
+        element->GetInputValue(&K,gauss,FrictionKEnum);
+
+        /*Get velocity magnitude*/
+        IssmDouble ub = VelMag(gauss);
+        IssmDouble Neff = EffectivePressure(gauss);
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(ub==0. && (m-1.)<0) {
+                alpha2=0.;
+        }
+        else {
+                /*Compute alpha^2*/
+                alpha2= (C*pow(ub,1./m-1.)) * Neff / pow((ub+pow(K*Neff,m)),1./m);
+        }
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
 IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
         /*Get effective pressure as a function of  flag */
@@ -1300 +1337 @@
                         iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
                         break;
+                case 15:
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.K",FrictionKEnum);
+                        break;
                 default:
                         _error_("friction law "<< frictionlaw <<" not supported");
@@ -1376 +1418 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.friction.u0",FrictionU0Enum));
                         break;
+                case 15:
+                        parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+                        break;
                 default: _error_("Friction law "<<frictionlaw<<" not implemented yet");
         }

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

@@ -57 +57 @@
                 void  GetAlpha2Schoof(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2RegCoulomb(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2RegCoulomb2(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Tsai(IssmDouble* palpha2,Gauss* gauss);
 

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

@@ -819 +819 @@
 syn keyword cConstant FrictionCoefficientcoulombEnum
 syn keyword cConstant FrictionEffectivePressureEnum
+syn keyword cConstant FrictionKEnum
 syn keyword cConstant FrictionMEnum
 syn keyword cConstant FrictionPEnum

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

@@ -815 +815 @@
         FrictionCoefficientcoulombEnum,
         FrictionEffectivePressureEnum,
+        FrictionKEnum,
         FrictionMEnum,
         FrictionPEnum,

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

@@ -821 +821 @@
                 case FrictionCoefficientcoulombEnum : return "FrictionCoefficientcoulomb";
                 case FrictionEffectivePressureEnum : return "FrictionEffectivePressure";
+                case FrictionKEnum : return "FrictionK";
                 case FrictionMEnum : return "FrictionM";
                 case FrictionPEnum : return "FrictionP";

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

@@ -812 +812 @@
 syn keyword juliaConstC FrictionCoefficientcoulombEnum
 syn keyword juliaConstC FrictionEffectivePressureEnum
+syn keyword juliaConstC FrictionKEnum
 syn keyword juliaConstC FrictionMEnum
 syn keyword juliaConstC FrictionPEnum

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

@@ -839 +839 @@
               else if (strcmp(name,"FrictionCoefficientcoulomb")==0) return FrictionCoefficientcoulombEnum;
               else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
+              else if (strcmp(name,"FrictionK")==0) return FrictionKEnum;
               else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
               else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
@@ -874 +875 @@
               else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum;
               else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
-              else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
+              if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
+              else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
               else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum;
               else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
@@ -997 +998 @@
               else if (strcmp(name,"BslcIce")==0) return BslcIceEnum;
               else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
-              else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
+              if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
+              else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
               else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
               else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"SmbPrecipitationsPresentday")==0) return SmbPrecipitationsPresentdayEnum;
               else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum;
-              else if (strcmp(name,"SmbRain")==0) return SmbRainEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"SmbRe")==0) return SmbReEnum;
+              if (strcmp(name,"SmbRain")==0) return SmbRainEnum;
+              else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
               else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
               else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"WaterfractionDrainage")==0) return WaterfractionDrainageEnum;
               else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
-              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
+              if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
+              else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
               else if (strcmp(name,"WaterPressureArmaPerturbation")==0) return WaterPressureArmaPerturbationEnum;
               else if (strcmp(name,"WaterPressureValuesAutoregression")==0) return WaterPressureValuesAutoregressionEnum;
@@ -1366 +1367 @@
               else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
               else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
-              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
          else stage=12;
    }
    if(stage==12){
-              if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
+              if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
+              else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
               else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
               else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
@@ -1489 +1490 @@
               else if (strcmp(name,"GroundedAreaScaled")==0) return GroundedAreaScaledEnum;
               else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
-              else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
          else stage=13;
    }
    if(stage==13){
-              if (strcmp(name,"Gset")==0) return GsetEnum;
+              if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
+              else if (strcmp(name,"Gset")==0) return GsetEnum;
               else if (strcmp(name,"Gsl")==0) return GslEnum;
               else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
@@ -1612 +1613 @@
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
               else if (strcmp(name,"None")==0) return NoneEnum;
-              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
          else stage=14;
    }
    if(stage==14){
-              if (strcmp(name,"NyeCO2")==0) return NyeCO2Enum;
+              if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
+              else if (strcmp(name,"NyeCO2")==0) return NyeCO2Enum;
               else if (strcmp(name,"NyeH2O")==0) return NyeH2OEnum;
               else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
@@ -1735 +1736 @@
               else if (strcmp(name,"TransientInput")==0) return TransientInputEnum;
               else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
-              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
          else stage=15;
    }
    if(stage==15){
-              if (strcmp(name,"Tria")==0) return TriaEnum;
+              if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
+              else if (strcmp(name,"Tria")==0) return TriaEnum;
               else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
               else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;