Download Citcom Tutorial: Convection Modeling with Finite Element Code and more Study notes Geology in PDF only on Docsity!
SIO 224
Convection modeling with citcom
- Introduction
This tutorial is adapted from a CIDER tutorial given by Peter van Keken, Carolina Lithgow-Bertelloni, and Louise Kellogg. This tutorial will allow you to investigate some aspects of mantle convection modeling using the finite element code Citcom. We will first explore the installation and basic use of the code and associated applications before investigating the relationships between the Nusselt and Rayleigh numbers for isoviscous and (for the brave hearted) temperature-dependent viscosity.
- Preliminaries
Citcom
Citcom is the name associated with a number of related finite element codes that found their origin at Caltech in the mid-90s. Louis Moresi is the principle author of the version we will be using. You are welcome to keep a version of this code but -as Louis requests in the source code- please be respectful of the time that went into creating this code. The code solves the incompressible Boussinesq equations for mantle convection which, from the notes are
∇ · v = 0 (1)
∇P ′^ = ∇ · T′^ − RaT ˆg (2)
∂T
∂t
where Q is a non-dimensional number representing internal heating.
Compiling and running the code
To run Citcom and analyze the models you’ll need to use the Unix shell. Locate the CITCOM directory and change into it. You should have a list of directories (such as Citcom.build.file, Perl scripts, tutorial). Change directories to tutorial. Here you should see a few example directories, including 00 test, a doc directory (which contains this document) and a number of executables (citcom, extract Nu, extract vrms). Your Mac should be able to run the citcom executable as is. Check this by changing the directory to 00 test and running
../citcom input-test
If you get messages like "2: divergence" you should be in business. This program should terminate in less than a minute with a message about a ppm file being written and some statements about cpu timing. If there are error messages you may need to recompile citcom. To recompile, change to the Shell- scripts directory in the CITCOM folder and type "build-citcom". When the build is complete, type "cp cll* tutorial/citcom" then retry the test.
- Basic model setup
In all cases we will assume a 2D Cartesian geometry. The initial model has aspect ratio 1 and a fluid of uniform viscosity. The box is heated from below (fixed temperature 1) and cooled from above (fixed temperature 0). The side boundaries have reflecting (symmetry) boundary conditions. The boundaries of the bottom are stress-free (i.e., not rigid). The equations are solved by a multigrid method, which requires the discretization into a fairly regular grid. In these examples we’ll limit ourselves to a uniform discretization, with the initial model at 33x33 nodal points. The initial condition for temperature is a conductive solution with a small (10%) harmonic perturba- tion that will initiate convection when the Rayleigh number is large enough. You can also specify the name of a file that contains the initial condition (which we will use to restart convection runs). In the first few examples we will try to get a steady-state, which in this code is obtained by letting time integration continue long enough. This works fine at moderate Ra, but becomes increasingly harder when the convection ’wants’ to be time-dependent.
- First full convection run
Find the example director 0 Ra1d4. Check out the input file ’input-2d Ra1e4’. You can do this with any editor through MacOS, or with vi or emacs from the shell. The top 10 lines set some variables which are defined as:
rayleigh: Thermal Rayleigh number maxstep: Maximum number of time steps finetuned: The multiplication factor for the self-selecting timestep. datafile: Common root of all output filenames storage spacing: Store output files every... write spacing: Store output files every... TDEPV: Indication whether temperature dependent viscosity is used rheology T type: Choice of viscosity function viscT1, viscN0: Parameters in viscosity function dimenx: Non-dimensional grid size of the mesh in x-direction levels: Number of multigrid levels. perturbk: Wave number of the initial perturbation previous temperature file: Name of file with initial condition.
The last variable is commented out (with a hash mark at the start of the line. Remove this hash mark when you want to do a restart from an existing file. The rest of the file is better left unchanged.... Run the citcom code with this input file
../citcom input-2d Ra1e
If you don’t like the output rolling on the screen or would like to check things in this shell while the code runs you can redirect the output and run it in the background with
nohup ../citcom input-2d Ra1e4 &
Use ’jobs’ to check whether it is still running, ’fg’ to bring it to the foreground, ’ctrl-C’ to kill it. The code will produce files ending in .horiz ave, .node data, and .ppm. The first contains, well, horizontal averages; the second contains the temperature in the nodes, and last one is a graphics file that can be converted to gif by a program called GraphicConvertor (if you don’t have this on your Mac
cp ra1e4-2d.00500.node data init.node data ../citcom input-etaT Ra1e
The viscosity in this case is prescribed by
η(T ) = exp (−bT ) (4)
with the parameter b represented by viscT1 in the input file. The input file in this directory is set up for the Blankenbach benchmark 2a (Ra 104 , b = ln 1000 , which has ’best’ values of Nu=10.066 and Vrms=480.43). Why are Nu and Vrms in this case higher than that for the isoviscous case with the same Ra? Can you get accurate results at 33x33? If you have time, try a few different Rayleigh numbers to see how the Nu/Ra relationship is affected Note that the code is quite versatile so you can explore any number of things including the effects of internal heating, etc. Feel free to play.
