Changeset 23897

Timestamp:

04/25/19 10:43:54 (6 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added Tsai's friction law

Location:

issm/trunk-jpl/src

Files:

• c/analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• c/classes/Loads/Friction.cpp (modified) (2 diffs)
• c/classes/Loads/Friction.h (modified) (1 diff)
• c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• c/shared/Enum/StringToEnumx.cpp (modified) (7 diffs)
• m/classes/frictiontsai.m (modified) (1 diff)

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

@@ -875 +875 @@
                         iomodel->FetchDataToInput(elements,"md.friction.Cmax",FrictionCmaxEnum);
                         break;
+                case 12:
+                        iomodel->FetchDataToInput(elements,"md.friction.m",FrictionMEnum);
+                        iomodel->FetchDataToInput(elements,"md.friction.C",FrictionCEnum);
+                        iomodel->FetchDataToInput(elements,"md.friction.f",FrictionfEnum);
+                        break;
                 default:
                         _error_("friction law "<< frictionlaw <<" not supported");
@@ -972 +977 @@
                         parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
                         break;
+                case 12:
+                        parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
+                        break;
                 default: _error_("Friction law "<<frictionlaw<<" not implemented yet");
         }

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

@@ -197 +197 @@
                 case 11:
                         GetAlpha2Schoof(palpha2,gauss);
+                        break;
+                case 12:
+                        GetAlpha2Tsai(palpha2,gauss);
                         break;
           default:
@@ -629 +632 @@
         else{
                 alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
+        }
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlpha2Tsai(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient 
+         *
+         * alpha2= min(C |ub|^m , f N ) / |ub|
+         *
+         * */
+
+        /*diverse: */
+        IssmDouble  C,f,m,alpha2;
+
+        /*Recover parameters: */
+        element->GetInputValue(&f,gauss,FrictionfEnum);
+        element->GetInputValue(&C,gauss,FrictionCEnum);
+        element->GetInputValue(&m,FrictionMEnum);
+
+        /*Get effective pressure and velocity magnitude*/
+        IssmDouble N  = EffectivePressure(gauss);
+        IssmDouble ub = VelMag(gauss);
+
+        /*Compute alpha^2*/
+        if(ub<1e-10){
+                alpha2 = 0.;
+        }
+        else{
+                alpha2= C*pow(ub,m);
+
+                if(alpha2>f*N) alpha2 = f*N;
+
+                alpha2 = alpha2/ub;
         }
 

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

@@ -42 +42 @@
                 void  GetAlpha2PISM(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Schoof(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2Tsai(IssmDouble* palpha2,Gauss* gauss);
 
                 IssmDouble EffectivePressure(Gauss* gauss);

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

@@ -517 +517 @@
         FrictionCoefficientEnum,
         FrictionEffectivePressureEnum,
+        FrictionfEnum,
         FrictionMEnum,
         FrictionPEnum,

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

@@ -523 +523 @@
                 case FrictionCoefficientEnum : return "FrictionCoefficient";
                 case FrictionEffectivePressureEnum : return "FrictionEffectivePressure";
+                case FrictionfEnum : return "Frictionf";
                 case FrictionMEnum : return "FrictionM";
                 case FrictionPEnum : return "FrictionP";

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

@@ -535 +535 @@
               else if (strcmp(name,"FrictionCoefficient")==0) return FrictionCoefficientEnum;
               else if (strcmp(name,"FrictionEffectivePressure")==0) return FrictionEffectivePressureEnum;
+              else if (strcmp(name,"Frictionf")==0) return FrictionfEnum;
               else if (strcmp(name,"FrictionM")==0) return FrictionMEnum;
               else if (strcmp(name,"FrictionP")==0) return FrictionPEnum;
@@ -628 +629 @@
               else if (strcmp(name,"SealevelUEsa")==0) return SealevelUEsaEnum;
               else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum;
-              else if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
+              if (strcmp(name,"SealevelriseSpcthickness")==0) return SealevelriseSpcthicknessEnum;
+              else if (strcmp(name,"SealevelriseStericRate")==0) return SealevelriseStericRateEnum;
               else if (strcmp(name,"SealevelNEsa")==0) return SealevelNEsaEnum;
               else if (strcmp(name,"SealevelUGia")==0) return SealevelUGiaEnum;
@@ -751 +752 @@
               else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum;
               else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum;
-              else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
          else stage=7;
    }
    if(stage==7){
-              if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
+              if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum;
+              else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum;
               else if (strcmp(name,"SurfaceAbsMisfit")==0) return SurfaceAbsMisfitEnum;
               else if (strcmp(name,"SurfaceAbsVelMisfit")==0) return SurfaceAbsVelMisfitEnum;
@@ -874 +875 @@
               else if (strcmp(name,"Outputdefinition76")==0) return Outputdefinition76Enum;
               else if (strcmp(name,"Outputdefinition77")==0) return Outputdefinition77Enum;
-              else if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
          else stage=8;
    }
    if(stage==8){
-              if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
+              if (strcmp(name,"Outputdefinition78")==0) return Outputdefinition78Enum;
+              else if (strcmp(name,"Outputdefinition79")==0) return Outputdefinition79Enum;
               else if (strcmp(name,"Outputdefinition7")==0) return Outputdefinition7Enum;
               else if (strcmp(name,"Outputdefinition80")==0) return Outputdefinition80Enum;
@@ -997 +998 @@
               else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum;
               else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum;
-              else if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
          else stage=9;
    }
    if(stage==9){
-              if (strcmp(name,"Fset")==0) return FsetEnum;
+              if (strcmp(name,"FSApproximation")==0) return FSApproximationEnum;
+              else if (strcmp(name,"Fset")==0) return FsetEnum;
               else if (strcmp(name,"FSpressure")==0) return FSpressureEnum;
               else if (strcmp(name,"FSSolver")==0) return FSSolverEnum;
@@ -1120 +1121 @@
               else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
               else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
-              else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
          else stage=10;
    }
    if(stage==10){
-              if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
+              if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum;
+              else if (strcmp(name,"NoMeltOnPartiallyFloating")==0) return NoMeltOnPartiallyFloatingEnum;
               else if (strcmp(name,"NoneApproximation")==0) return NoneApproximationEnum;
               else if (strcmp(name,"None")==0) return NoneEnum;
@@ -1243 +1244 @@
               else if (strcmp(name,"ViscousHeating")==0) return ViscousHeatingEnum;
               else if (strcmp(name,"Water")==0) return WaterEnum;
-              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
          else stage=11;
    }
    if(stage==11){
-              if (strcmp(name,"XY")==0) return XYEnum;
+              if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
+              else if (strcmp(name,"XY")==0) return XYEnum;
               else if (strcmp(name,"XYZ")==0) return XYZEnum;
               else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;

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

@@ -15 +15 @@
                                 case 0
                                         self=setdefaultparameters(self);
+                                case 1
+                                        self=structtoobj(frictiontsai(),varargin{1});
                                 otherwise
                                         error('constructor not supported');