Changeset 24063

Timestamp:

07/04/19 02:11:46 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: do not use effective pressure in thickness update, recompute it instead, and use implicit scheme for sheet thickness update

File:

• issm/trunk-jpl/src/c/analyses/HydrologyGlaDSAnalysis.cpp (modified) (10 diffs)

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

@@ -33 +33 @@
                 CreateFaces(iomodel);
                 for(int i=0;i<iomodel->numberoffaces;i++){
+
                         /*Get left and right elements*/
                         int element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
@@ -303 +304 @@
 
         /*Intermediaries */
+        bool        meltflag;
         IssmDouble  Jdet,w,v2,vx,vy,ub,h,h_r;
         IssmDouble  G,m,frictionheat,alpha2;
@@ -317 +319 @@
 
         /*Retrieve all inputs and parameters*/
-        bool meltflag;
+        element->GetVerticesCoordinates(&xyz_list);
         element->FindParam(&meltflag,HydrologyMeltFlagEnum);
-        element->GetVerticesCoordinates(&xyz_list);
         IssmDouble L         = element->FindParam(MaterialsLatentheatEnum);
         IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
@@ -362 +363 @@
                 b_input->GetInputValue(&b,gauss);
                 H_input->GetInputValue(&H,gauss);
-                m_input->GetInputValue(&m,gauss);
 
                 /*Get basal velocity*/
@@ -375 +375 @@
                 frictionheat=alpha2*ub*ub;
 
-                /*Compute melt*/
+                /*Compute melt (if necessary)*/
                 if(!meltflag){
                         m = (G + frictionheat)/(rho_ice*L);
                 }
+                else{
+                        m_input->GetInputValue(&m,gauss);
+                }
 
                 /*Compute closing rate*/
-                /*See Gagliardini and Werder 2018 eq. A2 (v = v2(phi_i) + v1*phi_{i+1})*/
                 phi_0 = rho_water*g*b + rho_ice*g*H;
                 A=pow(B,-n);
@@ -433 +435 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r;
-        IssmDouble  A,B,n;
+        IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r,H,b;
+        IssmDouble  A,B,n,phi,phi_0;
         IssmDouble  alpha,beta;
 
@@ -444 +446 @@
 
         /*Retrieve all inputs and parameters*/
-        IssmDouble  dt  = element->FindParam(TimesteppingTimeStepEnum);
-        IssmDouble  l_r = element->FindParam(HydrologyCavitySpacingEnum);
+        IssmDouble  dt       = element->FindParam(TimesteppingTimeStepEnum);
+        IssmDouble  l_r      = element->FindParam(HydrologyCavitySpacingEnum);
+        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
+        IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble g         = element->FindParam(ConstantsGEnum);
         Input* hr_input = element->GetInput(HydrologyBumpHeightEnum);_assert_(hr_input);
         Input* vx_input = element->GetInput(VxEnum);_assert_(vx_input);
         Input* vy_input = element->GetInput(VyEnum);_assert_(vy_input);
-        Input*  N_input = element->GetInput(EffectivePressureEnum); _assert_(N_input);
-        Input*  h_input = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
+        Input* H_input  = element->GetInput(ThicknessEnum); _assert_(H_input);
+        Input* b_input  = element->GetInput(BedEnum); _assert_(b_input);
+        Input* h_input  = element->GetInput(HydrologySheetThicknessEnum);_assert_(h_input);
         Input* B_input  = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
         Input* n_input  = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
+        Input* phi_input = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -460 +467 @@
 
                 /*Get input values at gauss points*/
+                phi_input->GetInputValue(&phi,gauss);
                 vx_input->GetInputValue(&vx,gauss);
                 vy_input->GetInputValue(&vy,gauss);
@@ -465 +473 @@
                 B_input->GetInputValue(&B,gauss);
                 n_input->GetInputValue(&n,gauss);
-                N_input->GetInputValue(&N,gauss);
                 hr_input->GetInputValue(&h_r,gauss);
+                b_input->GetInputValue(&b,gauss);
+                H_input->GetInputValue(&H,gauss);
+
+                /*Get values for a few potentials*/
+                phi_0 = rho_water*g*b + rho_ice*g*H;
+                N = phi_0 - phi;
 
                 /*Get basal velocity*/
@@ -485 +498 @@
 
                 /*Get new sheet thickness*/
-                h_new[iv] = ODE1(alpha,beta,h_old,dt,0);
+                h_new[iv] = ODE1(alpha,beta,h_old,dt,1);
 
                 /*Make sure it is positive*/