https://github.com/geodynamics/citcoms
Tip revision: cbfba74
citcoms-manual.tex
\documentclass[english]{book}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{geometry}
\geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
\setcounter{secnumdepth}{3}
\setcounter{tocdepth}{5}
\usepackage{color}
\usepackage{babel}
\usepackage{array}
\usepackage{longtable}
\usepackage{varioref}
\usepackage{float}
\usepackage{textcomp}
\usepackage{url}
\usepackage{graphicx}
\usepackage{times}
\usepackage{textpos}
\usepackage[unicode=true]{hyperref}

\makeatletter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
\newenvironment{lyxcode}
{\par\begin{list}{}{
\setlength{\rightmargin}{\leftmargin}
\setlength{\listparindent}{0pt}% needed for AMS classes
\raggedright
\setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt}
\normalfont\ttfamily}%
\item[]}
{\end{list}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\usepackage{hyperref}

\let\myUrl\url
\renewcommand{\url}[1]{(\myUrl{#1})}

\raggedbottom

\makeatother

\begin{document}

\definecolor{dark_grey}{gray}{0.3}
\definecolor{orange}{RGB}{255,177,20}

%LINE 1%
{
\renewcommand{\familydefault}{\sfdefault}

\pagenumbering{gobble}
\begin{center}
\resizebox{\textwidth}{!}{\textcolor{dark_grey}{\fontfamily{\sfdefault}\selectfont
COMPUTATIONAL INFRASTRUCTURE FOR GEODYNAMICS (CIG)
}}

%\hrule

\color{dark_grey}
\rule{\textwidth}{2pt}

\color{dark_grey}
{\Large }
\end{center}

\begin{center}
\resizebox{\textwidth}{!}{\colorbox
{orange}{\fontfamily{\rmdefault}\selectfont \textcolor{white} {
\hspace{0.05in}CitcomS\hspace{0.05in}
}}}
\end{center}

\begin{textblock*}{0in}(0in,-0.3in)
\begin{center}
\vspace{.2in}
\includegraphics[height=5.5in]
{graphics/cookbook8.png}
\end{center}
\end{textblock*}

\color{dark_grey}
\hfill{\Huge \fontfamily{\sfdefault}\selectfont User Manual \\
\raggedleft \huge \fontfamily{\sfdefault}\selectfont Version {3.2.0}\\}

\null
\vfill
\color{dark_grey}
\Large \hfill {\raggedleft \fontfamily{\sfdefault}\selectfont
Eh Tan\\Michael Gurnis\\Luis Armendariz\\Leif Strand\\Susan Kientz\\
}
{\raggedright\fontfamily{\sfdefault}\selectfont www.geodynamics.org}

%\hrule
\vspace*{-0.275in}
%LINE%
\begin{center}
\color{dark_grey}
\rule{\textwidth}{2pt}
\end{center}
}

\pagebreak
\pagenumbering{arabic}

%\noindent \begin{center}
%\begin{figure}[H]
%\noindent \centering{}\includegraphics[width=0.75\paperwidth]{graphics/citcoms_cover.pdf}
%\end{figure}
%\par\end{center}
%\noindent \begin{center}
%\thispagestyle{empty}
%\par\end{center}
%\title{CitcomS User Manual}
%\author{© California Institute of Technology\\Version 3.2.0}

\title{CitcomS User Manual}
\date{\noindent \today}

\maketitle
\tableofcontents{}

\listoffigures

\part{Preface}

\chapter*{Preface}

This document is organized into three parts. Part I consists of traditional
book front matter, including this preface. Part II begins with an
introduction to Pyre and the Pyre-compatible version of CitcomS and
their capabilities and proceeds to the details of implementation,
including several cookbooks'' of short tutorials. Part III provides
appendices and references.

The style of this publication is based on the Apple Publications Style
Guide \url{developer.apple.com/documentation/UserExperience/Conceptual/APStyleGuide/AppleStyleGuide2003.pdf},
as recommended by Python.org \url{www.python.org}. The documentation
was produced using \LaTeX{}, specifically the \texttt{pdflatex} tool from the
TeXLive distribution \url{www.tug.org/texlive}.

Errors and bug fixes in this manual should be directed to the CIG
Mantle Convection Mailing List \url{cig-mc@geodynamics.org}.

\section*{Who Will Use This Document}

This documentation is aimed at two categories of users: scientists
who prefer to use prepackaged and specialized analysis tools, and
experienced computational Earth scientists. Of the latter, there are
likely to be two classes of users: those who just run models, and
those who modify the source code. Users who modify the source are
likely to have familiarity with scripting, software installation,
and programming, but are not necessarily professional programmers.

The manual was written for the usage of CitcomS on a variety of different
platforms. CitcomS has run on shared memory computers (Sun, Hewlett-Packard,
SGI, and IBM), commercial distributed memory machines (Intel and Cray/SGI),
clusters (including machines on the NSF TeraGrid), Linux PCs and Mac
OS X desktops.

\section*{Citation}

Computational Infrastructure for Geodynamics (CIG) is making this
source code available to you in the hope that the software will enhance
your research in geophysics. The underlying C code for the finite
element package and the Python bindings for the framework were donated
to CIG in July of 2005. A number of individuals have contributed a
significant portion of their careers toward the development of CitcomS
and Pyre. It is essential that you recognize these individuals in
the normal scientific practice by citing the appropriate peer reviewed
papers and making appropriate acknowledgements.

The CitcomS development team asks that you cite both of the following:
\begin{itemize}
\item Zhong, S., M.T. Zuber, L.N. Moresi, and M. Gurnis (2000), The role
of temperature-dependent viscosity and surface plates in spherical
shell models of mantle convection. \emph{J. Geophys. Res.,} \emph{105,}
11,063-11,082
\item Tan, E., E. Choi, P. Thoutireddy, M. Gurnis, and M. Aivazis (2006),
GeoFramework: Coupling multiple models of mantle convection within
a computational framework,  \emph{Geochem., Geophys., Geosyst. 7,}
Q06001, doi:10.1029/2005GC001155.
\end{itemize}
following:
\begin{itemize}
\item McNamara, A.K., and S. Zhong (2004), Thermochemical structures within
a spherical mantle: Superplumes or Piles?, \emph{J. Geophys. Res.},
109, B07402, doi:10.1029/2003JB002847.
\end{itemize}
The developers also request that in your oral presentations and in
the authors of the code, and CIG \url{geodynamics.org}.

\section*{Support}

Pyre development is funded by the U.S. Dept. of Energy's Advanced
Simulation and Computing program \url{www.sandia.gov/NNSA/ASC} and
the National Science Foundation's \url{www.nsf.gov} Information Technology
Research (ITR) program (grant \#0205653). Continued support of CitcomS
is based upon work supported by the National Science Foundation under
Grant No. EAR-0406751. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the authors
and do not necessarily reflect the views of the National Science Foundation.

\section*{Conventions}

Throughout this documentation, any mention of username'' is meant
to indicate the user, meaning you should substitute your account name
in its place.

\part{Chapters}

\chapter{Introduction}

CitcomS is a finite element code designed to solve thermal convection
problems relevant to earth's mantle released under the GNU General
code runs on a variety of parallel processing computers, including
shared and distributed memory platforms. In an effort to increase
the functionality of CitcomS to include greater control during simulations
on large parallel systems, the software has been reengineered from
previous versions of CitcomS to work with a Python-based modeling
framework called Pyre. With Pyre, CitcomS can be dynamically coupled
with other CitcomS simulations or with other codes such as SNAC, which
solves crustal and lithospheric problems.

CitcomS is a finite element code written in C that solves for thermo-chemical
convection within a spherical shell. It can solve for problems within
either a full or a restricted spherical domain. Although the code
is capable of solving many different kinds of convection problems
using the flexibility of finite elements, there are aspects of CitcomS
which make it well-suited for solving problems in which the plate
tectonic history is incorporated. Variable viscosity, including temperature-,
pressure-, position-, composition-, and stress-dependent viscosity
are all possible.

The fundamental basis for the numerical solution of any time-dependent
convection problem is the sequential solution of an equation of motion
and an energy equation. Convection problems are initially valued with
boundary conditions, including all of the problems which are solved
with CitcomS. The normal sequence of steps for the solution of convection
problems starts with an initial temperature field. First, the momentum
equation is solved. The solution of this equation gives us the velocity
from which we then solve the advection-diffusion equation, giving
us a new temperature. CitcomS uses this interleaved strategy. It is
possible to run convection backward in time so as to guess an initial
condition for a normal forward running initial and boundary value
problem. However, users should be aware that even specialists in mantle
convection modeling are just now starting to explore methods in this
area and, as such, this is an emerging area of research.

This code uses an iterative solution scheme to solve for the velocity
and pressure and, as such, a converged solution cannot be guaranteed.
Nested inside the iterative scheme, the code uses either a conjugate
gradient solver or a full multi-grid solver to solve the discretized
matrix equations.

\section{History}

Citcom (for \textbf{C}alifornia \textbf{I}nstitute of \textbf{T}echnology
\textbf{Co}nvection in the \textbf{M}antle) was originally written
in the early 1990s by Louis Moresi. Although the code for three-dimensional
problems was incorporated from its inception, early versions of the
software only solved for time-dependent convection problems within
two-dimensional Cartesian domains. Moresi's original code turned out
to be incredibly modular and easily extensible. Consequently, the
fundamental finite element infrastructure which Louis wrote is still
in place and forms the basis for much of the code contained in the
present release.

In the mid-1990s Moresi wrote versions of the code that solved the
equations within three-dimensional Cartesian domains. Then Shijie
Zhong successfully parallelized Citcom using message passing routines
on a limited release Intel supercomputer. Zhong then created a spherical
version of the code which he named CitcomS. Lijie Han then created
a regional version of CitcomS as well as an alternate version of message
passing for an arbitrarily large number of processors. Clint Conrad
created the first Beowulf implementations of the code using the MPI
library, then Conrad and Eh Tan re-coded the message passing of the
fully spherical version so that problems run on arbitrarily large
numbers of processors could also be solved. A plethora of different
versions of Citcom exist both on computers at the California Institute
of Technology and around the world.

Consequently, by 2002, there were so many different versions of the
code that some rationalization was in order. The software was migrated
into a version control system and Eh Tan and Eun-seo Choi created
a version of CitcomS that generates either a fully spherical or regional
model. CitcomS was released to the community through the former GeoFramework
project as version 1.0 and 1.1.

By 2004, in order to increase the functionality of CitcomS, the developers
began to reengineer the code into an object-oriented environment specifically
so it could work with a Python-based modeling framework called Pyre.
This release of the software is essentially the product of those reengineering
efforts. Eh Tan was the principal developer of CitcomS, with considerable
help from Eun-seo Choi, Puru Thoutireddy, and Michael Aivazis.

CitcomS is one component of a larger collection of software encompassed
by the former GeoFramework project, a collaboration between the Center
for Advanced Computing Research (CACR) \url{www.cacr.caltech.edu}
and the Seismological Laboratory \url{www.gps.caltech.edu/seismo},
both at Caltech, and the Victorian Partnership for Advanced Computing
\url{www.vpac.org} in Australia. The GeoFramework project developed
a suite of tools to model multi-scale deformation for Earth science
problems. This effort was motivated by the need to understand interactions
between the long-term evolution of plate tectonics and shorter term
processes such as the evolution of faults during and between earthquakes.
During 2005 and 2006 much of the remaining software developed by GeoFramework
was released under a GPL license and made available from CIG \url{geodynamics.org}.

The second major release of CitcomS (2.0) incorporated the software
framework Pyre, free surface modeling methods, and stress boundary
conditions on the top and bottom surfaces. In the summer of 2005,
as part of the 2.0.1 release, CIG replaced the old build procedure
with the GNU Build System. A subsequent release, version 2.0.2, could
compile and run on 64-bit systems.

The release of CitcomS 2.1 incorporated new features and functionality,
the most important being the use of HDF5 (a parallel version of the
Hierarchical Data Format). The HDF5 format allows you to deal with
the massive data output created for production runs (see Chapter \vref{cha:Working-with-HDF5}).
This version accepted \texttt{.cfg} files on input, which are easier
to create and read. Other improvements included the incorporation
of geoid calculations that had been left out of earlier releases,
as well as new scripts to allow results to be visualized with MayaVi2
Mapping Tools (GMT) \url{gmt.soest.hawaii.edu} and OpenDX \url{www.opendx.org}.
Instructions were provided on using this version as a preinstalled
package on some of the NSF TeraGrid sites.

The release of CitcomS 2.2 incorporated the ability of tracing particles
in the flow. The tracer code was developed by Allen McNamara and Shijie
Zhong in 2004 and donated to CIG in early 2007. The tracer code has
a wide range of applications in the mantle convection. It can be used
in tracing the trajectory of passive particles, in delineating the
top boundary of subducted slabs to define the low viscosity wedges,
or in tracking the evolution of the chemical composition field.

The release of CitcomS 3.0 contains many new features, including two
implementations of compressible convection, one by Wei Leng and Shijie
Zhong and the other by Eh Tan; the ability to resume computation from
previous checkpoints; multi-component chemical convection; a fixed
non-Newtonian solver; an exchanger package for solver coupling; removing
the rigid body rotation component from the velocity by Shijie Zhong;
and an option to disable monitoring of maximum temperature. In addition,
a rheology option for pseudo-plasiticity, composition dependent viscosity
and heat generation, compressed ASCII output, and an easier way for

Additional changes that are not backward compatible include (1) the
viscosity field at element level is not smoothed (this might slow
down the convergence but will represent the viscosity field more accurately);
(2) the Lenardic filter on temperature is disabled by default; (3)
the rigid body rotation component is removed from the velocity by
default; (4) use of a better pseudo-random number generator to generate
the initial tracer; (5) the type of input parameter \texttt{coor}
is changed from a boolean to an integer; (6) setting \texttt{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 \texttt{tic\_method=-1}); and (7)
the input parameter \texttt{reset\_initial\_composition} becomes obsolete.
Among the seven backward-incompatible changes, the first four will
affect the results.

The release of CitcomS 3.1 added the ability to convert temperature
and composition to seismic velocities. The seismic velocities output
can be uploaded to the CIG Seismology Web Portal for SPECFEM3D simulation
to generate synthetic seismograms. Other enhancements include adding
a self-gravitational effect on geoid and adding the Consistent-Boundary-Flux
(CBF) method to compute dynamic topography, contributed by Shijie
Zhong. Also, GMT/NetCDF grd input was added for surface velocity boundary
conditions, initial temperature, material dependence and local Rayleigh
number in surface layers, contributed by Thorsten Becker. Mike Gurnis
top surface temperature boundary conditions. Finally, this version
allows multi-component chemical viscosity.

Some of the aforementioned features are backward-incompatible, including:
\begin{enumerate}
\item The global mesh is changed. The coordinates of all nodes are slightly
shifted laterally. The change is not visually discernible, but does
affect the solution compared to previous versions.
\item The convergence criterion of the Stokes solver is changed. A single
parameter accuracy'' is used to control the accuracy of the continuity
and momentum equations. The default value of accuracy'' was $10^{-6}$
and is now $10^{-4}$.
\item The input parameter mantle\_temp'' is moved from the {[}CitcomS.solver.param{]}
to the {[}CitcomS.solver.ic{]} section.
\end{enumerate}
This release of CitcomS 3.2 adds numerous changes since 3.1.1 which
was released in 2009. These include improved implementations of internal
stress, radial layer dependent viscosity, and velocity boundary conditions.
Improvements have also been made in solver convergence controls, VTK
output and several bug fixes. Finally, the Pyre configuration option
has been turned off by default, though it can be enabled during configuration.
Other enhancements include:
\begin{itemize}
\item If Python 2.7 is detected, Pyre is disabled due to conflicts.
\item Trial implementation of anisotropic viscosity.
\item Added support of multigrid solver in Exchanger.
\item Improved handling of grd-read velocity and traction boundary conditions.
\item Changed default parameters for z\_410, z\_lmantle, z\_lith to more
Earth-like parameters.
PDEPV psrw=on parameter.
\item Added Murnaghan's integrated linear EoS.
\item Added augmented Lagrangian component of stiffness matrix in cbf topo.
\item Performance optimizations for CUDA code. Contributed by Leif Strand.
\item Fixed several bugs in the code and manual, including:

\begin{itemize}
\item Fixed a bug in pseudo free surface, reported by Robert Moucha.
\item Correctly handle orphan tracers.
\item Added Pyre binding for write\_q\_files, z\_layer, r\_grid\_layers,
rr, nr, others.
\item Replaced fixed accuracy setting. Bug reported by Rob Moucha.
\item Fixed temperature perturbation bug for multiple processors.
\end{itemize}
\end{itemize}

Pyre is an object-oriented environment capable of specifying and launching
numerical simulations on multiple platforms, including Beowulf class
parallel computers and grid computing systems. Pyre allows the binding
of multiple components such as solid and fluid models used in Earth
science simulations, and different meshers. The Pyre framework enables
the elegant setup, modification and launching of massively parallel
three-dimensional solver applications.

The Pyre framework incorporates features aimed at enabling the scientific
non-expert to perform tasks easily without hindering the expert. Target
features for end users allow complete and intuitive simulation specification,
reasonable defaults, consistency checks of input, good diagnostics,
cycle, and good debugging support.

\section{Pyre and CitcomS}

Pyre provides a simulation framework that includes solver integration
The framework offers a way to add new solvers to CitcomS and to fine-tune
CitcomS simulations.

Developers have created Pyre classes for CitcomS to facilitate simulation
setup. However, they are not independent classes in a strict sense.
They still share the same underlying data structure and their functionality
is not divided clearly. CitcomS was not designed to be object-oriented
and to make it so would require significant investment of effort with
little return. However, the lack of object-oriented features does
not hinder the coupling of CitcomS with other solvers.

\section{Governing Equations\label{sec:Governing-Equations}}

With CitcomS, the mantle is treated as an anelastic, compressible,
viscous spherical shell under Truncated Anelastic Liquid Approximation.
With these assumptions, thermal convection is governed by the equations
for conservation of mass, momentum, and energy:

\begin{flushright}
$$\left(\rho u_{i}\right)_{,i}=0\label{eq:conservation of mass}$$

\par\end{flushright}

\begin{flushright}
$$-P_{,i}+\left(\eta(u_{i,j}+u_{j,i}-\frac{2}{3}u_{k,k}\delta_{ij})\right)_{,i}-\delta\rho g\delta_{ir}=0\label{eq:conservation of momentum}$$

\par\end{flushright}

\begin{flushright}
$$\rho c_{P}\left(T_{,t}+u_{i}T_{,i}\right)=\rho c_{P}\kappa T_{,ii}-\rho\alpha gu_{r}T+\Phi+\rho\left(Q_{L,t}+u_{i}Q_{L,i}\right)+\rho H\label{eq:conservation of energy}$$

\par\end{flushright}

\noindent where \emph{$\rho$} \textit{\emph{is the density,}} \emph{u}
is the velocity, \emph{P} is the dynamic pressure, $\eta$ is the
viscosity, $\delta_{ij}$ is the Kroneker delta tensor, $\delta\rho$
is the density anomaly, \emph{g} is the gravitational acceleration,
\emph{T} is the temperature, $c_{P}$ is the heat capacity, $\kappa$
is the thermal diffusivity, $\alpha$ is the thermal expansivity,
$\Phi$ is the viscous dissipation, $Q_{L}$ is the latent heat due
to phase transitions, and \emph{H} is the heat production rate. The
expression $X_{,y}$ represents the derivative of $X$ with respect
to $y$, where $i$ and $j$ are spatial indices, $r$ is the radial
direction, and $t$ is time. With phase transitions and temperature
and composition variations, the density anomalies are:

\begin{flushright}
$$\delta\rho=-\alpha\bar{\rho}(T-\bar{T_{a}})+\delta\rho_{ph}\Gamma+\delta\rho_{ch}C\label{eq:density anomalies}$$

\par\end{flushright}

\noindent where \emph{$\bar{\rho}$} \textit{\emph{is the radial profile
of density,}}  $\bar{T_{a}}$ is the radial profile of adiabatic temperature,
$\delta\rho_{ph}$ is the density jump across a phase change, $\delta\rho_{ch}$
is the density difference between the compositions\emph{, $\Gamma$}
\textit{\emph{is the phase function, and}} \emph{C} is the composition.
The phase function is defined as:

\begin{flushright}
$$\pi=\bar{\rho}g(1-r-d_{ph})-\gamma_{ph}(T-T_{ph})\label{eq:reduced pressure}$$

\par\end{flushright}

\begin{flushright}
$$\Gamma=\frac{1}{2}\left(1+\tanh\left(\frac{\pi}{\bar{\rho}gw_{ph}}\right)\right)\label{eq:phase function}$$

\par\end{flushright}

\noindent where $\pi$ is the reduced pressure, $d_{ph}$ and $T_{ph}$
are the ambient depth and temperature of a phase change, \emph{$\gamma_{ph}$}
is the Clapeyron slope of a phase change, and $w_{ph}$ is the width
of a phase transition.

These equations lead to the following normalization in which primed
quantities are nondimensional:

\begin{flushright}
$$\rho=\rho_{0}\rho^{'}\label{eq:rho dim}$$

\par\end{flushright}

\begin{flushright}
$$\alpha=\alpha_{0}\alpha^{'}\label{eq:alpha dim}$$

\par\end{flushright}

\begin{flushright}
$$g=g_{0}g^{'}\label{eq:g dim}$$

\par\end{flushright}

\begin{flushright}
$$\kappa=\kappa_{0}\kappa^{'}\label{eq:kappa dim}$$
$$\eta=\eta_{0}\eta^{'}\label{eq:eta dim}$$

\par\end{flushright}

\begin{flushright}
$$c_{P}=c_{P0}c_{P}^{'}\label{eq:cp dim}$$

\par\end{flushright}

\begin{flushright}
$$x_{i}=R_{0}x_{i}^{'}\label{eq:x dim}$$

\par\end{flushright}

$$u_{i}=\frac{\kappa_{0}}{R_{0}}u_{i}^{'}\label{eq:u dim}$$

\begin{flushright}
$$T_{0}=\Delta TT_{0}'\label{eq:T0 dim}$$

\par\end{flushright}

\begin{flushright}
$$T=\Delta T(T'+T_{0}')\label{eq:T dim}$$

\par\end{flushright}

\begin{flushright}
$$t=\frac{R_{0}^{2}}{\kappa_{0}}t^{'}\label{eq:t dim}$$

\par\end{flushright}

\begin{flushright}
$$H=\frac{\kappa_{0}}{R_{0}^{2}}c_{P0}\Delta TH^{'}\label{eq:H dim}$$

\par\end{flushright}

\begin{flushright}
$$P=\frac{\eta_{0}\kappa_{0}}{R_{0}^{2}}P^{'}\label{eq:p dim}$$
$$d_{ph}=R_{0}d_{ph}^{'}\label{eq:dph dim}$$

\par\end{flushright}

\begin{flushright}
$$\gamma_{ph}=\frac{\rho_{0}g_{0}R_{0}}{\Delta T}\gamma_{ph}^{'}\label{eq:gammaph dim}$$

\par\end{flushright}

\noindent where $\rho_{0}$ is the reference density, $\alpha_{0}$
is the reference thermal expansivity, $g_{0}$ is the reference gravity,
$\kappa_{0}$ is the reference thermal diffusivity, $\eta_{0}$ is
a reference viscosity, $c_{P0}$ is the reference heat capacity, $R_{0}$
is the radius of the Earth, $T_{0}$ is the temperature at the surface,
and $\Delta T$ is the temperature drop from the core-mantle boundary
(CMB) to the surface. Dropping the primes, the equations become:

$$u_{i,i}+\frac{1}{\bar{\rho}}\frac{d\bar{\rho}}{dr}u_{r}=0\label{eq:non-dimensional continuity eqn}$$

$$-P_{,i}+\left(\eta(u_{i,j}+u_{j,i}-\frac{2}{3}u_{k,k}\delta_{ij})\right)_{,i}+(Ra\bar{\rho}\alpha T-Ra_{b}\Gamma-Ra_{c}C)g\delta_{ir}=0\label{eq:non-dimensional momentum eqn}$$

\begin{flushright}
$$\begin{array}{c} \bar{\rho}c_{P}\left(T_{,t}+u_{i}T_{,i}\right)\left(1+2\Gamma\left(1-\Gamma\right)\frac{\gamma_{ph}^{2}}{d_{ph}}\frac{Ra_{b}}{Ra}Di\left(T+T_{0}\right)\right)=\bar{\rho}c_{P}\kappa T_{,ii}\\ -\bar{\rho}\alpha gu_{r}Di\left(T+T_{0}\right)\left(1+2\Gamma\left(1-\Gamma\right)\frac{\gamma_{ph}}{d_{ph}}\frac{Ra_{b}}{Ra}\right)+\frac{Di}{Ra}\Phi+\bar{\rho}H \end{array}\label{eq:non-dimensional energy eqn}$$

\par\end{flushright}

\noindent where \emph{Ra}, the thermal Rayleigh number, is defined
as:

\begin{flushright}
$$Ra=\frac{\rho_{0}g_{0}\alpha_{0}\Delta TR_{0}^{3}}{\eta_{0}\kappa_{0}}\label{eq:Ra, Rayleigh number}$$

\par\end{flushright}

\noindent This is not the usual definition of the Raleigh number that
is based on layer thickness; it is based on the radius of the Earth
$R_{0}.$ So for mantle convection problems where $R_{0}$ is slightly
more than twice the layer thickness, our $Ra$ is about a factor of
10 larger than by the usual definition. The phase-change Rayleigh
number, $Ra_{b}$, the chemical Rayleigh number, $Ra_{c}$, the internal
heating Rayleigh number, $Ra_{H}$, and the dissipation number \textit{Di}
, are defined as:

\begin{flushright}
$$Ra_{b}=Ra\frac{\delta\rho_{ph}}{\rho_{0}\alpha_{0}\Delta T}\label{eq:Rab, the phase-change Rayleigh number}$$

\par\end{flushright}

\begin{flushright}
$$Ra_{c}=Ra\frac{\delta\rho_{ch}}{\rho_{0}\alpha_{0}\Delta T}\label{eq:Rac, chemical Rayleigh number}$$

\par\end{flushright}

\begin{flushright}
$$Ra_{H}=RaH\frac{R_{0}^{3}-R_{CMB}^{3}}{3R_{0}^{3}}\label{eq:RaH, internal heating Rayleigh number}$$

\par\end{flushright}

$$Di=\frac{\alpha_{0}g_{0}R_{0}}{c_{P0}}\label{eq:Di, dissipation number}$$

\section{Numerical Methods}

The governing equations are solved with the finite element method
\cite{Hughes The Finite Element Method}. CitcomS employs an Uzawa
algorithm to solve the momentum equation coupled with the incompressibility
constraints \cite{Moresi/Gurnis Contraints,Ramage/Wathen Iterative solution}.
The energy equation is solved with a Steamline-Upwind Petrov-Galerkin
method \cite{Brooks A.N.}. Brick elements are used, such as eight
velocity nodes with trilinear shape functions and one constant pressure
node for each element. The use of brick elements in 3D (or rectangular
elements in 2D) is important for accurately determining the pressure,
which controls the dynamic topography, in incompressible Stokes flow
\cite{Hughes The Finite Element Method}. The discrete form of Equations
\ref{eq:non-dimensional continuity eqn} and \ref{eq:non-dimensional momentum eqn}
may be written in the following matrix form \cite{Zhong et al The role of temperature}:

\begin{flushright}
$$\mathrm{\left(\mathbf{B}^{\mathit{T}}+\mathbf{C}\right)\mathit{u=0}}\label{eq:discrete continuite eqn}$$

\par\end{flushright}

\begin{flushright}
$$\mathbf{A\mathit{u+}}\mathbf{B\mathit{p=f}}\label{eq:discrete momentum eqn}$$

\par\end{flushright}

\noindent where $\mathbf{A}$ is the stiffness'' matrix, \emph{u}
is a vector of unknown velocities, $\mathbf{B}$ is the discrete gradient
operator, $\mathbf{C}$ is the second term in Equation \ref{eq:non-dimensional continuity eqn},
\emph{p} is a vector of unknown pressures, and \emph{f} is a vector
composed of the body and boundary forces acting on the fluid. The
individual entries of $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$
and \emph{f} are obtained using a standard finite element formulation;
see \cite{Zhong et al The role of temperature} for the explicit entries.

In the incompressible case, $\mathbf{C}$ is zero. Equation \ref{eq:discrete momentum eqn}
can be transformed by premultiplying by $\mathrm{\mathbf{B}^{\mathit{T}}\mathit{\mathbf{A}^{\mathbf{\mathit{-1}}}}}$
and using Equation \ref{eq:discrete continuite eqn} to eliminate
the velocity unknowns:

\begin{flushright}
$$\mathbf{\mathbf{B}^{\mathit{T}}}\mathbf{A}^{\mathbf{\mathit{-1}}}\mathbf{B}\mathit{p=\mathbf{\mathbf{B}^{\mathit{T}}A^{\mathit{-1}}\mathit{f}}}\label{eq:cg}$$

\par\end{flushright}

\noindent This equation is solved using the Uzawa algorithm, an established
method for solving the minimization of a dual function \cite{Cahouet/Chabard Some fast 3D},
which simultaneously yields the velocity field. A conjugate gradient
scheme \cite{Ramage/Wathen Iterative solution,Atanga/Silvester Iterative methods}
is used for this iteration and forms the basis for the technique used
in CitcomS.

In the compressible case, there are two different strategies to solve
Equations \ref{eq:discrete continuite eqn} and \ref{eq:discrete momentum eqn}.
The first strategy is to add another layer of iterations when solving
Equation \ref{eq:cg}. The right-hand-side vector is updated by the
velocity solution of the previous iteration. This equation can be
solved using the same conjugate gradient scheme as the incompressible
case.

$$\mathbf{\mathbf{B}^{\mathit{T}}}\mathbf{A}^{\mathbf{\mathit{-1}}}\mathbf{B}\mathit{p^{(i)}=\mathbf{\mathbf{B}^{\mathit{T}}A^{\mathit{-1}}\mathit{f-\mathbf{C}u^{(i-1)}}}}\label{eq:iter-cg}$$

\noindent The second strategy is to transform Equation \ref{eq:discrete momentum eqn}
by premultiplying by $\mathrm{\left(\mathbf{B}^{\mathit{T}}+\mathbf{C}\right)\mathit{\mathbf{A}^{\mathbf{\mathit{-1}}}}}$
and using Equation \ref{eq:discrete continuite eqn} to eliminate
the velocity unknowns:
$$\left(\mathbf{B}^{\mathit{T}}+\mathbf{C}\right)\mathit{\mathbf{A}^{\mathbf{\mathit{-1}}}\mathbf{B}\mathit{p=\mathbf{\left(\mathbf{B}^{\mathit{T}}+\mathbf{C}\right)A^{\mathit{-1}}\mathit{f}}}}\label{eq:bicg}$$
This equation is solved using a bi-conjugate gradient stabilized scheme.

\section{Meshes and Geometry}

There are two forms of meshes and geometries for CitcomS. By default
CitcomS will produce a mesh within a regional geometry that is bound
by lines of constant latitude and longitude. There is an option to
generate a global mesh of a spherical shell.

For a regional mesh, CitcomS uses meshes that are regular, although
considerable flexibility exists for grid refinement in the regional
models. There is an option for mesh refinement in which the mesh is
refined as a function of latitude, longitude, or radius. Such refinement
is suitable for higher resolutions near boundary layers or within
the center of the map domain, but is incapable of increasing the resolution
near a curvilinear feature, as a plate boundary, unless that plate
boundary is orientated north-south or east-west.

In regional meshes, \emph{theta} (or $x$) is the colatitude measured
from the north pole, \emph{fi} (or $y$) is the east longitude, and
$z$ is the radius. \emph{theta} and \emph{fi} are in units of radians.
Figure \ref{fig:Global-Node-Numbering.} shows the organization of
the mesh in a regional problem. The numbering of the nodes is $z$-direction
first, then $x$-direction, then $y$-direction. This numbering convention
is used for the input and output data.

\noindent \begin{center}
\begin{figure}[H]
\begin{description}

\item [{\includegraphics[width=0.75\paperwidth]{graphics/c_fig3.jpg}}]~
\end{description}
\caption{\label{fig:Global-Node-Numbering.}Global Node Numbering. Left: Global
node numbering starts at the base of arrow A (theta\_min, fi\_min,
the tip. Upon reaching the tip, numbering continues from the base
of arrow B (\emph{nodez}+1) to its tip (2 \emph{nodez}), and so on
for all nodes on the plane fi = fi\_min. Right: After completing each
theta radius plane, the fi index is incremented and numbering commences
from (theta\_min, radius\_inner) as on the left.}
\end{figure}

\par\end{center}

For a global mesh, CitcomS is also capable of generating a mesh for
an entire spherical shell in which elements in map view are approximately
equal in area. In the full spherical mode, CitcomS has 12 caps numbered
0 to 11 (Figure \ref{fig:Topological-connectivity-of}).

\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics{graphics/cap_order.pdf}
\par\end{centering}

\caption{\label{fig:Topological-connectivity-of}Topological connectivity of
the 12 caps. N is the north pole and S is the south pole. The red
line marks the 0 degree meridian.}
\end{figure}

\par\end{center}

The caps are approximately square in map view so that the edges of
the square are oriented diagonally with respect to latitude and longitude.
The four corners of the domain are connected by great circles (Figure
\ref{fig:Orthographic-projection-of}). One would normally associate
at least one processor with one cap. However, CitcomS can further
decompose the domain such that additional processors are used to divide
caps uniformly along the two edges of the caps (Figure \ref{fig:Orthographic-projection-of})

\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[scale=0.75]{graphics/proc_mesh.eps}
\par\end{centering}

\caption{\label{fig:Orthographic-projection-of}Orthographic projection of
processors from a full CitcomS mesh in which there are 16 processors
in map view for each cap. The CitcomS cap is shown as distinct colors
while the processor domains within the caps are indicated by the intensity
of the color. This example was produced for a run with 2 processors
in radius such that the total number of processors was 12$\times$16$\times$2=384.}
\end{figure}

\par\end{center}

\chapter{Installation and Getting Help}

\section{Introduction}

CitcomS has been tested on Linux, Mac OS X, Sun Solaris and several
NSF TeraGrid platforms. You should have no trouble installing CitcomS
on most Unix-like systems.

Most users will install CitcomS from the released source package.
The following sections will lead you through the installation process.
the latest CitcomS source code directly from the CIG source code repository,
instead of using the prepared source package; see Section \ref{sec:Software-Repository}
later in this chapter.

\section{Getting Help}

For help, send e-mail to the CIG Mantle Convection Mailing List \url{cig-mc@geodynamics.org}.
You can subscribe to the Mailing List and view archived discussion
at the Geodynamics Mail Lists web page \url{geodynamics.org/cig/lists}.

\section{System Requirements}

Installation of CitcomS requires the following:
\begin{itemize}
\item A C compiler
\item An MPI library
\item If using Pyre: between Python 2.4 to 2.6, including header files (Python.h),
and the pythia package (see below).
\end{itemize}
You can install all the dependencies in one command on most Linux
machines. On Debian, Ubuntu or similar distributions, use this command
(as \texttt{root}):
\begin{lyxcode}
\$~apt-get~install~build-essential~python-dev~\textbackslash{}~\\ ~~~~~~~~~~openmpi-dev~openmpi-bin \end{lyxcode} On Red Hat, Fedora, CentOS, OpenSuSE or similar distributions, use this command (as \texttt{root}): \begin{lyxcode} \$~yum~install~make~automake~gcc~gcc-c++~kernel-devel~\textbackslash{}~\\
~~~~~~python~python-devel~openmpi~openmpi-libs~openmpi-devel
\end{lyxcode}
On Mac OS X, Python is installed by default. You will need to install
a C compiler and MPI library (see sections below).

MPI installations are typically configured for a particular compiler,
and provide a special wrapper command to invoke the right compiler.
Therefore, the choice of MPI implementation often determines which
C compiler to use.
\begin{quote}
\textbf{Note:} Users familiar with older versions of CitcomS may prefer
to install only the legacy CitcomS tools, \texttt{CitcomSFull} and
\texttt{CitcomSRegional}, and forgo use of Python and the Pyre framework.
By default in CitcomS 3.2.0 this is how CitcomS is configured. This
process requires only a C compiler and an MPI library. For more information,
see Section \ref{sec:Installing-without-Pyre} later in this chapter.
\end{quote}
If you are going to use solver coupling, you will also need a decent
C++ compiler (e.g., g++ 3.2 or newer) and the Exchanger package; see
Section \ref{sec:Exchanger-Configuration-(Optional)}. Optionally,
CitcomS can be configured to use a parallel HDF5 library; see Section
\vref{sec:HDF5-Configuration}.

\subsection{C Compiler}

On Unix or Linux systems, there is a high likelihood that a usable
C compiler is already installed. To check, type \texttt{cc} at the
shell prompt:
\begin{lyxcode}
\$~cc cc:~no~input~files \$
\end{lyxcode}
On Linux, if the \texttt{cc} command is not found, install GCC using
the package manager for your distribution.

The Mac OS X version of GCC is included in a software development
suite called Xcode. Xcode is available as a free download at the Apple
Developer Connection \url{developer.apple.com}.
\begin{quote}
\textcolor{red}{Warning:} \textcolor{black}{If you are using an Intel
compiler on an Itanium CPU, do not use the} \texttt{\textcolor{black}{-O3}}
\textcolor{black}{optimization flag as reports indicate that this
optimization level will generate incorrect codes. For any compiler,
you should always be careful about the correctness of the compiled
codes when using an} \texttt{\textcolor{black}{-O3}} \textcolor{black}{or
higher optimization level.}
\end{quote}

\subsection{MPI Library}

CitcomS requires a library which implements the MPI standard (either
version 1 or 2). Several free, open-source implementations of MPI
are available.

A popular choice is MPICH \url{www-unix.mcs.anl.gov/mpi/mpich}. Other
choices include LAM/MPI \url{www.lam-mpi.org} and Open MPI \url{www.open-mpi.org}.
Installing MPI from source involves walking through the standard GNU
build procedure (\texttt{configure \&\& make \&\& make install}) while
logged in as \texttt{root}.

Linux users may have a prebuilt MPI package available for their distribution.
On Mac OS X, the Fink package manager offers a prepackaged version
of LAM/MPI \url{www.lam-mpi.org}; so if you have Fink \url{fink.sourceforge.net}
installed, simply enter the following command from a Terminal window
to install LAM/MPI:
\begin{lyxcode}
\$~fink~install~lammpi~lammpi-dev \end{lyxcode} \subsubsection{MPI C Compiler Command} Once you have an MPI library installed, make sure its C complier command is on your PATH. Unfortunately, the name of this command varies from one MPI implementation to the next. The CitcomS configuration script searches for the following MPI C command names: \begin{lyxcode} mpicc~hcc~mpcc~mpcc\_r~mpxlc~cmpicc \end{lyxcode} \subsection{Python (Optional)} Your system may already have a suitable Python interpreter installed. To check, type the \texttt{python} command: \begin{lyxcode} \$~python~-V

Python~2.5.2
\end{lyxcode}
Mac OS X 10.3 and later ships with a suitable version of Python preinstalled.
If you're using an older version of Mac OS X, or for more information
at the Python web site.

On Linux, simply install the binary system package available for your
distribution. Be sure to select the Python development package (typically
called \texttt{python-dev}) in addition to the core Python package
-- even if you don't plan on doing any Python software development.
The development package contains the Python header files (Python.h),
which are necessary for building CitcomS.

If you are working on a cluster and \texttt{python} is too old, try
poking around a little. Sometimes multiple versions of Python are
installed on the same system:
\begin{lyxcode}
\$~python~-V Python~2.2.3 \$~which~python

/usr/bin/python

\$~ls~/usr/bin/python{*} /usr/bin/python~/usr/bin/python2~/usr/bin/python2.2 /usr/bin/python24~/usr/bin/python2.4 \$
\end{lyxcode}
In the above scenario, it is useful to create an alias to the newer
Python:
\begin{lyxcode}
\$~cd~\textasciitilde{}/bin \$~ln~-s~/usr/bin/python2.4~python

\$~export~PATH=\$HOME/bin:\$PATH \$~cd

\$~hash~-r \$~which~python

\textasciitilde{}/bin/python

\$~python~-V Python~2.4.1 \$
\end{lyxcode}
If absolutely necessary, one can easily build Python from source using
a C compiler. You can download Python from the Python website \url{www.python.org}.
Please note that Pyre in CitcomS does not work with Python 2.7 or
higher.

\subsection{Pythia (Optional)}

The current version of pythia package can be downloaded from the \href{http://www.geodynamics.org/cig/software/pythia}{CIG Pythia page}.

To obtain CitcomS, go to the Geodynamics Software Packages web page
the source archive and unpack it using the \texttt{tar} command:
\begin{lyxcode}
\$~tar~xzf~CitcomS-3.2.0.tar.gz \end{lyxcode} If you don't have GNU Tar, try the following command instead: \begin{lyxcode} \$~gunzip~-c~CitcomS-3.2.0.tar.gz~|~tar~xf~-
\end{lyxcode}

\section{Installation Procedure}

After unpacking the source, use the following procedure to install
CitcomS:
\begin{enumerate}
\item Navigate (i.e., \texttt{cd}) to the directory containing the CitcomS
source\texttt{.}~\\
\texttt{}~\\
\texttt{\$cd CitcomS-3.2.0} \item Type .\texttt{/configure} to configure the package for your system\texttt{.}~\\ \texttt{}~\\ \texttt{\$ ./configure}
\item Type \texttt{make} to build the package.\texttt{}~\\
\texttt{}~\\
\texttt{\$make} \end{enumerate} If you are content to run CitcomS from the build directory, then you are done. Upon successful completion, the \texttt{make} command creates a script called \texttt{citcoms} in the \texttt{bin} subdirectory; this is the script you will use to run CitcomS. You may wish to add the \texttt{bin} directory to your \texttt{PATH}. For more details about \texttt{configure}, see Section \ref{sec:Configuration} below. \subsection{Installing to a Secondary Location (Optional)} Optionally, after building CitcomS, you can install it in a secondary location using the \texttt{make install} command. This is not necessary for using CitcomS and is not recommended for most situations. It is documented only for the sake of completeness. By default, CitcomS is configured to install under \texttt{/usr/local}, which is not writable by normal users. To install as an ordinary user, give \texttt{make install} the \texttt{prefix} option, specifying a directory to which you have write access: \begin{lyxcode} \$~make~install~prefix=\$HOME/cig \end{lyxcode} The above command will install CitcomS under \texttt{\$HOME/cig}.
Afterwards, you may wish to add \texttt{PREFIX/bin} (\texttt{\$HOME/cig/bin}, in this example) to your \texttt{PATH}. After running \texttt{make install}, you may (if desired) run \texttt{make clean} in the build directory to save disk space. You are also free to delete the source/build directory altogether. (Note that \texttt{make install} installs the examples under \texttt{PREFIX/share/CitcomS/examples}.) \section{\label{sec:Configuration}Configuration} The \texttt{configure} script checks for various system features. As it runs, it prints messages informing you of which features it is checking for. Upon successful completion, it generates a \texttt{Makefile} in each source directory of the package. It also generates a \texttt{config.h} header file, which contains system-dependent definitions. The \texttt{configure} script will attempt to guess the correct values of various installation parameters. In the event that the default values used by \texttt{configure} are incorrect for your system, or \texttt{configure} is unable to guess the value of a certain parameter, you may have to specify the correct value by hand. \begin{quote} \textbf{Important:} If the \texttt{configure} script fails, and you don't know what went wrong, examine the log file \texttt{config.log}. This file contains a detailed transcript of all the checks \texttt{configure} performed. More importantly, it includes the error output (if any) from your compiler. When seeking help for \texttt{configure} failures on the CIG Mantle Convection Mailing List \url{cig-mc@geodynamics.org}, please send \texttt{config.log} as an attachment. \end{quote} Upon successful completion, \texttt{configure} will print a brief configuration summary. If Pyre is enabled, the \texttt{configure} script will automatically check for needed Python dependencies, including the Pythia package (which includes the Pyre framework). If necessary, \texttt{configure} will download missing Python packages from the Internet and store them under \texttt{\textit{\emph{deps/}}}. These are Python packages which are required (either directly or indirectly) by CitcomS. \subsection{Configure Usage} For a detailed list of \texttt{configure} variables and options, give \texttt{configure} the \texttt{-{}-help} option: \begin{lyxcode} \$~./configure~-{}-help
\end{lyxcode}
The following is a summary of the variables and options that are important
when installing CitcomS.

\subsection{Environment Variables}

Environment variables may be specified as arguments to \texttt{configure},
e.g.,
\begin{lyxcode}
\$~./configure~CC=icc~\#~use~the~Intel~compiler \end{lyxcode} \noindent \begin{center} \begin{table}[H] \noindent \centering{}% \begin{tabular}{|>{\centering}p{1in}|>{\raggedright}p{3.5in}|} \hline \textbf{Variable} & \textbf{Description}\tabularnewline \hline \hline PYTHON & Python interpreter. This is useful if you have Python installed in a non-standard location, e.g., \texttt{~~./configure \textbackslash{}} \texttt{~~~~~~PYTHON=/opt/python2.5/bin/python} By default, \texttt{configure} will search for a suitable Python interpreter using your \texttt{PATH} environment variable.\tabularnewline \hline CC & C compiler command. This is usually set to the name of an MPI wrapper command, e.g., \texttt{~~./configure \textbackslash{}} \texttt{~~~~~~CC=/opt/mpich-1.2.6/bin/mpicc} See Section \vref{sec:MPI-Configuration} for details and examples.\tabularnewline \hline CFLAGS & C compiler flags. This can be set for optimization flags.\tabularnewline \hline CPPFLAGS & C preprocessor flags; e.g., \texttt{-I<dir>} if you have headers in a nonstandard directory.\tabularnewline \hline LDFLAGS & Linker flags; e.g., \texttt{-L<dir>} if you have libraries in a nonstandard directory.\tabularnewline \hline LIBS & Library flags; e.g., -l<lib> if additional library is needed\tabularnewline \hline \end{tabular} \end{table} \par\end{center} \subsection{\label{sec:MPI-Configuration}MPI Configuration} By default, \texttt{configure} will search for a C compiler using your \texttt{PATH} environment variable. It prefers MPI wrapper commands (such as \texttt{mpicc}) to ordinary compiler commands (such as \texttt{cc} or \texttt{gcc}). You may specify the compiler command name manually using the \texttt{CC} variable: \begin{lyxcode} \$~./configure~CC=/opt/mpich-1.2.6/bin/mpicc
\end{lyxcode}
The \texttt{configure} script will test for the presence of the MPI
header (\texttt{mpi.h}) and an MPI library using the C compiler command.
If \texttt{CC} is set to an MPI wrapper command such as \texttt{mpicc},
and/or the MPI header and library files are installed in a standard
location (i.e., \texttt{/usr/include} and \texttt{/usr/lib}), these
\texttt{configure} tests should succeed without difficulty.

But if CC is set to an ordinary compiler command name (e.g., \texttt{cc}
or \texttt{gcc}) and MPI is installed in a non-standard location,
you must manually specify \texttt{CPPFLAGS} and \texttt{LDFLAGS},
so that the compiler can find the MPI header files and libraries.

\subsubsection{Manually Specifying MPI \texttt{include} and \texttt{lib} Directories}
\begin{lyxcode}
\$~./configure~\textbackslash{} CPPFLAGS=\textquotedbl{}-I/opt/mpich-1.2.6/include\textquotedbl{}~\textbackslash{} LDFLAGS=\textquotedbl{}-L/opt/mpich-1.2.6/lib~-lmpich\textquotedbl{} \end{lyxcode} \subsubsection{Manually Specifying MPI \texttt{include} and \texttt{lib} Directories and an Alternative Compiler} \begin{lyxcode} \$~./configure~\textbackslash{}

CC=icc~\textbackslash{}

CPPFLAGS=\textquotedbl{}-I/opt/mpich-1.2.6/include\textquotedbl{}~\textbackslash{}

LDFLAGS=\textquotedbl{}-L/opt/mpich-1.2.6/lib~-lmpich\textquotedbl{}
\end{lyxcode}
Note that it may be necessary to specify the name of the MPI library
itself in \texttt{LDFLAGS} using the \texttt{-l} compiler option.
If a library name is not given -- or if the given option doesn't work
-- \texttt{configure} will automatically try linking using \texttt{-lmpi}
and, if that fails, \texttt{-lmpich}.

\section{\label{sec:Exchanger-Configuration-(Optional)}Exchanger Configuration
(Optional)}

This version of CitcomS is capable of solver coupling. Two or more
instances of CitcomS solvers can be coupled together to solve a problem
with different length and time scales \cite{Tan et al GeoFramework Part I}.
An external Exchanger package is needed for solver coupling. You can
After you download and untar the package, you can configure and install
Exchanger using these commands:
\begin{lyxcode}
\$~cd~Exchanger-1.0.1/ \$~./configure~-{}-prefix=\$HOME/cig \$~make

\$~make~install \end{lyxcode} By default, CitcomS will attempt to auto-detect your Exchanger installation and will disable Exchanger support if it is not found. You may specify the location of your Exchanger installation by setting the PYTHONPATH environment variable to the appropriate installation prefix. \begin{lyxcode} \$~cd~CitcomS-3.2.0/

\$~export~PYTHONPATH=\$HOME/cig/lib/python2.5/site-packages:\$PYTHONPATH \$~./configure~-{}-with-exchanger
\end{lyxcode}

\section{\label{sec:HDF5-Configuration}HDF5 Configuration (Optional)}

For writing its output in binary format, CitcomS requires parallel
HDF5 (PHDF5). In turn, PHDF5 requires an MPI library with MPI-IO support
and a parallel file system. If an existing installation of the PHDF5
library is not available on your cluster, you can compile it from
source by following the instructions in the file \texttt{\small{release\_docs/INSTALL\_parallel}}
under the HDF5 source tree. Under Debian Linux, you may simply install
the \texttt{libhdf5-mpich}, \texttt{libhdf5-mpich-dev} and \texttt{hdf5-tools}
packages.
\begin{quote}
\textbf{NOTE:} Most post-processing and visualization scripts discussed
in Chapter \ref{cha:Postprocessing-and-Graphics} only understand
ASCII output. An additional step to convert the HDF5 output files
to ASCII files in necessary before using those scripts.
\end{quote}
By default, CitcomS will attempt to auto-detect your PHDF5 installation,
and will disable HDF5 support if it is not found. You may specify
the location of your PHDF5 installation by setting the PHDF5\_HOME
environment variable to the appropriate installation prefix.
\begin{lyxcode}
\$~export~PHDF5\_HOME=/opt/phdf5/1.6.5 \$~./configure~-{}-with-hdf5
\end{lyxcode}

While the following software is not necessary for the normal operation
of CitcomS, you may find it useful for accessing CitcomS data in HDF5
files.

\subsubsection{NumPy}

NumPy is an extension to Python which adds support for multi-dimensional
extracting it:
\begin{lyxcode}
\$~cd~numpy-1.0 \$~python~setup.py~install~-{}-prefix=\$HOME/cig \end{lyxcode} Alternatively, under Debian Linux you can install the \texttt{python-numpy} package. On Gentoo Linux, NumPy is available in the \texttt{dev-python/numpy} ebuild. \subsubsection{\label{sub:PyTables}PyTables} PyTables is an extension to Python and can expose HDF5 array datasets as Python NumPy arrays. It is available at PyTables \url{www.pytables.org}. To compile and install this extension, download the latest stable version and issue the following commands: \begin{lyxcode} \$~cd~pytables-1.3.3

\$~python~setup.py~install~-{}-prefix=\$HOME/cig
\end{lyxcode}
To install on Debian Linux, you may use the \texttt{python-tables}
package instead. On Gentoo Linux, it is available in the \texttt{dev-python/pytables}
ebuild.

\subsubsection{HDFView}

HDFView is a visual tool written in Java for browsing and editing

\subsubsection{\label{sub:OpenDXutils}OpenDXutils}

In order to import HDF5 files into OpenDX, you need the OpenDXutils
package from the Cactus project. Go to the OpenDXutils package website
to download and install the package. Note that you will need to set
both \texttt{DXMODULES} and \texttt{DXMDF} environment variables before
running OpenDX to load the package.

\section{GGRD Configuration (Optional)}

CitcomS can read data from files for initial temperature, boundary
conditions, material control and others. By default, these files are
in ASCII format, and the data in the files must be ordered in the
same way of the node ordering. As a result, whenever the mesh is changed,
these data files need to be regenerated to satisfy the ordering requirement.
This restriction can be lifted if using GRD data files. GRD files
are binary files in NetCDF format (Network Common Data Form). GRD
files are widely used in GMT (Generic Mapping Tools) as well.

The GRD file support can be enabled at configuration time. It requires
the installation of HC package \url{geodynamics.org/cig/software/packages/mc/hc},
a 1D global mantle circulation solver which can compute velocities,
tractions, and geoid for simple density distributions and plate velocities.
To install HC, you will also need the GMT and NetCDF packages (library
on installation. After HC is installed, note the directory (for example:
\texttt{\$HOME/cig/HC-1\_0}) of HC installation, and execute these commands: \begin{lyxcode} \$~export~HC\_HOME=\$HOME/cig/HC-1\_0 \$~./configure~-{}-with-ggrd
\end{lyxcode}
You may need to set \texttt{GMTHOME} and \texttt{NETCDFHOME} environment
variables if these packages is not installed in the system directory.

\section{PETSc Configuration (Optional)}

To use the PETSc functionality in CitcomS, CitcomS first needs to be built with the PETSc library. This
is done by passing \texttt{--with-petsc=yes} option to the \texttt{configure} script. If this option is
passed to \texttt{configure}, \texttt{configure} expects PETSc to be installed correctly on the system,
with the environment variables \texttt{PETSC\_DIR} and \texttt{PETSC\_ARCH} being set appropriately.
\texttt{configure} will check for the \texttt{PETSC\_DIR} and \texttt{PETSC\_ARCH} only if the
\texttt{--with-petsc=yes} option has been passed to the script.

The PETSc homepage is \url{http://www.mcs.anl.gov/petsc/} and information about installing PETSc can be found at
\url{http://www.mcs.anl.gov/petsc/documentation/installation.html}

\section{\label{sec:Batch-System-Configuration}Batch System Configuration}

If you are installing CitcomS on a cluster with a batch system, you
can configure Pyre such that the \texttt{citcoms} command automatically
submits jobs to the batch queue. Pyre contains support for the LSF,
PBS, and Globus batch systems.

The command to submit a batch job depends upon the particular batch
system used. Further, the command used in a batch script to launch
an MPI program varies from one cluster to the next. This command can
vary between two clusters, even if the clusters use the same batch
system! On some systems, \texttt{mpirun} is invoked directly from
the batch script. On others, a special wrapper is used instead. If
you need help to configure Pyre for your cluster, you can contact
CIG Mantle Convection Mailing List \url{cig-mc@geodynamics.org} for
assistance.

Properly configured, Pyre can handle job submissions automatically,
insulating users from the details of the batch system and the site
configuration. This feature has the most value when the system administrator
installs a global Pyre configuration file on the cluster (under \texttt{/etc/pythia-0.8}),
for the benefit of all users and all Pyre-based applications.

CIG's Pythia page \url{geodynamics.org/cig/software/packages/cs/pythia/docs/batch}.
to running CitcomS, see Section \vref{sub:Launchers-and-Schedulers}.

\section{\label{sec:Installing-without-Pyre}Installing with Pyre}

By default, CitcomS 3.2.0 does not install using Pyre. To build CitcomS
binaries with the Pyre wrappers (necessary for the example Cookbooks),
call \texttt{configure} using the \texttt{-{}-with-pyre} option:
\begin{lyxcode}
\$~./configure~-{}-with-pyre=yes \end{lyxcode} The only system requirements for this configuration are an MPI library and a C compiler. With Pyre the make command will build \texttt{pycitcoms} and \texttt{mpipycitcoms} as well as the script \texttt{citcoms} used for running CitcomS. Without Pyre the \texttt{make} command will build two command-line tools, \texttt{CitcomSFull} and \texttt{CitcomSRegional}, for running the full solver and the regional solver, respectively. \section{\label{sec:Software-Repository}Installing from the Software Repository} The CitcomS source code is available via a Subversion server at the Geodynamics website \url{geodynamics.org}. This allows users to view the revision history of the code and check out the most recent development version of the software. \begin{quote} \textbf{NOTE:} If you are content with the prepared source package, you may skip this section. \end{quote} \subsection{Tools You Will Need} In addition to the usual system requirements, you will need a handful of additional development tools installed in order to work with the source from the CIG software repository. First, you must have a Subversion client installed. To check, type \texttt{svn}; it should return a usage message. \begin{lyxcode} \$~svn

Type~'svn~help'~for~usage.
\end{lyxcode}

Second, you must have the GNU tools Autoconf, Automake, and Libtool
installed. To check, enter the following commands:
\begin{lyxcode}
\$~autoconf~-{}-version \$~automake~-{}-version

\$~libtoolize~-{}-version \end{lyxcode} For more information about these GNU tools, see the GNU website \url{www.gnu.org/software}. The CitcomS v3.2.0 source package was created with Autoconf 2.68, Automake 1.11.2, and Libtool 2.4.2. \subsection{Download Source from Subversion} To check out the latest version of the software, use the \texttt{svn checkout} command: \begin{lyxcode} \$~svn~checkout~http://geodynamics.org/svn/cig/mc/3D/CitcomS/trunk~CitcomS
\end{lyxcode}
This will create the local directory \texttt{CitcomS} (if it doesn't
already exist) and fill it with the latest CitcomS source from the
CIG software repository.

The \texttt{CitcomS} directory thus created is called a \emph{working
copy}. To merge the latest changes into an existing working copy,
use the \texttt{svn update} command:
\begin{lyxcode}
\$~cd~CitcomS \$~svn~update
\end{lyxcode}
This will preserve any local changes you have made to your working
copy.

\subsection{Generating the GNU Build System}

Your working directory should now contain a fresh checkout of CitcomS:
\begin{lyxcode}
\$~ls CitcomS \$
\end{lyxcode}
Before you can run \texttt{configure} or \texttt{make}, you must generate
the necessary files using the GNU tools. The easiest way to do this
is to run \texttt{autoreconf -i}:
\begin{lyxcode}
\$~cd~CitcomS \$~autoreconf~-i

\$~./configure \$~make
\end{lyxcode}
The \texttt{autoreconf} tool runs Autoconf to generate the \texttt{configure}
script from \texttt{configure.ac}. It also runs Automake to generate
\texttt{Makefile.in} from \texttt{Makefile.am} in each source directory.

The configure script stores dependency information for source files
in hidden subdirectories called \texttt{.deps}. If source files are
added, deleted, moved, or renamed -- which may happen if you run \texttt{svn
update} to merge the latest changes -- the \texttt{make} command may
report errors for source files which no longer exist. If this happens,
clean your existing configuration using the \texttt{make distclean}
command, and then re-run \texttt{configure} and \texttt{make}:
\begin{lyxcode}
\$~make~distclean \$~./configure

\$~make \end{lyxcode} The \texttt{make} \texttt{distclean} command deletes all files generated by \texttt{configure}, including the Makefiles themselves! Therefore, after running \texttt{make} \texttt{distclean}, you will not be able to run \texttt{make} again until you re-run \texttt{configure}. \chapter{Running CitcomS} \section{Using CitcomS without Pyre} Regardless of whether you build CitcomS with or without Pyre, two binary executables, \texttt{CitcomSRegional} and \texttt{CitcomSFull}, are always placed under the \texttt{bin} directory. These programs are compiled from pure C code and do not use Python or the Pyre framework. Each program has the same usage: \begin{lyxcode} \$~mpirun~{[}mpi\_options{]}~CitcomSRegional~inputfile~\\
\$~mpirun~{[}mpi\_options{]}~CitcomSFull~inputfile \end{lyxcode} Two input file examples, one for a regional spherical model and one for a full spherical model, are provided in the \texttt{examples/Regional} and \texttt{examples/Full} directories, respectively. The meaning of the input parameters is described in Appendix \vref{cha:Appendix-A:-Input}. The pure C version, compared to the Pyre version, shares the same input parameters and functionality, but is less flexible in changing parameters and in launching parallel jobs. Users are encouraged to use the Pyre version if possible. In the following sections and Chapter \ref{cha:Cookbooks}, we will concentrate on the Pyre version only. \section{Using CitcomS with Pyre} If you build CitcomS with the Pyre framework, an additional executable, \texttt{citcoms}, is placed under the \texttt{bin} directory. The \texttt{citcoms} executable is a Python script used for running both the regional and full spherical models using Pyre. Executed without any command line options, \texttt{citcoms} will run a regional model with default parameters. It can also run a full spherical model if the correct parameters are set (see Section \vref{sec:Cookbook-1:-Global} for an example). On input, CitcomS needs numerous parameters to be specified (see Appendix \ref{cha:Appendix-A:-Input} for a full list). All parameters have sensible default values. Since you will likely want to specify the parameters of your CitcomS runs, you will need to alter both computational details (such as the number of time steps) and controlling parameters specific to your problem (such as the Rayleigh number). These input parameters, or properties in the Pyre terminology, are grouped under several Pyre components. Most of the properties you will set using CitcomS have names which are identical to the parameters for CitcomS in pure C version, which are described in Appendix \ref{cha:Appendix-A:-Input}. \section{Changing Parameters } There are several methods to set the input parameters for CitcomS: via the command line, or by using a configuration file in \texttt{.cfg} format. \subsection{Using the Command Line} Pyre uses the following syntax to change properties from the command line. To change the value of a property of a component, use: \begin{lyxcode} -{}-{[}component{]}.{[}property{]}={[}value{]} \end{lyxcode} Each component is attached to a facility, so the option above can also be written as: \begin{lyxcode} -{}-{[}facility{]}.{[}property{]}={[}value{]} \end{lyxcode} Each facility has a default component attached to it. A different component can be attached to a facility by: \begin{lyxcode} -{}-{[}facility{]}={[}new\_component{]}~ \end{lyxcode} \subsection{Using a \texttt{.cfg} File} Entering all those parameters via the command line involves the risk of typographical errors, which can lead to undesired results. You may find it easier to write a brief \texttt{.cfg} input file that contains the parameters. This file has a format similar to a Windows INI file. The file is composed of one or more sections which are formatted as follows: \begin{lyxcode} {[}CitcomS.component1.component2{]}~\\ \#~this~is~a~comment~\\ property1~=~value1~\\ property2~=~value2~~;~this~is~another~comment \end{lyxcode} We strongly recommend that you use \texttt{.cfg} files for your work. The files are syntax-colored by the \texttt{vim} editor. (Upon termination of each run, all of the parameters are logged in a \texttt{.cfg} file.) \subsection{Using a \texttt{.pml} File} A \texttt{.pml} file is an XML file that specifies parameter values in a highly structured format. XML files are intended to be read and written by machines, not edited manually by humans. The \texttt{.pml} file format is intended for applications in which CitcomS input files are generated by another program, e.g., a GUI, web application, or a high-level structured editor. This file fomat will not be discussed or used further in the manual. It is composed of nested sections which are formatted as follows: \begin{lyxcode} <component~name='component1'>~\\ ~~~~<component~name='component2'>~\\ ~~~~~~~~<property~name='property1'>value1</property>~\\ ~~~~~~~~<property~name='property2'>value2</property>~\\ ~~~~</component>~\\ </component> \end{lyxcode} \subsection{Specification and Placement of Configuration Files} One or more configuration files and command line options may be specified on the command line: \begin{lyxcode} \$~citcoms~a.cfg~b.cfg~-{}-foo.bar=baz
\end{lyxcode}
In addition, the Pyre framework searches for configuration files named
\texttt{CitcomS.cfg} in several predefined locations. You may put
settings in any or all of these locations, depending on the scope
you want the settings to have:
\begin{enumerate}
\item \texttt{PREFIX/etc/CitcomS.cfg}, for system-wide settings;
\item \texttt{\$HOME/.pyre/CitcomS/CitcomS.cfg}, for user settings and preferences; \item the current directory (\texttt{./CitcomS.cfg}), for local overrides. \end{enumerate} Parameters given directly on the command line will override any input contained in a configuration file. If more than one configuration files are given on the command line, the later one overrides the former ones. Configuration files given on the command line override all other \texttt{CitcomS.cfg} files. The \texttt{CitcomS.cfg} files placed in (3) will override those in (2), (2) overrides (1), and (1) overrides only the built-in defaults. \section{Uniprocessor Example} CitcomS runs similarly in full spherical or regional modes. For the purpose of this example, you will perform a test run of the regional version on a workstation. Execute the following on the command line: \begin{lyxcode} \$~citcoms~-{}-steps=10~-{}-controller.monitoringFrequency=5~\textbackslash{}~\\
-{}-solver.datafile=example0~-{}-solver.mesher.nodex=17~-{}-solver.mesher.nodey=17
\end{lyxcode}
This runs a default convection problem in a regional domain for 10
time steps and with a mesh of 17 $\times$ 17 $\times$ 9 nodal points.
Since we did not provide the parameter \texttt{solver.mesher.nodez},
the default value 9 is used. The model results are written to files
\texttt{example0.{*}} with an interval of 5 time steps.

Instead of writing the input parameters on the command line, you can
put them in a \texttt{.cfg} file. The CitcomS source package contains
an \texttt{examples} directory (the \texttt{make install} command
installs the examples under \texttt{PREFIX/share/CitcomS/examples},
where \texttt{PREFIX} is the \texttt{CitcomS} installation directory).
In this directory, you will find a configuration file equivalent to
the previous example: \texttt{example0.cfg}. You can run the model
using:
\begin{lyxcode}
\$~citcoms~example0.cfg \end{lyxcode} \subsubsection{Example: Uniprocessor, \texttt{example0.cfg}} \begin{lyxcode} {[}CitcomS{]}~\\ steps~=~5~\\ ~\\ {[}CitcomS.controller{]}~\\ monitoringFrequency~=~1~\\ ~\\ {[}CitcomS.solver{]}~\\ datafile~=~example0~\\ ~\\ {[}CitcomS.solver.mesher{]}~\\ nodex~=~~17~\\ nodey~=~~17 \end{lyxcode} \section{\label{sec:Multiprocessor-Example}Multiprocessor Example} In order to run this example, you should be on a Beowulf cluster with four or more processors, or on a supercomputer; and you should be in the directory in which the input file is located, in this case, the \texttt{examples} directory. CitcomS has been extensively used on both environments, using up to several hundred processors. How to run a multiprocessor CitcomS model depends on your hardware and software settings, e.g., whether a batch system is used, what the names of the computers in a cluster are, and how the file system is organized. This section will lead you through the different settings of a parallel environment. \subsubsection{Example: Multiprocessor, \texttt{example1.cfg }} \begin{lyxcode} {[}CitcomS{]}~\\ steps~=~70~\\ ~\\ {[}CitcomS.controller{]}~\\ monitoringFrequency~=~10~\\ ~\\ {[}CitcomS.solver{]}~\\ datafile~=~example1~\\ ~\\ {[}CitcomS.solver.mesher{]}~\\ nprocx~=~~2~\\ nprocy~=~~2~\\ nodex~~=~17~\\ nodey~~=~17~\\ nodez~~=~~9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \end{lyxcode} This example uses 2 processors in colatitude ($x$-coordinate), 2 in longitude ($y$-direction), and 1 in the radial ($z$-direction), i.e., it uses 4 processors in total. In addition, there will be 17 nodes in$x$(theta), 17 nodes in$y$(fi), and 9 nodes in$z$(radius\_inner). The model will run for 70 time steps and the code will output the results every 10 time steps. It is important to realize that within the example script (and in finite element method, FEM) the term node'' refers to the mesh points defining the corners of the elements. In \texttt{example1.cfg}, this is indicated with: \begin{lyxcode} nodex~~=~17~\\ nodey~~=~17~\\ nodez~~=~~9~ \end{lyxcode} These quantities refer to the total number of FEM nodes in a given direction for the complete problem, and for the example it works out that within a given processor there will be 9$\times$9$\times$9 nodes. Note that in the$x$-direction (or$y$) that for the entire problem there are 17 nodes and there is one node shared between two processors. This shared node is duplicated in two adjacent processors. Unfortunately, nodes'' sometimes also refer to the individual computers which make up a Beowulf cluster or supercomputer. In the example scripts, this is indicated with: \begin{lyxcode} nprocx~=~~2~\\ nprocy~=~~2 \end{lyxcode} \begin{figure}[H] \begin{centering} \includegraphics[width=0.5\paperwidth]{graphics/cookbook2.pdf} \par\end{centering} \caption{Computational Domain. Map view on the configuration of the top layer of the computational nodes and the processors.} \end{figure} \subsection{Output Directories and Output Formats} CitcomS potentially generates a large number of ASCII files. This means that you will have to organize your directories carefully when running CitcomS so that you can manage these files as well as use a post-processing program contained in this distribution. How to best manage this large output depends on whether you will use a local file system or a parallel file system. For example, if you have a local hard disk on every machine (node) on a Beowulf cluster, with each hard disk mounted locally to the machine, this scenario is referred to as a local file system in this section. Or you might use some kind of parallel file system on your computer (e.g., NFS, GPFS, PVFS, to name a few), which is mounted on all of the nodes. Usually your home directory is mounted on the parallel file system. The local file system is usually more cost- and time-efficient than the parallel file system. If you want CitcomS to write its output to the local hard disks, you need to have a common directory structure on all of the local hard disks. For example, if the directory \texttt{/scratch} exists on all local hard disks, you can run the example script with: \begin{lyxcode} \$~citcoms~example1.cfg~-{}-solver.datadir=/scratch
\end{lyxcode}
property, which specifies the output directory. The output files are
then placed in \texttt{/scratch} on each individual machine with a
filename prefix \texttt{example1.}

However, if the output directory name on each local hard disk depends
on the machine hostname, you can run the example script with:
\begin{lyxcode}
\$~citcoms~example1.cfg~-{}-solver.datadir=/scratch\_\%HOSTNAME \end{lyxcode} The special string \texttt{\%HOSTNAME} will be substituted by the hostname of each machine. As the final example for a local file system, you can specify an arbitrary output directory for each machine. To do so, you must write a program to be executed on each machine which will print the output directory. The program must be named \texttt{citcoms\_datadir} and must reside on your path. An example of \texttt{citcoms\_datadir} can be found in the \texttt{visual/} directory. Then you can run the example script with: \begin{lyxcode} \$~citcoms~example1.cfg~-{}-solver.datadir=\%DATADIR
\end{lyxcode}
The special string \texttt{\%DATADIR} will be substituted by the output

If you want CitcomS to write its output to a parallel file system,
you have several choices. You can run the example script as follows
\begin{lyxcode}
\$~citcoms~example1.cfg~-{}-solver.datadir=/home/\emph{username} \end{lyxcode} The output files are then placed in your home directory with a filename prefix \texttt{example1}. A potential problem with this approach is that the directory \texttt{/home/}\texttt{\emph{username}} will be flooded with hundreds of files, perhaps even tens of thousands of files if you are running a model using several tens of processors for thousands of time steps. Alternatively you can have each machine write its output to its own directory, according to its MPI rank. You can run the example script with: \begin{lyxcode} \$~citcoms~example1.cfg~-{}-solver.datadir=/home/\emph{username}/\%RANK
\end{lyxcode}
The special string \texttt{\%RANK} will be substituted by the MPI
rank of each processor. You will see four new directories \texttt{/home/}\texttt{\emph{username}}\texttt{/0},
of MPI rank 0 will write its output in \texttt{/home/}\texttt{\emph{username}}\texttt{/0}
with a filename prefix \texttt{example1} (defined by the property
\texttt{datafile} inside \texttt{example1.cfg}) and so on.

The ASCII output can potentially take a lot of disk space. CitcomS
can write \texttt{gzip} compressed output directly. You can run the
example script with:
\begin{lyxcode}
\$~citcoms~example1.cfg~-{}-solver.datadir=/home/\emph{username}~\textbackslash{}~\\ -{}-solver.output.output\_format=ascii-gz \end{lyxcode} Be warned that the post-process scripts do not understand this output format yet. The last choice is the most powerful one. Instead of writing many ASCII files, CitcomS can write its results into a single HDF5 (Hierarchical Data Format) file per time step. These HDF5 files take less disk space than all the ASCII files combined and don't require additional post-processing to be visualized in OpenDX. In order to use this feature, you must compile CitcomS with the parallel HDF5 library if you haven't done so already (see Section \vref{sec:HDF5-Configuration}). You can run the example script with: \begin{lyxcode} \$~citcoms~example1.cfg~-{}-solver.datadir=/home/\emph{username}~\textbackslash{}~\\
-{}-solver.output.output\_format=hdf5
\end{lyxcode}
The output files will be stored in \texttt{/home/}\texttt{\emph{username}}\texttt{/}
with a filename prefix \texttt{example1} and a filename suffix \texttt{h5}.
how to work with the HDF5 output.

\subsection{\label{sub:Launchers-and-Schedulers}Launchers and Schedulers}

If you have used MPI before, you know that \texttt{mpirun} requires
several command-line options to launch a parallel job. Or if you have
used one of the batch systems, you will know that the batch system
requires you to write a script to launch a job. Fortunately, launching
a parallel CitcomS job is simplified by Pyre's \texttt{launcher} and
\texttt{scheduler} facilities. Many properties associated with \texttt{launcher}
and \texttt{scheduler} are pertinent to the cluster you are on, and
are best customized in a configuration file. Your personal CitcomS
configuration file (\texttt{\textasciitilde{}/.pyre/CitcomS/CitcomS.cfg})
is suitable for this purpose. On a cluster, the ideal setup is to
install a system-wide configuration file under \texttt{/etc/pythia-0.8},
for the benefit of all users.

Pyre's \texttt{scheduler} facility is used to specify the type of
batch system you are using (if any):
\begin{lyxcode}
{[}CitcomS{]}

scheduler~=~lsf
\end{lyxcode}
The valid values for \texttt{scheduler} are \texttt{lsf}, \texttt{pbs},
\texttt{globus}, and \texttt{none}.

Pyre's \texttt{launcher} facility is used to specify which MPI implementation
you are using:
\begin{lyxcode}
{[}CitcomS{]}

launcher~=~mpich
\end{lyxcode}
The valid values for \texttt{launcher} include \texttt{mpich} and
\texttt{lam-mpi}.

You may find the \texttt{dry} option useful while debugging the \texttt{launcher}
and \texttt{scheduler} configuration. To debug the scheduler configuration,
use the \texttt{-{}-scheduler.dry} option:
\begin{lyxcode}
citcoms~-{}-scheduler.dry
\end{lyxcode}
This option will cause CitcomS to perform a dry run,'' dumping the
batch script to the console, instead of actually submitting it for
execution (the output is only meaningful if you're using a batch system).
Likewise, to debug the launcher configuration, use the \texttt{-{}-launcher.dry}
option:
\begin{lyxcode}
citcoms~-{}-launcher.dry
\end{lyxcode}
This option will cause CitcomS to print the \texttt{mpirun} command,
instead of actually executing it. (If you're using a batch system,
a job will be submitted for execution; when it runs, CitcomS will
simply print the \texttt{mpirun} command, and the job will immediately
terminate.)

\subsubsection{Running without a Batch System}

On a cluster without a batch system, you need to specify on which
machines the job will run. Supposing the machines on your cluster
are named n001, n002, \ldots{}, etc., but you want to run the job
on machines n001, n003, n004, and n005 (maybe n002 is down for the
moment). To run the example, create a file named \texttt{mymachines.cfg}
which specifies the machines to use:
\begin{lyxcode}
{[}CitcomS.launcher{]}

nodegen~=~n\%03d

nodelist~=~{[}1,3-5{]}
\end{lyxcode}
The \texttt{nodegen} property is a printf-style format string, used
in conjunction with \texttt{nodelist} to generate the list of machine
names. The \texttt{nodelist} property is a comma-separated list of
machine names in square brackets.

Now, invoke the following:
\begin{lyxcode}
\$~citcoms~example1.cfg~mymachines.cfg \end{lyxcode} This strategy gives you the flexibility to create an assortment of \texttt{.cfg} files (with one \texttt{.cfg} file for each machine list) which can be easily paired with different parameter files. If your machine list does not change often, you may find it more convenient to specify default values for \texttt{nodegen} and \texttt{nodelist} in \texttt{\textasciitilde{}/.pyre/CitcomS/CitcomS.cfg} (which is read automatically). Then, you can run any simulation with no additional arguments: \begin{lyxcode} \$~citcoms~example1.cfg\end{lyxcode}
\begin{quote}
\textcolor{red}{Warning:} This assumes your machine list has enough
nodes for the simulation in question.
\end{quote}
You will notice that a machine file \texttt{mpirun.nodes} is generated.
It will contain a list of the nodes where CitcomS has run. Save the
machine file as it will be useful in the postprocessing step.

\subsubsection{Using a Batch System}

The settings which are important when using a batch system are summarized
in the sample configuration file which follows.
\begin{lyxcode}
{[}CitcomS{]}

scheduler~=~lsf~~~~;~the~type~of~the~installed~batch~system~\\
~\\
{[}CitcomS.lsf{]}

bsub-options~=~{[}-a~mpich\_gm{]}~~~~;~special~options~for~'bsub'~\\
~\\
{[}CitcomS.launcher{]}

command~=~mpirun.lsf~~~~;~'mpirun'~command~to~use~on~our~cluster~\\
~\\
{[}CitcomS.job{]}

queue~=~normal~~~~~~~~~~;~default~queue~for~jobs~\\
walltime~=~5{*}minute~~~~~;~run~time~limit~of~the~job
\end{lyxcode}
These settings are usually placed in \texttt{\textasciitilde{}/.pyre/CitcomS/CitcomS.cfg}
or in a system-wide configuration file. They can be overridden on
the command line, where one typically specifies the job name and the
allotted time for the job:
\begin{lyxcode}
\$~citcoms~example1.cfg~-{}-job.queue=debug~\textbackslash{} ~~~~-{}-job.name=example1~-{}-job.walltime=5{*}minute \end{lyxcode} The number of nodes to allocate for the job is determined automatically, based upon the simulation parameters. \subsection{Monitoring Your Jobs} Once launched, CitcomS will print the progress of the model to the standard error stream (stderr). Usually, the stderr is directed to your terminal so that you can monitor the progress. On some system, the stderr is redirected to a file. In any case, the progress is always saved in a log file (e.g., \texttt{example1.log}). The log file contains the convergence progress of the computation and, if an error occurs, debugging output. The time file (e.g., \texttt{example1.time}) contains the elapsed model time (in non-dimensional units) and CPU time (in seconds) of every time step. The format of the time file can be found in Appendix \vref{cha:Appendix-C:-CitComS,}. The log and time files are output by the rank-0 processor only. Following your successful run, you will want to retrieve the output files from all the nodes and process them so they can be visualized with the visualization program OpenDX (see Chapter\vref{cha:Postprocessing-and-Graphics}). \section{Using CitcomS on the TeraGrid} The TeraGrid is a set of parallel supercomputer facilities at eight partner sites in the U.S. which creates an integrated, persistent computational resource. Since TeraGrid software is based on commodity clusters, Linux/Unix, and Globus, it should be easier to scale from a laboratory development environment to a high-end environment in a straightforward manner which promotes application performance. Although the TeraGrid is a high-end resource, it was developed to be accessible to the general community of scientists and engineers as a production facility. TeraGrid accounts for small allocations are available directly from CIG for investigators in the U.S. CitcomS has already been installed and tested on several NSF TeraGrid platforms, including NCSA, SDSC, and TACC. To use CitcomS on these machines, please log in to your TeraGrid account and read the instructions at \texttt{\$TG\_COMMUNITY/CIG/CitcomS/TG\_README}. See CIG Community
Area Software on the TeraGrid \url{geodynamics.org/cig/software/csa/}
for additional information on access and to apply for allocation time.

\chapter{\label{cha:Working-with-HDF5}Working with CitcomS HDF5 Files}

\section{Introduction}

A typical run of CitcomS can create thousands if not millions of ASCII
output files. This situation is inefficient since it requires an extra
post-processing step for assembling the results from each processor
(see Chapter \vref{cha:Postprocessing-and-Graphics}). Since the v2.1
release, CitcomS solves this problem when running the software on
computers that have parallel file systems by assembling a binary HDF5
file in parallel I/O mode. You can skip this chapter if you are not
using the HDF5 output format.

The Hierarchical Data Format (HDF) is a portable file format developed
at the National Center for Supercomputing Applications (NCSA) \url{hdf.ncsa.uiuc.edu/HDF5}.
It is designed for storing, retrieving, analyzing, visualizing, and
converting scientific data. The current and most popular version is
HDF5, which stores multi-dimensional arrays together with ancillary
data in a portable self-describing format. It uses a hierarchical
structure that provides application programmers with a host of options
for organizing how data is stored in HDF5 files.

HDF5 files are organized in a hierarchical structure, similar to a
Unix file system. Two types of primary objects, \textit{groups} and
\textit{datasets}, are stored in this structure. A group contains
instances of zero or more groups or datasets, while a dataset stores
a multi-dimensional array of data elements. Both kinds of objects

A dataset is physically stored in two parts: a header and a data array.
as well as information that is needed to interpret the array portion
of the dataset. Essentially, it includes the name, datatype, dataspace,
and storage layout of the dataset. The name is a text string identifying
the dataset. The datatype describes the type of the data array elements.
The dataspace defines the dimensionality of the dataset, i.e., the
size and shape of the multi-dimensional array. The dimensions of a
dataset can be either fixed or unlimited (extensible). The storage
layout specifies how the data arrays are arranged in the file.

The data array contains the values of the array elements and can be
either stored together in a contiguous file space or split into smaller
\textit{chunks} stored at any allocated location. Chunks are defined
as equally-sized multi-dimensional subarrays (blocks) of the whole
data array and each chunk is stored in a separate contiguous file
space. Extensible datasets whose dimensions can grow are required
to be stored in chunks. One dimension is increased by allocating new
chunks at the end of the file to cover the extension.

by \textit{hyperslabs}, which consist of a subarray or strided subarray
of the multi-dimensional dataset. The selection is performed in the
file dataspace for the dataset. HDF5 also supports parallel I/O. Parallel
access is supported through MPI-IO. The file and datasets are collectively
created/opened by all participating processes. Each process accesses
part of a dataset by defining its own file dataspace for that dataset.
When accessing data, the data transfer property specifies whether
each process will perform independent I/O or all processes will perform
collective I/O.

\section{Input Parameters}

To enable HDF5 output in CitcomS, all you need to do is include the
following section in your \texttt{.cfg} input file.
\begin{lyxcode}
{[}CitcomS.solver.output{]}

output\_format~=~hdf5
\end{lyxcode}
Alternatively, you can specify the option \texttt{-{}-solver.output.output\_format=hdf5}
on the command line. The resulting filenames will start with the value
of your specified \texttt{datafile} input parameter and end with \texttt{.h5}.

\subsection{Optimizing Parallel I/O}

There are several platform-dependent parameters used by the HDF5 library
and the underlying MPI-IO routines to optimize the performance of
parallel I/O. The optimal values for these parameters may vary from
file system to file system. Ideally, before compiling CitcomS, the
build procedure would configure these parameters based on your platform,
but this is not implemented currently.

In order to facilitate the process of gathering I/O performance data
from a variety of parallel file systems such as GPFS, PVFS, IBRIX
FusionFS, etc., you can specify the following parameters in the CitcomS
input file. You may use these parameters to tune the performance of
the parallel I/O on your system.

All values are assumed to be specified in bytes, unless otherwise
indicated.
\begin{enumerate}
\item MPI file hints for collective buffering file access.

\begin{enumerate}
\item \texttt{\small{cb\_block\_size}}: Target nodes will access data in
chunks of this size.
\item \texttt{\small{cb\_buffer\_size}}: Specifies the total buffer space
on each target node. Set this parameter to a multiple of \texttt{\small{cb\_block\_size}}.
\end{enumerate}
\item HDF5 file access properties.

\begin{enumerate}
\item \texttt{\small{sieve\_buf\_size}}: Maximum size of data sieve buffer.
\item \texttt{\small{output\_alignment}}: Alignment interval size in bytes.
\item \texttt{\small{output\_alignment\_threshold}}: Allocation requests
above this size will be aligned on a memory address that is a multiple
of \texttt{\small{output\_alignment}}.
\item \texttt{\small{cache\_rdcc\_nelmts}}: Maximum number of chunks that
can be stored in the raw data chunk cache.
\item \texttt{\small{cache\_rdcc\_nbytes}}: Size of raw data chunk cache
in bytes.\\

\end{enumerate}
\end{enumerate}
For more details, you can refer to the following references:
\begin{itemize}
\item \textbf{MPI-2: Extensions to the Message-Passing Interface, section
9.2.8 \url{www-unix.mcs.anl.gov/mpi/mpi-standard/mpi-report-2.0/node182.htm}}
provides a list of MPI-IO reserved file hints, also available in Section
7.2.8 of the MPI-2 reference book.
\item \textbf{HDF5 documentation \url{hdf.ncsa.uiuc.edu/HDF5/doc/UG/08_TheFile.html}},
Chapter 2, Section 7.3, contains a list of HDF5 file access properties.
\item \textbf{HDF5 User's Guide: Data Caching \url{hdf.ncsa.uiuc.edu/HDF5/doc/Caching.html}}
offers a short explanation of raw data chunk caching.
\end{itemize}
You should also refer to the documentation for your particular parallel
file system and check if there are any other MPI or HDF5 hints that
could improve your parallel I/O performance.

Finally, here is an example section that would appear in a typical
CitcomS input file:
\begin{lyxcode}
{[}CitcomS.solver.output{]}

cb\_block\_size~=~1048576~~~~~~~~~~~~~~\#~1~MiB~\\
cb\_buffer\_size~=~4194304~~~~~~~~~~~~~\#~4~MiB~\\
sieve\_buf\_size~=~1048576~~~~~~~~~~~~~\#~1~MiB~\\
output\_alignment~=~262144~~~~~~~~~~~~\#~256~KiB~\\
output\_alignment\_threshold~=~524288~~\#~512~KiB~\\
cache\_rdcc\_nelmts~=~521~\\
cache\_rdcc\_nbytes~=~1048576
\end{lyxcode}

\section{\label{sec:Data-Layout}Data Layout}

Time independent data (e.g., node coordinates, grid connectivity)
are saved in a file (for example, \texttt{\small{test-case.h5}}) at
the first output stage. Subsequently, each output stage will save
time dependent data (e.g., velocity, and temperature) in a separate
file (for example, \texttt{\small{test-case.100.h5}} contains data
of time step 100). Most of the output data from CitcomS is specified
at the nodes of a logically Cartesian grid and is therefore well represented
by multi-dimensional arrays. A cap dimension is defined for addressing
each of the CitcomS caps, followed by three spatial indices $(i,j,k)$
in the case of 3D data, two spatial indices $(i,j)$ in the case of
2D surface data, or one spatial index $(k)$ in the case of 1D radial
data. An additional dimension is provided for storing the components
of vector and tensor data. These data arrays are stored in different
groups. Each group has a descriptive name. In addition, there is an
\texttt{\small{/input}} group archiving the input parameters of the
model. Two sample \texttt{\small{.h5}} files are provided in \texttt{\small{visual/samples/}}
directory.

\noindent \begin{center}
\begin{table}[H]
\noindent \begin{centering}
\begin{tabular}{|l|l|}
\hline
\textbf{Dataset} & \textbf{Shape}\tabularnewline
\hline
\hline
\texttt{/input} & \texttt{N/A}\tabularnewline
\hline
\texttt{/coord} & \texttt{(caps, nodex, nodey, nodez, 3)}\tabularnewline
\hline
\texttt{/connectivity} & \texttt{(cap\_elements, 8)}\tabularnewline
\hline
\end{tabular}
\par\end{centering}

\caption{Layout of the time-independent data file}
\end{table}

\par\end{center}

\noindent \begin{center}
\begin{table}[H]
\noindent \begin{centering}
\begin{tabular}{|l|l|}
\hline
\textbf{Dataset} & \textbf{Shape}\tabularnewline
\hline
\hline
\texttt{/velocity} & \texttt{(caps, nodex, nodey, nodez, 3)}\tabularnewline
\hline
\texttt{/temperature} & \texttt{(caps, nodex, nodey, nodez)}\tabularnewline
\hline
\texttt{/viscosity} & \texttt{(caps, nodex, nodey, nodez)}\tabularnewline
\hline
\texttt{/pressure} & \texttt{(caps, nodex, nodey, nodez)}\tabularnewline
\hline
\texttt{/stress} & \texttt{(caps, nodex, nodey, nodez, 6)}\tabularnewline
\hline
\texttt{/geoid} & \texttt{(8)} See Section \vref{sub:Geoid-Output-(test-case.geoid.10)}
for details\tabularnewline
\hline
\texttt{/surf/velocity} & \texttt{(caps, nodex, nodey, 2)}\tabularnewline
\hline
\texttt{/surf/heatflux} & \texttt{(caps, nodex, nodey)}\tabularnewline
\hline
\texttt{/surf/topography} & \texttt{(caps, nodex, nodey)}\tabularnewline
\hline
\texttt{/botm/velocity} & \texttt{(caps, nodex, nodey, 2)}\tabularnewline
\hline
\texttt{/botm/heatflux} & \texttt{(caps, nodex, nodey)}\tabularnewline
\hline
\texttt{/botm/topography} & \texttt{(caps, nodex, nodey)}\tabularnewline
\hline
\texttt{/horiz\_avg/temperature} & \texttt{(caps, nodez)}\tabularnewline
\hline
\texttt{/horiz\_avg/velocity\_xy} & \texttt{(caps, nodez)}\tabularnewline
\hline
\texttt{/horiz\_avg/velocity\_z} & \texttt{(caps, nodez)}\tabularnewline
\hline
\end{tabular}
\par\end{centering}

\caption{Layout of the time-dependent data file}
\end{table}

\par\end{center}

\section{Accessing Data}

As previously indicated, HDF5 is a self-describing binary file format.
As such we can use a variety of tools to inspect the structure of
an HDF5 file, as well as for retrieving data in any order.

\subsection{Inspecting Structures}

To quickly inspect the structure of an HDF5 file, you can use the
command \texttt{h5ls} which is included with the HDF5 software:
\begin{lyxcode}
\$~h5ls~-r~\emph{file}.h5 \end{lyxcode} \subsection{Converting to ASCII Files} You can convert the HDF5 files to the ASCII combined capfiles described in Appendix \ref{sec:ASCII-Output} for specific time steps by using the command included in CitcomS: \begin{lyxcode} \$~h5tocap~modelname~step1~{[}step2~{[}...{]}~{]}
\end{lyxcode}
You can also convert the HDF5 files to the velo files described in
Appendix \vref{sub:Velocity-and-Temperature} for restart purposes
by using the command included in CitcomS:
\begin{lyxcode}
\$~h5tovelo~modelname~step \end{lyxcode} \subsection{Accessing Data in Python} The small Python script \texttt{h5tocap.py} provides a good example of using the PyTables extension module to access the data contained in the CitcomS HDF5 files. Using PyTables, datasets can be retrieved from disk as NumPy arrays. The retrieval avoids unnecessary copying of data by using hyperslabs, which take advantage of Python's powerful array slice-indexing. For example, obtaining the node coordinates, temperature, and topography values over the entire surface of the sphere for time step 100 can be done easily with the following code snippet: \begin{lyxcode} import~tables~\\ ~\\ h5file~=~tables.openFile('samples/cookbook1.h5',~'r') surface\_coords~=~h5file.root.coord{[}0:12,:,:,-1,:{]}~\\ ~\\ data100~=~tables.openFile('samples/cookbook1.100.h5',~'r') surface\_temperature~=~data100.root.temperature{[}0:12,:,:,-1{]} surface\_topography~=~data100.root.surf.topography{[}0:12,:,:{]} \end{lyxcode} In this case, the slice \texttt{0:12} refers to all caps explicitly, while the empty slice \texttt{:}'' refers to the entire extent of the corresponding dimension. The values of \texttt{-1} above refer to the last$z$-index, which corresponds to the location of the surface nodes on each of the caps. Finally, note how both HDF5 datasets and groups are conveniently accessible as Python attributes on the PyTables file object. For more details, refer to the documentation for PyTables \url{www.pytables.org} and NumPy \url{numpy.scipy.org}. \subsection{Accessing Data Using HDFView} NCSA HDFView is a visual tool for accessing HDF files. You can use it for viewing the internal file hierarchy in a tree structure, creating new files, adding or deleting groups and datasets, and modifying existing datasets. HDFView is capable of displaying 2D slices of multi-dimensional datasets, with navigation arrow buttons that enable you to range over the entire extent of a third dimension. An example of using HDFView is found in Figure \ref{fig:A-screenshot-of} which uses data from Cookbook 1 in Section \vref{sec:Cookbook-1:-Global}. \noindent \begin{center} \begin{figure}[H] \noindent \begin{centering} \includegraphics[scale=0.53]{graphics/hdfview.png} \par\end{centering} \caption{\label{fig:A-screenshot-of}A screenshot of HDFView. The left panel shows the hierarchy of the groups and datasets. The right panel shows a 2D slice of a dataset. The bottom panel shows the metadata associated with the selected group or dataset.} \end{figure} \par\end{center} \chapter{\label{cha:Postprocessing-and-Graphics}Postprocessing and Graphics} \section{Introduction} Once you have run CitcomS, you would have a series of output files (potentially spread throughout the file systems of your Beowulf cluster or in a set of directories on your parallel file system). You now have to retrieve and combine the data for the time step (or age) of interest. To visualize your results, it is recommended that you use the open-source Open Visualization Data Explorer, better known as OpenDX. The software is available from the OpenDX website \url{www.opendx.org}. If you are using Linux, OpenDX is usually available as package \texttt{dx} or \texttt{opendx} in your distribution. If you are using Mac OS X, OpenDX is available via Fink \url{fink.sourceforge.net}. We provide experimental scripts for a 3D visualization program called MayaVi2; see Section \vref{sec:Using-MayaVi-for}. Scripts using GMT commands to plot 2D cross sections of temperature field are also provided; see Section \vref{sec:Using-GMT-Commands}. \section{Postprocessing on a Beowulf Cluster} Generally, the results from your CitcomS run will be distributed on disks attached to individual nodes of your Beowulf cluster. The output files are written in each node under the directory that you specified as the \texttt{datadir} property in the input file. To examine those files, log in to a node and change directories to the one you specified with a prefix. For example, if you set \texttt{datadir=/scratch/}\texttt{\emph{username}} and \texttt{datafile=test-case} in your input file, then your output files will be written to \texttt{/scratch/}\texttt{\emph{username}} and will have the prefix \texttt{test-case}. If you select HDF5 for the output format, the output files will have a \texttt{.h5} suffix. You won't need to postprocess the \texttt{.h5} files; you can visualize the results using OpenDX. See Section \vref{sec:Using-OpenDX-for-HDF5} for more details. If you select compressed ASCII \texttt{(ascii-gz)} for the output format, there is no post-processing script for it. You need to decompress and rename the output files before you can use the post-processing script described below. If you select ASCII for the output format, you will have many output files. An example of a filename for the velocity output is \texttt{test-case.velo.2.10} where \texttt{test-case} is the model prefix, \texttt{velo} means that this is a velocity data file, \texttt{2} corresponds to the processor number (i.e., it is output by the 2nd processor), and \texttt{10} corresponds to the time step. If your run used the time-dependent velocity boundary conditions (\texttt{solver.param.file\_vbcs=on}), the log file will also have a \texttt{current age} line that lists the time-step numbers and their corresponding times. To choose an age to export for postprocessing, you have to determine which time step corresponds to the age that interests you by looking at the log file. When you execute a CitcomS run, your input parameters will be saved in a file \texttt{pidxxxxx.cfg} where \texttt{xxxxx} is usually a five-digit number for the process ID. This pidfile contains most of the input parameters, which can be useful for archiving and postprocessing. The ASCII output files of CitcomS need to be postprocessed before you can perform the visualization. The script \texttt{autocombine.py} can postprocess (retrieve and combine) CitcomS output; it will retrieve CitcomS data to the current directory and combine the output into a few files. Using \texttt{autocombine.py}, retrieve and combine data of time-step 10: \begin{lyxcode} \$~autocombine.py~mpirun.nodes~pid12345.cfg~10~
\end{lyxcode}
This reads the MPI machinefile (\texttt{mpirun.nodes}) and the CitcomS
pidfile (\texttt{pid12345.cfg}), then calls other scripts to do the
actual job. The general usage for \texttt{autocombine.py} is:
\begin{lyxcode}
\$~autocombine.py~{[}machinefile{]}~{[}pidfile{]}~\textbackslash{}~\\ {[}step1{]}~{[}step2~or~more~...{]}~ \end{lyxcode} If your Beowulf cluster uses \texttt{ssh} (rather than \texttt{rsh}) to access the computation nodes, you must manually edit \texttt{visual/batchpaste.sh} and replace rsh' with ssh' in the script. Once \texttt{autocombine.py} has run, you will have 2 files (or 24 files for the full spherical version of CitcomS) formatted as follows: \begin{lyxcode} test-case.cap00.10~\\ test-case.cap00.10.general \end{lyxcode} The former file is the data file containing simulation results and is referred to as the capfile''; its format can be found in Appendix \ref{cha:Appendix-C:-CitComS,}. The latter file is the OpenDX header for the data. If you have enabled some \texttt{optional\_output}, and the optional fields are node-based (e.g, \texttt{comp\_nd}, \texttt{stress}, or \texttt{pressure}), these fields will be combined in separated files, with filenames: \begin{lyxcode} test-case.opt00.10~\\ test-case.opt00.10.general \end{lyxcode} The former file is the data file containing optional fields and is referred to as the optfile.'' Its format is column-based ASCII data. The latter file is the OpenDX header for the data. The header file contains information on the ordering of the optional fields. \section{Postprocessing in a Non-Cluster Environment} If you run CitcomS in a non-cluster environment or all of your data can be accessed from the local machine, you can still use \texttt{autocombine.py} to combine the data. In this case, the \texttt{machinefile} is not needed and can be replaced by \texttt{localhost}, such as: \begin{lyxcode} \$~autocombine.py~localhost~{[}pidfile{]}~\textbackslash{}~\\
{[}step1{]}~{[}step2~or~more~...{]}~
\end{lyxcode}

\section{Using OpenDX for Regional Sphere Visualization}

OpenDX modules designed for CitcomS can be found in the source directory
called \texttt{visual}. The optional \texttt{make install} command
installs the OpenDX modules under \texttt{PREFIX/share/CitcomS/visual}
(where \texttt{PREFIX} defaults to \texttt{/usr/local}).

In this example, you will use \texttt{visRegional.net} to visualize
the results of regional CitcomS.
\begin{enumerate}
\item Launch OpenDX by typing:

\begin{lyxcode}
\$~dx~ \end{lyxcode} You will see an OpenDX \textsf{Data Explorer} window. \item Click \textsf{Edit Visual Programs} and select \texttt{visRegional.net} from the file selector. \item You will see a \textsf{Visual Program Editor} window. \item Select the \textsf{import} tab in the \textsf{Visual Program Editor}'s main window and double-click on the \textsf{FileSelector} block, which will open a \textsf{Control Panel} window. \item In the \textsf{CitcomSImport filename} input box, select the header file of your postprocessed data, e.g., \begin{lyxcode} samples/regtest.cap00.100.general \end{lyxcode} \item In the pull-down menu, select \textsf{Execute$\triangleright$Execute on Change}. \item A new window will appear with the image of the model. \item If you want to zoom, rotate, change projection, and otherwise manipulate the view of the model, experiment with the menu \textsf{Options$\triangleright$View Control}. \end{enumerate} \noindent \begin{center} \begin{figure}[H] \begin{centering} \includegraphics{graphics/c_fig5.pdf} \par\end{centering} \caption{Regional Model Visualized with OpenDX. A snapshot of an upwelling (blue isosurface) with a slice of the temperature field (bisecting plane).} \end{figure} \par\end{center} Additional options in the control panel window for the regional model include: \begin{itemize} \item \textsf{\textbf{CitcomSImport reduced}} can increase or reduce the resolution of the velocity vectors and grids. Reducing resolution is often necessary for visualizing large dataset interactively. Note that the resolution of other fields (e.g., temperature and viscosity) is not reduced. \item \textsf{\textbf{Core radius}} is the size of the orange sphere which represents the core-mantle boundary. \item \textsf{\textbf{Isosurface value}} is the temperature isosurfaces. You can change the values of the isosurfaces, including adding additional isosurfaces or deleting existing ones. \item \textsf{\textbf{Slab cut dimension}} is the direction of the slab (an OpenDX term for cross-section). \item \textsf{\textbf{Slab cut position}} is the position of the slab. \end{itemize} You can change any of the parameters, visualization, or set-up by going back to the main window and clicking on each tab. If you click on each block, you will be able to change the settings for that function. \section{Using OpenDX for Full Sphere Visualization} \begin{enumerate} \item After launching OpenDX, you will see an OpenDX \textsf{Data Explorer} window. \item Click \textsf{Edit Visual Programs} and select \texttt{visFull.net} from the file selector. \item You will see a \textsf{Visual Program Editor} window. \item Select the \textsf{import} tab in the \textsf{Visual Program Editor}'s main window and double-click on the \textsf{FileSelector} block, which will open a \textsf{Control Panel} window. \item In the \textsf{Format String of CitcomSFullImport} input box, select one of the 12 header files of your postprocessed data, e.g., \begin{lyxcode} samples/fulltest.cap00.100.general \end{lyxcode} \item The results of a full spherical CitcomS consist of 12 cap files. In order to import the 12 cap files at the same time, edit the filename with the cap number replaced by printf-styled format string \texttt{\%02d}, e.g., \begin{lyxcode} samples/fulltest.cap\%02d.100.general \end{lyxcode} \item In the pull-down menu, select \textsf{Execute$\triangleright$Execute on Change}. \item A new window will appear with the image of the model. \item If you want to zoom, rotate, change projection, and otherwise manipulate the view of the model, experiment with the menu \textsf{Options$\triangleright$View Control}. \end{enumerate} Additional options in the control panel window for the spherical model include: \begin{itemize} \item \textsf{\textbf{CitcomSFullImport reduced}} can increase or reduce the resolution of the velocity vectors and grids. Reducing resolution is often necessary for visualizing large dataset interactively. Note that the resolution of other fields (e.g., temperature and viscosity) is not reduced. \item \textsf{\textbf{Core radius}} is the size of the orange sphere which represents the core-mantle boundary. \item \textsf{\textbf{Isosurface value}} is the temperature isosurfaces. You can change the values of the isosurfaces, including adding additional isosurfaces or deleting existing ones. \item \textsf{\textbf{Latitude of normal axis}} and \textsf{\textbf{Longitude of normal axis}} are the directions of the normal axis to the cross-section plane. \end{itemize} \section{\label{sec:Using-OpenDX-for-HDF5}Using OpenDX for HDF5 Visualization} If you use the HDF5 output format (\texttt{solver.output.output\_format=hdf5}), you can directly visualize the data without postprocessing. First, you need to install and set up the OpenDXutils package (see Section \vref{sub:OpenDXutils}). Then, open either \texttt{visRegional.net} or \texttt{visFull.net} in OpenDX. \begin{enumerate} \item In the pull-down menu, select \textsf{File$\triangleright$Load Macro,} and select the file \textsf{CitcomSImportHDF5.net} to load it. \item In the \textsf{Tools} panel of the main window, select the \textsf{CitcomSImportHDF5} module in the \textsf{Macros} category (highlighted in Figure \ref{fig:Import-CitComS.py-HDF5}). \item Place the module in the work space and rewire the network as shown in Figure \ref{fig:Import-CitComS.py-HDF5}. \begin{figure}[t] \noindent \begin{centering} \includegraphics[width=0.75\paperwidth,keepaspectratio]{graphics/CitcomSImportHDF5.png} \par\end{centering} \caption{\label{fig:Import-CitComS.py-HDF5} How to import CitcomS HDF5 data. The \textsf{CitcomSImportHDF5} module is highlighted in the Tools panel.} \end{figure} \item There are four input tabs in the \textsf{CitcomSImportHDF5} module. \begin{enumerate} \item The first tab (connected to the left \textsf{FileSelector} module) specifies the HDF5 file containing time-independent information (e.g., \texttt{samples/cookbook1.h5}). \item The second tab (connected to the right \textsf{FileSelector} module) specifies the HDF5 file containing time-dependent information (e.g., \textsf{samples}\texttt{/cookbook1.100.h5}). \item The third tab (connected to the \textsf{Integer} module) specifies the resolution reduction factor. \item The fourth tab (unconnected, default to 0) specifies which cap(s) to import. \begin{enumerate} \item To visualize a regional model, leave this tab unchanged. \item To visualize a full spherical model, double-click this module and change the value of this tab to \texttt{\{0,1,2,3,4,5,6,7,8,9,10,11\}}. \item You can also specify a subset of caps to visualize. \end{enumerate} \end{enumerate} \item After the data are imported, the visualization can be manipulated as described in previous sections. \end{enumerate} \section{\label{sec:Using-MayaVi-for}Using MayaVi for Visualization} This distribution also comes with scripts for visualizing CitcomS using MayaVi2. MayaVi2 is the successor of MayaVi for 2D/3D scientific data visualization. It is an interactive program allowing elaborate plots of scientific data. MayaVi2 relies on the Visualization Toolkit (VTK) and allows easy scripting in Python. To install MayaVi2, you should follow the instructions on the MayaVi2 website \url{https://svn.enthought.com/enthought/wiki/GrabbingAndBuilding} and on the SciPy website \url{www.scipy.org/Cookbook/MayaVi/Installation}. Besides MayaVi2, before using the script you will need to install PyTables (see Section \vref{sub:PyTables}) and PyVTK; see PyVTK website \url{cens.ioc.ee/projects/pyvtk} for instructions. The script acts as an extension module of MayaVi2. The installation of the script is still preliminary. You will need to include \texttt{PREFIX/share/CitcomS/visual/Mayavi2/citcoms\_plugins/} in your \texttt{PYTHONPATH} environment variable. You also need to create a directory called \texttt{\textasciitilde{}/.mayavi2} and copy into it the following: \begin{lyxcode} PREFIX/share/CitcomS/visual/Mayavi2/mayavi\_custom\_ui.py \end{lyxcode} To run the script, type this command: \begin{lyxcode} \$~mayavi2\_citcoms\_display.py~\textit{file.step.h5}
\end{lyxcode}
This will launch MayaVi2 and load the data from \texttt{file.step.h5}.
An example of using this script is found in Figure \ref{fig:Example-MayaVi2-visualization}.
to \textit{The Visualization Toolkit User's Guide} (ISBN 1-930934-13-0)
or see the VTK website \url{www.vtk.org}.

\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[scale=0.29]{graphics/mayavi2.png}
\par\end{centering}

\caption{\label{fig:Example-MayaVi2-visualization}Example MayaVi2 visualization
of \texttt{\small{cookbook1.100.h5}}. Here we display the velocity
field as vectors, as well a slice of the temperature field. Note the
visualization pipeline in the tree view on the left panel.}
\end{figure}

\par\end{center}

\section{\label{sec:Using-GMT-Commands}Using GMT Commands for Visualization}

The Generic Mapping Tools (GMT) is a collection of command-line tools
for manipulating geographic and Cartesian data that produces PostScript
(PS) or Encapsulated PostScript File (EPS) illustrations. GMT is widely
used in the geophysics community. Two scripts \texttt{plot\_layer.py}
and \texttt{plot\_annulus.py}, which can plot the temperature field
in horizontal and radial cross sections, respectively, are provided.
These scripts use GMT commands to generate EPS images. GMT is very
customizable. Users might wish to customize the scripts for their
need.

The usage of \texttt{plot\_layer.py} is:
\begin{lyxcode}
\$~plot\_layer.py~modelname~caps~step~layer \end{lyxcode} The script will look for the capfile(s) \texttt{modelname.cap00.step}, if \texttt{caps} is 1, or \texttt{modelname.cap00.step}, \texttt{modelname.cap01.step}, ... through \texttt{modelname.cap11.step}, if \texttt{caps} is 12. A slice of the capfile at the specified radial layer is then plotted in Mercator projection, if \texttt{caps} is 1, or in Hammer-Aitoff projection, if \texttt{caps} is 12. \begin{figure}[H] \noindent \begin{centering} \includegraphics[scale=0.75]{graphics/regtest.jpg} \par\end{centering} \caption{Temperature field generated by \texttt{plot\_layer.py}. This image is the mid-layer (layer 5) of the sample data \texttt{visual/samples/regtest.cap00.100}. The resolution of the image is poor because the data has only 9$\times$9$\times$9 nodes.} \end{figure} The usage of \texttt{plot\_annulus.py} is: \begin{lyxcode} \$~plot\_annulus.py~modelname~step
\end{lyxcode}
The script will ask you a few questions on how to specify the great
circle path. There are three options: (1) a starting point and an
azimuth, (2) a starting point and another point on the great circle,
and (3) a starting point and a rotation pole. The great circle path
is used to construct the radial cross section.

\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[scale=0.9]{graphics/fulltest.jpg}
\par\end{centering}

\caption{Temperature fields generated by \texttt{plot\_annulus.py}. The top
panel is the horizontal cross section at the mid layer in Hammer-Aitoff
projection. The thick line marks the surface track of the radial cross
section. The bottom panel is the radial cross section. The cross section
starts at $0^{\circ}$ latitude and $0^{\circ}$ longitude and has
an azimuth of ${45}^{\circ}$. The image is from the sample data \texttt{visual/samples/fulltest.cap{*}.100}.
There is always a gap on the left (starting point) of the cross section.
Users familiar with GMT are welcome to contribute a fix to the script.}
\end{figure}

\chapter{Cookbooks\label{cha:Cookbooks}}

\section{Introduction}

These cookbook examples are meant to serve as a guide to some of the
types of problems CitcomS can solve. Each of these cookbook problems
will only run with CitcomS if Pyre is enabled. Cookbook examples range
mantle convection problems. These cookbook examples are distributed
with the package under the \texttt{examples} directory. However, you
might need to edit these example scripts slightly to launch the job
on your cluster (see Section \vref{sec:Multiprocessor-Example} for

Cookbooks 1 to 4 introduce the basic parameters and are suitable for
all users. Cookbook 5 introduces time-dependent velocity boundary
conditions, reading in initial temperature field from files and tuning
of the advection solver. Cookbook 6 introduces the pseudo-free-surface
formulation for studying short wavelength dynamics topography. Cookbook
7 introduces thermo-chemical convection problem and tuning of the
Conjugate Gradient velocity solver. Cookbook 8 introduces compressible
convection problem, checkpointing/restarting, geoid and tuning of
the Multigrid and Uzawa velocity solver. Cookbook 9 introduces embedding
one CitcomS domain within another CitcomS domain, and further tuning
of the advection solver. Cookbook 10 introduces how to convert the
convection field to seismic velocities and how to generate synthetic
seismograms using SPECFEM3D\_GLOBE at the CIG Seismology Web Portal.\newpage{}

\section{\label{sec:Cookbook-1:-Global}Cookbook 1: Global Model}

\subsection{Problem}

This example solves for thermal convection within a full spherical
shell domain. The full spherical version of CitcomS is designed to
run on a cluster that decomposes the spherical shell into 12 equal
caps'' and then distributes the calculation for caps onto separate
processors. To run CitcomS with the full solver parameter set, it
is recommended that you have a minimum of 12 processors available

\begin{figure}[h]
\begin{centering}
\includegraphics{graphics/cookbook1.pdf}
\par\end{centering}

\caption{Global Model Caps.'' Left (A): Three-dimensional perspective image
showing seven of the 12 spherical caps used in a full CitcomS run.
Right (B): The temperature field at 1081 km depth from a Cookbook
1 run.}
\end{figure}

\subsection{Solution}

You will use \texttt{cookbook1.cfg}. The first set of parameters specifies
that you are going to use the full spherical version of the solver.
The default is to use the regional spherical version, so this must
be set.
\begin{lyxcode}
solver~=~full
\end{lyxcode}
The second set of parameters specifies the number of time steps (100),
how often full outputs of the computation are created (25), and the
prefix of output filenames (cookbook1).
\begin{lyxcode}
steps~=~100~\\
monitoringFrequency~=~25~\\
datafile~=~cookbook1
\end{lyxcode}
The \texttt{solver.ic} facility controls the temperature field for
the initial conditions. The last set of parameters includes the number
of perturbations to the initial temperature (1), the number of nodal
lines of the perturbation in the longitudinal direction, e.g., the
spherical harmonic order (3), the number of nodal lines in the latitudinal
direction, e.g., the spherical harmonic degree (2), which layer the
perturbation is on (5), and the amplitude of the perturbation (0.05).
Note that although the number of perturbations is assigned here as
\texttt{num\_perturbations=1}, that is actually the default value.
\begin{lyxcode}
num\_perturbations~=~1~\\
perturbl~=~3~\\
perturbm~=~2~\\
perturblayer~=~5~\\
perturbmag~=~0.05
\end{lyxcode}
This example is executed by typing
\begin{lyxcode}
\$~citcoms~cookbook1.cfg \end{lyxcode} \subsubsection{Example: Global Model, \texttt{cookbook1.cfg}} \begin{lyxcode} {[}CitcomS{]}~\\ solver~=~full~\\ steps~=~100~~~~~~~~~~~~~~~~~;~number~of~time~steps~\\ ~\\ {[}CitcomS.controller{]}~\\ monitoringFrequency~=~25~~~~;~how~often~outputs~are~created~\\ ~\\ {[}CitcomS.solver{]}~\\ datafile~=~cookbook1~~~~~~~~;~prefix~of~output~filenames~\\ ~\\ {[}CitcomS.solver.ic{]}~\\ num\_perturbations~=~1~\\ perturbl~=~3~\\ perturbm~=~2~\\ perturblayer~=~5~\\ perturbmag~=~0.05~ \end{lyxcode} Once you have run the model, you can visualize the results using OpenDX, described in the previous chapter. When you invoke Edit Visual Program,'' select \texttt{visFull.net}. \begin{figure}[H] \begin{centering} \includegraphics[scale=0.8]{graphics/cookbook1.jpg} \par\end{centering} \caption{Cookbook 1: Global Model. This image, created by OpenDX, depicts a slice through a spherical model of convection, with warmer colors indicating upwelling and the cooler colors showing downwelling.} \end{figure} \subsection{Discussion} You have generated a simple example of the full CitcomS model, with minimal parameter alterations. With a default Rayleigh number of$10^{5}$and perturbation on the initial temperature field which has a degree-3 and an order-2 pattern, after 100 time steps, the convection pattern remains relatively steady. As a side note, it is not required that 12 processors, with one spherical cap per processor, be used. As an end-member possibility, for example, 12 different jobs could be run on a single computer (\texttt{n001} in this example) by invoking: \begin{lyxcode} \$~citcoms~cookbook1.cfg~-{}-launcher.nodegen=\textquotedbl{}n\%03d\textquotedbl{}~\textbackslash{}~\\
-{}-launcher.nodelist={[}1,1,1,1,1,1,1,1,1,1,1,1{]}
\end{lyxcode}
This is not particularly efficient, but it does illustrate the flexibility
of both MPI and Pyre.

\newpage{}

\section{Cookbook 2: Domain Size and Velocity Boundary Conditions}

\subsection{Problem}

This example solves for thermal convection with velocity boundary
conditions imposed on the top surface within a given region of a sphere.
This requires using the regional version of CitcomS.

\subsection{Solution}

You will use \texttt{cookbook2.cfg}. The parameters specify the number
of time steps (60), the prefix of the output filenames (\texttt{cookbook2}),
and how often outputs will be created (30).
\begin{lyxcode}
steps~=~60~\\
monitoringFrequency~=~30~\\
datafile~=~cookbook2
\end{lyxcode}
The \texttt{solver.mesher} facility has several properties involved
in the generation of the computational mesh. Continue to use the default
values for the physical portion of the domain in which you are interested.
However, try modifying the layout of the mesh as shown. The parameters
specify two processors per cap in the $x$- and $y$-direction, one
processor per cap in the $z$-direction (by default), 17 nodes per
cap in the $x$- and $y$-direction, 9 nodes in the $z$-direction,
and the minimum/maximum extents in the colatitudal, longitudal, and
The unit for \texttt{theta} and \texttt{fi} is radian.
\begin{lyxcode}
nprocx~=~~2~\\
nprocy~=~~2~\\
nodex~~=~17~\\
nodey~~=~17~\\
nodez~~=~~9~\\
~\\
theta\_min~=~0.7854~\\
theta\_max~=~1.5708~\\
fi\_min~=~0.0~\\
fi\_max~=~0.7854~\\
\end{lyxcode}
The \texttt{solver.bc} facility allows you to impose a uniform velocity
across the top surface. You have a velocity which is purely in the
colatitude direction with a non-dimensional velocity of 100 (see Equation
\vref{eq:u dim} for how to dimensionalize it).
\begin{lyxcode}
topvbc~~~~=~~~1~\\
topvbxval~=~100~\\
topvbyval~=~~~0
\end{lyxcode}
In addition, the initial temperature field has a linear conductive
profile. The amplitude of initial temperature perturbation is set
to zero using the \texttt{solver.ic} facility.
\begin{lyxcode}
num\_perturbations~=~1~\\
perturbmag~~~~~~~~=~0.0~
\end{lyxcode}

\subsubsection{Example: Velocity Boundary Conditions, \texttt{cookbook2.cfg}}
\begin{lyxcode}
{[}CitcomS{]}~\\
steps~=~60~~~~~~~~~~~~~~~~~~;~number~of~time~steps~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~30~~~~;~how~often~outputs~are~created~\\
~\\
{[}CitcomS.solver{]}~\\
datafile~=~cookbook2~~~~~~~~;~prefix~of~output~filenames~\\
~\\
\#~Modify~the~layout~of~the~mesh.~\\
{[}CitcomS.solver.mesher{]}~\\
nprocx~=~~2~\\
nprocy~=~~2~\\
nodex~~=~17~\\
nodey~~=~17~\\
nodez~~=~~9~\\
~\\
theta\_min~=~0.7854~\\
theta\_max~=~1.5708~\\
fi\_min~=~0.0~\\
fi\_max~=~0.7854~\\
~\\
\#~Impose~a~uniform~velocity~across~the~top~surface.~\\
{[}CitcomS.solver.bc{]}~\\
topvbc~~~~=~~~1~\\
topvbxval~=~100~\\
topvbyval~=~~~0~\\
~\\
{[}CitcomS.solver.ic{]}~\\
num\_perturbations~=~1~\\
perturbmag~~~~~~~~=~0.0~
\end{lyxcode}

\subsection{Discussion}

Using OpenDX to visualize the results (Figure \ref{fig:Cookbook-2:-Velocity}),
this model allows you to create a plate-driven convection in which
there is a thermal upwelling on one wall, a thermal downwelling on
another, and uniform horizontal velocity across the top. The downwelling
is not exactly subduction because the sidewalls have reflecting boundary
conditions. This means that there is a symmetrical downwelling in
a vertical domain on the other side.

\begin{figure}[h]
\begin{centering}
\includegraphics{graphics/cookbook2dot2.pdf}
\par\end{centering}

\caption{\label{fig:Cookbook-2:-Velocity}Cookbook 2: Velocity Boundary Conditions.
This model, visualized with OpenDX, highlights a region of the sphere
and the heated upwellings (warm colors), downwellings (cool colors),
and the velocities (yellow arrows).}
\end{figure}

\newpage{}

\section{Cookbook 3: Temperature-Dependent Viscosity}

\subsection{Problem}

A common problem in geophysics is the exploration of natural convection
in the presence of variable viscosity, including temperature-dependent
or stress-dependent viscosity.

\subsection{Solution}

You will use \texttt{cookbook3.cfg}. The parameters specify the number
of time steps (200), how often outputs will be created (25), the prefix
of the output filenames (\texttt{cookbook3}), and the Rayleigh number
($10^{6}$).
\begin{lyxcode}
steps~=~200~\\
monitoringFrequency~=~25~\\
datafile~=~cookbook3~\\
rayleigh~=~1e6
\end{lyxcode}
The \texttt{solver.visc} facility assigns the viscosities. The parameters
specify whether the viscosity should be updated at every time step
(\texttt{on}), the number of viscous layers (4), the viscosity of
each layer (1,1,1,1), whether the viscosity is temperature dependent
(\texttt{on}), the type of viscosity law to use (4), the activation
energy of each layer (0.2,0.2,0.2,0.2), the temperature offset of
each layer (0,0,0,0), the activation volume of each layer (0,0,0,0),
whether to apply the minimum cutoff (\texttt{on}), the value of the
minimum cutoff (1.0), whether to apply the maximum cutoff (\texttt{on}),
and the value of the maximum cutoff (100.0).
\begin{lyxcode}
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
rheol~=~4~\\
viscE~=~0.2,0.2,0.2,0.2~\\
viscT~=~0,0,0,0~\\
viscZ~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0
\end{lyxcode}
The range of the layers are determined by the following parameters:
layer 1 extends from the top surface to a depth of \texttt{solver.const.z\_lith}
(default to 0.014, or 90 km dimensionally); layer 2 is below layer
1 and extends to a depth of \texttt{solver.const.z\_410} (default
to 0.06435, or 410 km dimensionally); layer 3 is below layer 2 and
extends to a depth of \texttt{solver.const.z\_lmantle} (default to
0.105, or 670 km dimensionally); layer 4 is below layer 3 and extends
to the bottom. These depth parameters also control the depth of various
phase transitions.

There are several viscosity laws coded in the program, which can be
selected by the \texttt{rheol} parameter. All available viscosity
laws can be found in Appendix \ref{sub:Viscosity-input}. For \texttt{rheol=4},
the temperature dependence of the viscosity of each layer is determined
by:

$$visc=visc0\times exp\left(\frac{viscE+viscZ\times(1-r)}{T^{*}+viscT}-\frac{viscE+viscZ\times(1-r_{inner})}{1+viscT}\right)\label{eq:non-dim visc law}$$
where $T^{*}=\min(\max(T,0),1)$. This equation is non-dimensionalized.
The last term on the right-hand side normalizes the viscosity such
that $visc=visc0$ when $T=1$ and $r=r_{inner}$. Comparing Equation
\ref{eq:non-dim visc law} with the dimensional Arrhenius activation
equation:
\begin{eqnarray}
\eta & = & \eta_{0}exp\left(\frac{E_{a}+PV_{a}}{RT}\right)=\eta_{0}exp\left(\frac{E_{a}+PV_{a}}{R\Delta T\left(T'+T_{0}'\right)}\right)\nonumber \\
& = & \eta_{0}exp\left(\frac{E_{a}+\rho g\left(1-r\right)V_{a}}{R\Delta T\left(T'+T_{0}'\right)}\right)=\eta_{0}exp\left(\frac{E_{a}/R\Delta T+\left(1-r\right)\rho gV_{a}/R\Delta T}{T'+T_{0}'}\right)\label{eq:arrhenius}
\end{eqnarray}
where $E_{a}$ and $V_{a}$ are the activation energy and activation
volume, $P$ is the hydrostatic pressure, $R$ is the universal gas
constant, and other symbols are as defined in Section \ref{sec:Governing-Equations},
it is clear that
$$viscE=\frac{E_{a}}{R\Delta T}\label{eq:viscE}$$
$$viscZ=\frac{\rho gV_{a}}{R\Delta T}\label{eq:viscZ}$$
$$viscT=T_{0}\label{eq:viscT}$$

\subsubsection{Example: Temperature-Dependent Viscosity, \texttt{cookbook3.cfg}}
\begin{lyxcode}
{[}CitcomS{]}~\\
steps~=~200~~~~~~~~~~~~~~~~~;~number~of~time~steps~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~25~~~~;~how~often~outputs~are~created~\\
~\\
{[}CitcomS.solver{]}~\\
datafile~=~cookbook3~~~~~~~~;~prefix~of~output~filenames~\\
rayleigh~=~1e6~~~~~~~~~~~~~~;~Rayleigh~number~\\
~\\
\#~Modify~the~layout~of~the~mesh.~\\
{[}CitcomS.solver.mesher{]}~\\
nprocx~=~~2~\\
nprocy~=~~2~\\
nodex~~=~17~\\
nodey~~=~17~\\
nodez~~=~~9~\\
~\\
\#~Assign~the~viscosities.~\\
{[}CitcomS.solver.visc{]}~\\
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
rheol~=~4~\\
viscE~=~0.2,0.2,0.2,0.2~\\
viscT~=~0,0,0,0~\\
viscZ~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0
\end{lyxcode}
\begin{figure}[H]
\begin{centering}
\includegraphics{graphics/cookbook3.pdf}
\par\end{centering}

\caption{Cookbook 3: Temperature-Dependent Viscosity. This model, visualized
with OpenDX, highlights a region that features both the heated upwelling
(warm colors) and the even distribution of the velocities (arrows).}
\end{figure}

\newpage{}

\section{Cookbook 4: Regionally Refined Meshes }

\subsection{Problem}

Frequently for convection problems, particularly those with variable
viscosity, there are features with strong gradients in temperature
or viscosity that can be better resolved with refined meshes. For
example, for the problem just studied in Cookbook 3, temperature-dependent
viscosity, a higher resolution is required for a narrow hot upwelling
while a lower resolution is sufficient for the wider cold downwelling.

\subsection{Solution}

The parameter that controls whether a mesh is to be uniform or refined
is \texttt{solver.mesher.coor}. Set it to \texttt{1} \texttt{(solver.mesher.coor=1)}
in order to read the uneven mesh point coordinates from an input file
(specified by \texttt{coor\_file}). The format of the coordinate input
file is described in Appendix \vref{sec:Coordinate-files}.
\begin{lyxcode}
coor~=~1~\\
coor\_file~=~coor.dat
\end{lyxcode}
The computational domain is bounded in colatitude between radian 1
and 2, in latitude between radian 0 and 1, and in radius between 0.55
and 1, as determined by the content of \texttt{coor\_file}.

\subsubsection{Example: Regionally Refined Meshes, \texttt{cookbook4.cfg}}
\begin{lyxcode}
{[}CitcomS{]}~\\
steps~=~250~~~~~~~~~~~~~~~~~;~number~of~time~steps~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~10~~~~;~how~often~outputs~are~created~\\
~\\
{[}CitcomS.solver{]}~\\
rayleigh~=~1e6~~~~~~~~~~~~~~;~Rayleigh~number~\\
datafile~=~cookbook4~~~~~~~~;~prefix~of~output~filenames~\\
~\\
{[}CitcomS.solver.ic{]}~\\
num\_perturbations~=~1~\\
perturbmag~=~0.05~\\
perturblayer~=~10~\\
~\\
\#~direction.~\\
perturbl~=~1~\\
perturbm~=~0~\\
~\\
{[}CitcomS.solver.mesher{]}~\\
coor~=~1~\\
coor\_file~=~coor.dat~\\
~\\
nprocx~=~2~\\
nprocy~=~1~\\
nprocz~=~1~\\
nodex~=~33~\\
nodey~=~17~\\
nodez~=~17~\\
~\\
{[}CitcomS.solver.visc{]}~\\
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
rheol~=~4~\\
viscE~=~0.2,0.2,0.2,0.2~\\
viscT~=~0,0,0,0~\\
viscZ~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0
\end{lyxcode}
\begin{figure}[H]
\begin{centering}
\includegraphics{graphics/cookbook4.pdf}
\par\end{centering}

\caption{Cookbook 4: Regionally Refined Mesher. This model shows the temperature
(upwelling -- warm colors, downwelling -- cool colors) and the uneven
distribution of the velocities (yellow arrows). Note the refined mesh
for the left and bottom regions of the model (inset).}
\end{figure}

\subsection{Discussion}

The resulting model is like a 2D model, but extends in the longitudinal
direction. To set up this type of model, the initial temperature gradient
was perturbed only in the longitudinal direction, by setting the parameters
\texttt{perturbl=1} and \texttt{perturbm=0}. The model results show
a thin thermal upwelling on one wall and a wide thermal downwelling
on the opposite wall. Also, at the bottom of the model a thin layer
of hot material can be shown. Note the higher resolution in the narrow
regions with hot material and the lower resolution for the wide thermal
downwelling region.

\newpage{}

\section{Cookbook 5: Subduction Models with Trench Rollback}

\subsection{Problem}

A common issue to address is the problem arising when the position
of the oceanic trench is not constant in time (trench rollback) for
a subduction zone. In addition, the trench rollback speed may vary
in time, which can be caused, for example, by a rotation pole jump.

\subsection{Solution}

In order to introduce in a convection model a trench rollback you
have to prepare the velocity files. The following scenario is proposed
for this cookbook: there are two plates centered along the equator,
with two different Euler poles (see Figure \ref{fig:Cookbook-5:-Left}).
Initially, both plates have Euler poles situated far apart, so the
velocities are approximately constant (about 5 cm/yr for the left
plate and about 2 cm/yr for the right plate). After 30 Ma, a pole
jump occurred for the right slab only, so the Euler pole moved closer.

Since we impose a top velocity boundary in our model, the following
parameters should be turned on:
\begin{lyxcode}
{[}CitcomS.solver.bc{]}~\\
topvbc~=~1~\\
~\\
{[}CitcomS.solver.param{]}~\\
file\_vbcs~=~on
\end{lyxcode}
The following two lines specify the location of the velocity files
and the starting age for the model:
\begin{lyxcode}
start\_age~=~55~\\
vel\_bound\_file~=~./velocity/bvel.dat~
\end{lyxcode}
Since the starting age is set to 55 Ma, there will be 57 velocity
files, one for each Ma (\texttt{vel.dat0, vel.dat1,} \texttt{...vel.dat55,
vel.dat56}). The format of the velocity file is described in Appendix
\vref{sec:Velocity-boundary-condition}. Note that the velocity file
uses a unit of cm/yr.

When \texttt{tic\_method=-1} is set, the initial temperature field
is read from setup files. The files have the same naming scheme and
format as the CitcomS velo output, as described in Appendix \vref{sub:Velocity-and-Temperature}.
The location and prefix of the file is specified in \texttt{datadir\_old}
and \texttt{datafile\_old}, and the time step of the velo file in
\texttt{solition\_cycles\_init}. In this example, a processor of MPI
rank \texttt{n} will read file \texttt{./ic/cookbook5.velo.n.0}
\begin{lyxcode}
{[}CitcomS.solver{]}~\\
datafile\_old~=~cookbook5~\\
~\\
{[}CitcomS.solver.ic{]}~\\
tic\_method~=~-1~~\\
solution\_cycles\_init~=~0
\end{lyxcode}
The descending slab will cause some numerical instability in the temperature
solver. To reduce the numerical instability, you will use a smaller
time-step size (0.75 of the maximum time-step size for a stable solution).
The change in maximum temperature will also be monitored. If the maximum
temperature increases too much (> 5\%) between time steps, the temperature
solver will rerun with half time-step size.
\begin{lyxcode}
{[}CitcomS.solver.tsolver{]}~\\
finetunedt~=~0.75~\\
monitor\_max\_T~=~on
\end{lyxcode}
\begin{figure}[H]
\begin{centering}
\includegraphics[scale=0.65]{graphics/cookbook5.pdf}
\par\end{centering}

\caption{\label{fig:Cookbook-5:-Left}Cookbook 5: Left (A): Initially, both
plates have Euler poles (magenta and green dots) situated far apart
from the plate, so the velocities are approximately constant (about
5 cm/yr for the left plate and about 2 cm/yr for the right plate).
Right (B): After 30 Ma a pole jump occurred for the right slab only,
so the Euler pole moved closer.}
\end{figure}

\subsubsection{Example: Subduction Models with Trench Rollback, \texttt{cookbook5.cfg}}
\begin{lyxcode}
{[}CitcomS{]}~\\
steps~=~1100~~~~~~~~~~~~~~~~;~number~of~time~steps~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~100~~~;~how~often~outputs~are~created~\\
~\\
{[}CitcomS.solver{]}~\\
datafile~=~cookbook5~\\
datafile\_old~=~cookbook5~\\
rayleigh~=~4.07e+08~\\
~\\
{[}CitcomS.solver.bc{]}~\\
topvbc~=~1~\\
~\\
{[}CitcomS.solver.param{]}~\\
file\_vbcs~=~on~\\
start\_age~=~55~\\
~\\
\#~Since~the~starting~age~is~set~to~55~Ma,~there~will~be~57~velocity~\\
\#~files,~one~for~each~Ma~(bvel.dat0,~bvel.dat1,~...~bvel.dat56).~\\
vel\_bound\_file~=~./velocity/bvel.dat~\\
~\\
{[}CitcomS.solver.ic{]}~\\
tic\_method~=~-1~\\
solution\_cycles\_init~=~0~\\
~\\
{[}CitcomS.solver.mesher{]}~\\
coor~=~1~\\
coor\_file~=~./coor.dat~\\
nprocx~=~1~\\
nprocy~=~2~\\
nprocz~=~1~\\
nodex~=~17~\\
nodey~=~65~\\
nodez~=~33~\\
~\\
{[}CitcomS.solver.tsolver{]}~\\
finetunedt~=~0.75~\\
monitor\_max\_T~=~on~\\
~\\
{[}CitcomS.solver.visc{]}~\\
num\_mat~=~4~\\
visc0~=~100,0.003,1,2~\\
TDEPV~=~on~\\
viscE~=~24,24,24,24~\\
viscT~=~0.182,0.182,0.182,0.182~\\
VMIN~=~on~\\
visc\_min~=~0.01~\\
VMAX~=~on~\\
visc\_max~=~100.0~
\end{lyxcode}
\begin{figure}[H]
\begin{centering}
\includegraphics{graphics/cookbook5dot2.pdf}
\par\end{centering}

\caption{\label{fig:Cookbook-5:-the}Cookbook 5: The pole jump after 30 Ma
produced a flat slab on one side (fore side) and a steep slab on the
other side (back side). }
\end{figure}

\subsection{Discussion}

The results for this problem are presented in Figure \ref{fig:Cookbook-5:-the}.
Since the two Euler poles are kept fixed for the first 30 Ma, the
shape of the subducting slab will be the same along the trench. At
25 Ma, the Euler pole for the right plate jumps toward the plate.
Therefore, the velocity along the trench varies from about 1 cm/yr
to about 2.5 cm/yr. This will produce in time a flat slab at one side
of the model, and a steep slab at the other side.

\newpage{}

\section{Cookbook 6: Pseudo-Free-Surface Formulation}

\subsection{Problem}

Free-slip boundary conditions are typically applied on the top surface
of mantle convection models, and the dynamic topography is obtained
by assuming that the normal stress on the top surface is instantaneously
compensated by the deformed surface. This type of boundary condition
works well for long-wavelength topography, but a free-surface formulation
becomes necessary in cases where intermediate to short wavelength
topography is of interest or the lithosphere has a very high effective
viscosity. Another situation where free-surface formulation is desired
is when two Pyre solvers, Snac and Citcom, are coupled.

The basic algorithm \cite{Zhong et al Free-surface formulation,Zhong et al Three-dimensional}
consists of four steps.
\begin{enumerate}
\item On the nodes of the top surface, compute the topography increment
by integrating normal velocity over time.
\item Calculate the normal traction on the top surface based on the accumulated
topography up to the current time step, and add it to the forcing
term in Equation \ref{eq:discrete momentum eqn}.
\item Update the velocity field with the changed forcing term.
\item If the velocity field has not converged yet, repeat steps 1 to 3.
\end{enumerate}

\subsection{Solution}

To verify that the above algorithm works, you will run two different
CitcomS models with each boundary condition (BC) (free-slip and pseudo-free-surface)
and compare the topography computed accordingly. The scripts in \texttt{cookbook6.cfg}
include the following parameters.
\begin{itemize}
\item Domain size: $45^{\circ}$$\times$ $45^{\circ}$ $\times$ $1200\mathrm{\, km}$
\item Mesh size: 61 $\times$ 61 $\times$ 25 nodes with mesh refinement
(\texttt{coor=1})
\item Boundary Conditions (BC): free-slip/free-surface for top; free-slip
for all the other surfaces
\item Initial Conditions (IC): spherical hot blob with diameter of $\frac{1}{4}$
of the radial dimension placed at the center of the domain
\end{itemize}
The pseudo-free-surface boundary condition is enabled by:
\begin{lyxcode}
{[}CitcomS.solver.bc{]}~\\
topvbc~=~2~\\
pseudo\_free\_surf~=~on
\end{lyxcode}
The time step size is fixed to $7.77\times10^{-10}$ (about 1000 yrs),
instead of determined dynamically by the velocity solution.
\begin{lyxcode}
{[}CitcomS.solver.tsolver{]}~\\
fixed\_timestep~=~7.77e-10
\end{lyxcode}

\subsection{Discussion}

The plots of the topography profile are consistent with results described
in \cite{Zhong et al Free-surface formulation}. In the graphs in
Figure \ref{fig:Cookbook-6.-Graphs}, the solid lines indicate free-slip
and the dashed lines indicate free-surface.

\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[scale=0.8]{graphics/cookbook6.pdf}
\par\end{centering}

\caption{\label{fig:Cookbook-6.-Graphs}Cookbook 6: Graphs of topography profiles. }
\end{figure}

\par\end{center}

\newpage{}

\section{Cookbook 7: Thermo-Chemical Convection}

\subsection{Problem}

This example solves for thermo-chemical convection within a full spherical
shell domain. Composition heterogeneity exists in the Earth mantle.
The density anomalies due to the composition heterogeneity, as well
as due to the thermal heterogeneity, drive the convection flow.

\subsection{Solution}

You will use \texttt{cookbook7.cfg}. Most of the parameters you have
encountered in previous cookbooks. The initial condition is the same
as in Cookbook 1, and the viscosity law is the same as in Cookbook
3.

A large number of ASCII files will be output. In order to keep the
current directory tidy, the output files can be put in a subdirectory
\texttt{output/} by specifying
\begin{lyxcode}
\end{lyxcode}
Since you are interested in the composition field, the compostion
output should be enabled and some output files disabled:
\begin{lyxcode}
output\_optional~=~tracer,comp\_nd
\end{lyxcode}
The most important parameters are in the \texttt{CitcomS.solver.tracer}
section. Turn on the tracer module:
\begin{lyxcode}
tracer~=~1
\end{lyxcode}
The number of tracers and their initial location must be set. Here,
you specify \texttt{tracer\_ic\_method=0}. The tracers will be generated
pseudo-randomly, with a total number equal to
\begin{lyxcode}
tracers\_per\_element$\times$\textrm{(total~number~of~finite~elements).}
\end{lyxcode}
If \texttt{tracer\_ic\_method=1}, all processors will read the location
(and flavors; see next paragraph) of the tracers from the same file
specified by \texttt{tracer\_file}. If \texttt{tracer\_ic\_method=2},
each processor will read the location (and flavors; see next paragraph)
of the tracers from the file specified by \texttt{datafile\_old} (in
the \texttt{CitcomS.solver} section) and \texttt{solution\_cycles\_init}
(in the \texttt{CitcomS.solver.ic} section).
\begin{lyxcode}
tracer\_ic\_method~=~0~\\
tracers\_per\_element~=~20~\\
tracer\_file~=~tracer.dat
\end{lyxcode}
Each tracer can have a flavor'' attached to it. The meaning of the
flavor depends on the application. Here, the flavor indicates the
chemical species of the tracer. There are two flavors of tracers (\texttt{tracer\_flavors=2}).
A tracer of flavor 0 is of the normal,'' or ambient, chemical composition,
while a tracer of flavor 1 is of the anomalous'' chemical composition.
Because the tracers are automatically generated, you need to specify
how to assign the flavor to each tracer. If \texttt{ic\_method\_for\_flavors=0},
tracers above \texttt{z\_interface} are assigned to flavor 0, and
tracers below \texttt{z\_interface} to flavor 1.
\begin{lyxcode}
tracer\_flavors~=~2~\\
ic\_method\_for\_flavors~=~0~\\
z\_interface~=~0.7
\end{lyxcode}
The thermo-chemical convection module is turned on by
\begin{lyxcode}
chemical\_buoyancy~=~1
\end{lyxcode}
The composition field is determined by the ratio method (\texttt{buoy\_type=1})
\cite{Tackely King tracer ratio method}. The density anomaly of the
anomalous chemical composition is defined by $B=\Delta\rho_{ch}/\left(\rho_{0}\alpha_{0}\Delta T\right)$,
the buoyancy ratio.
\begin{lyxcode}
buoy\_type~=~1~\\
buoyancy\_ratio~=~0.5
\end{lyxcode}
The code keeps track of the sum of the bulk composition. In a perfect
world, the sum of the bulk composition would not change with time.
But due to numerical issues, the sum of the bulk composition tends
to decrease with time. If \texttt{tracer\_ic\_method=2}, the code
will read in the location of the tracers, as well as the initial sum
of the bulk composition, from the previous run.

The regular grid is an auxiliary grid that helps locate the tracers
\cite{McNamara Zhong Thermochemical}. The optimal grid spacing of
the regular grid depends on the size of the CitcomS mesh. A general
rule of thumb is that the grid spacing of the regular grid should
be less than $18/(nodex\times nprocx)$. These parameters are used
only in the full spherical version.
\begin{lyxcode}
regular\_grid\_deltheta~=~1.0~\\
regular\_grid\_delphi~=~1.0
\end{lyxcode}
If you encounter any error, look at the end of the tracer log files
(\texttt{cookbook7.tracer\_log.{*}}) for error messages.

You will use the conjugate gradient solver to solve the matrix equation
$\mathbf{A}x=b$ for $x$. The conjugate gradient solver is more efficient
than the multigrid solver (\texttt{Solver=multigrid}) for small problems
(e.g., less than 17$\times$17$\times$17 nodes per processor). The
desired accuracy and maximum iterations are set by \texttt{accuracy}
and \texttt{vlowstep}. For the Uzawa algorithm, the maximum iterations
are set by \texttt{piterations}. \texttt{Both vlowstep }and \texttt{piterations}
should be large integers, but you can experiment with a larger value
for \texttt{accuracy} to speed up the computation.
\begin{lyxcode}
accuracy~=~1e-03~\\
vlowstep~=~1000~\\
piterations~=~1000
\end{lyxcode}

\subsubsection{Example: Thermo-Chemical Convection, \texttt{cookbook7.cfg}}
\begin{lyxcode}
{[}CitcomS{]}~\\
solver~=~full~\\
steps~=~15~~~~~~~~~~~~~~~~~;~number~of~time~steps~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~5~~~~;~how~often~outputs~are~created~\\
~\\
{[}CitcomS.solver{]}~\\
datafile~=~cookbook7~~~~~~~;~prefix~of~output~filenames~\\
rayleigh~=~1e7~\\
~\\
{[}CitcomS.solver.ic{]}~\\
num\_perturbations~=~1~\\
perturbl~=~3~\\
perturbm~=~2~\\
perturblayer~=~5~\\
perturbmag~=~0.05~\\
~\\
{[}CitcomS.solver.output{]}~\\
output\_optional~=~tracer,comp\_nd~\\
~\\
{[}CitcomS.solver.tracer{]}~\\
tracer~=~on~\\
tracer\_ic\_method~=~0~\\
tracers\_per\_element~=~20~\\
tracer\_file~=~tracer.dat~\\
~\\
tracer\_flavors~=~2~\\
ic\_method\_for\_flavors~=~0~\\
z\_interface~=~0.7~\\
~\\
chemical\_buoyancy~=~1~\\
buoy\_type~=~1~\\
buoyancy\_ratio~=~0.5~\\
~\\
regular\_grid\_deltheta~=~1.0~\\
regular\_grid\_delphi~=~1.0~\\
~\\
{[}CitcomS.solver.vsolver{]}~\\
accuracy~=~1e-04~\\
vlowstep~=~1000~\\
piterations~=~1000~\\
~\\
\#~Assign~the~viscosities.~\\
{[}CitcomS.solver.visc{]}~\\
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
rheol~=~4~\\
viscE~=~0.2,0.2,0.2,0.2~\\
viscT~=~0,0,0,0~\\
viscZ~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0
\end{lyxcode}
\noindent \begin{center}
\begin{figure}[H]
\noindent \begin{centering}
\includegraphics[width=0.3\paperwidth]{graphics/cookbook7.png}
\par\end{centering}

\caption{\label{fig:Cookbook-7:-The}Cookbook 7: The composition and velocity
field at the 15th step. The arrows are the velocity vectors. The composition
field is shown in an isosurface of 0.7 and in a cross section.}
\end{figure}

\par\end{center}

\subsection{Discussion}

When the model is running, it will output the progress of the run
to the screen. This line has information about the grid:
\begin{lyxcode}
Problem~has~9~x~9~x~9~nodes~per~cap,~6930~nodes~and~6144~elements~in~total
\end{lyxcode}
The following lines give the radial coordinate and reference density
of each nodes and their viscosity layer:
\begin{lyxcode}
~~~~~~1~~~~0.550000~~~~0.450000~1.000000e+00~~~~~4~\\
~~~~~~2~~~~0.606250~~~~0.393750~1.000000e+00~~~~~4~\\
~~~~~~3~~~~0.662500~~~~0.337500~1.000000e+00~~~~~4~\\
~~~~~~4~~~~0.718750~~~~0.281250~1.000000e+00~~~~~4~\\
~~~~~~5~~~~0.775000~~~~0.225000~1.000000e+00~~~~~4~\\
~~~~~~6~~~~0.831250~~~~0.168750~1.000000e+00~~~~~4~\\
~~~~~~7~~~~0.887500~~~~0.112500~1.000000e+00~~~~~4~\\
~~~~~~8~~~~0.943750~~~~0.056250~1.000000e+00~~~~~2~\\
~~~~~~9~~~~1.000000~~~~0.000000~1.000000e+00~~~~~1
\end{lyxcode}
This line gives the perturbation parameters used to construct the
initial temperature:
\begin{lyxcode}
Initial~temperature~perturbation:~~layer=5~~mag=0.05~~l=3~~m=2~
\end{lyxcode}
This line gives the magnitude the right-hand side vector in Equation
\ref{eq:discrete momentum eqn}:
\begin{lyxcode}
Momentum~equation~force~2.207120371e+03
\end{lyxcode}
The following lines give the convergence progress of the Stokes solver,
where \texttt{v} is the volume averaged norm of velocity, \texttt{p}
is the volume averaged norm of pressure, \texttt{div/v} is the volume
averaged norm of $\nabla\cdot(\bar{\rho}v)$ divided by \texttt{v},
\texttt{dv/v} is the volume averaged norm of velocity change divided
by \texttt{v}, and \texttt{dp/p} is the volume averaged norm of pressure
change divided by \texttt{p}. The convergence of the Stokes solver
is achieved if either (1) \texttt{div/v }is less than\texttt{ accuracy}
or (2) both \texttt{dv/v} and \texttt{dp/p} are less than \texttt{accuracy}
for two consecutive iterations. In this example, the second condition
is met.
\begin{lyxcode}
(000)~~~7.0~s~v=1.470010e+02~p=0.000000e+00~div/v=5.37e+00~dv/v=1.00e+00~dp/p=1.00e+00~step~0~\\
$\vdots$~\\
(090)~151.3~s~v=1.012752e+02~p=1.304499e+04~div/v=5.54e-04~dv/v=1.84e-05~dp/p=9.20e-05~step~0~\\
(091)~151.8~s~v=1.012752e+02~p=1.304512e+04~div/v=5.66e-04~dv/v=2.07e-05~dp/p=7.45e-05~step~0~
\end{lyxcode}
This line gives the rotation angle and the angular coordinates of
rotation pole of the removed net angular momentum:
\begin{lyxcode}
Angular~momentum:~rot=9.917258e-02~tr=8.924735e+01~fr=1.083630e+02
\end{lyxcode}
These lines give the heat flux across the top and bottom surfaces:
\begin{lyxcode}
surface~heat~flux=~2.227892~\\
bottom~heat~flux=~2.216542
\end{lyxcode}
The surface heat flux \texttt{$H_{surf}$} can be converted to the
Nusselt number \texttt{$Nu_{surf}$} by:

$$Nu_{surf}=H_{surf}\times(r_{outer}-r_{inner})\frac{r_{outer}}{r_{inner}}\label{eq:Nusselt number}$$

The results for this problem are presented in Figure \ref{fig:Cookbook-7:-The}.
The buoyancy ratio in this model is too low to stabilize the chemical
layer. A few thermo-chemical plumes are rising from the lower mantle,
especially the ones at the 4, 6, 9, 10, and 12 o'clock directions.
The resolution of this model is fairly low. The composition isosurface
is slightly discontinuous across the cap boundary. A model of higher
resolution will not have this kind of artifact.

\newpage{}

\subsection{Problem}

This example is a benchmark problem for compressible thermal convection.
The Stokes solver in CitcomS has been benchmarked and validated against
a semi-analytical solution. However, no analytical solution exists
for the benchmark on the energy equation solver, which is nonlinear.
The steady-state solution is usually used for the comparison with
other numerical solutions.

\subsection{Solution}

This cookbook example will run for 10,000 time steps to reach steady
state. It will use 12 processors and take 1 to 2 days to finish on
a modern computer. At every 1,000th time-step interval, a checkpoint
for the internal state of the solver is saved.
\begin{lyxcode}
checkpointFrequency~=~1000
\end{lyxcode}
If the solver is interrupted before finishing the computation, one
can resume the computation from the checkpointed state. To shorten
the computation time, a checkpoint at the 9,000th time step (\texttt{CitcomS-cookbook8-
restart.tar.gz}) is provided at CIG's CitcomS web page \url{geodynamics.org/cig/software/packages/mc/citcoms}.
(Note that the checkpoint files are produced by an x86 machine and
may not be usable by other types of machines, e.g., PowerPC.) To resume
the computation from the 9,000th time-step checkpoint, untar the package
in the \texttt{examples/Cookbook8/} directory and set these parameters:
\begin{lyxcode}
{[}CitcomS.solver.ic{]}~\\
restart~=~on~\\
solution\_cycles\_init~=~9000
\end{lyxcode}
If you restart from the checkpoint, the computation will resume from
the 9,000th time-step. Note that the \texttt{coord} files are only
output at the 0th time-step, so you will need to run the model without
restarting for one time-step to get the \texttt{coord} files.

A compressible convection model has four dimensionless numbers: the
Rayleigh number, the dissipation number, the Gruneisen parameter $\gamma=\alpha_{0}K_{S0}/\rho_{0}c_{P0}$,
and the non-dimensionalized absolute surface temperature, which are
defined in Section \ref{sec:Governing-Equations}. The Rayleigh number
of CitcomS is scaled by the radius of the Earth. If scaled to the
thickness of the mantle, the effective Rayleigh number is $7\times10{}^{3}$
($=7.68175583\times10{}^{4}\times(r_{outer}-r_{inner})^{3}$). Similarily,
the effective dissipation number is 0.275 ($=0.5\times(r_{outer}-r_{inner})$).
Under these non-dimensional numbers, the convection is of low vigor
and low compressibility.
\begin{lyxcode}
rayleigh~=~7.68175583e4~\\
dissipation\_number~=~0.5~\\
gruneisen~=~1.25~\\
surfaceT~=~0.1
\end{lyxcode}
Since we are going to use the multigrid solver, the grid size is specified
by:
\begin{lyxcode}
nodex~=~33~\\
nodey~=~33~\\
nodez~=~33~\\
levels~=~5
\end{lyxcode}
The additional parameter \texttt{levels} specifies the nested levels
of multigrid units and is subjected to the following constraint:
$$\mathrm{nodex}=1+\mathrm{nprocx}\times\mathrm{mgunitx}\times2^{levels-1}\label{eq:mgunit}$$
where \texttt{mgunitx} must be an integer. In the pure C version,
\texttt{mgunitx} is an input parameter.

If \texttt{reference\_state=}1, then constant gravity, heat capacity,
thermal expansivity, and $\rho_{r}=exp\left(\frac{D_{i}}{\gamma}(1-r)\right)$
are used as the reference state. If \texttt{reference\_state=}0, the
reference state is read from a file \texttt{refstate\_file}. See Appendix
\ref{cha:Appendix-A:-Input} for the file format.
\begin{lyxcode}
reference\_state~=~1~\\
refstate\_file~=~ref.dat
\end{lyxcode}
We are interested in the geoid, dynamic topography, and heat flux.
The dynamic topography and heat flux are computed on the surface grids.
The geoid is computed in the spherical harmonics, with a maximal degree
of \texttt{20}.
\begin{lyxcode}
output\_optional~=~geoid,surf,botm~\\
output\_ll\_max~=~20
\end{lyxcode}
The initial temperature is a conductive profile with a single spherical
harmonic perturbation. The perturbation is located at mid-depth and
is defined as:
$mag\times\sin\left(\frac{(r-r_{in})\pi}{r_{out}-r_{in}}\right)\left(\sin(m\phi)+\cos(m\phi)\right)P_{lm}(\cos\theta)$

\begin{lyxcode}
tic\_method~=~3~\\
num\_perturbations~=~1~\\
perturbl~=~3~\\
perturbm~=~2~\\
perturblayer~=~17~\\
perturbmag~=~0.01
\end{lyxcode}
You will need the output of dynamic topography and the geoid. The
dynamics topography will be computed by the Consistent Boundary Flux
(CBF) method, and the effect of self-gravitation is included in the
geoid. The maximum spherical harmonics degree for the geoid is 20.
\begin{lyxcode}
output\_optional~=~geoid,surf,botm~\\
use\_cbf\_topo~=~on~\\
self\_gravitation~=~on~\\
output\_ll\_max~=~20
\end{lyxcode}
Various parameters tune the performance of the solver. The maximum
size of each time step is determined dynamically by the Courant criterion.
To enhance the stability of the energy equation solver, you will only
use three quarters of the maximum Courant time-step size.
\begin{lyxcode}
finetunedt~=~0.75
\end{lyxcode}
You will use the multigrid solver to solve the matrix equation $\mathbf{A}x=b$
for $x$. The multigrid solver is more efficient than the conjugate
more than 17$\times$17$\times$17 nodes per processor). Several parameters
control the behavior of the multigrid solver: \texttt{mg\_cycle=1}
for the V cycle and \texttt{2} for the W cycle; \texttt{down\_heavy}
and \texttt{up\_heavy} are the number of smoothing cycles for downward/upward
smoothing; \texttt{vlowstep} and \texttt{vhighstep} are the number
of smoothing passes at lowest/highest levels; and \texttt{max\_mg\_cycles}
is the maximum number of multigrid cycles per solve. All these parameters
should be small integers.
\begin{lyxcode}
Solver~=~multigrid~\\
mg\_cycle~=~1~\\
down\_heavy~=~2~\\
up\_heavy~=~2~\\
vlowstep~=~20~\\
vhighstep~=~2~\\
max\_mg\_cycles~=~50
\end{lyxcode}
The following parameter turn on the pre-conditioner for the matrix
equation solver (either multigrid or conjugate gradient).
\begin{lyxcode}
precond~=~on
\end{lyxcode}
The stiffness matrix uses augmented Lagrangian formulation to improve
the convergence for large viscosity variations \cite{Moresi/Zhong/Gurnis The accuracy}.
These parameters specify whether to enable the formulation and how
much weight to use for the formulation.
\begin{lyxcode}
aug\_lagr~=~on~\\
aug\_number~=~2.0e3
\end{lyxcode}
The discrete Stokes equations \ref{eq:discrete continuite eqn} and
\ref{eq:discrete momentum eqn} are solved using the Uzawa algorithm,
which iteratively updates the pressure and velocity solutions. Three
variations of the Uzawa algorithm are used in CitcomS, one for the
incompressible case, and the other two for the compressible case.
One parameter, \texttt{piterations,} common to the three variations,
specifies the maximum number of iterations and the desired residual
level for the continuity equation \ref{eq:discrete continuite eqn}.
Sometimes, larger value of \texttt{piterations} is required for convergence
if complicated velocity boundary conditions are used.
\begin{lyxcode}
piterations~=~375
\end{lyxcode}
For the compressible case, two choices of the Uzawa algorithm are
available. If \texttt{uzawa=cg}, the algorithm described in Equation
\ref{eq:iter-cg} is used. In this case, an additional parameter controls
the maximum number of outer iterations.
\begin{lyxcode}
uzawa~=~cg~\\
compress\_iter\_maxstep~=~100
\end{lyxcode}
If \texttt{uzawa=bicg}, the algorithm described in Equation \ref{eq:bicg}
is used, and no additional parameter is needed.

The overall accuracy of the velocity solver is controlled by a single
parameter. The solver is converged if the residuals of Equation \ref{eq:discrete continuite eqn}
and \ref{eq:discrete momentum eqn} both are smaller than \texttt{accuracy},
or if the changes in the velocity and pressure both are smaller than
\texttt{accuracy} for the last two iterations.
\begin{lyxcode}
accuracy~=~0.001
\end{lyxcode}
Finally, the net angular momentum of the velocity solution is removed.
The net angular momentum and rigid body rotation are unconstrained
by the Stokes equations, if free-slip boundary conditions are used
for the top and bottom boundaries in a full spherical model. That
is, you can add an arbitrary amount of rotation to the velocity solution,
and the resultant velocity is still a valid solution of the Stokes
equations. Since no external torque is applied to the mantle in free-slip
boundary conditions, the angular momentum must be constant in time.
Enforcing the angular momentum to be zero by removing it from the
velocity solution is often desirable.
\begin{lyxcode}
remove\_rigid\_rotation~=~off~\\
remove\_angular\_momentum~=~on
\end{lyxcode}
However, for models with imposed plate velocity, it is advisable to
turn off both \texttt{remove\_rigid\_rotation} and \texttt{remove\_angular\_momentum}.

\begin{lyxcode}
\#~Cookbook~8:~Compressible~Model~and~Checkpointing~\\
~\\
{[}CitcomS{]}~\\
solver~=~full~\\
steps~=~10000~\\
~\\
{[}CitcomS.controller{]}~\\
monitoringFrequency~=~1000~\\
checkpointFrequency~=~1000~\\
~\\
{[}CitcomS.solver{]}~\\
datafile~=~cookbook8~\\
datafile\_old~=~cookbook8~\\
~\\
rayleigh~=~7.68175583e4~\\
dissipation\_number~=~0.5~\\
gruneisen~=~1.25~\\
surfaceT~=~0.1~\\
~\\
{[}CitcomS.solver.mesher{]}~\\
nodex~=~33~\\
nodey~=~33~\\
nodez~=~33~\\
levels~=~5~\\
~\\
coor~=~1~\\
coor\_file~=~coord.dat~\\
~\\
{[}CitcomS.solver.param{]}~\\
reference\_state~=~1~\\
refstate\_file~=~ref.dat~\\
~\\
{[}CitcomS.solver.output{]}~\\
output\_optional~=~geoid,surf,botm~\\
use\_cbf\_topo~=~on~\\
self\_gravitation~=~on~\\
output\_ll\_max~=~20~\\
~\\
{[}CitcomS.solver.ic{]}~\\
tic\_method~=~3~\\
num\_perturbations~=~1~\\
perturbl~=~3~\\
perturbm~=~2~\\
perturblayer~=~17~\\
perturbmag~=~0.01~\\
~\\
restart~=~off~\\
solution\_cycles\_init~=~9000~\\
~\\
{[}CitcomS.solver.visc{]}~\\
VISC\_UPDATE~=~on~\\
visc\_smooth\_method~=~1~\\
~\\
TDEPV~=~on~\\
rheol~=~1~\\
num\_mat~=~4~\\
viscE~=~2.99573,2.99573,2.99573,2.99573~\\
viscT~=~0.5,0.5,0.5,0.5~\\
visc0~=~1,1,1,1~\\
~\\
VMIN~=~on~\\
visc\_min~=~0.001~\\
VMAX~=~on~\\
visc\_max~=~1e+06~\\
~\\
{[}CitcomS.solver.tsolver{]}~\\
finetunedt~=~0.75~\\
~\\
{[}CitcomS.solver.vsolver{]}~\\
uzawa~=~cg~\\
Solver~=~multigrid~\\
mg\_cycle~=~1~\\
down\_heavy~=~2~\\
up\_heavy~=~2~\\
vlowstep~=~20~\\
vhighstep~=~2~\\
max\_mg\_cycles~=~50~\\
~\\
aug\_lagr~=~on~\\
aug\_number~=~2.0e3~\\
precond~=~on~\\
~\\
piterations~=~375~\\
accuracy~=~0.001~\\
compress\_iter\_maxstep~=~100~\\
~\\
remove\_rigid\_rotation~=~off~\\
remove\_angular\_momentum~=~on
\end{lyxcode}

\subsection{Discussion}

The results for this problem are presented in Figure \ref{fig:Cookbook-8}.
A tetrahedral symmetric pattern is developed for the convection. The
surface heat flux $H_{surf}$ at the steady state is 3.892, and the
bottom heat flux $H_{botm}$ is 12.817. The heat flux imbalance ($H_{botm}r_{inner}^{2}/H_{surf}r_{outer}^{2}-1$)
is -0.38\%. The Nussult number, converted from the surface heat flux
by Equation \ref{eq:Nusselt number}, is 3.182. The total viscous
heating is 7.68991, and the total adiabatic cooling is 7.71719. Under
a steady state, these two terms should be exactly balanced.

The geoid is output as spherical harmonic coefficients. A post-processing
program is provided to project the geoid coefficients to a regular
(longitude, latitude) grid. To use a regular grid of $1^{\circ}$
spacing, run this command:
\begin{lyxcode}
\$~visual/project\_geoid~cookbook8.geoid.0.10000~geoid.xyz~361~181 \end{lyxcode} It will generate a file \texttt{geoid.xyz}. There are 3 columns in the file, which are (longitude, latitude, geoid). The geoid is in a unit of meters. Several dimensional constants are required for the geoid computation. These constants have sensible default values (in SI units) for the Earth. Note that the temperature drop from the core-mantle boundary to the surface ($\Delta T$in Equation \ref{eq:T dim}) is derived from these constants and the Rayleigh number. \begin{lyxcode} {[}CitcomS.solver.const{]}~\\ radius=6.371e+06~\\ density=3340~\\ thermdiff=1e-06~\\ gravacc=9.81~\\ thermexp=3e-05~\\ refvisc=1e+21~\\ cp=1200~\\ density\_above=1030~\\ density\_below=6600~ \end{lyxcode} \begin{figure}[H] \begin{centering} \includegraphics[scale=0.75]{graphics/cookbook8.png} \par\end{centering} \caption{\label{fig:Cookbook-8}Cookbook 8: The steady state temperature field at the 10,000th time step. A tetrahedral symmetric convection pattern is developed. Two temperature isosurfaces of 0.4 and 0.8 are shown. } \end{figure} \newpage{} \section{Cookbook 9: Nested Solver Coupling for Ridge-Plume Interaction} To run this cookbook, you need to install the Exchanger package. See Section \vref{sec:Exchanger-Configuration-(Optional)} for installation instructions. \subsection{Problem} In order to model the interaction between a segment of a mid-ocean ridge and a plume, a small grid spacing is required to resolve the thin plume conduit and the flow near the ridge. In addition, the plume must be placed sufficiently away from the sidewalls to avoid possible boundary effects. Satisfying both of these requirements in a model entails a long computation time. Using solver coupling, a high-resolution model with a smaller domain can be nested within a low-resolution model, and the computation time significantly reduced. \subsection{Solution} You will use two solvers in the model. A special command-line option is required for a coupled model. Type the following command to run this cookbook example: \begin{lyxcode} \$~citcoms~-{}-coupled~cookbook9.cfg
\end{lyxcode}
The embedded solver (esolver) is nested within the domain of the containing
solver (csolver). The velocity boundary conditions of the esolver
are fixed normal velocity and shear stress, whose values are received
from the esolver. The containing solver can be either a full CitcomS
solver or a regional one. The embedded solver can only be a regional
one. Both solvers are instances of a regional CitcomS solver in the
cookbook example.
\begin{lyxcode}
csolver~=~regional~~~~~~;~containing~solver~\\
esolver~=~regional~~~~~~;~embedded~solver
\end{lyxcode}
You will use four processors for this model. Each solver will use
two processors. The rank-0 and rank-1 processors will be assigned
to the csolver, and the rank-2 and rank-3 processors will be assigned
to the esolver.
\begin{lyxcode}
{[}CitcomS.layout{]}~\\
containing\_group~=~{[}0-1{]}~\\
embedded\_group~=~{[}2-3{]}
\end{lyxcode}
The procedure for the nested solver coupling is described elsewhere
\cite{Tan et al GeoFramework Part I} and only briefly described below:
\begin{enumerate}
\item The containing solver solves the velocity and computes its time-step
size $dt_{c}$.
\item The containing solver interpolates the velocity, stress and temperature
onto the boundary of the embedded model.
\item The embedded model solves the velocity and temperature, using the
boundary conditions from step 2, with a time-step size $dt_{e}$;
$dt_{e}$ is always smaller than $dt_{c}$.
\item Step 3 is repeated until the sum of $dt_{e}$ is equal to $dt_{c}$.
\item Optionally (if \texttt{two\_way\_communication} is \texttt{on}), the
containing solver updates its temperature according to the embedded
solver's temperature.
\item The containing solver solves the temperature with a time-step size
$dt_{c}$.
\item The entire process repeats beginning at step 1.
\end{enumerate}
The two solver processes are separated and only communicate through
the couplers, which in turn use the exchanger package to pass messages.
The containing coupler and controller (ccoupler and ccontroller) are
associated with the csolver, and the embedded coupler and controller
(ecoupler and econtroller) with the esolver. Each solver will track
its own number of time steps. The esolver has a smaller time-step
size and, hence, a larger number of time steps. The model will finish
when either of the solvers reaches the 200th time step. The csolver
will output for every 2 steps and the esolver for every 10 steps.
Additional output will occur when the embedded solver is at procedure
4 above.
\begin{lyxcode}
{[}CitcomS{]}~\\
steps~=~200~\\
~\\
{[}CitcomS.ccontroller{]}~\\
monitoringFrequency~=~2~\\
~\\
{[}CitcomS.econtroller{]}~\\
monitoringFrequency~=~10
\end{lyxcode}
A few parameters must be identical for the ccoupler and ecoupler.
You will use two-way communication, which enables the ecoupler to
send temperature information to the ccoupler. Otherwise the communication
is one way, ccoupler to ecoupler only, and the csolver is not affected
by the esolver.
\begin{lyxcode}
two\_way\_communication~=~on
\end{lyxcode}
There is an option to exchange initial temperature, which could ensure
that the initial temperature field of both solvers is consistent.
contain consistent temperature fields and don't need to exchange that
information again.
\begin{lyxcode}
exchange\_initial\_temperature~=~off
\end{lyxcode}
A good initial guess of pressure field can speed up the convergence
of the iterative Stokes solver. We will use the pressure solution
of the containing solver as the initial guess of pressure of the embedded
solver.
\begin{lyxcode}
exchange\_pressure~=~on
\end{lyxcode}
Three additional parameters are for the ecoupler only. The following
two parameters specify whether the top and bottom boundary of the
esolver will be imposed with the velocity and temperature from the
csolver. For models with prescribed top or bottom boundary conditions,
you can turn on these two parameters:
\begin{lyxcode}
exclude\_top~=~off~\\
exclude\_bottom~=~off
\end{lyxcode}
The next parameter deserves a little more explanation. The continuity
equation \ref{eq:conservation of mass} can be rearranged according
to Gauss's divergence theorem to become:
\begin{eqnarray}
\ointop_{surface\ area}\rho_{0}\overrightarrow{u}\cdot\overrightarrow{n}dA & = & 0\label{eq:gauss div}
\end{eqnarray}
This equation requires that the normal velocity $\overrightarrow{u}$
imposed on the embedded solver must have zero outflow. Due to finite
accuracy of the containing solver and the interpolation process, the
normal velocity received by the embedded solver might not have zero
outflow, which could cause a convergence problem in the embedded solver.
One can amend the received normal velocity slightly to satisfy Equation
\ref{eq:gauss div} by setting:
\begin{lyxcode}
amending\_outflow~=~on
\end{lyxcode}
Without turning on \texttt{side\_sbcs}, the model will fail to run.
\begin{lyxcode}
{[}CitcomS.esolver.bc{]}~\\
side\_sbcs~=~on
\end{lyxcode}
Plate motion is imposed on top of the esolver, which has a mid-ocean
ridge with a 5 cm/yr half-spreading rate. A transform fault cuts through
the ridge. Setting \texttt{start\_age} to 0, the plate motion is read
from the file \texttt{velocity/bvel.dat0} and will be constant in
time.
\begin{lyxcode}
{[}CitcomS.csolver.bc{]}~\\
topvbc~=~1~\\
~\\
{[}CitcomS.csolver.param{]}~\\
file\_vbcs~=~on~\\
vel\_bound\_file~=~velocity/bvel.dat~\\
start\_age~=~0
\end{lyxcode}
The domain of the csolver is bigger than that of the esolver. The
radial dimension of the esolver will shrink slightly with respect
to the csolver so that the domain of the esolver is completely inside
the domain of the csolver. The radial coordinates of the csolver is
refined near the top and bottom boundaries (\texttt{coor=2}), with
the lower 10\% of the mesh divided by 15\% of the nodes, the upper
10\% of the mesh divided by 20\% of the nodes, and the remaining 80\%
of the mesh divided by 75\% of the nodes.

Therefore, when setting the values of \texttt{coor\_refine} (see below),
be careful not to set the first or third value (in example below,
\texttt{0.1}) too low, or the second or fourth value (in example below,
\texttt{0.15} and \texttt{0.2}) too high. This can cause the mesh
size to change too rapidly across an element, which will generate
pronounced numerical artifacts in the temperature solver.
\begin{lyxcode}
{[}CitcomS.csolver.mesher{]}~\\
nprocz~=~2~\\
nodex~=~17~\\
nodey~=~17~\\
nodez~=~25~\\
theta\_min~=~1.171~\\
theta\_max~=~1.971~\\
fi\_min~=~0~\\
fi\_max~=~0.8~\\
coor~=~2~\\
coor\_refine~=~0.1,0.15,0.1,0.2~\\
~\\
{[}CitcomS.esolver.mesher{]}~\\
nprocz~=~2~\\
nodex~=~21~\\
nodey~=~21~\\
nodez~=~49~\\
theta\_min~=~1.421~\\
theta\_max~=~1.721~\\
fi\_min~=~0.25~\\
fi\_max~=~0.55
\end{lyxcode}
The temperature solver has several parameters to tune its behavior.
The temperature solver is turned on (\texttt{ADV}). A Lenardic-type
filter, which removes numerical artifacts while keeping total energy
conserved, is disabled (\texttt{filter\_temp}). The maximum temperature
can be monitored between time steps (\texttt{monitor\_max\_T}). If
the maximum temperature increases too much (> 5\%) between time steps,
the temperature solver will rerun with half time-step size. The time-step
size is usually determined dynamically according to Courant criterion
and is reduced by a fraction (\texttt{finetunedt}) to improve accuracy.
The time-step size can also be specified statically in \texttt{fixed\_timestep},
if \texttt{fixed\_timestep} is non-zero. The temperature solver uses
an explicit predictor-corrector algorithm. Using 0.5 for the predictor
this algorithm is second-order accurate.
\begin{lyxcode}
filter\_temp~=~off~\\
monitor\_max\_T~=~off~\\
finetunedt~=~0.9~\\
fixed\_timestep~=~0~\\
\end{lyxcode}
The purpose of the csolver is to provide the proper boundary condition
for the esolver, and we don't need high accuracy for that. So we will
use low accuracy for the csolver and moderate accuracy for the esolver.
\begin{lyxcode}
{[}CitcomS.csolver.vsolver{]}~\\
accuracy~=~2e-2~\\
~\\
{[}CitcomS.esolver.vsolver{]}~\\
accuracy~=~1e-3~\\

\end{lyxcode}

\subsubsection{Example: Nested Solver Coupling for Ridge-Plume Interaction, cookbook9.cfg}
\begin{lyxcode}
\#~Cookbook~9:~Nested~Solver~Coupling~for~Ridge-Plume~Interaction~\\
~\\
{[}CitcomS{]}~\\
steps~=~200~\\
csolver~=~regional~~~~~~;~containing~solver~\\
esolver~=~regional~~~~~~;~embedded~solver~\\
~\\
{[}CitcomS.ccoupler{]}~\\
two\_way\_communication~=~on~\\
exchange\_initial\_temperature~=~off~\\
exchange\_pressure~=~on~\\
~\\
{[}CitcomS.ecoupler{]}~\\
two\_way\_communication~=~on~\\
exchange\_initial\_temperature~=~off~\\
exchange\_pressure~=~on~\\
amending\_outflow~=~on~\\
exclude\_top~=~off~\\
exclude\_bottom~=~off~\\
~\\
\#~assigning~processors~to~csolver~and~esolver~\\
{[}CitcomS.layout{]}~\\
containing\_group~=~{[}0-1{]}~\\
embedded\_group~=~{[}2-3{]}~\\
~\\
{[}CitcomS.ccontroller{]}~\\
monitoringFrequency~=~2~\\
~\\
{[}CitcomS.econtroller{]}~\\
monitoringFrequency~=~10~\\
~\\
{[}CitcomS.csolver{]}~\\
datafile\_old~=~cntn~\\
datafile~=~cntn~\\
rayleigh~=~1e8~\\
~\\
{[}CitcomS.csolver.mesher{]}~\\
nprocz~=~2~\\
nodex~=~17~\\
nodey~=~17~\\
nodez~=~25~\\
theta\_min~=~1.171~\\
theta\_max~=~1.971~\\
fi\_min~=~0~\\
fi\_max~=~0.8~\\
coor~=~2~\\
coor\_refine~=~0.1,0.15,0.1,0.2~\\
~\\
{[}CitcomS.csolver.bc{]}~\\
topvbc~=~1~\\
~\\
{[}CitcomS.csolver.param{]}~\\
file\_vbcs~=~on~\\
vel\_bound\_file~=~velocity/bvel.dat~\\
start\_age~=~0~\\
~\\
{[}CitcomS.csolver.ic{]}~\\
tic\_method~=~-1~\\
solution\_cycles\_init~=~0~\\
~\\
{[}CitcomS.csolver.tsolver{]}~\\
filter\_temp~=~off~\\
monitor\_max\_T~=~off~\\
finetunedt~=~0.9~\\
fixed\_timestep~=~0~\\
~\\
{[}CitcomS.csolver.vsolver{]}~\\
accuracy~=~2e-2~\\
~\\
{[}CitcomS.csolver.visc{]}~\\
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
viscE~=~1.6,1.6,1.6,1.6~\\
viscT~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0~\\
~\\
\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#\#~\\
~\\
{[}CitcomS.esolver{]}~\\
datafile\_old~=~embd~\\
datafile~=~embd~\\
rayleigh~=~1e8~\\
~\\
{[}CitcomS.esolver.mesher{]}~\\
nprocz~=~2~\\
nodex~=~21~\\
nodey~=~21~\\
nodez~=~49~\\
theta\_min~=~1.421~\\
theta\_max~=~1.721~\\
fi\_min~=~0.25~\\
fi\_max~=~0.55~\\
~\\
{[}CitcomS.esolver.ic{]}~\\
tic\_method~=~-1~\\
solution\_cycles\_init~=~0~\\
~\\
{[}CitcomS.esolver.vsolver{]}~\\
accuracy~=~1e-3~\\
~\\
{[}CitcomS.esolver.tsolver{]}~\\
filter\_temp~=~off~\\
monitor\_max\_T~=~off~\\
finetunedt~=~0.9~\\
fixed\_timestep~=~0~\\
~\\
{[}CitcomS.esolver.visc{]}~\\
VISC\_UPDATE~=~on~\\
num\_mat~=~4~\\
visc0~=~1,1,1,1~\\
TDEPV~=~on~\\
viscE~=~1.6,1.6,1.6,1.6~\\
viscT~=~0,0,0,0~\\
VMIN~=~on~\\
visc\_min~=~1.0~\\
VMAX~=~on~\\
visc\_max~=~100.0~\\
~\\
{[}CitcomS.esolver.bc{]}~\\
side\_sbcs~=~on
\end{lyxcode}

\subsection{Discussion}

The solution of this problem is presented in Figure \ref{fig:Cookbook-9}.
is elongated in the ridge-parallel direction. The lithosphere subducts
at the left and right edges of the csolver. If the domain of the csolver
were bigger, the plume would be further away from the subducted slabs.

\begin{figure}[H]
\begin{centering}
\includegraphics[width=0.6\paperwidth]{graphics/cookbook9.png}
\par\end{centering}

lithosphere, and the plume conduit is elongated in the ridge-parallel
direction. The temperature isosurface is at 0.8. The grid spacings
of both meshes are reduced approximately threefold for better visualization.}
\end{figure}

\newpage{}

\section{Cookbook 10: Synthetic Seismograms from Mantle Convection Models}

\subsection{Problem}

This example is a step-up of Cookbook \ref{fig:Cookbook-7:-The},
thermo-chemical convection within a full spherical shell domain. Composition
heterogeneity exists in the Earth mantle. The density anomalies due
to the composition heterogeneity, as well as due to the thermal heterogeneity,
drive the convection flow. The same composition heterogeneity causes
seismic velocity anomalies. The goal is to convert the mantle convection
models, with the aid of a mineral physics model, to seismic anomalies,
and then generate synthetic seismograms using SPECFEM3D\_GLOBE.

\subsection{Solution}

This cookbook uses CitcomS and SPECFEM3D\_GLOBE. SPECFEM3D\_GLOBE
is a spectral element code to simulate seismic wave propogation in
the code from this link \url{geodynamics.org/cig/software/packages/seismo/specfem3d-globe/}
and install it on your machine, which requires a Fortran 90 compiler
and a cluster, or you can go to CIG Seismology Web Portal \url{https://crust.geodynamics.org/portals/seismo/}
to launch a job, which only requires a browser. We will use the later
approach in the cookbook.

\subsubsection{Running CitcomS}

Most of the parameters are copied from Cookbook \ref{fig:Cookbook-7:-The}.
We will highlight a few differences. We use a higher Rayleigh number.
\begin{lyxcode}
rayleigh~=~1e9
\end{lyxcode}
The radius of the bottom surface is set to 0.546 in the cookbook,
which must be smaller than the CMB radius in SPECFEM3D\_GLOBE. (If
PREM is used as the 1D reference model in SPECFEM3D\_GLOBE, the CMB
\begin{lyxcode}
\end{lyxcode}
We convert the temperature and composition fields to seismic velocities
using a mineral phyiscs model. Currently, there is only one model
implemented, based on the model of \textit{Trampert, Vacher, and Vlaar}\cite{Trampert mineral physics model 2001}.
In their paper, the (depth dependent) temperature and composition
derivatives of seismic velocities are given as polynomial coefficients.
Note that the paper has its own reference profile, and is only valid
between 1000 km < depth < 2600 km. We will extend the model to whole
mantle and use PREM as the reference profile. This clearly is over-stretching
the model, but we are using it for demonstration purposes.
\begin{lyxcode}
mineral\_physics\_model~=~3
\end{lyxcode}
We will output the converted density and seismic velocities for SPECFEM3D\_GLOBE
consumption. This option enables the binary seismic output (\texttt{{*}.seismic.{*}}
files), as well as other binary output (\texttt{{*}.domain} and \texttt{{*}.coord\_bin.{*}}
files) that SPECFEM3D\_GLOBE will read. The formats of these binary
files are provided in Appendix \ref{sec:Misc.-Binary-Output}.
\begin{lyxcode}
output\_optional~=~tracer,~comp\_nd,~seismic
\end{lyxcode}

\subsubsection{Example: Synthetic Seismograms from Mantle Convection Model, cookbook10.cfg}
\begin{lyxcode}
\#~Cookbook~10:~Synthetic~Seismograms~from~Mantle~Convection~Models

{[}CitcomS{]}

solver~=~full

steps~=~15~\\
~\\
{[}CitcomS.controller{]}

monitoringFrequency~=~5~\\
~\\
{[}CitcomS.solver{]}

datafile~=~cookbook10

rayleigh~=~1e9~\\
~\\
{[}CitcomS.solver.mesher{]}

~\\
{[}CitcomS.solver.ic{]}

num\_perturbations~=~1

perturbl~=~3

perturbm~=~2

perturblayer~=~5

perturbmag~=~0.05~\\
~\\
{[}CitcomS.solver.output{]}

output\_optional~=~tracer,comp\_nd,seismic~\\
~\\
{[}CitcomS.solver.param{]}

mineral\_physics\_model~=~3~\\
~\\
{[}CitcomS.solver.tracer{]}

tracer~=~on

tracer\_ic\_method~=~0

tracers\_per\_element~=~20

tracer\_file~=~tracer.dat~\\
~\\
tracer\_flavors~=~2

ic\_method\_for\_flavors~=~0

z\_interface~=~0.7~\\
~\\
chemical\_buoyancy~=~1

buoy\_type~=~1

buoyancy\_ratio~=~0.5~\\
~\\
regular\_grid\_deltheta~=~1.0

regular\_grid\_delphi~=~1.0~\\
~\\
{[}CitcomS.solver.vsolver{]}

accuracy~=~1e-04

vlowstep~=~1000

piterations~=~1000~\\
~\\
\#~Assign~the~viscosities.

{[}CitcomS.solver.visc{]}

VISC\_UPDATE~=~on

num\_mat~=~4

visc0~=~1,1,1,1

TDEPV~=~on

rheol~=~4

viscE~=~0.2,0.2,0.2,0.2

viscT~=~0,0,0,0

viscZ~=~0,0,0,0

VMIN~=~on

visc\_min~=~1.0

VMAX~=~on

visc\_max~=~100.0
\end{lyxcode}

\subsubsection{Running SPECFEM3D\_GLOBE from CIG Seismology Web Portal}

Using your favorite browser, go to the CIG Seismology Web Portal \url{https://crust.geodynamics.org/portals/seismo/}.
Free registration is required to use the portal. Your initial account
has 10,000 SUs (CPU-hours) allocated. A short (20 minutes of seismograms),
low-resolution SPECFEM3D\_GLOBE run consumes about 450 SUs. A short,
high-resolution SPECFEM3D\_GLOBE run consumes about 1,000 SUs. Further
SUs can be allocated upon request.

After you login to the portal, take a few minutes to read the instructions
and familiarize yourself with the interface. You will see a toolbar
with big icons on the top and several grayed-out tabs below the toolbar.
The toolbar allows you to specify the earthquake source (\texttt{Events}),
the receivers (\texttt{Stations}), and simulation parameters (\texttt{3D
Parameters}), including the resolution (\texttt{Mesh}) and the seismic
velocity model for the crust and mantle (\texttt{3D Models}). In the
cookbook, we will upload the result of CitcomS as a 3D Model and upload
our own fictitious stations and event.

\begin{center}
\par\end{center}

Click on the \texttt{3D Models} icon (1), and the \texttt{Upload}
tab (2), then click on the link \texttt{citcoms\_isotropic\_no\_crust.tgz}
seismic velocities to SPECFEM3D\_GLOBE. However, the tarball does
not contains the actual CitcomS output, which will be provided by
you later. As you might have guessed from the filename, the code will
generate an isotropic mantle without an overlay of a 3D crust, such
as the CRUST2.0 model. The 3D crust model obeys the continent-ocean
distribution of current Earth. Since our CitcomS model has no notation
on where is the Pacific and where is Asia, we do not want to impose
a 3D crust model. You can untar the file by running:
\begin{lyxcode}
\$~gunzip~-c~citcoms\_isotropic\_no\_crust.tgz~|~tar~xf~- \end{lyxcode} You will supply the CitcomS output to the code. Copy these CitcomS output files \texttt{cookbook.domain}, \texttt{cookbook.coord\_bin.{*}}, and \texttt{cookbook.seismic.{*}.15} to the directory \texttt{citcoms\_isotropic\_no\_crust/shared/}, so that the content of the directory looks like: \begin{lyxcode} \$~ls~citcoms\_isotropic\_no\_crust/shared/

cookbook10.coord\_bin.0~~~cookbook10.seismic.0.15~

cookbook10.coord\_bin.1~~~cookbook10.seismic.10.15~

cookbook10.coord\_bin.10~~cookbook10.seismic.11.15~

cookbook10.coord\_bin.11~~cookbook10.seismic.1.15~

cookbook10.coord\_bin.2~~~cookbook10.seismic.2.15~

cookbook10.coord\_bin.3~~~cookbook10.seismic.3.15~

cookbook10.coord\_bin.4~~~cookbook10.seismic.4.15~

cookbook10.coord\_bin.5~~~cookbook10.seismic.5.15~

cookbook10.coord\_bin.6~~~cookbook10.seismic.6.15~

cookbook10.coord\_bin.7~~~cookbook10.seismic.7.15~

cookbook10.coord\_bin.8~~~cookbook10.seismic.8.15~

cookbook10.coord\_bin.9~~~cookbook10.seismic.9.15~

cookbook10.domain
\end{lyxcode}
These files contain the coordinate, density, and seismic velocities
data. We will need to tell SPECFEM3D\_GLOBE where to find the CitcomS
output. Edit the header file \texttt{citcoms\_isotropic\_no\_crust/citcoms\_parm.h}
to become:
\begin{lyxcode}
/{*}~filename~prefix~of~citcoms~output~{*}/

const~char~citcoms\_model\_filename\_base{[}{]}~=~\textquotedbl{}@THIS\_DIR@/shared/cookbook10\textquotedbl{};

/{*}~time~step~of~citcoms~output~{*}/

const~int~citcoms\_step~=~15;
\end{lyxcode}
The portal has a few conventions on where to put the data. First,
the magic string \texttt{@THIS\_DIR@} will be replaced by the absolute
path when the tar file is expanded by the portal. Second, the tar
file will be distributed to all computer nodes when portal launches
the job. Third, the subdirectory \texttt{shared/} is special such
that files in this directory will not be distributed to all computer
nodes. Data files of huge size should be put in this directory. You
might want to edit \texttt{citcoms\_isotropic\_no\_crust/description.txt}
to provide a description of your model. The text in this file will
be displayed on the portal. The final step before uploading is to
create a new tar file containing the code and data.
\begin{lyxcode}
\$~tar~cf~cookbook10.tar~citcoms\_isotropic\_no\_crust/ \$~gzip~cookbook10.tar
\end{lyxcode}
Now, you can upload the data by clicking on the \texttt{Browse} button
(4), choosing \texttt{cookbook10.tar.gz} from the resulting dialog
box, then clicking on the \texttt{Upload} button. After the portal
confirms that the upload is successful, click on the \texttt{3D Models}
icon again; you will see \texttt{cookbook10} is among the available
3D models.

\begin{center}
\includegraphics[scale=0.6]{graphics/portal-param-a.png}
\par\end{center}

You have to select a few parameters for a SPECFEM3D\_GLOBE run. Click
on the \texttt{3D Parameters} icon (5), the \texttt{New} tab (6),
select the parameters as shown in the figure above, and save the parameters
for later use. We will use a global mesh that is capable of resolving
seismic waves of 27 second period. Since CitcomS mesh is a perfect
sphere, you need to disable \texttt{topography} and \texttt{ellipticity}
to have a matching mesh in SPECFEM3D\_GLOBE. The \texttt{oceans} is
disabled because we don't have a 3D crust.

We will provide our own fictitious earthquake source and seismic stations.
The earthquake source is defined in \texttt{cmt-solution.cookbook10}
in the cookbook directory. The source is located at 15 km depth of
$0^{\circ}$N, $0^{\circ}$E. The half duration of the source is 0
second, which means a step function is used. Click on the \texttt{Events}
The stations are defined in \texttt{stations.cookbook10} in the cookbook
directory. The stations are on the Equator and are 70 to 100 degrees
away from the source. Click on the \texttt{Stations} icon and the

\begin{center}
\includegraphics[scale=0.75]{graphics/portal-run-a.png}
\par\end{center}

Finally, we are ready to launch a SPECFEM3D\_GLOBE simulation. Click
on the \texttt{Runs} icon (7), then the \texttt{New} tab (8), and
select our recently uploaded event, stations and parameters. Set the
record length to 30 minutes, which will give us synthetic seismograms
of 30 minutes. Click the \texttt{Save} button to save the settings.

\begin{center}
\includegraphics[scale=0.6]{graphics/portal-submit2-a.png}
\par\end{center}

You will be presented with a summary of the simulation settings, including
a map showing the location of the events and stations. After you have
reviewed all settings, scroll to the bottom of the page. Click on
the \texttt{Start} button. The request for the simulation is sent
to the TeraGrid supercomputing center. Your simulation will be waiting
in the queue and will be launched within several hours. You can log
off the portal for now.

\begin{center}
\includegraphics[scale=0.6]{graphics/portal-result.png}
\par\end{center}

When your simulation is finished, the portal will send you an e-mail
where there is a list of finished and unfinished runs. Click on the
finished run, and you will see a list of output. Click on the \texttt{specfem3dglobe.tar.gz}
\begin{lyxcode}
\$~gunzip~-c~specfem3dglobe.tar.gz~|~tar~xf~- \end{lyxcode} You will find the synthetic seismograms \texttt{{*}.sem.sac} in the directory \texttt{OUTPUT\_FILES/}. \subsection{Discussion} The dimensional value of total temperature contrast across the mantle is derived from the Rayleigh number in Equation \ref{eq:Ra, Rayleigh number} and the dimensional parameters in the \texttt{{[}CitcomS.solver.const{]}} section. In the very beginning of the simulation, the total temperature contrast across the mantle is displayed on the screen. This value will be used to scale the non-dimensional temperature when converting the temperature anomaly to seismic velocities. \begin{lyxcode} Total~temperature~contrast~=~3934.055176~K \end{lyxcode} The event and stations are all on the equatorial plane. Figure \ref{fig:cookbook10} shows the temperature field of the equatorial slice. A large warm upwelling is beneath$60^{\circ}$E. We should be able to observe its effect on delaying teleseismic waves in the$70^{\circ}$to$90^{\circ}$range. \noindent \begin{center} \begin{figure}[H] \begin{centering} \includegraphics[scale=0.6]{graphics/cookbook10.png} \par\end{centering} \caption{\label{fig:cookbook10}The equatorial slice of temperature from the CitcomS result. Three S wave paths probing the mantle are shown. } \end{figure} \par\end{center} The synthetic seismograms need to be post-processed. The post-processing programs we will use can be downloaded from the portal. Follow the instructions on the portal to build the programs. We will band-pass filter within the period of 30 and 500 seconds, and convolve the synthetics with a gaussian source time function of 50 seconds half duration, and append \texttt{conv} to the output filenames. \begin{lyxcode} \$~process\_syn.pl~-S~-t~30/500~-h~50~-x~conv~{*}.sac
\end{lyxcode}
Optionally, you can rotate the synthetics from East-North components
\begin{lyxcode}
\$~rotate.pl~{*}LHE{*}.conv \end{lyxcode} \part{Appendices} \appendix \chapter{\label{cha:Appendix-A:-Input}Input Parameters for CitcomS} \section{Input Parameters Grouped by Functionality} This section explains the meaning of the input parameters for CitcomS. These parameters are grouped by their functionality. Parameters are given with their default values. \subsection{Parameters that Control Input Files} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{reference\_state=1}}~\\ \texttt{\small{refstate\_file=\textquotedbl{}refstate.dat\textquotedbl{}}} & If \texttt{\small{reference\_state}} is set to \texttt{\small{1}}, a simple reference state of$\rho_{r}=\exp\left((1-r)D_{i}/\gamma\right)$is used, with constant gravity, thermal expansivity, and heat capacity. If \texttt{\small{reference\_state}} is set to \texttt{\small{0}}, the reference state is read from a file \texttt{\small{refstate\_file}}\texttt{.} \tabularnewline \hline \texttt{\small{mineral\_physics\_model=3}} & Using which mineral physics model to convert the seismic velocities.\tabularnewline \hline \texttt{\small{file\_vbcs=off}}~\\ \texttt{\small{vel\_bound\_file=\textquotedbl{}bvel.dat\textquotedbl{}}} & If \texttt{\small{file\_vbcs}} is set to \texttt{\small{on}}, the top surface velocity boundary conditions are read in from files which have location and name specified by \texttt{\small{vel\_bound\_file}}\texttt{.} Requires setting \texttt{\small{topvbc=1}} to take effect. If you wish to have a uniform top surface velocity boundary condition or some simple geometric pattern, then \texttt{\small{file\_vbcs}} should be set to zero. \tabularnewline \hline \texttt{\small{mat\_control=off}}\\ \texttt{\small{mat\_file=\textquotedbl{}mat.dat\textquotedbl{}}} & If \texttt{\small{mat\_control}} is set to \texttt{\small{on}}, then the time- and positional-dependent viscosity factor is defined from the files specified by \texttt{\small{mat\_file}}. These parameters allow you to define the material group of each element (such as a moving weak zone). Not working in this version.\tabularnewline \hline \texttt{\small{lith\_age=off}}\\ \texttt{\small{lith\_age\_file=\textquotedbl{}age.dat\textquotedbl{}}}\\ \texttt{\small{lith\_age\_time=off}} & If \texttt{\small{lith\_age}} is set to \texttt{\small{on}}, then the age of each surface nodes is read from the files specified by \texttt{\small{lith\_age\_file}}. These parameters control the thermal age of the top thermal boundary condition. If \texttt{\small{lith\_age\_time}} is \texttt{\small{on}}, the files are time-dependent.\tabularnewline \hline \end{tabular} \subsection{Parameters that Control Output Files} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{output\_format=ascii}} & Choose the format and layout of the output files. Can be either \texttt{ascii}, \texttt{ascii-gz} or \texttt{hdf5}. If \texttt{ascii-gz} is chosen, the code places gzipped files into \texttt{data\_dir}, and will put all time-step output into subdirectories of \texttt{data\_dir}. The same naming logic holds for reading old velo files.\tabularnewline \hline \texttt{\small{output\_optional=\textquotedbl{}surf,}}~\\ \texttt{\small{botm,tracer\textquotedbl{}}} & Choose additional output, including \texttt{surf}, \texttt{botm}, \texttt{geoid}, \texttt{seismic}, \texttt{stress}, \texttt{pressure}, \texttt{connectivity}, \texttt{horiz\_avg, tracer, heating, comp\_el} and\texttt{ comp\_nd}.\tabularnewline \hline \texttt{\small{datadir=\textquotedbl{}.\textquotedbl{}}} & Controls the location of output files. \tabularnewline \hline \texttt{\small{datafile=\textquotedbl{}regtest\textquotedbl{}}} & Controls the prefix of output file names such as \texttt{regtest.xxx}. Cannot contain the \texttt{/}'' character if \texttt{output\_format=ascii}.\tabularnewline \hline \texttt{\small{storage\_spacing=10}}\\ (in non-Pyre version)\\ or\\ \texttt{\small{monitoringFrequency=100}}\\ (in Pyre version) & Controls the interval between output files. CitcomS dynamically determines the size of the time step; this means that you might not get an output at the exact time required, but you can always get close depending on how small this number is. Do not make this number too small since outputs slow the code down and you may end up with an unmanageable number of output files. \tabularnewline \hline \texttt{\small{checkpointFrequency=100}} & The time-step interval between checkpoint output, which can be used later to resume the computation.\tabularnewline \hline \texttt{\small{output\_ll\_max=20}} & This parameter controls the maximum degree of spherical harmonics coefficients for geoid output.\tabularnewline \hline \texttt{self\_gravitation=off} & Considering the effect the self gravitation on the geoid or not.\tabularnewline \hline \texttt{use\_cbf\_topo=off} & Using the Consistent Boundary Flux (CBF) method to compute the dynamic topography or not.\tabularnewline \hline \end{tabular} \subsection{Mesh and Processors Setup } \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{nproc\_surf=1}} & This specifies the number of spherical caps of the mesh; must be \texttt{\small{1}} for regional spherical model and \texttt{\small{12}} for full spherical model.\tabularnewline \hline \texttt{\small{nprocx=1}}~\\ \texttt{\small{nprocy=1}}~\\ \texttt{\small{nprocz=1}} & These specify the number of processors in each spherical cap. \medskip{} For a full spherical model, \texttt{nprocx} must be equal to \texttt{nprocy}\tabularnewline \hline \texttt{\small{nodex=9}}~\\ \texttt{\small{nodey=9}}~\\ \texttt{\small{nodez=9}} & These specify the number of FEM nodes in each spherical cap. These parameters are not used in the C version if multigrid solver is used. \medskip{} For a full spherical model, \texttt{nodex} must be equal to \texttt{nodey}\tabularnewline \hline \texttt{\small{mgunitx=8}}~\\ \texttt{\small{mgunity=8}}~\\ \texttt{\small{mgunitz=8}}~\\ \texttt{\small{levels=1}} & These specify the nested level of multigrid units. Used by multigrid solver only. These parameters are not completely independent to each other. The constraint of Equation \ref{eq:mgunit} must be satisfied. These parameters are not used in the C version if conjgrad solver is used. \medskip{} For a full spherical model, \texttt{mgunitx} must be equal to \texttt{mgunity}\tabularnewline \hline \end{tabular} \subsection{\noindent Domain Size} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{theta\_min=1.0708}}~\\ \texttt{\small{theta\_max=2.0708}}~\\ \texttt{\small{fi\_min=0}}~\\ \texttt{\small{fi\_max=1}} & These parameters specify the horizontal extent of the computational domain. \texttt{theta\_min} and \texttt{theta\_max} are the colatitude measured in radians from the north pole. \texttt{fi\_min} and \texttt{fi\_max} are the longitudes measured from the prime meridian eastward in radians. Only in regional CitcomS.\tabularnewline \hline \texttt{\small{radius\_inner=0.55}}~\\ \texttt{\small{radius\_outer=1.0 }} & These parameters specify the radial extent of the computational domain. \texttt{radius\_inner} and \texttt{radius\_outer} are the inner and outer radii in non-dimensional units. It is probably more convenient to normalize lengths by the radius of the Earth. If you do this, then the Rayleigh number must be calculated with the radius of the Earth, not the thickness of the mantle. The core mantle boundary is close to having a non-dimensional radius of 0.55.\tabularnewline \hline \texttt{coor}\texttt{\small{=0}}{\small \par} \medskip{} \texttt{coor\_file=\textquotedbl{}coor.dat\textquotedbl{}}~\\ \texttt{coor\_refine=0.1,0.15,0.1,0.2} & If \texttt{coor=0}, there will be uniform mesh in the latitudinal, longitudinal, and radial directions. If \texttt{coor=1}, the coordinate is reading from the file specified by \texttt{coor\_file}. This is used to have a regular but uneven spacing between elements. If \texttt{coor=}2, there will be uniform mesh in the latitudinal and longitudinal directions. The mesh in the radial direction is generated according to \texttt{coor\_refine}, which is a vector of 4 numbers. The 1st value of \texttt{coor\_refine} specifies the radius fraction of the bottom layer, the 2nd value specifies the fraction of the nodes in the bottom layer, the 3rd value specifies the top layer fraction, and the last value specifies the top layer node fraction.\tabularnewline \hline \end{tabular} \subsection{Restarting the Code} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{restart=off}} & If \texttt{restart} is \texttt{\small{on}}, each processor will resume the computation from the checkpoint files.\tabularnewline \hline \texttt{\small{post\_p=off}} & Similar to \texttt{restart}, except that the model will then only run for 1 time step, which can be useful to regenerate the flow field and calculate the associated observables.\tabularnewline \hline \texttt{\small{datadir\_old=\textquotedbl{}.\textquotedbl{}}}~\\ \texttt{\small{datafile\_old=\textquotedbl{}regtest\textquotedbl{}}}~\\ \texttt{\small{solution\_cycles\_init=0}} & If \texttt{restart} is \texttt{\small{on}}, for example, processor \#5 will read its initial conditions from checkpoint file \texttt{regtest.chkpt.5.0} in this case. \tabularnewline \hline \end{tabular} \subsection{Run Length} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{minstep=1}}~\\ \texttt{\small{maxtotstep=1000000}}~\\ (only in pure C version)\texttt{\small{}}~\\ or\texttt{\small{}}~\\ \texttt{\small{steps=1}}~\\ (only in Pyre version) & The maximum and minimum number of time steps for the model, including the 0th time step.\tabularnewline \hline \texttt{\small{cpu\_limits\_in\_seconds=}}~\\ \texttt{\small{~~~360000000}} & Controls the termination of the code based on total wall clock time used. Available only in pure C version.\tabularnewline \hline \end{tabular} \subsection{Initial Conditions} \begin{flushleft} \begin{longtable}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{tic\_method=0}} & Which method to use to generate the initial temperature field.\tabularnewline \hline \texttt{\small{datadir\_old=\textquotedbl{}.\textquotedbl{}}}~\\ \texttt{\small{datafile\_old=\textquotedbl{}regtest\textquotedbl{}}}~\\ \texttt{\small{solution\_cycles\_init=0}}~\\ \texttt{\small{zero\_elapsed\_time=on}} & If \texttt{\small{tic\_method=-1}}, the initial temperature is read from files. For example, processor \#5 will read its initial temperature from old velo file \texttt{\small{regtest.velo.5.0}} in this case. If \texttt{\small{zero\_elapsed\_time}} is \texttt{\small{on}}, the initial time is set to zero. If it is \texttt{\small{off}} and \texttt{\small{tic\_method=-1}}, the initial time is read in from the old velo files. Note that this option has no effect when \texttt{\small{restart=on}}.\tabularnewline \hline \texttt{\small{num\_perturbations=1}}~\\ \texttt{\small{perturbmag=0.05}}~\\ \texttt{\small{perturbl=1}}~\\ \texttt{\small{perturbm=1}}~\\ \texttt{\small{perturblayer=5}}~\\ \texttt{\small{}}~\\ \texttt{\small{half\_space\_age=40}}~\\ \texttt{\small{mantle\_temp=1.0}}~\\ \texttt{\small{}}~\\ \texttt{\small{blob\_center={[}-999,-999,-999{]}}}~\\ \texttt{\small{blob\_radius=0.063}}~\\ \texttt{\small{blob\_dT=0.18}} & If \texttt{tic\_method}\texttt{\small{=0}}. The initial temperature is a linear temperature gradient with perturbations at specific layers, where \texttt{num\_perturbations} specifies the number of perturbations, and \texttt{perturblayer} specifies the layers to be perturbed, representing the number of the mesh node in radial direction. There must be as many entries as \texttt{num\_perturbations} in a comma-separated list. The perturbation added to each layer is given by: $mag\times\cos(m\phi)\times P_{lm}(\cos\theta)$ for the full sphere, and by: $mag\times\cos(\frac{\left(\theta-\theta_{min}\right)l\pi}{\theta_{max}-\theta_{min}})\times\cos(\frac{\left(\phi-\phi_{min}\right)m\pi}{\phi_{max}-\phi_{min}})$ for the regional sphere.\\ If \texttt{tic\_method}\texttt{\small{=}}1, T is \texttt{1} everywhere, except a cold thermal boundary layer at the top, whose temperature is determined by the half-space cooling model and \texttt{\small{half\_space\_age}} (in million of years, Myrs). \\ If \texttt{tic\_method}\texttt{\small{=2,}} T is \texttt{\small{mantle\_temp}} everywhere, except for a warm spherical blob and a cold thermal boundary layer at the top, whose temperature is determined by the half-space cooling model and \texttt{\small{half\_space\_age}} (in Myrs). The location of the blob is default to the center of the computational domain.\\ If \texttt{tic\_method}\texttt{\small{=}}3, the initial temperature is a conductive profile with perturbations to all layers. The perturbation is given by: $mag\times\sin\left(\frac{(r-r_{in})\pi}{r_{out}-r_{in}}\right)\left(\sin(m\phi)+\cos(m\phi)\right)P_{lm}(\cos\theta)$ for the full sphere, and by: $mag\times\sin\left(\frac{(r-r_{in})\pi}{r_{out}-r_{in}}\right)\times\cos(\frac{\left(\theta-\theta_{min}\right)l\pi}{\theta_{max}-\theta_{min}})\times\cos(\frac{\left(\phi-\phi_{min}\right)m\pi}{\phi_{max}-\phi_{min}})$ for the regional sphere.\\ If \texttt{tic\_method}\texttt{\small{=}}4, the initial temperature is read from grd files.\\ If \texttt{tic\_method}\texttt{\small{=10,}} T is \texttt{\small{mantle\_temp}} everywhere, except for a cold thermal boundary layer at the top and perturbations at all layers, similar to \texttt{tic\_method=3}.\\ If \texttt{tic\_method}\texttt{\small{=11,}} T is \texttt{\small{mantle\_temp}} everywhere, except for a hot thermal boundary layer at the bottom and perturbations at all layers, similar to \texttt{tic\_method=3}.\\ If \texttt{tic\_method}\texttt{\small{=12,}} T is \texttt{\small{mantle\_temp}} everywhere, except for a cold thermal boundary layer, a hot thermal boundary layer at the bottom and perturbations at all layers, similar to \texttt{tic\_method=3}. If \texttt{tic\_method}\texttt{\small{=90,}} T is \texttt{\small{0}} everywhere, except for a single perturbation at the middle layer. This initial temperature is good for comparison with analytical solutions.\tabularnewline \hline \end{longtable} \par\end{flushleft} \subsection{Boundary Conditions} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{topvbc=0}}~\\ \texttt{\small{topvbxval=0.0}}~\\ \texttt{\small{topvbyval=0.0}} & Surface velocity boundary condition parameters. If \texttt{topvbc} is 0, \texttt{topvbxval} and \texttt{topvbyval} specify the tangential surface stress (stress BC). If \texttt{topvbc} is 1, \texttt{topvbxval} and \texttt{topvbyval} specify the tangential surface velocity (velocity BC). The surface normal velocity is zero in these two cases (impermeable BC).\tabularnewline \hline \texttt{\small{botvbc=0}}~\\ \texttt{\small{botvbxval=0.0}}~\\ \texttt{\small{botvbyval=0.0}} & As above, but for bottom velocity boundary conditions.\tabularnewline \hline \texttt{\small{side\_sbcs=off}} & Enable traction boundary condition for the sidewalls or not. Must be \texttt{\small{on}} for coupled model.\tabularnewline \hline \texttt{\small{pseudo\_free\_surf=off}} & Enable pseudo free surface or not.\tabularnewline \hline t\texttt{\small{optbc=1}}~\\ \texttt{\small{toptbcval=0.0 }} & Surface temperature boundary conditions. If \texttt{toptbc} is 0, \texttt{toptbcval} specifies the surface heat flux (not working in this version). If \texttt{toptbc} is 1, the \texttt{toptbcval} specifies the surface temperature. \tabularnewline \hline \texttt{\small{bottbc=1}}~\\ \texttt{\small{bottbcval=1.0}} & As above, but for bottom temperature boundary conditions. \tabularnewline \hline \texttt{\small{temperature\_bound\_adj=off}}~\\ \texttt{\small{depth\_bound\_adj=0.157}}~\\ \texttt{\small{width\_bound\_adj=0.08727}}~\\ \texttt{\small{lith\_age\_depth=0.0314}} & Additional parameters for temperature boundary conditions when \texttt{lith\_age} is \texttt{\small{on}}.\tabularnewline \hline \end{tabular} \subsection{Non-Dimensional Numbers} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{rayleigh=1.0e+5}} & This specifies the Rayleigh number, which is one of the most important parameters you may want to change.\tabularnewline \hline \texttt{\small{dissipation\_number=0.0}} & The dissipation number.\tabularnewline \hline \texttt{\small{gruneisen=0.0}} & The Gruneisen parameter. When this parameter is 0, the code treats it as infinity (i.e., incompressible case). \tabularnewline \hline \texttt{\small{surfaceT=0.1}} & The non-dimensional value of surface temperature.\tabularnewline \hline \texttt{\small{Q0=0.0 }} & This specifies the internal heating number. \tabularnewline \hline \end{tabular} \subsection{Depth Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{z\_lith=0.014}}~\\ \texttt{\small{z\_410=0.06435}}~\\ \texttt{\small{z\_lmantle=0.105}}~\\ \texttt{\small{z\_cmb=0.439}} & These specify the non-dimensional depth of the Moho, 410km discontinuity, 660km discontinuity and D\textquotedbl{}. These parameters are used to determine the depth of viscosity layers and phase changes (see next two sections). The names are only suggestive. You are free to refer \texttt{z\_lith} to an arbitrary depth, for example.\tabularnewline \hline \end{tabular} \subsection{Viscosity\label{sub:Viscosity-input}} \begin{longtable}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Viscosity=system}} & This parameter must be set as indicated.\tabularnewline \hline \texttt{\small{visc\_smooth\_method=3 }} & This specifies which method to use to smooth viscosity projection for the multigrid solver. \tabularnewline \hline \texttt{\small{VISC\_UPDATE=on }} & If \texttt{VISC\_UPDATE} is \texttt{\small{on}}, viscosity will be updated every time step. Otherwise, viscosity will be time-independent.\tabularnewline \hline \texttt{\small{num\_mat=4}} & This specifies the number of material layers. Material 1 is at depth between 0 and \texttt{z\_lith}, material 2 between \texttt{z\_lith} and \texttt{z\_410}, material 3 between \texttt{z\_410} and \texttt{z\_lmantle}, and material 4 between \texttt{z\_lmantle} and the bottom. If \texttt{mat\_control} is on, then a multiplicative factor is applied to the viscosity, as defined below. \tabularnewline \hline \texttt{\small{visc0=1,1,1,1}} & The pre-exponent factor of layered viscosity structure ($\eta_{0}$in the equations below). There must be as many entries as \texttt{num\_mat} in a comma-separated list. \tabularnewline \hline \texttt{\small{TDEPV=off}} & Enable temperature dependence or not.\tabularnewline \hline \texttt{\small{rheol=3}} & When \texttt{rheol}=1, temperature-dependent viscosity is computed by\texttt{\footnotesize{$\eta=\eta_{0}\times\exp(E{}_{\eta}(T_{\eta}-T))$}}~\\ When \texttt{rheol}=2, temperature-dependent viscosity is computed by \texttt{\footnotesize{$\eta=\eta_{0}\times\exp(-T/T_{\eta})$}}~\\ When \texttt{rheol}=3, temperature-dependent viscosity is computed by \texttt{\footnotesize{$\eta=\eta_{0}\times\exp(\frac{E_{\eta}}{T^{*}+T_{\eta}}-\frac{E_{\eta}}{1+T_{\eta}})$}}, where$T^{*}=\min(\max(T,0),1)$.\texttt{\footnotesize{}}~\\ When \texttt{rheol}=4, temperature-dependent viscosity is computed by \texttt{\footnotesize{$\eta=\eta_{0}\times\exp(\frac{E_{\eta}+Z_{\eta}(1-r)}{T^{*}+T_{\eta}})$}}~\\ When \texttt{rheol}=5, same as \texttt{rheol}=3, except the viscosity cut-off is applied before \texttt{mat\_control}.\\ When \texttt{rheol}=6, temperature-dependent viscosity is computed by \texttt{\footnotesize{$\eta=\eta_{0}\times\exp(E_{\eta}(T_{\eta}-T^{*})+Z_{\eta}(1-r))$}}~\\ When \texttt{rheol}=7, temperature-dependent viscosity is computed by \texttt{\footnotesize{$\eta=\eta_{0}\times\exp(\frac{E_{\eta}+Z_{\eta}(1-r)}{T+T_{\eta}}-\frac{E_{\eta}+Z_{\eta}(1-r_{inner})}{1+T_{\eta}})$}}, where$r_{inner}$is the inner radius of the mesh.\\ When \texttt{rheol}=8, same as \texttt{rheol}=3, except viscosity reduction is applied when the temperature exceeds the solidus.\\ When \texttt{rheol}=9, same as \texttt{rheol}=3, except using$T$, not$T^{*}.$\\ When \texttt{rheol}=10, same as \texttt{rheol}=8, except using$T$, not$T^{*}.$\tabularnewline \hline \texttt{\small{viscE=1,1,1,1}}~\\ \texttt{\small{viscT=1,1,1,1}}~\\ \texttt{\small{viscZ=1,1,1,1}} & Parameters defining viscosity law ($E_{\eta}$,$T_{\eta}$, and$Z_{\eta}$in the equations above, respectively). There must be as many entries as \texttt{num\_mat} in a comma-separated list. \tabularnewline \hline \texttt{\small{SDEPV=off}} & Enable stress dependence (non-Newtonian) or not.\tabularnewline \hline \texttt{\small{sdepv\_expt=1,1,1,1}}~\\ \texttt{\small{sdepv\_misfit=0.001}} & If \texttt{\small{SDEPV}} is on, these specify the exponent in the viscosity law and the criterion of convergence test. There must be as many entries as \texttt{\small{num\_mat}} in a comma-separated list.\tabularnewline \hline \texttt{\small{PDEPV=off}} & Pseudo-plastic rheology, implemented by adding a plastic viscosity$\eta_{p}=\frac{\sigma_{y}}{\left(2\epsilon_{II}+10^{-7}\right)+\eta_{p}^{0}}$where the yield stress is defined as$\sigma_{y}=\min\left(a+b\left(1-z\right),y\right)$and$\epsilon_{II}$is the second shear strain rate invariant (all non-dimensional). The effective viscosity$\eta'$is computed from$\eta$and$\eta_{p}$as$\eta'=\frac{\eta_{p}\eta}{\left(\eta+\eta_{p}\right)}$(pdevp\_eff=on) or$\eta'=\min\left(\eta,\eta_{p}\right)$(pdepv\_eff=off).\tabularnewline \hline \texttt{\small{pdepv\_a=1e20,1e20,1e20,1e20}}~\\ \texttt{\small{pdepv\_b=0,0,0,0}}~\\ \texttt{\small{pdepv\_y=1e20,1e20,1e20,1e20}} & Parameters for$\sigma_{y}$.\tabularnewline \hline \texttt{\small{pdepv\_eff=on}} & Effective or minimum viscosity plasticity'' (see above).\tabularnewline \hline \texttt{\small{pdepv\_offset=0}} & Offset for plastic viscosity$\eta_{p}^{0}$.\tabularnewline \hline \texttt{\small{CDEPV=off}} & Compositionally dependent viscosity pre-factor (only for \texttt{\small{tracer=on}}). Will assign two pre-multipliers for$\mathrm{C}=0$and$\mathrm{C}=1$and compute a geometric mean for all elements depending on the nodal composition. Only works for two flavors and assuming C varies between$[0;1]$.\tabularnewline \hline \texttt{\small{cdepv\_ff=1,1}} & Flavor-wise viscosity pre-factors for$\mathrm{C}=0$and$\mathrm{C}=1, respectively.\tabularnewline \hline \texttt{\small{low\_visc\_channel=off}}~\\ \texttt{\small{low\_visc\_wedge=off}} & Used in conjunction with tracers. The tracers define the upper boundary of the subducted slab. \tabularnewline \hline \texttt{\small{lv\_min\_radius=0.9764}}~\\ \texttt{\small{lv\_max\_radius=0.9921}}\\ \texttt{\small{lv\_channel\_thickness=0.0047}}~\\ \texttt{\small{lv\_reduction=0.5}} & If \texttt{\small{low\_visc\_channel}}{\small{ }}or\texttt{\small{ low\_visc\_wedge}} is on, this specifies the radial extents of the low viscosity zones and the viscosity reduction factor, respectively.\tabularnewline \hline \texttt{\small{VMIN=off}}~\\ \texttt{\small{visc\_min=0.001}} & If \texttt{VMIN} is on, minimum viscosity is cut off at \texttt{visc\_min}. \tabularnewline \hline \texttt{\small{VMAX=off}}~\\ \texttt{\small{visc\_max=1000}} & If \texttt{VMAX} is on, maximum viscosity is cut off at \texttt{visc\_max.} \tabularnewline \hline \end{longtable} \subsection{Phase Change Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Ra\_410=0.0}}~\\ \texttt{\small{clapeyron410=0.0235}}~\\ \texttt{\small{transT410=0.78}}~\\ \texttt{\small{width410=0.0058 }} & These specify the phase change parameters (phase change Rayleigh number, Clapeyron slope, ambient temperature, and phase change width, respectively). The depth of this phase change is specified by \texttt{z\_410}.\tabularnewline \hline \texttt{\small{Ra\_670=0.0}}~\\ \texttt{\small{clapeyron670=-0.0235}}~\\ \texttt{\small{transT670=0.875}}~\\ \texttt{\small{width670=0.0058}} & As above. The depth of this phase change is specified by \texttt{z\_lmantle}.\tabularnewline \hline \texttt{\small{Ra\_cmb=0.0}}~\\ \texttt{\small{clapeyroncmb=-0.0235}}~\\ \texttt{\small{transTcmb=0.875}}~\\ \texttt{\small{widthcmb=0.0058}} & As above. The depth of this phase change is specified by \texttt{z\_cmb}.\tabularnewline \hline \end{tabular} \subsection{Momentum Equation Solver Parameters} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{stokes\_flow\_only=off}} & If you wish only to solve for the velocity once (e.g., Stokes flow) then change this parameter to \texttt{\small{on}}. However, if you want to do a convection problem, this should be \texttt{\small{off}}.\tabularnewline \hline \texttt{\small{remove\_rigid\_rotation=on}}~\\ \texttt{\small{remove\_angular\_momentum=on}} & Whether to remove rigid body rotation component or net angular momentum from the velocity solution. Only effective for global models.\tabularnewline \hline \texttt{\small{Solver=cgrad}} & Can be either \texttt{\small{cgrad}} for conjugate gradient solver or \texttt{\small{multigrid}} for multigrid solver for the outer loop of the momentum solver.\tabularnewline \hline \texttt{\small{node\_assemble=on}} & Whether to assemble stiffness matrix at the node level or not.\tabularnewline \hline \texttt{\small{mg\_cycle=1}}~\\ \texttt{\small{down\_heavy=3}}~\\ \texttt{\small{up\_heavy=3}}~\\ \texttt{\small{vlowstep=1000}}~\\ \texttt{\small{vhighstep=3}}~\\ \texttt{\small{max\_mg\_cycles=50}} & Parameters for \texttt{\small{multigrid}} and \texttt{\small{cgrad}} solvers. For \texttt{\small{multigrid}} solver: \texttt{\small{mg\_cycle=1}} for V cycle and 2 for W cycle; \texttt{\small{down\_heavy}} and \texttt{\small{up\_heavy}} are the number of smoothing passes for downward/upward smoothing; \texttt{\small{vlowstep}} and \texttt{\small{vhighstep}} are the number of smoothing passes at lowest/highest levels; \texttt{\small{max\_mg\_cycles}} is the maximum number of multigrid cycles used per solve; all these parameters should be small integers. \\ For \texttt{\small{cgrad}} solver, \texttt{\small{vlowstep}} is the maximum iterations of the conjugate gradient solver and should be a large integer. \tabularnewline \hline \texttt{\small{piterations=1000}} & Maximum iterations of the outer loop for the momentum solver.\tabularnewline \hline \texttt{\small{accuracy=1.0e-4}} & Convergence criterion for the momentum solver. \tabularnewline \hline \texttt{\small{uzawa=cg}}~\\ \texttt{\small{compress\_iter\_maxstep=100}} & Can be either \texttt{\small{cg}} or \texttt{\small{bicg}}. If set to \texttt{\small{cg}}, additional parameters control the maximum iterations.\tabularnewline \hline \texttt{\small{precond=on}} & Whether to use the preconditioner.\tabularnewline \hline \texttt{\small{aug\_lagr=on}}~\\ \texttt{\small{aug\_number=2.0e3}} & Whether to use augmented stiff matrix and the weight of the augmented stiff matrix.\tabularnewline \hline \end{tabular} \subsection{Energy Equation Solver Parameters} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{ADV=on}} & If on, solves the energy equation.\tabularnewline \hline \texttt{\small{fixed\_timestep=0.0}} & If it is equal to 0, the size of the time step is variable and is determined dynamically. Otherwise, the size of the time step is fixed at the specified value.\tabularnewline \hline \texttt{\small{finetunedt=0.9}} & Set the size of the time step to the specified fraction of a maximum stable advection time step. Must be between 0 and 1.\tabularnewline \hline \texttt{\small{adv\_sub\_iterations=2}}~\\ \texttt{\small{adv\_gamma=0.5}} & The number of iterations and the weight of each iteration for the energy predictor-corrector solver.\tabularnewline \hline \texttt{\small{filter\_temp=off}} & Using a Lenardic filter to remove the overshoots and undershoots of the temperature field or not. The filter conserves the total energy.\tabularnewline \hline \texttt{\small{monitor\_max\_T=on}} & If on, the maximum value of the current and previous temperature fields are compared. If the maximum temperature increases more than 5\%, the energy equation solver is rerun with a smaller time-step size. Keep this parameter on to prevent numerical instability.\tabularnewline \hline \texttt{\small{inputdiffusivity=1}} & Currently, don't change this parameter. It is used only in problems which are integrated backward in time.\tabularnewline \hline \end{tabular} \subsection{Age Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{start\_age=40.0}} & Set initial age (in Myrs). This age determines which files of various time-dependent input to read in. \tabularnewline \hline \texttt{\small{reset\_startage=off}} & If \texttt{\small{on}}, the initial age is set to \texttt{start\_age}. If it is \texttt{\small{off}} and \texttt{tic\_method=-1}, the initial age is read in from previous output. \tabularnewline \hline \end{tabular} \subsection{Debugging Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{DESCRIBE=off}}~\\ \texttt{\small{BEGINNER=off }}~\\ \texttt{\small{VERBOSE=off}} & These parameters affect the echo behavior of the CitcomS parser. Only in non-Pyre version.\tabularnewline \hline \texttt{\small{verbose=off }} & This is used for debugging. If verbose is \texttt{\small{on}}, additional information is output to a .\texttt{\small{info}} file. \tabularnewline \hline \texttt{\small{see\_convergence=on}} & If \texttt{\small{on}}, the velocity residual will be output on the screen for every iteration of the momentum solver.\tabularnewline \hline \end{tabular} \subsection{HDF5 Output Parameters} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{cb\_block\_size=1048576}}~\\ \texttt{\small{cb\_buffer\_size=4194304}} & Size for collective buffer in MPI-IO.\tabularnewline \hline \texttt{\small{sieve\_buf\_size=1048576}} & Size of data sieve buffer.\tabularnewline \hline \texttt{\small{output\_alignment=262144}}~\\ \texttt{\small{output\_alignment\_threshold=}}~\\ \texttt{\small{~~~524288}} & Memory alignment.\tabularnewline \hline \texttt{\small{cache\_mdc\_nelmts=10330}}~\\ \texttt{\small{cache\_rdcc\_nelmts=521}}~\\ \texttt{\small{cache\_rdcc\_nbytes=1048576}} & Cache size for chunked dataset.\tabularnewline \hline \end{tabular} \subsection{Tracer Parameters} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{tracer=off}} & If on, enables the tracers. \tabularnewline \hline \texttt{\small{tracer\_ic\_method=0}}~\\ \texttt{\small{tracers\_per\_element=10}}~\\ \texttt{\small{tracer\_file=\textquotedbl{}tracer.dat\textquotedbl{}}} & Specify the initial location of the tracers. If \texttt{\small{tracer\_ic\_method=0}}, the tracers are generated randomly, with the number of tracers per element specified by \texttt{\small{tracers\_per\_element}}. If \texttt{\small{tracer\_ic\_method=1}}, the location of the tracers is read from a file specified in \texttt{\small{tracer\_file}}. If \texttt{\small{tracer\_ic\_method=}}2, the location of the tracers is read from the old tracer output, similar to \texttt{\small{tic\_method=-1}}.\tabularnewline \hline \texttt{\small{tracer\_flavors=0}} & The number of flavors among the tracers. If set to 0, flavor counting is turned off.\tabularnewline \hline \texttt{\small{ic\_method\_for\_flavors=0}}~\\ \texttt{\small{z\_interface=0.7}} & Specify the initial flavors of the tracers. Used only when \texttt{\small{tracer\_ic\_method=0}}. Currently only \texttt{\small{ic\_method\_for\_flavors=0}} is supported. A layered structure is generated. \texttt{\small{z\_interface}} must be a comma-separated list with \texttt{\small{(tracer\_flavors-1)}} entries. Tracers above the 1st radius in\texttt{\small{ z\_interface}} are of flavor 0, between the 1st and 2nd radii are of flavor 1, ... etc.\tabularnewline \hline \texttt{\small{regular\_grid\_deltheta=1.0}}~\\ \texttt{\small{regular\_grid\_delphi=1.0}} & The grid spacing of the regular grid. The regular grid is an auxiliary grid to help locate the tracers.\tabularnewline \hline \texttt{\small{chemical\_buoyancy=on}} & If on, enables thermo-chemical convection.\tabularnewline \hline \texttt{\small{buoy\_type=1}} & If \texttt{\small{buoy\_type=1}}, the composition field is determined by the ratio method \cite{Tackely King tracer ratio method}. (No other methods are implemented yet.) \tabularnewline \hline \texttt{\small{buoyancy\_ratio=1.0}} & The ratio of chemical density anomaly to the reference density, must be a comma-separated list with \texttt{\small{(tracer\_flavors-1)}} entries.\tabularnewline \hline \texttt{\small{itracer\_warnings=on}} & The warning level of the tracer module.\tabularnewline \hline \texttt{\small{tracer\_enriched=off}}~\\ \texttt{\small{Q0\_enriched=0.0}} & Whether the composition anomaly is associated with radioactive heating anomaly. If \texttt{\small{on}}, specifies the internal heating number inside the composition anomaly.\tabularnewline \hline \end{tabular} \subsection{Dimensional Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline r\texttt{\small{adius=6371e3}}~\\ \texttt{\small{density=3340.0}}~\\ \texttt{\small{thermdiff=1.0e-6}}~\\ \texttt{\small{gravacc=9.81}}~\\ \texttt{\small{thermexp=3.0e-5}}~\\ \texttt{\small{refvisc=1.0e+21}}~\\ \texttt{\small{cp=1200}}~\\ \texttt{\small{density\_above=1030.0}}~\\ \texttt{\small{density\_below=6600.0}} & Various dimensional information in SI units.\tabularnewline \hline \end{tabular} \subsection{Required Information} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Problem=convection}}~\\ \texttt{\small{Geometry=sphere}} & For this version of CitcomS, all of these parameters must be set as indicated.\tabularnewline \hline \end{tabular} \subsection{PETSc Related Flags} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{use\_petsc=off}} & Set to \texttt{on} to use PETSc functionality. Turning on this flag only will\\ cause the PETSc version of the Uzawa algorithm to be used in the Stokes solver \tabularnewline\hline \texttt{\small{petsc\_schur=off}} & Set to \texttt{on} to use the Schur complement reduction'' approach in\\ the Stokes solver\tabularnewline\hline \end{tabular} \section{CitcomS Facilities and Properties} This section lists the facilities and properties in the Pyre version of CitcomS. Most of the properties have names which are identical to the parameters used in the non-Pyre version of CitcomS and are explained in the section above. This section highlights those which have changed and those which are entirely new. Parameters are given with their default values. \subsection{Top-Level Facilities and Properties} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{steps=1}} & How many time steps to run.\tabularnewline \hline \texttt{\small{launcher=mpich}} & A facility specifying which launcher to use. Choices include \texttt{\small{mpich}} and \texttt{\small{lam-mpi}}.\tabularnewline \hline \texttt{\small{scheduler=none}} & A facility specifying which scheduler to use. Choices are \texttt{\small{lsf}}, \texttt{\small{pbs}}, and \texttt{\small{globus}}.\tabularnewline \hline \texttt{\small{controller}} & A facility. User cannot change it.\tabularnewline \hline \texttt{\small{solver=regional}} & A facility specifying which solver to use. Must be either \texttt{\small{regional}} or \texttt{\small{full}}.\tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{launcher}}} The launcher facility controls how CitcomS, an MPI application, is executed on multiple processors. It is the equivalent to \texttt{mpirun} in MPI. \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{dry=off}} & If on, prints the \texttt{mpirun} command line and exits without launching the job.\tabularnewline \hline \texttt{\small{nodegen}} & A printf-styled format string, used in conjunction with \texttt{nodelist} to generate a list of machine names. Example: \texttt{\small{n\%03d}}\tabularnewline \hline \texttt{\small{nodelist}} & A comma-separated list of machine names in square brackets. Example: \texttt{ {[}101-103,105,107{]} }\tabularnewline \hline \texttt{\small{command}} & \texttt{The command used to launch the job. The default value should be customized if using any job scheduler, otherwise the default value should work well.}\tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{scheduler}}} The scheduler facility controls how CitcomS submits jobs to a batch scheduler. Available schedulers include: \texttt{none} (default, not using a batch scheduler), \texttt{lsf} (for Load Sharing Facility batch scheduler), \texttt{pbs} (for Portable Batch System batch scheduler), \texttt{sge} (for Sun Grid Engine batch scheduler), and tacc-ranger (a variant of PBS, customized for the Ranger cluster at Texas Advanced Computing Center). Only the command options for these schedulers are listed below. For more information on configuring Pyre for your batch system, see CIG's Pythia page \url{geodynamics.org/cig/software/packages/cs/pythia/docs/batch}. \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{dry=off}} & If \texttt{\small{on}}, prints the batch script and exits without scheduling the job.\tabularnewline \hline \texttt{\small{wait=off}} & If \texttt{\small{on}}, waits for batch jobs to finish before exiting.\tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{job}}} The \texttt{\small{job}} facility controls options for individual batch jobs. \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{name}} & The name of the job.\tabularnewline \hline \texttt{\small{queue}} & The name of the queue to which the job is scheduled.\tabularnewline \hline \texttt{\small{walltime}} & Time limit for the job, specified using units, e.g., \texttt{5{*}minute}, \texttt{2{*}hour}.\tabularnewline \hline \texttt{\small{stdin=/dev/null}} & File to read as the input stream.\tabularnewline \hline \texttt{\small{stdout=stdout.txt}} & File to write for the output stream.\tabularnewline \hline \texttt{\small{stderr=stderr.txt}} & File to write for the error stream.\tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{controller}}} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{monitoringFrequency=100}} & The time-step interval between disk output. It replaces the \texttt{storage\_spacing} parameter in the (old) CitcomS input file.\tabularnewline \hline \texttt{\small{checkpointFrequency=100}} & The time-step interval between checkpoint output, which can be used later to resume the computation.\tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver }}} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{mesher }} & A facility. Must be either \texttt{\small{full-sphere}} or \texttt{\small{regional-sphere}}. User does not need to specify it. It is set by the \texttt{solver} facility automatically. \tabularnewline \hline \texttt{\small{tsolver }} & A facility for the temperature solver using Petrov-Galerkin time integration. \tabularnewline \hline \texttt{\small{vsolver }} & A facility for the velocity solver using Boussinesq approximation and the Uzawa algorithm. \tabularnewline \hline \texttt{\small{bc }} & A facility of boundary conditions.\tabularnewline \hline \texttt{\small{const}} & A facility for dimensional constants.\tabularnewline \hline \texttt{\small{ic}} & A facility for initial conditions.\tabularnewline \hline \texttt{\small{output}} & A facility for output options.\tabularnewline \hline \texttt{\small{param}} & A facility for some additional input parameter files.\tabularnewline \hline \texttt{\small{phase }} & A facility for phase change parameters.\tabularnewline \hline \texttt{\small{tracer }} & A facility for tracer parameters.\tabularnewline \hline \texttt{\small{visc}} & A facility for viscosity parameters.\tabularnewline \hline \texttt{\small{datadir=\textquotedbl{}.\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{datafile=\textquotedbl{}regtest\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{datadir\_old=\textquotedbl{}.\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{datafile\_old=\textquotedbl{}regtest\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{rayleigh=100000}} & \tabularnewline \hline \texttt{\small{dissipation\_number=0.0}} & \tabularnewline \hline \texttt{\small{gruneisen=0.0}} & \tabularnewline \hline \texttt{\small{surfaceT=0.1}} & \tabularnewline \hline \texttt{\small{Q0=0}} & \tabularnewline \hline \texttt{\small{stokes\_flow\_only=off}} & \tabularnewline \hline \texttt{\small{verbose=off}} & \tabularnewline \hline \texttt{\small{see\_convergence=on}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.mesher}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{nprocx=1}} & \tabularnewline \hline \texttt{\small{nprocy=1}} & \tabularnewline \hline \texttt{\small{nprocz=1 }} & \tabularnewline \hline \texttt{\small{coor=0}} & \tabularnewline \hline \texttt{\small{coor\_file=\textquotedbl{}coor.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{coor\_refine=0.1,0.15,0.1,0.2}} & \tabularnewline \hline \texttt{\small{nodex=9}} & \tabularnewline \hline \texttt{\small{nodey=9}} & \tabularnewline \hline \texttt{\small{nodez=9}} & \tabularnewline \hline \texttt{\small{levels=1 }} & \tabularnewline \hline \texttt{\small{radius\_outer=1}} & \tabularnewline \hline \texttt{\small{radius\_inner=0.55}} & \tabularnewline \hline \texttt{\small{theta\_min=1.0708}} & \tabularnewline \hline \texttt{\small{theta\_max=2.0708}} & \tabularnewline \hline \texttt{\small{fi\_min=0 }} & \tabularnewline \hline \texttt{\small{fi\_max=1}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.tsolver}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{ADV=on}} & \tabularnewline \hline \texttt{\small{filter\_temp=off}} & \tabularnewline \hline \texttt{\small{monitor\_max\_T=on}} & \tabularnewline \hline \texttt{\small{finetunedt=0.9}} & \tabularnewline \hline \texttt{\small{fixed\_timestep=0.0}} & \tabularnewline \hline \texttt{\small{inputdiffusivity=1}} & \tabularnewline \hline \texttt{\small{adv\_gamma=0.5}} & \tabularnewline \hline \texttt{\small{adv\_sub\_iterations=2}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.vsolver}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Solver=cgrad}} & \tabularnewline \hline \texttt{\small{node\_assemble=on}} & \tabularnewline \hline \texttt{\small{precond=on}} & \tabularnewline \hline \texttt{\small{accuracy=1.0e-4}} & \tabularnewline \hline \texttt{\small{uzawa=cg}} & \tabularnewline \hline \texttt{\small{compress\_iter\_maxstep=100}} & \tabularnewline \hline \texttt{\small{mg\_cycle=1 }} & \tabularnewline \hline \texttt{\small{down\_heavy=3}} & \tabularnewline \hline \texttt{\small{up\_heavy=3}} & \tabularnewline \hline \texttt{\small{vlowstep=1000}} & \tabularnewline \hline \texttt{\small{vhighstep=3}} & \tabularnewline \hline \texttt{\small{max\_mg\_cycles=50}} & \tabularnewline \hline \texttt{\small{piterations=1000}} & \tabularnewline \hline \texttt{\small{aug\_lagr=on}} & \tabularnewline \hline \texttt{\small{aug\_number=2000}} & \tabularnewline \hline \texttt{\small{remove\_rigid\_rotation=on}} & \tabularnewline \hline \texttt{\small{remove\_angular\_momentum=on}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.bc }}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{side\_sbcs=off}} & \tabularnewline \hline \texttt{\small{pseudo\_free\_surf=off}} & \tabularnewline \hline \texttt{\small{topvbc=0}} & \tabularnewline \hline \texttt{\small{topvbxval=0}} & \tabularnewline \hline \texttt{\small{topvbyval=0}} & \tabularnewline \hline \texttt{\small{botvbc=0}} & \tabularnewline \hline \texttt{\small{botvbxval=0}} & \tabularnewline \hline \texttt{\small{botvbyval=0}} & \tabularnewline \hline \texttt{\small{toptbc=1}} & \tabularnewline \hline \texttt{\small{toptbcval=0}} & \tabularnewline \hline \texttt{\small{bottbc=1}} & \tabularnewline \hline \texttt{\small{bottbcval=1}} & \tabularnewline \hline \texttt{\small{temperature\_bound\_adj=off}} & \tabularnewline \hline \texttt{\small{depth\_bound\_adj=0.157}} & \tabularnewline \hline \texttt{\small{width\_bound\_adj=0.08727}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.const }}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{radius=6.371e+06}} & \tabularnewline \hline \texttt{\small{density=3340.0}} & \tabularnewline \hline \texttt{\small{thermdiff=1e-06}} & \tabularnewline \hline \texttt{\small{gravacc=9.81}} & \tabularnewline \hline \texttt{\small{thermexp=3e-05}} & \tabularnewline \hline \texttt{\small{refvisc=1e+21}} & \tabularnewline \hline \texttt{\small{cp=1200}} & \tabularnewline \hline \texttt{\small{density\_above=1030.0}} & \tabularnewline \hline \texttt{\small{density\_below=6600.0}} & \tabularnewline \hline \texttt{\small{z\_lith=0.014}} & \tabularnewline \hline \texttt{\small{z\_410=0.06435}} & \tabularnewline \hline \texttt{\small{z\_lmantle=0.105}} & \tabularnewline \hline \texttt{\small{z\_cmb=0.439}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.ic }}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{restart=off}} & \tabularnewline \hline \texttt{\small{post\_p=off}} & \tabularnewline \hline \texttt{\small{solution\_cycles\_init=0}} & \tabularnewline \hline \texttt{\small{zero\_elapsed\_time=on}} & \tabularnewline \hline \texttt{\small{tic\_method=0}} & \tabularnewline \hline \texttt{\small{num\_perturbations=1}} & \tabularnewline \hline \texttt{\small{perturbl=1}} & \tabularnewline \hline \texttt{\small{perturbm=1}} & \tabularnewline \hline \texttt{\small{perturblayer=5}} & \tabularnewline \hline \texttt{\small{perturbmag=0.05}} & \tabularnewline \hline \texttt{\small{half\_space\_age=40}} & \tabularnewline \hline \texttt{\small{mantle\_temp=1.0}} & \tabularnewline \hline \texttt{\small{blob\_center={[}-999,-999,-999{]} }} & \tabularnewline \hline \texttt{\small{blob\_radius=0.063 }} & \tabularnewline \hline \texttt{\small{blob\_dT=0.18}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.output}}} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{output\_format=\textquotedbl{}ascii\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{output\_optional=\textquotedbl{}surf,botm,}}~\\ \texttt{\small{tracer\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{output\_ll\_max=20}} & \tabularnewline \hline \texttt{self\_gravitation=off} & \tabularnewline \hline \texttt{use\_cbf\_topo=off} & \tabularnewline \hline \texttt{\small{cb\_block\_size=1048576}} & \tabularnewline \hline \texttt{\small{cb\_buffer\_size=4194304}} & \tabularnewline \hline \texttt{\small{sieve\_buf\_size=1048576}} & \tabularnewline \hline \texttt{\small{output\_alignment=262144}} & \tabularnewline \hline \texttt{\small{output\_alignment\_threshold=}}~\\ \texttt{\small{~~~~524288}} & \tabularnewline \hline \texttt{\small{cache\_mdc\_nelmts=10330}} & \tabularnewline \hline \texttt{\small{cache\_rdcc\_nelmts=521}} & \tabularnewline \hline \texttt{\small{cache\_rdcc\_nbytes=1048576}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.param}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{reference\_state=1}} & \tabularnewline \hline \texttt{\small{refstate\_file=\textquotedbl{}refstate.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{mineral\_physics\_model=3}} & \tabularnewline \hline \texttt{\small{file\_vbcs=off}} & \tabularnewline \hline \texttt{\small{vel\_bound\_file=\textquotedbl{}bvel.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{mat\_control=off}} & \tabularnewline \hline \texttt{\small{mat\_file=\textquotedbl{}mat.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{lith\_age=off}} & \tabularnewline \hline \texttt{\small{lith\_age\_file=\textquotedbl{}age.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{lith\_age\_time=off}} & \tabularnewline \hline \texttt{\small{lith\_age\_depth=0.0314}} & \tabularnewline \hline \texttt{\small{start\_age=40 }} & \tabularnewline \hline \texttt{\small{reset\_startage=off}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.phase}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Ra\_410=0 }} & \tabularnewline \hline \texttt{\small{clapeyron410=0.0235}} & \tabularnewline \hline \texttt{\small{transT410=0.78 }} & \tabularnewline \hline \texttt{\small{width410=0.0058}} & \tabularnewline \hline \texttt{\small{Ra\_670=0}} & \tabularnewline \hline \texttt{\small{clapeyron670=-0.0235}} & \tabularnewline \hline \texttt{\small{transT670=0.78}} & \tabularnewline \hline \texttt{\small{width670=0.0058}} & \tabularnewline \hline \texttt{\small{Ra\_cmb=0}} & \tabularnewline \hline \texttt{\small{clapeyroncmb=-0.0235}} & \tabularnewline \hline \texttt{\small{transTcmb=0.875}} & \tabularnewline \hline \texttt{\small{widthcmb=0.0058}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.tracer}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{tracer=off}} & \tabularnewline \hline \texttt{\small{tracer\_ic\_method=0}} & \tabularnewline \hline \texttt{\small{tracers\_per\_element=10}} & \tabularnewline \hline \texttt{\small{tracer\_file=\textquotedbl{}tracer.dat\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{tracer\_flavors=0}} & \tabularnewline \hline \texttt{\small{ic\_method\_for\_flavors=0}} & \tabularnewline \hline \texttt{\small{z\_interface=0.7}} & \tabularnewline \hline \texttt{\small{itracer\_warnings=on}} & \tabularnewline \hline \texttt{\small{regular\_grid\_deltheta=1.0}} & \tabularnewline \hline \texttt{\small{regular\_grid\_delphi=1.0}} & \tabularnewline \hline \texttt{\small{chemical\_buoyancy=on}} & \tabularnewline \hline \texttt{\small{buoy\_type=1}} & \tabularnewline \hline \texttt{\small{buoyancy\_ratio=1.0}} & \tabularnewline \hline \texttt{\small{tracer\_enriched=off}} & \tabularnewline \hline \texttt{\small{Q0\_enriched=0.0}} & \tabularnewline \hline \end{tabular} \subsection{\texttt{\normalsize{solver.visc}}} \noindent % \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{Viscosity=\textquotedbl{}system\textquotedbl{}}} & \tabularnewline \hline \texttt{\small{visc\_smooth\_method=3}} & \tabularnewline \hline \texttt{\small{VISC\_UPDATE=on}} & \tabularnewline \hline \texttt{\small{num\_mat=4}} & \tabularnewline \hline \texttt{\small{visc0=1,1,1,1,}} & \tabularnewline \hline \texttt{\small{TDEPV=off}} & \tabularnewline \hline \texttt{\small{rheol=3}} & \tabularnewline \hline \texttt{\small{viscE=1,1,1,1,}} & \tabularnewline \hline \texttt{\small{viscT=1,1,1,1,}} & \tabularnewline \hline \texttt{\small{viscZ=1,1,1,1,}} & \tabularnewline \hline \texttt{\small{SDEPV=off}} & \tabularnewline \hline \texttt{\small{sdepv\_misfit=0.001}} & \tabularnewline \hline \texttt{\small{sdepv\_expt=1,1,1,1,}} & \tabularnewline \hline \texttt{\small{PDEPV=off}} & \tabularnewline \hline \texttt{\small{pdepv\_a=1e20,1e20,1e20,1e20}} & \tabularnewline \hline \texttt{\small{pdepv\_b=0,0,0,0}} & \tabularnewline \hline \texttt{\small{pdepv\_y=1e20,1e20,1e20,1e20}} & \tabularnewline \hline \texttt{\small{pdepv\_eff=on}} & \tabularnewline \hline \texttt{\small{pdepv\_offset=0}} & \tabularnewline \hline \texttt{\small{CDEPV=off}} & \tabularnewline \hline \texttt{\small{cdepv\_ff=1,1}} & \tabularnewline \hline \texttt{\small{low\_visc\_channel=0}} & \tabularnewline \hline \texttt{\small{low\_visc\_wedge=0}} & \tabularnewline \hline \texttt{\small{lv\_min\_radius=0.9764}} & \tabularnewline \hline \texttt{\small{lv\_max\_radius=0.9921}} & \tabularnewline \hline \texttt{\small{lv\_channel\_thickness=0.0047}} & \tabularnewline \hline \texttt{\small{lv\_reduction=0.5}} & \tabularnewline \hline \texttt{\small{VMIN=off}} & \tabularnewline \hline \texttt{\small{visc\_min=0.001}} & \tabularnewline \hline \texttt{\small{VMAX=off}} & \tabularnewline \hline \texttt{\small{visc\_max=1000}} & \tabularnewline \hline \end{tabular} \subsection{\noindent \texttt{\normalsize{layout}}} The following facilities are only used in the coupled CitcomS solver. \\ \\ \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{containing\_group={[}0-11{]}}} & The processors to be assigned to the csolver.\tabularnewline \hline \texttt{\small{embedded\_group={[}12{]}}} & The processors to be assigned to the esolver.\tabularnewline \hline \end{tabular} \subsection{\noindent \texttt{\normalsize{ccoupler}}} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{si\_unit=off}} & Whether to exchange information in SI units.\tabularnewline \hline \texttt{\small{cartesian=off}} & Whether to exchange information in Cartesian coordinate.\tabularnewline \hline \texttt{\small{two\_way\_communication=on}} & Whether to transfer information from esolver to csolver.\tabularnewline \hline \texttt{\small{exchange\_initial\_temperature}}~\\ \texttt{\small{=on}} & Whether to ensure consistent initial temperature between csolver or esolver.\tabularnewline \hline \texttt{\small{exchange\_pressure=on}} & Whether to send the pressure solution to esolver.\tabularnewline \hline \end{tabular} \subsection{\noindent \texttt{\normalsize{ecoupler}}} \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{si\_unit=off}} & Whether to exchange information in SI units.\tabularnewline \hline \texttt{\small{cartesian=off}} & Whether to exchange information in Cartesian coordinate.\tabularnewline \hline \texttt{\small{two\_way\_communication=on}} & Whether to transfer information from esolver to csolver.\tabularnewline \hline \texttt{\small{exchange\_initial\_temperature=on}} & Whether to ensure consistent initial temperature between csolver or esolver.\tabularnewline \hline \texttt{\small{exchange\_pressure=on}} & Whether to receive the initial guess of pressure from csolver.\tabularnewline \hline \texttt{\small{exclude\_top=off}} & Whether to exclude the top boundary when receiving boundary conditions from the csolver. If the esolver has imposed plate velocity and constant temperature at the surface, set this to on.\tabularnewline \hline \texttt{\small{exclude\_bottom=off}} & Whether to exclude the bottom boundary when receiving boundary conditions from the csolver.\tabularnewline \hline \texttt{\small{amending\_outflow=off}} & Whether to amend the velocity boundary conditions to satisfy the continuity equation.\tabularnewline \hline \end{tabular} \subsection{\noindent \texttt{\normalsize{journal}}} The Pyre facility \texttt{journal} provides five types of debugging output for conceptually different types of information. The journal output stream can be instantiated as separated channels multiple times, and each channel can be individually activated or deactivated, and sent to different locations (the terminal, sockets, files, devices, etc.). The streams and their default states of the journal streams are:\\ \\ \begin{tabular}{|>{\raggedright}p{1.85in}|>{\raggedright}p{4.25in}|} \hline \texttt{\small{debug}} & Debugging information. Default off.\tabularnewline \hline \texttt{\small{error}} & Unrecoverable runtime error. Default on. \tabularnewline \hline \texttt{\small{firewall}} & Fatal programming error. Default on. \tabularnewline \hline \texttt{\small{info}} & Descriptive information. Default off.\tabularnewline \hline \texttt{\small{warning}} & Recoverable runtime error. Default off.\tabularnewline \hline \end{tabular}\\ \\ For coupled solvers, there is a large amount of debugging information that outputs through the journal facility. That output can be turned on/off with command line options. The most important ones are \texttt{journal.debug.Exchanger} and \texttt{journal.debug.CitcomS-Exchanger}. Other less important options can be found in \texttt{tests/coupled.cfg}. To enable the journal output, try running the scripts with the following options: \begin{lyxcode} {[}CitcomS.journal.debug{]}~\\ Exchanger~=~on~\\ CitcomS-Exchanger~=~on \end{lyxcode} \chapter{CitcomS Input File Format} \section{Introduction} CitcomS expects Unix-styled ASCII files (i.e., no carriage character following new line character) for all input files. This can be a nuisance in DOS/Windows systems. You may want to find a text editor that can write Unix-style ASCII files. In the following, words in normal \texttt{courier} must be input exactly as shown, while \texttt{\textbf{\emph{italicized}}} words should be substituted by your values. All parameters are in non-dimensional units unless specified. \section{\label{sec:Coordinate-files}Coordinate Files} For regional version of CitcomS, the mesh must be regular, but the mesh spacing may be unequal. The \texttt{coor\_file} has the format: \begin{lyxcode} nsd=~1~\\ 1~\textbf{~~~~}\textbf{\emph{theta1}}~\\ 2~\textbf{~~~~}\textbf{\emph{theta2}}~\\ ...~~~...\textbf{}~\\ \textbf{\emph{nodex~theta\_nodex}}~\\ nsd=~2~\\ 1~\textbf{~~~~}\textbf{\emph{phi1}}~\\ 2~\textbf{~~~~}\textbf{\emph{phi2}}~\\ ...~~~...\textbf{}~\\ \textbf{\emph{nodey~phi\_nodey}}~\\ nsd=~3~\\ 1~\textbf{~~~~}\textbf{\emph{r1}}~\\ 2~\textbf{~~~~}\textbf{\emph{r2}}~\\ ...~~~...\textbf{}~\\ \textbf{\emph{nodez~r\_nodez}} \end{lyxcode} For the full spherical version of CitcomS, the mesh of each cap must be regular and equidistant in the horizontal dimension. Only the vertical dimension is specified by \texttt{coor\_file}. The \texttt{coor\_file} has the format: \begin{lyxcode} nsd=~3~\\ 1~\textbf{~~~~}\textbf{\emph{r1}}~\\ 2~\textbf{~~~~}\textbf{\emph{r2}}~\\ ...~~~...\textbf{}~\\ \textbf{\emph{nodez~r\_nodez}} \end{lyxcode} \section{\label{sec:Velocity-boundary-condition}Velocity Boundary Condition Files} If \texttt{vel\_bound\_file} is set to \texttt{bvel.dat}, then it is necessary to have one \texttt{bvel.dat} file per cap for each million-year interval. For example, if \texttt{start\_age=83}, then the following files are needed for the regional mesh: \begin{lyxcode} bvel.dat84~\\ bvel.dat83~\\ bvel.dat82~\\ ...~\\ bvel.dat0\end{lyxcode} \begin{quote} \textbf{TIP:} Even though the model starts at 83 million years before the present, by default it will attempt to read in both a file for 84 and a file for 83 million years ago. To deal with this, create a duplicate copy of \texttt{bvel.dat83} and call it \texttt{bvel.dat84} \end{quote} In this example, the model age will start at 83 Ma and decrease toward 0 Ma. At any particular age, the velocity boundary conditions are interpolated in time according to the nearest \texttt{bvel.dat} files. Once the age becomes negative, the velocity boundary conditions are specified by the content of \texttt{bvel.dat0.} \begin{quote} \textbf{TIP:} By setting \texttt{start\_age=0} and providing identical \texttt{bvel.dat0} and \texttt{bvel.dat1}, you can effectively get time-invariant velocity boundary conditions. \end{quote} For the global mesh, each of the 12 caps requires one file for each million-year interval. For example, these files are needed for an 82-million-year age: \begin{lyxcode} bvel.dat82.0~\\ bvel.dat82.1~\\ bvel.dat82.2~\\ ...~\\ bvel.dat82.11 \end{lyxcode} Each velocity boundary condition file has the format: \begin{lyxcode} \textbf{\emph{Vx~Vy}} \end{lyxcode} The units of \texttt{Vx} and \texttt{Vy} are cm/yr. \section{Material Files} In this version of CitcomS, the implementation of material support is not working. The material output has been disabled and is documented here for completion. If \texttt{mat\_file} is set to \texttt{mat.dat}, then it is necessary to have one \texttt{\small{mat.dat}} file for each million-year interval. For example, if \texttt{start\_age=83}, then the following files are needed: \begin{lyxcode} mat.dat84~\\ mat.dat83~\\ mat.dat82~\\ ...~\\ mat.dat0 \end{lyxcode} The same note about \texttt{vel\_bound\_file} (see Section \ref{sec:Velocity-boundary-condition} above) applies. The format of \texttt{mat\_file} is: \begin{lyxcode} \textbf{\emph{n~viscosity\_factor}} \end{lyxcode} \section{Lithosphere Age Files} If \texttt{lith\_age\_file} is set to \texttt{lith.dat}, then it is necessary to have one \texttt{bvel.dat} file for each million-year interval. For example, if \texttt{start\_age=83}, then the following files are needed for the regional mesh: \begin{lyxcode} lith.dat84~\\ lith.dat83~\\ lith.dat82~\\ ...~\\ lith.dat0 \end{lyxcode} The same note about \texttt{vel\_bound\_file} (see Section \ref{sec:Velocity-boundary-condition} above) still applies. The input age is in millions of years. For the global mesh, each of the 12 caps requires one file for each million-year interval. For example, these files are needed for an 82-million-year age: \begin{lyxcode} lith.dat82.0~\\ lith.dat82.1~\\ lith.dat82.2~\\ ...~\\ lith.dat82.11 \end{lyxcode} The format of \texttt{lith\_age\_file} is: \begin{lyxcode} \textbf{\emph{n~age}} \end{lyxcode} \section{Tracer Files} This file contains the initial location of all tracers. The first line is the number of tracers and the number of columns in the file. The first three columns are the coordinates of the tracers. The fourth column, if present, indicates the flavor of the tracers. \begin{lyxcode} \textbf{\emph{num\_tracers~num\_columns}}~\\ \textbf{\emph{theta~phi~radius~{[}flavor{]}}} \end{lyxcode} \section{Reference State Files} This file contains the profiles of the reference state. This file must contain at least \texttt{\small{nodez}} lines, and each line must contain 7 columns, where the meaning of each column is: \begin{lyxcode} \textbf{\rho_{r}$~$g$~$\alpha$~$c_{P}$~}\textbf{\emph{reserved~reserved~reserved}} \end{lyxcode} \chapter{\label{cha:Appendix-C:-CitComS,}CitcomS Output File Format} \section{Introduction} The format of the output files of CitcomS is described here. In the following sections, the model prefix is assumed as \texttt{test-case}, the processor number as 0, and the time step as 10. All outputs are in non-dimensional units unless specified. \section{Postprocessed Cap Output (\texttt{\large{test-case.cap00.10}})} The command \texttt{autocombine.py} produces 1 cap file for regional CitcomS and 12 cap files for the full CitcomS. The first line of the cap file is a comment describing the node geometry (\texttt{nodex}$\times$\texttt{nodey}$\times$\texttt{nodez}). The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode} colatitude~longitude~radius~vel\_colat~vel\_lon~vel\_r~temperature~viscosity \end{lyxcode} \section{Postprocessed Opt Output (\texttt{\large{test-case.opt00.10}})} If the optional output contains node-based fields (e.g., \texttt{comp\_nd}, \texttt{pressure}, and \texttt{stress}), the command \texttt{autocombine.py} produces 1 opt file for regional CitcomS and 12 opt files for the full CitcomS. The first line of the opt file is a comment describing the node geometry (\texttt{nodex}$\times$\texttt{nodey}$\times$\texttt{nodez}). The rest of the file is column-based, where the meaning of each column is defined in the corresponding \texttt{.general} file. For example, if the \texttt{.general} file contains these two lines, then the first column in the opt file is the composition, the 2nd column is the pressure, and the last six columns are the stress. \begin{lyxcode} field~=~composition,~pressure,~stress~\\ structure~=~scalar,~scalar,~6-vector \end{lyxcode} \section{Postprocessed Surf Output (\texttt{\large{test-case.surf0.10}})} An undocumented (and unmaintained) post-processing script \texttt{batchsurf.py} can combine the \texttt{coord}, \texttt{surf}, and \texttt{botm} output. The format of the combined surf file is listed for completion. The first line of the surf file is a comment describing the node geometry (\texttt{nodex}$\times$\texttt{nodey}). The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode} colatitude~longitude~topography~heat\_flux~vel\_colat~vel\_lon~ \end{lyxcode} \section{Time Output (\texttt{\large{test-case.time}})} This file reports non-dimensional elapsed model time and spent CPU time (in seconds), and it is only outputted on computer node 0. The meaning of each column is: \begin{lyxcode} step~total\_t~delta\_t~total\_cpu\_time~step\_cpu\_time~ \end{lyxcode} \section{\label{sec:ASCII-Output}ASCII Output} \subsection{Coordinate Output (\texttt{\normalsize{test-case.coord.0}})} This file is only outputted at the 0th time step. The first line is a header. The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode} colatitude~longitude~radius \end{lyxcode} \subsection{\label{sub:Velocity-and-Temperature}Velocity and Temperature Output (\texttt{\normalsize{test-case.velo.0.10}})} The first two lines of this file are headers. The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode} vel\_colat~vel\_lon~vel\_r~temperature \end{lyxcode} \subsection{Viscosity Output (\texttt{\normalsize{test-case.visc.0.10}})} The first line of this file is a header. The rest of the file is column-based, where the meaning of the only column is: \begin{lyxcode} viscosity \end{lyxcode} \subsection{Material Output (\texttt{\normalsize{test-case.mat.0}})} In this version of CitcomS, the implementation of material support is not working. The material output has been disabled and is documented here for completion. This file is only outputted at the 0th time step. There is no header. The meaning of each column is: \begin{lyxcode} element\_number~layer~material \end{lyxcode} \subsection{Surface Variables Output (\texttt{\normalsize{test-case.surf.0.10}} and \texttt{\normalsize{test-case.botm.0.10}})} The first line of each file is a header. The rest of each file is column-based, where the meaning of each column is: \begin{lyxcode} topography~heat\_flux~vel\_colat~vel\_lon \end{lyxcode} \subsection{Stress Output (\texttt{\normalsize{test-case.stress.0.10}})} The first two lines of the file are headers. The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode}$S_{\theta\theta}$~$S_{\phi\phi}$~$S_{rr}$~$S_{\theta\phi}$~$S_{\theta r}$~$S_{\phi r}$~ \end{lyxcode} You can use the above values to form the following symmetric stress tensor: $$\left[\begin{array}{ccc} S_{\theta\theta} & S_{\theta\phi} & S_{\theta r}\\ S_{\phi\theta} & S_{\phi\phi} & S_{\phi r}\\ S_{r\theta} & S_{r\phi} & S_{rr} \end{array}\right]\label{eq:symmetric stress tensor}$$ \subsection{Pressure Output (\texttt{\normalsize{test-case.pressure.0.10}})} The first line of the file is a header. The rest of the file is column-based, where the meaning of the only column is: \begin{lyxcode} pressure \end{lyxcode} \subsection{Horizontal Average Output (\texttt{\normalsize{test-case.horiz\_avg.0.10}})} The first line of the file is a header. The rest of the file is column-based, where the meaning of each column is: \begin{lyxcode} radius~temperature~RMS(V\_horizontal)~RMS(V\_vertical) \end{lyxcode} \subsection{\label{sub:Geoid-Output-(test-case.geoid.10)}Geoid Output (\texttt{\normalsize{test-case.geoid.10}})} The first line of the file is a header. The rest of the file is column-based, where the meaning of each column is: \noindent \begin{center} \texttt{\footnotesize{}} \begin{table}[H] \noindent \centering{}% \begin{tabular}{|c|l|} \hline col 1 & spherical harmonic degrees (\texttt{l})\tabularnewline \hline col 2 & spherical harmonic orders (m)\tabularnewline \hline col 3 & cosine coefficients of total geoid (S$_{\mathrm{lm}}$)\tabularnewline \hline col 4 & sine coefficients of total geoid (C$_{\mathrm{lm}}$)\tabularnewline \hline col 5 & cosine coefficients of the geoid due to surface topography\tabularnewline \hline col 6 & sine coefficents of the geoid due to surface topography\tabularnewline \hline col 7 & cosine coefficients of the geoid due to internal density variation\tabularnewline \hline col 8 & sine coefficients of the geoid due to internal density variation\tabularnewline \hline \end{tabular} \end{table} \par\end{center}{\footnotesize \par} The units of the geoid coefficients are meters. The geoid field can be reconstructed by ${\displaystyle \sum}\left(S_{lm}\sin(m\phi)+C_{lm}\cos(m\phi)\right)P_{lm}(\cos\theta)$ where$P_{\mathrm{lm}}$is the associated Legendre polynomials. \subsection{Tracer Output (\texttt{\normalsize{test-case.tracer.0.10}})} The first line of the file is a header. The second field in the header is the number of tracers in the file. The third field in the header is the number of columns in the file. The rest of the file is column-based, where the meaning of the columns is: \begin{lyxcode} theta~phi~radius~{[}flavor{]} \end{lyxcode} \subsection{Composition Output (\texttt{\normalsize{test-case.comp\_el.0.10}} and \texttt{\normalsize{test-case.comp\_nd.0.10}})} These files contain the composition field either on the element (\texttt{comp\_el}) or on the node (\texttt{comp\_nd}). The format of the files is the same. The first two lines of the file are headers. The fourth field in the first line is the number of composition components. In the second line, the first field is the current bulk sum of composition 1. The second field is the initial bulk sum of composition 1. The third and fourth fields (if they exist) are the current and initial bulk sum of composition 2, etc. The rest of the file is column-based, where the meaning of the columns is: \begin{lyxcode} composition1~{[}composition2~...~compositionN{]} \end{lyxcode} \subsection{Heating Output (\texttt{\normalsize{test-case.heating.0.10}})} These files contain various heating terms on the element. The first two lines of the file are the headers. The rest of the file is column-based, where the meaning of the columns is: \begin{lyxcode} adiabatic\_heating~viscous\_heating~latent\_heating \end{lyxcode} \section{HDF5 Output (\texttt{\large{test-case.h5}})} The format and layout of HDF5 output is described in Section \vref{sec:Data-Layout}. \section{Misc. Binary Output\label{sec:Misc.-Binary-Output}} \subsection{Checkpoint Output (\texttt{\normalsize{test-case.checkpoint.0.10}})} These files are used for restarting the run. These files are not portable. You cannot run a model on machine A, copy the checkpoint files to machine B, and expect you can always restart the simulation on machine B. Their format is undocumented on purpose and will remain so. \subsection{Domain Output (\texttt{\normalsize{test-case.domain}})} This file is only outputted at the 0th time step. This file is output by rank-0 processor only. It contains the domain bounds of all processors in binary format. The file begins with four intergers: \texttt{nproc}, \texttt{ncolumns}, \texttt{const1}, \texttt{const2}, where \texttt{nproc} is the number of processors, \texttt{ncolumns} is 10, \texttt{const1} and \texttt{const2} are used to detect binary incompatibility. The rest of the file has \texttt{nproc}$\times$\texttt{ncolumns} doubles, the columns contains the min/max of radius and the$\theta,\phi$coordinates of the four corners: \begin{lyxcode} r\_min~r\_max~theta0~phi0~theta1~phi1~theta2~phi2~theta3~phi3 \end{lyxcode} \subsection{Coordinate Binary Output (\texttt{\normalsize{test-case.coord\_bin.0}})} These files are only outputted at the 0th time step. These files contain the coordinates of nodes in binary format. The file begins with four intergers: \texttt{nodex}, \texttt{nodey}, \texttt{nodez}, \texttt{const}, then followed by the$\theta,\phi,r\$ coordinates of all nodes.

\subsection{Seismic Output (\texttt{\normalsize{test-case.seismic.0.10}})}

These files contain the density, P wave velocity and S wave velocity
in binary format. The density of all nodes is stored in double, then
followed by P wave velocity of all nodes, and finally S wave velocity
of all nodes.

(C) 1989, 1991 Free Software Foundation, Inc. 59 Temple Place, Suite
330, Boston, MA 02111-1307 USA} \\
Everyone is permitted to copy and distribute verbatim copies of this
license document, but changing it is not allowed.

\section*{Preamble}

The licenses for most software are designed to take away your freedom
to share and change it. By contrast, the GNU General Public License
is intended to guarantee your freedom to share and change free software
-- to make sure the software is free for all its users. This General
Public License applies to most of the Free Software Foundation's software
and to any other program whose authors commit to using it. (Some other
Free Software Foundation software is covered by the GNU Library General

When we speak of free software, we are referring to freedom, not price.
Our General Public Licenses are designed to make sure that you have
the freedom to distribute copies of free software (and charge for
this service if you wish), that you receive source code or can get
it if you want it, that you can change the software or use pieces
of it in new free programs; and that you know you can do these things.

To protect your rights, we need to make restrictions that forbid anyone
to deny you these rights or to ask you to surrender the rights. These
restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.

For example, if you distribute copies of such a program, whether gratis
or for a fee, you must give the recipients all the rights that you
have. You must make sure that they, too, receive or can get the source
code. And you must show them these terms so they know their rights.

We protect your rights with two steps:
\begin{enumerate}
\item Offer you this license which gives you legal permission to copy, distribute
and/or modify the software.
\end{enumerate}
Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software. If the software is modified by someone else and passed on,
we want its recipients to know that what they have is not the original,
so that any problems introduced by others will not reflect on the
original authors' reputations.

Finally, any free program is threatened constantly by software patents.
We wish to avoid the danger that redistributors of a free program
will individually obtain patent licenses, in effect making the program
proprietary. To prevent this, we have made it clear that any patent
must be licensed for everyone's free use or not licensed at all.

The precise terms and conditions for copying, distribution and modification
follow.

\section*{GNU GENERAL PUBLIC LICENSE TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION
AND MODIFICATION }

\begin{itemize}

\item[0.]This License applies to any program or other work which
contains a notice placed by the copyright holder saying it may be
distributed under the terms of this General Public License. The Program''
below refers to any such program or work, and a work based on the
Program'' means either the Program or any derivative work under copyright
law: that is to say, a work containing the Program or a portion of
it, either verbatim or with modifications and/or translated into another
language. (Hereinafter, translation is included without limitation
in the term modification.'') Each licensee is addressed as you.''\\
\\
Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope. The act of running
the Program is not restricted, and the output from the Program is
covered only if its contents constitute a work based on the Program
(independent of having been made by running the Program). Whether
that is true depends on what the Program does.

\end{itemize}
\begin{enumerate}
\item You may copy and distribute verbatim copies of the Program's source
code as you receive it, in any medium, provided that you conspicuously
and appropriately publish on each copy an appropriate copyright notice
and disclaimer of warranty; keep intact all the notices that refer
to this License and to the absence of any warranty; and give any other
recipients of the Program a copy of this License along with the Program.

You may charge a fee for the physical act of transferring a copy,
and you may at your option offer warranty protection in exchange for
a fee.

\item You may modify your copy or copies of the Program or any portion of
it, thus forming a work based on the Program, and copy and distribute
such modifications or work under the terms of Section 1 above, provided
that you also meet all of these conditions:

\begin{enumerate}
\item You must cause the modified files to carry prominent notices stating
that you changed the files and the date of any change.
\item You must cause any work that you distribute or publish, that in whole
or in part contains or is derived from the Program or any part thereof,
to be licensed as a whole at no charge to all third parties under
\item If the modified program normally reads commands interactively when
run, you must cause it, when started running for such interactive
use in the most ordinary way, to print or display an announcement
including an appropriate copyright notice and a notice that there
is no warranty (or else, saying that you provide a warranty) and that
users may redistribute the program under these conditions, and telling
the user how to view a copy of this License. (Exception: if the Program
itself is interactive but does not normally print such an announcement,
your work based on the Program is not required to print an announcement.)
\end{enumerate}

These requirements apply to the modified work as a whole. If identifiable
sections of that work are not derived from the Program, and can be
reasonably considered independent and separate works in themselves,
then this License, and its terms, do not apply to those sections when
you distribute them as separate works. But when you distribute the
same sections as part of a whole which is a work based on the Program,
the distribution of the whole must be on the terms of this License,
whose permissions for other licensees extend to the entire whole,
and thus to each and every part regardless of who wrote it.

Thus, it is not the intent of this section to claim rights or contest
your rights to work written entirely by you; rather, the intent is
to exercise the right to control the distribution of derivative or
collective works based on the Program.

In addition, mere aggregation of another work not based on the Program
with the Program (or with a work based on the Program) on a volume
of a storage or distribution medium does not bring the other work
under the scope of this License.

\item You may copy and distribute the Program (or a work based on it, under
Section 2) in object code or executable form under the terms of Sections
1 and 2 above provided that you also do one of the following:

\begin{enumerate}
\item Accompany it with the complete corresponding machine-readable source
code, which must be distributed under the terms of Sections 1 and
2 above on a medium customarily used for software interchange; or,
\item Accompany it with a written offer, valid for at least three years,
to give any third party, for a charge no more than your cost of physically
performing source distribution, a complete machine-readable copy of
Sections 1 and 2 above on a medium customarily used for software interchange;
or,
\item Accompany it with the information you received as to the offer to
distribute corresponding source code. (This alternative is allowed
only for noncommercial distribution and only if you received the program
in object code or executable form with such an offer, in accord with
Subsection b above.)
\end{enumerate}

The source code for a work means the preferred form of the work for
making modifications to it. For an executable work, complete source
code means all the source code for all modules it contains, plus any
associated interface definition files, plus the scripts used to control
compilation and installation of the executable. However, as a special
exception, the source code distributed need not include anything that
is normally distributed (in either source or binary form) with the
major components (compiler, kernel, and so on) of the operating system
on which the executable runs, unless that component itself accompanies
the executable.

If distribution of executable or object code is made by offering access
copy the source code from the same place counts as distribution of
the source code, even though third parties are not compelled to copy
the source along with the object code.

\item You may not copy, modify, sublicense, or distribute the Program except
as expressly provided under this License. Any attempt otherwise to
copy, modify, sublicense or distribute the Program is void, and will
who have received copies, or rights, from you under this License will
not have their licenses terminated so long as such parties remain
in full compliance.
\item You are not required to accept this License, since you have not signed
it. However, nothing else grants you permission to modify or distribute
the Program or its derivative works. These actions are prohibited
by law if you do not accept this License. Therefore, by modifying
or distributing the Program (or any work based on the Program), you
indicate your acceptance of this License to do so, and all its terms
and conditions for copying, distributing or modifying the Program
or works based on it.
\item Each time you redistribute the Program (or any work based on the Program),
to copy, distribute or modify the Program subject to these terms and
conditions. You may not impose any further restrictions on the recipients'
exercise of the rights granted herein. You are not responsible for
enforcing compliance by third parties to this License.
\item If, as a consequence of a court judgment or allegation of patent infringement
or for any other reason (not limited to patent issues), conditions
are imposed on you (whether by court order, agreement or otherwise)
you from the conditions of this License. If you cannot distribute
and any other pertinent obligations, then as a consequence you may
not distribute the Program at all. For example, if a patent license
would not permit royalty-free redistribution of the Program by all
those who receive copies directly or indirectly through you, then
the only way you could satisfy both it and this License would be to
refrain entirely from distribution of the Program.

If any portion of this section is held invalid or unenforceable under
any particular circumstance, the balance of the section is intended
to apply and the section as a whole is intended to apply in other
circumstances.

It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any
such claims; this section has the sole purpose of protecting the integrity
of the free software distribution system, which is implemented by
to the wide range of software distributed through that system in reliance
on consistent application of that system; it is up to the author/donor
to decide if he or she is willing to distribute software through any
other system and a licensee cannot impose that choice.

This section is intended to make thoroughly clear what is believed
to be a consequence of the rest of this License.

\item If the distribution and/or use of the Program is restricted in certain
countries either by patents or by copyrighted interfaces, the original
an explicit geographical distribution limitation excluding those countries,
so that distribution is permitted only in or among countries not thus
excluded. In such case, this License incorporates the limitation as
if written in the body of this License.
\item The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time. Such new versions
will be similar in spirit to the present version, but may differ in
detail to address new problems or concerns.

Each version is given a distinguishing version number. If the Program
specifies a version number of this License which applies to it and
any later version,'' you have the option of following the terms
and conditions either of that version or of any later version published
by the Free Software Foundation. If the Program does not specify a
version number of this License, you may choose any version ever published
by the Free Software Foundation.

\item If you wish to incorporate parts of the Program into other free programs
whose distribution conditions are different, write to the author to
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this. Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software
and of promoting the sharing and reuse of software generally.
\end{enumerate}

\subsection*{NO WARRANTY }

\begin{itemize}

\item[11.]BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS
NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE
LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS
AND/OR OTHER PARTIES PROVIDE THE PROGRAM AS IS'' WITHOUT WARRANTY
OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED
TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
PROGRAM IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME
THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.

\item[12.]IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED
TO IN WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY
MODIFY AND/OR REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE
TO YOU FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED
INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE
OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER
OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

\end{itemize}

\section*{END OF TERMS AND CONDITIONS }

\subsection*{How to Apply These Terms to Your New Programs}

If you develop a new program, and you want it to be of the greatest
possible use to the public, the best way to achieve this is to make
it free software which everyone can redistribute and change under
these terms.

To do so, attach the following notices to the program. It is safest
to attach them to the start of each source file to most effectively
convey the exclusion of warranty; and each file should have at least
the copyright'' line and a pointer to where the full notice is found.
For example:
\begin{quote}
One line to give the program's name and a brief idea of what it does.

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published
by the Free Software Foundation; either version 2 of the License,
or (at your option) any later version.

This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
\end{quote}
Also add information on how to contact you by electronic and paper
mail.

If the program is interactive, make it output a short notice like
this when it starts in an interactive mode:
\begin{quote}
comes with ABSOLUTELY NO WARRANTY; for details type show w'. This
is free software, and you are welcome to redistribute it under certain
conditions; type show c' for details.
\end{quote}
The hypothetical commands show w' and show c' should show the appropriate
parts of the General Public License. Of course, the commands you use
may be called something other than show w' and show c'; they could

You should also get your employer (if you work as a programmer) or
your school, if any, to sign a copyright disclaimer'' for the program,
if necessary. Here is a sample; alter the names:
\begin{quote}
Yoyodyne, Inc., hereby disclaims all copyright interest in the program
Gnomovision' (which makes passes at compilers) written by James Hacker.

(signature of Ty Coon)\\
1 April 1989 \\
Ty Coon, President of Vice
\end{quote}
into proprietary programs. If your program is a subroutine library,
you may consider it more useful to permit linking proprietary applications
with the library. If this is what you want to do, use the GNU Library
\begin{thebibliography}{10}
\bibitem[1]{Moresi et al Plate tectonics}Moresi, L., M. Gurnis, and
S. Zhong (2000), Plate tectonics and convection in the Earth's mantle:
Toward a numerical simulation, \emph{Comp. Sci. Engin., 2}, 22-33.

\bibitem[2]{Moresi/Solomatov Numerical}Moresi, L.N., and V.S. Solomatov
(1995), Numerical investigation of 2D convection with extremely large
viscosity variations, \emph{Phys. Fluid,} \emph{7}, 2,154-2,162.

\bibitem[3]{Moresi/Zhong/Gurnis The accuracy}Moresi, L., S. Zhong,
and M.Gurnis (1996), The accuracy of finite element solutions of Stokes'
flow with strongly varying viscosity, \textit{Phys. Earth Planet.
Inter.}, 97, 83-94.

\bibitem[4]{Moresi/Gurnis Contraints}Moresi, L.N., and M. Gurnis
(1996), Constraints on the lateral strength of slabs from three-dimensional
dynamic flow models, \emph{Earth Planet. Sci. Lett.,} \emph{138},
15-28.

\bibitem[5]{Zhong/Moresi The role of faults}Zhong, S., M. Gurnis,
and L. Moresi (1998), The role of faults, nonlinear rheology, and
viscosity structure in generating plates from instantaneous mantle
flow models, \emph{J. Geophys. Res.,} \emph{103}, 15,255-15,268.

\bibitem[6]{Zhong et al The role of temperature}Zhong, S., M.T. Zuber,
L.N. Moresi, and M. Gurnis (2000), The role of temperature-dependent
viscosity and surface plates in spherical shell models of mantle convection,
\emph{J. Geophys. Res., 105}, 11,063-11,082.

\bibitem[7]{Tan et al Slabs in the lower}Tan, E., M. Gurnis, and
L. Han (2002), Slabs in the lower mantle and their modulation of plume
formation, \emph{Geochem. Geophys. Geosys.,} \emph{3}, 1067.

Gurnis (2003), Seismic tomography, surface uplift, and the breakup
of Gondwanaland: Integrating mantle convection backwards in time,
\emph{Geochem. Geophys. Geosys.,} \emph{4(3)}, 1031, doi:10.1029/2001GC000299.

\bibitem[9]{Hughes The Finite Element Method}Hughes, T.J.R. \emph{The
Finite Element Method: Linear Static and Dynamic Finite Element Analysis.}
Englewood Cliffs, New Jersey: Prentice-Hall, Inc.; 1987. 672 p.

\bibitem[10]{Ramage/Wathen Iterative solution}Ramage, A., and A.J.
Wathen (1994), Iterative solution techniques for the Stokes and Navier-Stokes
equations, \emph{Int. J. Numer. Methods. Fluids,} \emph{19}, 67-83.

\bibitem[11]{Brooks A.N.}Brooks, A.N. \emph{A Petrov-Galerkin Finite
Element Formulation for Convecton Dominated Flows.} Unpublished doctoral
thesis, California Institute of Technology, Pasadena, CA, 1981.

\bibitem[12]{Cahouet/Chabard Some fast 3D}Cahouet, J., and J.-P.
Chabard (1988), Some fast 3D finite element solvers for the generalized
Stokes problem\emph{, Int. J. Numer. Methods. Fluids,} \emph{8}, 869-895.

\bibitem[13]{Atanga/Silvester Iterative methods}Atanga, J., and D.
Silvester (1992), Iterative methods for stabilized mixed velocity-pressure
finite elements, \emph{Int. J. Numer. Methods. Fluids,} \emph{14},
71-81.

\bibitem[14]{Hager/OConnell A simple global}Hager, B.H., and R.J.
O'Connell (1981), A simple global model of plate dynamics and mantle
convection, \emph{J. Geophys. Res.,} \emph{86}, 4,843-4,867.

\bibitem[15]{Tan et al GeoFramework Part I}Tan, E., E. Choi, P. Thoutireddy,
M. Gurnis, and M. Aivazis (2006), GeoFramework: Coupling multiple
models of mantle convection within a computational framework,  \emph{Geochem.,
Geophys., Geosyst. 7,} Q06001, doi:10.1029/2005GC001155.

\bibitem[16]{Zhong et al Free-surface formulation}Zhong, S., M. Gurnis,
and L. Moresi (1996), Free-surface formulation of mantle convection--I.
Basic theory and appication to plumes\emph{, Geophys. J. Int.}, \emph{127},
708-718.

\bibitem[17]{Zhong et al Three-dimensional}Zhong, S., A. Paulson,
and J. Wahr (2003), Three-dimensional finite-element modeling of Earth's
viscoelastic deformation: effects of lateral variations in lithospheric
thickness, \emph{Geophys. J. Int.,} \emph{155}, 679-695.

\bibitem[18]{McNamara Zhong Thermochemical}McNamara, A.K., and S.
Zhong (2004), Thermochemical structures within a spherical mantle:
Superplumes or Piles? \emph{J. Geophys. Res.}, \textit{109}, B07402,
doi:10.1029/2003JB002847.

\bibitem[19]{Tackely King tracer ratio method}Tackley, P.J., and
S.D. King (2003), Testing the tracer ratio method for modeling active
compositional fields in mantle convection simulations, \emph{Geochem.
Geophys. Geosyst., 4,} 8302, doi:10.1029/2001GC000214.

\bibitem[20]{Trampert mineral physics model 2001}Trampert, J., P.
Vacher, and N. Vlaar (2001), Sensitivities of seismic velocities to
temperature, pressure and composition in the lower mantle, \emph{Phys.
Earth Planet. Inter.}, \textit{124}, 255-267.\end{thebibliography}

\end{document}
`