[cig-commits] r5758 - in mc/3D/CitcomS/branches/compressible: lib
module
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Jan 11 13:01:07 PST 2007
Author: tan2
Date: 2007-01-11 13:01:06 -0800 (Thu, 11 Jan 2007)
New Revision: 5758
Modified:
mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
mc/3D/CitcomS/branches/compressible/lib/Instructions.c
mc/3D/CitcomS/branches/compressible/lib/Material_properties.c
mc/3D/CitcomS/branches/compressible/lib/Output_h5.c
mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
mc/3D/CitcomS/branches/compressible/lib/global_defs.h
mc/3D/CitcomS/branches/compressible/lib/material_properties.h
mc/3D/CitcomS/branches/compressible/module/mesher.c
mc/3D/CitcomS/branches/compressible/module/setProperties.c
Log:
1. Stored the Gruneisen parameter in its inverse
2. Calculated reference state
3. Removed contribution of ref. temperature from buoyancy
Modified: mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -78,8 +78,9 @@
const int nel=E->lmesh.nel;
const int lev=E->mesh.levmax;
- density(E, E->rho);
thermal_buoyancy(E,E->buoyancy);
+ //chemical_buoyancy(E,E->buoyancy);
+ //buoyancy_to_density(E, E->buoyancy, E->rho);
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -291,7 +292,7 @@
bdbmu[2][1]=bdbmu[2][2]=bdbmu[2][3]=
bdbmu[3][1]=bdbmu[3][2]=bdbmu[3][3]=0.0;
- if(E->control.gruneisen > 0)
+ if(E->control.inv_gruneisen > 0)
for(i=1;i<=dims;i++)
for(j=1;j<=dims;j++)
for(i1=1;i1<=6;i1++)
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -162,6 +162,8 @@
construct_sub_element(E);
construct_shape_functions(E);
+ reference_state(E);
+
/* construct_c3x3matrix(E); */ /* this matrix results from spherical geometry*/
mass_matrix(E);
@@ -366,7 +368,9 @@
input_float("rayleigh",&(E->control.Atemp),"essential",m);
input_float("dissipation",&(E->control.Di),"0.0",m);
- input_float("gruneisen",&(E->control.gruneisen),"0.0",m);
+ input_float("gruneisen",&(tmp),"0.0",m);
+ if(abs(tmp) > 1e-6)
+ E->control.inv_gruneisen = 1/tmp;
/* data section */
input_float("Q0",&(E->control.Q0),"0.0",m);
Modified: mc/3D/CitcomS/branches/compressible/lib/Material_properties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -30,6 +30,8 @@
#include "config.h"
#endif
+#include <math.h>
+
#include "global_defs.h"
#include "material_properties.h"
@@ -39,32 +41,47 @@
int noz = E->lmesh.noz;
int nno = E->lmesh.nno;
- /* reference density profile */
- E->rhoref = (double *) malloc((noz+1)*sizeof(double));
+ /* reference profile of density */
+ E->rho_ref = (double *) malloc((noz+1)*sizeof(double));
- /* coefficient of thermal expansion */
- E->thermexp = (double *) malloc((noz+1)*sizeof(double));
+ /* reference profile of coefficient of thermal expansion */
+ E->thermexp_ref = (double *) malloc((noz+1)*sizeof(double));
-
+ /* reference profile of temperature */
+ E->T_ref = (double *) malloc((noz+1)*sizeof(double));
}
-void density(struct All_variables *E, double *rho)
+void reference_state(struct All_variables *E)
{
+ int noz = E->lmesh.noz;
int i;
- for(i=1; i<=E->lmesh.nno; i++) {
- rho[i] = 1;
+ double r, z, tmp, T0;
+
+ tmp = E->control.Di * E->control.inv_gruneisen;
+ T0 = E->data.surf_temp / E->data.ref_temperature;
+
+ for(i=1; i<=noz; i++) {
+ r = E->sx[1][3][i];
+ z = 1 - r;
+ E->rho_ref[i] = exp(tmp*z);
+ E->thermexp_ref[i] = 1;
+ E->T_ref[i] = T0 * (exp(E->control.Di * z) - 1);
}
+ for(i=1; i<=noz; i++) {
+ fprintf(stderr, "%d %f %f %f %f\n",
+ i, E->sx[1][3][i], 1-E->sx[1][3][i],
+ E->rho_ref[i], E->thermexp_ref[i]);
+ }
}
-void thermal_expansion(struct All_variables *E, double *thermexp)
+void density(struct All_variables *E, double *rho)
{
int i;
- for(i=1; i<=E->lmesh.noz; i++) {
- thermexp[i] = 1;
+ for(i=1; i<=E->lmesh.nno; i++) {
+ rho[i] = 1;
}
+
}
-
-
Modified: mc/3D/CitcomS/branches/compressible/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -1310,6 +1310,7 @@
int rank;
hsize_t *dims;
double *data;
+ float tmp;
input = h5create_group(E->hdf5.file_id, "input", (size_t)0);
@@ -1459,7 +1460,12 @@
status = set_attribute_float(input, "rayleigh", E->control.Atemp);
status = set_attribute_float(input, "dissipation", E->control.Di);
- status = set_attribute_float(input, "gruneisen", E->control.gruneisen);
+ status = set_attribute_float(input, "gruneisen",
+ ((abs(E->control.inv_gruneisen) > 1e-6)?
+ 1/E->control.inv_gruneisen :
+ E->control.inv_gruneisen
+ )
+ );
status = set_attribute_float(input, "Q0", E->control.Q0);
status = set_attribute_int(input, "stokes_flow_only", E->control.stokes);
Modified: mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -137,7 +137,8 @@
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++) {
int nz = ((i-1) % E->lmesh.noz) + 1;
- buoy[m][i] = temp * E->rho[i] * E->thermexp[nz] * E->T[m][i];
+ buoy[m][i] = temp * E->rho_ref[nz] * E->thermexp_ref[nz]
+ * (E->T[m][i] - E->T_ref[nz]);
}
phase_change_apply_410(E, buoy);
Modified: mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -266,7 +266,7 @@
}
}
- if(E->control.gruneisen > 0) {
+ if(E->control.inv_gruneisen > 0) {
for(j=1;j<=vpts;j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) * 2.0 / 3.0;
}
Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-11 21:01:06 UTC (rev 5758)
@@ -515,9 +515,12 @@
float sob_tolerance;
+ /* Rayleigh # */
float Atemp;
+ /* Dissipation # */
float Di;
- float gruneisen;
+ /* inverse of Gruneisen parameter */
+ float inv_gruneisen;
float inputdiff;
float VBXtopval;
@@ -724,8 +727,10 @@
double *BI[MAX_LEVELS][NCS],*BPI[MAX_LEVELS][NCS];
+ double *rho;
+ double *rho_ref, *thermexp_ref, *T_ref;
+
double *P[NCS],*F[NCS],*H[NCS],*S[NCS],*U[NCS];
- double *rho, *rhoref, *thermexp;
double *T[NCS],*Tdot[NCS],*buoyancy[NCS];
double *u1[NCS];
double *temp[NCS],*temp1[NCS];
Modified: mc/3D/CitcomS/branches/compressible/lib/material_properties.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/material_properties.h 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/lib/material_properties.h 2007-01-11 21:01:06 UTC (rev 5758)
@@ -28,5 +28,8 @@
void mat_prop_allocate(struct All_variables *E);
+void reference_state(struct All_variables *E);
+
+
void density(struct All_variables *E, double *rho);
Modified: mc/3D/CitcomS/branches/compressible/module/mesher.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/mesher.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/module/mesher.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -32,6 +32,7 @@
#include "global_defs.h"
#include "parallel_related.h"
+#include "material_properties.h"
void allocate_common_vars(struct All_variables*);
@@ -116,6 +117,7 @@
construct_sub_element(E);
construct_shape_functions(E);
+ reference_state(E);
mass_matrix(E);
construct_surf_det (E);
Modified: mc/3D/CitcomS/branches/compressible/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-11 16:31:31 UTC (rev 5757)
+++ mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-01-11 21:01:06 UTC (rev 5758)
@@ -434,6 +434,7 @@
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
+ float tmp;
if (!PyArg_ParseTuple(args, "OOO:Solver_set_properties",
&obj, &properties, &out))
@@ -451,7 +452,10 @@
getFloatProperty(properties, "rayleigh", E->control.Atemp, fp);
getFloatProperty(properties, "dissipation", E->control.Di, fp);
- getFloatProperty(properties, "gruneisen", E->control.gruneisen, fp);
+ getFloatProperty(properties, "gruneisen", tmp, fp);
+ if(abs(tmp) > 1e-6)
+ E->control.inv_gruneisen = 1/tmp;
+
getFloatProperty(properties, "Q0", E->control.Q0, fp);
getIntProperty(properties, "stokes_flow_only", E->control.stokes, fp);
More information about the cig-commits
mailing list