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

tan2 at geodynamics.org tan2 at geodynamics.org
Wed Sep 12 12:24:30 PDT 2007


Author: tan2
Date: 2007-09-12 12:24:29 -0700 (Wed, 12 Sep 2007)
New Revision: 7957

Modified:
   mc/3D/CitcomS/trunk/lib/Instructions.c
   mc/3D/CitcomS/trunk/lib/Output.c
   mc/3D/CitcomS/trunk/lib/Output_gzdir.c
   mc/3D/CitcomS/trunk/lib/Tracer_setup.c
   mc/3D/CitcomS/trunk/lib/global_defs.h
Log:
Added heating file to optional output

Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c	2007-09-12 19:24:15 UTC (rev 7956)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c	2007-09-12 19:24:29 UTC (rev 7957)
@@ -1145,6 +1145,7 @@
     E->output.tracer = 0;
     E->output.comp_el = 0;
     E->output.comp_nd = 0;
+    E->output.heating = 0;
 
     while(1) {
         /* get next field */
@@ -1190,6 +1191,8 @@
             E->output.comp_el = 1;
         else if(strcmp(prev, "comp_nd")==0)
             E->output.comp_nd = 1;
+        else if(strcmp(prev, "heating")==0)
+            E->output.heating = 1;
         else
             if(E->parallel.me == 0)
                 fprintf(stderr, "Warning: unknown field for output_optional: %s\n", prev);
@@ -1325,7 +1328,7 @@
     open_log(E);
     open_time(E);
     open_info(E);
-   
+
     if (strcmp(E->output.format, "ascii") == 0) {
         E->problem_output = output;
     }

Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c	2007-09-12 19:24:15 UTC (rev 7956)
+++ mc/3D/CitcomS/trunk/lib/Output.c	2007-09-12 19:24:29 UTC (rev 7957)
@@ -98,32 +98,33 @@
 
   output_surf_botm(E, cycles);
 
-  if(E->control.disptn_number != 0)
-    output_heating(E, cycles);
+  /* optiotnal output below */
 
-  /* optiotnal output below */
   /* compute and output geoid (in spherical harmonics coeff) */
-  if (E->output.geoid == 1)
+  if (E->output.geoid)
       output_geoid(E, cycles);
 
-  if (E->output.stress == 1)
+  if (E->output.stress)
     output_stress(E, cycles);
 
-  if (E->output.pressure == 1)
+  if (E->output.pressure)
     output_pressure(E, cycles);
 
-  if (E->output.horiz_avg == 1)
+  if (E->output.horiz_avg)
       output_horiz_avg(E, cycles);
 
-  if(E->output.tracer == 1 && E->control.tracer == 1)
+  if(E->output.tracer && E->control.tracer)
       output_tracer(E, cycles);
 
-  if (E->output.comp_nd == 1 && E->composition.on)
+  if (E->output.comp_nd && E->composition.on)
       output_comp_nd(E, cycles);
 
-  if (E->output.comp_el == 1 && E->composition.on)
-          output_comp_el(E, cycles);
+  if (E->output.comp_el && E->composition.on)
+      output_comp_el(E, cycles);
 
+  if(E->output.heating && E->control.disptn_number != 0)
+      output_heating(E, cycles);
+
   return;
 }
 

Modified: mc/3D/CitcomS/trunk/lib/Output_gzdir.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_gzdir.c	2007-09-12 19:24:15 UTC (rev 7956)
+++ mc/3D/CitcomS/trunk/lib/Output_gzdir.c	2007-09-12 19:24:29 UTC (rev 7957)
@@ -95,6 +95,7 @@
 void gzdir_output_horiz_avg(struct All_variables *, int);
 void gzdir_output_tracer(struct All_variables *, int);
 void gzdir_output_pressure(struct All_variables *, int);
+void gzdir_output_heating(struct All_variables *, int);
 
 
 void sub_netr(float, float, float, float *, float *, double *);
@@ -154,16 +155,16 @@
 
   /* optiotnal output below */
   /* compute and output geoid (in spherical harmonics coeff) */
-  if (E->output.geoid == 1)
+  if (E->output.geoid)
       gzdir_output_geoid(E, cycles);
 
-  if (E->output.stress == 1)
+  if (E->output.stress)
     gzdir_output_stress(E, cycles);
 
-  if (E->output.pressure == 1)
+  if (E->output.pressure)
     gzdir_output_pressure(E, cycles);
 
-  if (E->output.horiz_avg == 1)
+  if (E->output.horiz_avg)
       gzdir_output_horiz_avg(E, cycles);
 
   if(E->control.tracer){
@@ -176,8 +177,11 @@
       gzdir_output_comp_nd(E, cycles);
 
   if (E->output.comp_el && E->composition.on)
-          gzdir_output_comp_el(E, cycles);
+      gzdir_output_comp_el(E, cycles);
 
+  if(E->output.heating && E->control.disptn_number != 0)
+      gzdir_output_heating(E, cycles);
+
   return;
 }
 
@@ -1083,6 +1087,31 @@
     return;
 }
 
+
+void output_heating(struct All_variables *E, int cycles)
+{
+    int j, e;
+    char output_file[255];
+    FILE *fp1;
+
+    snprintf(output_file,255,"%s/%d/heating.%d.%d.gz", E->control.data_dir,
+	    cycles,E->parallel.me, cycles);
+    fp1 = gzdir_output_open(output_file,"w");
+
+    gzfprintf(fp1,"%.5e\n",E->monitor.elapsed_time);
+
+    for(j=1;j<=E->sphere.caps_per_proc;j++) {
+        gzfprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
+        for(e=1; e<=E->lmesh.nel; e++)
+            gzfprintf(fp1, "%.4e %.4e %.4e\n", E->heating_adi[j][e],
+                      E->heating_visc[j][e], E->heating_latent[j][e]);
+    }
+    gzclose(fp1);
+
+    return;
+}
+
+
 /*
 
 restart facility for zipped/VTK style , will init temperature

Modified: mc/3D/CitcomS/trunk/lib/Tracer_setup.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Tracer_setup.c	2007-09-12 19:24:15 UTC (rev 7956)
+++ mc/3D/CitcomS/trunk/lib/Tracer_setup.c	2007-09-12 19:24:29 UTC (rev 7957)
@@ -742,6 +742,7 @@
     double theta,phi,rad;
     double xmin,xmax,ymin,ymax,zmin,zmax;
     double random1,random2,random3;
+    double u, v, s, r;
 
 
     allocate_tracer_arrays(E,j,tracers_cap);
@@ -769,30 +770,25 @@
 
     while (E->trace.ntracers[j]<tracers_cap) {
 
-        number_of_tries++;
-        max_tries=100*tracers_cap;
+        do {
+            /* pick two uniformly distributed random numbers in [-1;1] */
+            u = -1 + drand48()*2;
+            v = -1 + drand48()*2;
+            s = u*u + v*v;
+            /* length has to be <= 1 */
+        } while(s > 1);
+        r = 2.0 * sqrt(1.0-s);
 
-        if (number_of_tries>max_tries) {
-            fprintf(E->trace.fpt,"Error(make_tracer_array)-too many tries?\n");
-            fprintf(E->trace.fpt,"%d %d %d\n",max_tries,number_of_tries,RAND_MAX);
-            fflush(E->trace.fpt);
-            exit(10);
-        }
+        /* cartesian coordinates */
+        x = u * r;
+        y = v * r;
+        z = 2.0*s -1 ;
 
-#if 1
-        random1=drand48();
-        random2=drand48();
-        random3=drand48();
-#else
-        random1=(1.0*rand())/(1.0*RAND_MAX);
-        random2=(1.0*rand())/(1.0*RAND_MAX);
-        random3=(1.0*rand())/(1.0*RAND_MAX);
-#endif
+        /* skip if outside the bounding box */
+        if ((x < xmin) || (x > xmax)) continue;
+        if ((y < ymin) || (y > ymax)) continue;
+        if ((z < zmin) || (z > zmax)) continue;
 
-        x=xmin+random1*(xmax-xmin);
-        y=ymin+random2*(ymax-ymin);
-        z=zmin+random3*(zmax-zmin);
-
         /* first check if within shell */
 
         cart_to_sphere(E,x,y,z,&theta,&phi,&rad);

Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h	2007-09-12 19:24:15 UTC (rev 7956)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h	2007-09-12 19:24:29 UTC (rev 7957)
@@ -390,9 +390,9 @@
 
 struct CONTROL {
     int PID;
-  
+
     char output_written_external_command[500];   /* a unix command to run when output files have been created */
-  
+
   /* this clashed with a GMT definition of ORTHO TWB */
   int CITCOM_ORTHO,CITCOM_ORTHOZ;   /* indicates levels of mesh symmetry */
 
@@ -621,6 +621,7 @@
     int tracer;       /* whether to output tracer coordinate */
     int comp_el;      /* whether to output composition at elements */
     int comp_nd;      /* whether to output composition at nodes */
+    int heating;      /* whether to output heating terms at elements */
 
 
   /* flags used by GZDIR */



More information about the cig-commits mailing list