Changeset 18992

Timestamp:

01/06/15 16:42:52 (10 years ago)

Author:

bdef

Message:

NEW: adding the complement to invert for the sliding parameter in Hydro friction law

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

@@ -42 +42 @@
 void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
 
+        switch(this->law){
+                case 1:
+                        GetAlphaViscousComplement(palpha_complement,gauss);
+                        break;
+                case 3:
+                        GetAlphaHydroComplement(palpha_complement,gauss);
+                        break;
+          default:
+                        _error_("not supported");
+        }
+
+}/*}}}*/
+
+void Friction::GetAlphaViscousComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
+
         /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. 
          * FrictionGetAlphaComplement is used in control methods on drag, and it computes: 
          * alpha_complement= Neff ^r * vel ^s*/
-
-        if(this->law!=1)_error_("not supported");
 
         /*diverse: */
@@ -104 +117 @@
 }
 /*}}}*/
+void Friction::GetAlphaHydroComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/
+
+        /*diverse: */
+        IssmDouble  q_exp;
+        IssmDouble  C_param;
+        IssmDouble  As;
+        IssmDouble  Neff;
+        IssmDouble  n;
+        IssmDouble  alpha;
+        IssmDouble  Chi;
+        IssmDouble  vx,vy,vz,vmag;
+        IssmDouble  alpha_complement;
+
+        /*Recover parameters: */
+        element->GetInputValue(&q_exp,FrictionQEnum);
+        element->GetInputValue(&C_param,FrictionCEnum);
+        element->GetInputValue(&As,FrictionAsEnum);
+
+        element->GetInputValue(&Neff,gauss,FrictionEffectivePressureEnum);
+        element->GetInputValue(&n,gauss,MaterialsRheologyNEnum);
+
+        if(Neff<0)Neff=0;
+
+        //We need the velocity magnitude to evaluate the basal stress:
+        switch(dim){
+                case 1:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        vmag=sqrt(vx*vx);
+                        break;
+                case 2:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        element->GetInputValue(&vy,gauss,VyEnum);
+                        vmag=sqrt(vx*vx+vy*vy);
+                        break;
+                case 3:
+                        element->GetInputValue(&vx,gauss,VxEnum);
+                        element->GetInputValue(&vy,gauss,VyEnum);
+                        element->GetInputValue(&vz,gauss,VzEnum);
+                        vmag=sqrt(vx*vx+vy*vy+vz*vz);
+                        break;
+                default:
+                        _error_("not supported");
+        }
+        if (q_exp==1){
+                alpha=1;
+        }
+        else{
+                alpha=(pow(q_exp-1,q_exp-1))/pow(q_exp,q_exp);
+        }
+        Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(vmag==0.) alpha_complement=0.;
+        else if(Neff==0.) alpha_complement=0.;
+        else alpha_complement=-(C_param*Neff/(vmag*n)) * pow((Chi/(alpha*pow(Chi,q_exp)+1)),((1-n)/n)) *(Chi/(As*(alpha*pow(Chi,q_exp)+1)))-(alpha*q_exp*pow(Chi,q_exp+1)/(As*(alpha*pow(Chi,q_exp)+1)));
+        _assert_(!xIsNan<IssmDouble>(alpha_complement));
+        /*Assign output pointers:*/
+        *palpha_complement=alpha_complement;
+}
+/*}}}*/
 void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
 
@@ -191 +264 @@
         }
         Chi=vmag/(pow(C_param,n)*pow(Neff,n)*As);
-
 
         /*Check to prevent dividing by zero if vmag==0*/

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

@@ -30 +30 @@
                 void  Echo(void);
                 void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
+                void  GetAlphaViscousComplement(IssmDouble* alpha_complement,Gauss* gauss);
+                void  GetAlphaHydroComplement(IssmDouble* alpha_complement,Gauss* gauss);
                 void  GetAlpha2(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);