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

Timestamp:

07/18/22 12:42:48 (3 years ago)

Author:

vverjans

Message:

CHG: allowing stochastic water pressure for frictionschoof

File:

• issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (modified) (2 diffs)

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

@@ -749 +749 @@
          *               C |u_b|^(m-1)
          * alpha2= __________________________
-         *          (1+(C/(Cmax N))^1/m |u_b| )^m
+         *          (1+(C/(Cmax Neff))^1/m |u_b| )^m
          *
          * */
@@ -764 +764 @@
         C = coeff*coeff;
 
-        /*Get effective pressure and velocity magnitude*/
-        IssmDouble N  = EffectivePressure(gauss);
+        /*Get effective pressure*/
+        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,FrictionSchoofWaterPressureEnum);
+                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 ub = VelMag(gauss);
 
         /*Compute alpha^2*/
-        if((ub<1e-10) ||(N==0.0)){
+        if((ub<1e-10) ||(Neff==0.0)){
                 alpha2 = 0.;
         }
         else{
-                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*N),1./m)*ub,m);
+                alpha2= (C*pow(ub,m-1.)) / pow(1.+  pow(C/(Cmax*Neff),1./m)*ub,m);
         }