Changeset 26601

Timestamp:

11/10/21 18:30:41 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: better way if fixing friction singularity when v=0

File:

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

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

@@ -197 +197 @@
         IssmDouble vmag = VelMag(gauss);
 
-        /*Check to prevent dividing by zero if vmag==0*/
-        if(vmag==0. && (s-1.)<=0.) alpha_complement=0.;
-        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
+        if(s==1.){
+                /*This is to make AD happy and avoid 0^0*/
+                alpha_complement=pow(Neff,r);
+        }
+        else{
+                /*Check to prevent dividing by zero if vmag==0*/
+                if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
+                else alpha_complement=pow(Neff,r)*pow(vmag,(s-1.));
+        }
 
         /*Assign output pointers:*/
@@ -344 +350 @@
         IssmDouble vmag = VelMag(gauss);
 
-        /*Check to prevent dividing by zero if vmag==0*/
-        if(vmag==0. && (s-1.)<=0.){
-                alpha2=0.;
+        if(s==1.){
+                /*This is to make AD happy and avoid 0^0*/
+                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
         }
         else{
-                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+                /*Check to prevent dividing by zero if vmag==0*/
+                if(vmag==0. && (s-1.)<0.) alpha2=0.;
+                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
         }
 
@@ -541 +549 @@
 
         /*Check to prevent dividing by zero if vmag==0*/
-        if(vmag==0. && (s-1.)<=0.) alpha2=0.;
-        else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+        if(s==1.){
+                /*This is to make AD happy and avoid 0^0*/
+                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
+        }
+        else{
+                if(vmag==0. && (s-1.)<0.) alpha2=0.;
+                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+        }
 
         /*Assign output pointers:*/
@@ -586 +600 @@
         IssmDouble vmag = VelMag(gauss);
 
-        /*Check to prevent dividing by zero if vmag==0*/
-        if(vmag==0. && (s-1.)<=0.) alpha2=0.;
-        else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+        if(s==1.){
+                /*This is to make AD happy and avoid 0^0*/
+                alpha2=drag_coefficient*drag_coefficient*pow(Neff,r);
+        }
+        else{
+                /*Check to prevent dividing by zero if vmag==0*/
+                if(vmag==0. && (s-1.)<0.) alpha2=0.;
+                else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.));
+        }
 
         /*Assign output pointers:*/