[cig-commits] r7932 - mc/3D/CitcomS/trunk/lib

[tan2 at geodynamics.org]

Author: tan2
Date: 2007-09-05 15:20:46 -0700 (Wed, 05 Sep 2007)
New Revision: 7932

Modified:
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
Set the min of 2nd strain rate invariant to 1e-16 to prevent infinite viscosity when SDEPV=on

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-09-05 22:20:08 UTC (rev 7931)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-09-05 22:20:46 UTC (rev 7932)
@@ -83,6 +83,7 @@
     E->viscosity.sdepv_misfit = 1.0;
     input_boolean("SDEPV",&(E->viscosity.SDEPV),"off",m);
     if (E->viscosity.SDEPV) {
+      E->viscosity.pdepv_visited = 0;
       input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
     }
 
@@ -526,8 +527,14 @@
     two = 2.0;
 
     for(m=1;m<=E->sphere.caps_per_proc;m++)  {
-        strain_rate_2_inv(E,m,eedot,1);	/* should there be a check if velocities have been computed here? */
 
+        /* get second invariant for all elements */
+        strain_rate_2_inv(E,m,eedot,1);
+
+        /* eedot cannot be too small, or the viscosity will go to inf */
+	for(e=1;e<=nel;e++)
+            eedot[e] = max(eedot[e], 1.0e-16);
+
         for(e=1;e<=nel;e++)   {
             exponent1= one/E->viscosity.sdepv_expt[E->mat[m][e]-1];
             scale=pow(eedot[e],exponent1-one);