Changeset 25973

Timestamp:

02/03/21 19:54:22 (4 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added coulomb limited friction

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)
• m/classes/frictioncoulomb2.m (modified) (1 diff)

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

@@ -916 +916 @@
                         iomodel->FetchDataToInput(inputs,elements,"md.friction.f",FrictionfEnum);
                         break;
+                case 13:
+                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.coefficient",FrictionCEnum);
+                        iomodel->FetchDataToInput(inputs,elements,"md.friction.m",FrictionMEnum);
+                        if(FrictionCoupling==3){
+                                iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);}
+                        else if(FrictionCoupling==4){
+                                iomodel->FetchDataToInput(inputs,elements,"md.friction.effective_pressure",EffectivePressureEnum);
+
+                        }
+                        break;
                 default:
                         _error_("friction law "<< frictionlaw <<" not supported");
@@ -1027 +1038 @@
                 case 12:
                         parameters->AddObject(new IntParam(FrictionCouplingEnum,2));
+                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
+                        break;
+                case 13:
+                        parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
                         parameters->AddObject(iomodel->CopyConstantObject("md.friction.effective_pressure_limit",FrictionEffectivePressureLimitEnum));
                         break;

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

@@ -261 +261 @@
                 case 12:
                         GetAlpha2Tsai(palpha2,gauss);
+                        break;
+                case 13:
+                        GetAlpha2Coulomb2(palpha2,gauss);
                         break;
           default:
@@ -735 +738 @@
 
                 alpha2 = alpha2/ub;
+        }
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlpha2Coulomb2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient
+         *
+         *               C^2 |u_b|^(m-1) * (.5*N)
+         * alpha2= ___________________________________
+         *          (C^(2/m) |u_b| + (0.5*N)^(1/m) )^m
+         *
+         * */
+
+        /*diverse: */
+        IssmDouble  C,m,alpha2;
+
+        /*Recover parameters: */
+        element->GetInputValue(&C,gauss,FrictionCEnum);
+        element->GetInputValue(&m,gauss,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= (pow(C,2)*pow(ub,m-1.)*(0.5*N)) / pow(pow(C,2./m)*ub + pow(0.5*N,1./m),m);
         }
 

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

@@ -32 +32 @@
                 void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Coulomb(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2Coulomb2(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Josh(IssmDouble* palpha2,Gauss* gauss);

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

@@ -4 +4 @@
 %      frictioncoulomb=frictioncoulomb2();
 
-classdef frictioncoulomb
+classdef frictioncoulomb2
         properties (SetAccess=public) 
                 coefficient              = NaN;