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

gurnis at geodynamics.org gurnis at geodynamics.org
Mon Jan 29 12:32:28 PST 2007


Author: gurnis
Date: 2007-01-29 12:32:27 -0800 (Mon, 29 Jan 2007)
New Revision: 5921

Modified:
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
Small change for rheol=5 to use harmonic mean when element materials are read in. Also, this involved a minor conflict resolution with Eh s recent change to this file

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-01-27 20:28:12 UTC (rev 5920)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2007-01-29 20:32:27 UTC (rev 5921)
@@ -408,18 +408,22 @@
                                  - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
 
                     if(E->control.mat_control==1) {
-                        visc1 = E->VIP[m][i];
+                     /* visc1 = E->VIP[m][i];
                         visc2 = 2.0/(1./visc1 + 1.);
                         tempa_exp = tempa*
-                            exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
-                                 - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+	                exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+		           - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
                         visc1 = tempa*E->viscosity.max_value;
                         if(tempa_exp > visc1) tempa_exp=visc1;
                         EEta[m][ (i-1)*vpts + jj ] = visc2*tempa_exp;
-                        /* if(E->parallel.me == 0 && visc1 < 1.0e-03)
-                           fprintf(stderr,"%f  %f   %e  %e  %e\n",zzz,temp,visc1,visc2,
-                           EEta[m][ (i-1)*vpts + jj ]); */
-                    }
+                     */
+                       visc1 = E->VIP[m][i];
+                       visc2 = tempa*
+	               exp( E->viscosity.E[l-1]/(temp+E->viscosity.T[l-1])
+		          - E->viscosity.E[l-1]/(one +E->viscosity.T[l-1]) );
+                       if(visc1 <= 0.95) visc2=visc1;
+                       EEta[m][ (i-1)*vpts + jj ] = visc2;
+                      }
 
                 }
             }



More information about the cig-commits mailing list