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

[becker at geodynamics.org]

Author: becker
Date: 2011-12-20 10:04:26 -0800 (Tue, 20 Dec 2011)
New Revision: 19305

Modified:
   mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
   mc/3D/CitcomS/trunk/lib/Instructions.c
   mc/3D/CitcomS/trunk/lib/Topo_gravity.c
   mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Merged with recent check in, minor changes. 



Modified: mc/3D/CitcomS/trunk/lib/Ggrd_handling.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Ggrd_handling.c	2011-12-20 01:44:07 UTC (rev 19304)
+++ mc/3D/CitcomS/trunk/lib/Ggrd_handling.c	2011-12-20 18:04:26 UTC (rev 19305)
@@ -420,6 +420,9 @@
   /* no interpolation */
   use_nearneighbor = TRUE;
 #endif
+  if(E->control.ggrd_mat_is_code)
+    use_nearneighbor = TRUE;
+    
   /*
      if we have not initialized the time history structure, do it now
   */
@@ -510,7 +513,9 @@
       i1 = 0;
     }
     /*
+      
        loop through all elements and assign
+       
     */
     for (m=1;m <= E->sphere.caps_per_proc;m++) {
       for (j=1;j <= elz;j++)  {	/* this assumes a regular grid sorted as in (1)!!! */
@@ -591,17 +596,56 @@
 	    }
 	  }
 	}else{
-	  /* outside the lithosphere */
+	  /* outside the lithosphere/layers under consideration */
+	  if(E->control.ggrd_mat_is_code){
+	    for (k=1;k <= ely;k++){
+	      for (i=1;i <= elx;i++)   {
+		el = j + (i-1) * elz + (k-1)*elxlz;
+		E->VIP[m][el] = 0; /* zero code --> unity scale */
+	      }
+	    }
+	  }else{
+	    for (k=1;k <= ely;k++){
+	      for (i=1;i <= elx;i++)   {
+		el = j + (i-1) * elz + (k-1)*elxlz;
+		E->VIP[m][el] = 1.0;
+	      }
+	    }
+	  }
+	}
+      }	/* end elz loop */
+    } /* end m loop */
+    if(E->control.ggrd_mat_is_code){
+      /* 
+	 code assignment 
+      */
+    
+      if(E->parallel.me==0){
+	for(i=0;i < E->control.ggrd_mat_is_code;i++)
+	  fprintf(stderr,"ggrd_read_mat_from_file: code %03i: viscosity: %12.3e\n",
+		  i+1,E->control.ggrd_mat_code_viscosities[i]);
+      }
+      
+      /* the grids were actually code flags from 1...cmax */
+      for (m=1;m <= E->sphere.caps_per_proc;m++) {
+	for (j=1;j <= elz;j++) {
 	  for (k=1;k <= ely;k++){
 	    for (i=1;i <= elx;i++)   {
 	      el = j + (i-1) * elz + (k-1)*elxlz;
-	      /* no scaling else */
-	      E->VIP[m][el] = 1.0;
+	      if((int)E->VIP[m][el] < 1){ /* background */
+		E->VIP[m][el] = 1.0;
+	      }else{
+		if((((int)E->VIP[m][el]) > E->control.ggrd_mat_is_code)||(((int)E->VIP[m][el]) < 1)){
+		  fprintf(stderr,"%i\n",(int)E->VIP[m][el]);
+		  myerror(E,"ggrd_mat_code_viscosities: input code out of bounds");
+		}
+		E->VIP[m][el] = E->control.ggrd_mat_code_viscosities[(int)(E->VIP[m][el]-1)];
+	      }
 	    }
 	  }
 	}
-      }	/* end elz loop */
-    } /* end m loop */
+      }
+    }
   } /* end assignment loop */
   if((!timedep) && (!E->control.ggrd.mat_control_init)){
     /* forget the grid */

Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c	2011-12-20 01:44:07 UTC (rev 19304)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c	2011-12-20 18:04:26 UTC (rev 19305)
@@ -278,7 +278,7 @@
   void set_mg_defaults();
   float tmp;
   double ell_tmp;
-  int m=E->parallel.me;
+  int m=E->parallel.me,i;
   double levmax;
   int tmp_int_in;
 
@@ -542,6 +542,21 @@
   /* if < 0, will assign only to layer == -ggrd_mat_control */
   input_int("ggrd_mat_control",&(E->control.ggrd.mat_control),"0",m); 
   input_boolean("ggrd_mat_limit_prefactor",&(E->control.ggrd_mat_limit_prefactor),"on",m); /* limit prefactor to with 1e+/-5 */
+  input_int("ggrd_mat_is_code",&(E->control.ggrd_mat_is_code),"0",m); /* the viscosity grids are
+									   actually codes for
+									   different types of
+									   rheologies, from 1 .... cmax
+									   
+									   
+									*/
+  if(E->control.ggrd_mat_is_code){
+    /* we need code assignments */
+    E->control.ggrd_mat_code_viscosities = (float *)malloc(sizeof(float)*E->control.ggrd_mat_is_code);
+    for(i=0;i < E->control.ggrd_mat_is_code;i++)
+      E->control.ggrd_mat_code_viscosities[i] = 1;
+    input_float_vector("ggrd_mat_code_viscosities",
+		       E->control.ggrd_mat_is_code,(E->control.ggrd_mat_code_viscosities),m);
+  }
   input_string("ggrd_mat_file",E->control.ggrd.mat_file,"",m); /* file to read prefactors from */
   input_string("ggrd_mat_depth_file",
 	       E->control.ggrd_mat_depth_file,"_i_do_not_exist_",m); 

Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c	2011-12-20 01:44:07 UTC (rev 19304)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c	2011-12-20 18:04:26 UTC (rev 19305)
@@ -304,7 +304,7 @@
 	if ((e-1)%E->lmesh.elz==0) {
 	  construct_c3x3matrix_el(E,e,&E->element_Cc,&E->element_Ccx,lev,m,0);
 	}
-	/* get b at velocity Gauss points */
+	/* get B at velocity Gauss points */
 	get_ba(&(E->N), GNx, &E->element_Cc, &E->element_Ccx,rtf, dims, ba);
 	/* regular stress tensor */
 	for(p=1;p <= 6;p++)

Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h	2011-12-20 01:44:07 UTC (rev 19304)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h	2011-12-20 18:04:26 UTC (rev 19305)
@@ -534,7 +534,8 @@
   float ggrd_vtop_omega[4];
   char ggrd_mat_depth_file[1000];
   ggrd_boolean ggrd_mat_is_3d;
-  int ggrd_mat_limit_prefactor;
+  int ggrd_mat_limit_prefactor,ggrd_mat_is_code;
+  float *ggrd_mat_code_viscosities;
   float  ggrd_lower_depth_km,ggrd_lower_scale,ggrd_lower_offset;
 #endif
     double accuracy,inner_accuracy_scale;