[cig-commits] r5880 - mc/3D/CitcomS/branches/compressible/lib
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Jan 23 17:10:13 PST 2007
Author: tan2
Date: 2007-01-23 17:10:13 -0800 (Tue, 23 Jan 2007)
New Revision: 5880
Modified:
mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c
Log:
Fixed a bug when computing pressure correction.
Using full spherical solver, Di=0, init T perturbation (l,m)=(3,2), and
default values for other parameters, the differences in CG and BiCGstab
solvers are 0.5% in velocity and stress, 0.2% in pressure.
Modified: mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-24 00:44:09 UTC (rev 5879)
+++ mc/3D/CitcomS/branches/compressible/lib/Pan_problem_misc_functions.c 2007-01-24 01:10:13 UTC (rev 5880)
@@ -138,7 +138,7 @@
for(i=1;i<=E->lmesh.nno;i++) {
int nz = ((i-1) % E->lmesh.noz) + 1;
buoy[m][i] = temp * E->rho_ref[nz] * E->thermexp_ref[nz]
- * (E->T[m][i] - E->T_ref[nz]);
+ * E->T[m][i];
}
phase_change_apply_410(E, buoy);
Modified: mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-01-24 00:44:09 UTC (rev 5879)
+++ mc/3D/CitcomS/branches/compressible/lib/Stokes_flow_Incomp.c 2007-01-24 01:10:13 UTC (rev 5880)
@@ -156,8 +156,7 @@
}
time0 = CPU_time0();
- valid = 1;
- r0dotz0 = 0;
+ count = 0;
/* calculate the initial velocity residual */
v_res = initial_vel_residual(E, V, P, F, imp);
@@ -167,8 +166,6 @@
assemble_div_u(E, V, r1, lev);
residual = incompressibility_residual(E, V, r1);
- 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,
@@ -183,6 +180,9 @@
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) ) {
@@ -325,11 +325,11 @@
int m, i, j, count, lev;
int valid;
- double alpha, delta, omega;
+ double alpha, beta, omega;
double r0dotrt, r1dotrt;
double residual, dpressure, dvelocity;
- double *r1[NCS], *r2[NCS], *z0[NCS], *p1[NCS], *p2[NCS];
+ 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];
@@ -345,10 +345,10 @@
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));
- z0[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));
@@ -358,7 +358,7 @@
}
time0 = CPU_time0();
- valid = 1;
+ count = 0;
/* calculate the initial velocity residual */
v_res = initial_vel_residual(E, V, P, F, imp);
@@ -375,9 +375,6 @@
rt[m][j] = r1[m][j];
- count = 0;
- r0dotrt = alpha = omega = 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,
@@ -392,15 +389,20 @@
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) ) {
+ ((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();
@@ -414,22 +416,22 @@
p2[m][j] = r1[m][j];
else {
/* p2 = r1 + <r1,rt>/<r0,rt> * alpha/omega * (p1 - omega*v0) */
- delta = (r1dotrt / r0dotrt) * (alpha / omega);
+ 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] + delta
+ p2[m][j] = r1[m][j] + beta
* (p1[m][j] - omega * v0[m][j]);
}
- /* preconditioner BPI ~= inv(K), z0 = BPI*p2 */
+ /* preconditioner BPI ~= inv(K), pt = BPI*p2 */
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
- z0[m][j] = E->BPI[lev][m][j] * p2[m][j];
+ pt[m][j] = E->BPI[lev][m][j] * p2[m][j];
- /* solve K*u0 = grad(p2) for u1 */
- assemble_grad_p(E, p2, F, lev);
+ /* 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);
@@ -484,30 +486,37 @@
r2[m][j] = s0[m][j] - omega * t0[m][j];
- /* P = P + alpha * z0 + omega * st */
+ /* P = P + alpha * pt + omega * st */
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
- P[m][j] += alpha * z0[m][j] + omega * st[m][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++)
- V[m][j] -= alpha * u0[m][j] + omega * E->u1[m][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 = alpha * sqrt( (global_pdot(E, z0, z0, lev) +
- global_pdot(E, st, st, lev))
- / (1.0e-32 + global_pdot(E, P, P, lev)));
- dvelocity = alpha * sqrt( (global_vdot(E, u0, u0, lev) +
- global_vdot(E, E->u1, E->u1, lev))
- / (1.0e-32 + E->monitor.vdotv));
+ 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) {
@@ -542,7 +551,7 @@
for(m=1; m<=E->sphere.caps_per_proc; m++) {
free((void *) r1[m]);
free((void *) r2[m]);
- free((void *) z0[m]);
+ free((void *) pt[m]);
free((void *) p1[m]);
free((void *) p2[m]);
free((void *) rt[m]);
More information about the cig-commits
mailing list