Changeset 23114

Timestamp:

08/17/18 13:45:39 (7 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on PISM friction law

Location:

issm/trunk-jpl/src/c

Files:

• analyses/StressbalanceAnalysis.cpp (modified) (2 diffs)
• classes/Loads/Friction.cpp (modified) (1 diff)

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

@@ -845 +845 @@
                         InputUpdateFromConstantx(elements,1.,FrictionQEnum);
                         break;
+                case 10:
+                        iomodel->FetchDataToInput(elements,"md.friction.till_friction_angle",FrictionTillFrictionAngleEnum);
+                        iomodel->FetchDataToInput(elements,"md.friction.sediment_compressibility_coefficient",FrictionSedimentCompressibilityCoefficientEnum);
+                        iomodel->FetchDataToInput(elements,"md.friction.watercolumn_max",FrictionWatercolumnMaxEnum);
+                        break;
                 default:
-                        _error_("not supported");
+                        _error_("friction law "<< frictionlaw <<" not supported");
         }
 
@@ -910 +915 @@
         if(frictionlaw==5) parameters->AddObject(iomodel->CopyConstantObject("md.friction.f",FrictionFEnum));
         if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
+        if(frictionlaw==10){
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.pseudoplasticity_exponent",FrictionPseudoplasticityExponentEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.threshold_speed",FrictionThresholdSpeedEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.delta",FrictionDeltaEnum));
+                parameters->AddObject(iomodel->CopyConstantObject("md.friction.void_ratio",FrictionVoidRatioEnum));
+        }
 
 }/*}}}*/

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

@@ -681 +681 @@
          */
 
-        /*Compute effective pressure first */
-        IssmDouble  N,delta,W,Wmax,e0,Cc,P0;
+        /*compute ice pressure P0*/
+        IssmDouble thickness,base,P0;
+        element->GetInputValue(&thickness, gauss,ThicknessEnum);
+        element->GetInputValue(&base, gauss,BaseEnum);
+        //element->GetInputValue(&sealevel, gauss,SealevelEnum);
+        IssmDouble rho_ice   = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble gravity   = element->GetMaterialParameter(ConstantsGEnum);
+        P0 = gravity*rho_ice*thickness;
+
+        /*Compute effective pressure*/
+        IssmDouble  N,delta,W,Wmax,e0,Cc;
         element->parameters->FindParam(&delta,FrictionDeltaEnum);
         element->parameters->FindParam(&e0,FrictionVoidRatioEnum);
-        //element->parameters->FindParam(&P0,XXEnum);
-        _error_("don't know what P0 is");
         element->GetInputValue(&Cc,gauss,FrictionSedimentCompressibilityCoefficientEnum);
         element->GetInputValue(&W,gauss,WatercolumnEnum);
-        Wmax = 0.; _error_("Don't know what to do with Wmax yet");
+        element->GetInputValue(&Wmax,gauss,FrictionWatercolumnMaxEnum);
         N = delta*P0*pow(10.,(e0/Cc)*(1.-W/Wmax));