[cig-commits] r8021 - in mc/3D/CitcomS/trunk: CitcomS/Components
lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Sep 24 18:31:31 PDT 2007
Author: tan2
Date: 2007-09-24 18:31:30 -0700 (Mon, 24 Sep 2007)
New Revision: 8021
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Param.py
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Material_properties.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Introducing choice of reference state:
When solver.param.reference_state=0, the reference state is read from a file (solver.param.refstate_file)
When solver.param.reference_state=1 (default), the reference state is calculated using Adams-Williamson equation of state.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Param.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Param.py 2007-09-25 01:21:03 UTC (rev 8020)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Param.py 2007-09-25 01:31:30 UTC (rev 8021)
@@ -50,6 +50,9 @@
import pyre.inventory
+ reference_state = pyre.inventory.int("reference_state", default=1)
+ refstate_file = pyre.inventory.str("refstate_file", default="refstate.dat")
+
file_vbcs = pyre.inventory.bool("file_vbcs", default=False)
vel_bound_file = pyre.inventory.str("vel_bound_file", default="bvel.dat")
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-09-25 01:21:03 UTC (rev 8020)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-09-25 01:31:30 UTC (rev 8021)
@@ -401,6 +401,11 @@
input_int("file_vbcs",&(E->control.vbcs_file),"0",m);
input_string("vel_bound_file",E->control.velocity_boundary_file,"",m);
+ input_int("reference_state",&(E->refstate.choice),"1",m);
+ if(E->refstate.choice == 0) {
+ input_string("refstate_file",E->refstate.filename,"refstate.dat",m);
+ }
+
input_int("mat_control",&(E->control.mat_control),"0",m);
input_string("mat_file",E->control.mat_file,"",m);
Modified: mc/3D/CitcomS/trunk/lib/Material_properties.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Material_properties.c 2007-09-25 01:21:03 UTC (rev 8020)
+++ mc/3D/CitcomS/trunk/lib/Material_properties.c 2007-09-25 01:31:30 UTC (rev 8021)
@@ -34,8 +34,12 @@
#include "global_defs.h"
#include "material_properties.h"
+#include "parallel_related.h"
+static void read_refstate(struct All_variables *E);
+static void adams_williamson_eos(struct All_variables *E);
+
void mat_prop_allocate(struct All_variables *E)
{
int noz = E->lmesh.noz;
@@ -45,6 +49,9 @@
/* reference profile of density */
E->refstate.rho = (double *) malloc((noz+1)*sizeof(double));
+ /* reference profile of gravity */
+ E->refstate.gravity = (double *) malloc((noz+1)*sizeof(double));
+
/* reference profile of coefficient of thermal expansion */
E->refstate.thermal_expansivity = (double *) malloc((noz+1)*sizeof(double));
@@ -54,9 +61,6 @@
/* reference profile of thermal conductivity */
/*E->refstate.thermal_conductivity = (double *) malloc((noz+1)*sizeof(double));*/
- /* reference profile of gravity */
- E->refstate.gravity = (double *) malloc((noz+1)*sizeof(double));
-
/* reference profile of temperature */
/*E->refstate.Tadi = (double *) malloc((noz+1)*sizeof(double));*/
@@ -65,38 +69,108 @@
void reference_state(struct All_variables *E)
{
- int noz = E->lmesh.noz;
- int nel = E->lmesh.nel;
int i;
- double r, z, beta;
- beta = E->control.disptn_number * E->control.inv_gruneisen;
-
/* All refstate variables (except Tadi) must be 1 at the surface.
* Otherwise, the scaling of eqns in the code might not be correct. */
- /* Adams-Williamson EoS */
- for(i=1; i<=noz; i++) {
- r = E->sx[1][3][i];
- z = 1 - r;
- E->refstate.rho[i] = exp(beta*z);
- E->refstate.thermal_expansivity[i] = 1;
- E->refstate.heat_capacity[i] = 1;
- /*E->refstate.thermal_conductivity[i] = 1;*/
- E->refstate.gravity[i] = 1;
- /*E->refstate.Tadi[i] = (E->control.adiabaticT0 + E->control.surface_temp) * exp(E->control.disptn_number * z) - E->control.surface_temp;*/
+ /* select the choice of reference state */
+ switch(E->refstate.choice) {
+ case 0:
+ /* read from a file */
+ read_refstate(E);
+ break;
+ case 1:
+ /* Adams-Williamson EoS */
+ adams_williamson_eos(E);
+ break;
+ default:
+ if (E->parallel.me) {
+ fprintf(stderr, "Unknown option for reference state\n");
+ fprintf(E->fp, "Unknown option for reference state\n");
+ fflush(E->fp);
+ }
+ parallel_process_termination();
}
if(E->parallel.me == 0) {
fprintf(stderr, "nz radius depth rho\n");
}
if(E->parallel.me < E->parallel.nprocz)
- for(i=1; i<=noz; i++) {
+ for(i=1; i<=E->lmesh.noz; i++) {
fprintf(stderr, "%d %f %f %e\n",
i+E->lmesh.nzs-1, E->sx[1][3][i], 1-E->sx[1][3][i],
E->refstate.rho[i]);
}
+ return;
}
+static void read_refstate(struct All_variables *E)
+{
+ FILE *fp;
+ int i;
+ char buffer[255];
+ double not_used1, not_used2, not_used3;
+
+ fp = fopen(E->refstate.filename, "r");
+ if(fp == NULL) {
+ fprintf(stderr, "Cannot open reference state file: %s\n",
+ E->refstate.filename);
+ parallel_process_termination();
+ }
+
+ /* skip these lines, which belong to other processors */
+ for(i=1; i<E->lmesh.nzs; i++) {
+ fgets(buffer, 255, fp);
+ }
+
+ for(i=1; i<=E->lmesh.noz; i++) {
+ fgets(buffer, 255, fp);
+ sscanf(buffer, "%lf %lf %lf %lf %lf %lf %lf\n",
+ &(E->refstate.rho[i]),
+ &(E->refstate.gravity[i]),
+ &(E->refstate.thermal_expansivity[i]),
+ &(E->refstate.heat_capacity[i]),
+ ¬_used1,
+ ¬_used2,
+ ¬_used3);
+
+ /**** debug ****
+ fprintf(stderr, "%d %f %f %f %f\n",
+ i,
+ E->refstate.rho[i],
+ E->refstate.gravity[i],
+ E->refstate.thermal_expansivity[i],
+ E->refstate.heat_capacity[i]);
+ /* end of debug */
+ }
+
+ fclose(fp);
+ return;
+}
+
+
+static void adams_williamson_eos(struct All_variables *E)
+{
+ int i;
+ double r, z, beta;
+
+ beta = E->control.disptn_number * E->control.inv_gruneisen;
+
+ for(i=1; i<=E->lmesh.noz; i++) {
+ r = E->sx[1][3][i];
+ z = 1 - r;
+ E->refstate.rho[i] = exp(beta*z);
+ E->refstate.gravity[i] = 1;
+ E->refstate.thermal_expansivity[i] = 1;
+ E->refstate.heat_capacity[i] = 1;
+ /*E->refstate.thermal_conductivity[i] = 1;*/
+ /*E->refstate.Tadi[i] = (E->control.adiabaticT0 + E->control.surface_temp) * exp(E->control.disptn_number * z) - E->control.surface_temp;*/
+ }
+
+ return;
+}
+
+
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-09-25 01:21:03 UTC (rev 8020)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-09-25 01:31:30 UTC (rev 8021)
@@ -524,6 +524,8 @@
struct REF_STATE {
+ int choice;
+ char filename[200];
double *rho;
double *thermal_expansivity;
double *heat_capacity;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-09-25 01:21:03 UTC (rev 8020)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-09-25 01:31:30 UTC (rev 8021)
@@ -358,6 +358,11 @@
PUTS(("[CitcomS.solver.param]\n"));
+ getIntProperty(properties, "reference_state", E->refstate.choice, fp);
+ if(E->refstate.choice == 0) {
+ getStringProperty(properties, "refstate_file", E->refstate.filename, fp);
+ }
+
getIntProperty(properties, "file_vbcs", E->control.vbcs_file, fp);
getStringProperty(properties, "vel_bound_file", E->control.velocity_boundary_file, fp);
More information about the cig-commits
mailing list