[cig-commits] r13196 - in mc/3D/CitcomS/trunk: CitcomS/Components/Stokes_solver examples/Cookbook6 examples/Cookbook7 examples/Cookbook8 examples/Cookbook9 examples/Full examples/Regional lib module module/Exchanger tests
tan2 at geodynamics.org
tan2 at geodynamics.org
Wed Oct 29 16:17:11 PDT 2008
Author: tan2
Date: 2008-10-29 16:17:11 -0700 (Wed, 29 Oct 2008)
New Revision: 13196
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg
mc/3D/CitcomS/trunk/examples/Cookbook7/cookbook7.cfg
mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
mc/3D/CitcomS/trunk/examples/Cookbook9/cookbook9.cfg
mc/3D/CitcomS/trunk/examples/Full/input.sample
mc/3D/CitcomS/trunk/examples/Regional/input.sample
mc/3D/CitcomS/trunk/lib/Checkpoints.c
mc/3D/CitcomS/trunk/lib/Element_calculations.c
mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
mc/3D/CitcomS/trunk/lib/Global_operations.c
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output_h5.c
mc/3D/CitcomS/trunk/lib/Output_vtk.c
mc/3D/CitcomS/trunk/lib/Size_does_matter.c
mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/Exchanger/AreaWeightedNormal.cc
mc/3D/CitcomS/trunk/module/Exchanger/misc.cc
mc/3D/CitcomS/trunk/module/setProperties.c
mc/3D/CitcomS/trunk/tests/coupled.cfg
mc/3D/CitcomS/trunk/tests/multicoupled.cfg
mc/3D/CitcomS/trunk/tests/test2.sh
mc/3D/CitcomS/trunk/tests/test5.sh
Log:
Fixed two bugs in vtk velocity output
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Stokes_solver/Incompressible.py 2008-10-29 23:17:11 UTC (rev 13196)
@@ -89,8 +89,7 @@
node_assemble = prop.bool("node_assemble", default=True)
precond = prop.bool("precond", default=True)
- accuracy = prop.float("accuracy", default=1.0e-6)
- tole_compressibility = prop.float("tole_compressibility", default=1.0e-7)
+ accuracy = prop.float("accuracy", default=1.0e-4)
mg_cycle = prop.int("mg_cycle", default=1)
down_heavy = prop.int("down_heavy", default=3)
up_heavy = prop.int("up_heavy", default=3)
@@ -105,9 +104,11 @@
uzawa = prop.str("uzawa", default="cg",
validator=prop.choice(["cg", "bicg"]))
compress_iter_maxstep = prop.int("compress_iter_maxstep", default=100)
- relative_err_accuracy = prop.float("relative_err_accuracy", default=0.001)
remove_rigid_rotation = prop.bool("remove_rigid_rotation", default=True)
+ # Not used. Retained here for backward compatibility.
+ tole_compressibility = prop.float("tole_compressibility", default=1.0e-7)
+ relative_err_accuracy = prop.float("relative_err_accuracy", default=0.001)
# version
__id__ = "$Id$"
Modified: mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Cookbook6/cookbook6.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -30,8 +30,7 @@
[CitcomS.solver.vsolver]
precond = on
-accuracy = 1e-10
-tole_compressibility = 1e-06
+accuracy = 1e-6
vlowstep = 100000
piterations = 100000
Modified: mc/3D/CitcomS/trunk/examples/Cookbook7/cookbook7.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook7/cookbook7.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Cookbook7/cookbook7.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -47,9 +47,8 @@
[CitcomS.solver.vsolver]
Solver = cgrad
-accuracy = 1e-06
+accuracy = 1e-04
vlowstep = 1000
-tole_compressibility = 1e-07
piterations = 1000
Modified: mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Cookbook8/cookbook8.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -91,9 +91,7 @@
piterations = 375
accuracy = 0.001
-tole_compressibility = 1e-08
compress_iter_maxstep = 100
-relative_err_accuracy = 0.001
remove_rigid_rotation = on
Modified: mc/3D/CitcomS/trunk/examples/Cookbook9/cookbook9.cfg
===================================================================
--- mc/3D/CitcomS/trunk/examples/Cookbook9/cookbook9.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Cookbook9/cookbook9.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -80,10 +80,6 @@
adv_sub_iterations = 2
-[CitcomS.csolver.vsolver]
-tole_compressibility = 8e-6
-
-
[CitcomS.csolver.visc]
VISC_UPDATE = on
num_mat = 4
@@ -125,10 +121,6 @@
solution_cycles_init = 0
-[CitcomS.esolver.vsolver]
-tole_compressibility = 6e-6
-
-
[CitcomS.esolver.tsolver]
ADV = on
filter_temp = off
Modified: mc/3D/CitcomS/trunk/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Full/input.sample 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Full/input.sample 2008-10-29 23:17:11 UTC (rev 13196)
@@ -211,7 +211,6 @@
piterations=375
accuracy=1.0e-6
-tole_compressibility=1.0e-7
ADV=on
fixed_timestep=0.0
Modified: mc/3D/CitcomS/trunk/examples/Regional/input.sample
===================================================================
--- mc/3D/CitcomS/trunk/examples/Regional/input.sample 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/examples/Regional/input.sample 2008-10-29 23:17:11 UTC (rev 13196)
@@ -214,7 +214,6 @@
piterations=375
accuracy=1.0e-6
-tole_compressibility=1.0e-7
ADV=on
fixed_timestep=0.0
Modified: mc/3D/CitcomS/trunk/lib/Checkpoints.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Checkpoints.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Checkpoints.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -443,13 +443,14 @@
static void momentum_checkpoint(struct All_variables *E, FILE *fp)
{
- int m, i;
- int lev = E->mesh.levmax;
+ int m;
+ float junk[2];
+ junk[0] = junk[1] = 0;
write_sentinel(fp);
- fwrite(&(E->monitor.vdotv), sizeof(float), 1, fp);
- fwrite(&(E->monitor.incompressibility), sizeof(float), 1, fp);
+ /* for backward compatibility */
+ fwrite(junk, sizeof(float), 2, fp);
/* the 0-th element of P/NP/EVI/VI is not init'd
* and won't be used when read it. */
@@ -469,16 +470,18 @@
{
void v_from_vector();
void p_to_nodes();
+ double global_v_norm2(), global_p_norm2();
- int m, i;
+ int m;
int lev = E->mesh.levmax;
+ float junk[2];
read_sentinel(fp, E->parallel.me);
- if(fread(&(E->monitor.vdotv), sizeof(float), 1, fp)!=1)
+ /* for backward compatibility */
+ if(fread(junk, sizeof(float), 2, fp)!=2)
myerror("read_momentum_checkpoint: error at vdotv",E);
- if(fread(&(E->monitor.incompressibility), sizeof(float), 1, fp)!=1)
- myerror("read_momentum_checkpoint: error at incomp",E);
+
for(m=1; m<=E->sphere.caps_per_proc; m++) {
/* Pressure at equation points */
if(fread(E->P[m], sizeof(double), E->lmesh.npno+1, fp) != E->lmesh.npno+1)
@@ -488,6 +491,9 @@
myerror("read_momentum_checkpoint: error at U",E);
}
+ E->monitor.vdotv = global_v_norm2(E, E->U);
+ E->monitor.pdotp = global_p_norm2(E, E->P);
+
/* update velocity array */
v_from_vector(E);
Modified: mc/3D/CitcomS/trunk/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Element_calculations.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Element_calculations.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -84,6 +84,7 @@
void get_elt_f();
void get_elt_tr();
void strip_bcs_from_residual();
+ double global_vdot();
const int neq=E->lmesh.neq;
const int nel=E->lmesh.nel;
@@ -115,6 +116,17 @@
(E->solver.exchange_id_d)(E, E->F, lev);
strip_bcs_from_residual(E,E->F,lev);
+
+ /* compute the norm of E->F */
+ E->monitor.fdotf = sqrt(global_vdot(E, E->F, E->F, lev));
+
+ if(E->parallel.me==0) {
+ fprintf(stderr, "Momentum equation force %.9e\n",
+ E->monitor.fdotf);
+ fprintf(E->fp, "Momentum equation force %.9e\n",
+ E->monitor.fdotf);
+ }
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/General_matrix_functions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/General_matrix_functions.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -65,7 +65,7 @@
void report();
int count,counts,cycles,convergent,valid;
- int i, neq, gneq,m;
+ int i, neq, m;
char message[200];
@@ -73,7 +73,6 @@
double residual,prior_residual,r0;
double *D1[NCS], *r[NCS], *Au[NCS];
- gneq = E->mesh.NEQ[high_lev];
neq = E->lmesh.NEQ[high_lev];
for (m=1;m<=E->sphere.caps_per_proc;m++)
@@ -81,7 +80,7 @@
d0[m][i] = 0.0;
}
- r0=residual=sqrt(global_vdot(E,F,F,high_lev)/gneq);
+ r0=residual=sqrt(global_vdot(E,F,F,high_lev));
prior_residual=2*residual;
count = 0;
@@ -140,7 +139,7 @@
count++;
-
+ E->monitor.momentum_residual = residual;
E->control.total_iteration_cycles += count;
E->control.total_v_solver_calls += 1;
@@ -282,7 +281,7 @@
d1[m][j]+=vel[levmax][m][j];
}
- residual = sqrt(global_vdot(E,F,F,hl)/E->mesh.NEQ[hl]);
+ residual = sqrt(global_vdot(E,F,F,hl));
for(i=E->mesh.levmin;i<=E->mesh.levmax;i++)
for(m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -352,7 +351,7 @@
d0[m][i] = 0.0;
}
- residual = sqrt(global_vdot(E,r1,r1,level)/E->mesh.NEQ[level]);
+ residual = sqrt(global_vdot(E,r1,r1,level));
assert(residual != 0.0 /* initial residual for CG = 0.0 */);
count = 0;
@@ -394,7 +393,7 @@
r2[m][i] = r1[m][i] - alpha * Ap[m][i];
}
- residual = sqrt(global_vdot(E,r2,r2,level)/E->mesh.NEQ[level]);
+ residual = sqrt(global_vdot(E,r2,r2,level));
for(m=1;m<=E->sphere.caps_per_proc;m++) {
shuffle[m] = r0[m]; r0[m] = r1[m]; r1[m] = r2[m]; r2[m] = shuffle[m];
Modified: mc/3D/CitcomS/trunk/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Global_operations.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Global_operations.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -583,9 +583,78 @@
MPI_Allreduce(&temp, &prod,1,MPI_DOUBLE,MPI_SUM,E->parallel.world);
return (prod);
- }
+}
+/* return ||V||^2 */
+double global_v_norm2(struct All_variables *E, double **V)
+{
+ int i, m, d;
+ int eqn1, eqn2, eqn3;
+ double prod, temp;
+
+ temp = 0.0;
+ prod = 0.0;
+ for (m=1; m<=E->sphere.caps_per_proc; m++)
+ for (i=1; i<=E->lmesh.nno; i++) {
+ eqn1 = E->id[m][i].doff[1];
+ eqn2 = E->id[m][i].doff[2];
+ eqn3 = E->id[m][i].doff[3];
+ /* L2 norm */
+ temp += (V[m][eqn1] * V[m][eqn1] +
+ V[m][eqn2] * V[m][eqn2] +
+ V[m][eqn3] * V[m][eqn3]) * E->NMass[m][i];
+ }
+
+ MPI_Allreduce(&temp, &prod, 1, MPI_DOUBLE, MPI_SUM, E->parallel.world);
+
+ return (prod/E->mesh.volume);
+}
+
+
+/* return ||P||^2 */
+double global_p_norm2(struct All_variables *E, double **P)
+{
+ int i, m;
+ double prod, temp;
+
+ temp = 0.0;
+ prod = 0.0;
+ for (m=1; m<=E->sphere.caps_per_proc; m++)
+ for (i=1; i<=E->lmesh.npno; i++) {
+ /* L2 norm */
+ temp += P[m][i] * P[m][i] * E->eco[m][i].area;
+ }
+
+ MPI_Allreduce(&temp, &prod, 1, MPI_DOUBLE, MPI_SUM, E->parallel.world);
+
+ return (prod/E->mesh.volume);
+}
+
+
+/* return ||A||^2, where A_i is \int{div(u) d\Omega_i} */
+double global_div_norm2(struct All_variables *E, double **A)
+{
+ int i, m;
+ double prod, temp;
+
+ temp = 0.0;
+ prod = 0.0;
+ for (m=1; m<=E->sphere.caps_per_proc; m++)
+ for (i=1; i<=E->lmesh.npno; i++) {
+ /* L2 norm of div(u) */
+ temp += A[m][i] * A[m][i] / E->eco[m][i].area;
+
+ /* L1 norm */
+ /*temp += fabs(A[m][i]);*/
+ }
+
+ MPI_Allreduce(&temp, &prod, 1, MPI_DOUBLE, MPI_SUM, E->parallel.world);
+
+ return (prod/E->mesh.volume);
+}
+
+
double global_tdot_d(E,A,B,lev)
struct All_variables *E;
double **A,**B;
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -605,8 +605,6 @@
input_boolean("remove_rigid_rotation",&(E->control.remove_rigid_rotation),"on",m);
- input_float("tole_compressibility",&(E->control.tole_comp),"0.0",m);
-
input_boolean("self_gravitation",&(E->control.self_gravitation),"off",m);
input_boolean("use_cbf_topo",&(E->control.use_cbf_topo),"off",m); /* make default on later XXX TWB */
@@ -649,7 +647,6 @@
if(strcmp(E->control.uzawa, "cg") == 0) {
/* more convergence parameters for "cg" */
input_int("compress_iter_maxstep",&(E->control.compress_iter_maxstep),"100",m);
- input_float("relative_err_accuracy",&(E->control.relative_err_accuracy),"0.001",m);
}
else if(strcmp(E->control.uzawa, "bicg") == 0) {
}
@@ -809,6 +806,9 @@
/* lump mass matrix for the energy eqn */
E->TMass[j] = (double *) malloc((nno+1)*sizeof(double));
+ /* nodal mass */
+ E->NMass[j] = (double *) malloc((nno+1)*sizeof(double));
+
nxyz = max(nox*noz,nox*noy);
nxyz = 2*max(nxyz,noz*noy);
@@ -856,7 +856,7 @@
E->GNX[i][j] = (struct Shape_function_dx *)malloc((nel+1)*sizeof(struct Shape_function_dx));
E->GDA[i][j] = (struct Shape_function_dA *)malloc((nel+1)*sizeof(struct Shape_function_dA));
- E->MASS[i][j] = (float *) malloc((nno+1)*sizeof(float));
+ E->MASS[i][j] = (double *) malloc((nno+1)*sizeof(double));
E->ECO[i][j] = (struct COORD *) malloc((nno+2)*sizeof(struct COORD));
E->TWW[i][j] = (struct FNODE *) malloc((nel+2)*sizeof(struct FNODE));
@@ -964,6 +964,11 @@
{
int m,n,i,j,k,l;
+ E->monitor.incompressibility = 0;
+ E->monitor.fdotf = 0;
+ E->monitor.vdotv = 0;
+ E->monitor.pdotp = 0;
+
m=0;
n=1;
for (j=1;j<=E->sphere.caps_per_proc;j++) {
@@ -1023,8 +1028,7 @@
E->control.v_steps_low = 10;
E->control.v_steps_upper = 1;
- E->control.accuracy = 1.0e-6;
- E->control.vaccuracy = 1.0e-8;
+ E->control.accuracy = 1.0e-4;
E->control.verbose=0; /* debugging/profiles */
/* SECOND: values for which an obvious default setting is useful */
Modified: mc/3D/CitcomS/trunk/lib/Output_h5.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_h5.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Output_h5.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -1610,7 +1610,6 @@
status = set_attribute_int(input, "precond", E->control.precondition);
status = set_attribute_double(input, "accuracy", E->control.accuracy);
- status = set_attribute_float(input, "tole_compressibility", E->control.tole_comp);
status = set_attribute_int(input, "mg_cycle", E->control.mg_cycle);
status = set_attribute_int(input, "down_heavy", E->control.down_heavy);
Modified: mc/3D/CitcomS/trunk/lib/Output_vtk.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_vtk.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Output_vtk.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -123,7 +123,7 @@
{
int i, j;
double sint, sinf, cost, cosf;
- float *V[3];
+ float *V[4];
const int lev = E->mesh.levmax;
fputs(" <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n", fp);
@@ -139,7 +139,7 @@
cost = E->SinCos[lev][j][2][i];
cosf = E->SinCos[lev][j][3][i];
- fprintf(fp, "%.6e %.6e %.6e %.6e\n",
+ fprintf(fp, "%.6e %.6e %.6e\n",
V[1][i]*cost*cosf - V[2][i]*sinf + V[3][i]*sint*cosf,
V[1][i]*cost*sinf + V[2][i]*cosf + V[3][i]*sint*sinf,
-V[1][i]*sint + V[3][i]*cost);
Modified: mc/3D/CitcomS/trunk/lib/Size_does_matter.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Size_does_matter.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -1102,8 +1102,13 @@
} /* end of for m */
+ if(lev == E->mesh.levmax)
+ for (m=1;m<=E->sphere.caps_per_proc;m++)
+ for(node=1;node<=E->lmesh.NNO[lev];node++)
+ E->NMass[m][node] = E->MASS[lev][m][node];
+
if (E->control.NMULTIGRID||E->control.EMULTIGRID||E->mesh.levmax==lev)
- (E->exchange_node_f)(E,E->MASS[lev],lev);
+ (E->exchange_node_d)(E,E->MASS[lev],lev);
for (m=1;m<=E->sphere.caps_per_proc;m++)
for(node=1;node<=E->lmesh.NNO[lev];node++)
Modified: mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/Stokes_flow_Incomp.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -30,6 +30,7 @@
/* Functions which solve for the velocity and pressure fields using Uzawa-type iteration loop. */
#include <math.h>
+#include <string.h>
#include <sys/types.h>
#include "element_definitions.h"
#include "global_defs.h"
@@ -37,23 +38,21 @@
void myerror(struct All_variables *,char *);
-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_CG(struct All_variables *E,
- double **V, double **P, double **F,
- double imp, int *steps_max);
-static float solve_Ahat_p_fhat_BiCG(struct All_variables *E,
+static void solve_Ahat_p_fhat(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max);
+static void solve_Ahat_p_fhat_CG(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max);
+static void 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,
+static void 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);
-static double incompressibility_residual(struct All_variables *E,
- double **V, double **r);
+static void initial_vel_residual(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp);
/* Master loop for pressure and (hence) velocity field */
@@ -103,45 +102,85 @@
/* ========================================================================= */
+static double momentum_eqn_residual(struct All_variables *E,
+ double **V, double **P, double **F)
+{
+ /* Compute the norm of (F - grad(P) - K*V)
+ * This norm is ~= E->monitor.momentum_residual */
+ void assemble_del2_u();
+ void assemble_grad_p();
+ void strip_bcs_from_residual();
+ double global_v_norm2();
+
+ int i, m;
+ double *r1[NCS], *r2[NCS];
+ double res;
+ const int lev = E->mesh.levmax;
+ const int neq = E->lmesh.neq;
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ r1[m] = malloc((neq+1)*sizeof(double));
+ r2[m] = malloc((neq+1)*sizeof(double));
+ }
+
+ /* r2 = F - grad(P) - K*V */
+ assemble_grad_p(E, P, E->u1, lev);
+ assemble_del2_u(E, V, r1, lev, 1);
+ for(m=1; m<=E->sphere.caps_per_proc; m++)
+ for(i=0; i<neq; i++)
+ r2[m][i] = F[m][i] - E->u1[m][i] - r1[m][i];
+
+ strip_bcs_from_residual(E, r2, lev);
+
+ res = sqrt(global_v_norm2(E, r2));
+
+ for(m=1; m<=E->sphere.caps_per_proc; m++) {
+ free(r1[m]);
+ free(r2[m]);
+ }
+ return(res);
+}
+
+
static void print_convergence_progress(struct All_variables *E,
int count, double time0,
- double sq_vdotv,
- double dvelocity, double dpressure)
+ double v_norm, double p_norm,
+ double dv, double dp,
+ double div)
{
- double CPU_time0();
+ double CPU_time0(), t;
+ t = CPU_time0() - time0;
- fprintf(E->fp, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
- "dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv,E->monitor.incompressibility,
- dvelocity, dpressure, E->monitor.solution_cycles);
- fprintf(stderr, "AhatP (%03d) after %6.2f s v=%.3e div/v=%.3e "
- "dv/v=%.3e and dp/p=%.3e for step %d\n",
- count, CPU_time0()-time0, sq_vdotv, E->monitor.incompressibility,
- dvelocity, dpressure, E->monitor.solution_cycles);
+ fprintf(E->fp, "(%03d) %5.1f s v=%e p=%e "
+ "div/v=%.2e dv/v=%.2e dp/p=%.2e step %d\n",
+ count, t, v_norm, p_norm, div, dv, dp,
+ E->monitor.solution_cycles);
+ fprintf(stderr, "(%03d) %5.1f s v=%e p=%e "
+ "div/v=%.2e dv/v=%.2e dp/p=%.2e step %d\n",
+ count, t, v_norm, p_norm, div, dv, dp,
+ E->monitor.solution_cycles);
return;
}
-static float solve_Ahat_p_fhat(struct All_variables *E,
+static void 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 == 0)
- residual = solve_Ahat_p_fhat_CG(E, V, P, F, imp, steps_max);
+ solve_Ahat_p_fhat_CG(E, V, P, F, imp, steps_max);
else {
if(strcmp(E->control.uzawa, "cg") == 0)
- residual = solve_Ahat_p_fhat_iterCG(E, V, P, F, imp, steps_max);
+ solve_Ahat_p_fhat_iterCG(E, V, P, F, imp, steps_max);
else if(strcmp(E->control.uzawa, "bicg") == 0)
- residual = solve_Ahat_p_fhat_BiCG(E, V, P, F, imp, steps_max);
+ solve_Ahat_p_fhat_BiCG(E, V, P, F, imp, steps_max);
else
myerror(E, "Error: unknown Uzawa iteration\n");
}
- return(residual);
+ return;
}
@@ -149,22 +188,25 @@
* conjugate gradient (CG) iterations
*/
-static float solve_Ahat_p_fhat_CG(struct All_variables *E,
- double **V, double **P, double **FF,
- double imp, int *steps_max)
+static void solve_Ahat_p_fhat_CG(struct All_variables *E,
+ double **V, double **P, double **FF,
+ double imp, int *steps_max)
{
int m, j, count, valid, lev, npno, neq;
- int gnpno;
- double *r1[NCS], *r2[NCS], *z1[NCS], *s1[NCS], *s2[NCS], *F[NCS];
+ double *r1[NCS], *r2[NCS], *z1[NCS], *s1[NCS], *s2[NCS], *cu[NCS];
+ double *F[NCS];
double *shuffle[NCS];
- double alpha, delta, r0dotz0, r1dotz1,sq_vdotv;
- double residual, v_res;
+ double alpha, delta, r0dotz0, r1dotz1;
+ double v_res;
- double global_vdot(), global_pdot();
+ double global_pdot();
+ double global_v_norm2(), global_p_norm2(), global_div_norm2();
double time0, CPU_time0();
- double dpressure, dvelocity;
+ double v_norm, p_norm;
+ double dvelocity, dpressure;
+ int converging;
void assemble_c_u();
void assemble_div_u();
@@ -174,7 +216,6 @@
int solve_del2_u();
void parallel_process_termination();
- gnpno = E->mesh.npno;
npno = E->lmesh.npno;
neq = E->lmesh.neq;
lev = E->mesh.levmax;
@@ -186,12 +227,13 @@
z1[m] = (double *)malloc((npno+1)*sizeof(double));
s1[m] = (double *)malloc((npno+1)*sizeof(double));
s2[m] = (double *)malloc((npno+1)*sizeof(double));
+ cu[m] = (double *)malloc((npno+1)*sizeof(double));
}
time0 = CPU_time0();
count = 0;
+ v_res = E->monitor.fdotf;
-
/* copy the original force vector since we need to keep it intact
between iterations */
for(m=1;m<=E->sphere.caps_per_proc;m++)
@@ -203,14 +245,17 @@
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;
+ cu[m][j] = 0.0;
- assemble_c_u(E, V, r2, lev);
+ assemble_c_u(E, V, cu, lev);
}
/* calculate the initial velocity residual */
- v_res = initial_vel_residual(E, V, P, F, imp);
+ /* In the compressible case, the initial guess of P might be bad.
+ * Do not correct V with it. */
+ if(E->control.inv_gruneisen == 0)
+ initial_vel_residual(E, V, P, F, imp*v_res);
/* initial residual r1 = div(V) */
@@ -221,32 +266,34 @@
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];
+ r1[m][j] += cu[m][j];
}
- residual = incompressibility_residual(E, V, r1);
+ E->monitor.vdotv = global_v_norm2(E, V);
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, r1)
+ / (1e-32 + E->monitor.vdotv));
-
- sq_vdotv = sqrt(E->monitor.vdotv);
-
- /* pressure and velocity corrections */
+ v_norm = sqrt(E->monitor.vdotv);
+ p_norm = sqrt(E->monitor.pdotp);
+ dvelocity = 1.0;
dpressure = 1.0;
- dvelocity = 1.0;
+ converging = 0;
-
if (E->control.print_convergence && E->parallel.me==0) {
- print_convergence_progress(E, count, time0, sq_vdotv,
- dvelocity, dpressure);
+ print_convergence_progress(E, count, time0,
+ v_norm, p_norm,
+ dvelocity, dpressure,
+ E->monitor.incompressibility);
}
r0dotz0 = 0;
while( (count < *steps_max) &&
- (E->monitor.incompressibility >= E->control.tole_comp) &&
- (dpressure >= imp) && (dvelocity >= imp) ) {
+ (E->monitor.incompressibility > imp) &&
+ (converging < 2) ) {
+ /* require two consecutive converging iterations to quit the while-loop */
-
/* preconditioner BPI ~= inv(K), z1 = BPI*r1 */
for(m=1; m<=E->sphere.caps_per_proc; m++)
for(j=1; j<=npno; j++)
@@ -287,11 +334,7 @@
/* alpha = <r1, z1> / <s2, F> */
- if(valid)
- /* alpha defined this way is the same as R&W */
- alpha = r1dotz1 / global_pdot(E, s2, F, lev);
- else
- alpha = 0.0;
+ alpha = r1dotz1 / global_pdot(E, s2, F, lev);
/* r2 = r1 - alpha * div(u1) */
@@ -313,22 +356,36 @@
/* compute velocity and incompressibility residual */
- assemble_div_u(E, V, F, lev);
- incompressibility_residual(E, V, F);
+ E->monitor.vdotv = global_v_norm2(E, V);
+ E->monitor.pdotp = global_p_norm2(E, P);
+ v_norm = sqrt(E->monitor.vdotv);
+ p_norm = sqrt(E->monitor.pdotp);
+ dvelocity = alpha * sqrt(global_v_norm2(E, E->u1) / (1e-32 + E->monitor.vdotv));
+ dpressure = alpha * sqrt(global_v_norm2(E, s2) / (1e-32 + E->monitor.pdotp));
- /* compute velocity and pressure corrections */
- dpressure = alpha * sqrt(global_pdot(E, s2, s2, lev)
- / (1.0e-32 + global_pdot(E, P, P, lev)));
- dvelocity = alpha * sqrt(global_vdot(E, E->u1, E->u1, lev)
- / (1.0e-32 + E->monitor.vdotv));
+ /* how many consecutive converging iterations? */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+ assemble_div_u(E, V, z1, lev);
+ if(E->control.inv_gruneisen != 0)
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(j=1;j<=npno;j++) {
+ z1[m][j] += cu[m][j];
+ }
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, z1)
+ / (1e-32 + E->monitor.vdotv));
+
count++;
- sq_vdotv = sqrt(E->monitor.vdotv);
if (E->control.print_convergence && E->parallel.me==0) {
- print_convergence_progress(E, count, time0, sq_vdotv,
- dvelocity, dpressure);
+ print_convergence_progress(E, count, time0,
+ v_norm, p_norm,
+ dvelocity, dpressure,
+ E->monitor.incompressibility);
}
@@ -348,6 +405,14 @@
} /* end loop for conjugate gradient */
+ assemble_div_u(E, V, z1, lev);
+ if(E->control.inv_gruneisen != 0)
+ for(m=1;m<=E->sphere.caps_per_proc;m++)
+ for(j=1;j<=npno;j++) {
+ z1[m][j] += cu[m][j];
+ }
+
+
for(m=1; m<=E->sphere.caps_per_proc; m++) {
free((void *) F[m]);
free((void *) r1[m]);
@@ -355,20 +420,21 @@
free((void *) z1[m]);
free((void *) s1[m]);
free((void *) s2[m]);
+ free((void *) cu[m]);
}
*steps_max=count;
- return(residual);
+ return;
}
/* 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 **FF,
- double imp, int *steps_max)
+static void solve_Ahat_p_fhat_BiCG(struct All_variables *E,
+ double **V, double **P, double **FF,
+ double imp, int *steps_max)
{
void assemble_div_rho_u();
void assemble_del2_u();
@@ -377,17 +443,19 @@
int solve_del2_u();
void parallel_process_termination();
- double global_vdot(), global_pdot();
+ double global_pdot();
+ double global_v_norm2(), global_p_norm2(), global_div_norm2();
double CPU_time0();
- int gnpno;
int npno, neq;
int m, j, count, lev;
int valid;
- double alpha, beta, omega,sq_vdotv;
+ double alpha, beta, omega;
double r0dotrt, r1dotrt;
- double residual, dpressure, dvelocity;
+ double v_norm, p_norm;
+ double dvelocity, dpressure;
+ int converging;
double *F[NCS];
double *r1[NCS], *r2[NCS], *pt[NCS], *p1[NCS], *p2[NCS];
@@ -397,7 +465,6 @@
double time0, v_res;
- gnpno = E->mesh.npno;
npno = E->lmesh.npno;
neq = E->lmesh.neq;
lev = E->mesh.levmax;
@@ -420,6 +487,7 @@
time0 = CPU_time0();
count = 0;
+ v_res = E->monitor.fdotf;
/* copy the original force vector since we need to keep it intact
between iterations */
@@ -429,42 +497,45 @@
/* calculate the initial velocity residual */
- v_res = initial_vel_residual(E, V, P, F, imp);
+ initial_vel_residual(E, V, P, F, imp*v_res);
/* initial residual r1 = div(rho_ref*V) */
assemble_div_rho_u(E, V, r1, lev);
- residual = incompressibility_residual(E, V, r1);
+ E->monitor.vdotv = global_v_norm2(E, V);
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, r1)
+ / (1e-32 + E->monitor.vdotv));
- /* 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];
-
-
- sq_vdotv = sqrt(E->monitor.vdotv);
-
- /* pressure and velocity corrections */
+ v_norm = sqrt(E->monitor.vdotv);
+ p_norm = sqrt(E->monitor.pdotp);
+ dvelocity = 1.0;
dpressure = 1.0;
- dvelocity = 1.0;
+ converging = 0;
if (E->control.print_convergence && E->parallel.me==0) {
- print_convergence_progress(E, count, time0, sq_vdotv,
- dvelocity, dpressure);
+ print_convergence_progress(E, count, time0,
+ v_norm, p_norm,
+ dvelocity, dpressure,
+ E->monitor.incompressibility);
}
+ /* 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];
+
+
valid = 1;
r0dotrt = alpha = omega = 0;
while( (count < *steps_max) &&
- ((E->monitor.incompressibility >= E->control.tole_comp) &&
- (dpressure >= imp) && (dvelocity >= imp)) ) {
+ (E->monitor.incompressibility > imp) &&
+ (converging < 2) ) {
+ /* require two consecutive converging iterations to quit the while-loop */
-
-
/* r1dotrt = <r1, rt> */
r1dotrt = global_pdot(E, r1, rt, lev);
if(r1dotrt == 0.0) {
@@ -571,23 +642,31 @@
/* compute velocity and incompressibility residual */
+ E->monitor.vdotv = global_v_norm2(E, V);
+ E->monitor.pdotp = global_p_norm2(E, P);
+ v_norm = sqrt(E->monitor.vdotv);
+ p_norm = sqrt(E->monitor.pdotp);
+ dvelocity = sqrt(global_v_norm2(E, F) / (1e-32 + E->monitor.vdotv));
+ dpressure = sqrt(global_v_norm2(E, s0) / (1e-32 + E->monitor.pdotp));
+
+ /* how many consecutive converging iterations? */
+ if(dvelocity < imp && dpressure < imp)
+ converging++;
+ else
+ converging = 0;
+
assemble_div_rho_u(E, V, t0, lev);
- incompressibility_residual(E, V, t0);
+ E->monitor.incompressibility = sqrt(global_div_norm2(E, t0)
+ / (1e-32 + E->monitor.vdotv));
- /* 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++;
- sq_vdotv = sqrt(E->monitor.vdotv);
-
if(E->control.print_convergence && E->parallel.me==0) {
- print_convergence_progress(E, count, time0, sq_vdotv,
- dvelocity, dpressure);
+ print_convergence_progress(E, count, time0,
+ v_norm, p_norm,
+ dvelocity, dpressure,
+ E->monitor.incompressibility);
}
@@ -626,7 +705,7 @@
*steps_max=count;
- return(residual);
+ return;
}
@@ -635,21 +714,23 @@
* conjugate gradient (CG) iterations with an outer iteration
*/
-static float solve_Ahat_p_fhat_iterCG(struct All_variables *E,
- double **V, double **P, double **F,
- double imp, int *steps_max)
+static void solve_Ahat_p_fhat_iterCG(struct All_variables *E,
+ double **V, double **P, double **F,
+ double imp, int *steps_max)
{
int m, i;
int cycles, num_of_loop;
- double residual;
double relative_err_v, relative_err_p;
double *old_v[NCS], *old_p[NCS],*diff_v[NCS],*diff_p[NCS];
+ double div_res;
const int npno = E->lmesh.npno;
const int neq = E->lmesh.neq;
const int lev = E->mesh.levmax;
- double global_vdot(),global_pdot();
+ double global_v_norm2(),global_p_norm2();
+ double global_div_norm2();
+ void assemble_div_rho_u();
for (m=1;m<=E->sphere.caps_per_proc;m++) {
old_v[m] = (double *)malloc(neq*sizeof(double));
@@ -660,41 +741,48 @@
cycles = E->control.p_iterations;
- residual = 1.0;
+ initial_vel_residual(E, V, P, F,
+ imp * E->monitor.fdotf);
+
+ div_res = 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) {
+ while((relative_err_v >= imp || relative_err_p >= imp) &&
+ (div_res > imp) &&
+ (num_of_loop <= E->control.compress_iter_maxstep)) {
for (m=1;m<=E->sphere.caps_per_proc;m++) {
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);
+ solve_Ahat_p_fhat_CG(E, V, P, F, imp, &cycles);
+ /* compute norm of div(rho*V) */
+ assemble_div_rho_u(E, V, E->u1, lev);
+ div_res = sqrt(global_div_norm2(E, E->u1) / (1e-32 + E->monitor.vdotv));
+
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)) );
+ relative_err_v = sqrt( global_v_norm2(E,diff_v) /
+ (1.0e-32 + E->monitor.vdotv) );
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)) );
+ relative_err_p = sqrt( global_p_norm2(E,diff_p) /
+ (1.0e-32 + E->monitor.pdotp) );
- 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);
+ fprintf(stderr, "itercg -- div(rho*v)/v=%.2e dv/v=%.2e and dp/p=%.2e loop %d\n\n", div_res, relative_err_v, relative_err_p, num_of_loop);
+ fprintf(E->fp, "itercg -- div(rho*v)/v=%.2e dv/v=%.2e and dp/p=%.2e loop %d\n\n", div_res, relative_err_v, relative_err_p, num_of_loop);
}
+ num_of_loop++;
+
} /* end of while */
for (m=1;m<=E->sphere.caps_per_proc;m++) {
@@ -704,36 +792,23 @@
free((void *) diff_p[m]);
}
- return(residual);
+ return;
}
-static double initial_vel_residual(struct All_variables *E,
- double **V, double **P, double **F,
- double imp)
+static void initial_vel_residual(struct All_variables *E,
+ double **V, double **P, double **F,
+ double acc)
{
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, valid;
- 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++)
@@ -749,7 +824,7 @@
/* solve K*u1 = F for u1 */
- valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ valid = solve_del2_u(E, E->u1, F, acc, lev);
if(!valid && (E->parallel.me==0)) {
fputs("Warning: solver not converging! 0\n", stderr);
fputs("Warning: solver not converging! 0\n", E->fp);
@@ -762,31 +837,5 @@
for(i=0; i<neq; i++)
V[m][i] += E->u1[m][i];
- return(v_res);
+ return;
}
-
-
-
-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(r) / 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/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2008-10-29 23:17:11 UTC (rev 13196)
@@ -105,8 +105,12 @@
#define MAX_LEVELS 12 /* max. number of multigrid levels */
#define NCS 14 /* max. number of sphere caps */
-typedef float higher_precision; /* matrix coeffs etc */
-typedef double higher_precision1; /* intermediate calculations for finding above coeffs */
+/* type of elt_del and elt_c arrays */
+#if 1
+ typedef float higher_precision;
+#else
+ typedef double higher_precision;
+#endif
/* Common structures */
@@ -372,17 +376,19 @@
int stop_topo_loop;
int topo_loop;
- float elapsed_time;
- float incompressibility;
- float vdotv;
+ double momentum_residual;
+ double incompressibility;
+ double fdotf;
+ double vdotv;
+ double pdotp;
+
double cpu_time_at_start;
double cpu_time_at_last_cycle;
+ float elapsed_time;
- float T_interior;
- float T_maxvaried;
-
- float T_interior_max_for_exit;
-
+ float T_interior;
+ float T_maxvaried;
+ float T_interior_max_for_exit;
};
struct CONTROL {
@@ -438,8 +444,6 @@
double augmented;
int NASSEMBLE;
- float tole_comp;
-
float sob_tolerance;
/* Rayleigh # */
@@ -458,9 +462,6 @@
/* float adiabaticT0; */
/**/
- float relative_err_accuracy;
-
- /**/
int compress_iter_maxstep;
int self_gravitation; /* self gravitation */
@@ -522,7 +523,6 @@
float ggrd_vtop_omega[4];
#endif
double accuracy;
- double vaccuracy;
char velocity_boundary_file[1000];
char temperature_boundary_file[1000];
char mat_file[1000];
@@ -754,14 +754,14 @@
double *T[NCS],*Tdot[NCS],*buoyancy[NCS];
double *u1[NCS];
double *temp[NCS],*temp1[NCS];
- float *NP[NCS],*Mass[NCS];
- float *MASS[MAX_LEVELS][NCS];
- double *TMass[NCS];
+ double *Mass[NCS], *MASS[MAX_LEVELS][NCS];
+ double *TMass[NCS], *NMass[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 *NP[NCS];
//float *stress[NCS];
float *gstress[NCS];
float *Fas670[NCS],*Fas410[NCS],*Fas670_b[NCS],*Fas410_b[NCS];
Modified: mc/3D/CitcomS/trunk/module/Exchanger/AreaWeightedNormal.cc
===================================================================
--- mc/3D/CitcomS/trunk/module/Exchanger/AreaWeightedNormal.cc 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/module/Exchanger/AreaWeightedNormal.cc 2008-10-29 23:17:11 UTC (rev 13196)
@@ -52,7 +52,7 @@
const Sink& sink,
const All_variables* E) :
size_(boundary.size()),
- toleranceOutflow_(E->control.tole_comp),
+ toleranceOutflow_(E->control.accuracy),
nwght(size_, 0)
{
computeWeightedNormal(boundary, E);
Modified: mc/3D/CitcomS/trunk/module/Exchanger/misc.cc
===================================================================
--- mc/3D/CitcomS/trunk/module/Exchanger/misc.cc 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/module/Exchanger/misc.cc 2008-10-29 23:17:11 UTC (rev 13196)
@@ -158,8 +158,7 @@
void commonE(All_variables *E)
{
- E->control.accuracy = 1e-6;
- E->control.tole_comp = 1e-7;
+ E->control.accuracy = 1e-4;
E->parallel.nprocxy = 1;
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2008-10-29 23:17:11 UTC (rev 13196)
@@ -846,7 +846,6 @@
getIntProperty(properties, "precond", E->control.precondition, fp);
getDoubleProperty(properties, "accuracy", E->control.accuracy, fp);
- getFloatProperty(properties, "tole_compressibility", E->control.tole_comp, fp);
getIntProperty(properties, "mg_cycle", E->control.mg_cycle, fp);
getIntProperty(properties, "down_heavy", E->control.down_heavy, fp);
@@ -867,7 +866,6 @@
if(strcmp(E->control.uzawa, "cg") == 0) {
/* more convergence parameters for "cg" */
getIntProperty(properties, "compress_iter_maxstep", E->control.compress_iter_maxstep, fp);
- getFloatProperty(properties, "relative_err_accuracy", E->control.relative_err_accuracy, fp);
}
}
Modified: mc/3D/CitcomS/trunk/tests/coupled.cfg
===================================================================
--- mc/3D/CitcomS/trunk/tests/coupled.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/tests/coupled.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -77,9 +77,6 @@
datafile = embd
-[CitcomS.esolver.vsolver]
-tole_compressibility = 3e-7
-
[CitcomS.esolver.mesher]
nodex = 17
nodey = 17
Modified: mc/3D/CitcomS/trunk/tests/multicoupled.cfg
===================================================================
--- mc/3D/CitcomS/trunk/tests/multicoupled.cfg 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/tests/multicoupled.cfg 2008-10-29 23:17:11 UTC (rev 13196)
@@ -90,10 +90,6 @@
datafile = embd2
-[CitcomS.esolver1.vsolver]
-tole_compressibility = 6e-7
-
-
[CitcomS.esolver1.mesher]
nodex = 9
nodey = 9
@@ -148,4 +144,12 @@
monitoringFrequency = 1
+[CitcomS.csolver.tsolver]
+monitor_max_T = off
+[CitcomS.esolver1.tsolver]
+monitor_max_T = off
+
+[CitcomS.esolver2.tsolver]
+monitor_max_T = off
+
Modified: mc/3D/CitcomS/trunk/tests/test2.sh
===================================================================
--- mc/3D/CitcomS/trunk/tests/test2.sh 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/tests/test2.sh 2008-10-29 23:17:11 UTC (rev 13196)
@@ -38,11 +38,7 @@
rm $FINE.* $COARSE.*
-
-coupledcitcoms.py \
---launcher.nodegen="n%03d" \
---launcher.nodelist=[101-118,120-129,131-170] \
---launcher.nodes=5 \
+../bin/citcoms --coupled \
--layout.coarse=[0-3] \
--layout.fine=[4] \
--coarse=regional \
@@ -75,8 +71,7 @@
--fine.bc.bottbc=0 \
--fine.bc.bottbcval=0 \
--fine.vsolver.precond=False \
---fine.vsolver.accuracy=1e-9 \
---fine.vsolver.tole_compressibility=1e-6 \
+--fine.vsolver.accuracy=1e-4 \
--fine.vsolver.piterations=8000 \
--fine.vsolver.vlowstep=5000 \
--fine.vsolver.vhighstep=10 \
Modified: mc/3D/CitcomS/trunk/tests/test5.sh
===================================================================
--- mc/3D/CitcomS/trunk/tests/test5.sh 2008-10-29 22:22:02 UTC (rev 13195)
+++ mc/3D/CitcomS/trunk/tests/test5.sh 2008-10-29 23:17:11 UTC (rev 13196)
@@ -37,11 +37,7 @@
rm $FINE.* $COARSE.*
-./coupledcitcoms.sh \
---launcher.nodes=4 \
---launcher.nodegen="n%03d" \
---launcher.nodelist=[171-172,171-172] \
-\
+../bin//citcoms --coupled \
--layout.coarse=[0-1] \
--layout.fine=[2-3] \
\
@@ -76,9 +72,7 @@
--fine.bc.side_sbcs=on \
\
--coarse.vsolver.accuracy=1e-3 \
---coarse.vsolver.tole_compressibility=1e-3 \
--fine.vsolver.accuracy=1e-3 \
---fine.vsolver.tole_compressibility=1e-3 \
\
--fge.excludeTop=on \
--fge.excludeBottom=on \
More information about the CIG-COMMITS
mailing list