[cig-commits] r7898 - in mc/3D/CitcomS/trunk: CitcomS/Components
CitcomS/Solver lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Aug 27 18:09:15 PDT 2007
Author: tan2
Date: 2007-08-27 18:09:14 -0700 (Mon, 27 Aug 2007)
New Revision: 7898
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Const.py
mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Material_properties.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Added various heating term to the energy equation.
A new input parameter: solver.surfaceT for non-dimensional surface temperature.
Removed old (and unused) input parameter: solver.const.surf_temp.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Const.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Const.py 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Const.py 2007-08-28 01:09:14 UTC (rev 7898)
@@ -81,7 +81,6 @@
cp = pyre.inventory.float("cp", default=1200.0)
density_above = pyre.inventory.float("density_above", default=1030.0)
density_below = pyre.inventory.float("density_below", default=6600.0)
- surftemp = pyre.inventory.float("surftemp", default=273.0)
z_lith = pyre.inventory.float("z_lith", default=0.014)
z_410 = pyre.inventory.float("z_410", default=0.06435)
Modified: mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2007-08-28 01:09:14 UTC (rev 7898)
@@ -306,6 +306,7 @@
rayleigh = inv.float("rayleigh", default=1e+05)
dissipation_number = inv.float("dissipation_number", default=0.0)
gruneisen = inv.float("gruneisen", default=0.0)
+ surfaceT = inv.float("surfaceT", default=0.1)
Q0 = inv.float("Q0", default=0.0)
stokes_flow_only = inv.bool("stokes_flow_only", default=False)
Modified: mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-28 01:09:14 UTC (rev 7898)
@@ -540,7 +540,8 @@
double prod,sfn;
struct Shape_function1 GM;
struct Shape_function1_dA dGamma;
- double temp;
+ double temp,rho,heating;
+ int nz;
void get_global_1d_shape_fn();
@@ -600,7 +601,15 @@
Q += E->composition.comp_el[m][0][el] * E->control.Q0ER;
}
+ nz = ((el-1) % E->lmesh.noz) + 1;
+ rho = 0.5 * (E->refstate.rho[nz] + E->refstate.rho[nz+1]);
+ if(E->control.disptn_number == 0)
+ heating = rho * Q;
+ else
+ heating = (rho * Q - E->heating_adi[m][el] + E->heating_visc[m][el])
+ * E->heating_latent[m][el];
+
/* construct residual from this information */
@@ -610,7 +619,7 @@
for(i=1;i<=vpts;i++)
Eres[j] -=
PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i]
- * (dT[i] - Q + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])
+ * (dT[i] - heating + v1[i]*tx1[i] + v2[i]*tx2[i] + v3[i]*tx3[i])
+ diff*dOmega.vpt[i] * (GNx.vpt[GNVXINDEX(0,j,i)]*tx1[i]*rtf[3][i] +
GNx.vpt[GNVXINDEX(1,j,i)]*tx2[i]*sint[i] +
GNx.vpt[GNVXINDEX(2,j,i)]*tx3[i] );
@@ -621,7 +630,7 @@
for(j=1;j<=ends;j++) {
Eres[j]=0.0;
for(i=1;i<=vpts;i++)
- Eres[j] -= PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i] * (dT[i] - Q + v1[i] * tx1[i] + v2[i] * tx2[i] + v3[i] * tx3[i]);
+ Eres[j] -= PG.vpt[GNVINDEX(j,i)] * dOmega.vpt[i] * (dT[i] - heating + v1[i] * tx1[i] + v2[i] * tx2[i] + v3[i] * tx3[i]);
}
}
@@ -731,7 +740,7 @@
static void process_visc_heating(struct All_variables *E, int m,
- double *heating, double *total_heating)
+ double *heating)
{
void strain_rate_2_inv();
int e, i;
@@ -754,18 +763,12 @@
free(strain_sqr);
-
- /* sum up */
- *total_heating = 0;
- for(e=1; e<=E->lmesh.nel; e++)
- *total_heating += heating[e] * E->eco[m][e].area;
-
return;
}
static void process_adi_heating(struct All_variables *E, int m,
- double *heating, double *total_heating)
+ double *heating)
{
int e, ez, i, j;
double temp, temp2;
@@ -779,20 +782,15 @@
for(i=1; i<=ends; i++) {
j = E->ien[m][e].node[i];
temp += E->sphere.cap[m].V[3][j]
- * (E->T[m][j] + E->data.surf_temp);
+ * (E->T[m][j] + E->control.surface_temp);
}
- /* XXX: missing gravity */
- heating[e] = temp * temp2 * 0.25
+ heating[e] = temp * temp2 * 0.125
* (E->refstate.thermal_expansivity[ez] + E->refstate.thermal_expansivity[ez + 1])
- * (E->refstate.rho[ez] + E->refstate.rho[ez + 1]);
+ * (E->refstate.rho[ez] + E->refstate.rho[ez + 1])
+ * (E->refstate.gravity[ez] + E->refstate.gravity[ez + 1]);
}
- /* sum up */
- *total_heating = 0;
- for(e=1; e<=E->lmesh.nel; e++)
- *total_heating += heating[e] * E->eco[m][e].area;
-
return;
}
@@ -814,10 +812,10 @@
for(i=1; i<=ends; i++) {
j = E->ien[m][e].node[i];
temp2 += (1.0 - B[m][j]) * B[m][j]
- * E->sphere.cap[m].V[3][j] * (E->T[m][j] + E->data.surf_temp)
+ * E->sphere.cap[m].V[3][j] * (E->T[m][j] + E->control.surface_temp)
* E->control.disptn_number;
temp3 += clapeyron * (1.0 - B[m][j])
- * B[m][j] * (E->T[m][j] + E->data.surf_temp)
+ * B[m][j] * (E->T[m][j] + E->control.surface_temp)
* E->control.disptn_number;
}
heating_adi[e] += temp2 * temp1;
@@ -832,11 +830,9 @@
{
int e;
- if(E->control.Ra_410 != 0 || E->control.Ra_670 != 0.0 ||
- E->control.Ra_cmb != 0) {
- for(e=1; e<=E->lmesh.nel; e++)
- heating_latent[e] = 1.0;
- }
+ /* reset */
+ for(e=1; e<=E->lmesh.nel; e++)
+ heating_latent[e] = 1.0;
if(E->control.Ra_410 != 0.0) {
latent_heating(E, m, heating_latent, heating_adi,
@@ -871,36 +867,48 @@
}
+static double total_heating(struct All_variables *E, double **heating)
+{
+ int m, e;
+ double sum, total;
+ /* sum up within each processor */
+ sum = 0;
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ for(e=1; e<=E->lmesh.nel; e++)
+ sum += heating[m][e] * E->eco[m][e].area;
+ }
+
+ /* sum up for all processors */
+ MPI_Allreduce(&sum, &total, 1,
+ MPI_DOUBLE, MPI_SUM, E->parallel.world);
+
+ return total;
+}
+
+
static void process_heating(struct All_variables *E)
{
int m;
- double total_visc_heating[NCS], total_adi_heating[NCS];
- double temp1, temp2, temp3, temp4;
+ double total_visc_heating, total_adi_heating;
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- process_visc_heating(E, m, E->heating_visc[m], &(total_visc_heating[m]));
- process_adi_heating(E, m, E->heating_adi[m], &(total_adi_heating[m]));
+ process_visc_heating(E, m, E->heating_visc[m]);
+ process_adi_heating(E, m, E->heating_adi[m]);
process_latent_heating(E, m, E->heating_latent[m], E->heating_adi[m]);
}
/* compute total amount of visc/adi heating over all processors */
- temp3 = temp4 = 0;
- for(m=1; m<=E->sphere.caps_per_proc; m++) {
- MPI_Allreduce(&(total_visc_heating[m]), &temp1,
- E->sphere.caps_per_proc,
- MPI_DOUBLE, MPI_SUM, E->parallel.world);
- MPI_Allreduce(&(total_adi_heating[m]), &temp2,
- E->sphere.caps_per_proc,
- MPI_DOUBLE, MPI_SUM, E->parallel.world);
+ total_visc_heating = total_heating(E, E->heating_visc);
+ total_adi_heating = total_heating(E, E->heating_adi);
- temp3 += temp1;
- temp4 += temp2;
- }
-
if(E->parallel.me == 0)
fprintf(E->fp, "Step: %d, Total_heating(visc, adi): %g %g\n",
- E->monitor.solution_cycles, temp3, temp4);
+ E->monitor.solution_cycles,
+ total_visc_heating, total_adi_heating);
+ fprintf(stderr, "Step: %d, Total_heating(visc, adi): %g %g\n",
+ E->monitor.solution_cycles,
+ total_visc_heating, total_adi_heating);
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-28 01:09:14 UTC (rev 7898)
@@ -438,11 +438,11 @@
myerror(E, "Error: unknown Uzawa iteration\n");
}
-
- /* data section */
+ input_float("surfaceT",&(E->control.surface_temp),"0.1",m);
input_float("Q0",&(E->control.Q0),"0.0",m);
/* Q0_enriched gets read in Tracer_setup.c */
+ /* data section */
input_float("gravacc",&(E->data.grav_acc),"9.81",m);
input_float("thermexp",&(E->data.therm_exp),"3.0e-5",m);
input_float("cp",&(E->data.Cp),"1200.0",m);
@@ -451,7 +451,6 @@
input_float("density_above",&(E->data.density_above),"1030.0",m);
input_float("density_below",&(E->data.density_below),"6600.0",m);
input_float("refvisc",&(E->data.ref_viscosity),"1.0e21",m);
- input_float("surftemp",&(E->data.surf_temp),"273.0",m);
input_float("radius",&tmp,"6371e3.0",m);
E->data.radius_km = tmp / 1e3;
@@ -674,7 +673,7 @@
E->heating_adi[j][i] = 0;
E->heating_visc[j][i] = 0;
- E->heating_latent[j][i] = 1.0; //TODO: why 1?
+ E->heating_latent[j][i] = 1.0;
}
for(i=1;i<=E->lmesh.npno;i++)
@@ -839,7 +838,6 @@
E->data.gas_const = 8.3;
E->data.surf_heat_flux = 4.4e-2;
E->data.grav_const = 6.673e-11;
- E->data.surf_temp = 0.0;
E->data.youngs_mod = 1.0e11;
E->data.Te = 0.0;
E->data.T_sol0 = 1373.0; /* Dave's values 1991 (for the earth) */
Modified: mc/3D/CitcomS/trunk/lib/Material_properties.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Material_properties.c 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/lib/Material_properties.c 2007-08-28 01:09:14 UTC (rev 7898)
@@ -71,11 +71,7 @@
double r, z, beta, T0;
beta = E->control.disptn_number * E->control.inv_gruneisen;
- T0 = E->data.surf_temp / E->data.ref_temperature;
- if(E->parallel.me == 0)
- fprintf(stderr, "Di=%e, gamma=%e, surf_temp=%e, dT=%e\n",
- E->control.disptn_number, 1.0/E->control.inv_gruneisen,
- E->data.surf_temp, E->data.ref_temperature);
+ T0 = E->control.surface_temp / E->data.ref_temperature;
/* All refstate variables (except Tadi) must be 1 at the surface.
* Otherwise, the scaling of eqns in the code might not be correct. */
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-28 01:09:14 UTC (rev 7898)
@@ -1368,7 +1368,6 @@
status = set_attribute_float(input, "cp", E->data.Cp);
status = set_attribute_float(input, "density_above", E->data.density_above);
status = set_attribute_float(input, "density_below", E->data.density_below);
- status = set_attribute_float(input, "surftemp", E->data.surf_temp);
status = set_attribute_float(input, "z_lith", E->viscosity.zlith);
status = set_attribute_float(input, "z_410", E->viscosity.z410);
@@ -1464,6 +1463,7 @@
(E->control.inv_gruneisen == 0)?
1.0/E->control.inv_gruneisen :
E->control.inv_gruneisen);
+ status = set_attribute_float(input, "surfaceT", E->control.surface_temp);
status = set_attribute_float(input, "Q0", E->control.Q0);
status = set_attribute_int(input, "stokes_flow_only", E->control.stokes);
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-28 01:09:14 UTC (rev 7898)
@@ -455,6 +455,9 @@
/* inverse of Gruneisen parameter */
float inv_gruneisen;
+ /* surface temperature */
+ float surface_temp;
+
/**/
float relative_err_accuracy;
@@ -552,7 +555,6 @@
float melt_viscosity;
float permeability;
float grav_const;
- float surf_temp;
float youngs_mod;
float Te;
float ref_temperature;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-27 20:09:06 UTC (rev 7897)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-28 01:09:14 UTC (rev 7898)
@@ -197,7 +197,6 @@
getFloatProperty(properties, "cp", E->data.Cp, fp);
getFloatProperty(properties, "density_above", E->data.density_above, fp);
getFloatProperty(properties, "density_below", E->data.density_below, fp);
- getFloatProperty(properties, "surftemp", E->data.surf_temp, fp);
E->data.therm_cond = E->data.therm_diff * E->data.density * E->data.Cp;
E->data.ref_temperature = E->control.Atemp * E->data.therm_diff
@@ -467,6 +466,7 @@
else
E->control.inv_gruneisen = 0;
+ getFloatProperty(properties, "surfaceT", E->control.surface_temp, fp);
getFloatProperty(properties, "Q0", E->control.Q0, fp);
getIntProperty(properties, "stokes_flow_only", E->control.stokes, fp);
More information about the cig-commits
mailing list