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

[gurnis at geodynamics.org]

Author: gurnis
Date: 2006-10-08 17:14:42 -0700 (Sun, 08 Oct 2006)
New Revision: 4744

Modified:
   mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
   mc/3D/CitcomS/trunk/lib/Problem_related.c
   mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
Log:
For Viscosity_structures.c, added a new rheology, rheol=5, that is a variant of rheol=3. The difference lay in the way that plate boundaries are weakened. Must set mat_control=1 and enter materials into the VIP array. Removed debug statements from Problem_related.c in preparation of new release.

Modified: mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c	2006-10-08 13:56:34 UTC (rev 4743)
+++ mc/3D/CitcomS/trunk/lib/Full_read_input_from_files.c	2006-10-09 00:14:42 UTC (rev 4744)
@@ -73,11 +73,8 @@
       pos_age = 1;
     }
 
-    fprintf(stderr,"inside full... caps_per_proc=%d\n",E->sphere.caps_per_proc);
-
     for (m=1;m<=E->sphere.caps_per_proc;m++)  {
       cap = E->sphere.capid[m] - 1;  /* capid: 1-12 */
-      fprintf(stderr,"inside full... cap=%d\n",cap);
 
       switch (action) { /* set up files to open */
 

Modified: mc/3D/CitcomS/trunk/lib/Problem_related.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Problem_related.c	2006-10-08 13:56:34 UTC (rev 4743)
+++ mc/3D/CitcomS/trunk/lib/Problem_related.c	2006-10-09 00:14:42 UTC (rev 4744)
@@ -65,8 +65,6 @@
   const int dims=E->mesh.nsd,dofs=E->mesh.dof;
   const int ends=enodes[dims];
 
-  fprintf(stderr,"inside read_mat_from_files\n");
-
   elx=E->lmesh.elx;
   elz=E->lmesh.elz;
   ely=E->lmesh.ely;
@@ -98,7 +96,6 @@
 
     for(m=1;m<=E->sphere.caps_per_proc;m++) {
        cap = E->sphere.capid[m] - 1; /* capid: 1-12 */
-       fprintf(stderr,"Problem_related cap=%d\n",cap);
 
        sprintf(output_file,"%s%0.0f.%d",E->control.mat_file,newage1,cap);
        fprintf(E->fp,"%s %f %s\n","open mat file newage1",newage1,output_file);
@@ -126,9 +123,6 @@
     fclose(fp1);
     fclose(fp2);
 
-    fprintf(stderr,"lmesh.exs=%d\n",E->lmesh.exs);
-    fprintf(stderr,"lmesh.eys=%d\n",E->lmesh.eys);
-
     for (m=1;m<=E->sphere.caps_per_proc;m++)
       for (k=1;k<=ely;k++)
 	for (i=1;i<=elx;i++)   {

Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2006-10-08 13:56:34 UTC (rev 4743)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c	2006-10-09 00:14:42 UTC (rev 4744)
@@ -212,6 +212,7 @@
     int m,i,j,k,l,z,jj,kk,imark;
     float zero,e_6,one,eta0,Tave,depth,temp,tempa,temp1,TT[9];
     float zzz,zz[9];
+    float visc1, visc2, tempa_exp;
     const int vpts = vpoints[E->mesh.nsd];
     const int ends = enodes[E->mesh.nsd];
     const int nel = E->lmesh.nel;
@@ -314,8 +315,55 @@
       break;
 
 
+    case 5:
 
+      /* same as rheol 3, except alternative margin, VIP, formulation */
+      for(m=1;m<=E->sphere.caps_per_proc;m++)
+        for(i=1;i<=nel;i++)   {
+          l = E->mat[m][i];
+          tempa = E->viscosity.N0[l-1];
+          j = 0;
 
+          for(kk=1;kk<=ends;kk++) {
+            TT[kk] = E->T[m][E->ien[m][i].node[kk]];
+            zz[kk] = (1.-E->sx[m][3][E->ien[m][i].node[kk]]);
+          }
+
+          for(jj=1;jj<=vpts;jj++) {
+            temp=0.0;
+            zzz=0.0;
+            for(kk=1;kk<=ends;kk++)   {
+              TT[kk]=max(TT[kk],zero);
+              temp += min(TT[kk],one) * E->N.vpt[GNVINDEX(kk,jj)];
+              zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
+            }
+
+            if(E->control.mat_control==0)
+              EEta[m][ (i-1)*vpts + jj ] = 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(E->control.mat_control==1) {
+               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]) );
+               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 ]); */
+              }
+
+	    }
+        }
+      break;
+
+
+
+
     }
 
     return;