[cig-commits] r9191 - cs/avm/trunk

[tan2 at geodynamics.org]

Author: tan2
Date: 2008-01-29 13:49:29 -0800 (Tue, 29 Jan 2008)
New Revision: 9191

Modified:
   cs/avm/trunk/avm-partition.c
Log:
Fixed bugs and non-dimensionalize [rho, ks, mu]

Modified: cs/avm/trunk/avm-partition.c
===================================================================
--- cs/avm/trunk/avm-partition.c	2008-01-29 21:48:57 UTC (rev 9190)
+++ cs/avm/trunk/avm-partition.c	2008-01-29 21:49:29 UTC (rev 9191)
@@ -30,9 +30,29 @@
 #include <stdio.h>
 #include <stdlib.h>
 
-#define IREGION_CRUST_MANTLE 1
+/* at least five columns in the input files */
+/* which are [rank, order, temperature, pressure, density] */
+#define NUM_REQUIRED_COLUMNS 2
+
+/* max. number of compositions components */
+/* increase it if more components are needed */
 #define MAX_NUM_COMPOSITIONS 10
 
+/* reference density (SI unit) */
+/* corresponding to const.density in CitcomS */
+#define REFERENCE_DENSITY 3340.0
+
+/* gravititional constant from 2006 CODATA (SI unit) */
+#define GRAV_CONST 6.67428e-11
+
+/* a constant for specfem mesh data */
+#define IREGION_CRUST_MANTLE 1
+
+
+/****************************************************************************/
+/****************************************************************************/
+
+
 /*
  * This function takes the temperatur (in Kelvin), pressure (in Pascal), and
  * composition (a vector of concentration) at a point, and returns the
@@ -47,7 +67,23 @@
     return;
 }
 
+/*
+ * Take dimensional density and elastic moduli (SI units) and return
+ * non-dimensionalized results.
+ * See SPECFEM3D manaul for the dimensional scheme.
+ */
+void non_dimensionalize(double *rho, double *ks, double *mu)
+{
+    const double inv_ref_rho = 1.0 / REFERENCE_DENSITY;
+    const double inv_ref_moduli = 1.0;
 
+    *rho *= inv_ref_rho;
+    *ks *= inv_ref_moduli;
+    *mu *= inv_ref_moduli;
+    return;
+}
+
+
 /*
  * Scan input str to get a double vector *values. The vector length is from
  * input len. The input str contains white-space seperated numbers.
@@ -79,7 +115,7 @@
                               double *rho, double *compositions)
 {
     char buffer[256], *p;
-    double temp[MAX_NUM_COMPOSITIONS + 5];
+    double temp[MAX_NUM_COMPOSITIONS + NUM_REQUIRED_COLUMNS];
     int c, i, j;
 
     p = fgets(buffer, 255, in);
@@ -87,7 +123,7 @@
         return 0;
     }
 
-    c = scan_double_vector(buffer, num_compositions, temp);
+    c = scan_double_vector(buffer, num_compositions+NUM_REQUIRED_COLUMNS, temp);
 
     *rank = (int) temp[0];
     *glob = (int) temp[1];
@@ -95,9 +131,9 @@
     *pressure = temp[3];
     *rho = temp[4];
 
-    for (i = 0, j = 5; i < num_compositions; ++i, ++j) {
+    for (i = 0, j = NUM_REQUIRED_COLUMNS; i < num_compositions; ++i, ++j) {
         compositions[i] = temp[j];
-    }    
+    }
     return c;
 }
 
@@ -106,12 +142,12 @@
 {
     FILE **out; int nout;
     int i;
-    
+
     if (argc < 2) {
         fprintf(stderr, "usage: %s FILE...\n", argv[0]);
         return 1;
     }
-    
+
     nout = 0;
     out = 0;
     for (i = 1; i < argc; ++i) {
@@ -119,6 +155,7 @@
         char buffer[256], *p;
         FILE *in;
         int c, num_columns;
+        int junk1, junk2;
         int num_compositions;
         int rank, glob;
         double temperature, pressure, rho;
@@ -127,13 +164,14 @@
 
         /* open input file(s) */
         ifilename = argv[i];
+        fprintf(stderr, "\rReading %s ...                  ", ifilename);
         in = fopen(ifilename, "r");
         if (!in) {
             fprintf(stderr, "%s: fopen: `%s': ", argv[0], ifilename);
             perror("");
             return 1;
         }
-        
+
         /* read header */
         p = fgets(buffer, 255, in);
         if (!p) {
@@ -141,24 +179,28 @@
             return 1;
         }
 
-        c = scanf(buffer, "%* %d", &num_columns);
-        if (c != 2) {
+        c = sscanf(buffer, "%d %d %d", &junk1, &junk2, &num_columns);
+        if (c != 3) {
             fprintf(stderr, "%s: bad input file: `%s'\n", argv[0], ifilename);
             return 1;
         }
 
         /* except the first five columns, the rest columns are compositions */
-        num_compositions = num_columns - 5;
+        num_compositions = num_columns - NUM_REQUIRED_COLUMNS;
         if (num_compositions > MAX_NUM_COMPOSITIONS) {
             fprintf(stderr, "%s: too many compositions in input file: `%s'\n", argv[0], ifilename);
             return 1;
         }
+        if (num_compositions < 0) {
+            fprintf(stderr, "%s: too few columns in input file: `%s'\n", argv[0], ifilename);
+            return 1;
+        }
 
         while (read_convection_variables(in, num_compositions,
                                          &rank, &glob,
                                          &temperature, &pressure,
                                          &rho, compositions)
-               >= (5 + num_compositions) ) {
+               >= (NUM_REQUIRED_COLUMNS + num_compositions) ) {
             int nproc;
             
             /* open more output files as needed */
@@ -186,6 +228,7 @@
         fclose(out[nout]);
     }
     free(out);
+    fprintf(stderr, "\nDone\n");
 
     return 0;
 }