[cig-commits] r8250 - in mc/3D/CitcomS/branches/v3.0: .
CitcomS/Components examples examples/Cookbook4
examples/Cookbook5 examples/Cookbook6 examples/Cookbook7
examples/Full examples/Regional lib visual
tan2 at geodynamics.org
tan2 at geodynamics.org
Thu Nov 8 15:28:47 PST 2007
Author: tan2
Date: 2007-11-08 15:28:46 -0800 (Thu, 08 Nov 2007)
New Revision: 8250
Added:
mc/3D/CitcomS/branches/v3.0/visual/project_geoid.c
Modified:
mc/3D/CitcomS/branches/v3.0/
mc/3D/CitcomS/branches/v3.0/ChangeLog
mc/3D/CitcomS/branches/v3.0/CitcomS/Components/Output.py
mc/3D/CitcomS/branches/v3.0/NEWS
mc/3D/CitcomS/branches/v3.0/configure.ac
mc/3D/CitcomS/branches/v3.0/examples/Cookbook4/cookbook4.cfg
mc/3D/CitcomS/branches/v3.0/examples/Cookbook5/cookbook5.cfg
mc/3D/CitcomS/branches/v3.0/examples/Cookbook6/cookbook6.cfg
mc/3D/CitcomS/branches/v3.0/examples/Cookbook7/cookbook7.cfg
mc/3D/CitcomS/branches/v3.0/examples/Full/input.sample
mc/3D/CitcomS/branches/v3.0/examples/Makefile.am
mc/3D/CitcomS/branches/v3.0/examples/Regional/input.sample
mc/3D/CitcomS/branches/v3.0/lib/Element_calculations.c
mc/3D/CitcomS/branches/v3.0/lib/Full_obsolete.c
mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c
mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c
mc/3D/CitcomS/branches/v3.0/lib/Instructions.c
mc/3D/CitcomS/branches/v3.0/lib/Obsolete.c
mc/3D/CitcomS/branches/v3.0/lib/Output.c
mc/3D/CitcomS/branches/v3.0/lib/Pan_problem_misc_functions.c
mc/3D/CitcomS/branches/v3.0/lib/Regional_obsolete.c
mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c
mc/3D/CitcomS/branches/v3.0/lib/Stokes_flow_Incomp.c
mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c
mc/3D/CitcomS/branches/v3.0/visual/Makefile.am
Log:
Merged revisions 8194-8245 via svnmerge from
svn+ssh://svn@geodynamics.org/cig/mc/3D/CitcomS/trunk
........
r8194 | tan2 | 2007-10-30 14:49:58 -0700 (Tue, 30 Oct 2007) | 1 line
Compute d(rho)/dr/rho from rho(r)
........
r8195 | tan2 | 2007-10-30 14:50:52 -0700 (Tue, 30 Oct 2007) | 1 line
Fixed a bug in dimensionalizing density. Provided the formula of geoid calculation in the comments. Rearranged the order of functions.
........
r8196 | tan2 | 2007-10-30 14:53:50 -0700 (Tue, 30 Oct 2007) | 1 line
A post-processing program to project geoid coefficents onto a regular (longitude, latitude) mesh
........
r8197 | tan2 | 2007-10-30 14:54:14 -0700 (Tue, 30 Oct 2007) | 1 line
Added the C program project_geoid to the makefile
........
r8199 | tan2 | 2007-10-30 15:29:44 -0700 (Tue, 30 Oct 2007) | 1 line
Minor modification
........
r8201 | tan2 | 2007-11-01 16:33:30 -0700 (Thu, 01 Nov 2007) | 1 line
Print dv/v=dp/p=1.0 for the 1st Uzawa iteraion
........
r8202 | tan2 | 2007-11-01 16:33:50 -0700 (Thu, 01 Nov 2007) | 1 line
Fixed an error in comment
........
r8204 | tan2 | 2007-11-05 17:03:35 -0800 (Mon, 05 Nov 2007) | 1 line
Scaled topo with variable gravity. Fixed an error in comment. Rearranged computation.
........
r8205 | tan2 | 2007-11-05 17:03:55 -0800 (Mon, 05 Nov 2007) | 1 line
Removed functions related sph. harm in lib/Regional_obsolete.c
........
r8206 | tan2 | 2007-11-05 17:04:20 -0800 (Mon, 05 Nov 2007) | 1 line
Shrank the size of sph. harm arrays
........
r8207 | tan2 | 2007-11-05 17:04:43 -0800 (Mon, 05 Nov 2007) | 1 line
Init'd some variables about vtk_io, which might be accessed with uninit'd values in output_finalize()
........
r8212 | tan2 | 2007-11-06 15:17:54 -0800 (Tue, 06 Nov 2007) | 1 line
Fixed a few memory errors
........
r8213 | tan2 | 2007-11-06 15:18:12 -0800 (Tue, 06 Nov 2007) | 1 line
Increase vlowstep to match the default value in pyre
........
r8214 | tan2 | 2007-11-06 15:18:35 -0800 (Tue, 06 Nov 2007) | 1 line
Removed unused multigrid parameters
........
r8215 | tan2 | 2007-11-06 15:18:54 -0800 (Tue, 06 Nov 2007) | 1 line
Added cgrad solver convergence parameters, increased buoyancy_ratio and lower the # of steps
........
r8226 | tan2 | 2007-11-07 11:51:56 -0800 (Wed, 07 Nov 2007) | 1 line
Print a warning when matrix eqn solver not converging
........
r8227 | tan2 | 2007-11-07 11:52:17 -0800 (Wed, 07 Nov 2007) | 1 line
Removed comp_el from default output, since it is not required for restart anymore.
........
r8228 | tan2 | 2007-11-07 11:52:39 -0800 (Wed, 07 Nov 2007) | 1 line
Decreased the # of processors. This is the only way I can reproduce single-cell convection as in the manual.
........
r8235 | tan2 | 2007-11-08 11:18:26 -0800 (Thu, 08 Nov 2007) | 1 line
Dereased the timestep size to reduce artifacts in advection
........
r8236 | tan2 | 2007-11-08 11:18:52 -0800 (Thu, 08 Nov 2007) | 1 line
Update NEWS
........
r8237 | tan2 | 2007-11-08 11:19:12 -0800 (Thu, 08 Nov 2007) | 1 line
Update the version number
........
r8241 | tan2 | 2007-11-08 13:17:14 -0800 (Thu, 08 Nov 2007) | 1 line
Updated file ChangeLog to r8240
........
r8242 | tan2 | 2007-11-08 13:36:55 -0800 (Thu, 08 Nov 2007) | 1 line
Removed binary checkpoint files from makefile, as the file size is too big for distribution.
........
r8243 | tan2 | 2007-11-08 13:38:09 -0800 (Thu, 08 Nov 2007) | 1 line
Updated file ChangeLog to r8242
........
r8244 | tan2 | 2007-11-08 14:31:21 -0800 (Thu, 08 Nov 2007) | 1 line
Replaced a system call by std C library remove() and disabled another system call (backup input file). Partially fixed issue130. All remaining system calls are in lib/Output_gzdir.c.
........
r8245 | tan2 | 2007-11-08 14:41:31 -0800 (Thu, 08 Nov 2007) | 1 line
Updated file ChangeLog to r8244
........
Property changes on: mc/3D/CitcomS/branches/v3.0
___________________________________________________________________
Name: svnmerge-integrated
- /mc/3D/CitcomS/trunk:1-8172
+ /mc/3D/CitcomS/trunk:1-8172,8194-8245
Modified: mc/3D/CitcomS/branches/v3.0/ChangeLog
===================================================================
--- mc/3D/CitcomS/branches/v3.0/ChangeLog 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/ChangeLog 2007-11-08 23:28:46 UTC (rev 8250)
@@ -1,3 +1,130 @@
+2007-11-08 22:31 tan2
+
+ * lib/Instructions.c, lib/Pan_problem_misc_functions.c: Replaced a
+ system call by std C library remove() and disabled another system
+ call (backup input file). Partially fixed issue130. All remaining
+ system calls are in lib/Output_gzdir.c.
+
+2007-11-08 21:38 tan2
+
+ * ChangeLog: Updated file ChangeLog to r8242
+
+2007-11-08 21:36 tan2
+
+ * examples/Makefile.am: Removed binary checkpoint files from
+ makefile, as the file size is too big for distribution.
+
+2007-11-08 21:17 tan2
+
+ * ChangeLog: Updated file ChangeLog to r8240
+
+2007-11-08 19:19 tan2
+
+ * configure.ac: Update the version number
+
+2007-11-08 19:18 tan2
+
+ * NEWS: Update NEWS
+
+2007-11-08 19:18 tan2
+
+ * examples/Cookbook5/cookbook5.cfg: Dereased the timestep size to
+ reduce artifacts in advection
+
+2007-11-07 19:52 tan2
+
+ * examples/Cookbook4/cookbook4.cfg: Decreased the # of processors.
+ This is the only way I can reproduce single-cell convection as in
+ the manual.
+
+2007-11-07 19:52 tan2
+
+ * CitcomS/Components/Output.py, lib/Output.c: Removed comp_el from
+ default output, since it is not required for restart anymore.
+
+2007-11-07 19:51 tan2
+
+ * lib/Stokes_flow_Incomp.c: Print a warning when matrix eqn solver
+ not converging
+
+2007-11-06 23:18 tan2
+
+ * examples/Cookbook7/cookbook7.cfg: Added cgrad solver convergence
+ parameters, increased buoyancy_ratio and lower the # of steps
+
+2007-11-06 23:18 tan2
+
+ * examples/Cookbook6/cookbook6.cfg: Removed unused multigrid
+ parameters
+
+2007-11-06 23:18 tan2
+
+ * examples/Full/input.sample, examples/Regional/input.sample:
+ Increase vlowstep to match the default value in pyre
+
+2007-11-06 23:17 tan2
+
+ * visual/project_geoid.c: Fixed a few memory errors
+
+2007-11-06 01:04 tan2
+
+ * lib/Output.c: Init'd some variables about vtk_io, which might be
+ accessed with uninit'd values in output_finalize()
+
+2007-11-06 01:04 tan2
+
+ * lib/Full_obsolete.c, lib/Global_operations.c, lib/Obsolete.c,
+ lib/Sphere_harmonics.c, lib/Topo_gravity.c: Shrank the size of
+ sph. harm arrays
+
+2007-11-06 01:03 tan2
+
+ * lib/Regional_obsolete.c: Removed functions related sph. harm in
+ lib/Regional_obsolete.c
+
+2007-11-06 01:03 tan2
+
+ * lib/Topo_gravity.c: Scaled topo with variable gravity. Fixed an
+ error in comment. Rearranged computation.
+
+2007-11-01 23:33 tan2
+
+ * lib/Full_version_dependent.c: Fixed an error in comment
+
+2007-11-01 23:33 tan2
+
+ * lib/Stokes_flow_Incomp.c: Print dv/v=dp/p=1.0 for the 1st Uzawa
+ iteraion
+
+2007-10-30 22:29 tan2
+
+ * visual/Makefile.am: Minor modification
+
+2007-10-30 21:54 tan2
+
+ * visual/Makefile.am: Added the C program project_geoid to the
+ makefile
+
+2007-10-30 21:53 tan2
+
+ * visual/project_geoid.c: A post-processing program to project
+ geoid coefficents onto a regular (longitude, latitude) mesh
+
+2007-10-30 21:50 tan2
+
+ * lib/Topo_gravity.c: Fixed a bug in dimensionalizing density.
+ Provided the formula of geoid calculation in the comments.
+ Rearranged the order of functions.
+
+2007-10-30 21:49 tan2
+
+ * lib/Element_calculations.c, lib/Instructions.c: Compute
+ d(rho)/dr/rho from rho(r)
+
+2007-10-22 20:57 tan2
+
+ * ChangeLog: Updated file ChangeLog to r8170
+
2007-10-22 20:26 tan2
* NEWS: Updated NEWS
Modified: mc/3D/CitcomS/branches/v3.0/CitcomS/Components/Output.py
===================================================================
--- mc/3D/CitcomS/branches/v3.0/CitcomS/Components/Output.py 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/CitcomS/Components/Output.py 2007-11-08 23:28:46 UTC (rev 8250)
@@ -53,7 +53,7 @@
"ascii-gz",
"hdf5"]))
output_optional = inv.str("output_optional",
- default="surf,botm,tracer,comp_el")
+ default="surf,botm,tracer")
# experimental vtk output
gzdir_vtkio = inv.int("gzdir_vtkio", default=0)
Modified: mc/3D/CitcomS/branches/v3.0/NEWS
===================================================================
--- mc/3D/CitcomS/branches/v3.0/NEWS 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/NEWS 2007-11-08 23:28:46 UTC (rev 8250)
@@ -2,41 +2,53 @@
Please send CitcomS.py bug reports to cig-mc at geodynamics.org.
-Version 3.0.0b
+Version 3.0.0
-* Two implementations of compressible convection under Truncated An-elastic
- Liquid Approximation are available, one by Wei Leng and Shijie Zhong and
- the other by Eh Tan. See Cookbook 8 in the manual for more details.
+This release of CitcomS (3.0) contains many new features, including:
-* The computation can be resumed from previous checkpoints. See Cookbook 8
- in the manual for more details.
+* two implementations of compressible convection, one by Wei Leng and Shijie
+ Zhong and the other by Eh Tan (Cookbook 8);
-* The chemical convection can handle multi-component composition.
+* the ability to resume computation from previous checkpoints (Cookbook 8);
-* The non-Newtonian Stokes' solver was not working properly and is fixed.
+* multi-component chemical convection;
-* The viscosity can be composition dependent or pseudo-plastic (contributed
- by Thorsten Becker).
+* a fixed non-Newtonian solver;
-* The ASCII output can be compressed on-the-fly (contributed by Thorsten
- Becker).
+* an exchanger package for solver coupling (Cookbook 9);
-* An easier way for radial mesh refinement is available (contributed by
- Thorsten Becker). See Cookbook 9 in the manual for more details.
+* removing the rigid body rotation component from the velocity by Shijie
+ Zhong (Cookbook 8);
-* An exchanger package for nested solver coupling is available. See Cookbook
- 9 in the manual for more details.
+* an option to disable monitoring of maximum temperature (Cookbook 9);
-* The rigid body rotation component is removed from the velocity (contributed
- by Shijie Zhong).
+* a rheology option for pseudo-plasiticity, composition dependent viscosity
+ and heat generation, compressed ASCII output, and an easier way for mesh
+ refinement for the radial coordinate by Thorsten Becker.
-* A new parameter can disable monitoring of maximum temperature. See Cookbook
- 9 in the manual for more details.
-* Several other changes are not backward compatible. A few of them will
- affect the solution, and others are changes in input parameters only. See
- Section 1.2 in the manual for the details.
+There are additional seven backward-incompatible changes. Among them, the
+first four will affect the results. The same input file will produce slightly
+different results in v3.0.0 than in v2.2.2.
+* the viscosity field at element level is not smoothed (this might slow down
+ the convergence but will represent the viscosity field more accurately);
+
+* the Lenardic filter on temperature is disabled by default;
+
+* the rigid body rotation component is removed from the velocity by default;
+
+* use of a better pseudo-random number generator to generate the initial
+ tracer;
+
+* the type of input parameter coor is changed from a boolean to an integer;
+
+* setting restart=on will resume the computation from the checkpoint files and
+ will not need the tracer files (the old way of reading initial temperature
+ from velo files can be achieved by tic_method=-1);
+
+* the input parameter reset_initial_composition becomes obsolete.
+
Version 2.2.2
Modified: mc/3D/CitcomS/branches/v3.0/configure.ac
===================================================================
--- mc/3D/CitcomS/branches/v3.0/configure.ac 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/configure.ac 2007-11-08 23:28:46 UTC (rev 8250)
@@ -25,7 +25,7 @@
# $Id$
AC_PREREQ(2.59)
-AC_INIT([CitcomS], [3.0.0b], [cig-mc at geodynamics.org], [CitcomS])
+AC_INIT([CitcomS], [3.0.0], [cig-mc at geodynamics.org], [CitcomS])
AC_CONFIG_AUX_DIR([./aux-config])
AC_CONFIG_SRCDIR([bin/Citcom.c])
AC_CONFIG_HEADER([config.h])
Modified: mc/3D/CitcomS/branches/v3.0/examples/Cookbook4/cookbook4.cfg
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Cookbook4/cookbook4.cfg 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Cookbook4/cookbook4.cfg 2007-11-08 23:28:46 UTC (rev 8250)
@@ -5,7 +5,7 @@
[CitcomS.controller]
-monitoringFrequency = 10 ; how often outputs are created
+monitoringFrequency = 50 ; how often outputs are created
[CitcomS.solver]
@@ -29,9 +29,9 @@
coor = 1
coor_file = coor.dat
-nprocx = 4
-nprocy = 2
-nprocz = 2
+nprocx = 1
+nprocy = 1
+nprocz = 1
nodex = 33
nodey = 17
nodez = 17
Modified: mc/3D/CitcomS/branches/v3.0/examples/Cookbook5/cookbook5.cfg
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Cookbook5/cookbook5.cfg 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Cookbook5/cookbook5.cfg 2007-11-08 23:28:46 UTC (rev 8250)
@@ -1,7 +1,7 @@
# Cookbook 5: Subduction Models with Trench Rollback
[CitcomS]
-steps = 1000 ; number of time steps
+steps = 1100 ; number of time steps
[CitcomS.controller]
monitoringFrequency = 100 ; how often outputs are created
@@ -42,6 +42,10 @@
fi_max = 0.5
radius_inner = 0.7
+[CitcomS.solver.tsolver]
+finetunedt = 0.75
+monitor_max_T = on
+
[CitcomS.solver.visc]
num_mat = 4
visc0 = 100,0.003,1,2
Modified: mc/3D/CitcomS/branches/v3.0/examples/Cookbook6/cookbook6.cfg
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Cookbook6/cookbook6.cfg 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Cookbook6/cookbook6.cfg 2007-11-08 23:28:46 UTC (rev 8250)
@@ -32,11 +32,7 @@
precond = on
accuracy = 1e-10
tole_compressibility = 1e-06
-mg_cycle = 1
-down_heavy = 5
-up_heavy = 5
vlowstep = 100000
-vhighstep = 3
piterations = 100000
Modified: mc/3D/CitcomS/branches/v3.0/examples/Cookbook7/cookbook7.cfg
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Cookbook7/cookbook7.cfg 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Cookbook7/cookbook7.cfg 2007-11-08 23:28:46 UTC (rev 8250)
@@ -2,11 +2,11 @@
[CitcomS]
solver = full
-steps = 30 ; number of time steps
+steps = 15 ; number of time steps
[CitcomS.controller]
-monitoringFrequency = 10 ; how often outputs are created
+monitoringFrequency = 5 ; how often outputs are created
[CitcomS.solver]
@@ -39,12 +39,20 @@
chemical_buoyancy = 1
buoy_type = 1
-buoyancy_ratio = 0.4
+buoyancy_ratio = 0.5
regular_grid_deltheta = 1.0
regular_grid_delphi = 1.0
+[CitcomS.solver.vsolver]
+Solver = cgrad
+accuracy = 1e-06
+vlowstep = 1000
+tole_compressibility = 1e-07
+piterations = 1000
+
+
# Assign the viscosities.
[CitcomS.solver.visc]
VISC_UPDATE = on
Modified: mc/3D/CitcomS/branches/v3.0/examples/Full/input.sample
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Full/input.sample 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Full/input.sample 2007-11-08 23:28:46 UTC (rev 8250)
@@ -206,7 +206,7 @@
mg_cycle=1
down_heavy=3
up_heavy=3
-vlowstep=500
+vlowstep=1000
vhighstep=3
piterations=375
Modified: mc/3D/CitcomS/branches/v3.0/examples/Makefile.am
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Makefile.am 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Makefile.am 2007-11-08 23:28:46 UTC (rev 8250)
@@ -98,18 +98,6 @@
Cookbook7/cookbook7.cfg \
Cookbook8/cookbook8.cfg \
Cookbook8/coord.dat \
- Cookbook8/restart/cookbook8.chkpt.0.9000 \
- Cookbook8/restart/cookbook8.chkpt.1.9000 \
- Cookbook8/restart/cookbook8.chkpt.2.9000 \
- Cookbook8/restart/cookbook8.chkpt.3.9000 \
- Cookbook8/restart/cookbook8.chkpt.4.9000 \
- Cookbook8/restart/cookbook8.chkpt.5.9000 \
- Cookbook8/restart/cookbook8.chkpt.6.9000 \
- Cookbook8/restart/cookbook8.chkpt.7.9000 \
- Cookbook8/restart/cookbook8.chkpt.8.9000 \
- Cookbook8/restart/cookbook8.chkpt.9.9000 \
- Cookbook8/restart/cookbook8.chkpt.10.9000 \
- Cookbook8/restart/cookbook8.chkpt.11.9000 \
Cookbook9/cookbook9.cfg \
Cookbook9/ic/cntn.velo.0.0 \
Cookbook9/ic/cntn.velo.1.0 \
Modified: mc/3D/CitcomS/branches/v3.0/examples/Regional/input.sample
===================================================================
--- mc/3D/CitcomS/branches/v3.0/examples/Regional/input.sample 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/examples/Regional/input.sample 2007-11-08 23:28:46 UTC (rev 8250)
@@ -209,7 +209,7 @@
mg_cycle=1
down_heavy=3
up_heavy=3
-vlowstep=500
+vlowstep=1000
vhighstep=3
piterations=375
Modified: mc/3D/CitcomS/branches/v3.0/lib/Element_calculations.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Element_calculations.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Element_calculations.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -828,13 +828,13 @@
{
void get_global_shape_fn();
void construct_c3x3matrix_el();
- int p, a, i;
- double temp, beta, x[4];
+ int p, a, i, j, nz;
+ double temp, beta, rho_avg, x[4];
struct Shape_function GN;
struct Shape_function_dx GNx;
struct Shape_function_dA dOmega;
- double rtf[4][9];
+ double rtf[4][9], rho[9];
const int dims = E->mesh.nsd;
const int ends = enodes[dims];
@@ -846,25 +846,54 @@
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 */
+ switch (E->refstate.choice) {
+ case 1:
+ /* the reference state is computed by rho=exp((1-r)Di/gamma) */
+ /* so d(rho)/dr/rho == -Di/gamma */
- 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
+ beta = - E->control.disptn_number * E->control.inv_gruneisen;
+
+ for(a=1;a<=ends;a++) {
+ for (i=1;i<=dims;i++) {
+ x[i] = E->N.ppt[GNPINDEX(a,1)]
+ * E->element_Cc.ppt[BPINDEX(3,i,a,1)];
+ }
+ p=dims*(a-1);
+ elt_c[p ][0] = -x[1] * temp * beta;
+ elt_c[p+1][0] = -x[2] * temp * beta;
+ elt_c[p+2][0] = -x[3] * temp * beta;
}
- 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;
+ break;
+ default:
+ /* compute d(rho)/dr/rho from rho(r) */
+
+ for(a=1;a<=ends;a++) {
+ j = E->IEN[lev][m][el].node[a];
+ nz = (j - 1) % E->lmesh.noz + 1;
+ rho[a] = E->refstate.rho[nz];
+ }
+
+ rho_avg = 0;
+ for(a=1;a<=ends;a++) {
+ rho_avg += rho[a];
+ }
+ rho_avg /= ends;
+
+ for(a=1;a<=ends;a++) {
+ for (i=1;i<=dims;i++) {
+ x[i] = rho[a] * GNx.ppt[GNPXINDEX(2,a,1)]
+ * E->N.ppt[GNPINDEX(a,1)]
+ * E->element_Cc.ppt[BPINDEX(3,i,a,1)];
+ }
+ p=dims*(a-1);
+ elt_c[p ][0] = -x[1] * temp / rho_avg;
+ elt_c[p+1][0] = -x[2] * temp / rho_avg;
+ elt_c[p+2][0] = -x[3] * temp / rho_avg;
+ }
+
}
+
return;
}
Modified: mc/3D/CitcomS/branches/v3.0/lib/Full_obsolete.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Full_obsolete.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Full_obsolete.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -423,12 +423,12 @@
if (E->parallel.nprocz==1) return;
jumpp = E->sphere.hindice;
- nsl = E->sphere.hindice*2+3;
+ nsl = E->sphere.hindice*2;
me = E->parallel.me;
- TG = ( float *)malloc((nsl+1)*sizeof(float));
+ TG = ( float *)malloc(nsl*sizeof(float));
if (E->parallel.me_loc[3]==dest_proc)
- RG = ( float *)malloc((nsl+1)*sizeof(float));
+ RG = ( float *)malloc(nsl*sizeof(float));
for (i=0;i<E->sphere.hindice;i++) {
TG[i] = sphc[i];
@@ -481,11 +481,11 @@
if (E->parallel.nprocxy==1) return;
- nsl = E->sphere.hindice*2+2;
+ nsl = E->sphere.hindice*2;
me = E->parallel.me;
for (i=1;i<E->parallel.nprocxy*E->parallel.surf_proc_per_cap;i++)
- RG[i] = ( float *)malloc((nsl+1)*sizeof(float));
+ RG[i] = ( float *)malloc(nsl*sizeof(float));
idb=0;
@@ -524,7 +524,6 @@
return;
}
-#endif
Modified: mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Full_version_dependent.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -153,7 +153,7 @@
case 0:
/* generate uniform mesh in radial direction */
dr = (E->sphere.ro-E->sphere.ri)/(E->mesh.noz-1);
-
+
for (k=1;k <= E->mesh.noz;k++) {
rr[k] = E->sphere.ri + (k-1)*dr;
}
@@ -170,7 +170,7 @@
fscanf(fp1,"%d %f",&nn,&tt1);
rr[k]=tt1;
}
-
+
fclose(fp1);
break;
case 2:
@@ -185,7 +185,7 @@
myerror(E,"coor flag undefined in Full_version_dependent");
break;
}
-
+
for (i=1;i<=E->lmesh.noz;i++) {
k = E->lmesh.nzs+i-1;
RR[i] = rr[k];
@@ -325,20 +325,20 @@
break;
case 3:
-
- /* set up a linear temperature profile first */
+
+ /* set up a conductive temperature profile first */
for(m=1;m<=E->sphere.caps_per_proc;m++)
for(i=1;i<=noy;i++)
for(j=1;j<=nox;j++)
for(k=1;k<=noz;k++) {
node=k+(j-1)*noz+(i-1)*nox*noz;
r1=E->sx[m][3][node];
- E->T[m][node] = (E->control.TBCtopval*E->sphere.ro
- - E->control.TBCbotval*E->sphere.ri)
- / (E->sphere.ro - E->sphere.ri)
- + (E->control.TBCbotval - E->control.TBCtopval)
- * E->sphere.ro * E->sphere.ri / r1
- / (E->sphere.ro - E->sphere.ri);
+ E->T[m][node] = (E->control.TBCtopval*E->sphere.ro
+ - E->control.TBCbotval*E->sphere.ri)
+ / (E->sphere.ro - E->sphere.ri)
+ + (E->control.TBCbotval - E->control.TBCtopval)
+ * E->sphere.ro * E->sphere.ri / r1
+ / (E->sphere.ro - E->sphere.ri);
}
/* This part put a temperature anomaly for whole mantle. The horizontal
pattern of the anomaly is given by spherical harmonic ll & mm. */
Modified: mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Global_operations.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -452,12 +452,12 @@
int jumpp,total,j,d;
float *sphcs,*temp;
- temp = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
- sphcs = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
+ temp = (float *) malloc((E->sphere.hindice*2)*sizeof(float));
+ sphcs = (float *) malloc((E->sphere.hindice*2)*sizeof(float));
/* pack */
jumpp = E->sphere.hindice;
- total = E->sphere.hindice*2+3;
+ total = E->sphere.hindice*2;
for (j=0;j<E->sphere.hindice;j++) {
sphcs[j] = sphc[j];
sphcs[j+jumpp] = sphs[j];
@@ -758,7 +758,7 @@
float *sphcs,*temp;
if (E->parallel.nprocz > 1) {
- total = E->sphere.hindice*2+3;
+ total = E->sphere.hindice*2;
temp = (float *) malloc(total*sizeof(float));
sphcs = (float *) malloc(total*sizeof(float));
Modified: mc/3D/CitcomS/branches/v3.0/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Instructions.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Instructions.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -137,14 +137,14 @@
(E->solver.parallel_communication_routs_v)(E);
(E->solver.parallel_communication_routs_s)(E);
+ reference_state(E);
+
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);
-
/* this matrix results from spherical geometry */
/* construct_c3x3matrix(E); */
@@ -1426,7 +1426,7 @@
if((E->output.gzdir.vtk_io == 3)||(E->parallel.me == 0)){
/* delete the geo files */
get_vtk_filename(files,1,E,0);
- sprintf(message,"rm %s",files);system(message);
+ remove(files);
if(E->parallel.me == 0){
/* close the log */
fclose(E->output.gzdir.vtk_fp);
Modified: mc/3D/CitcomS/branches/v3.0/lib/Obsolete.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Obsolete.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Obsolete.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -333,7 +333,7 @@
double t,f;
double dth,dfi,sqrt_multis();
- E->sphere.con = (double *)malloc((E->sphere.hindice+3)*sizeof(double));
+ E->sphere.con = (double *)malloc(E->sphere.hindice*sizeof(double));
for (ll=0;ll<=E->output.llmax;ll++)
for (mm=0;mm<=ll;mm++) {
E->sphere.con[E->sphere.hindex[ll][mm]] =
@@ -344,7 +344,7 @@
E->sphere.tablenplm = (double **) malloc((E->sphere.nox+1)
*sizeof(double*));
for (i=1;i<=E->sphere.nox;i++)
- E->sphere.tablenplm[i]= (double *)malloc((E->sphere.hindice+3)
+ E->sphere.tablenplm[i]= (double *)malloc(E->sphere.hindice
*sizeof(double));
E->sphere.tablencosf = (double **) malloc((E->sphere.noy+1)
Modified: mc/3D/CitcomS/branches/v3.0/lib/Output.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Output.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Output.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -63,9 +63,11 @@
int m = E->parallel.me;
input_string("output_format", E->output.format, "ascii",m);
- input_string("output_optional", E->output.optional, "surf,botm,tracer,comp_el",m);
+ input_string("output_optional", E->output.optional, "surf,botm,tracer",m);
/* gzdir type of I/O */
+ E->output.gzdir.vtk_io = 0;
+ E->output.gzdir.rnr = 0;
if(strcmp(E->output.format, "ascii-gz") == 0){
/*
vtk_io = 1: write files for post-processing into VTK
Modified: mc/3D/CitcomS/branches/v3.0/lib/Pan_problem_misc_functions.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Pan_problem_misc_functions.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Pan_problem_misc_functions.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -80,7 +80,10 @@
if (E->parallel.me==0) {
sprintf(unique_name,"%06d.%s-%s",E->control.PID,comment,name);
sprintf(command,"cp -f %s %s\n",name,unique_name);
+#if 0
+ /* disable copying file, since some MPI implementation doesn't support it */
system(command);
+#endif
}
}
@@ -459,7 +462,7 @@
{
int k,klim,nb,nt,nm;
double drb,dr0,drt,dr,drm,range,r,mrange, brange,bfrac,trange, tfrac;
-
+
brange = (double)E->control.coor_refine[0];
bfrac = (double)E->control.coor_refine[1];
trange = (double)E->control.coor_refine[2];
@@ -496,7 +499,7 @@
rr[k] = r;
}
}
-/*
+/*
get r spacing at radial locations and node numbers as specified
@@ -530,10 +533,10 @@
rr[1] = E->control.rrlayer[0];
for(j = 1; j < E->control.rlayers; j++){
- ddr = (E->control.rrlayer[j] - E->control.rrlayer[j - 1]) /
+ ddr = (E->control.rrlayer[j] - E->control.rrlayer[j - 1]) /
(E->control.nrlayer[j] - E->control.nrlayer[j - 1]);
for(k = E->control.nrlayer[j-1]+1;k <= E->control.nrlayer[j];k++)
rr[k] = rr[k-1]+ddr;
}
-
+
}
Modified: mc/3D/CitcomS/branches/v3.0/lib/Regional_obsolete.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Regional_obsolete.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Regional_obsolete.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -410,159 +410,7 @@
}
-void sum_across_depth_sph(E,sphc,sphs,dest_proc)
- struct All_variables *E;
-int dest_proc;
-float *sphc,*sphs;
-{
- int jumpp,i,j,nsl,idb,to_proc,from_proc,mst,me;
-
- float *RG,*TG;
-
- MPI_Status status[100];
- MPI_Status status1;
- MPI_Request request[100];
-
- if (E->parallel.nprocz==1) return;
-
- jumpp = E->sphere.hindice;
- nsl = E->sphere.hindice*2+3;
- me = E->parallel.me;
-
- TG = ( float *)malloc((nsl+1)*sizeof(float));
- if (E->parallel.me_loc[3]==dest_proc)
- RG = ( float *)malloc((nsl+1)*sizeof(float));
-
- for (i=0;i<E->sphere.hindice;i++) {
- TG[i] = sphc[i];
- TG[i+jumpp] = sphs[i];
- }
-
-
- if (E->parallel.me_loc[3]!=dest_proc) { /* send TG */
- to_proc = E->parallel.me_sph*E->parallel.nprocz+E->parallel.nprocz-1;
- mst = me;
- MPI_Send(TG,nsl,MPI_FLOAT,to_proc,mst,E->parallel.world);
- }
-
- parallel_process_sync(E);
-
- if (E->parallel.me_loc[3]==dest_proc) {
- for (i=1;i<E->parallel.nprocz;i++) {
- from_proc = me - i;
- mst = from_proc;
- MPI_Recv(RG,nsl,MPI_FLOAT,from_proc,mst,E->parallel.world,&status1);
-
- for (j=0;j<E->sphere.hindice;j++) {
- sphc[j] += RG[j];
- sphs[j] += RG[j+jumpp];
- }
- }
- }
-
- free((void *) TG);
- if (E->parallel.me_loc[3]==dest_proc)
- free((void *) RG);
-
- return;
-}
-
-
-void sum_across_surf_sph(E,TG,loc_proc)
- struct All_variables *E;
-int loc_proc;
-float *TG;
-{
-
- int i,j,nsl,idb,to_everyone,from_proc,mst,me;
-
- float *RG[20];
-
- MPI_Status status[100];
- MPI_Status status1;
- MPI_Request request[100];
-
- if (E->parallel.nprocxy==1) return;
-
- nsl = E->sphere.hindice*2+2;
- me = E->parallel.me;
-
- for (i=1;i<E->parallel.nprocxy;i++)
- RG[i] = ( float *)malloc((nsl+1)*sizeof(float));
-
-
- idb=0;
- for (i=1;i<=E->parallel.nprocxy;i++) {
- to_everyone = E->parallel.nprocz*(i-1) + loc_proc;
-
- if (me!=to_everyone) { /* send TG */
- idb++;
- mst = me;
- MPI_Isend(TG,nsl,MPI_FLOAT,to_everyone,mst,E->parallel.world,&request[idb-1]);
- }
- }
-
- /* parallel_process_sync(E); */
-
- idb=0;
- for (i=1;i<=E->parallel.nprocxy;i++) {
- from_proc = E->parallel.nprocz*(i-1) + loc_proc;
- if (me!=from_proc) { /* me==0 receive all TG and add them up */
- mst = from_proc;
- idb++;
- MPI_Irecv(RG[idb],nsl,MPI_FLOAT,from_proc,mst,E->parallel.world,&request[idb-1]);
- }
- }
-
- MPI_Waitall(idb,request,status);
-
- for (i=1;i<E->parallel.nprocxy;i++)
- for (j=0;j<nsl; j++) {
- TG[j] += RG[i][j];
- }
-
- /* parallel_process_sync(E); */
-
- for (i=1;i<E->parallel.nprocxy;i++)
- free((void *) RG[i]);
-
- return;
-}
-
-
-void set_communication_sphereh(E)
- struct All_variables *E;
-{
- int i;
-
- i = cases[E->sphere.caps_per_proc];
-
- E->parallel.nproc_sph[1] = incases3[i].xy[0];
- E->parallel.nproc_sph[2] = incases3[i].xy[1];
-
- E->sphere.lelx = E->sphere.elx/E->parallel.nproc_sph[1];
- E->sphere.lely = E->sphere.ely/E->parallel.nproc_sph[2];
- E->sphere.lsnel = E->sphere.lely*E->sphere.lelx;
- E->sphere.lnox = E->sphere.lelx + 1;
- E->sphere.lnoy = E->sphere.lely + 1;
- E->sphere.lnsf = E->sphere.lnox*E->sphere.lnoy;
-
- for (i=0;i<=E->parallel.nprocz-1;i++)
- if (E->parallel.me_loc[3] == i) {
- E->parallel.me_sph = (E->parallel.me-i)/E->parallel.nprocz;
- E->parallel.me_loc_sph[1] = E->parallel.me_sph%E->parallel.nproc_sph[1];
- E->parallel.me_loc_sph[2] = E->parallel.me_sph/E->parallel.nproc_sph[1];
- }
-
- E->sphere.lexs = E->sphere.lelx * E->parallel.me_loc_sph[1];
- E->sphere.leys = E->sphere.lely * E->parallel.me_loc_sph[2];
-
- return;
-}
-
-
-
/* ========================================================== */
/* from Boundary_conditions.c */
/* =========================================================== */
@@ -967,43 +815,6 @@
/* ========================================================== */
-/* from Global_operations.c */
-/* =========================================================== */
-
-void sum_across_depth_sph1(E,sphc,sphs)
-struct All_variables *E;
-float *sphc,*sphs;
-{
- int jumpp,total,j,d;
- float *sphcs,*temp;
-
- temp = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
- sphcs = (float *) malloc((E->sphere.hindice*2+3)*sizeof(float));
-
- total = E->sphere.hindice*2+3;
- jumpp = E->sphere.hindice;
- for (j=0;j<E->sphere.hindice;j++) {
- sphcs[j] = sphc[j];
- sphcs[j+jumpp] = sphs[j];
- }
-
- MPI_Allreduce(sphcs,temp,total,MPI_FLOAT,MPI_SUM,E->parallel.vertical_comm);
-
- for (j=0;j<E->sphere.hindice;j++) {
- sphc[j] = temp[j];
- sphs[j] = temp[j+jumpp];
- }
-
- MPI_Allreduce(sphcs,temp,total,MPI_FLOAT,MPI_SUM,E->parallel.vertical_comm);
-
- free((void*) temp);
- free((void*) sphcs);
-
-return;
-}
-
-
-/* ========================================================== */
/* from */
/* =========================================================== */
Modified: mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Sphere_harmonics.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -30,10 +30,10 @@
/* spherical harmonic coeff (0=cos, 1=sin)
for surface topo, cmb topo and geoid */
for (i=0;i<=1;i++) {
- E->sphere.harm_geoid[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_geoid_from_bncy[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_geoid_from_tpgt[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
- E->sphere.harm_geoid_from_tpgb[i]=(float*)malloc((E->sphere.hindice+2)*sizeof(float));
+ E->sphere.harm_geoid[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
+ E->sphere.harm_geoid_from_bncy[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
+ E->sphere.harm_geoid_from_tpgt[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
+ E->sphere.harm_geoid_from_tpgb[i]=(float*)malloc(E->sphere.hindice*sizeof(float));
}
compute_sphereh_table(E);
@@ -114,9 +114,9 @@
E->sphere.tablessinf[m] = (double **) malloc((E->lmesh.nsf+1)*sizeof(double*));
for (i=1;i<=E->lmesh.nsf;i++) {
- E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice+3)*sizeof(double));
- E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
- E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+3)*sizeof(double));
+ E->sphere.tablesplm[m][i]= (double *)malloc((E->sphere.hindice)*sizeof(double));
+ E->sphere.tablescosf[m][i]= (double *)malloc((E->output.llmax+1)*sizeof(double));
+ E->sphere.tablessinf[m][i]= (double *)malloc((E->output.llmax+1)*sizeof(double));
}
}
Modified: mc/3D/CitcomS/branches/v3.0/lib/Stokes_flow_Incomp.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Stokes_flow_Incomp.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Stokes_flow_Incomp.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -103,6 +103,27 @@
/* ========================================================================= */
+static void print_convergence_progress(struct All_variables *E,
+ int count, double time0,
+ double sq_vdotv,
+ double dvelocity, double dpressure)
+{
+ double CPU_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);
+
+ return;
+}
+
+
+
static float solve_Ahat_p_fhat(struct All_variables *E,
double **V, double **P, double **F,
double imp, int *steps_max)
@@ -209,27 +230,20 @@
sq_vdotv = sqrt(E->monitor.vdotv);
+ /* pressure and velocity corrections */
+ dpressure = 1.0;
+ dvelocity = 1.0;
+
+
if (E->control.print_convergence && E->parallel.me==0) {
- 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,
- 0.0, 0.0, E->monitor.solution_cycles);
- fflush(E->fp);
- 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,
- 0.0, 0.0, E->monitor.solution_cycles);
+ print_convergence_progress(E, count, time0, sq_vdotv,
+ dvelocity, dpressure);
}
- /* pressure and velocity corrections */
- dpressure = 1.0;
- dvelocity = 1.0;
-
- valid = 1;
r0dotz0 = 0;
- while( (valid) && (count < *steps_max) &&
+ while( (count < *steps_max) &&
(E->monitor.incompressibility >= E->control.tole_comp) &&
(dpressure >= imp) && (dvelocity >= imp) ) {
@@ -262,6 +276,10 @@
/* solve K*u1 = grad(s2) for u1 */
assemble_grad_p(E, s2, F, lev);
valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ if(!valid && (E->parallel.me==0)) {
+ fputs("Warning: solver not converging! 1\n", stderr);
+ fputs("Warning: solver not converging! 1\n", E->fp);
+ }
strip_bcs_from_residual(E, E->u1, lev);
@@ -309,15 +327,9 @@
sq_vdotv = sqrt(E->monitor.vdotv);
- if(E->control.print_convergence && E->parallel.me==0) {
- 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);
+ if (E->control.print_convergence && E->parallel.me==0) {
+ print_convergence_progress(E, count, time0, sq_vdotv,
+ dvelocity, dpressure);
}
@@ -435,28 +447,21 @@
sq_vdotv = sqrt(E->monitor.vdotv);
+ /* pressure and velocity corrections */
+ dpressure = 1.0;
+ dvelocity = 1.0;
+
if (E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp,"AhatP (%03d) after %g 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,
- 0.0, 0.0, E->monitor.solution_cycles);
- fflush(E->fp);
- fprintf(stderr,"AhatP (%03d) after %g 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,
- 0.0, 0.0, E->monitor.solution_cycles);
+ print_convergence_progress(E, count, time0, sq_vdotv,
+ dvelocity, dpressure);
}
- /* pressure and velocity corrections */
- dpressure = 1.0;
- dvelocity = 1.0;
-
valid = 1;
r0dotrt = alpha = omega = 0;
- while( (valid) && (count < *steps_max) &&
+ while( (count < *steps_max) &&
((E->monitor.incompressibility >= E->control.tole_comp) &&
(dpressure >= imp) && (dvelocity >= imp)) ) {
@@ -496,7 +501,10 @@
/* 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");
+ if(!valid && (E->parallel.me==0)) {
+ fputs("Warning: solver not converging! 1\n", stderr);
+ fputs("Warning: solver not converging! 1\n", E->fp);
+ }
strip_bcs_from_residual(E, u0, lev);
@@ -523,7 +531,10 @@
/* solve K*u1 = grad(st) for u1 */
assemble_grad_p(E, st, F, lev);
valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
- if(!valid) fprintf(stderr, "not valid 2\n");
+ if(!valid && (E->parallel.me==0)) {
+ fputs("Warning: solver not converging! 2\n", stderr);
+ fputs("Warning: solver not converging! 2\n", E->fp);
+ }
strip_bcs_from_residual(E, E->u1, lev);
@@ -573,18 +584,12 @@
count++;
-
+
sq_vdotv = sqrt(E->monitor.vdotv);
if(E->control.print_convergence && E->parallel.me==0) {
- fprintf(E->fp, "AhatP (%03d) after %g 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 %g 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);
+ print_convergence_progress(E, count, time0, sq_vdotv,
+ dvelocity, dpressure);
}
@@ -718,7 +723,7 @@
int neq = E->lmesh.neq;
int gneq = E->mesh.neq;
int lev = E->mesh.levmax;
- int i, m;
+ int i, m, valid;
double v_res;
v_res = sqrt(global_vdot(E, F, F, lev) / gneq);
@@ -746,7 +751,11 @@
/* solve K*u1 = F for u1 */
- solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ valid = solve_del2_u(E, E->u1, F, imp*v_res, lev);
+ if(!valid && (E->parallel.me==0)) {
+ fputs("Warning: solver not converging! 0\n", stderr);
+ fputs("Warning: solver not converging! 0\n", E->fp);
+ }
strip_bcs_from_residual(E, E->u1, lev);
Modified: mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c
===================================================================
--- mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/lib/Topo_gravity.c 2007-11-08 23:28:46 UTC (rev 8250)
@@ -30,10 +30,6 @@
#include "element_definitions.h"
#include "global_defs.h"
-static void geoid_from_buoyancy();
-static void geoid_from_topography();
-
-
void get_STD_topo(E,tpg,tpgb,divg,vort,ii)
struct All_variables *E;
float **tpg,**tpgb;
@@ -417,64 +413,51 @@
/* ===================================================================
=================================================================== */
-void compute_geoid(E, harm_geoid, harm_geoid_from_bncy,
- harm_geoid_from_tpgt, harm_geoid_from_tpgb)
- struct All_variables *E;
- float *harm_geoid[2], *harm_geoid_from_bncy[2];
- float *harm_geoid_from_tpgt[2], *harm_geoid_from_tpgb[2];
+static void geoid_from_buoyancy(struct All_variables *E,
+ float *harm_geoid[2])
{
- int i, p;
+ /* Compute the geoid due to internal density distribution.
+ *
+ * geoid(ll,mm) = 4*pi*G*R*(r/R)^(ll+2)*dlayer*rho(ll,mm)/g/(2*ll+1)
+ *
+ * E->buoyancy needs to be converted to density (-therm_exp*ref_T/Ra/g)
+ * and dimensionalized (data.density). dlayer needs to be dimensionalized.
+ */
- geoid_from_buoyancy(E, harm_geoid_from_bncy);
- geoid_from_topography(E, harm_geoid_from_tpgt, harm_geoid_from_tpgb);
-
- if (E->parallel.me == (E->parallel.nprocz-1)) {
- for (i = 0; i < 2; i++)
- for (p = 0; p <= E->sphere.hindice; p++) {
- harm_geoid[i][p] = harm_geoid_from_bncy[i][p]
- + harm_geoid_from_tpgt[i][p]
- + harm_geoid_from_tpgb[i][p];
- }
- }
-
-}
-
-
-static void geoid_from_buoyancy(E, harm_geoid)
- struct All_variables *E;
- float *harm_geoid[2];
-{
int m,k,ll,mm,node,i,j,p,noz,snode;
- float *TT[NCS],radius,*geoid[2],dlayer,con1,con2,scaling2,scaling;
+ float *TT[NCS],radius,*geoid[2],dlayer,con1,grav,scaling2,scaling;
double buoy2rho;
void sphere_expansion();
void sum_across_depth_sph1();
/* scale for buoyancy */
- scaling2 = -E->data.therm_exp*E->data.ref_temperature*E->data.density/E->control.Atemp;
+ scaling2 = -E->data.therm_exp*E->data.ref_temperature*E->data.density
+ / E->control.Atemp;
/* scale for geoid */
- scaling = 1.0e6*4.0*M_PI*E->data.radius_km*E->data.radius_km*
- E->data.grav_const/E->data.grav_acc;
+ scaling = 4.0 * M_PI * 1.0e3 * E->data.radius_km * E->data.grav_const
+ / E->data.grav_acc;
/* density of one layer */
for(m=1;m<=E->sphere.caps_per_proc;m++)
TT[m] = (float *) malloc ((E->lmesh.nsf+1)*sizeof(float));
/* sin coeff */
- geoid[0] = (float*)malloc((E->sphere.hindice+1)*sizeof(float));
+ geoid[0] = (float*)malloc(E->sphere.hindice*sizeof(float));
/* cos coeff */
- geoid[1] = (float*)malloc((E->sphere.hindice+1)*sizeof(float));
+ geoid[1] = (float*)malloc(E->sphere.hindice*sizeof(float));
/* reset arrays */
- for (p = 0; p <= E->sphere.hindice; p++) {
+ for (p = 0; p < E->sphere.hindice; p++) {
harm_geoid[0][p] = 0;
harm_geoid[1][p] = 0;
}
- /* loop over each layer */
+ /* loop over each layer, notice the range is [1,noz) */
for(k=1;k<E->lmesh.noz;k++) {
- buoy2rho = scaling2 / E->refstate.gravity[k];
+ /* correction for variable gravity */
+ grav = 0.5 * (E->refstate.gravity[k] + E->refstate.gravity[k+1]);
+ buoy2rho = scaling2 / grav;
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++) {
@@ -488,21 +471,21 @@
sphere_expansion(E,TT,geoid[0],geoid[1]);
/* thickness of the layer */
- dlayer = (E->sx[1][3][k+1]-E->sx[1][3][k]);
+ dlayer = (E->sx[1][3][k+1]-E->sx[1][3][k])*E->data.radius_km*1e3;
- /* radius of the layer */
+ /* mean radius of the layer */
radius = (E->sx[1][3][k+1]+E->sx[1][3][k])*0.5;
- /* geoid contribution of density at this layer */
- for (ll=1;ll<=E->output.llmax;ll++)
+ /* geoid contribution of density at this layer, ignore degree-0 term */
+ for (ll=1;ll<=E->output.llmax;ll++) {
+ con1 = scaling * dlayer * pow(radius,((double)(ll+2)))
+ / (2.0*ll+1.0);
for (mm=0;mm<=ll;mm++) {
- con1 = scaling/(2.0*ll+1.0)*dlayer;
- con2 = pow(radius,((double)(ll+2)));
-
p = E->sphere.hindex[ll][mm];
- harm_geoid[0][p] += con1*con2*geoid[0][p];
- harm_geoid[1][p] += con1*con2*geoid[1][p];
+ harm_geoid[0][p] += con1*geoid[0][p];
+ harm_geoid[1][p] += con1*geoid[1][p];
}
+ }
//if(E->parallel.me==0) fprintf(stderr,"layer %d %.5e %g %g %g\n",k,radius,dlayer,con1,con2);
}
@@ -520,33 +503,53 @@
-static void geoid_from_topography(E, geoid_tpgt, geoid_tpgb)
- struct All_variables *E;
- float *geoid_tpgt[2], *geoid_tpgb[2];
+static void geoid_from_topography(struct All_variables *E,
+ float *geoid_tpgt[2],
+ float *geoid_tpgb[2])
{
- float con1,con2,scaling,den_contrast1,den_contrast2,stress_scaling,topo_scaling1,topo_scaling2;
+ /* Compute the geoid due to surface and CMB dynamic topography.
+ *
+ * geoid(ll,mm) = 4*pi*G*R*delta_rho*topo(ll,mm)/g/(2*ll+1)
+ *
+ * E->slice.tpg is essentailly non-dimensional stress(rr) and need
+ * to be dimensionalized by stress_scaling/(delta_rho*g).
+ *
+ * In theory, the degree-0 and 1 coefficients of topography must be 0.
+ * The geoid coefficents for these degrees are ingnored as a result.
+ */
+
+ float con1,con2,scaling,den_contrast1,den_contrast2,stress_scaling,topo_scaling1,topo_scaling2,grav1,grav2;
int i,j,k,ll,mm,s;
float *tpgt[2], *tpgb[2];
void sum_across_depth_sph1();
void sphere_expansion();
for (i=0; i<2; i++) {
- tpgt[i] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
- tpgb[i] = (float *)malloc((E->sphere.hindice+2)*sizeof(float));
+ tpgt[i] = (float *)malloc(E->sphere.hindice*sizeof(float));
+ tpgb[i] = (float *)malloc(E->sphere.hindice*sizeof(float));
}
stress_scaling = E->data.ref_viscosity*E->data.therm_diff/
(E->data.radius_km*E->data.radius_km*1e6);
- /* density contrast across surface */
- 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)*E->refstate.rho[1];
+ /* density contrast across surface, need to dimensionalize reference density */
+ den_contrast1 = E->data.density*E->refstate.rho[E->lmesh.noz] - E->data.density_above;
+ /* density contrast across CMB, need to dimensionalize reference density */
+ den_contrast2 = E->data.density_below - E->data.density*E->refstate.rho[1];
+ /* gravity at surface */
+ grav1 = E->refstate.gravity[E->lmesh.noz] * E->data.grav_acc;
+ /* gravity at CMB */
+ grav2 = E->refstate.gravity[1] * E->data.grav_acc;
+
/* 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 * grav1);
+ topo_scaling2 = stress_scaling / (den_contrast2 * grav2);
+ /* scale for geoid */
+ scaling = 4.0 * M_PI * 1.0e3 * 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]);
@@ -557,10 +560,6 @@
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;
-
/* zero degree-0 and 1 term */
geoid_tpgt[0][E->sphere.hindex[0][0]]
= geoid_tpgt[0][E->sphere.hindex[1][0]]
@@ -632,6 +631,30 @@
+void compute_geoid(E, harm_geoid, harm_geoid_from_bncy,
+ harm_geoid_from_tpgt, harm_geoid_from_tpgb)
+
+ struct All_variables *E;
+ float *harm_geoid[2], *harm_geoid_from_bncy[2];
+ float *harm_geoid_from_tpgt[2], *harm_geoid_from_tpgb[2];
+{
+ int i, p;
+
+ geoid_from_buoyancy(E, harm_geoid_from_bncy);
+ geoid_from_topography(E, harm_geoid_from_tpgt, harm_geoid_from_tpgb);
+
+ if (E->parallel.me == (E->parallel.nprocz-1)) {
+ for (i = 0; i < 2; i++)
+ for (p = 0; p < E->sphere.hindice; p++) {
+ harm_geoid[i][p] = harm_geoid_from_bncy[i][p]
+ + harm_geoid_from_tpgt[i][p]
+ + harm_geoid_from_tpgb[i][p];
+ }
+ }
+
+}
+
+
/* ===================================================================
Consistent boundary flux method for stress ... Zhong,Gurnis,Hulbert
Modified: mc/3D/CitcomS/branches/v3.0/visual/Makefile.am
===================================================================
--- mc/3D/CitcomS/branches/v3.0/visual/Makefile.am 2007-11-08 23:28:38 UTC (rev 8249)
+++ mc/3D/CitcomS/branches/v3.0/visual/Makefile.am 2007-11-08 23:28:46 UTC (rev 8250)
@@ -131,8 +131,12 @@
chmod 644 "$(DESTDIR)$(visualdir)/$$f"; \
done
+
+bin_PROGRAMS = project_geoid
+project_geoid_SOURCES = project_geoid.c
+
if COND_HDF5
- bin_PROGRAMS = h5tocap h5tovelo
+ bin_PROGRAMS += h5tocap h5tovelo
h5tocap_SOURCES = h5tocap.c h5util.c h5util.h
h5tovelo_SOURCES = h5tovelo.c h5util.c h5util.h
h5tovelo_LDADD = $(LIBHDF5)
Copied: mc/3D/CitcomS/branches/v3.0/visual/project_geoid.c (from rev 8245, mc/3D/CitcomS/trunk/visual/project_geoid.c)
More information about the cig-commits
mailing list