[cig-commits] r5925 - mc/3D/CitcomS/branches/compressible/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Mon Jan 29 18:07:31 PST 2007
Author: tan2
Date: 2007-01-29 18:07:30 -0800 (Mon, 29 Jan 2007)
New Revision: 5925
Modified:
mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
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.c
mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
mc/3D/CitcomS/branches/compressible/lib/global_defs.h
Log:
First cut on computing adibatic/viscous/latent heating. Most of the code is
copied from CitcomCU r2546 and refactored later. It is not tested yet.
Modified: mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Advection_diffusion.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -286,6 +286,7 @@
int bc;
unsigned int **FLAGS;
{
+ void process_heating();
void get_global_shape_fn();
void pg_shape_fn();
void element_residual();
@@ -305,36 +306,39 @@
const int ends=enodes[dims];
const int sphere_key = 0;
+ /* adiabatic and dissipative heating*/
+ process_heating(E);
+
for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++)
- DTdot[m][i] = 0.0;
+ for(i=1;i<=E->lmesh.nno;i++)
+ DTdot[m][i] = 0.0;
for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(el=1;el<=E->lmesh.nel;el++) {
+ for(el=1;el<=E->lmesh.nel;el++) {
- velo_from_element(E,VV,m,el,sphere_key);
+ velo_from_element(E,VV,m,el,sphere_key);
- get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
+ get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,E->mesh.levmax,m);
- pg_shape_fn(E,el,&PG,&GNx,VV,diff,m);
- element_residual(E,el,PG,GNx,dOmega,VV,T,Tdot,Q0,Eres,diff,E->sphere.cap[m].TB,FLAGS,m);
+ pg_shape_fn(E,el,&PG,&GNx,VV,diff,m);
+ element_residual(E,el,PG,GNx,dOmega,VV,T,Tdot,Q0,Eres,diff,E->sphere.cap[m].TB,FLAGS,m);
- for(a=1;a<=ends;a++) {
- a1 = E->ien[m][el].node[a];
- DTdot[m][a1] += Eres[a];
- }
+ for(a=1;a<=ends;a++) {
+ a1 = E->ien[m][el].node[a];
+ DTdot[m][a1] += Eres[a];
+ }
} /* next element */
(E->exchange_node_d)(E,DTdot,E->mesh.levmax);
for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++) {
- if(!(E->node[m][i] & (TBX | TBY | TBZ)))
- DTdot[m][i] *= E->Mass[m][i]; /* lumped mass matrix */
- else
- DTdot[m][i] = 0.0; /* lumped mass matrix */
- }
+ for(i=1;i<=E->lmesh.nno;i++) {
+ if(!(E->node[m][i] & (TBX | TBY | TBZ)))
+ DTdot[m][i] *= E->Mass[m][i]; /* lumped mass matrix */
+ else
+ DTdot[m][i] = 0.0; /* lumped mass matrix */
+ }
return;
}
@@ -730,3 +734,165 @@
return;
}
+
+
+static void process_visc_heating(struct All_variables *E, int m,
+ double *heating, double *total_heating)
+{
+ int e, i;
+ double visc, temp;
+ const int vpts = vpoints[E->mesh.nsd];
+
+ temp = E->control.disptn_number / E->control.Atemp;
+
+ strain_rate_2_inv(E, m, heating, 0);
+
+ *total_heating = 0;
+
+ for(e=1; e<=E->lmesh.nel; e++) {
+ visc = 0.0;
+ for(i = 1; i <= vpts; i++)
+ visc += E->EVi[m][(e-1)*vpts + i];
+ visc = visc / vpts;
+
+ heating[e] *= temp * visc;
+ *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)
+{
+ int e, ee, i, j;
+ double temp;
+ const int ends = enodes[E->mesh.nsd];
+
+ *total_heating = 0;
+
+ for(e=1; e<=E->lmesh.nel; e++) {
+ ee = (e - 1) % E->lmesh.elz + 1;
+
+ temp = 0.0;
+ for(i=1; i<=ends; i++) {
+ j = E->ien[m][e].node[i];
+ temp += E->V[m][3][j] * (E->T[m][j] + E->data.surf_temp);
+ }
+ temp = temp * E->control.disptn_number / ends;
+
+ heating[e] = temp * (E->refstate.expansivity[ee] + E->refstate.expansivity[ee + 1]) * 0.5;
+ *total_heating += heating[e] * E->eco[m][e].area;
+ }
+
+ return;
+}
+
+
+static void latent_heating(struct All_variables *E, int m,
+ double *heating_latent, double *heating_adi,
+ float **B, float Ra, float clapeyron,
+ float depth, float transT, float width)
+{
+ double temp1, temp2, temp3;
+ int e, i, j;
+ const int ends = enodes[E->mesh.nsd];
+
+ temp1 = 2.0 * width * clapeyron * Ra / E->control.Atemp;
+
+ for(e=1; e<=E->lmesh.nel; e++) {
+ temp2 = 0;
+ temp3 = 0;
+ for(i=1; i<=ends; i++) {
+ j = E->ien[m][e].node[i];
+ temp2 += temp1 * (1.0 - B[m][j]) * B[m][j]
+ * E->V[m][3][j] * (E->T[m][j] + E->data.surf_temp)
+ * E->control.disptn_number;
+ temp3 += temp1 * clapeyron * (1.0 - B[m][j])
+ * B[m][j] * (E->T[m][j] + E->data.surf_temp)
+ * E->control.disptn_number;
+ }
+ temp2 = temp2 / ends;
+ temp3 = temp3 / ends;
+ heating_adi[e] += temp2;
+ heating_latent[e] += temp3;
+ }
+ return;
+}
+
+
+static void process_latent_heating(struct All_variables *E, int m,
+ double *heating_latent, double *heating_adi)
+{
+ 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;
+ }
+
+ if(E->control.Ra_410 != 0.0) {
+ latent_heating(E, m, heating_latent, heating_adi,
+ E->Fas410, E->control.Ra_410,
+ E->control.clapeyron410, E->viscosity.z410,
+ E->control.transT410, E->control.width410);
+
+ }
+
+ if(E->control.Ra_670 != 0.0) {
+ latent_heating(E, m, heating_latent, heating_adi,
+ E->Fas670, E->control.Ra_670,
+ E->control.clapeyron670, E->viscosity.zlm,
+ E->control.transT670, E->control.width670);
+ }
+
+ if(E->control.Ra_cmb != 0.0) {
+ latent_heating(E, m, heating_latent, heating_adi,
+ E->Fascmb, E->control.Ra_cmb,
+ E->control.clapeyroncmb, E->viscosity.zcmb,
+ E->control.transTcmb, E->control.widthcmb);
+ }
+
+
+ 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 / heating_latent[e];
+ }
+
+ return;
+}
+
+
+
+void process_heating(struct All_variables *E)
+{
+ int m;
+ double temp1, temp2;
+ double total_visc_heating[NCS], total_adi_heating[NCS];
+ FILE *fp;
+ char filename[250];
+
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ const int vpts = vpoints[dims];
+
+
+ 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_latent_heating(E, m, E->heating_latent[m], E->heating_adi[m]);
+ }
+
+ /* compute total amount of visc/adi heating over all processors */
+ MPI_Allreduce(&(total_visc_heating[1]), &temp1, E->sphere.caps_per_proc,
+ MPI_DOUBLE, MPI_SUM, E->parallel.world);
+ MPI_Allreduce(&(total_adi_heating[1]), &temp2, E->sphere.caps_per_proc,
+ MPI_DOUBLE, MPI_SUM, E->parallel.world);
+
+ if(E->parallel.me == 0)
+ fprintf(E->fp, "Total_heating(adi, visc): %g %g\n", temp2, temp1);
+
+ return;
+}
Modified: mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -614,7 +614,7 @@
{
void assemble_div_u();
assemble_div_u(E, U, result, level);
- assemble_dlnrho(E, E->dlnrhodr, U, result, level);
+ assemble_dlnrho(E, E->refstate.dlnrhodr, U, result, level);
return;
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -468,6 +468,10 @@
E->mat[j] = (int *) malloc((nel+2)*sizeof(int));
E->VIP[j] = (float *) malloc((nel+2)*sizeof(float));
+ E->heating_adi[j] = (double *) malloc((nel+1)*sizeof(double));
+ E->heating_visc[j] = (double *) malloc((nel+1)*sizeof(double));
+ E->heating_latent[j] = (double *) malloc((nel+1)*sizeof(double));
+
nxyz = max(nox*noz,nox*noy);
nxyz = 2*max(nxyz,noz*noy);
@@ -594,6 +598,10 @@
for(i=1;i<=E->lmesh.nel;i++) {
E->mat[j][i]=1;
E->VIP[j][i]=1.0;
+
+ E->heating_adi[j][i] = 0;
+ E->heating_visc[j][i] = 0;
+ E->heating_latent[j][i] = 1.0; //TODO: why 1?
}
for(i=1;i<=E->lmesh.npno;i++)
Modified: mc/3D/CitcomS/branches/compressible/lib/Material_properties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -43,16 +43,16 @@
int nel = E->lmesh.nel;
/* reference profile of density */
- E->rho_ref = (double *) malloc((noz+1)*sizeof(double));
+ E->refstate.rho = (double *) malloc((noz+1)*sizeof(double));
/* reference profile of coefficient of thermal expansion */
- E->thermexp_ref = (double *) malloc((noz+1)*sizeof(double));
+ E->refstate.expansivity = (double *) malloc((noz+1)*sizeof(double));
/* reference profile of temperature */
- E->T_ref = (double *) malloc((noz+1)*sizeof(double));
+ E->refstate.Tadi = (double *) malloc((noz+1)*sizeof(double));
/* reference profile of d(ln(rho_ref))/dr */
- E->dlnrhodr = (double *) malloc((nel+1)*sizeof(double));
+ E->refstate.dlnrhodr = (double *) malloc((nel+1)*sizeof(double));
}
@@ -65,25 +65,29 @@
tmp = 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=%f, gamma=%f, surf_temp=%f, dT=%f\n",
+ E->control.disptn_number, 1/E->control.inv_gruneisen,
+ 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.disptn_number * z) - 1);
+ E->refstate.rho[i] = exp(tmp*z);
+ E->refstate.expansivity[i] = 1;
+ E->refstate.Tadi[i] = T0 * (exp(E->control.disptn_number * z) - 1);
}
for(i=1; i<=nel; i++) {
// TODO: dln(rho)/dr
- E->dlnrhodr[i] = - tmp;
+ E->refstate.dlnrhodr[i] = - tmp;
}
if(E->parallel.me < E->parallel.nprocz)
for(i=1; i<=noz; i++) {
fprintf(stderr, "%d %f %f %f %f\n",
i+E->lmesh.nzs-1, E->sx[1][3][i], 1-E->sx[1][3][i],
- E->rho_ref[i], E->thermexp_ref[i]);
+ E->refstate.rho[i], E->refstate.Tadi[i]);
}
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Output.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Output.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -47,6 +47,7 @@
void output_horiz_avg(struct All_variables *, int);
void output_tracer(struct All_variables *, int);
void output_pressure(struct All_variables *, int);
+void output_heating(struct All_variables *, int);
extern void parallel_process_termination();
extern void heat_flux(struct All_variables *);
@@ -83,6 +84,9 @@
if(E->control.tracer==1)
output_tracer(E, cycles);
+ if(E->control.disptn_number != 0)
+ output_heating(E, cycles);
+
/* optiotnal output below */
/* compute and output geoid (in spherical harmonics coeff) */
if (E->output.geoid == 1)
@@ -406,6 +410,30 @@
}
+void output_heating(struct All_variables *E, int cycles)
+{
+ int j, e;
+ char output_file[255];
+ FILE *fp1;
+
+ sprintf(output_file,"%s.heating.%d.%d", E->control.data_file,
+ E->parallel.me, cycles);
+ fp1 = output_open(output_file);
+
+ fprintf(fp1,"%.5e\n",E->monitor.elapsed_time);
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++) {
+ fprintf(fp1,"%3d %7d\n", j, E->lmesh.nel);
+ for(e=1; e<=E->lmesh.nel; e++)
+ fprintf(fp1, "%.4e %.4e %.4e\n", E->heating_adi[j][e],
+ E->heating_visc[j][e], E->heating_latent[j][e]);
+ }
+ fclose(fp1);
+
+ return;
+}
+
+
void output_time(struct All_variables *E, int cycles)
{
double CPU_time0();
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-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -140,8 +140,8 @@
/* We don't need to substract adiabatic T profile from T here,
* since the horizontal average of buoy will be removed.
*/
- buoy[m][i] = temp * E->rho_ref[nz] * E->thermexp_ref[nz]
- * E->T[m][i];
+ buoy[m][i] = temp * E->refstate.rho[nz]
+ * E->refstate.expansivity[nz] * E->T[m][i];
}
phase_change_apply_410(E, buoy);
Modified: mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -393,8 +393,8 @@
r0dotrt = alpha = omega = 0;
while( (valid) && (count < *steps_max) &&
- ((E->monitor.incompressibility >= E->control.tole_comp) ||
- (dpressure >= imp) || (dvelocity >= imp)) ) {
+ ((E->monitor.incompressibility >= E->control.tole_comp) &&
+ (dpressure >= imp) && (dvelocity >= imp)) ) {
Modified: mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -219,6 +219,16 @@
velo_from_element(E,VV,m,e,sphere_key);
+ /* Vxyz is the strain rate vector, whose relationship with
+ * the strain rate tensor (e) is that:
+ * Vxyz[1] = e11
+ * Vxyz[2] = e22
+ * Vxyz[3] = e33
+ * Vxyz[4] = 2*e12
+ * Vxyz[5] = 2*e13
+ * Vxyz[6] = 2*e23
+ * where 1 is theta, 2 is phi, and 3 is r
+ */
for(j=1;j<=vpts;j++) {
pre[j] = E->EVi[m][(e-1)*vpts+j]*dOmega.vpt[j];
dilation[j] = 0.0;
@@ -268,7 +278,7 @@
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;
+ dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
}
for(j=1;j<=vpts;j++) {
@@ -291,9 +301,10 @@
div /= E->eco[m][e].area;
vor /= E->eco[m][e].area;
- Szz -= E->P[m][e]; /* add the pressure term */
- Sxx -= E->P[m][e]; /* add the pressure term */
- Syy -= E->P[m][e]; /* add the pressure term */
+ /* add the pressure term */
+ Szz -= E->P[m][e];
+ Sxx -= E->P[m][e];
+ Syy -= E->P[m][e];
for(i=1;i<=ends;i++) {
node = E->ien[m][e].node[i];
Modified: mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-01-30 02:07:30 UTC (rev 5925)
@@ -481,48 +481,104 @@
struct Shape_function_dA dOmega;
struct Shape_function_dx GNx;
- double edot[4][4],dudx[4][4],rtf[4][9];
- float VV[4][9];
- int e,i,p,q,n,nel,k;
+ double edot[4][4], rtf[4][9];
+ float VV[4][9], Vxyz[7][9], dilation[9];
+
+ int e, i, j, p, q, n;
+
+ const int nel = E->lmesh.nel;
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
const int lev = E->mesh.levmax;
- const int nno = E->lmesh.nno;
- const int vpts = vpoints[dims];
- const int sphere_key = 0;
+ const int ppts = ppoints[dims];
+ const int sphere_key = 1;
- nel = E->lmesh.nel;
- for(e=1;e<=nel;e++) {
+ for(e=1; e<=nel; e++) {
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
+ /* get shape function on presure gauss points */
+ get_global_shape_fn(E, e, &GN, &GNx, &dOmega, 2,
+ sphere_key, rtf, lev, m);
- velo_from_element(E,VV,m,e,sphere_key);
+ velo_from_element(E, VV, m, e, sphere_key);
- for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- dudx[p][q] = 0.0;
- for(i=1;i<=ends;i++)
- for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- dudx[p][q] += VV[p][i] * GNx.ppt[GNPXINDEX(q-1,i,1)];
+ /* Vxyz is the strain rate vector, whose relationship with
+ * the strain rate tensor (e) is that:
+ * Vxyz[1] = e11
+ * Vxyz[2] = e22
+ * Vxyz[3] = e33
+ * Vxyz[4] = 2*e12
+ * Vxyz[5] = 2*e13
+ * Vxyz[6] = 2*e23
+ * where 1 is theta, 2 is phi, and 3 is r
+ */
+ for(j=1; j<=ppts; j++) {
+ Vxyz[1][j] = 0.0;
+ Vxyz[2][j] = 0.0;
+ Vxyz[3][j] = 0.0;
+ Vxyz[4][j] = 0.0;
+ Vxyz[5][j] = 0.0;
+ Vxyz[6][j] = 0.0;
+ }
- for(p=1;p<=dims;p++)
- for(q=1;q<=dims;q++)
- edot[p][q] = dudx[p][q] + dudx[q][p];
+ for(j=1; j<=ppts; j++) {
+ for(i=1; i<=ends; i++) {
+ Vxyz[1][j] += (VV[1][i] * GNx.ppt[GNPXINDEX(0, i, j)]
+ + VV[3][i] * E->N.ppt[GNPINDEX(i, j)])
+ * rtf[3][j];
+ Vxyz[2][j] += ((VV[2][i] * GNx.ppt[GNPXINDEX(1, i, j)]
+ + VV[1][i] * E->N.ppt[GNPINDEX(i, j)]
+ * cos(rtf[1][j])) / sin(rtf[1][j])
+ + VV[3][i] * E->N.ppt[GNPINDEX(i, j)])
+ * rtf[3][j];
+ Vxyz[3][j] += VV[3][i] * GNx.ppt[GNPXINDEX(2, i, j)];
- if (dims==2)
- EEDOT[e] = edot[1][1]*edot[1][1] + edot[2][2]*edot[2][2]
- + edot[1][2]*edot[1][2]*2.0;
+ Vxyz[4][j] += ((VV[1][i] * GNx.ppt[GNPXINDEX(1, i, j)]
+ - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]
+ * cos(rtf[1][j])) / sin(rtf[1][j])
+ + VV[2][i] * GNx.ppt[GNPXINDEX(0, i, j)])
+ * rtf[3][j];
+ Vxyz[5][j] += VV[1][i] * GNx.ppt[GNPXINDEX(2, i, j)]
+ + rtf[3][j] * (VV[3][i] * GNx.ppt[GNPXINDEX(0, i, j)]
+ - VV[1][i] * E->N.ppt[GNPINDEX(i, j)]);
+ Vxyz[6][j] += VV[2][i] * GNx.ppt[GNPXINDEX(2, i, j)]
+ + rtf[3][j] * (VV[3][i]
+ * GNx.ppt[GNPXINDEX(1, i, j)]
+ / sin(rtf[1][j])
+ - VV[2][i] * E->N.ppt[GNPINDEX(i, j)]);
+ }
+ }
- else if (dims==3)
- EEDOT[e] = edot[1][1]*edot[1][1] + edot[1][2]*edot[1][2]*2.0
- + edot[2][2]*edot[2][2] + edot[2][3]*edot[2][3]*2.0
- + edot[3][3]*edot[3][3] + edot[1][3]*edot[1][3]*2.0;
+ if(E->control.inv_gruneisen > 0)
+ for(j=1; j<=ppts; j++)
+ dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
+ else
+ for(j=1; j<=ppts; j++)
+ dilation[j] = 0;
+ edot[1][1] = edot[2][2] = edot[3][3] = 0;
+ edot[1][2] = edot[1][3] = edot[2][3] = 0;
+
+ /* edot is 2 * (the deviatoric strain rate tensor) */
+ for(j=1; j<=ppts; j++) {
+ edot[1][1] += 2.0 * (Vxyz[1][j] - dilation[j]);
+ edot[2][2] += 2.0 * (Vxyz[2][j] - dilation[j]);
+ edot[3][3] += 2.0 * (Vxyz[3][j] - dilation[j]);
+ edot[1][2] += Vxyz[4][j];
+ edot[1][3] += Vxyz[5][j];
+ edot[2][3] += Vxyz[6][j];
+ }
+
+ /* TODO: figure out why 2.0 factor is here */
+ EEDOT[e] = edot[1][1] * edot[1][1]
+ + edot[1][2] * edot[1][2] * 2.0
+ + edot[2][2] * edot[2][2]
+ + edot[2][3] * edot[2][3] * 2.0
+ + edot[3][3] * edot[3][3]
+ + edot[1][3] * edot[1][3] * 2.0;
}
if(SQRT)
Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-29 20:51:24 UTC (rev 5924)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-01-30 02:07:30 UTC (rev 5925)
@@ -539,6 +539,15 @@
};
+
+struct REF_STATE {
+ double *rho;
+ double *expansivity;
+ double *Tadi;
+ double *dlnrhodr;
+};
+
+
struct DATA {
float layer_km;
float radius_km;
@@ -639,8 +648,6 @@
struct SBC sbc;
struct Output output;
- int filed[20];
-
struct COORD *eco[NCS];
struct IEN *ien[NCS]; /* global */
struct SIEN *sien[NCS];
@@ -662,6 +669,9 @@
struct CC element_Cc;
struct CCX element_Ccx;
+ struct REF_STATE refstate;
+
+
higher_precision *Eqn_k1[MAX_LEVELS][NCS],*Eqn_k2[MAX_LEVELS][NCS],*Eqn_k3[MAX_LEVELS][NCS];
int *Node_map [MAX_LEVELS][NCS];
int *Node_eqn [MAX_LEVELS][NCS];
@@ -670,8 +680,9 @@
double *BI[MAX_LEVELS][NCS],*BPI[MAX_LEVELS][NCS];
double *rho;
- double *rho_ref, *thermexp_ref, *T_ref;
- double *dlnrhodr;
+ double *heating_adi[NCS];
+ double *heating_visc[NCS];
+ double *heating_latent[NCS];
double *P[NCS],*F[NCS],*H[NCS],*S[NCS],*U[NCS];
double *T[NCS],*Tdot[NCS],*buoyancy[NCS];
More information about the cig-commits
mailing list