[cig-commits] r7831 - in mc/3D/CitcomS/trunk: CitcomS/Solver lib
module
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Aug 15 17:27:44 PDT 2007
Author: tan2
Date: 2007-08-15 17:27:43 -0700 (Wed, 15 Aug 2007)
New Revision: 7831
Modified:
mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
mc/3D/CitcomS/trunk/lib/Advection_diffusion.c
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Convection.c
mc/3D/CitcomS/trunk/lib/Drive_solvers.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Makefile.am
mc/3D/CitcomS/trunk/lib/Obsolete.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
mc/3D/CitcomS/trunk/lib/Size_does_matter.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/Topo_gravity.c
mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/mesher.c
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
Merging compressible branch (up to r6584) to trunk.
Modified: mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/CitcomS/Solver/Solver.py 2007-08-16 00:27:43 UTC (rev 7831)
@@ -290,6 +290,8 @@
datadir_old = inv.str("datadir_old", default="")
rayleigh = inv.float("rayleigh", default=1e+05)
+ dissipation_number = inv.float("dissipation_number", default=0.0)
+ gruneisen = inv.float("gruneisen", default=0.0)
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-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Advection_diffusion.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -127,6 +127,7 @@
predictor(E,E->T,E->Tdot);
for(psc_pass=0;psc_pass<E->advection.temp_iterations;psc_pass++) {
+ /* XXX: replace inputdiff with refstate.conductivity */
pg_solver(E,E->T,E->Tdot,DTdot,E->convection.heat_sources,E->control.inputdiff,1,E->node);
corrector(E,E->T,E->Tdot,DTdot);
temperatures_conform_bcs(E);
@@ -151,7 +152,7 @@
if(E->control.filter_temperature)
filter(E);
-
+
/* time0= CPU_time0()-time1; */
/* if(E->control.verbose) */
/* fprintf(E->fp_out,"time=%f\n",time0); */
@@ -294,6 +295,7 @@
int bc;
unsigned int **FLAGS;
{
+ void process_heating();
void get_global_shape_fn();
void pg_shape_fn();
void element_residual();
@@ -313,6 +315,9 @@
const int ends=enodes[dims];
const int sphere_key = 1;
+ /* 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;
@@ -325,6 +330,7 @@
get_global_shape_fn(E, el, &GN, &GNx, &dOmega, 0,
sphere_key, rtf, E->mesh.levmax, m);
+ /* XXX: replace diff with refstate.conductivity */
pg_shape_fn(E, el, &PG, &GNx, VV,
rtf, diff, m);
element_residual(E, el, PG, GNx, dOmega, VV, T, Tdot,
@@ -343,7 +349,7 @@
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 */
+ DTdot[m][i] *= E->TMass[m][i]; /* lumped mass matrix */
else
DTdot[m][i] = 0.0; /* lumped mass matrix */
}
@@ -709,3 +715,178 @@
return;
}
+
+
+static void process_visc_heating(struct All_variables *E, int m,
+ double *heating, double *total_heating)
+{
+ int e, i;
+ double visc, temp;
+ float *strain_sqr;
+ const int vpts = vpoints[E->mesh.nsd];
+
+ strain_sqr = (float*) malloc((E->lmesh.nel+1)*sizeof(float));
+ temp = E->control.disptn_number / E->control.Atemp / vpts;
+
+ strain_rate_2_inv(E, m, strain_sqr, 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];
+
+ heating[e] = temp * visc * strain_sqr[e];
+ }
+
+ 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)
+{
+ int e, ez, i, j;
+ double temp, temp2;
+ const int ends = enodes[E->mesh.nsd];
+
+ temp2 = E->control.disptn_number / ends;
+ for(e=1; e<=E->lmesh.nel; e++) {
+ ez = (e - 1) % E->lmesh.elz + 1;
+
+ temp = 0.0;
+ 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);
+ }
+
+ /* XXX: missing gravity */
+ heating[e] = temp * temp2 * 0.25
+ * (E->refstate.expansivity[ez] + E->refstate.expansivity[ez + 1])
+ * (E->refstate.rho[ez] + E->refstate.rho[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;
+}
+
+
+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->sphere.cap[m].V[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/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -1,6 +1,6 @@
/*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- *
+ *
*<LicenseText>
*
* CitcomS by Louis Moresi, Shijie Zhong, Lijie Han, Eh Tan,
@@ -22,7 +22,7 @@
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*</LicenseText>
- *
+ *
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
#include <math.h>
@@ -343,7 +343,7 @@
for(element=1;element<=nel;element++) {
- get_elt_k(E,element,elt_K,level,m);
+ get_elt_k(E,element,elt_K,level,m,0);
if (E->control.augmented_Lagr)
get_aug_k(E,element,elt_K,level,m);
@@ -635,7 +635,7 @@
for(el=1;el<=E->lmesh.NEL[lev];el++) {
- get_elt_k(E,el,E->elt_k[lev][m][el].k,lev,m); /* not for penalty */
+ get_elt_k(E,el,E->elt_k[lev][m][el].k,lev,m,0);
if (E->control.augmented_Lagr)
get_aug_k(E,el,E->elt_k[lev][m][el].k,lev,m);
@@ -696,14 +696,11 @@
void construct_node_maps();
void construct_node_ks();
void construct_elt_ks();
- void construct_elt_gs();
void rebuild_BI_on_boundary();
if (E->control.NMULTIGRID)
project_viscosity(E);
- construct_elt_gs(E);
-
if (E->control.NMULTIGRID || E->control.NASSEMBLE) {
construct_node_ks(E);
}
@@ -759,18 +756,18 @@
float r;
{
float zlith, z410, zlm;
-
+
int llayers = 0;
zlith = E->viscosity.zlith;
z410 = E->viscosity.z410;
zlm = E->viscosity.zlm;
-
+
if (r > (E->sphere.ro-zlith))
llayers = 1;
- else if ((r > (E->sphere.ro-z410))&&
+ else if ((r > (E->sphere.ro-z410))&&
(r <= (E->sphere.ro-zlith)))
llayers = 2;
- else if ((r > (E->sphere.ro-zlm)) &&
+ else if ((r > (E->sphere.ro-zlm)) &&
(r <= (E->sphere.ro-z410)))
llayers = 3;
else
Modified: mc/3D/CitcomS/trunk/lib/Convection.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Convection.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Convection.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -48,11 +48,11 @@
void convection_initial_fields();
void twiddle_thumbs();
- E->advection.timestep = 0.0;
- E->advection.timesteps = 0;
+ E->advection.temp_iterations = 2; /* petrov-galerkin iterations: minimum value. */
E->advection.total_timesteps = 1;
E->advection.sub_iterations = 1;
E->advection.last_sub_iterations = 1;
+ E->advection.gamma = 0.5;
E->advection.dt_reduced = 1.0;
E->monitor.T_maxvaried = 1.05;
Modified: mc/3D/CitcomS/trunk/lib/Drive_solvers.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Drive_solvers.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -79,8 +79,6 @@
const int dims = E->mesh.nsd;
const int addi_dof = additional_dof[dims];
- E->monitor.elapsed_time_vsoln = E->monitor.elapsed_time;
-
velocities_conform_bcs(E,E->U);
assemble_forces(E,0);
@@ -169,8 +167,6 @@
const int dims = E->mesh.nsd;
const int addi_dof = additional_dof[dims];
- E->monitor.elapsed_time_vsoln = E->monitor.elapsed_time;
-
velocities_conform_bcs(E,E->U);
E->monitor.stop_topo_loop = 0;
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -32,6 +32,7 @@
#include <math.h>
#include "element_definitions.h"
#include "global_defs.h"
+#include "material_properties.h"
@@ -152,21 +153,90 @@
return;
}
+
/*==============================================================
+ Function to supply the element strain-displacement matrix Ba,
+ which is used to compute element stiffness matrix
+ ============================================================== */
+
+void get_ba(struct Shape_function *N, struct Shape_function_dx *GNx,
+ struct CC *cc, struct CCX *ccx, double rtf[4][9],
+ int dims, double ba[9][9][4][7])
+{
+ int k, a, n;
+ const int vpts = VPOINTS3D;
+ const int ends = ENODES3D;
+
+ double ra[9], si[9], ct[9];
+
+ for(k=1;k<=vpts;k++) {
+ ra[k] = rtf[3][k];
+ si[k] = sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k])/si[k];
+ }
+
+ for(n=1;n<=dims;n++)
+ for(k=1;k<=vpts;k++)
+ for(a=1;a<=ends;a++) {
+ ba[a][k][n][1]=
+ (GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(1,n,1,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ )*ra[k];
+
+ ba[a][k][n][2]=
+ (N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(1,n,a,k)]*ct[k]
+ + N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(2,n,2,a,k)]
+ )/si[k]
+ )*ra[k];
+
+ ba[a][k][n][3]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(3,n,a,k)];
+
+ ba[a][k][n][4]=
+ (GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(2,n,1,a,k)]
+ - N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(2,n,a,k)]*ct[k]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(1,n,2,a,k)]
+ )/si[k]
+ )*ra[k];
+
+ ba[a][k][n][5]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(1,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(0,a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*(ccx->vpt[BVXINDEX(3,n,1,a,k)]
+ - cc->vpt[BVINDEX(1,n,a,k)])
+ )*ra[k];
+
+ ba[a][k][n][6]=
+ GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ -ra[k]*N->vpt[GNVINDEX(a,k)]*cc->vpt[BVINDEX(2,n,a,k)]
+ +(GNx->vpt[GNVXINDEX(1,a,k)]*cc->vpt[BVINDEX(3,n,a,k)]
+ + N->vpt[GNVINDEX(a,k)]*ccx->vpt[BVXINDEX(3,n,2,a,k)]
+ )/si[k]*ra[k];
+ }
+
+ return;
+}
+
+/*==============================================================
Function to supply the element k matrix for a given element e.
============================================================== */
-void get_elt_k(E,el,elt_k,lev,m)
+void get_elt_k(E,el,elt_k,lev,m,iconv)
struct All_variables *E;
int el,m;
double elt_k[24*24];
- int lev;
+ int lev, iconv;
{
double bdbmu[4][4];
int pn,qn,ad,bd;
- int a,b,i,j,j1,k,n;
- double rtf[4][9],W[9],ra[9],si[9],ct[9];
+ int a,b,i,j,i1,j1,k;
+ double rtf[4][9],W[9];
void get_global_shape_fn();
void construct_c3x3matrix_el();
struct Shape_function GN;
@@ -174,16 +244,34 @@
struct Shape_function_dx GNx;
double ba[9][9][4][7]; /* integration points,node,3x6 matrix */
+ /* d1 and d2 are the elastic coefficient matrix for incompressible
+ and compressible creeping flow, respectively */
+ const double d1[7][7] =
+ { {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} };
+ const double d2[7][7] =
+ { {0., 0., 0., 0., 0., 0., 0.},
+ {0., 2.-2./3., -2./3., -2./3., 0., 0., 0.},
+ {0., -2./3., 2.-2./3., -2./3., 0., 0., 0.},
+ {0., -2./3., -2./3., 2.-2./3., 0., 0., 0.},
+ {0., 0., 0., 0., 1., 0., 0.},
+ {0., 0., 0., 0., 0., 1., 0.},
+ {0., 0., 0., 0., 0., 0., 1.} };
const int nn=loc_mat_size[E->mesh.nsd];
- const int vpts=vpoints[E->mesh.nsd];
- const int ends=enodes[E->mesh.nsd];
+ const int vpts = VPOINTS3D;
+ const int ends = ENODES3D;
const int dims=E->mesh.nsd;
const int sphere_key = 1;
get_global_shape_fn(E,el,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
- if ((el-1)%E->lmesh.ELZ[lev]==0)
+ if (iconv || (el-1)%E->lmesh.ELZ[lev]==0)
construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,0);
/* Note N[a].gauss_pt[n] is the value of shape fn a at the nth gaussian
@@ -191,70 +279,39 @@
for(k=1;k<=vpts;k++) {
W[k]=g_point[k].weight[dims-1]*dOmega.vpt[k]*E->EVI[lev][m][(el-1)*vpts+k];
+ }
- ra[k] = rtf[3][k];
- si[k] = sin(rtf[1][k]);
- ct[k] = cos(rtf[1][k])/si[k];
- }
+ get_ba(&(E->N), &GNx, &E->element_Cc, &E->element_Ccx,
+ rtf, E->mesh.nsd, ba);
-
- for(k=1;k<=VPOINTS3D;k++)
- for(a=1;a<=ends;a++)
- for(n=1;n<=dims;n++) {
- ba[a][k][n][1]=
- (GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(1,n,1,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)])*ra[k];
-
- ba[a][k][n][2]=
- (E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]*ct[k]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(2,n,2,a,k)])
- /si[k])*ra[k];
-
- ba[a][k][n][3]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)];
-
- ba[a][k][n][4]=
- (GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(2,n,1,a,k)]
- - E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]*ct[k]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(1,n,2,a,k)])
- /si[k])*ra[k];
-
- ba[a][k][n][5]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(1,n,a,k)]
- +(GNx.vpt[GNVXINDEX(0,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*(E->element_Ccx.vpt[BVXINDEX(3,n,1,a,k)]
- - E->element_Cc.vpt[BVINDEX(1,n,a,k)]))*ra[k];
-
- ba[a][k][n][6]=
- GNx.vpt[GNVXINDEX(2,a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- -ra[k]*E->N.vpt[GNVINDEX(a,k)]*E->element_Cc.vpt[BVINDEX(2,n,a,k)]
- +(GNx.vpt[GNVXINDEX(1,a,k)]*E->element_Cc.vpt[BVINDEX(3,n,a,k)]
- + E->N.vpt[GNVINDEX(a,k)]*E->element_Ccx.vpt[BVXINDEX(3,n,2,a,k)])
- /si[k]*ra[k];
- }
-
for(a=1;a<=ends;a++)
for(b=a;b<=ends;b++) {
- ad=dims*(a-1);
- bd=dims*(b-1);
-
bdbmu[1][1]=bdbmu[1][2]=bdbmu[1][3]=
bdbmu[2][1]=bdbmu[2][2]=bdbmu[2][3]=
bdbmu[3][1]=bdbmu[3][2]=bdbmu[3][3]=0.0;
+ if(E->control.inv_gruneisen > 0)
for(i=1;i<=dims;i++)
for(j=1;j<=dims;j++)
- for(j1=1;j1<=6;j1++)
- for(k=1;k<=VPOINTS3D;k++)
- bdbmu[i][j] +=
- W[k]*((j1>3)?1.0:2.0)*ba[a][k][i][j1]*ba[b][k][j][j1];
+ for(i1=1;i1<=6;i1++)
+ for(j1=1;j1<=6;j1++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] +=
+ W[k]*d2[i1][j1]*ba[a][k][i][i1]*ba[b][k][j][j1];
+ else
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(i1=1;i1<=6;i1++)
+ for(j1=1;j1<=6;j1++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] +=
+ W[k]*d1[i1][j1]*ba[a][k][i][i1]*ba[b][k][j][j1];
- /**/
+
+ /**/
+ ad=dims*(a-1);
+ bd=dims*(b-1);
+
pn=ad*nn+bd;
qn=bd*nn+ad;
@@ -518,14 +575,55 @@
return;
}
+
+
+
+/* compute div(rho_ref*V) = div(V) + Vz*d(ln(rho_ref))/dz */
+
+static void assemble_dlnrho(struct All_variables *E, double *dlnrhodr,
+ double **U, double **result, int level)
+{
+ int e, nz, j3, a, b, m;
+
+ const int nel = E->lmesh.NEL[level];
+ const int ends = enodes[E->mesh.nsd];
+ const int dims = E->mesh.nsd;
+ double tmp;
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(e=1; e<=nel; e++) {
+ //nz = ((e-1) % E->lmesh.elz) + 1;
+ //tmp = dlnrhodr[e] / ends;
+ for(a=1; a<=ends; a++) {
+ b = E->IEN[level][m][e].node[a];
+ j3 = E->ID[level][m][b].doff[3];
+ result[m][e] += dlnrhodr[e] * U[m][j3];
+ }
+ }
+
+ return;
+}
+
+
+
+
+void assemble_div_rho_u(struct All_variables *E,
+ double **U, double **result, int level)
+{
+ void assemble_div_u();
+ assemble_div_u(E, U, result, level);
+ assemble_dlnrho(E, E->refstate.dlnrhodr, U, result, level);
+
+ return;
+}
+
+
/* ==========================================
Assemble a div_u vector element by element
========================================== */
-void assemble_div_u(E,U,divU,level)
- struct All_variables *E;
- double **U,**divU;
- int level;
+void assemble_div_u(struct All_variables *E,
+ double **U, double **divU, int level)
{
int e,j1,j2,j3,p,a,b,m;
@@ -552,7 +650,6 @@
}
}
-
return;
}
@@ -793,7 +890,7 @@
nodeb=E->ien[m][el].node[b];
if ((E->node[m][nodeb]&type)&&(E->sphere.cap[m].VB[j][nodeb]!=0.0)){
if(!got_elt_k) {
- get_elt_k(E,el,elt_k,E->mesh.levmax,m);
+ get_elt_k(E,el,elt_k,E->mesh.levmax,m,1);
got_elt_k = 1;
}
q = dims*(b-1)+j-1;
Modified: mc/3D/CitcomS/trunk/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Full_version_dependent.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -63,16 +63,11 @@
if(E->mesh.nsd != 3)
E->mesh.noy = 1;
- E->mesh.nnx[1] = E->mesh.nox;
- E->mesh.nnx[2] = E->mesh.noy;
- E->mesh.nnx[3] = E->mesh.noz;
E->mesh.elx = E->mesh.nox-1;
E->mesh.ely = E->mesh.noy-1;
E->mesh.elz = E->mesh.noz-1;
- E->mesh.nno = E->sphere.caps;
- for(d=1;d<=E->mesh.nsd;d++)
- E->mesh.nno *= E->mesh.nnx[d];
+ E->mesh.nno = E->sphere.caps*E->mesh.nox*E->mesh.noy*E->mesh.noz;
E->mesh.nel = E->sphere.caps*E->mesh.elx*E->mesh.elz*E->mesh.ely;
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -43,6 +43,7 @@
#include "citcom_init.h"
#include "initial_temperature.h"
#include "lith_age.h"
+#include "material_properties.h"
#include "output.h"
#include "output_h5.h"
#include "parallel_related.h"
@@ -54,6 +55,7 @@
void allocate_common_vars(struct All_variables*);
void allocate_velocity_vars(struct All_variables*);
void check_bc_consistency(struct All_variables*);
+void construct_elt_gs(struct All_variables*);
void construct_id(struct All_variables*);
void construct_ien(struct All_variables*);
void construct_lm(struct All_variables*);
@@ -124,7 +126,10 @@
construct_sub_element(E);
construct_shape_functions(E);
+ construct_elt_gs(E);
+ reference_state(E);
+
/* this matrix results from spherical geometry */
/* construct_c3x3matrix(E); */
@@ -215,7 +220,6 @@
else if ( strcmp(E->control.PROBLEM_TYPE,"convection-chemical") == 0) {
E->control.CONVECTION = 1;
- E->control.CHEMISTRY_MODULE=1;
set_convection_defaults(E);
}
@@ -378,6 +382,10 @@
input_int("piterations",&(E->control.p_iterations),"100,0,nomax",m);
input_float("rayleigh",&(E->control.Atemp),"essential",m);
+ input_float("dissipation_number",&(E->control.disptn_number),"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);
@@ -475,6 +483,13 @@
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));
+
+ /* lump mass matrix for the energy eqn */
+ E->TMass[j] = (double *) malloc((nno+1)*sizeof(double));
+
nxyz = max(nox*noz,nox*noy);
nxyz = 2*max(nxyz,noz*noy);
@@ -484,6 +499,10 @@
} /* end for cap j */
+ /* density field */
+ E->rho = (double *) malloc((nno+1)*sizeof(double));
+
+ /* horizontal average */
E->Have.T = (float *)malloc((E->lmesh.noz+2)*sizeof(float));
E->Have.V[1] = (float *)malloc((E->lmesh.noz+2)*sizeof(float));
E->Have.V[2] = (float *)malloc((E->lmesh.noz+2)*sizeof(float));
@@ -597,11 +616,16 @@
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++)
E->P[j][i] = 0.0;
+ mat_prop_allocate(E);
phase_change_allocate(E);
set_up_nonmg_aliases(E,j);
@@ -692,12 +716,8 @@
E->control.v_steps_low = 10;
E->control.v_steps_upper = 1;
- E->control.max_res_red_each_p_mg = 1.0e-3;
E->control.accuracy = 1.0e-6;
E->control.vaccuracy = 1.0e-8;
- E->control.true_vcycle=0;
- E->control.depth_dominated=0;
- E->control.eqn_zigzag=0;
E->control.verbose=0; /* debugging/profiles */
/* SECOND: values for which an obvious default setting is useful */
@@ -706,11 +726,9 @@
E->control.ORTHOZ = 1; /* for orthogonal meshes by default */
- E->control.KERNEL = 0;
E->control.stokes=0;
E->control.restart=0;
E->control.CONVECTION = 0;
- E->control.SLAB = 0;
E->control.CART2D = 0;
E->control.CART3D = 0;
E->control.CART2pt5D = 0;
@@ -722,14 +740,7 @@
E->control.augmented_Lagr = 0;
E->control.augmented = 0.0;
- /* Default: all optional modules set to `off' */
- E->control.MELTING_MODULE = 0;
- E->control.CHEMISTRY_MODULE = 0;
-
E->control.GRID_TYPE=1;
- E->mesh.hwidth[1]=E->mesh.hwidth[2]=E->mesh.hwidth[3]=1.0; /* divide by this one ! */
- E->mesh.magnitude[1]=E->mesh.magnitude[2]=E->mesh.magnitude[3]=0.0;
- E->mesh.offset[1]=E->mesh.offset[2]=E->mesh.offset[3]=0.0;
E->parallel.nprocx=1; E->parallel.nprocz=1; E->parallel.nprocy=1;
@@ -745,7 +756,6 @@
E->sphere.ri = 0.5;
E->control.precondition = 0; /* for larger visc contrasts turn this back on */
- E->control.vprecondition = 1;
E->mesh.toptbc = 1; /* fixed t */
E->mesh.bottbc = 1;
Modified: mc/3D/CitcomS/trunk/lib/Makefile.am
===================================================================
--- mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Makefile.am 2007-08-16 00:27:43 UTC (rev 7831)
@@ -83,6 +83,8 @@
interuption.h \
lith_age.h \
Lith_age.c \
+ Material_properties.c \
+ material_properties.h \
Nodal_mesh.c \
Output.c \
output.h \
Modified: mc/3D/CitcomS/trunk/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Obsolete.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -896,6 +896,75 @@
}
+
+/* ========================================================== */
+/* From Pan_problem_misc_functions.c */
+/* =========================================================== */
+
+double SIN_D(x)
+ double x;
+{
+#if defined(__osf__)
+ return sind(x);
+#else
+ return sin((x/180.0) * M_PI);
+#endif
+
+}
+
+double COT_D(x)
+ double x;
+{
+#if defined(__osf__)
+ return cotd(x);
+#else
+ return tan(((90.0-x)/180.0) * M_PI);
+#endif
+
+}
+
+
+/* non-runaway malloc */
+
+void * Malloc1(bytes,file,line)
+ int bytes;
+ char *file;
+ int line;
+{
+ void *ptr;
+
+ ptr = malloc((size_t)bytes);
+ if (ptr == (void *)NULL) {
+ fprintf(stderr,"Memory: cannot allocate another %d bytes \n(line %d of file %s)\n",bytes,line,file);
+ parallel_process_termination();
+ }
+
+ return(ptr);
+}
+
+
+/* returns the out of plane component of the cross product of
+ the two vectors assuming that one is looking AGAINST the
+ direction of the axis of D, anti-clockwise angles
+ are positive (are you sure ?), and the axes are ordered 2,3 or 1,3 or 1,2 */
+
+
+float cross2d(x11,x12,x21,x22,D)
+ float x11,x12,x21,x22;
+ int D;
+{
+ float temp;
+ if(1==D)
+ temp = ( x11*x22-x12*x21);
+ if(2==D)
+ temp = (-x11*x22+x12*x21);
+ if(3==D)
+ temp = ( x11*x22-x12*x21);
+
+ return(temp);
+}
+
+
/* version */
/* $Id$ */
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -49,6 +49,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);
FILE* output_open_mode(char *, char *);
extern void parallel_process_termination();
@@ -90,12 +91,15 @@
/*output_mat(E);*/
}
-
+
output_velo(E, cycles);
output_visc(E, cycles);
output_surf_botm(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)
@@ -503,6 +507,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/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -1310,6 +1310,7 @@
int rank;
hsize_t *dims;
double *data;
+ float tmp;
input = h5create_group(E->hdf5.file_id, "input", (size_t)0);
@@ -1458,6 +1459,13 @@
status = set_attribute_string(input, "datafile_old", E->control.old_P_file);
status = set_attribute_float(input, "rayleigh", E->control.Atemp);
+ status = set_attribute_float(input, "dissipation_number", E->control.disptn_number);
+ 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/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -137,15 +137,20 @@
void get_buoyancy(struct All_variables *E, double **buoy)
{
- int i,m;
+ int i,j,m;
double temp,temp2;
void remove_horiz_ave2(struct All_variables*, double**);
temp = E->control.Atemp;
- for(m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=1;i<=E->lmesh.nno;i++)
- buoy[m][i] = temp * E->T[m][i];
+ for(i=1;i<=E->lmesh.nno;i++) {
+ int nz = ((i-1) % E->lmesh.noz) + 1;
+ /* 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->refstate.rho[nz]
+ * E->refstate.expansivity[nz] * E->T[m][i];
+ }
/* chemical buoyancy */
if(E->control.tracer &&
@@ -160,53 +165,20 @@
phase_change_apply_670(E, buoy);
phase_change_apply_cmb(E, buoy);
+ /* convert density to buoyancy */
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=E->lmesh.noz;i++)
+ for(j=0;j<E->lmesh.nox*E->lmesh.noy;j++) {
+ int n = j*E->lmesh.noz + i;
+ buoy[m][n] *= E->refstate.gravity[i];
+ }
+
remove_horiz_ave2(E,buoy);
return;
}
-double SIN_D(x)
- double x;
-{
-#if defined(__osf__)
- return sind(x);
-#else
- return sin((x/180.0) * M_PI);
-#endif
-}
-
-double COT_D(x)
- double x;
-{
-#if defined(__osf__)
- return cotd(x);
-#else
- return tan(((90.0-x)/180.0) * M_PI);
-#endif
-
-}
-
-
-/* non-runaway malloc */
-
-void * Malloc1(bytes,file,line)
- int bytes;
- char *file;
- int line;
-{
- void *ptr;
-
- ptr = malloc((size_t)bytes);
- if (ptr == (void *)NULL) {
- fprintf(stderr,"Memory: cannot allocate another %d bytes \n(line %d of file %s)\n",bytes,line,file);
- parallel_process_termination();
- }
-
- return(ptr);
-}
-
-
/* Read in a file containing previous values of a field. The input in the parameter
file for this should look like: `previous_name_file=string' and `previous_name_column=int'
where `name' is substituted by the argument of the function.
@@ -230,8 +202,6 @@
{
int input_string();
- float cross2d();
-
char discard[5001];
char *token;
char *filename;
@@ -349,27 +319,6 @@
}
-/* returns the out of plane component of the cross product of
- the two vectors assuming that one is looking AGAINST the
- direction of the axis of D, anti-clockwise angles
- are positive (are you sure ?), and the axes are ordered 2,3 or 1,3 or 1,2 */
-
-
-float cross2d(x11,x12,x21,x22,D)
- float x11,x12,x21,x22;
- int D;
-{
- float temp;
- if(1==D)
- temp = ( x11*x22-x12*x21);
- if(2==D)
- temp = (-x11*x22+x12*x21);
- if(3==D)
- temp = ( x11*x22-x12*x21);
-
- return(temp);
-}
-
/* =================================================
my version of arc tan
=================================================*/
Modified: mc/3D/CitcomS/trunk/lib/Process_buoyancy.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Process_buoyancy.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -107,6 +107,7 @@
uT = 0.0;
area = 0.0;
for(i=1;i<=vpts;i++) {
+ /* XXX: missing unit conversion, heat capacity and conductivity */
uT += u[i]*T[i]*dOmega.vpt[i] + dTdz[i]*dOmega.vpt[i];
}
Modified: mc/3D/CitcomS/trunk/lib/Regional_obsolete.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Regional_obsolete.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -926,10 +926,6 @@
int m,i,it;
- E->monitor.length_scale = E->data.radius_km/E->mesh.layer[2]; /* km */
- E->monitor.time_scale = pow(E->data.radius_km*1000.0,2.0)/ /* Million years */
- (E->data.therm_diff*3600.0*24.0*365.25*1.0e6);
-
if ( (ii == 0) || ((ii % E->control.record_every) == 0)
|| E->control.DIRECTII) {
get_STD_topo(E,E->slice.tpg,E->slice.tpgb,E->slice.divg,E->slice.vort,ii);
Modified: mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Regional_version_dependent.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -64,21 +64,13 @@
if(E->mesh.nsd != 3)
E->mesh.noy = 1;
- E->mesh.nnx[1] = E->mesh.nox;
- E->mesh.nnx[2] = E->mesh.noy;
- E->mesh.nnx[3] = E->mesh.noz;
E->mesh.elx = E->mesh.nox-1;
E->mesh.ely = E->mesh.noy-1;
E->mesh.elz = E->mesh.noz-1;
- E->mesh.nno = E->sphere.caps;
- for(d=1;d<=E->mesh.nsd;d++)
- E->mesh.nno *= E->mesh.nnx[d];
+ E->mesh.nno = E->sphere.caps*E->mesh.nox*E->mesh.noy*E->mesh.noz;
-/*
E->mesh.nel = E->sphere.caps*E->mesh.elx*E->mesh.elz*E->mesh.ely;
-*/
- E->mesh.nel = E->mesh.elx*E->mesh.elz*E->mesh.ely;
E->mesh.nnov = E->mesh.nno;
Modified: mc/3D/CitcomS/trunk/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -867,141 +867,189 @@
Na be delta(a,b)
========================================== */
-void mass_matrix(E)
- struct All_variables *E;
+void mass_matrix(struct All_variables *E)
+{
+ int m,node,i,nint,e,lev;
+ int n[9], nz;
+ void get_global_shape_fn();
+ double myatan(),rtf[4][9],area,centre[4],temp[9],dx1,dx2,dx3;
+ struct Shape_function GN;
+ struct Shape_function_dA dOmega;
+ struct Shape_function_dx GNx;
-{ int m,node,i,nint,e,lev;
- int n[9];
- void get_global_shape_fn();
- double myatan(),rtf[4][9],area,centre[4],temp[9],dx1,dx2,dx3;
- struct Shape_function GN;
- struct Shape_function_dA dOmega;
- struct Shape_function_dx GNx;
+ const int vpts=vpoints[E->mesh.nsd];
+ const int sphere_key=1;
- const int vpts=vpoints[E->mesh.nsd];
- const int sphere_key=1;
+ /* ECO .size can also be defined here */
- /* ECO .size can also be defined here */
+ for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ for(node=1;node<=E->lmesh.NNO[lev];node++)
+ E->MASS[lev][m][node] = 0.0;
- for(node=1;node<=E->lmesh.NNO[lev];node++)
- E->MASS[lev][m][node] = 0.0;
+ for(e=1;e<=E->lmesh.NEL[lev];e++) {
- for(e=1;e<=E->lmesh.NEL[lev];e++) {
+ get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
- get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,sphere_key,rtf,lev,m);
+ area = centre[1] = centre[2] = centre[3] = 0.0;
- area = centre[1] = centre[2] = centre[3] = 0.0;
+ for(node=1;node<=enodes[E->mesh.nsd];node++)
+ n[node] = E->IEN[lev][m][e].node[node];
- for(node=1;node<=enodes[E->mesh.nsd];node++)
- n[node] = E->IEN[lev][m][e].node[node];
+ for(i=1;i<=E->mesh.nsd;i++) {
+ for(node=1;node<=enodes[E->mesh.nsd];node++)
+ centre[i] += E->X[lev][m][i][n[node]];
- for(i=1;i<=E->mesh.nsd;i++) {
- for(node=1;node<=enodes[E->mesh.nsd];node++)
- centre[i] += E->X[lev][m][i][E->IEN[lev][m][e].node[node]];
+ centre[i] = centre[i]/enodes[E->mesh.nsd];
+ } /* end for i */
- centre[i] = centre[i]/enodes[E->mesh.nsd];
- } /* end for i */
+ /* dx3=radius, dx1=theta, dx2=phi */
+ dx3 = sqrt(centre[1]*centre[1]+centre[2]*centre[2]+centre[3]*centre[3]);
+ dx1 = acos( centre[3]/dx3 );
+ dx2 = myatan(centre[2],centre[1]);
- dx3 = sqrt(centre[1]*centre[1]+centre[2]*centre[2]+centre[3]*centre[3]);
- dx1 = acos( centre[3]/dx3 );
- dx2 = myatan(centre[2],centre[1]);
+ /* center of this element in the spherical coordinate */
+ E->ECO[lev][m][e].centre[1] = dx1;
+ E->ECO[lev][m][e].centre[2] = dx2;
+ E->ECO[lev][m][e].centre[3] = dx3;
- E->ECO[lev][m][e].centre[1] = dx1;
- E->ECO[lev][m][e].centre[2] = dx2;
- E->ECO[lev][m][e].centre[3] = dx3;
+ /* delta(theta) of this element */
+ dx1 = max( fabs(E->SX[lev][m][1][n[3]]-E->SX[lev][m][1][n[1]]),
+ fabs(E->SX[lev][m][1][n[2]]-E->SX[lev][m][1][n[4]]) );
- dx1 = max( fabs(E->SX[lev][m][1][n[3]]-E->SX[lev][m][1][n[1]]),
- fabs(E->SX[lev][m][1][n[2]]-E->SX[lev][m][1][n[4]]) );
- E->ECO[lev][m][e].size[1] = dx1*E->ECO[lev][m][e].centre[3];
+ /* length of this element in the theta-direction */
+ E->ECO[lev][m][e].size[1] = dx1*E->ECO[lev][m][e].centre[3];
- dx1 = fabs(E->SX[lev][m][2][n[3]]-E->SX[lev][m][2][n[1]]);
- if (dx1>M_PI)
- dx1 = min(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) + 2.0*M_PI -
- max(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) ;
+ /* delta(phi) of this element */
+ dx1 = fabs(E->SX[lev][m][2][n[3]]-E->SX[lev][m][2][n[1]]);
+ if (dx1>M_PI)
+ dx1 = min(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) + 2.0*M_PI -
+ max(E->SX[lev][m][2][n[3]],E->SX[lev][m][2][n[1]]) ;
- dx2 = fabs(E->SX[lev][m][2][n[2]]-E->SX[lev][m][2][n[4]]);
- if (dx2>M_PI)
- dx2 = min(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) + 2.0*M_PI -
- max(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) ;
+ dx2 = fabs(E->SX[lev][m][2][n[2]]-E->SX[lev][m][2][n[4]]);
+ if (dx2>M_PI)
+ dx2 = min(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) + 2.0*M_PI -
+ max(E->SX[lev][m][2][n[2]],E->SX[lev][m][2][n[4]]) ;
- dx2 = max(dx1,dx2);
+ dx2 = max(dx1,dx2);
- E->ECO[lev][m][e].size[2] = dx2*E->ECO[lev][m][e].centre[3]
- *sin(E->ECO[lev][m][e].centre[1]);
+ /* length of this element in the phi-direction */
+ E->ECO[lev][m][e].size[2] = dx2*E->ECO[lev][m][e].centre[3]
+ *sin(E->ECO[lev][m][e].centre[1]);
- dx3 = 0.25*(
- fabs(E->SX[lev][m][3][n[5]]+E->SX[lev][m][3][n[6]]
- +E->SX[lev][m][3][n[7]]+E->SX[lev][m][3][n[8]]
- -E->SX[lev][m][3][n[1]]-E->SX[lev][m][3][n[2]]
- -E->SX[lev][m][3][n[3]]-E->SX[lev][m][3][n[4]]));
+ /* delta(radius) of this element */
+ dx3 = 0.25*(fabs(E->SX[lev][m][3][n[5]]+E->SX[lev][m][3][n[6]]
+ +E->SX[lev][m][3][n[7]]+E->SX[lev][m][3][n[8]]
+ -E->SX[lev][m][3][n[1]]-E->SX[lev][m][3][n[2]]
+ -E->SX[lev][m][3][n[3]]-E->SX[lev][m][3][n[4]]));
- E->ECO[lev][m][e].size[3] = dx3;
+ /* length of this element in the radius-direction */
+ E->ECO[lev][m][e].size[3] = dx3;
- for(nint=1;nint<=vpts;nint++)
- area += g_point[nint].weight[E->mesh.nsd-1] * dOmega.vpt[nint];
- E->ECO[lev][m][e].area = area;
+ /* volume (area in 2D) of this element */
+ for(nint=1;nint<=vpts;nint++)
+ area += g_point[nint].weight[E->mesh.nsd-1] * dOmega.vpt[nint];
+ E->ECO[lev][m][e].area = area;
- for(node=1;node<=enodes[E->mesh.nsd];node++) {
- temp[node] = 0.0;
- for(nint=1;nint<=vpts;nint++)
- temp[node] += dOmega.vpt[nint]*g_point[nint].weight[E->mesh.nsd-1]
- *E->N.vpt[GNVINDEX(node,nint)]; /* int Na dV */
- }
+ for(node=1;node<=enodes[E->mesh.nsd];node++) {
+ temp[node] = 0.0;
+ for(nint=1;nint<=vpts;nint++)
+ temp[node] += dOmega.vpt[nint]*g_point[nint].weight[E->mesh.nsd-1]
+ *E->N.vpt[GNVINDEX(node,nint)]; /* int Na dV */
+ }
- for(node=1;node<=enodes[E->mesh.nsd];node++)
- E->MASS[lev][m][E->IEN[lev][m][e].node[node]] += temp[node];
+ for(node=1;node<=enodes[E->mesh.nsd];node++)
+ E->MASS[lev][m][E->IEN[lev][m][e].node[node]] += temp[node];
- for(node=1;node<=enodes[E->mesh.nsd];node++)
- E->TWW[lev][m][e].node[node] = temp[node];
+ /* weight of each node, equivalent to pmass in ConMan */
+ for(node=1;node<=enodes[E->mesh.nsd];node++)
+ E->TWW[lev][m][e].node[node] = temp[node];
- } /* end of ele*/
+ } /* end of ele*/
- } /* m */
+ } /* end of for m */
- if (E->control.NMULTIGRID||E->control.EMULTIGRID||E->mesh.levmax==lev)
- (E->exchange_node_f)(E,E->MASS[lev],lev);
+ if (E->control.NMULTIGRID||E->control.EMULTIGRID||E->mesh.levmax==lev)
+ (E->exchange_node_f)(E,E->MASS[lev],lev);
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(node=1;node<=E->lmesh.NNO[lev];node++)
- E->MASS[lev][m][node] = 1.0/E->MASS[lev][m][node];
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(node=1;node<=E->lmesh.NNO[lev];node++)
+ E->MASS[lev][m][node] = 1.0/E->MASS[lev][m][node];
- } /* lev */
+ } /* end of for lev */
- /* compute volume of this processor mesh and the whole mesh */
- E->lmesh.volume = 0;
- E->mesh.volume = 0;
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(e=1;e<=E->lmesh.nel;e++)
- E->lmesh.volume += E->eco[m][e].area;
+ for(node=1;node<=E->lmesh.nno;node++)
+ E->TMass[m][node] = 0.0;
- MPI_Allreduce(&E->lmesh.volume, &E->mesh.volume, 1, MPI_DOUBLE, MPI_SUM, E->parallel.world);
+ for(e=1;e<=E->lmesh.nel;e++) {
+ get_global_shape_fn(E,e,&GN,&GNx,&dOmega,0,
+ sphere_key,rtf,E->mesh.levmax,m);
+ for(node=1;node<=enodes[E->mesh.nsd];node++) {
+ temp[node] = 0.0;
+ nz = ((E->ien[m][e].node[node]-1) % E->lmesh.noz) + 1;
+ for(nint=1;nint<=vpts;nint++) {
+ temp[node] += E->refstate.rho[nz]
+ * E->refstate.capacity[nz]
+ * dOmega.vpt[nint]
+ * g_point[nint].weight[E->mesh.nsd-1]
+ * E->N.vpt[GNVINDEX(node,nint)]; /* int Na dV */
+ }
+ }
+ /* lumped mass matrix, equivalent to tmass in ConMan */
+ for(node=1;node<=enodes[E->mesh.nsd];node++)
+ E->TMass[m][E->ien[m][e].node[node]] += temp[node];
- if (E->control.verbose) {
- fprintf(E->fp_out, "rank=%d my_volume=%e total_volume=%e\n",
- E->parallel.me, E->lmesh.volume, E->mesh.volume);
+ } /* end of for e */
+ } /* end of for m */
- for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
- fprintf(E->fp_out,"output_mass lev=%d\n",lev);
- for (m=1;m<=E->sphere.caps_per_proc;m++) {
- fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
- for(e=1;e<=E->lmesh.NEL[lev];e++)
- fprintf(E->fp_out,"%d %g \n",e,E->ECO[lev][m][e].area);
- for (node=1;node<=E->lmesh.NNO[lev];node++)
- fprintf(E->fp_out,"Mass[%d]= %g \n",node,E->MASS[lev][m][node]);
- }
- }
- fflush(E->fp_out);
- }
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(node=1;node<=E->lmesh.nno;node++)
+ E->TMass[m][node] = 1.0 / E->TMass[m][node];
- return;
+
+ /* compute volume of this processor mesh and the whole mesh */
+ E->lmesh.volume = 0;
+ E->mesh.volume = 0;
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(e=1;e<=E->lmesh.nel;e++)
+ E->lmesh.volume += E->eco[m][e].area;
+
+ MPI_Allreduce(&E->lmesh.volume, &E->mesh.volume, 1, MPI_DOUBLE,
+ MPI_SUM, E->parallel.world);
+
+
+ if (E->control.verbose) {
+ fprintf(E->fp_out, "rank=%d my_volume=%e total_volume=%e\n",
+ E->parallel.me, E->lmesh.volume, E->mesh.volume);
+
+ for(lev=E->mesh.levmin;lev<=E->mesh.levmax;lev++) {
+ fprintf(E->fp_out,"output_mass lev=%d\n",lev);
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
+ for(e=1;e<=E->lmesh.NEL[lev];e++)
+ fprintf(E->fp_out,"%d %g \n",e,E->ECO[lev][m][e].area);
+ for (node=1;node<=E->lmesh.NNO[lev];node++)
+ fprintf(E->fp_out,"Mass[%d]= %g \n",node,E->MASS[lev][m][node]);
+ }
+ }
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ fprintf(E->fp_out,"m=%d %d \n",E->sphere.capid[m],m);
+ for (node=1;node<=E->lmesh.nno;node++)
+ fprintf(E->fp_out,"TMass[%d]= %g \n",node,E->TMass[m][node]);
+ }
+ fflush(E->fp_out);
+ }
+
+ return;
}
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -35,13 +35,28 @@
#include "global_defs.h"
#include <stdlib.h>
+static float solve_Ahat_p_fhat(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max);
+static float solve_Ahat_p_fhat_BA(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max);
+static float solve_Ahat_p_fhat_TALA(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max);
+static double initial_vel_residual(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp);
+static double incompressibility_residual(struct All_variables *E,
+ double **V, double **r);
+
+
/* Master loop for pressure and (hence) velocity field */
void solve_constrained_flow_iterative(E)
struct All_variables *E;
{
- float solve_Ahat_p_fhat();
void v_from_vector();
void p_to_nodes();
@@ -63,7 +78,6 @@
struct All_variables *E;
{
- float solve_Ahat_p_fhat();
void v_from_vector_pseudo_surf();
void p_to_nodes();
@@ -84,21 +98,38 @@
/* ========================================================================= */
-float solve_Ahat_p_fhat(struct All_variables *E,
- double **V, double **P, double **FF,
- double imp, int *steps_max)
+static float solve_Ahat_p_fhat(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max)
{
+ float residual;
+
+ if(E->control.inv_gruneisen < 1e-6)
+ residual = solve_Ahat_p_fhat_BA(E, V, P, F, imp, steps_max);
+ else
+ residual = solve_Ahat_p_fhat_TALA(E, V, P, F, imp, steps_max);
+
+ return(residual);
+}
+
+
+/* Solve incompressible Stokes flow (Boussinesq Approximation) using
+ * conjugate gradient (CG) iterations
+ */
+
+static float solve_Ahat_p_fhat_BA(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max)
+{
int m, i, j, count, valid, lev, npno, neq;
int gnpno, gneq;
- double *r1[NCS], *F[NCS];
- double *r0[NCS], *r2[NCS], *z0[NCS], *z1[NCS], *s1[NCS], *s2[NCS];
+ double *r1[NCS], *r2[NCS], *z1[NCS], *s1[NCS], *s2[NCS];
double *shuffle[NCS];
- double alpha, delta, s2dotAhat, r0dotr0, r1dotz1;
+ double alpha, delta, r0dotz0, r1dotz1;
double residual, v_res;
double global_vdot(), global_pdot();
- double *dvector();
double time0, CPU_time0();
float dpressure, dvelocity;
@@ -114,81 +145,37 @@
gneq = E->mesh.neq;
npno = E->lmesh.npno;
neq = E->lmesh.neq;
+ lev = E->mesh.levmax;
- for (m=1; m<=E->sphere.caps_per_proc; m++) {
- F[m] = (double *)malloc((neq+1)*sizeof(double));
-
- r0[m] = (double *)malloc((npno+1)*sizeof(double));
+ for (m=1; m<=E->sphere.caps_per_proc; m++) {
r1[m] = (double *)malloc((npno+1)*sizeof(double));
r2[m] = (double *)malloc((npno+1)*sizeof(double));
- z0[m] = (double *)malloc((npno+1)*sizeof(double));
z1[m] = (double *)malloc((npno+1)*sizeof(double));
s1[m] = (double *)malloc((npno+1)*sizeof(double));
s2[m] = (double *)malloc((npno+1)*sizeof(double));
}
- /* Copy the original force vector. FF shouldn't be modified. */
- for (m=1;m<=E->sphere.caps_per_proc;m++)
- for(i=0;i<neq;i++)
- F[m][i] = FF[m][i];
-
time0 = CPU_time0();
+ count = 0;
-
/* calculate the initial velocity residual */
- lev = E->mesh.levmax;
- v_res = sqrt(global_vdot(E, F, F, lev)/gneq);
+ v_res = initial_vel_residual(E, V, P, F, imp);
- if (E->parallel.me==0) {
- fprintf(E->fp, "initial residue of momentum equation F %.9e %d\n",
- v_res, gneq);
- fprintf(stderr, "initial residue of momentum equation F %.9e %d\n",
- v_res, gneq);
- }
-
- /* F = F - grad(P) - K*V */
- assemble_grad_p(E, P, E->u1, lev);
- for(m=1; m<=E->sphere.caps_per_proc; m++)
- for(i=0; i<neq; i++)
- F[m][i] = F[m][i] - E->u1[m][i];
-
- assemble_del2_u(E, V, E->u1, lev, 1);
- for(m=1; m<=E->sphere.caps_per_proc; m++)
- for(i=0; i<neq; i++)
- F[m][i] = F[m][i] - E->u1[m][i];
-
- strip_bcs_from_residual(E, F, lev);
-
-
- /* solve K*u1 = F for u1 */
- valid=solve_del2_u(E, E->u1, F, imp*v_res, E->mesh.levmax);
- strip_bcs_from_residual(E, E->u1, lev);
-
-
- /* V = V + u1 */
- for(m=1; m<=E->sphere.caps_per_proc; m++)
- for(i=0; i<neq; i++)
- V[m][i] += E->u1[m][i];
-
-
- /* r1 = div(V) */
+ /* initial residual r1 = div(V) */
assemble_div_u(E, V, r1, lev);
+ residual = incompressibility_residual(E, V, r1);
- /* incompressiblity residual = norm(r1) / norm(V) */
- residual = sqrt(global_pdot(E, r1, r1, lev)/gnpno);
- E->monitor.vdotv = sqrt(global_vdot(E, V, V, lev)/gneq);
- E->monitor.incompressibility = residual / E->monitor.vdotv;
-
- count = 0;
-
- if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
- "for step %d\n", count, CPU_time0()-time0,
- E->monitor.incompressibility, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
- "for step %d\n", count, CPU_time0()-time0,
- E->monitor.incompressibility, E->monitor.solution_cycles);
+ if (E->control.print_convergence && E->parallel.me==0) {
+ fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ 0.0, 0.0, E->monitor.solution_cycles);
+ fflush(E->fp);
+ fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ 0.0, 0.0, E->monitor.solution_cycles);
}
@@ -196,12 +183,15 @@
dpressure = 1.0;
dvelocity = 1.0;
+ valid = 1;
+ r0dotz0 = 0;
+
while( (valid) && (count < *steps_max) &&
(E->monitor.incompressibility >= E->control.tole_comp) &&
(dpressure >= imp) && (dvelocity >= imp) ) {
- /* preconditioner B, z1 = B*r1 */
+ /* preconditioner BPI ~= inv(K), z1 = BPI*r1 */
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
z1[m][j] = E->BPI[lev][m][j] * r1[m][j];
@@ -209,48 +199,54 @@
/* r1dotz1 = <r1, z1> */
r1dotz1 = global_pdot(E, r1, z1, lev);
+ assert(r1dotz1 != 0.0 /* Division by zero in head of incompressibility iteration */);
- if ((count == 0))
+ /* update search direction */
+ if(count == 0)
for (m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
s2[m][j] = z1[m][j];
else {
/* s2 = z1 + s1 * <r1,z1>/<r0,z0> */
- r0dotr0 = global_pdot(E, r0, z0, lev);
- assert(r0dotr0 != 0.0 /* Division by zero in head of incompressibility iteration */);
- delta = r1dotz1 / r0dotr0;
+ delta = r1dotz1 / r0dotz0;
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
s2[m][j] = z1[m][j] + delta * s1[m][j];
}
+
/* solve K*u1 = grad(s2) for u1 */
assemble_grad_p(E, s2, F, lev);
valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
strip_bcs_from_residual(E, E->u1, lev);
- /* alpha = <r1, z1> / <s2, div(u1)> */
+ /* F = div(u1) */
assemble_div_u(E, E->u1, F, lev);
- s2dotAhat = global_pdot(E, s2, F, lev);
+
+ /* alpha = <r1, z1> / <s2, F> */
if(valid)
/* alpha defined this way is the same as R&W */
- alpha = r1dotz1 / s2dotAhat;
+ alpha = r1dotz1 / global_pdot(E, s2, F, lev);
else
alpha = 0.0;
/* r2 = r1 - alpha * div(u1) */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ r2[m][j] = r1[m][j] - alpha * F[m][j];
+
+
/* P = P + alpha * s2 */
- /* V = V - alpha * u1 */
for(m=1; m<=E->sphere.caps_per_proc; m++)
- for(j=1; j<=npno; j++) {
- r2[m][j] = r1[m][j] - alpha * F[m][j];
+ for(j=1; j<=npno; j++)
P[m][j] += alpha * s2[m][j];
- }
+
+ /* V = V - alpha * u1 */
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=0; j<neq; j++)
V[m][j] -= alpha * E->u1[m][j];
@@ -258,13 +254,8 @@
/* compute velocity and incompressibility residual */
assemble_div_u(E, V, F, lev);
- E->monitor.vdotv = global_vdot(E, V, V, E->mesh.levmax);
- E->monitor.incompressibility = sqrt((gneq/gnpno)
- *(1.0e-32
- + global_pdot(E, F, F, lev)
- / (1.0e-32+E->monitor.vdotv)));
+ incompressibility_residual(E, V, F);
-
/* compute velocity and pressure corrections */
dpressure = alpha * sqrt(global_pdot(E, s2, s2, lev)
/ (1.0e-32 + global_pdot(E, P, P, lev)));
@@ -285,29 +276,25 @@
}
- /* swap array pointers */
+ /* shift array pointers */
for(m=1; m<=E->sphere.caps_per_proc; m++) {
shuffle[m] = s1[m];
s1[m] = s2[m];
s2[m] = shuffle[m];
- shuffle[m] = r0[m];
- r0[m] = r1[m];
+ shuffle[m] = r1[m];
r1[m] = r2[m];
r2[m] = shuffle[m];
-
- shuffle[m] = z0[m];
- z0[m] = z1[m];
- z1[m] = shuffle[m];
}
- } /* end loop for conjugate gradient */
+ /* shift <r0, z0> = <r1, z1> */
+ r0dotz0 = r1dotz1;
+ } /* end loop for conjugate gradient */
+
for(m=1; m<=E->sphere.caps_per_proc; m++) {
- free((void *) r0[m]);
free((void *) r1[m]);
free((void *) r2[m]);
- free((void *) z0[m]);
free((void *) z1[m]);
free((void *) s1[m]);
free((void *) s2[m]);
@@ -318,8 +305,353 @@
return(residual);
}
-/* ========================================================================== */
+/* Solve incompressible Stokes flow (Boussinesq Approximation) using
+ * bi-conjugate gradient stablized (BiCG-stab)iterations
+ */
+static float solve_Ahat_p_fhat_TALA(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max)
+{
+ void assemble_div_u();
+ void assemble_del2_u();
+ void assemble_grad_p();
+ void strip_bcs_from_residual();
+ int solve_del2_u();
+ void parallel_process_termination();
+ double global_vdot(), global_pdot();
+ double CPU_time0();
+ int gnpno, gneq;
+ int npno, neq;
+ int m, i, j, count, lev;
+ int valid;
+ double alpha, beta, omega;
+ double r0dotrt, r1dotrt;
+ double residual, dpressure, dvelocity;
+
+ double *r1[NCS], *r2[NCS], *pt[NCS], *p1[NCS], *p2[NCS];
+ double *rt[NCS], *v0[NCS], *s0[NCS], *st[NCS], *t0[NCS];
+ double *u0[NCS];
+ double *shuffle[NCS];
+
+ double time0, v_res;
+
+ gnpno = E->mesh.npno;
+ gneq = E->mesh.neq;
+ npno = E->lmesh.npno;
+ neq = E->lmesh.neq;
+ lev = E->mesh.levmax;
+
+ for (m=1; m<=E->sphere.caps_per_proc; m++) {
+ r1[m] = (double *)malloc((npno+1)*sizeof(double));
+ r2[m] = (double *)malloc((npno+1)*sizeof(double));
+ pt[m] = (double *)malloc((npno+1)*sizeof(double));
+ p1[m] = (double *)malloc((npno+1)*sizeof(double));
+ p2[m] = (double *)malloc((npno+1)*sizeof(double));
+ rt[m] = (double *)malloc((npno+1)*sizeof(double));
+ v0[m] = (double *)malloc((npno+1)*sizeof(double));
+ s0[m] = (double *)malloc((npno+1)*sizeof(double));
+ st[m] = (double *)malloc((npno+1)*sizeof(double));
+ t0[m] = (double *)malloc((npno+1)*sizeof(double));
+
+ u0[m] = (double *)malloc((neq+1)*sizeof(double));
+ }
+
+ time0 = CPU_time0();
+ count = 0;
+
+ /* calculate the initial velocity residual */
+ v_res = initial_vel_residual(E, V, P, F, imp);
+
+
+ /* initial residual r1 = div(rho_ref*V) */
+ assemble_div_rho_u(E, V, r1, lev);
+ residual = incompressibility_residual(E, V, r1);
+
+
+ /* initial conjugate residual rt = r1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ rt[m][j] = r1[m][j];
+
+
+ if (E->control.print_convergence && E->parallel.me==0) {
+ fprintf(E->fp,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ 0.0, 0.0, E->monitor.solution_cycles);
+ fflush(E->fp);
+ fprintf(stderr,"AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ 0.0, 0.0, E->monitor.solution_cycles);
+ }
+
+
+ /* pressure and velocity corrections */
+ dpressure = 1.0;
+ dvelocity = 1.0;
+
+ valid = 1;
+ r0dotrt = alpha = omega = 0;
+
+ while( (valid) && (count < *steps_max) &&
+ ((E->monitor.incompressibility >= E->control.tole_comp) &&
+ (dpressure >= imp) && (dvelocity >= imp)) ) {
+
+
+
+ /* r1dotrt = <r1, rt> */
+ r1dotrt = global_pdot(E, r1, rt, lev);
+ if(r1dotrt == 0.0) {
+ /* TODO: can we resume the computation even when BiCGstab failed?
+ */
+ fprintf(E->fp, "BiCGstab method failed!!\n");
+ fprintf(stderr, "BiCGstab method failed!!\n");
+ parallel_process_termination();
+ }
+
+
+ /* update search direction */
+ if(count == 0)
+ for (m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ p2[m][j] = r1[m][j];
+ else {
+ /* p2 = r1 + <r1,rt>/<r0,rt> * alpha/omega * (p1 - omega*v0) */
+ beta = (r1dotrt / r0dotrt) * (alpha / omega);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ p2[m][j] = r1[m][j] + beta
+ * (p1[m][j] - omega * v0[m][j]);
+ }
+
+
+ /* preconditioner BPI ~= inv(K), pt = BPI*p2 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ pt[m][j] = E->BPI[lev][m][j] * p2[m][j];
+
+
+ /* solve K*u0 = grad(pt) for u1 */
+ assemble_grad_p(E, pt, F, lev);
+ valid = solve_del2_u(E, u0, F, imp*v_res, lev);
+ if(!valid) fprintf(stderr, "not valid 1\n");
+ strip_bcs_from_residual(E, u0, lev);
+
+
+ /* v0 = div(rho_ref*u0) */
+ assemble_div_rho_u(E, u0, v0, lev);
+
+
+ /* alpha = r1dotrt / <rt, v0> */
+ alpha = r1dotrt / global_pdot(E, rt, v0, lev);
+
+
+ /* s0 = r1 - alpha * v0 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ s0[m][j] = r1[m][j] - alpha * v0[m][j];
+
+
+ /* stop iteration if norm(s) is small enough */
+ if(global_pdot(E, s0, s0, lev) < imp*gnpno) {
+ // is the check correct?
+ // update solution, TODO
+ //break;
+ }
+
+
+ /* preconditioner BPI ~= inv(K), st = BPI*s0 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ st[m][j] = E->BPI[lev][m][j] * s0[m][j];
+
+
+ /* solve K*u1 = grad(st) for u1 */
+ assemble_grad_p(E, st, F, lev);
+ valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ if(!valid) fprintf(stderr, "not valid 2\n");
+ strip_bcs_from_residual(E, E->u1, lev);
+
+
+ /* t0 = div(rho_ref * u1) */
+ assemble_div_rho_u(E, E->u1, t0, lev);
+
+
+ /* omega = <t0, s0> / <t0, t0> */
+ omega = global_pdot(E, t0, s0, lev) / global_pdot(E, t0, t0, lev);
+
+
+ /* r2 = s0 - omega * t0 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ r2[m][j] = s0[m][j] - omega * t0[m][j];
+
+
+ /* P = P + alpha * pt + omega * st */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ s0[m][j] = alpha * pt[m][j] + omega * st[m][j];
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=1; j<=npno; j++)
+ P[m][j] += s0[m][j];
+
+
+ /* V = V - alpha * u0 - omega * u1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=0; j<neq; j++)
+ F[m][j] = alpha * u0[m][j] + omega * E->u1[m][j];
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(j=0; j<neq; j++)
+ V[m][j] -= F[m][j];
+
+
+ /* compute velocity and incompressibility residual */
+ assemble_div_rho_u(E, V, t0, lev);
+ incompressibility_residual(E, V, t0);
+
+ /* compute velocity and pressure corrections */
+ dpressure = sqrt( global_pdot(E, s0, s0, lev)
+ / (1.0e-32 + global_pdot(E, P, P, lev)) );
+ dvelocity = sqrt( global_vdot(E, F, F, lev)
+ / (1.0e-32 + E->monitor.vdotv) );
+
+
+ count++;
+
+ if(E->control.print_convergence && E->parallel.me==0) {
+ fprintf(E->fp, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ dvelocity, dpressure, E->monitor.solution_cycles);
+ fprintf(stderr, "AhatP (%03d) after %g seconds with div/v=%.3e "
+ "dv/v=%.3e and dp/p=%.3e for step %d\n",
+ count, CPU_time0()-time0, E->monitor.incompressibility,
+ dvelocity, dpressure, E->monitor.solution_cycles);
+ }
+
+
+ /* shift array pointers */
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ shuffle[m] = p1[m];
+ p1[m] = p2[m];
+ p2[m] = shuffle[m];
+
+ shuffle[m] = r1[m];
+ r1[m] = r2[m];
+ r2[m] = shuffle[m];
+ }
+
+ /* shift <r0, rt> = <r1, rt> */
+ r0dotrt = r1dotrt;
+
+ } /* end loop for conjugate gradient */
+
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ free((void *) r1[m]);
+ free((void *) r2[m]);
+ free((void *) pt[m]);
+ free((void *) p1[m]);
+ free((void *) p2[m]);
+ free((void *) rt[m]);
+ free((void *) v0[m]);
+ free((void *) s0[m]);
+ free((void *) st[m]);
+ free((void *) t0[m]);
+
+ free((void *) u0[m]);
+ }
+
+ *steps_max=count;
+
+ return(residual);
+
+}
+
+
+
+static double initial_vel_residual(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp)
+{
+ void assemble_del2_u();
+ void assemble_grad_p();
+ void strip_bcs_from_residual();
+ int solve_del2_u();
+ double global_vdot();
+
+ int neq = E->lmesh.neq;
+ int gneq = E->mesh.neq;
+ int lev = E->mesh.levmax;
+ int i, m;
+ double v_res;
+
+ v_res = sqrt(global_vdot(E, F, F, lev) / gneq);
+
+ if (E->parallel.me==0) {
+ fprintf(E->fp, "initial residue of momentum equation F %.9e %d\n",
+ v_res, gneq);
+ fprintf(stderr, "initial residue of momentum equation F %.9e %d\n",
+ v_res, gneq);
+ }
+
+
+ /* F = F - grad(P) - K*V */
+ assemble_grad_p(E, P, E->u1, lev);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ F[m][i] = F[m][i] - E->u1[m][i];
+
+ assemble_del2_u(E, V, E->u1, lev, 1);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ F[m][i] = F[m][i] - E->u1[m][i];
+
+ strip_bcs_from_residual(E, F, lev);
+
+
+ /* solve K*u1 = F for u1 */
+ solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ strip_bcs_from_residual(E, E->u1, lev);
+
+
+ /* V = V + u1 */
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ V[m][i] += E->u1[m][i];
+
+ return(v_res);
+}
+
+
+
+static double incompressibility_residual(struct All_variables *E,
+ double **V, double **r)
+{
+ double global_pdot();
+ double global_vdot();
+
+ int gnpno = E->mesh.npno;
+ int gneq = E->mesh.neq;
+ int lev = E->mesh.levmax;
+ double tmp1, tmp2;
+
+ /* incompressiblity residual = norm(F) / norm(V) */
+
+ tmp1 = global_vdot(E, V, V, lev);
+ tmp2 = global_pdot(E, r, r, lev);
+ E->monitor.incompressibility = sqrt((gneq / gnpno)
+ *( (1.0e-32 + tmp2)
+ / (1.0e-32 + tmp1) ));
+
+ E->monitor.vdotv = tmp1;
+
+ return(sqrt(tmp2/gnpno));;
+}
Modified: mc/3D/CitcomS/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -191,6 +191,7 @@
int i,j,e,node,m;
float VV[4][9],Vxyz[9][9],Szz,Sxx,Syy,Sxy,Sxz,Szy,div,vor;
+ double dilation[9];
double pre[9],tww[9],rtf[4][9];
double velo_scaling, stress_scaling;
@@ -218,8 +219,19 @@
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;
Vxyz[1][j] = 0.0;
Vxyz[2][j] = 0.0;
Vxyz[3][j] = 0.0;
@@ -262,16 +274,24 @@
+ VV[2][i]*E->N.vpt[GNVINDEX(i,j)]*cos(rtf[1][j])
- VV[1][i]*GNx.vpt[GNVXINDEX(1,i,j)] );
}
- Sxx += 2.0 * pre[j] * Vxyz[1][j];
- Syy += 2.0 * pre[j] * Vxyz[2][j];
- Szz += 2.0 * pre[j] * Vxyz[3][j];
- Sxy += pre[j] * Vxyz[4][j];
- Sxz += pre[j] * Vxyz[5][j];
- Szy += pre[j] * Vxyz[6][j];
- div += Vxyz[7][j]*dOmega.vpt[j];
- vor += Vxyz[8][j]*dOmega.vpt[j];
}
+ if(E->control.inv_gruneisen > 0) {
+ for(j=1;j<=vpts;j++)
+ dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
+ }
+
+ for(j=1;j<=vpts;j++) {
+ Sxx += 2.0 * pre[j] * (Vxyz[1][j] - dilation[j]);
+ Syy += 2.0 * pre[j] * (Vxyz[2][j] - dilation[j]);
+ Szz += 2.0 * pre[j] * (Vxyz[3][j] - dilation[j]);
+ Sxy += pre[j] * Vxyz[4][j];
+ Sxz += pre[j] * Vxyz[5][j];
+ Szy += pre[j] * Vxyz[6][j];
+ div += Vxyz[7][j]*dOmega.vpt[j];
+ vor += Vxyz[8][j]*dOmega.vpt[j];
+ }
+
Sxx /= E->eco[m][e].area;
Syy /= E->eco[m][e].area;
Szz /= E->eco[m][e].area;
@@ -281,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];
@@ -682,8 +703,8 @@
/* eu [m*dims+2] = VV[3][m+1]; */
/* } */
-/* get_elt_k(E,el,eltk,lev,j); */
-/* get_elt_k(E,elb,eltkb,lev,j); */
+/* get_elt_k(E,el,eltk,lev,j,1); */
+/* get_elt_k(E,elb,eltkb,lev,j,1); */
/* get_elt_f(E,el,eltf,0,j); */
/* get_elt_f(E,elb,eltfb,0,j); */
/* get_elt_g(E,el,eltg,lev,j); */
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -82,23 +82,23 @@
if (E->viscosity.SDEPV) {
input_float_vector("sdepv_expt",E->viscosity.num_mat,(E->viscosity.sdepv_expt),m);
}
-
+
input_boolean("PDEPV",&(E->viscosity.PDEPV),"off",m); /* plasticity addition by TWB */
if (E->viscosity.PDEPV) {
E->viscosity.pdepv_visited = 0;
-
+
input_boolean("pdepv_eff",&(E->viscosity.pdepv_eff),"on",m);
input_float_vector("pdepv_a",E->viscosity.num_mat,(E->viscosity.pdepv_a),m);
input_float_vector("pdepv_b",E->viscosity.num_mat,(E->viscosity.pdepv_b),m);
input_float_vector("pdepv_y",E->viscosity.num_mat,(E->viscosity.pdepv_y),m);
-
+
input_float("pdepv_offset",&(E->viscosity.pdepv_offset),"0.0",m);
}
if(E->viscosity.PDEPV || E->viscosity.SDEPV)
input_float("sdepv_misfit",&(E->viscosity.sdepv_misfit),"0.001",m);
-
+
input_boolean("CDEPV",&(E->viscosity.CDEPV),"off",m);
if(E->viscosity.CDEPV){ /* compositional viscosity */
if(!E->control.tracer)
@@ -191,16 +191,16 @@
if(E->viscosity.SDEPV)
visc_from_S(E,evisc,propogate);
-
+
if(E->viscosity.PDEPV) /* "plasticity" */
visc_from_P(E,evisc);
/* i think this should me placed differently i.e. before the
- stress dependence but I won't change it because it's by
+ stress dependence but I won't change it because it's by
someone else
- TWB
+ TWB
*/
if(E->viscosity.channel || E->viscosity.wedge)
apply_low_visc_wedge_channel(E, evisc);
@@ -420,7 +420,7 @@
EEta[m][ (i-1)*vpts + jj ] = tempa*
exp( (E->viscosity.E[l-1] + E->viscosity.Z[l-1]*zzz )
/ (E->viscosity.T[l-1]+temp) );
-
+
}
}
break;
@@ -478,21 +478,21 @@
case 6: /* eta = N_0 exp(E(T_0-T) + (1-z) Z_0 ) */
-
+
for(m=1;m <= E->sphere.caps_per_proc;m++)
for(i=1;i <= nel;i++) {
l = E->mat[m][i];
- if(E->control.mat_control)
+ if(E->control.mat_control)
tempa = E->viscosity.N0[l-1] * E->VIP[m][i];
else
tempa = E->viscosity.N0[l-1];
j = 0;
-
+
for(kk=1;kk<=ends;kk++) {
TT[kk] = E->T[m][E->ien[m][i].node[kk]];
zz[kk] = (1.0 - E->sx[m][3][E->ien[m][i].node[kk]]);
}
-
+
for(jj=1;jj <= vpts;jj++) {
temp=0.0;zzz=0.0;
for(kk=1;kk <= ends;kk++) {
@@ -508,10 +508,10 @@
}
}
break;
-
+
}
return;
@@ -551,35 +551,35 @@
return;
}
-void visc_from_P(E,EEta) /* "plasticity" implementation
-
- viscosity will be limited by a yield stress
-
+void visc_from_P(E,EEta) /* "plasticity" implementation
+
+ viscosity will be limited by a yield stress
+
\sigma_y = min(a + b * (1-r), y)
-
+
where a,b,y are parameters input via pdepv_a,b,y
-
- and
-
- \eta_y = \sigma_y / (2 \eps_II)
-
+
+ and
+
+ \eta_y = \sigma_y / (2 \eps_II)
+
where \eps_II is the second invariant. Then
-
+
\eta_eff = (\eta_0 \eta_y)/(\eta_0 + \eta_y)
-
+
for pdepv_eff = 1
-
- or
+ or
+
\eta_eff = min(\eta_0,\eta_y)
-
+
for pdepv_eff = 0
-
+
where \eta_0 is the regular viscosity
-
-
+
+
TWB
-
+
*/
struct All_variables *E;
float **EEta;
@@ -596,14 +596,14 @@
eedot = (float *) malloc((2+nel)*sizeof(float));
for(m=1;m<=E->sphere.caps_per_proc;m++) {
-
+
if(E->viscosity.pdepv_visited){
strain_rate_2_inv(E,m,eedot,1); /* get second invariant for all elements */
}else{
for(e=1;e<=nel;e++) /* initialize with unity if no velocities around */
- eedot[e] = 1.0e-5;
+ eedot[e] = 1.0e-5;
if(m == E->sphere.caps_per_proc)
E->viscosity.pdepv_visited = 1;
}
@@ -614,13 +614,13 @@
for(kk=1;kk <= ends;kk++) /* nodal depths */
zz[kk] = (1.0 - E->sx[m][3][E->ien[m][e].node[kk]]); /* for depth, zz = 1 - r */
-
+
for(jj=1;jj <= vpts;jj++){ /* loop through integration points */
zzz = 0.0; /* get mean depth of integration point */
- for(kk=1;kk<=ends;kk++)
+ for(kk=1;kk<=ends;kk++)
zzz += zz[kk] * E->N.vpt[GNVINDEX(kk,jj)];
-
+
/* depth dependent yield stress */
tau = E->viscosity.pdepv_a[l] + zzz * E->viscosity.pdepv_b[l];
@@ -629,12 +629,12 @@
/* yield viscosity */
eta_p = tau/(2.0 * eedot[e] + 1e-7) + E->viscosity.pdepv_offset;
-
- if(E->viscosity.pdepv_eff){
+
+ if(E->viscosity.pdepv_eff){
/* two dashpots in series */
eta_new = 1.0/(1.0/EEta[m][ (e-1)*vpts + jj ] + 1.0/eta_p);
- }else{
+ }else{
/* min viscosities*/
eta_new = min(EEta[m][ (e-1)*vpts + jj ], eta_p);
}
@@ -648,7 +648,7 @@
return;
}
-/*
+/*
multiply with compositional factor which is determined by a geometric
mean average from the tracer composition, assuming two flavors and
@@ -662,7 +662,7 @@
float comp,comp_fac,CC[9],tcomp;
double vmean,cc_loc;
int m,l,z,jj,kk,i;
-
+
const int vpts = vpoints[E->mesh.nsd];
const int nel = E->lmesh.nel;
const int ends = enodes[E->mesh.nsd];
@@ -684,9 +684,9 @@
cc_loc = 0.0;
for(kk = 1; kk <= ends; kk++)
cc_loc += CC[kk] * E->N.vpt[GNVINDEX(kk, jj)];
-
+
/* geometric mean of viscosity */
- vmean = exp(cc_loc * E->viscosity.cdepv_ff[1] +
+ vmean = exp(cc_loc * E->viscosity.cdepv_ff[1] +
(1.0-cc_loc) * E->viscosity.cdepv_ff[0]);
/* multiply the viscosity with this prefactor */
EEta[m][ (i-1)*vpts + jj ] *= vmean;
@@ -778,14 +778,21 @@
}
}
+ 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];
- edot[2][2] += 2.0 * Vxyz[2][j];
- edot[3][3] += 2.0 * Vxyz[3][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];
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-16 00:27:43 UTC (rev 7831)
@@ -268,7 +268,7 @@
double *angle1[MAX_LEVELS][NCS][5];
double dircos[4][4];
- double *R[MAX_LEVELS],*R_redundant;
+ double *R[MAX_LEVELS];
double ro,ri;
struct CAP cap[NCS];
@@ -289,33 +289,24 @@
int NEQ[MAX_LEVELS];
int NNO[MAX_LEVELS];
int NNOV[MAX_LEVELS];
- int NLNO[MAX_LEVELS];
int NPNO[MAX_LEVELS];
int NEL[MAX_LEVELS];
int NOX[MAX_LEVELS];
int NOZ[MAX_LEVELS];
int NOY[MAX_LEVELS];
- int NMX[MAX_LEVELS];
int ELX[MAX_LEVELS];
int ELZ[MAX_LEVELS];
int ELY[MAX_LEVELS];
- int LNDS[MAX_LEVELS];
- int LELS[MAX_LEVELS];
int SNEL[MAX_LEVELS];
- int neqd;
int neq;
int nno;
int nnov;
- int nlno;
int npno;
int nel;
int snel;
int elx;
int elz;
int ely;
- int nnx[4]; /* general form of ... */
- int gnox;
- int gelx;
int nox;
int noz;
int noy;
@@ -331,7 +322,6 @@
int NXS[MAX_LEVELS];
int NYS[MAX_LEVELS];
int NZS[MAX_LEVELS];
- int nmx;
int nsf; /* nodes for surface observables */
int toptbc;
int bottbc;
@@ -339,26 +329,11 @@
int botvbc;
int sidevbc;
- char topvbc_file[100];
- char botvbc_file[100];
- char sidevbc_file[100];
- char gridfile[4][100];
-
int periodic_x;
int periodic_y;
float layer[4]; /* dimensionless dimensions */
- float lidz;
- float bl1width[4],bl2width[4],bl1mag[4],bl2mag[4];
- float hwidth[4],magnitude[4],offset[4],width[4]; /* grid compression information */
double volume;
- int fnodal_malloc_size;
- int dnodal_malloc_size;
- int feqn_malloc_size;
- int deqn_malloc_size;
- int bandwidth;
- int null_source;
- int null_sink;
int matrix_size[MAX_LEVELS];
} ;
@@ -385,53 +360,20 @@
struct MONITOR {
- char node_output[100][6]; /* recording the format of the output data */
- char sobs_output[100][6]; /* recording the format of the output data */
- int node_output_cols;
- int sobs_output_cols;
-
int solution_cycles;
int solution_cycles_init;
int stop_topo_loop;
int topo_loop;
- float time_scale;
- float length_scale;
- float viscosity_scale;
- float geoscale;
- float tpgscale;
- float grvscale;
-
- float delta_v_last_soln;
float elapsed_time;
- float elapsed_time_vsoln;
- float elapsed_time_vsoln1;
- float reference_stress;
float incompressibility;
float vdotv;
- float nond_av_heat_fl;
- float nond_av_adv_hfl;
double cpu_time_at_start;
double cpu_time_at_last_cycle;
- float tpgkmag;
- float grvkmag;
- float Nusselt;
- float Vmax;
- float Vsrms;
- float Vrms;
- float Vrms_surface;
- float Vrms_base;
- float F_surface;
- float F_base;
- float Frat_surface;
- float Frat_base;
float T_interior;
float T_maxvaried;
- float Sigma_max;
- float Sigma_interior;
- float Vi_average;
};
@@ -441,8 +383,6 @@
char output_written_external_command[500]; /* a unix command to run when output files have been created */
int ORTHO,ORTHOZ; /* indicates levels of mesh symmetry */
- char B_is_good[MAX_LEVELS]; /* general information controlling program flow */
- char Ahat_is_good[MAX_LEVELS]; /* general information controlling program flow */
char data_prefix[50];
char data_prefix_old[50];
@@ -453,22 +393,13 @@
char data_file[200];
char old_P_file[200];
- char post_topo_file[100];
- char slabgeoid_file[100];
-
- char which_data_files[1000];
- char which_horiz_averages[1000];
- char which_running_data[1000];
- char which_observable_data[1000];
-
char PROBLEM_TYPE[20]; /* one of ... */
- int KERNEL;
int CONVECTION;
int stokes;
int restart;
int post_p;
int post_topo;
- int SLAB;
+ int compressible;
char GEOMETRY[20]; /* one of ... */
int CART2D;
@@ -489,10 +420,6 @@
int AVS;
int CONMAN;
- int read_density;
- int read_slab;
- int read_slabgeoid;
-
int pseudo_free_surf;
int tracer;
@@ -517,7 +444,13 @@
float sob_tolerance;
+ /* Rayleigh # */
float Atemp;
+ /* Dissipation # */
+ float disptn_number;
+ /* inverse of Gruneisen parameter */
+ float inv_gruneisen;
+
float inputdiff;
float VBXtopval;
float VBXbotval;
@@ -545,26 +478,15 @@
float Q0;
float Q0ER;
- float jrelax;
-
int precondition;
- int vprecondition;
int keep_going;
int v_steps_low;
int v_steps_high;
int v_steps_upper;
- int max_vel_iterations;
int p_iterations;
- int max_same_visc;
- float max_res_red_each_p_mg;
- float sub_stepping_factor;
int mg_cycle;
- int true_vcycle;
int down_heavy;
int up_heavy;
- int depth_dominated;
- int eqn_viscosity;
- int eqn_zigzag;
int verbose;
int side_sbcs;
@@ -585,15 +507,20 @@
int print_convergence;
int sdepv_print_convergence;
- /* modules */
- int MELTING_MODULE;
- int CHEMISTRY_MODULE;
+};
- /* for embedded setting */
- int embedded;
+struct REF_STATE {
+ double *rho;
+ double *expansivity;
+ double *capacity;
+ double *conductivity;
+ double *gravity;
+ double *Tadi;
+ double *dlnrhodr;
};
+
struct DATA {
float layer_km;
float radius_km;
@@ -742,8 +669,6 @@
struct SIEN *sien[NCS];
struct ID *id[NCS];
struct COORD *ECO[MAX_LEVELS][NCS];
- struct IEN *IEN_redundant[NCS];
- struct ID *ID_redundant[NCS];
struct IEN *IEN[MAX_LEVELS][NCS]; /* global at each level */
struct FNODE *TWW[MAX_LEVELS][NCS]; /* for nodal averages */
struct ID *ID[MAX_LEVELS][NCS];
@@ -758,6 +683,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];
@@ -765,18 +693,22 @@
double *BI[MAX_LEVELS][NCS],*BPI[MAX_LEVELS][NCS];
+ double *rho;
+ 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];
double *u1[NCS];
double *temp[NCS],*temp1[NCS];
- float *NP[NCS],*edot[NCS],*Mass[NCS],*tw[NCS];
+ float *NP[NCS],*edot[NCS],*Mass[NCS];
float *MASS[MAX_LEVELS][NCS];
- double *ZZ;
+ double *TMass[NCS];
double *SX[MAX_LEVELS][NCS][4],*X[MAX_LEVELS][NCS][4];
double *sx[NCS][4],*x[NCS][4];
double *surf_det[NCS][5];
double *SinCos[MAX_LEVELS][NCS][4];
- float *TT;
float *V[NCS][4],*GV[NCS][4],*GV1[NCS][4];
float *stress[NCS];
@@ -786,7 +718,6 @@
float *Vi[NCS],*EVi[NCS];
float *VI[MAX_LEVELS][NCS],*EVI[MAX_LEVELS][NCS];
- float *TW[MAX_LEVELS][NCS]; /* nodal weightings */
int num_zero_resid[MAX_LEVELS][NCS];
int *zero_resid[MAX_LEVELS][NCS];
Modified: mc/3D/CitcomS/trunk/module/mesher.c
===================================================================
--- mc/3D/CitcomS/trunk/module/mesher.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/module/mesher.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -32,6 +32,7 @@
#include "global_defs.h"
#include "parallel_related.h"
+#include "material_properties.h"
extern void initial_mesh_solver_setup(struct All_variables *);
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-16 00:03:05 UTC (rev 7830)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-16 00:27:43 UTC (rev 7831)
@@ -428,6 +428,7 @@
PyObject *obj, *properties, *out;
struct All_variables *E;
FILE *fp;
+ float tmp;
if (!PyArg_ParseTuple(args, "OOO:Solver_set_properties",
&obj, &properties, &out))
@@ -444,6 +445,11 @@
getStringProperty(properties, "datafile_old", E->control.data_prefix_old, fp);
getFloatProperty(properties, "rayleigh", E->control.Atemp, fp);
+ getFloatProperty(properties, "dissipation_number", E->control.disptn_number, 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