[cig-commits] r7834 - in mc/3D/CitcomS/trunk:
CitcomS/Components/Stokes_solver lib module visual
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Aug 16 12:28:12 PDT 2007
Author: tan2
Date: 2007-08-16 12:28:10 -0700 (Thu, 16 Aug 2007)
New Revision: 7834
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/trunk/lib/Construct_arrays.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.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/setProperties.c
mc/3D/CitcomS/trunk/visual/pasteCitcomData.py
Log:
Merging compressible branch to trunk.
The TALA solver is working. Two new input parameters: dissipation_number and
gruneisen.
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2007-08-16 19:28:10 UTC (rev 7834)
@@ -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/trunk/lib/Construct_arrays.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Construct_arrays.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -683,6 +683,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/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -155,8 +155,8 @@
/*==============================================================
- Function to supply the element strain-displacement matrix Ba,
- which is used to compute element stiffness matrix
+ Function to supply the element strain-displacement matrix Ba at velocity
+ quadrature points, which is used to compute element stiffness matrix
============================================================== */
void get_ba(struct Shape_function *N, struct Shape_function_dx *GNx,
@@ -167,61 +167,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];
+ ra[k] = rtf[3][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)];
+ for(n=1;n<=dims;n++) {
+ 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];
- return;
+ 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;
}
+
/*==============================================================
Function to supply the element k matrix for a given element e.
============================================================== */
@@ -237,6 +238,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;
@@ -246,22 +251,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;
@@ -290,25 +279,26 @@
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(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];
+ 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] );
+ 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] );
- /**/
+
+ /**/
ad=dims*(a-1);
bd=dims*(b-1);
@@ -334,7 +324,7 @@
elt_k[qn+2*nn ] = bdbmu[1][3] ;
elt_k[qn+2*nn+1] = bdbmu[2][3] ;
elt_k[qn+2*nn+2] = bdbmu[3][3] ;
- /**/
+ /**/
} /* Sum over all the a,b's to obtain full elt_k matrix */
@@ -576,29 +566,33 @@
}
+/* =====================================================
+ 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];
+ }
}
return;
@@ -606,13 +600,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;
}
@@ -758,6 +756,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/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -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*);
@@ -127,6 +128,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);
@@ -384,9 +387,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(abs(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);
/* Q0_enriched gets read in Tracer_setup.c */
@@ -575,11 +584,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/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -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/trunk/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Pan_problem_misc_functions.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -143,6 +143,7 @@
temp = E->control.Atemp;
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.nno;i++) {
int nz = ((i-1) % E->lmesh.noz) + 1;
/* We don't need to substract adiabatic T profile from T here,
@@ -425,7 +426,7 @@
parallel_process_termination();
}
-/*
+/*
@@ -448,17 +449,17 @@
brange *= range; /* bottom */
trange *= range; /* top */
mrange *= range; /* middle */
-
+
nb = nz * bfrac;
nt = nz * tfrac;
- nm = nz-nb-nt;
+ nm = nz-nb-nt;
if((nm < 1)||(nt < 2)||(nb < 2))
myerror(E,"get_r_spacing_fine: refinement out of bounds");
-
+
drb = brange/(nb-1);
drt = trange/(nt-1);
drm = mrange / (nm + 1);
-
+
for(r=ri,k=1;k<=nb;k++,r+=drb){
rr[k] = r;
}
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -38,12 +38,15 @@
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,
+static float solve_Ahat_p_fhat_CG(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,
+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,24 +107,26 @@
{
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);
+ 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);
+ //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
*/
-static float solve_Ahat_p_fhat_BA(struct All_variables *E,
+static float solve_Ahat_p_fhat_CG(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 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_TALA(struct All_variables *E,
+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/trunk/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Topo_gravity.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -56,8 +56,6 @@
if (E->parallel.me_loc[3]==E->parallel.nprocz-1)
get_surf_stress(E,SXX,SYY,SZZ,SXY,SXZ,SZY);
-/* topo_scaling1=1.0/((E->data.density-E->data.density_above)*E->data.grav_acc); */
-/* topo_scaling2=1.0/((E->data.density_below-E->data.density)*E->data.grav_acc); */
topo_scaling1 = topo_scaling2 = 1.0;
@@ -276,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;
}
@@ -330,12 +328,6 @@
(E->exchange_node_f)(E,divv,lev);
(E->exchange_node_f)(E,vorv,lev);
- /* stress_scaling = 1.0e-6*E->data.ref_viscosity*E->data.therm_diff/ */
- /* (E->data.radius_km*E->data.radius_km); */
-
- /* velo_scaling = 100.*365.*24.*3600.*1.0e-3*E->data.therm_diff/E->data.radius_km; */
- /* cm/yr */
-
stress_scaling = velo_scaling = 1.0;
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -455,6 +447,7 @@
{
int m,k,ll,mm,node,i,j,p,noz,snode;
float *TT[NCS],radius,*geoid[2],dlayer,con1,con2,scaling2,scaling;
+ double buoy2rho;
void sphere_expansion();
void sum_across_depth_sph1();
@@ -481,13 +474,15 @@
/* loop over each layer */
for(k=1;k<E->lmesh.noz;k++) {
-
+ buoy2rho = scaling2 * E->refstate.rho[k]
+ * E->refstate.expansivity[k] / E->refstate.gravity[k];
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=E->lmesh.noy;i++)
for(j=1;j<=E->lmesh.nox;j++) {
node=k+(j-1)*E->lmesh.noz+(i-1)*E->lmesh.nox*E->lmesh.noz;
p = j + (i-1)*E->lmesh.nox;
- TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])*0.5*scaling2;
+ TT[m][p] = (E->buoyancy[m][node]+E->buoyancy[m][node+1])
+ * 0.5 * buoy2rho;
}
/* expand TT into spherical harmonics */
@@ -545,18 +540,14 @@
(E->data.radius_km*E->data.radius_km*1e6);
/* density contrast across surface */
- den_contrast1 = (E->data.density-E->data.density_above);
+ den_contrast1 = (E->data.density-E->data.density_above)*E->refstate.rho[E->lmesh.noz];
/* density contrast across CMB */
- den_contrast2 = (E->data.density_below-E->data.density);
+ den_contrast2 = (E->data.density_below-E->data.density)*E->refstate.rho[1];
/* scale for surface and CMB topo */
- topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc);
- topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc);
+ topo_scaling1 = stress_scaling / (den_contrast1*E->data.grav_acc*E->refstate.gravity[E->lmesh.noz]);
+ topo_scaling2 = stress_scaling / (den_contrast2*E->data.grav_acc*E->refstate.gravity[1]);
- /* scale for geoid */
- scaling = 1.0e3*4.0*M_PI*E->data.radius_km*
- E->data.grav_const/E->data.grav_acc;
-
if (E->parallel.me_loc[3] == E->parallel.nprocz-1) {
/* expand surface topography into sph. harm. */
sphere_expansion(E, E->slice.tpg, tpgt[0], tpgt[1]);
@@ -567,6 +558,10 @@
tpgt[j][i] *= topo_scaling1;
}
+ /* scale for geoid */
+ scaling = 1.0e3 * 4.0 * M_PI * E->data.radius_km * E->data.grav_const
+ / (E->data.grav_acc * E->refstate.gravity[E->lmesh.noz]);
+
/* compute geoid due to surface topo, skip degree-0 and 1 term */
for (j=0; j<2; j++)
for (ll=2; ll<=E->output.llmax; ll++) {
@@ -588,6 +583,11 @@
for (i=0; i<E->sphere.hindice; i++) {
tpgb[j][i] *= topo_scaling2;
}
+
+ /* scale for geoid */
+ scaling = 1.0e3 * 4.0 * M_PI * E->data.radius_km * E->data.grav_const
+ / (E->data.grav_acc * E->refstate.gravity[1]);
+
/* compute geoid due to bottom topo, skip degree-0 and 1 term */
for (j=0; j<2; j++)
for (ll=2; ll<=E->output.llmax; ll++) {
Modified: mc/3D/CitcomS/trunk/lib/Viscosity_structures.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/Viscosity_structures.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -748,6 +748,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++) {
@@ -778,12 +779,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/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2007-08-16 19:28:10 UTC (rev 7834)
@@ -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;
@@ -446,11 +448,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;
@@ -517,7 +527,6 @@
double *conductivity;
double *gravity;
double *Tadi;
- double *dlnrhodr;
};
@@ -674,6 +683,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/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2007-08-16 19:28:10 UTC (rev 7834)
@@ -447,8 +447,11 @@
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)
+ /* special case: if tmp==0, set gruneisen as inf */
+ if(tmp != 0)
E->control.inv_gruneisen = 1/tmp;
+ else
+ E->control.inv_gruneisen = 0;
getFloatProperty(properties, "Q0", E->control.Q0, fp);
@@ -753,6 +756,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);
Modified: mc/3D/CitcomS/trunk/visual/pasteCitcomData.py
===================================================================
--- mc/3D/CitcomS/trunk/visual/pasteCitcomData.py 2007-08-16 09:40:16 UTC (rev 7833)
+++ mc/3D/CitcomS/trunk/visual/pasteCitcomData.py 2007-08-16 19:28:10 UTC (rev 7834)
@@ -122,7 +122,7 @@
# how many header lines for each infix
headers = {'coord': 1,
'botm': 1,
- 'comp_nd': 1,
+ 'comp_nd': 2,
'pressure': 2,
'stress': 2,
'surf': 1,
More information about the cig-commits
mailing list