[cig-commits] r9206 - cs/avm/trunk
tan2 at geodynamics.org
tan2 at geodynamics.org
Fri Feb 1 11:47:33 PST 2008
Author: tan2
Date: 2008-02-01 11:47:33 -0800 (Fri, 01 Feb 2008)
New Revision: 9206
Added:
cs/avm/trunk/avm-util.c
cs/avm/trunk/avm-util.h
cs/avm/trunk/elastic-models.c
cs/avm/trunk/elastic-models.h
cs/avm/trunk/prem-derived-models.c
cs/avm/trunk/prem-derived-models.h
Modified:
cs/avm/trunk/Makefile
cs/avm/trunk/avm-partition.c
Log:
Support of multiple ways to convert a convection model to a elastic model.
* An integer "model_type" from the tracer input indicates the type of convection model.
* Based on "model_type", a function pointer for the elastic model and the appropriate model struct are returned from select_elastic_model().
* The elastic model interpretes the convection model, ie the meaning of each column, and converts it to density and bulk and shear moduli.
* Avaiable models: given (T, C) and (dA/dT, dA/dC), computes deviation of A from PREM. A can be elastic moduli (model_type=0) or seismic velocities (model_type=1).
Modified: cs/avm/trunk/Makefile
===================================================================
--- cs/avm/trunk/Makefile 2008-02-01 19:28:55 UTC (rev 9205)
+++ cs/avm/trunk/Makefile 2008-02-01 19:47:33 UTC (rev 9206)
@@ -48,8 +48,8 @@
@rm -f avm-partition avm write_CitcomS_tracers
-avm-partition: avm-partition.c
- $(CC) $(CFLAGS) $(PFLAGS) -o $@ $@.c
+avm-partition: avm-partition.c elastic-models.c elastic-models.h avm-util.c avm-util.h prem-derived-models.c prem-derived-models.h
+ $(CC) $(CFLAGS) -o $@ avm-partition.c elastic-models.c avm-util.c prem-derived-models.c
avm: avm.f90
Modified: cs/avm/trunk/avm-partition.c
===================================================================
--- cs/avm/trunk/avm-partition.c 2008-02-01 19:28:55 UTC (rev 9205)
+++ cs/avm/trunk/avm-partition.c 2008-02-01 19:47:33 UTC (rev 9206)
@@ -29,15 +29,9 @@
#include <stdio.h>
#include <stdlib.h>
+#include "elastic-models.h"
+#include "avm-util.h"
-/* 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
-
/* a constant for specfem mesh data */
#define IREGION_CRUST_MANTLE 1
@@ -46,101 +40,33 @@
/****************************************************************************/
-/*
- * This function takes the temperatur (in Kelvin), pressure (in Pascal), and
- * composition (a vector of concentration) at a point, and returns the
- * density, bulk and shear moduli (in Pascal).
- */
-void elastic_model(double temperature, double pressure, double rho,
- double *compositions, int num_compositions,
- double *ks, double *mu)
-{
- /* XXX: user modification here */
- *ks = *mu = 1.0;
- 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.
- */
-int scan_double_vector(const char *str, int len, double *values)
-{
- char *nptr, *endptr;
- int i;
-
- /* cast to avoid compiler warning */
- nptr = endptr = (char *) str;
-
- for (i = 0; i < len; ++i) {
- values[i] = strtod(nptr, &endptr);
- if (nptr == endptr) {
- /* error: no conversion is performed */
- return -1;
- }
- nptr = endptr;
- }
-
- return len;
-}
-
-
-int read_convection_variables(FILE *in, int num_compositions,
- int *rank, int *glob,
- double *temperature, double *pressure,
- double *rho, double *compositions)
-{
- char buffer[256], *p;
- double temp[MAX_NUM_COMPOSITIONS + NUM_REQUIRED_COLUMNS];
- int c, i, j;
-
- p = fgets(buffer, 255, in);
- if (!p) {
- return 0;
- }
-
- c = scan_double_vector(buffer, num_compositions+NUM_REQUIRED_COLUMNS, temp);
-
- *rank = (int) temp[0];
- *glob = (int) temp[1];
- *temperature = temp[2];
- *pressure = temp[3];
- *rho = temp[4];
-
- for (i = 0, j = NUM_REQUIRED_COLUMNS; i < num_compositions; ++i, ++j) {
- compositions[i] = temp[j];
- }
- return c;
-}
-
-
int main(int argc, char **argv)
{
FILE **out; int nout;
- int i;
+ int i, old_type, old_num_columns, old_num_compositions;
+ elastic_model_fn elastic_model = NULL;
+ void *context = NULL;
if (argc < 2) {
fprintf(stderr, "usage: %s FILE...\n", argv[0]);
return 1;
}
+ old_type = old_num_columns = old_num_compositions = -1;
nout = 0;
- out = 0;
+ out = NULL;
for (i = 1; i < argc; ++i) {
char *ifilename;
char buffer[256], *p;
FILE *in;
- int c, num_columns;
- int junk1, junk2;
- int num_compositions;
- int rank, glob;
- double temperature, pressure, rho;
- double compositions[MAX_NUM_COMPOSITIONS];
- double ks, mu;
+ int c, junk;
+ int model_type, num_columns, num_compositions;
+ double *fields;
+ double rho, ks, mu;
/* open input file(s) */
ifilename = argv[i];
- fprintf(stderr, "\rReading %s ... ", ifilename);
+ fprintf(stderr, "\rReading %s ... ", ifilename);
in = fopen(ifilename, "r");
if (!in) {
fprintf(stderr, "%s: fopen: `%s': ", argv[0], ifilename);
@@ -155,30 +81,46 @@
return 1;
}
- c = sscanf(buffer, "%d %d %d", &junk1, &junk2, &num_columns);
- if (c != 3) {
+ c = sscanf(buffer, "%d %d %d %d", &junk, &model_type,
+ &num_columns, &num_compositions);
+
+ /* validate header */
+ if (c != 4) {
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 - 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) {
+ if (num_compositions < 0 || num_compositions < 0) {
fprintf(stderr, "%s: too few columns in input file: `%s'\n", argv[0], ifilename);
return 1;
}
+ if (old_type == -1) {
+ /* first time in loop */
+ old_type = model_type;
+ old_num_compositions = num_compositions;
- while (read_convection_variables(in, num_compositions,
- &rank, &glob,
- &temperature, &pressure,
- &rho, compositions)
- >= (NUM_REQUIRED_COLUMNS + num_compositions) ) {
- int nproc;
-
+ /* select and setup elastic model for later use */
+ select_elastic_model(model_type, num_columns, num_compositions,
+ elastic_model, context);
+ }
+ else {
+ if ((old_type != model_type) ||
+ (old_num_columns != num_columns) ||
+ (old_num_compositions != num_compositions)) {
+ fprintf(stderr, "%s: inconsistent header with previous input in input file: `%s'\n", argv[0], ifilename);
+ return 1;
+ }
+ }
+ fields = malloc(num_columns * sizeof(double));
+
+ /* read content */
+ while (read_double_vector(in, num_columns, fields)
+ == num_columns) {
+ int nproc, rank, iglob, error;
+
+ /* the 1st and 2nd columns are always rank and iglob */
+ rank = (int) fields[0];
+ iglob = (int) fields[1];
+
/* open more output files as needed */
nproc = rank + 1;
if (nproc > nout) {
@@ -190,13 +132,17 @@
}
}
- /* convert convection variables to elastic moduli */
- elastic_model(temperature, pressure, rho,
- compositions, num_compositions,
- &ks, &mu);
- fprintf(out[rank], "%d %lf %lf %lf\n", glob, rho, ks, mu);
+ /* convert convection fields to density and elastic moduli */
+ error = elastic_model(context, fields,
+ &rho, &ks, &mu);
+ if (error) {
+ fprintf(stderr, "%s: error while converting `%s' to elastic moduli\n", argv[0], ifilename);
+ exit(1);
+ }
+ fprintf(out[rank], "%d %lf %lf %lf\n", iglob, rho, ks, mu);
}
-
+
+ free(fields);
fclose(in);
}
Added: cs/avm/trunk/avm-util.c
===================================================================
--- cs/avm/trunk/avm-util.c (rev 0)
+++ cs/avm/trunk/avm-util.c 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,174 @@
+
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * avm-util.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "avm-util.h"
+
+/*
+ * Scan input str to get a double vector *values. The vector length is from
+ * input len. The input str contains white-space seperated numbers.
+ */
+static int scan_double_vector(const char *str, int len, double *values)
+{
+ char *nptr, *endptr;
+ int i;
+
+ /* cast to avoid compiler warning */
+ nptr = endptr = (char *) str;
+
+ for (i = 0; i < len; ++i) {
+ values[i] = strtod(nptr, &endptr);
+ if (nptr == endptr) {
+ /* error: no conversion is performed */
+ return -1;
+ }
+ nptr = endptr;
+ }
+
+ /** debug **
+ for (i = 0; i < len; ++i) fprintf(stderr, "%e, ", values[i]);
+ fprintf(stderr, "\n");
+ /**/
+ return len;
+}
+
+
+/*
+ * From input file, read a line, which contains white-space seperated numbers
+ * of lenght num_columns, store the numbers in a double array.
+ */
+int read_double_vector(FILE *in, int num_columns, double *fields)
+{
+ char buffer[256], *p;
+
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ return 0;
+ }
+
+ return scan_double_vector(buffer, num_columns, fields);
+}
+
+
+/*
+ * Given a radius (in meter), returns Vp/Vs/rho of PREM.
+ */
+void get_prem(double radius, double *vp, double *vs, double *rho)
+{
+#define NUM_PREM_LAYERS 11
+
+ /* radius of various layers */
+ const double prem_radius[NUM_PREM_LAYERS] =
+ {1221.5,
+ 3480.0,
+ 3630.0,
+ 5600.0,
+ 5701.0,
+ 5771.0,
+ 5971.0,
+ 6151.0,
+ 6345.6,
+ 6356.0,
+ 6371.0};
+
+ /* polynomial coefficients of PREM */
+ const double prem_vs[NUM_PREM_LAYERS][4] =
+ {{ 3.6678, 0.0000, -4.4475, 0.0000},
+ { 0.0010, 0.0000, 0.0000, 0.0000},
+ { 6.9254, 1.4672, -2.0834, 0.9783},
+ {11.1671, -13.7818, 17.4575, -9.2777},
+ {22.3459, -17.2473, -2.0834, 0.9783},
+ { 9.9839, -4.9324, 0.0000, 0.0000},
+ {22.3512, -18.5856, 0.0000, 0.0000},
+ { 8.9496, -4.4597, 0.0000, 0.0000},
+ { 2.1519, 2.3481, 0.0000, 0.0000},
+ { 3.9000, 0.0000, 0.0000, 0.0000},
+ { 3.2000, 0.0000, 0.0000, 0.0000}};
+
+ const double prem_vp[NUM_PREM_LAYERS][4] =
+ {{11.2622, 0.0000, -6.3640, 0.0000},
+ {11.0487, -4.0362, 4.8023,-13.5732},
+ {15.3891, -5.3181, 5.5242, -2.5514},
+ {24.9520, -40.4673, 51.4832,-26.6419},
+ {29.2766, -23.6027, 5.5242, -2.5514},
+ {19.0957, -9.8672, 0.0000, 0.0000},
+ {39.7027, -32.6166, 0.0000, 0.0000},
+ {20.3926, -12.2569, 0.0000, 0.0000},
+ { 4.1875, 3.9382, 0.0000, 0.0000},
+ { 6.8000, 0.0000, 0.0000, 0.0000},
+ { 5.8000, 0.0000, 0.0000, 0.0000}};
+
+ const double prem_rho[NUM_PREM_LAYERS][4] =
+ {{13.0885, 0.0000, -8.8381, 0.0000},
+ {12.5815, -1.2638, -3.6426, -5.5281},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 7.9565, -6.4761, 5.5283, -3.0807},
+ { 5.3197, -1.4836, 0.0000, 0.0000},
+ {11.2494, -8.0298, 0.0000, 0.0000},
+ { 7.1089, -3.8045, 0.0000, 0.0000},
+ { 2.6910, 0.6924, 0.0000, 0.0000},
+ { 2.9000, 0.0000, 0.0000, 0.0000},
+ { 2.6000, 0.0000, 0.0000, 0.0000}};
+
+
+ /* normalized radius */
+ double r = radius / prem_radius[NUM_PREM_LAYERS - 1];
+ double r2 = r * r;
+ double r3 = r2 * r;
+ int j;
+
+
+ /* find layer */
+ for (j = 0; j < NUM_PREM_LAYERS; ++j)
+ if (radius < prem_radius[j]) break;
+
+ if (j < 0) j = 0;
+ if (j >= NUM_PREM_LAYERS) j = NUM_PREM_LAYERS - 1;
+
+ /* expand polynomials */
+ *vp = prem_vp[j][0]
+ + prem_vp[j][1] * r
+ + prem_vp[j][2] * r2
+ + prem_vp[j][3] * r3;
+ *vs = prem_vs[j][0]
+ + prem_vs[j][1] * r
+ + prem_vs[j][2] * r2
+ + prem_vs[j][3] * r3;
+ *rho = prem_rho[j][0]
+ + prem_rho[j][1] * r
+ + prem_rho[j][2] * r2
+ + prem_rho[j][3] * r3;
+
+#undef NUM_PREM_LAYERS
+
+ return;
+}
+
+
Added: cs/avm/trunk/avm-util.h
===================================================================
--- cs/avm/trunk/avm-util.h (rev 0)
+++ cs/avm/trunk/avm-util.h 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,31 @@
+
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * avm-util.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+
+int read_double_vector(FILE *in, int num_columns, double *fields);
+void get_prem(double radius, double *vp, double *vs, double *rho);
Added: cs/avm/trunk/elastic-models.c
===================================================================
--- cs/avm/trunk/elastic-models.c (rev 0)
+++ cs/avm/trunk/elastic-models.c 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,87 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * elastic-models.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "elastic-models.h"
+#include "prem-derived-models.h"
+
+/****************************************************************************/
+
+/*
+ * These functions takes a vector fileds, whis is num_columns in length and
+ * the last num_compositions entries are compositions, and returns the
+ * density, bulk and shear moduli. All quantities are in SI units.
+ */
+
+/****************************************************************************/
+
+int user_defined_em(void *context, double *fields,
+ double *rho, double *ks, double *mu)
+{
+ /* define your own elastic model here */
+ *rho = *ks = *mu = 1.0;
+
+ return 0;
+}
+
+/****************************************************************************/
+/*
+ * choose which elastic model to use
+ */
+void select_elastic_model(int model_type,
+ int num_columns, int num_compositions,
+ elastic_model_fn fn, void *context)
+{
+ /*
+ * Select which elastic model to use. The convection is:
+ *
+ * model_type 0-9 are based on PREM.
+ * model_type 100 and later are reserved for user modification.
+ * other model_type are not defined yet.
+ */
+ switch (model_type) {
+ case 0:
+ /* using modulus derivatives */
+ prem_mod_derv_em_init(num_columns, num_compositions, &context);
+ fn = prem_mod_derv_em;
+ return;
+ case 1:
+ /* using velocity derivatives */
+ prem_vel_derv_em_init(num_columns, num_compositions, &context);
+ fn = prem_vel_derv_em;
+ return;
+ case 100:
+ /* user defined */
+ fn = user_defined_em;
+ return;
+ default:
+ fprintf(stderr, "Error: unknown type of elastic model\n");
+ exit(1);
+ }
+ return;
+}
Added: cs/avm/trunk/elastic-models.h
===================================================================
--- cs/avm/trunk/elastic-models.h (rev 0)
+++ cs/avm/trunk/elastic-models.h 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,34 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * elastic-models.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+/* a function pointer for various elastic models */
+typedef int (*elastic_model_fn)(void*, double*, double*, double*, double*);
+
+
+void select_elastic_model(int model_type,
+ int num_columns, int num_compositions,
+ elastic_model_fn fn, void *context);
Added: cs/avm/trunk/prem-derived-models.c
===================================================================
--- cs/avm/trunk/prem-derived-models.c (rev 0)
+++ cs/avm/trunk/prem-derived-models.c 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,394 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * prem-derived-models.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "avm-util.h"
+#include "prem-derived-models.h"
+
+
+/****************************************************************************/
+/*
+ * Helper functions
+ */
+
+static void read_header(FILE *in, const char *ifilename,
+ int num_compositions, int *num_layers)
+{
+ char buffer[256], *p;
+ int c, itmp ;
+
+ /* read header, the format is:
+ * [num_layers, num_compositions]
+ */
+ p = fgets(buffer, 255, in);
+ if (!p) {
+ fprintf(stderr, "Bad input file: `%s'\n", ifilename);
+ exit(2);
+ }
+
+ c = sscanf(buffer, "%d %d", num_layers, &itmp);
+
+ /* validate header */
+ if (c != 2) {
+ fprintf(stderr, "Bad input file: `%s'\n", ifilename);
+ exit(2);
+ }
+ if (itmp < num_compositions) {
+ fprintf(stderr, "Insufficient # of compositions in input file: `%s'\n", ifilename);
+ exit(2);
+ }
+
+ return;
+}
+
+
+/****************************************************************************/
+/*
+ * This model uses temperature and composition derivatives of elastic moduli
+ * to compute elastic moduli from PREM.
+ *
+ * The derivatives are read from file "modulus-derivatives.dat".
+ *
+ * The input fields are:
+ * [rank, iglob, radius, temperature, composition(s)]
+ * The temperature and compositions fields are the deviations from the radial
+ * profile.
+ */
+
+typedef struct {
+ int num_layers;
+ int num_compositions;
+
+ double *layer_radius;
+
+ /* temperature derivatives */
+ double *layer_drhodt;
+ double *layer_dksdt;
+ double *layer_dmudt;
+
+ /* composition derivatives */
+ double **layer_drhodc;
+ double **layer_dksdc;
+ double **layer_dmudc;
+
+} prem_modulus_derv_t;
+
+
+void prem_mod_derv_em_init(int num_columns, int num_compositions,
+ void **context)
+{
+ /* see the comment of fields format in the next function */
+ const int num_required_fields = num_compositions + 5;
+
+ int i, j, k;
+ int total_columns, num_layers;
+ int itmp;
+ double *tmp;
+ const char ifilename[] = "modulus-derivatives.dat";
+ FILE *in;
+ prem_modulus_derv_t *p;
+
+ /* if *p is NULL, we need to set up the context */
+ if (!(*context)) return;
+ *context = p = malloc(sizeof(prem_modulus_derv_t));
+ p->num_compositions = num_compositions;
+
+ /* check the # of input fields */
+ if (num_columns != num_required_fields) {
+ fprintf(stderr, "Error: %d columns are required in field input\n",
+ num_required_fields);
+ exit(2);
+ }
+
+ /* read d*dc, d*dt from file, the filename is hard coded */
+ in = fopen(ifilename, "r");
+ if (!in) {
+ fprintf(stderr, "Cannot open input `%s': ", ifilename);
+ perror("");
+ exit(2);
+ }
+
+ read_header(in, ifilename, num_compositions, &num_layers);
+ p->num_layers = num_layers;
+
+ /* allocate memory */
+ p->layer_radius = malloc(num_layers * sizeof(double));
+ p->layer_dksdt = malloc(num_layers * sizeof(double));
+ p->layer_dmudt = malloc(num_layers * sizeof(double));
+ p->layer_drhodt = malloc(num_layers * sizeof(double));
+
+ p->layer_dksdc = malloc(num_layers * sizeof(double*));
+ p->layer_dmudc = malloc(num_layers * sizeof(double*));
+ p->layer_drhodc = malloc(num_layers * sizeof(double*));
+ for (i = 0; i < num_layers; ++i) {
+ p->layer_dksdc[i] = malloc(num_compositions * sizeof(double));
+ p->layer_dmudc[i] = malloc(num_compositions * sizeof(double));
+ p->layer_drhodc[i] = malloc(num_compositions * sizeof(double));
+ }
+
+
+ /* read content, the file format is:
+ * [radius, dks/dt, dmu/dt, drho/dt,
+ * dks/dc1, dmu/dc1, drho/dc1,
+ * dks/dc2, dmu/dc2, drho/dc2, ...]
+ */
+ total_columns = 4 + 3 * num_compositions;
+ tmp = malloc(total_columns * sizeof(double));
+
+ for (i = 0; i < num_layers; ++i) {
+ itmp = read_double_vector(in, total_columns, tmp);
+ if (itmp != total_columns) {
+ fprintf(stderr, "Bad line in input file: `%s'\n", ifilename);
+ exit(2);
+ }
+
+ p->layer_radius[i] = tmp[0];
+ p->layer_dksdt[i] = tmp[1];
+ p->layer_dmudt[i] = tmp[2];
+ p->layer_drhodt[i] = tmp[3];
+
+ for (j = 0, k = 4; j < num_compositions; ++j, k += 3) {
+ p->layer_dksdc[i][j] = tmp[k];
+ p->layer_dmudc[i][j] = tmp[k+1];
+ p->layer_drhodc[i][j] = tmp[k+2];
+ }
+ }
+ fclose(in);
+
+ free(tmp);
+ return;
+}
+
+
+int prem_mod_derv_em(void *context, double *fields,
+ double *rho, double *ks, double *mu)
+{
+ int i, j;
+ double radius, temperature, *compositions;
+ double vp0, vs0, rho0;
+ double ks0, mu0;
+ prem_modulus_derv_t *p;
+
+ p = context;
+
+ /* 5 fields are required, the 1st two are not used here */
+ radius = fields[2];
+ temperature = fields[3];
+ compositions = &(fields[4]);
+
+ /* get PREM value */
+ get_prem(radius, &vp0, &vs0, &rho0);
+
+ ks0 = rho0 * (vp0*vp0 - 4.0/3.0*vs0*vs0);
+ mu0 = rho0 * vs0 * vs0;
+
+ /* find layer */
+ for (j = 0; j < p->num_layers; ++j)
+ if (radius < p->layer_radius[j]) break;
+
+ if (j < 0) j = 0;
+ if (j >= p->num_layers) j = p->num_layers - 1;
+
+ /* contribution from temperature */
+ *ks = ks0 + temperature * p->layer_dksdt[j];
+ *mu = mu0 + temperature * p->layer_dmudt[j];
+ *rho = rho0 + temperature * p->layer_drhodt[j];
+
+ /* contribution from compositions */
+ for (i = 0; i < p->num_compositions; ++i) {
+ *ks += compositions[i] * p->layer_dksdc[j][i];
+ *mu += compositions[i] * p->layer_dmudc[j][i];
+ *rho += compositions[i] * p->layer_drhodc[j][i];
+ }
+
+ return 0;
+}
+
+
+/****************************************************************************/
+/*
+ * This model uses temperature and composition derivatives of seismic velocity
+ * to compute elastic moduli from PREM.
+ *
+ * The derivatives are read from file "velocity-derivatives.dat".
+ *
+ * The input fields are:
+ * [rank, iglob, radius, temperature, composition(s)]
+ * The temperature and compositions fields are the deviations from the radial
+ * profile.
+ */
+
+typedef struct {
+ int num_layers;
+ int num_compositions;
+
+ double *layer_radius;
+
+ /* temperature derivatives */
+ double *layer_drhodt;
+ double *layer_dvpdt;
+ double *layer_dvsdt;
+
+ /* composition derivatives */
+ double **layer_drhodc;
+ double **layer_dvpdc;
+ double **layer_dvsdc;
+
+} prem_velocity_derv_t;
+
+
+void prem_vel_derv_em_init(int num_columns, int num_compositions,
+ void **context)
+{
+ /* see the comment of fields format in the next function */
+ const int num_required_fields = num_compositions + 5;
+
+ int i, j, k;
+ int total_columns, num_layers;
+ int itmp;
+ double *tmp;
+ const char ifilename[] = "velocity-derivatives.dat";
+ FILE *in;
+ prem_velocity_derv_t *p;
+
+ /* if *p is NULL, we need to set up the context */
+ if (!(*context)) return;
+ *context = p = malloc(sizeof(prem_velocity_derv_t));
+ p->num_compositions = num_compositions;
+
+ /* check the # of input fields */
+ if (num_columns != num_required_fields) {
+ fprintf(stderr, "Error: %d columns are required in field input\n",
+ num_required_fields);
+ exit(2);
+ }
+
+ /* read d*dc, d*dt from file, the filename is hard coded */
+ in = fopen(ifilename, "r");
+ if (!in) {
+ fprintf(stderr, "Cannot open input `%s': ", ifilename);
+ perror("");
+ exit(2);
+ }
+
+ read_header(in, ifilename, num_compositions, &num_layers);
+ p->num_layers = num_layers;
+
+ /* allocate memory */
+ p->layer_radius = malloc(num_layers * sizeof(double));
+ p->layer_dvpdt = malloc(num_layers * sizeof(double));
+ p->layer_dvsdt = malloc(num_layers * sizeof(double));
+ p->layer_drhodt = malloc(num_layers * sizeof(double));
+
+ p->layer_dvpdc = malloc(num_layers * sizeof(double*));
+ p->layer_dvsdc = malloc(num_layers * sizeof(double*));
+ p->layer_drhodc = malloc(num_layers * sizeof(double*));
+ for (i = 0; i < num_layers; ++i) {
+ p->layer_dvpdc[i] = malloc(num_compositions * sizeof(double));
+ p->layer_dvsdc[i] = malloc(num_compositions * sizeof(double));
+ p->layer_drhodc[i] = malloc(num_compositions * sizeof(double));
+ }
+
+
+ /* read content, the file format is:
+ * [radius, dvp/dt, dvs/dt, drho/dt,
+ * dvp/dc1, dvs/dc1, drho/dc1,
+ * dvp/dc2, dvs/dc2, drho/dc2, ...]
+ */
+ total_columns = 4 + 3 * num_compositions;
+ tmp = malloc(total_columns * sizeof(double));
+
+ for (i = 0; i < num_layers; ++i) {
+ itmp = read_double_vector(in, total_columns, tmp);
+ if (itmp != total_columns) {
+ fprintf(stderr, "Bad line in input file: `%s'\n", ifilename);
+ exit(2);
+ }
+
+ p->layer_radius[i] = tmp[0];
+ p->layer_dvpdt[i] = tmp[1];
+ p->layer_dvsdt[i] = tmp[2];
+ p->layer_drhodt[i] = tmp[3];
+
+ for (j = 0, k = 4; j < num_compositions; ++j, k += 3) {
+ p->layer_dvpdc[i][j] = tmp[k];
+ p->layer_dvsdc[i][j] = tmp[k+1];
+ p->layer_drhodc[i][j] = tmp[k+2];
+ }
+ }
+ fclose(in);
+
+ free(tmp);
+ return;
+}
+
+
+int prem_vel_derv_em(void *context, double *fields,
+ double *rho, double *ks, double *mu)
+{
+ int i, j;
+ double radius, temperature, *compositions;
+ double vp0, vs0, rho0;
+ double vp, vs;
+ prem_velocity_derv_t *p;
+
+ p = context;
+
+ radius = fields[2];
+ temperature = fields[3];
+ compositions = &(fields[4]);
+
+ /* get PREM value */
+ get_prem(radius, &vp0, &vs0, &rho0);
+
+ /* find layer */
+ for (j = 0; j < p->num_layers; ++j)
+ if (radius < p->layer_radius[j]) break;
+
+ if (j < 0) j = 0;
+ if (j >= p->num_layers) j = p->num_layers - 1;
+
+ /* contribution from temperature */
+ vp = vp0 + temperature * p->layer_dvpdt[j];
+ vs = vs0 + temperature * p->layer_dvsdt[j];
+ *rho = rho0 + temperature * p->layer_drhodt[j];
+
+ /* contribution from compositions */
+ for (i = 0; i < p->num_compositions; ++i) {
+ *ks += compositions[i] * p->layer_dvpdc[j][i];
+ *mu += compositions[i] * p->layer_dvsdc[j][i];
+ *rho += compositions[i] * p->layer_drhodc[j][i];
+ }
+
+ *ks = *rho * (vp*vp - 4.0/3.0*vs*vs);
+ *mu = *rho * vs * vs;
+
+ return 0;
+}
+
+
Added: cs/avm/trunk/prem-derived-models.h
===================================================================
--- cs/avm/trunk/prem-derived-models.h (rev 0)
+++ cs/avm/trunk/prem-derived-models.h 2008-02-01 19:47:33 UTC (rev 9206)
@@ -0,0 +1,36 @@
+/*
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ *<LicenseText>
+ *
+ * prem-derived-models.c by Leif Strand and Eh Tan
+ * Copyright (C) 2007, California Institute of Technology.
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ *
+ *</LicenseText>
+ *
+ *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ */
+
+
+void prem_mod_derv_em_init(int num_columns, int num_compositions,
+ void **context);
+int prem_mod_derv_em(void *context, double *fields,
+ double *rho, double *ks, double *mu);
+
+void prem_vel_derv_em_init(int num_columns, int num_compositions,
+ void **context);
+int prem_vel_derv_em(void *context, double *fields,
+ double *rho, double *ks, double *mu);
More information about the cig-commits
mailing list