Changeset 22995

Timestamp:

07/22/18 12:06:19 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on new friction law for Henning

Location:

issm/trunk-jpl/src/c/classes/Loads

Files:

• Friction.cpp (modified) (3 diffs)
• Friction.h (modified) (1 diff)

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

@@ -195 +195 @@
         /*Check to prevent dividing by zero if vmag==0*/
         if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
-        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan<IssmDouble>(alpha_complement));
+        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
+        
+        _assert_(!xIsNan<IssmDouble>(alpha_complement));
+        _assert_(!xIsInf<IssmDouble>(alpha_complement));
 
         /*Assign output pointers:*/
@@ -230 +233 @@
                 case 9:
                         GetAlpha2Josh(palpha2,gauss);
+                        break;
+                case 10:
+                        GetAlpha2PISM(palpha2,gauss);
                         break;
           default:
@@ -666 +672 @@
         element->parameters->FindParam(&gamma,FrictionGammaEnum);
         alpha2 = alpha2 / exp((T-Tpmp)/gamma);
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlpha2PISM(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+        /*Here, we want to parameterize the friction as
+         *
+         * ... more comments to come
+         */
+
+        /*Compute effective pressure first */
+        IssmDouble  N,delta,W,Wmax,e0,Cc,P0;
+        //element->parameters->FindParam(&delta,TimeEnum);
+        //element->parameters->FindParam(&e0,TimeEnum);
+        //element->parameters->FindParam(&P0,TimeEnum);
+        //element->GetInputValue(&Cc,gauss,FrictionPressureAdjustedTemperatureEnum);
+        //element->GetInputValue(&W,gauss,FrictionPressureAdjustedTemperatureEnum);
+        N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));
+
+        /*Compute yield stress following a Mohr-Colomb criterion*/
+        IssmDouble phi;
+        //element->GetInputValue(&phi,gauss,FrictionPressureAdjustedTemperatureEnum);
+        IssmDouble tau_c = N*tan(phi);
+
+        /*Compute basal speed*/
+        IssmDouble ub;
+        element->GetInputValue(&ub,gauss,VelEnum);
+
+        /*now compute alpha^2*/
+        IssmDouble u0,q;
+        //element->parameters->FindParam(&u0,TimeEnum);
+        //element->parameters->FindParam(&q,TimeEnum);
+        IssmDouble alpha2 = tau_c/(pow(ub,1.-q)*pow(u0,q));
+
+        /*Final checks in debuging mode*/
+        _assert_(!xIsNan<IssmDouble>(alpha2));
+        _assert_(!xIsInf<IssmDouble>(alpha2));
 
         /*Assign output pointers:*/

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

@@ -43 +43 @@
                 void  GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2PISM(IssmDouble* palpha2,Gauss* gauss);
 
                 IssmDouble EffectivePressure(Gauss* gauss);