Changeset 23961

Timestamp:

05/30/19 09:54:32 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: in non linear iteration, do not refer to N but use its linearization on phi

File:

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

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

@@ -169 +169 @@
         /*Intermediaries */
         IssmDouble  Jdet,dphi[3],h,k;
+        IssmDouble  A,B,n,phi_old,phi,phi_0,H,b,v1;
         IssmDouble* xyz_list = NULL;
 
@@ -189 +190 @@
         IssmDouble dt        = element->FindParam(TimesteppingTimeStepEnum);
         IssmDouble rho_water = element->FindParam(MaterialsRhoFreshwaterEnum);
+        IssmDouble rho_ice   = element->FindParam(MaterialsRhoIceEnum);
         IssmDouble g         = element->FindParam(ConstantsGEnum);
         IssmDouble e_v       = element->FindParam(HydrologyEnglacialVoidRatioEnum);
@@ -194 +196 @@
         Input* phi_input = element->GetInput(HydraulicPotentialEnum);      _assert_(phi_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* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
+        Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
 
         /* Start  looping on the number of gaussian points: */
@@ -205 +211 @@
 
                 phi_input->GetInputDerivativeValue(&dphi[0],xyz_list,gauss);
+                phi_input->GetInputValue(&phi,gauss);
                 h_input->GetInputValue(&h,gauss);
                 k_input->GetInputValue(&k,gauss);
+                B_input->GetInputValue(&B,gauss);
+                n_input->GetInputValue(&n,gauss);
+                b_input->GetInputValue(&b,gauss);
+                H_input->GetInputValue(&H,gauss);
 
                 /*Get norm of gradient of hydraulic potential and make sure it is >0*/
@@ -220 +231 @@
                                                         coeff*dbasis[0*numnodes+i]*dbasis[0*numnodes+j]
                                                         + coeff*dbasis[1*numnodes+i]*dbasis[1*numnodes+j]);
+                        }
+                }
+
+                /*Closing rate term, see Gagliardini and Werder 2018 eq. A2 (v = v1*phi_i + v2(phi_{i+1}))*/
+                phi_0 = rho_water*g*b + rho_ice*g*H;
+                A=pow(B,-n);
+                v1 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*( - n));
+                for(int i=0;i<numnodes;i++){
+                        for(int j=0;j<numnodes;j++){
+                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(-v1)*basis[i]*basis[j];
                         }
                 }
@@ -248 +269 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,w,v,vx,vy,ub,h,N,h_r;
+        IssmDouble  Jdet,w,v2,vx,vy,ub,h,h_r;
         IssmDouble  G,m,frictionheat,alpha2;
-        IssmDouble  A,B,n,phi_old;
+        IssmDouble  A,B,n,phi_old,phi,phi_0;
+        IssmDouble  H,b;
         IssmDouble* xyz_list = NULL;
 
@@ -272 +294 @@
         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* G_input      = element->GetInput(BasalforcingsGeothermalfluxEnum);_assert_(G_input);
         Input* B_input      = element->GetInput(MaterialsRheologyBEnum);         _assert_(B_input);
         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
         Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
+        Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
         /*Build friction element, needed later: */
@@ -297 +321 @@
                 B_input->GetInputValue(&B,gauss);
                 n_input->GetInputValue(&n,gauss);
-                N_input->GetInputValue(&N,gauss);
                 hr_input->GetInputValue(&h_r,gauss);
+                phi_input->GetInputValue(&phi,gauss);
+                b_input->GetInputValue(&b,gauss);
+                H_input->GetInputValue(&H,gauss);
 
                 /*Get basal velocity*/
@@ -315 +341 @@
 
                 /*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);
-                v = 2./pow(n,n)*A*h*pow(fabs(N),n-1.)*N;
-
-                for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v-m)*basis[i];
+                v2 = 2./pow(n,n)*A*h*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+= - Jdet*gauss->weight*(w-v2-m)*basis[i];
 
                 /*Transient term if dt>0*/
@@ -413 +441 @@
                 /*Get new sheet thickness*/
                 h_new[iv] = h_old + dt*(w-v);
+                if(h_new[iv]<1.e-12) h_new[iv] = 1.e-12;
         }