Changeset 17348

Timestamp:

02/25/14 08:11:21 (11 years ago)

Author:

Mathieu Morlighem

Message:

NEW: added friction law for Henning:

File:

• issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp (modified) (6 diffs)

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

@@ -257 +257 @@
 
         /*Intermediaries */
-        IssmDouble ub,vb,slope2,drag,thickness,connectivity;
+        int        frictionlaw = 1;
+        IssmDouble ub,vb,slope2,drag,thickness,surface,connectivity;
         IssmDouble slope[2];
 
@@ -274 +275 @@
         Input* slopey_input    = element->GetInput(SurfaceSlopeYEnum);        _assert_(slopey_input);
         Input* thickness_input = element->GetInput(ThicknessEnum);            _assert_(thickness_input);
+        Input* surface_input   = element->GetInput(SurfaceEnum);              _assert_(surface_input);
         Input* drag_input      = element->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
 
@@ -283 +285 @@
 
                 thickness_input->GetInputValue(&thickness,gauss);
+                surface_input->GetInputValue(&surface,gauss);
                 drag_input->GetInputValue(&drag,gauss);
                 slopex_input->GetInputValue(&slope[0],gauss);
@@ -288 +291 @@
                 slope2=slope[0]*slope[0]+slope[1]*slope[1];
 
-                /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
-                ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
-                vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
-                ///*Ritz et al. 1996*/
-                //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
-                //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
-                ///*Rutt et al. 2009*/
-                //ub=-drag*rho_ice*gravity*thickness*slope[0];
-                //vb=-drag*rho_ice*gravity*thickness*slope[1];
+                switch(frictionlaw){
+                        case 1:
+                                /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
+                                ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+                                vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+                                break;
+                        case 2:
+                                /*Ritz et al. 1996*/
+                                ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+                                vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+                                break;
+                        case 3:
+                                /*Rutt et al. 2009*/
+                                ub=-drag*rho_ice*gravity*thickness*slope[0];
+                                vb=-drag*rho_ice*gravity*thickness*slope[1];
+                                break;
+                        case 4:
+                                /*Henning Akesson*/
+                                drag = -4e-15 * surface + 8.6e-12;
+                                ub=-drag*rho_ice*gravity*thickness*slope[0];
+                                vb=-drag*rho_ice*gravity*thickness*slope[1];
+                                break;
+                        default:
+                                _error_("Not supported yet");
+                }
 
                 pe->values[2*iv+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
@@ -309 +328 @@
 
         /*Intermediaries */
+        int         frictionlaw = 1;
         int         nodeup,nodedown,numsegments;
         IssmDouble  ub,vb,slope2,drag,surface,thickness,constant_part,z,Jdet;
@@ -374 +394 @@
                         drag_input->GetInputValue(&drag,gauss);
 
-                        /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
-                        ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
-                        vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
-                        ///*Ritz et al. 1996*/
-                        //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
-                        //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
-                        ///*Rutt et al. 2009*/
-                        //ub=-drag*rho_ice*gravity*thickness*slope[0];
-                        //vb=-drag*rho_ice*gravity*thickness*slope[1];
+                        switch(frictionlaw){
+                                case 1:
+                                        /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
+                                        ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
+                                        vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
+                                        break;
+                                case 2:
+                                        /*Ritz et al. 1996*/
+                                        ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
+                                        vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
+                                        break;
+                                case 3:
+                                        /*Rutt et al. 2009*/
+                                        ub=-drag*rho_ice*gravity*thickness*slope[0];
+                                        vb=-drag*rho_ice*gravity*thickness*slope[1];
+                                        break;
+                                case 4:
+                                        /*Henning Akesson*/
+                                        drag = -4e-15 * surface + 8.6e-12;
+                                        ub=-drag*rho_ice*gravity*thickness*slope[0];
+                                        vb=-drag*rho_ice*gravity*thickness*slope[1];
+                                        break;
+                                default:
+                                        _error_("Not supported yet");
+                        }
 
                         pe->values[2*nodedown+0]+=ub/connectivity[0];