[cig-commits] r7798 - in mc/3D/CitcomS/branches/compressible:
CitcomS/Components/Stokes_solver lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Aug 9 15:57:29 PDT 2007
Author: tan2
Date: 2007-08-09 15:57:28 -0700 (Thu, 09 Aug 2007)
New Revision: 7798
Modified:
mc/3D/CitcomS/branches/compressible/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.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_h5.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
mc/3D/CitcomS/branches/compressible/module/setProperties.c
Log:
Finished the compressible Stokes solver for TALA.
Two non-dimensional parameters are added: "dissipation_number" and "gruneisen"
under the Solver component. One can use the original incompressible solver by
setting "gruneisen=0". The code will treat this as "gruneisen=infinity".
Setting non-zero value to "gruneisen" will switch to compressible solver.
One can use the TALA solver for incompressible case by setting "gruneisen" to
a non-zero value while setting "dissipation_number=0". This is useful when
debugging the compressible solver.
Two implementations are available: one by Wei Leng (U. Colorado) and one by
Eh Tan (CIG). Leng's version uses the original conjugate gradient method for
the Uzawa iteration and moves the contribution of compressibility to the RHS,
similar to the method of Ita and King, JGR, 1994. Tan's version uses the
bi-conjugate gradient stablized method for the Uzawa iteration, similar to the
method of Tan and Gurnis, JGR, 2007. Both versions agree very well. In the
benchmark case, 33x33x33 nodes per cap, Di/gamma=1.0, Ra=1.0, delta function
of load at the mid mantle, the peak velocity differs by only 0.007%. Leng's
version is enabled by default. Edit function solve_Ahat_p_fhat() in
lib/Stokes_flow_Incomp.c to switch to Tan's version.
Modified: mc/3D/CitcomS/branches/compressible/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/branches/compressible/CitcomS/Components/Stokes_solver/Incompressible.py 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/CitcomS/Components/Stokes_solver/Incompressible.py 2007-08-09 22:57:28 UTC (rev 7798)
@@ -102,6 +102,9 @@
aug_lagr = prop.bool("aug_lagr", default=True)
aug_number = prop.float("aug_number", default=2.0e3)
+ compress_iter_maxstep = prop.int("compress_iter_maxstep", default=100)
+ relative_err_accuracy = prop.float("relative_err_accuracy", default=0.01)
+
# version
__id__ = "$Id$"
Modified: mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Construct_arrays.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -681,6 +681,30 @@
}
+/*==============================================
+ For compressible cases, construct c matrix,
+ where c = \frac{d rho_r}{dr} / rho_r * u_r
+ ==============================================*/
+
+void construct_elt_cs(struct All_variables *E)
+{
+ int m, el, lev;
+ void get_elt_c();
+
+/* if(E->control.verbose && E->parallel.me==0) */
+/* fprintf(stderr,"storing elt c matrices\n"); */
+
+ for(lev=E->mesh.gridmin;lev<=E->mesh.gridmax;lev++)
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(el=1;el<=E->lmesh.NEL[lev];el++) {
+ get_elt_c(E,el,E->elt_c[lev][m][el].c,lev,m);
+ }
+
+
+ return;
+}
+
+
/* ==============================================================
routine for constructing stiffness and node_maps
============================================================== */
Modified: mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Element_calculations.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -168,61 +168,62 @@
const int vpts = VPOINTS3D;
const int ends = ENODES3D;
- double ra[9], si[9], ct[9];
+ double ra[9], isi[9], ct[9];
+ double gnx0, gnx1, gnx2, shp, cc1, cc2, cc3;
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];
+ isi[k] = 1.0 / sin(rtf[1][k]);
+ ct[k] = cos(rtf[1][k]) * isi[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];
+ for(a=1;a<=ends;a++)
+ for(k=1;k<=vpts;k++) {
+ gnx0 = GNx->vpt[GNVXINDEX(0,a,k)];
+ gnx1 = GNx->vpt[GNVXINDEX(1,a,k)];
+ gnx2 = GNx->vpt[GNVXINDEX(2,a,k)];
+ shp = N->vpt[GNVINDEX(a,k)];
+ cc1 = cc->vpt[BVINDEX(1,n,a,k)];
+ cc2 = cc->vpt[BVINDEX(2,n,a,k)];
+ cc3 = cc->vpt[BVINDEX(3,n,a,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][1] = ( gnx0 * cc1
+ + shp * ccx->vpt[BVXINDEX(1,n,1,a,k)]
+ + shp * cc3 ) * ra[k];
- ba[a][k][n][3]=
- GNx->vpt[GNVXINDEX(2,a,k)]*cc->vpt[BVINDEX(3,n,a,k)];
+ ba[a][k][n][2] = ( shp * cc1 * ct[k]
+ + shp * cc3
+ + ( gnx1 * cc2
+ + shp * ccx->vpt[BVXINDEX(2,n,2,a,k)] )
+ * isi[k] ) * ra[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][3] = gnx2 * cc3;
- 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][4] = ( gnx0 * cc2
+ + shp * ccx->vpt[BVXINDEX(2,n,1,a,k)]
+ - shp * cc2 * ct[k]
+ + ( gnx1 * cc1
+ + shp * ccx->vpt[BVXINDEX(1,n,2,a,k)] )
+ * isi[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];
+ ba[a][k][n][5] = gnx2 * cc1
+ + ( gnx0 * cc3
+ + shp * ( ccx->vpt[BVXINDEX(3,n,1,a,k)]
+ - cc1 ) ) * ra[k];
+
+ ba[a][k][n][6] = gnx2 * cc2
+ - ra[k] * shp * cc2
+ + ( gnx1 * cc3
+ + shp * ccx->vpt[BVXINDEX(3,n,2,a,k)] )
+ * isi[k] * ra[k];
}
- return;
+ return;
}
+
+
/*==============================================================
Function to supply the element k matrix for a given element e.
============================================================== */
@@ -238,6 +239,10 @@
int a,b,i,j,i1,j1,k;
double rtf[4][9],W[9];
+
+ const double two = 2.0;
+ const double two_thirds = 2.0/3.0;
+
void get_global_shape_fn();
void construct_c3x3matrix_el();
struct Shape_function GN;
@@ -247,22 +252,6 @@
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 = VPOINTS3D;
@@ -291,23 +280,23 @@
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 < 1e-6)
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];
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] += W[k] * ( two * ( ba[a][k][i][1]*ba[b][k][j][1] +
+ ba[a][k][i][2]*ba[b][k][j][2] +
+ ba[a][k][i][3]*ba[b][k][j][3] ) +
+ ba[a][k][i][4]*ba[b][k][j][4] +
+ ba[a][k][i][5]*ba[b][k][j][5] +
+ ba[a][k][i][6]*ba[b][k][j][6] );
- 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]*d2[i1][j1]*ba[a][k][i][i1]*ba[b][k][j][j1];
+ if(E->control.inv_gruneisen != 0)
+ for(i=1;i<=dims;i++)
+ for(j=1;j<=dims;j++)
+ for(k=1;k<=VPOINTS3D;k++)
+ bdbmu[i][j] -= W[k] * two_thirds *
+ ( ba[a][k][i][1] + ba[a][k][i][2] + ba[a][k][i][3] ) *
+ ( ba[b][k][j][1] + ba[b][k][j][2] + ba[b][k][j][3] );
/**/
@@ -578,28 +567,32 @@
}
+/* =====================================================
+ Assemble grad(rho_ref*ez)*V element by element.
+ Note that the storage is not zero'd before assembling.
+ ===================================================== */
-
-/* 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)
+void assemble_c_u(struct All_variables *E,
+ double **U, double **result, int level)
{
- int e, nz, j3, a, b, m;
+ int e,j1,j2,j3,p,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;
+ const int npno = E->lmesh.NPNO[level];
- 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++) {
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(a=1;a<=ends;a++) {
+ p = (a-1)*dims;
+ for(e=1;e<=nel;e++) {
b = E->IEN[level][m][e].node[a];
- j3 = E->ID[level][m][b].doff[3];
- result[m][e] += dlnrhodr[e] * U[m][j3];
+ j1= E->ID[level][m][b].doff[1];
+ j2= E->ID[level][m][b].doff[2];
+ j3= E->ID[level][m][b].doff[3];
+ result[m][e] += E->elt_c[level][m][e].c[p ][0] * U[m][j1]
+ + E->elt_c[level][m][e].c[p+1][0] * U[m][j2]
+ + E->elt_c[level][m][e].c[p+2][0] * U[m][j3];
}
}
@@ -608,13 +601,17 @@
+/* =====================================================
+ Assemble div(rho_ref*V) = div(V) + grad(rho_ref*ez)*V
+ element by element
+ ===================================================== */
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);
+ assemble_c_u(E, U, result, level);
return;
}
@@ -761,6 +758,56 @@
/*==============================================================
+ Function to supply the element c matrix for a given element e.
+ ============================================================== */
+
+void get_elt_c(struct All_variables *E, int el,
+ higher_precision elt_c[24][1], int lev, int m)
+{
+ void get_global_shape_fn();
+ void construct_c3x3matrix_el();
+ int p, a, i;
+ double temp, beta, x[4];
+
+ struct Shape_function GN;
+ struct Shape_function_dx GNx;
+ struct Shape_function_dA dOmega;
+ double rtf[4][9];
+
+ const int dims = E->mesh.nsd;
+ const int ends = enodes[dims];
+ const int sphere_key = 1;
+
+ get_global_shape_fn(E,el,&GN,&GNx,&dOmega,2,sphere_key,rtf,lev,m);
+
+ if ((el-1)%E->lmesh.ELZ[lev]==0)
+ construct_c3x3matrix_el(E,el,&E->element_Cc,&E->element_Ccx,lev,m,1);
+
+ temp = p_point[1].weight[dims-1] * dOmega.ppt[1];
+ beta = E->control.disptn_number * E->control.inv_gruneisen;
+
+ for(a=1;a<=ends;a++) {
+ for (i=1;i<=dims;i++) {
+#if 1
+ /* hard coded dln(rho)/dr here */
+
+ x[i] = - beta * E->N.ppt[GNPINDEX(a,1)]
+ * E->element_Cc.ppt[BPINDEX(3,i,a,1)];
+#else
+ /* compute dln(rho)/dr from rho(r) here */
+ /* XXX */
+#endif
+ }
+ p=dims*(a-1);
+ elt_c[p ][0] = -x[1] * temp;
+ elt_c[p+1][0] = -x[2] * temp;
+ elt_c[p+2][0] = -x[3] * temp;
+ }
+ return;
+}
+
+
+/*==============================================================
Function to supply the element g matrix for a given element e.
============================================================== */
Modified: mc/3D/CitcomS/branches/compressible/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Instructions.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -56,6 +56,7 @@
void allocate_velocity_vars(struct All_variables*);
void check_bc_consistency(struct All_variables*);
void construct_elt_gs(struct All_variables*);
+void construct_elt_cs(struct All_variables*);
void construct_id(struct All_variables*);
void construct_ien(struct All_variables*);
void construct_lm(struct All_variables*);
@@ -125,6 +126,8 @@
construct_sub_element(E);
construct_shape_functions(E);
construct_elt_gs(E);
+ if(E->control.inv_gruneisen != 0)
+ construct_elt_cs(E);
reference_state(E);
@@ -407,11 +410,15 @@
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(fabs(tmp) > 1e-6)
+ /* special case: if tmp==0, set gruneisen as inf */
+ if(tmp != 0)
E->control.inv_gruneisen = 1/tmp;
else
E->control.inv_gruneisen = 0;
+ input_int("compress_iter_maxstep",&(E->control.compress_iter_maxstep),"100",m);
+ input_float("relative_err_accuracy",&(E->control.relative_err_accuracy),"0.01",m);
+
/* data section */
input_float("Q0",&(E->control.Q0),"0.0",m);
input_float("gravacc",&(E->data.grav_acc),"9.81",m);
@@ -595,11 +602,14 @@
E->ELEMENT[i][j][k] = 0;
/*ccccc*/
- E->elt_del[i][j]=(struct EG *) malloc((nel+1)*sizeof(struct EG));
+ E->elt_del[i][j] = (struct EG *) malloc((nel+1)*sizeof(struct EG));
- E->EVI[i][j] = (float *) malloc((nel+2)*vpoints[E->mesh.nsd]*sizeof(float));
- E->BPI[i][j] = (double *) malloc((npno+1)*sizeof(double));
+ if(E->control.inv_gruneisen != 0)
+ E->elt_c[i][j] = (struct EC *) malloc((nel+1)*sizeof(struct EC));
+ E->EVI[i][j] = (float *) malloc((nel+2)*vpoints[E->mesh.nsd]*sizeof(float));
+ E->BPI[i][j] = (double *) malloc((npno+1)*sizeof(double));
+
E->ID[i][j] = (struct ID *) malloc((nno+2)*sizeof(struct ID));
E->VI[i][j] = (float *) malloc((nno+2)*sizeof(float));
E->NODE[i][j] = (unsigned int *)malloc((nno+2)*sizeof(unsigned int));
Modified: mc/3D/CitcomS/branches/compressible/lib/Material_properties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Material_properties.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -60,8 +60,6 @@
/* reference profile of temperature */
E->refstate.Tadi = (double *) malloc((noz+1)*sizeof(double));
- /* reference profile of d(ln(rho_ref))/dr */
- E->refstate.dlnrhodr = (double *) malloc((nel+1)*sizeof(double));
}
@@ -70,19 +68,19 @@
int noz = E->lmesh.noz;
int nel = E->lmesh.nel;
int i;
- double r, z, tmp, T0;
+ double r, z, beta, T0;
- tmp = E->control.disptn_number * E->control.inv_gruneisen;
+ beta = E->control.disptn_number * E->control.inv_gruneisen;
T0 = E->data.surf_temp / E->data.ref_temperature;
if(E->parallel.me == 0)
- fprintf(stderr, "Di=%f, gamma=%f, surf_temp=%f, dT=%f\n",
- E->control.disptn_number, 1/E->control.inv_gruneisen,
+ fprintf(stderr, "Di=%e, gamma=%e, surf_temp=%e, dT=%e\n",
+ E->control.disptn_number, 1.0/E->control.inv_gruneisen,
E->data.surf_temp, E->data.ref_temperature);
for(i=1; i<=noz; i++) {
r = E->sx[1][3][i];
z = 1 - r;
- E->refstate.rho[i] = exp(tmp*z);
+ E->refstate.rho[i] = exp(beta*z);
E->refstate.expansivity[i] = 1;
E->refstate.capacity[i] = 1;
E->refstate.conductivity[i] = 1;
@@ -90,49 +88,20 @@
E->refstate.Tadi[i] = T0 * (exp(E->control.disptn_number * z) - 1);
}
- if (tmp) {
- for(i=1; i<=nel; i++) {
- // TODO: dln(rho)/dr
- //E->refstate.dlnrhodr[i] = - tmp;
- int a, p, j3;
- int m = 1;
- int lev = E->mesh.levmax;
- int dims = E->mesh.nsd;
- int nz = ((i-1) % E->lmesh.elz) + 1;
- double rho_mean = 0.5 * (E->refstate.rho[nz] + E->refstate.rho[nz+1]);
- double rho[9], tmp;
-
- rho[1] = rho[2] = rho[3] = rho[4] = E->sx[1][3][nz]*E->refstate.rho[nz];
- rho[5] = rho[6] = rho[7] = rho[8] = E->sx[1][3][nz+1]*E->refstate.rho[nz+1];
-
- tmp = 0.0;
- for(a=1; a<=8; a++) {
- p = (a-1)*dims;
- tmp += E->elt_del[lev][m][i].g[p+2][0] * rho[a];
- }
- E->refstate.dlnrhodr[i] = tmp / rho_mean;
-
- }
+ if(E->parallel.me == 0) {
+ fprintf(stderr, "nz radius depth rho Tadi\n");
+ fprintf(E->fp, "nz radius depth rho Tadi\n");
}
- else {
- for(i=1; i<=nel; i++) E->refstate.dlnrhodr[i] = 0;
- }
-
-
if(E->parallel.me < E->parallel.nprocz)
for(i=1; i<=noz; i++) {
- fprintf(stderr, "%d %f %f %f %f\n",
+ fprintf(stderr, "%d %f %f %e %e\n",
i+E->lmesh.nzs-1, E->sx[1][3][i], 1-E->sx[1][3][i],
E->refstate.rho[i], E->refstate.Tadi[i]);
+ fprintf(E->fp, "%d %f %f %e %e\n",
+ i+E->lmesh.nzs-1, E->sx[1][3][i], 1-E->sx[1][3][i],
+ E->refstate.rho[i], E->refstate.Tadi[i]);
}
- if(E->parallel.me < E->parallel.nprocz)
- for(i=1; i<=E->lmesh.elz; i++) {
- fprintf(stderr, "%d %f\n",
- i+E->lmesh.nzs-1, E->refstate.dlnrhodr[i]);
- }
-
- //parallel_process_termination();
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Output_h5.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -1461,11 +1461,9 @@
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
- )
- );
+ (E->control.inv_gruneisen == 0)?
+ 1.0/E->control.inv_gruneisen :
+ E->control.inv_gruneisen));
status = set_attribute_float(input, "Q0", E->control.Q0);
status = set_attribute_int(input, "stokes_flow_only", E->control.stokes);
Modified: mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -44,6 +44,9 @@
static float solve_Ahat_p_fhat_BiCG(struct All_variables *E,
double **V, double **P, double **F,
double imp, int *steps_max);
+static float solve_Ahat_p_fhat_iterCG(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);
@@ -104,16 +107,18 @@
{
float residual;
- if(E->control.inv_gruneisen < 1e-6)
+ if(E->control.inv_gruneisen == 0)
residual = solve_Ahat_p_fhat_CG(E, V, P, F, imp, steps_max);
- else
- residual = solve_Ahat_p_fhat_BiCG(E, V, P, F, imp, steps_max);
+ else {
+ //residual = solve_Ahat_p_fhat_BiCG(E, V, P, F, imp, steps_max);
+ residual = solve_Ahat_p_fhat_iterCG(E, V, P, F, imp, steps_max);
+ }
return(residual);
}
-/* Solve incompressible Stokes flow (Boussinesq Approximation) using
+/* Solve incompressible Stokes flow using
* conjugate gradient (CG) iterations
*/
@@ -121,7 +126,7 @@
double **V, double **P, double **F,
double imp, int *steps_max)
{
- int m, i, j, count, valid, lev, npno, neq;
+ int m, j, count, valid, lev, npno, neq;
int gnpno, gneq;
double *r1[NCS], *r2[NCS], *z1[NCS], *s1[NCS], *s2[NCS];
@@ -134,6 +139,7 @@
double time0, CPU_time0();
float dpressure, dvelocity;
+ void assemble_c_u();
void assemble_div_u();
void assemble_del2_u();
void assemble_grad_p();
@@ -158,12 +164,27 @@
time0 = CPU_time0();
count = 0;
+ if(E->control.inv_gruneisen != 0) {
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(j=1;j<=npno;j++)
+ r2[m][j] = 0.0;
+
+ assemble_c_u(E, V, r2, lev);
+ }
+
/* calculate the initial velocity residual */
v_res = initial_vel_residual(E, V, P, F, imp);
- /* initial residual r1 = div(V) */
+ /* initial residual r1 = div(V) or r1 = div(rho*V) if compressible */
assemble_div_u(E, V, r1, lev);
+
+ if(E->control.inv_gruneisen != 0)
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(j=1;j<=npno;j++) {
+ r1[m][j] += r2[m][j];
+ }
+
residual = incompressibility_residual(E, V, r1);
if (E->control.print_convergence && E->parallel.me==0) {
@@ -305,15 +326,15 @@
return(residual);
}
-/* Solve incompressible Stokes flow (Boussinesq Approximation) using
- * bi-conjugate gradient stablized (BiCG-stab)iterations
+/* Solve compressible Stokes flow using
+ * bi-conjugate gradient stablized (BiCG-stab) iterations
*/
static float solve_Ahat_p_fhat_BiCG(struct All_variables *E,
double **V, double **P, double **F,
double imp, int *steps_max)
{
- void assemble_div_u();
+ void assemble_div_rho_u();
void assemble_del2_u();
void assemble_grad_p();
void strip_bcs_from_residual();
@@ -325,7 +346,7 @@
int gnpno, gneq;
int npno, neq;
- int m, i, j, count, lev;
+ int m, j, count, lev;
int valid;
double alpha, beta, omega;
@@ -576,7 +597,87 @@
}
+/* Solve compressible Stokes flow using
+ * conjugate gradient (CG) iterations with an outer iteration
+ */
+static float solve_Ahat_p_fhat_iterCG(struct All_variables *E,
+ double **V, double **P, double **FF,
+ double imp, int *steps_max)
+{
+ int m, i;
+ int cycles, num_of_loop;
+ double residual;
+ double relative_err_v, relative_err_p;
+ double *F[NCS],*old_v[NCS], *old_p[NCS],*diff_v[NCS],*diff_p[NCS];
+
+ const int npno = E->lmesh.npno;
+ const int neq = E->lmesh.neq;
+ const int lev = E->mesh.levmax;
+
+ double global_vdot(),global_pdot();
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ F[m] = (double *)malloc((neq+1)*sizeof(double));
+ old_v[m] = (double *)malloc((neq+1)*sizeof(double));
+ diff_v[m] = (double *)malloc((neq+1)*sizeof(double));
+ old_p[m] = (double *)malloc((npno+1)*sizeof(double));
+ diff_p[m] = (double *)malloc((npno+1)*sizeof(double));
+ }
+
+ cycles = E->control.p_iterations;
+
+ residual = 1.0;
+ relative_err_v = 1.0;
+ relative_err_p = 1.0;
+ num_of_loop = 0;
+
+ while((relative_err_v >= E->control.relative_err_accuracy ||
+ relative_err_p >= E->control.relative_err_accuracy) &&
+ num_of_loop <= E->control.compress_iter_maxstep) {
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ /* copy the original force vector since we need to keep it intact between iterations */
+ for(i=0;i<neq;i++) F[m][i] = FF[m][i];
+ for(i=0;i<neq;i++) old_v[m][i] = V[m][i];
+ for(i=1;i<=npno;i++) old_p[m][i] = P[m][i];
+ }
+
+ residual = solve_Ahat_p_fhat_CG(E,V,P,F,E->control.accuracy,&cycles);
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=0;i<neq;i++) diff_v[m][i] = V[m][i] - old_v[m][i];
+
+ relative_err_v = sqrt( global_vdot(E,diff_v,diff_v,lev) /
+ (1.0e-32 + global_vdot(E,V,V,lev)) );
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(i=1;i<=npno;i++) diff_p[m][i] = P[m][i] - old_p[m][i];
+
+ relative_err_p = sqrt( global_pdot(E,diff_p,diff_p,lev) /
+ (1.0e-32 + global_pdot(E,P,P,lev)) );
+
+ num_of_loop++;
+
+ if(E->parallel.me == 0) {
+ fprintf(stderr, "Relative error err_v / v = %e and err_p / p = %e after %d loops\n\n", relative_err_v, relative_err_p, num_of_loop);
+ fprintf(E->fp, "Relative error err_v / v = %e and err_p / p = %e after %d loops\n\n", relative_err_v, relative_err_p, num_of_loop);
+ }
+
+ } /* end of while */
+
+ for (m=1;m<=E->sphere.caps_per_proc;m++) {
+ free((void *) F[m]);
+ free((void *) old_v[m]);
+ free((void *) old_p[m]);
+ free((void *) diff_v[m]);
+ free((void *) diff_p[m]);
+ }
+
+ return(residual);
+}
+
+
static double initial_vel_residual(struct All_variables *E,
double **V, double **P, double **F,
double imp)
Modified: mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Topo_gravity.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -274,7 +274,7 @@
}
}
- if(E->control.inv_gruneisen > 0) {
+ if(E->control.inv_gruneisen != 0) {
for(j=1;j<=vpts;j++)
dilation[j] = (Vxyz[1][j] + Vxyz[2][j] + Vxyz[3][j]) / 3.0;
}
Modified: mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/Viscosity_structures.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -514,6 +514,7 @@
Vxyz[4][j] = 0.0;
Vxyz[5][j] = 0.0;
Vxyz[6][j] = 0.0;
+ dilation[j] = 0.0;
}
for(j=1; j<=ppts; j++) {
@@ -544,12 +545,10 @@
}
}
- if(E->control.inv_gruneisen > 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;
Modified: mc/3D/CitcomS/branches/compressible/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/lib/global_defs.h 2007-08-09 22:57:28 UTC (rev 7798)
@@ -157,6 +157,9 @@
struct EG {
higher_precision g[24][1]; };
+struct EC {
+ higher_precision c[24][1]; };
+
struct EK {
double k[24*24]; };
@@ -399,7 +402,6 @@
int restart;
int post_p;
int post_topo;
- int compressible;
char GEOMETRY[20]; /* one of ... */
int CART2D;
@@ -445,11 +447,19 @@
/* Rayleigh # */
float Atemp;
+
/* Dissipation # */
float disptn_number;
+
/* inverse of Gruneisen parameter */
float inv_gruneisen;
+ /**/
+ float relative_err_accuracy;
+
+ /**/
+ int compress_iter_maxstep;
+
float inputdiff;
float VBXtopval;
float VBXbotval;
@@ -515,7 +525,6 @@
double *conductivity;
double *gravity;
double *Tadi;
- double *dlnrhodr;
};
@@ -660,6 +669,7 @@
struct ID *ID[MAX_LEVELS][NCS];
struct SUBEL *EL[MAX_LEVELS][NCS];
struct EG *elt_del[MAX_LEVELS][NCS];
+ struct EC *elt_c[MAX_LEVELS][NCS];
struct EK *elt_k[MAX_LEVELS][NCS];
struct CC *cc[NCS];
struct CCX *ccx[NCS];
Modified: mc/3D/CitcomS/branches/compressible/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-08-09 22:06:35 UTC (rev 7797)
+++ mc/3D/CitcomS/branches/compressible/module/setProperties.c 2007-08-09 22:57:28 UTC (rev 7798)
@@ -451,7 +451,8 @@
getFloatProperty(properties, "rayleigh", E->control.Atemp, fp);
getFloatProperty(properties, "dissipation_number", E->control.disptn_number, fp);
getFloatProperty(properties, "gruneisen", tmp, fp);
- if(fabs(tmp) > 1e-6)
+ /* special case: if tmp==0, set gruneisen as inf */
+ if(tmp != 0)
E->control.inv_gruneisen = 1.0/tmp;
else
E->control.inv_gruneisen = 0;
@@ -772,6 +773,9 @@
getIntProperty(properties, "aug_lagr", E->control.augmented_Lagr, fp);
getDoubleProperty(properties, "aug_number", E->control.augmented, fp);
+ getIntProperty(properties, "compress_iter_maxstep", E->control.compress_iter_maxstep, fp);
+ getFloatProperty(properties, "relative_err_accuracy", E->control.relative_err_accuracy, fp);
+
PUTS(("\n"));
Py_INCREF(Py_None);
More information about the cig-commits
mailing list