Changeset 24064

Timestamp:

07/04/19 02:12:40 (6 years ago)

Author:

Mathieu Morlighem

Message:

CHG: do not grow channels on boundary

Location:

issm/trunk-jpl/src/c/classes/Loads

Files:

• Channel.cpp (modified) (4 diffs)
• Channel.h (modified) (1 diff)

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

@@ -54 +54 @@
         int e1 = iomodel->faces[4*index+2];
         int e2 = iomodel->faces[4*index+3];
+
+        if(e2==-1){
+                this->boundary = true;
+        }
+        else{
+                this->boundary = false;
+        }
 
         /*Set Element hook (4th column may be -1 for boundary edges)*/
@@ -545 +552 @@
                 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.);
+                /*Compute the effective conductivity Ks = k h^alpha |grad Phi|^{beta-2} (same for sheet)*/
+                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*/
@@ -585 +591 @@
         /*Initialize Element matrix and return if necessary*/
         Tria*  tria=(Tria*)element;
-        if(!tria->IsIceInElement()){
+        if(!tria->IsIceInElement() || this->boundary){
                 this->S = 0.;
                 return;
@@ -674 +680 @@
 
         /*Solve ODE*/
-        this->S = ODE1(alpha,beta,this->S,dt,0);
+        this->S = ODE1(alpha,beta,this->S,dt,2);
+
+        /*Make sure Area > 0*/
+        if(this->S<0.) this->S = 0.;
 
         /*Clean up and return*/

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

@@ -20 +20 @@
         private: 
                 IssmDouble S;
+                bool       boundary;
 
         public: