Changeset 23977

Timestamp:

05/31/19 20:03:07 (6 years ago)

Author:

Mathieu Morlighem

Message:

BUG: fixed a few things...

File:

• issm/trunk-jpl/src/c/classes/Loads/Channel.cpp (modified) (10 diffs)

issm/trunk-jpl/src/c/classes/Loads/Channel.cpp ¶

@@ -381 +381 @@
         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
         Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -402 +401 @@
 
                 /*Get input values at gauss points*/
+                phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
+                b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
                 phi_input->GetInputValue(&phi,gauss);
-                phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
                 h_input->GetInputValue(&h,gauss);
                 ks_input->GetInputValue(&ks,gauss);
@@ -409 +409 @@
                 n_input->GetInputValue(&n,gauss);
                 b_input->GetInputValue(&b,gauss);
-                b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
                 H_input->GetInputValue(&H,gauss);
 
@@ -484 +483 @@
         /*Intermediaries */
         IssmDouble  Jdet,v2,Afactor,Bfactor,fFactor;
-        IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds;
-        IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks;
+        IssmDouble  A,B,n,phi_old,phi,phi_0,dphimds,dphi[2];
+        IssmDouble  H,h,b,db[2],dphids,qc,dPw,ks,Ngrad;
         IssmDouble  xyz_list[NUMVERTICES][3];
         IssmDouble  xyz_list_tria[3][3];
@@ -512 +511 @@
         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
         Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -531 +529 @@
                 /*Get input values at gauss points*/
                 b_input->GetInputDerivativeValue(&db[0],&xyz_list_tria[0][0],gauss);
+                phi_input->GetInputDerivativeValue(&dphi[0],&xyz_list_tria[0][0],gauss);
                 h_input->GetInputValue(&h,gauss);
                 ks_input->GetInputValue(&ks,gauss);
@@ -540 +539 @@
 
                 /*Get values for a few potentials*/
+                phi_0   = rho_water*g*b + rho_ice*g*H;
+                dphids  = dphi[0]*tx + dphi[1]*ty;
                 dphimds = rho_water*g*(db[0]*tx + db[1]*ty);
+                Ngrad   = fabs(dphids);
+                if(Ngrad<1.e-12) Ngrad = 1.e-12;
 
                 /*Compute the effective conductivity Kc = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
-                IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(fabs(dphids),BETA_C-2.);
-                IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(fabs(dphids),BETA_S-2.);
+                IssmDouble Kc = kc * pow(this->S,ALPHA_C) * pow(Ngrad,BETA_C-2.);
+                IssmDouble Ks = ks * pow(h      ,ALPHA_S) * pow(Ngrad,BETA_S-2.);
 
                 /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
@@ -564 +567 @@
                 /*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);
                 v2 = 2./pow(n,n)*A*this->S*(pow(fabs(phi_0 - phi),n-1.)*(phi_0 +(n-1.)*phi));
@@ -616 +618 @@
         Input* n_input      = element->GetInput(MaterialsRheologyNEnum);         _assert_(n_input);
         Input* ks_input     = element->GetInput(HydrologySheetConductivityEnum); _assert_(ks_input);
-        Input* phiold_input = element->GetInput(HydraulicPotentialOldEnum);      _assert_(phiold_input);
         Input* phi_input    = element->GetInput(HydraulicPotentialEnum);         _assert_(phi_input);
 
@@ -676 +677 @@
         /*Compute new S based on Forward Euler (explicit)*/
         this->S = this->S + dt*( (Xi - Pi)/(rho_ice*L) - vc);
-        if(this->S<0); this->S = 0.;
+        if(this->S<0.); this->S = 0.;
 
         /*Clean up and return*/