/
README
323 lines (245 loc) · 14.9 KB
/
README
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
A Modified 3D Cartesian/Regional Spherical Citcom
Shijie Zhong
Dept of Physics
University of Colorado
Boulder, Colorado 80309
szhong@anquetil.colorado.edu
1. Background information.
This parallel 3D Cartesian/Regional spherical Citcom code was modified
from the original Citcom code that was originally developed by L. Moresi
[Moresi and Gurnis, 1996]. The current version has the following
functionalities: modeling thermal convection and thermochemical convection
in either 3D Cartesian or 3D Regional spherical geometry. The current
version also considers the Boussinesq and extended-Boussinesq approximations
(i.e., with phase changes, adiabatic heating, viscous heating, and latent
heating [Christensen and Yuen [1985]). Compared with the original muiltigrid
Citcom, the major new additions are: 1) parallel computing with MPI, 2) both
Cartesian and Regional Spherical geometry, 3) thermochemical convection
with particles, 4) full-multigrid solver with consistent projection,
and 5) extended Boussinesq approximations. For more information on the code
and benchmarks, see Zhong [2005a,b].
It is worthwhile to point out that there is an earlier version of Citcom that
also computes isochemical thermal convection for regional spherical geometry
[Billen et al., 2003] but with different implementation and formulation.
2. How to compile and run the code.
2.1. Compiling the code
In the src/ directory of this distribution, you will find 31 C source
files (*.c), 7 header files (*.h), and a makefile (Makefile).
Before compiling the code, please make sure that your system is installed
with MPI. You may need to edit the Makefile to specify the location of mpicc
and set various compiler flags for your specific machine. Then, to compile
the code and produce an executable citcom.mpi, type
$ make citcom.mpi
2.2. Running the code
In the examples/ directory of this distribution, you will find several
input files. One of them (input1) is a sample input file. There are also
some input files (in the examples/Busse1993 directory) for benchmarks
against Busse et al. [1993].
You may first create directories for storing your output files, before
running a calculation. Suppose that you use the sample input file input1,
you need to create directory CASE1 (see input1) in the directory from
which your code is going to run (often on a disk associated with the
head node of your cluster). You also need to create directory CASE1
on disks that are local to each compute node for this calculcation
(e.g., compute_node01, compute_node02, ..., compute_nodeXX). Most of
the output files are to be written to disks associated with compute
nodes.
Then to run a calculation using citcom.mpi with input file input1, the
number of CPU NP, and machine file mc, type
$ mpirun -np NP -nolocal -machinefile mc citcom.mpi input1
where mc is a file listing processors to be used for the calculation. The
above command running the code may be system-dependent, in particular,
the options -nolocal and -machinefile. You may need to discuss it with
your system administrator.
3. Getting to know the input file
The sample input file input1 provides some explanations for the various
input parameters. You may have to start with some simple calculations that
may only need a few essential parameters and then gradually learn other
functionalities of the code and their associated parameters.
The input file is divided into 10 sections which are briefly discussed
here.
Sect. 1. Input and Output Files Information.
The parameters in this section specify where the output and input files
are. For example, the first parameter line
datafile="CASE1/caseA"
specifies the directory and file head for all the output files.
Some output files, including log files, are stored in the directory
from which you run the calculation (often on the head node), and in
this example case, in subdirectory CASE1 as caseA.log0 and so on.
Most other output files (velocity or temperature) are stored in
disks on compute nodes in directory CASE1 with file head caseA
(e.g., caseA.ave.0.10000 and so on). Refer to Output.c to see
exactly what data is printed in the output files.
The second parameter "use_scratch" specifies whether output files
are stored to compute nodes (recommended mode).
The next three parameters are about restarting a calculation for which
previous output, e.g., temperature, may be needed.
The next two important parameters "maxstep" and "storage_spacing" specify
the number of timesteps for the model run and output frequency. Here
output frequency for temperature and averaged properties (e.g., heat
flux) is often different. For details, check Output.c. You can modify
output in anyway you want.
Sect. 2. Geometry, Ra numbers, Internal heating, Thermochemical/Purely thermal convection.
The first two parameters specify the multigrid solver, which may never
need to change. The next two parameters are Rayleigh number and compositional
Rayleigh number (their ratio is the buoyancy number). The compositional
Rayleigh number is only relevant if the next parameter "composition=1"
which turns on thermochemical convection.
The next two parameters "Q0" and "Q0_enriched" are nondimensional internal
heat generation rate for normal mantle and dense component of the mantle,
the latter is only relevant if "composition=1".
Parameters "markers_per_ele" and "comp_depth" specify the average number of
particles per element and initial depth of the density interface, again
both are only relevant if "composition=1".
The last two parameters "visc_heating" and "adi_heating" are switches for
viscous heating and adiabatic heating that are implemented similar to
Christensen and Yuen [1985] for extended Boussinesq approximation.
Sect. 3. Grid And Multiprocessor Information
The first three parameters specify the number of processors used in
x (theta), z (r) and y (fi) directions. Of course, we assume that
the number of elements in each direction is fixed. The product of
three numbers is the total number of processors used for the calculation.
The next three parameters "nodex, nodey, and nodez" are numbers of
nodes in each direction but ONLY relevant for conjugate gradient
solver.
The next four parameters specify the total number of elements for
the calculation with multigrid solver. The first three give the
base level grid in each direction, and the fourth parameter "levels"
indicates how many times it gets doubled. For example, in the sample
input file input1, mgunitx=6,mgunitz=6,mgunity=6,levels=4, and the
total number of elements in each direction are 48, i.e., grid size
of 48x48x48.
Note that number of processors in each direction is not entirely
independent to the base level grid. More specifically, mgunitx must
be divisible by nprocx, and the same rule applies to y and z directions.
Sect. 4. Coordinate Information
The first parameter "Geometry" specifies which geometry you want for
the calculation, either 3D Cartesian or 3D regional spherical.
If 3D Cartesian is chosen, then the parameters for regional spherical
are irrelevant, and vice verse.
Here I will explain the setup for Cartesian cases, and regional spherical
parameters are similar.
For Cartesian setup, parameters "dimenx", "dimeny", and "dimenz" give
the box size in each direction. The next three parameters "z_grid_layers",
"zz" and "nz" are for coordinate information for z direction.
"z_grid_layers" minus 1 is the number of layer in which the grid spacing
is uniform. "zz" and "nz" give the starting and ending z coordinates and
nodal index. For the sample input, z_grid_layers=4, zz=0.0,0.1,0.9,1.0,
nz=1,7,43,49, are for three grid layers that are bounded by coordinates
and nodal indices given in zz and nz, in z direction. The next six
parameters are for x and y directions, respectively. In the sample file,
they give uniform spacing in those two directions.
The next three parameters "z_lmantle", "z_410" and "z_lith" are the
nondimensional depth for 670-km, 410-km, and lithosphere. This information
may be used later in rheology definition.
Sect. 5. Rheology
The first parameter "rheol" specifies the type of rheological equation
that one may use. You need to check Viscosity_structures.c to see what are
available and it is very easy to implement your own. There are three
rheological equations (options) available with this routine: rheol=0,1,or 2.
* For rheol=0, eta = N0*exp[E(1-T)], where N0 is the pre-exponential
factor, and all are nondimensional.
* For rheol=1, eta = N0*exp[E/(T+T_offset)], where T_offset is the
nondimensional surface temperature caused by nondimensionalization of
absolute temperature.
* For rheol=2, eta = N0*exp{[E+(1-z)*V]/(T+T_offset)}, where (1-z) is the
nondimensional depth and V is nondimensional activation volume.
* For rheol=3, eta = N0*exp[E(T_offset-T)].
* For rheol=4, eta = N0*exp[E(T_offset-T)+(1-z)*V].
The second parameter "TDEPV" is a switch for temperature-dependent viscosity.
The parameter "VISC_UPDATE" is a switch for whether viscosity is to be updated.
It's clear that if "TDEPV" is on, then "VISC_UPDATE" should also be on.
The next parameter "update_every_steps" indicates how often viscosity is
to be updated (admittedly, this is a short-cut, but our experience seems
to suggest that viscosity and stiffness matrix do not have to be updated
every time step).
The parameter "num_mat" specifies the number of material group that may
have different rheology (e.g., upper mantle, lower mantle, ...). In the
sample input file, num_mat=4 i.e., there can be four different material
groups. To see how they are defined, check routine construct_mat_group in
Construct_arrays.c.
The next two parameters "visc0" and "viscE" are pre-exponential constants and
activation energy for each material group. In the sample input file, there
are four numbers for each parameter, and each number is for a material group.
The next two parameters "viscT" and "viscZ" can be defined for your own
purposes. And check Viscosity_structures.c to see how you may use them.
Parameter "SDEPV" is a switch for non-Newtonian rheology, and "sdepv_misfit"
is the accuracy level for non-Newtonian iterations. Parameters "sdepv_expt"
and "sdepv_trns" are stress exponent and transition stress for each material
group. Non-Newtonian rheology is implemented in the code, but is rarely used.
So be careful with non-Newtonian rheology, if you decide to use it.
Parameters "VMIN" and "visc_min" specify whether or not a minimum viscosity
is imposed and what it is, if imposed. Likewise, "VMAX" and "visc_max" are
for the maximum viscosity.
Parameters "visc_smooth_cycles" and "Viscosity" do not need to vary.
Sect. 6. Dimensional information and depth-dependence.
The dimensional information is not really very useful except when 410-km
and 670-km phase changes are used (next section).
Parameters "visc_factor" specifies a factor of viscosity increase with depth
as a linear function. This is in addition to whatever layered viscosity
structure we impose in section 5. In the sample input, "visc_factor=1" implies
that this linear increase with depth does not exist.
Parameters "thermexp_factor" and "thermdiff_factor" are the factors of thermal
expansion decrease and thermal diffusivity increase with depths as linear
functions.
Parameters "dispation_number" and "surf_temp" are dissipation number and non-
dimensional surface temperature.
Sect. 7. Phase changes
Parameters "Ra_410" and "Ra_670" are the density drops at 410-km and 670-km
phase changes in SI unit. The next two parameters are the Clapeyron slopes.
Parameters "width410" and "width670" are the phase change widths (see
Christensen and Yuen [1985] for details).
You may need to check Phase_change.c to see how these dimensional numbers
get nondimensionalized.
Sect. 8. Boundary conditions and initial perturbation
The default boundary conditions are free-slip on all sides. Periodic boundary
conditions are also implemented, however, you may want to talk to me before
using it.
Parameter "num_perturbations" specifies the number of perturbations of
different wavelengths. Parameters "perturbmag" and "perturbk" are the
magnitude and wavenumber of the perturbation. If you want, you may specify
perturbation for spherical harmonic order and degree for regional spherical
calculations with parameters "perturbl" and "perturbm". Initial fields are
specified in Convection.c (convection_initial_temperature). You should take
a look at this section of code, as they are often modified.
Sect. 9. Solver related matters.
Parameters in this section are rarely changed. Two parameters you may need to
be aware of are "accuracy" and "tole_compressibility", and they determine the
accuracy of the solutions.
4. Utilities
In the util/ direction of this distribution, there is a shell script to
convert the output file to paralllel VTK files, which can be visualized
by several visualization packages, such as Paraview and LLNL VisIt. To
invoke the script, type
$ citcomcu_write_vtk
to see the usage and arguments.
The script has the following constraints:
* The script can only convert files stored on local file system.
* The script does not convert spherical coordinates to Cartesian
coordinates.
* The converted VTK files have their (x,y,z) coordinate axes mapped to the
(z,x,y) axes, respectively, in CitcomCU data.
5. Acknowledgement and citation issues
The developers of this code have spent considerable amount of time in developing
and testing the code. We would appreciate if you can reference the following
two papers for publications that result from using this code.
Zhong, S., Constraints on thermochemical convection of the mantle from plume
heat flux, plume excess temperature and upper mantle temperature,
J. Geophys. Res., 111, B04409, doi:10.1029/2005JB003972, 2006.
Moresi L.N. and M. Gurnis, Constraints on lateral strength of slabs
from 3-D dynamic flow models, Earth Planet. Sci. Lett., 138, 15-28, 1996.
6. References
Billen, M.I., Gurnis, M., and Simons, M., Multiscale dynamic models of the
Tonga-Kermadec subduction zone, Geophys. J. Int., 153, 359-388, 2003.
Busse, F.H. et al., 3D convection at infinite Prandtl number in Cartesian
Geometry -- a benchmark comparison, Geophys. Astrophys. Fluid Dynamics,
75, 39-59, 1993.
Christensen, U.R., and Yuen, D.A., Layered convection induced by phase changes,
J. Geophys. Res., 90, 10,291-10,300, 1985.
Moresi L.N. and M. Gurnis, Constraints on lateral strength of slabs
from 3-D dynamic flow models, Earth Planet. Sci. Lett., 138, 15-28, 1996.
Zhong, S., Constraints on thermochemical convection of the mantle from plume
heat flux, plume excess temperature and upper mantle temperature,
J. Geophys. Res., 2005a, in review.
Zhong, S., Dynamics of thermal plumes in 3D isoviscous thermal convection,
Geophys. J. Int., 162, 289-300, 2005.