Changeset 24057

Timestamp:

07/01/19 06:45:47 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: use alpha and beta in thickness update and channel area update

Location:

issm/trunk-jpl/src/c

Files:

• analyses/HydrologyGlaDSAnalysis.cpp (modified) (2 diffs)
• classes/Loads/Channel.cpp (modified) (2 diffs)

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

@@ -425 +425 @@
 
         /*Intermediaries */
-        IssmDouble  Jdet,w,v,vx,vy,ub,h_old,N,h_r;
+        IssmDouble  Jdet,vx,vy,ub,h_old,N,h_r;
         IssmDouble  A,B,n;
+        IssmDouble  alpha,beta;
 
         /*Fetch number vertices for this element*/
@@ -462 +463 @@
                 ub = sqrt(vx*vx + vy*vy);
 
-                /*Compute cavity opening w*/
-                w  = 0.;
-                if(h_old<h_r) w = ub*(h_r-h_old)/l_r;
-
-                /*Compute closing rate*/
+                /*Get A from B and n*/
                 A=pow(B,-n);
-                v = 2./pow(n,n)*A*h_old*pow(fabs(N),n-1.)*N;
+
+                /*Define alpha and beta*/
+                if(h_old<h_r){
+                        alpha = -ub/l_r - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+                        beta  = ub*h_r/l_r;
+                }
+                else{
+                        alpha = - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+                        beta  = 0.;
+                }
 
                 /*Get new sheet thickness*/
-                h_new[iv] = h_old + dt*(w-v);
+                h_new[iv] = ODE1(alpha,beta,h_old,dt,0);
+
+                /*Make sure it is positive*/
                 if(h_new[iv]<1.e-12) h_new[iv] = 1.e-12;
         }

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

@@ -648 +648 @@
         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(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*/
-        IssmDouble qc = - Ks * dphids;
-
         /*d(phi - phi_m)/ds*/
         IssmDouble dPw = dphids - dphimds;
+
+        /*Approx. discharge in the sheet flowing folwing in the direction of the channel ofver a width lc*/
+        IssmDouble qc = - ks * pow(h,ALPHA_S) * pow(Ngrad,BETA_S-2.) * dphids;
 
         /*Compute f factor*/
@@ -664 +660 @@
         }
 
-        /*Compute total discharge*/
-        IssmDouble Q = -Kc*dphids;
-
-        /*Compute Pi and Xi*/
-        IssmDouble Pi = -C_W*c_t*rho_water*(Q+fFactor)*dPw;
-        IssmDouble Xi = fabs(Q*dphids) + fabs(lc * qc * dphids);
-
-        /*Compute closing rate*/
         A=pow(B,-n);
-        IssmDouble vc = 2./pow(n,n)*A*this->S*pow(fabs(phi_0 - phi),n-1.)*(phi_0 - phi);
-
-        /*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.;
+
+        IssmDouble C = C_W*c_t*rho_water;
+        IssmDouble Qprime = -kc * pow(Ngrad,BETA_C-2.)*dphids;
+        IssmDouble N = phi_0 - phi;
+
+        IssmDouble alpha = 1./(rho_ice*L)*(
+                                fabs(Qprime*pow(this->S,ALPHA_C-1.)*dphids)
+                                + C*Qprime*pow(this->S,ALPHA_C-1.)*dPw
+                                ) - 2./pow(n,n)*A*pow(fabs(N),n-1.)*N;
+
+        IssmDouble beta = 1./(rho_ice*L)*( fabs(lc*qc*dphids) + C*fFactor*dPw );
+
+        /*Solve ODE*/
+        this->S = ODE1(alpha,beta,this->S,dt,0);
 
         /*Clean up and return*/