Changeset 27476

Timestamp:

12/21/22 13:46:29 (2 years ago)

Author:

eyhli

Message:

ADD: Regularised Coulomb friction law (Joughin et al. 2019)

Location:

issm/trunk-jpl

Files:

• src/c/analyses/AdjointHorizAnalysis.cpp (modified) (5 diffs)
• src/c/classes/Loads/Friction.cpp (modified) (6 diffs)
• src/c/classes/Loads/Friction.h (modified) (2 diffs)
• src/c/modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.cpp (modified) (2 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• test/Archives/Archive358.arch (added)
• test/NightlyRun/test358.m (added)

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

@@ -2391 +2391 @@
                 case 2:
                 case 11:
+                case 14:
                         dragcoefficient_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoefficient_input);
                         break;
@@ -2484 +2485 @@
                 case 2:
                 case 11:
+                case 14:
                         dragcoeff_input = element->GetInput(FrictionCEnum); _assert_(dragcoeff_input);
                         break;
@@ -2593 +2595 @@
                 case 2:
                 case 11:
+                case 14:
                         dragcoeff_input = element->GetInput(FrictionCEnum); _assert_(dragcoeff_input);
                         break;
@@ -2714 +2717 @@
                 case 2:
                 case 11:
+                case 14:
                         dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input);
                         break;
@@ -2832 +2836 @@
                 case 2:
                 case 11:
+                case 14:
                         dragcoeff_input = basalelement->GetInput(FrictionCEnum); _assert_(dragcoeff_input);
                         break;

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

@@ -125 +125 @@
                         GetAlphaSchoofComplement(palpha_complement,gauss);
                         break;
+                case 14:
+                        GetAlphaRegCoulombComplement(palpha_complement,gauss);
+                        break;
           default:
                         _error_("not supported");
@@ -295 +298 @@
         *palpha_complement=alpha_complement;
 }/*}}}*/
+void Friction::GetAlphaRegCoulombComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
+
+        /* Compute the complement of Schoof's law for inversion
+         * d alpha2                       
+         * -------- = |u_b|^(1/m-1) * (|u_b + u_0|)^(-1/m)
+         *  dC                           
+         */
+
+        /*diverse: */
+        IssmDouble  m, u0;
+        IssmDouble  alpha_complement;
+
+        /*Recover parameters: */
+        element->GetInputValue(&m,gauss,FrictionMEnum);
+        element->parameters->FindParam(&u0,FrictionU0Enum);
+
+        /*Get velocity magnitude*/
+        IssmDouble ub = VelMag(gauss);
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(ub==0.) {
+                alpha_complement=0.;
+        }
+        else {
+                /*Compute friction complement*/
+                alpha_complement= (pow(ub,1./m-1.)) / pow(ub/u0 + 1.,1./m);     
+        }
+
+        /*Assign output pointers:*/
+        *palpha_complement=alpha_complement;
+}/*}}}*/
 void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
@@ -341 +375 @@
                 case 13:
                         GetAlpha2Coulomb2(palpha2,gauss);
+                        break;
+                case 14:
+                        GetAlpha2RegCoulomb(palpha2,gauss);
                         break;
           default:
@@ -910 +947 @@
         *palpha2=alpha2;
 }/*}}}*/
-
+void Friction::GetAlpha2RegCoulomb(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient
+         *
+         *               C |u_b|^(1/m-1)
+         * alpha2= __________________________
+         *          (|u_b|/u0 + 1 )^(1/m)
+         *
+         * */
+
+        /*diverse: */
+        IssmDouble  C,coeff,u0,m,alpha2;
+
+        /*Recover parameters: */
+        element->GetInputValue(&coeff,gauss,FrictionCEnum);
+        element->GetInputValue(&m,gauss,FrictionMEnum);
+        element->parameters->FindParam(&u0,FrictionU0Enum);
+
+        /* scale C for a better inversion */
+        C = coeff*coeff;
+
+        /*Get velocity magnitude*/
+        IssmDouble ub = VelMag(gauss);
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(ub==0.) {
+                alpha2=0.;
+        }
+        else {
+                /*Compute alpha^2*/
+                alpha2= (C*pow(ub,1./m-1.)) / pow(ub/u0 + 1.,1./m);
+        }
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
 IssmDouble Friction::EffectivePressure(Gauss* gauss){/*{{{*/
         /*Get effective pressure as a function of  flag */
@@ -1192 +1264 @@
                         }
                         break;
+                case 14:
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.C",FrictionCEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
+                        break;
                 default:
                         _error_("friction law "<< frictionlaw <<" not supported");
@@ -1264 +1340 @@
                         parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
                         break;
+                case 14:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.u0",FrictionU0Enum));
+                        break;
                 default: _error_("Friction law "<<frictionlaw<<" not implemented yet");
         }

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

@@ -41 +41 @@
                 void  GetAlphaViscousComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlphaSchoofComplement(IssmDouble* alpha_complement,Gauss* gauss);
+                void  GetAlphaRegCoulombComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlphaWeertmanComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
@@ -55 +56 @@
                 void  GetAlpha2PISM(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Schoof(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2RegCoulomb(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Tsai(IssmDouble* palpha2,Gauss* gauss);
 

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

@@ -63 +63 @@
         DatasetInput* weights_input=basalelement->GetDatasetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
 
-        /* get the friction law: if 11-Schoof, which has a different name of C */
+        /* get the friction law: if 2-Weertman, 11-Schoof or 14-RegularizedCoulomb, which has a different names of C */
         element->FindParam(&frictionlaw, FrictionLawEnum);
         Input* drag_input;
@@ -69 +69 @@
                 case 2:
                 case 11:
+                case 14:
                         drag_input = basalelement->GetInput(FrictionCEnum); _assert_(drag_input);
                         break;

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

@@ -205 +205 @@
         FrictionLinearizeEnum,
         FrictionPseudoplasticityExponentEnum,
+        FrictionU0Enum,
         FrictionThresholdSpeedEnum,
         FrictionVoidRatioEnum,