[CIG-MC] Re: CitcomCU in 2D?
Shijie Zhong
Shijie.Zhong at Colorado.Edu
Tue Nov 6 14:49:46 PST 2007
I just sent Tobias an email about 2D citcom. Perhaps I should have just sent the email to cig-mc.
Basically, for 2d cartesian, if you do not want to have the parallel computing features, it is better to use other versions of citcom including Louis' original copy if you have one. Louis' original citcom works for both 2D/3D cartesian but for 2D, one does not need to have the third dimension (i.e., real 2D).
If there are enough interests in only 2D cartesian, which is still very useful, I probably can send in a version of my 2D cartesian, which is similar to Louis' original one, except for the full multigrid solvers (probably 2-3 times faster than the original). It should also have the particle stuff. Wei Leng, one of my students, has also recently implemented fully compressible convection in the code. Perhaps after AGU, we can work on it.
Shijie
Shijie Zhong
Department of Physics
University of Colorado at Boulder
Boulder, CO 80309
Tel: 303-735-5095; Fax: 303-492-7935
Web: http://anquetil.colorado.edu/szhong
---- Original message ----
>Date: Tue, 6 Nov 2007 14:29:08 -0800
>From: Magali Billen <billen at geology.ucdavis.edu>
>Subject: Re: [CIG-MC] Re: CitcomCU in 2D?
>To: Eh Tan <tan2 at geodynamics.org>
>Cc: Tobias Hoeink <Tobias.Hoeink at rice.edu>,cig-mc at geodynamics.org
>
> Hello Eh and Tobias,
> I've been using CitcomCU in 2D using the conjugate
> gradient solver.
> However, I found that there was a line missing from
> the function solve_del2_u
> when choosing the conjugate gradient solver. In
> the version I have (1.0.2),
> the parameter called valid was not set if you choose
> the conjugate gradient solver.
> This parameter is returned by the call to this
> function and is needed for the program
> to run correctly.
> I added the line below marked by /* MIB, FEB 07 */.
> This should be line 416 in
> General_matrix_functions.c
> if(!(E->control.NMULTIGRID ||
> E->control.EMULTIGRID))
> {
> valid = (residual < acc) ? 0 : 1; /* MIB, FEB
> 07 */
> cycles = E->control.v_steps_low;
> time = CPU_time0();
> residual = conj_grad(E, D1, r, Au, acc, &cycles,
> high_lev);
> for(i = 0; i < neq; i++)
> {
> d0[i] += D1[i];
> D1[i] = 0.0;
> }
> - Magali
> On Nov 6, 2007, at 2:04 PM, Eh Tan wrote:
>
> Hi Tobias,
> You can run "2D" model by using only 2 nodes in y
> direction. A small
> (but non-zero) dimeny might be required to ensure
> the element has an
> aspect ratio close to 1. You also have to use the
> conjugate gradient
> solver, not the multigrid solver. The multigrid
> solver has special
> requirement for the number of nodes in each
> direction.
> A working input file is attached. I only confirmed
> that the solver is
> working and converging, but didn't check the
> output, so the model may
> not make sense. Customize the input file for your
> need, and let me know
> if you have any problem running it.
> Cheers,
> --
> Eh Tan
> Staff Scientist
> Computational Infrastructure for Geodynamics
> 2750 E. Washington Blvd. Suite 210
> Pasadena, CA 91107
> (626) 395-1693
> http://www.geodynamics.org
> Tobias Hoeink wrote:
>
> Hello from Houston,
> Is there a way to run CitcomCU (cartesian
> geometry) in 2 dimensions? I
> have looked into the source a little and it
> looks like there is a
> 'cart2d' option. However, it appears disabled.
> Would it work to use 'cart3d' and to play with
> the geometry
> parameters, e.g. by setting dimeny=0.0 or the
> like?
> I would appreciate your help.
> Thanks,
> Tobias
> Tobias Höink
> Keith-Wiess Postdoctoral Fellow
> Department of Earth Science, MS-126
> Rice University
> 6100 Main St, Houston, Texas 77005
> Mailing Address:
> P.O. Box 1892, Houston, Texas 77251-1892
> Tel. (713) 348 4497
> Fax. (713) 348 5214
> Tobias.Hoeink at rice.edu
> <mailto:Tobias.Hoeink at rice.edu>
> http://terra.rice.edu/department/staff/tobias.hoeink/
>
> # Start up file for citcom using getpar.
> # 1. Input and Output Files Information
> datafile="CASE1/caseA" # directory and file head
> for output files.
> use_scratch="local" # output files under current
> directory.
> # use_scratch="szhong" # output files onto local
> disk to each compute node.
> oldfile="CASE1/caseA" # directory and file head
> for restart files.
> restart=0 # check Convection.c for different
> options
> restart_timesteps=20000 # timestep to restart
> stokes_flow_only=0 # 1: only solve velocity once;
> 0: time-dependent problem.
> maxstep=20000 # max time steps
> storage_spacing=10 # write data every ...
> # 2. Geometry, Ra numbers, Internal heating,
> Thermochemical/Purely thermal convection
> Solver=cgrad node_assemble=1 # conjugate
> gradient: unsure if it still works
> rayleigh=10.97394e5 # Rayleigh number
> rayleigh_comp=1e6 # Compositional Rayleigh
> number. relevant for composition=1
> composition=0 # 0: purely thermal convection; 1:
> thermochem
> Q0=0 # Dimensionless internal heating rate
> Q0_enriched=0 # Q0 for C=1 layer, when relevant
> markers_per_ele=15 # number of particles per
> element
> comp_depth=0.605 # initial depth of compositional
> interface
> visc_heating=0 # 1: visc heating on; or 0: visc
> heating off.
> adi_heating=0 # 1: adiabatic heating on; or 0:
> adiabatic heating off.
>
> # 3. Grid And Multiprocessor Information
> nprocx=1 # number of processors in x or theta
> dir
> nprocz=1 # number of processors in z or r dir
> nprocy=1 # number of processors in y or fi dir
> nodex=15 nodez=15 nodey=2 # only relevant
> with conj-grad
> mgunitx=14 # multigrid base level in x or t
> mgunitz=14 # multigrid base level in z or r
> mgunity=1 # multigrid base level in y or f
> levels=1 # and how many times it gets doubled
> # 4. Coordinate Information
> Geometry=cart3d # cart3d or Rsphere
> #CART: irrelevant if Geometry=Rsphere.
> dimenx=1.0 # box size in x direction
> dimenz=1.0 # box size in z direction
> dimeny=0.01 # box size in y direction
> z_grid_layers=4 # minus 1 is number of layers
> with uniform grid in z.
> zz=0.0,0.1,0.9,1.0 # starting and ending z
> coodinates
> nz=1,3,13,15 # starting and ending node in
> z direction
> x_grid_layers=2
> xx=0,1
> nx=1,15
> y_grid_layers=2
> yy=0,1
> ny=1,2
> z_lmantle=0.76655052 # 2870 km-670 km
> z_410=0.857143 # 2870 km-410 km
> z_lith=0.9651568 # 2870-100 km
> # Regional SPHERICAL: irrelevant if
> Geometry=cart3d
> radius_inner=0.55 radius_outer=1.0 # inner and
> outer radius
> theta_north=73.5 theta_south=106.5 # colatitude
> in deg
> fi_west=0 fi_east=36.0 # longitude in deg
> r_grid_layers=4
> rr=0.55,0.59,0.96,1.0
> nr=1,6,44,49
> t_grid_layers=2
> tt=73.5,106.5
> nt=1,49
> f_grid_layers=2
> ff=0,36
> nf=1,49
> r_lmantle=0.89482 # (Ro-670 km)/Ro
> r_410=0.9356358
> r_lith=0.984301 # (Ro-100 km)/Ro
> # 5. Rheology
> rheol=0 # 0,1,2,... diff. option for diff.
> rheology. Check Visco....c
> TDEPV=off # on/off for temperature-dependent
> viscosity
> VISC_UPDATE=off # on/off for updating viscosity
> or not.
> update_every_steps=2 # update viscosity every n
> timesteps
> num_mat=4 # number of material group with possibly
> different rheology
> visc0=1.0e0,1.0e0,1.0e0,1.0e0 # pre-exponential
> constant
> viscE=6.9077553,6.9077553,6.9077553,6.9077553
> #activation energy
> viscT=273,273,273,273 #surface temperature:
> viscZ=5e-6,5e-6,5e-6,5e-6 #activation volume
> SDEPV=off # on/off for non-Newtonian
> sdepv_misfit=0.010 # accuracy for non-Newtonian
> iteration
> sdepv_expt=1,1,1,1 # for each mat group, stress
> exponent
> sdepv_trns=1.e0,1.e0,1.e0,1.e0 # transition
> stress
> VMIN=on visc_min=5.0e-2 # viscosity lower cutoff
> VMAX=on visc_max=2.0e04 # viscosity upper cutoff
> visc_smooth_cycles=1 # how viscosity is smoothed
> in multigrid
> Viscosity=system # always
> # 6. DIMENSIONAL INFORMATION and Depth-dependence
> layerd=2870000.0 #meter
> radius=6370000.0
> ReferenceT=3800.0
> refvisc=1.0e20
> density=3300.0
> thermdiff=1.0e-6
> gravacc=9.8
> thermexp=5e-5
> cp=1250
> wdensity=0.0
> visc_factor=1.0
> thermexp_factor=1.0
> thermdiff_factor=1.00
> dissipation_number=2.601
> surf_temp=0.078947
> # 7. phase changes: to turn off any of the phase
> changes, let Ra_XXX=0
> Ra_410=0.0 #kg/m^3
> Ra_670=0.0
> clapeyron410=3.0e6 #Pa K-1
> clapeyron670=-3.0e6
> width410=3.5e4 #meter
> width670=3.5e4
> # 8. BOUNDARY CONDITIONS and Initial perturbations
>
> topvbc=0 # velocity boundary conditions top and
> bottom
> topvbxval=0.0
> topvbyval=0.0
> botvbc=0
> botvbxval=0.0
> botvbyval=0.0
> #
> toptbc=1 bottbc=1 # temperature bc's top and
> bottom
> toptbcval=0.0 bottbcval=1.0 #
> periodicx=off #
> periodicy=off #
> flowthroughx=off #
> flowthroughy=off #
> num_perturbations=1 # N, Number of
> perturbations
> perturbmag=0.001 # A list of N magnitudes
> perturbk=1.0 # A list of N wavenumbers (/PI)
> perturbl=6.0 # A list of N wavenumbers (/PI)
> perturbm=0.0 # A list of N wavenumbers (/PI)
> # 9. SOLVER RELATED MATTERS
> Problem=convection # always, almost
> aug_lagr=on
> aug_number=1.0e3
> precond=on
> orthogonal=off
> maxsub=1
> viterations=2 # Uzawa iteration loops.
> vhighstep=3 # Smoothing passes at highest level
> (finest grid).
> piterations=375 # Uzawa iteration loops.
> accuracy=1.0e-4 # Desired accuracy of Uzawa
> algorithm.
> tole_compressibility=1e-7
> # Tuning of energy equation
> adv_sub_iterations=2
> finetunedt=0.75
> ll_max=20
> nlong=180
> nlati=90
> # Data input and program debugging
> DESCRIBE=off #
> BEGINNER=off #
> VERBOSE=off #
> verbose=off #
> COMPRESS=off #
> see_convergence=1
> # vim:ts=8:sw=8
> _______________________________________________
> CIG-MC mailing list
> CIG-MC at geodynamics.org
> http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc
>
> ---------------------------------------------
> Assistant Professor
> U.C. Davis, Geology
> (530) 754-5696
> billen at geology.ucdavis.edu
> ----------------------------------------------
> "Life is not measured by the breaths you take, but
> by the moments that take your breath away."
>________________
>_______________________________________________
>CIG-MC mailing list
>CIG-MC at geodynamics.org
>http://geodynamics.org/cgi-bin/mailman/listinfo/cig-mc
More information about the CIG-MC
mailing list