Changeset 18770

Timestamp:

11/10/14 14:32:06 (10 years ago)

Author:

seroussi

Message:

NEW: starting to add new friction law with water layer

Location:

issm/trunk-jpl/src

Files:

• c/classes/Loads/Friction.cpp (modified) (3 diffs)
• c/classes/Loads/Friction.h (modified) (1 diff)
• m/classes/frictionwaterlayer.m (added)

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

@@ -54 +54 @@
                 case 4:
                         GetAlpha2Temp(palpha2,gauss);
+                        break;
+                case 5:
+                        GetAlpha2WaterLayer(palpha2,gauss);
                         break;
           default:
@@ -250 +253 @@
         *palpha2=alpha2;
 }/*}}}*/
-void Friction::GetAlphaComplement(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");
+void Friction::GetAlpha2WaterLayer(IssmDouble* palpha2, Gauss* gauss){/*{{{*/
+
+        /*This routine calculates the basal friction coefficient 
+          alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
 
         /*diverse: */
         IssmDouble  r,s;
-        IssmDouble  vx,vy,vz,vmag;
-        IssmDouble  drag_p,drag_q;
+        IssmDouble  drag_p, drag_q;
         IssmDouble  Neff;
+        IssmDouble  thickness,bed;
+        IssmDouble  vx,vy,vz,vmag;
         IssmDouble  drag_coefficient;
-        IssmDouble  bed,thickness;
-        IssmDouble  alpha_complement;
+        IssmDouble  alpha2;
 
         /*Recover parameters: */
@@ -285 +285 @@
         if(Neff<0)Neff=0;
 
+        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");
+        }
+
+        /*Check to prevent dividing by zero if vmag==0*/
+        if(vmag==0. && (s-1.)<0.) alpha2=0.;
+        else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+        _assert_(!xIsNan<IssmDouble>(alpha2));
+
+        /*Assign output pointers:*/
+        *palpha2=alpha2;
+}/*}}}*/
+void Friction::GetAlphaComplement(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: */
+        IssmDouble  r,s;
+        IssmDouble  vx,vy,vz,vmag;
+        IssmDouble  drag_p,drag_q;
+        IssmDouble  Neff;
+        IssmDouble  drag_coefficient;
+        IssmDouble  bed,thickness;
+        IssmDouble  alpha_complement;
+
+        /*Recover parameters: */
+        element->GetInputValue(&drag_p,FrictionPEnum);
+        element->GetInputValue(&drag_q,FrictionQEnum);
+        element->GetInputValue(&thickness, gauss,ThicknessEnum);
+        element->GetInputValue(&bed, gauss,BaseEnum);
+        element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum);
+        IssmDouble rho_water   = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+        IssmDouble rho_ice     = element->GetMaterialParameter(MaterialsRhoIceEnum);
+        IssmDouble gravity     = element->GetMaterialParameter(ConstantsGEnum);
+
+        //compute r and q coefficients: */
+        r=drag_q/drag_p;
+        s=1./drag_p;
+
+        //From bed and thickness, compute effective pressure when drag is viscous:
+        Neff=gravity*(rho_ice*thickness+rho_water*bed);
+        if(Neff<0)Neff=0;
+
         //We need the velocity magnitude to evaluate the basal stress:
         switch(dim){

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

@@ -34 +34 @@
                 void  GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss);
+                void  GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss);
                 void  GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss);
 };