\documentclass[english]{book}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{geometry}
\geometry{verbose,tmargin=1in,bmargin=1in,lmargin=1in,rmargin=1in}
\pagestyle{headings}
\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.3.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}


\section*{About This Document}

This document is organized into three parts. Part I consists of traditional
book front matter, including this preface. Part II describes the capabilities of 
CitcomS and 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}
Additionally, if you are using tracers in CitcomS, please cite the
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
your paper acknowledgements that you indicate your use of this code,
the authors of the code, and CIG \url{geodynamics.org}.


\section*{Support}

Pyre development was 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
Public License (see Appendix \vref{cha:License}). Written in C, the
code runs on a variety of parallel processing computers, including
shared and distributed memory platforms. 


\section{About CitcomS}
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
\url{svn.enthought.com/enthought/wiki/MayaVi} in addition to Generic
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
mesh refinement for the radial coordinate were added by Thorsten Becker. 

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
contributed added capability to read in time- and geographic-dependent,
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}
The release of CitcomS 3.2 added 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. Other enhancements include:
\begin{itemize}
\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.
\item Added test implementation of steady state strain-rate weakening via
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 Replaced fixed accuracy setting. Bug reported by Rob Moucha.
\item Fixed temperature perturbation bug for multiple processors. 
\end{itemize}
\end{itemize}

This release (3.3) removes the Python dependencies from CitcomS and provides
a tool to convert Python style parameter files to the original CitcomS
parameter files.  Python support was proving problematic because of
incompatibility with versions of Python newer than 2.6 (final release 2008)
and other internal issues.  The removal of Python should
simplify much of the code while retaining almost all functionality.
A full list of the changes includes:

\begin{itemize}
\item Removed Python sections of CitcomS
\item Wrote Python to standard parameter file conversion tool (Py2C)
\item Converted Python cookbooks parameter files to original style parameter files
\item Added output of parameters to a pidXXXXXX style file
\item Updated the manual to correspond to the new changes
\end{itemize}

\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}
\begin{equation}
\left(\rho u_{i}\right)_{,i}=0\label{eq:conservation of mass}
\end{equation}

\par\end{flushright}

\begin{flushright}
\begin{equation}
-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}
\end{equation}
 
\par\end{flushright}

\begin{flushright}
\begin{equation}
\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}
\end{equation}

\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}
\begin{equation}
\delta\rho=-\alpha\bar{\rho}(T-\bar{T_{a}})+\delta\rho_{ph}\Gamma+\delta\rho_{ch}C\label{eq:density anomalies}
\end{equation}

\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}
\begin{equation}
\pi=\bar{\rho}g(1-r-d_{ph})-\gamma_{ph}(T-T_{ph})\label{eq:reduced pressure}
\end{equation}

\par\end{flushright}

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

\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}
\begin{equation}
\rho=\rho_{0}\rho^{'}\label{eq:rho dim}
\end{equation}

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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


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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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

\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:

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


\begin{equation}
-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}
\end{equation}


\begin{flushright}
\begin{equation}
\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}
\end{equation}

\par\end{flushright}

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

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

\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}
\begin{equation}
Ra_{b}=Ra\frac{\delta\rho_{ph}}{\rho_{0}\alpha_{0}\Delta T}\label{eq:Rab, the phase-change Rayleigh number}
\end{equation}

\par\end{flushright}

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

\par\end{flushright}

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

\par\end{flushright}

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



\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}
\begin{equation}
\mathrm{\left(\mathbf{B}^{\mathit{T}}+\mathbf{C}\right)\mathit{u=0}}\label{eq:discrete continuite eqn}
\end{equation}

\par\end{flushright}

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

\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}
\begin{equation}
\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}
\end{equation}

\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.

\begin{equation}
\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}
\end{equation}


\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:
\begin{equation}
\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}
\end{equation}
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,
radius\_inner), and advances from 1 at the base to \emph{nodez} at
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})
as well as in radius. 

\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.
Advanced users and software developers may be interested in downloading
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
\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.

\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}

\section{Downloading and Unpacking Source}

To obtain CitcomS, go to the Geodynamics Software Packages web page
\url{geodynamics.org/cig/software/citcoms}, download
the source archive and unpack it using the \texttt{tar} command: 
\begin{lyxcode}
\$~tar~xzf~CitcomS-3.3.0.tar.gz
\end{lyxcode}

If you don't have GNU Tar, try the following command instead: 
\begin{lyxcode}
\$~gunzip~-c~CitcomS-3.3.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.3.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.

\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 
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: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}

\subsection{Additional Tools}

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
arrays for use in scientific computing. You may download NumPy from
the NumPy home page \url{numpy.scipy.org}. To compile and install
this extension, download it and issue the following commands after
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
HDF5 files. You may download it from the HDFView home page \url{hdf.ncsa.uiuc.edu/hdf-java-html/hdfview}.


\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
\url{www.cactuscode.org/Visualization/openDX} and follow the instructions
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
and header files). See the \texttt{README.TXT} file in HC for instructions
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{\label{sec:Software-Repository}Installing from the Software Repository}

The CitcomS source code is available via Git 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 Git client installed. To check, type
\texttt{git}; it should return a usage message.
\begin{lyxcode}
\$~git

usage: git [--version] [--help] [-C <path>] [-c name=value]
           [--exec-path[=<path>]] [--html-path] [--man-path] [--info-path]
           [-p|--paginate|--no-pager] [--no-replace-objects] [--bare]
           [--git-dir=<path>] [--work-tree=<path>] [--namespace=<name>]
           <command> [<args>]
\end{lyxcode}
For more information on Git, we recommend one of the following websites:

http://github.com/
http://git-scm.com/
https://www.atlassian.com/git/

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.3.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{git 
clone} command:
\begin{lyxcode}
\$~git~clone~https://github.com/geodynamics/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{git pull update} command:
\begin{lyxcode}
\$~cd~citcoms

\$~git~pull
\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}

Upon the completion of a successful build, 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}.

\section{Changing Parameters }

A configuration file should be used to set the input parameters for CitcomS.
A configuration file is a plain text file whose format is based on the
Windows INI file. The format of a configuration file is as follows:

\begin{lyxcode}
\#CitcomS.component1.component2\\
\#~this~is~a~comment\\
property1=value1\\
property2=value2~~;~this~is~another~comment
property3=val1,val2,val3~~;~this~property~is~specified~by~a~list~of~values
\end{lyxcode}

Upon termination of each run, all of the parameters are logged in a 
\texttt{pidXXXXXX} file, where \texttt{XXXXXX} is the process id of the
CitcomS application. Section names such as 
\texttt{\#CitcomS.component1.component2} are optional, but are recommended for 
the purpose of organizing the contents of the configuration file

\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. 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 the \texttt{example0} configuration file. Switch to the \texttt{examples}
directory and execute the following on the command line:
\begin{lyxcode}
\$~CitcomSRegional example0
\end{lyxcode}
or, equivalently,
\begin{lyxcode}
\$~mpirun~-np~1~CitcomSRegional example0
\end{lyxcode}

\subsubsection{Example: Uniprocessor, \texttt{example0}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=5

# CitcomS.controller
storage_spacing=1

# CitcomS.solver
datafile=example0
rayleigh=100000

# CitcomS.solver.mesher
fi_max=1
fi_min=0
nodex=17
nodey=17
nodez=9
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=5
perturbm=1
perturbmag=0.05

# CitcomS.solver.visc
num_mat=4
\end{verbatim}

\section{\label{sec:Multiprocessor-Example}Multiprocessor Example}

In the \texttt{examples} directory, type the following at the command line:
\begin{lyxcode}
\$~mpirun~-np~4~CitcomSRegional example1
\end{lyxcode}

\subsubsection{Example: Multiprocessor, \texttt{example1}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=70

# CitcomS.controller
storage_spacing=10

# CitcomS.solver
datafile=example1
rayleigh=100000

# CitcomS.solver.mesher
fi_max=1
fi_min=0
nodex=17
nodey=17
nodez=9
nprocx=2
nprocy=2
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=5
perturbm=1
perturbmag=0.05

# CitcomS.solver.visc
num_mat=4
\end{verbatim}

This example has:
\begin{itemize}
  \item 2 processors in colatitude ($x$-coordinate), specified by 
    the \texttt{nprocx=2} line in the configuration file
  \item 2 processors in longitude ($y$-direction), specified by 
    the \texttt{nprocy=2} line
  \item 1 processor in the radial ($z$-direction). The default value for
    the \texttt{nprocx,nprocy,nprocz} parameters is 1. Hence the 
    \texttt{nprocz=1} line is omitted from the configuration file.
  \item The total number of processors used is given by
  $\mathtt{nprocx}\times\mathtt{nprocy}\times\mathtt{nprocz}$. Hence this
    model uses 4 processors in total.
  \item There are 17 nodes in $x$ (theta) direction.
  \item There are 17 nodes in $y$ (fi) direction.
  \item There are 9 nodes in $z$ (radius\_inner) direction.
\end{itemize}
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}, 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 cluster or supercomputer or the number of CPUs in a given 
computer. We run one ``process'' on each node. In the example script 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 Formats}
The possible output formats in CitcomS are \texttt{ascii,ascii-gz,hdf5,vtk} with
\texttt{ascii} being the default. The output format is specified by the
\texttt{output\_format} configuration parameter. 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}
  output\_format=ascii-gz
\end{lyxcode}
Be warned that the post-process scripts do not understand this output
format yet.

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}
  output\_format=hdf5
\end{lyxcode}
The output files will be stored in 
\texttt{./}
with a filename prefix \texttt{example1} and a filename suffix \texttt{h5}.
See Chapter \vref{cha:Working-with-HDF5} for more information on
how to work with the HDF5 output.


\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.


\section{About HDF5}

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
are accompanied by supporting metadata.

A dataset is physically stored in two parts: a header and a data array.
The header contains miscellaneous metadata describing the dataset
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.

HDF5 also supports access to portions (or selections) of a dataset
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}.
You may adjust your visualization pipeline by adding any appropriate
VTK filters and modules. For more details about this process, refer
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. Cookbook examples range
from regional to full spherical shell problems that address traditional
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
more information).

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.


\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
on your cluster. 

\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}. 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~\\
storage\_spacing=25\\
datafile=cookbook1
\end{lyxcode}
The following parameters control 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}
\$mpirun -np 12 CitcomSFull cookbook1
\end{lyxcode}

\subsubsection{Example: Global Model, \texttt{cookbook1}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=100
solver=full

# CitcomS.controller
storage_spacing=25

# CitcomS.solver
datafile=cookbook1
rayleigh=100000

# CitcomS.solver.mesher
fi_max=1
fi_min=0
nodex=9
nodey=9
nodez=9
nproc_surf=12
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.ic
num_perturbations=1
perturbl=3
perturblayer=5
perturbm=2
perturbmag=0.05

# CitcomS.solver.visc
num_mat=4
\end{verbatim}
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.

\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}. 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\\
storage\_spacing=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
radial (\texttt{theta}, \texttt{fi}, \texttt{radius}) directions.
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\\
radius\_inner=0.55\\
radius\_outer=1.0
\end{lyxcode}
The following parameters allow 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 following parameters. 
\begin{lyxcode}
num\_perturbations=1\\
perturbmag=0.0
\end{lyxcode}

\subsubsection{Example: Velocity Boundary Conditions, \texttt{cookbook2}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=60

# CitcomS.controller
storage_spacing=30      ; how often outputs are created

# CitcomS.solver
datafile=cookbook2      ; prefix of output filenames
rayleigh=100000

# Modify the layout of the mesh
# CitcomS.solver.mesher
nodex=17
nodey=17
nodez=9
nprocx=2
nprocy=2

theta_max=1.5708
theta_min=0.7854
fi_max=0.7854
fi_min=0.0
radius_inner=0.55
radius_outer=1.0

# Import a uniform velocity across the top surface.
# CitcomS.solver.bc
topvbc=1
topvbxval=100
topvbyval=0

# In addition, set the initial temperature perturbation to zero.
# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=5
perturbm=1
perturbmag=0.0

# CitcomS.solver.visc
num_mat=4
\end{verbatim}

\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}. 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 parameters specify whether the viscosity should be updated at every time
step (\texttt{on}), the number of viscous layers (4), the viscosity of
>>>>>>> python-removal
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
=======
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
>>>>>>> python-removal
\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: 

\begin{equation}
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}
\end{equation}
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
\begin{equation}
viscE=\frac{E_{a}}{R\Delta T}\label{eq:viscE}
\end{equation}
\begin{equation}
viscZ=\frac{\rho gV_{a}}{R\Delta T}\label{eq:viscZ}
\end{equation}
\begin{equation}
viscT=T_{0}\label{eq:viscT}
\end{equation}

\subsubsection{Example: Temperature-Dependent Viscosity, \texttt{cookbook3}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=200

# CitcomS.controller
storage_spacing=25

# CitcomS.solver
datafile=cookbook3
rayleigh=1e6

# CitcomS.solver.mesher
fi_max=1
fi_min=0
nodex=17
nodey=17
nodez=9
nprocx=2
nprocy=2
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=5
perturbm=1
perturbmag=0.05

# CitcomS.solver.visc
TDEPV=on
VISC_UPDATE=on
VMAX=on
VMIN=on
num_mat=4
rheol=4
visc0=1,1,1,1
viscE=0.2,0.2,0.2,0.2
viscT=0,0,0,0
viscZ=0,0,0,0
visc_max=100.0
visc_min=1.0
\end{verbatim}
\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{coor}. Set it to \texttt{1} \texttt{(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}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=250

# CitcomS.controller
storage_spacing=50

# CitcomS.solver
datafile=cookbook4
rayleigh=1e6

# CitcomS.solver.mesher
coor=1
coor_file=coor.dat
fi_max=1
fi_min=0
nodex=33
nodey=17
nodez=17
nprocx=1
nprocy=1
nprocz=1
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=10
perturbm=0
perturbmag=0.05

# CitcomS.solver.visc
TDEPV=on
VISC_UPDATE=on
VMAX=on
VMIN=on
num_mat=4
rheol=4
visc0=1,1,1,1
viscE=0.2,0.2,0.2,0.2
viscT=0,0,0,0
viscZ=0,0,0,0
visc_max=100.0
visc_min=1.0
\end{verbatim}
\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\\
datadir\_old=./ic\\
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}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=1100

# CitcomS.controller
storage_spacing=100

# CitcomS.solver
datadir_old=./ic
datafile=cookbook5
datafile_old=cookbook5
rayleigh=4.07e+08

# CitcomS.solver.mesher
coor=1
coor_file=./coor.dat
fi_max=1
fi_min=0
nodex=17
nodey=65
nodez=33
nprocx=1
nprocy=2
nprocz=1
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.tsolver
finetunedt=0.75
monitor_max_T=on

# CitcomS.solver.bc
topvbc=1

# CitcomS.solver.ic
num_perturbations=1
perturbl=1
perturblayer=5
perturbm=1
perturbmag=0.05
solution_cycles_init=0
tic_method=-1

# CitcomS.solver.param
file_vbcs=on
start_age=55
vel_bound_file=./velocity/bvel.dat

# CitcomS.solver.visc
TDEPV=on
VMAX=on
VMIN=on
num_mat=4
visc0=100,0.003,1,2
viscE=24,24,24,24
viscT=0.182,0.182,0.182,0.182
visc_max=100.0
visc_min=0.01
\end{verbatim}
\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.

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}. 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}
datadir=output
\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}
Solver=cgrad\\
accuracy=1e-03\\
vlowstep=1000\\
piterations=1000
\end{lyxcode}

\subsubsection{Example: Thermo-Chemical Convection, \texttt{cookbook7}}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=15
solver=full

# CitcomS.controller
storage_spacing=5

# CitcomS.solver
datadir=output
datafile=cookbook7
rayleigh=1e7

# CitcomS.solver.mesher
fi_max=1
fi_min=0
nodex=9
nodey=9
nodez=9
nproc_surf=12
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.vsolver
Solver=cgrad
accuracy=1e-04
piterations=1000
vlowstep=1000

# CitcomS.solver.ic
num_perturbations=1
perturbl=3
perturblayer=5
perturbm=2
perturbmag=0.05

# CitcomS.solver.output
output_optional=tracer,comp_nd

# CitcomS.solver.tracer

# CitcomS.solver.visc
TDEPV=on
VISC_UPDATE=on
VMAX=on
VMIN=on
num_mat=4
rheol=4
visc0=1,1,1,1
viscE=0.2,0.2,0.2,0.2
viscT=0,0,0,0
viscZ=0,0,0,0
visc_max=100.0
visc_min=1.0
\end{verbatim}
\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}
~~~~nz~~~~~radius~~~~~~depth~~~~rho~~~~~~~~~~~~~~layer~\\
~~~~~~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:

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


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{}


\section{Cookbook 8: Compressible Steady-State Convection}


\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 
is provided in the \texttt{examples/Cookbook8/} directory. To restart the
computation from this checkpoint, set these parameters:
\begin{lyxcode}
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: 
\begin{equation}
\mathrm{nodex}=1+\mathrm{nprocx}\times\mathrm{mgunitx}\times2^{levels-1}\label{eq:mgunit}
\end{equation}
where \texttt{mgunitx}, which \textbf{must be an integer}, 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
gradient solver (\texttt{Solver=cgrad}) for larger problems (e.g.,
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}.


\subsubsection{Example: Compressible Steady-State Convection, cookbook8}
\begin{verbatim}
# CitcomS
cpu_limits_in_seconds=360000000
minstep=10000
solver=full

# CitcomS.controller
checkpointFrequency=1000
storage_spacing=1000

# CitcomS.solver
datadir=output
datadir_old=restart
datafile=cookbook8
datafile_old=cookbook8
dissipation_number=0.5
gruneisen=1.25
rayleigh=7.68175583e4
surfaceT=0.1

# CitcomS.solver.mesher
coor=1
coor_file=coord.dat
fi_max=1
fi_min=0
levels=5
mgunitx=2
mgunity=2
mgunitz=2
nodex=33
nodey=33
nodez=33
nproc_surf=12
theta_max=2.0708
theta_min=1.0708

# CitcomS.solver.tsolver
finetunedt=0.75

# CitcomS.solver.vsolver
Solver=multigrid
accuracy=0.001
aug_lagr=on
aug_number=2.0e3
compress_iter_maxstep=100
down_heavy=2
max_mg_cycles=50
mg_cycle=1
piterations=375
precond=on
remove_angular_momentum=on
remove_rigid_rotation=off
up_heavy=2
uzawa=cg
vhighstep=2
vlowstep=20

# CitcomS.solver.ic
num_perturbations=1
perturbl=3
perturblayer=17
perturbm=2
perturbmag=0.01
restart=off
solution_cycles_init=9000
tic_method=3

# CitcomS.solver.output
output_ll_max=20
output_optional=geoid,surf,botm
self_gravitation=on
use_cbf_topo=on

# CitcomS.solver.param
reference_state=1
refstate_file=ref.dat

# CitcomS.solver.visc
TDEPV=on
VISC_UPDATE=on
VMAX=on
VMIN=on
num_mat=4
rheol=1
visc0=1,1,1,1
viscE=2.99573,2.99573,2.99573,2.99573
viscT=0.5,0.5,0.5,0.5
visc_max=1e+06
visc_min=0.001
visc_smooth_method=1
\end{verbatim}

\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}
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}

\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}

\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:

\begin{equation}
\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}
\end{equation}



\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.


\chapter{\label{cha:License}License }

\textbf{GNU GENERAL PUBLIC LICENSE Version 2, June 1991. Copyright
(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
Public License instead.) You can apply it to your programs, too.

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 Copyright the software, and 
\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
the terms of this License. 
\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
the corresponding source code, to be distributed under the terms 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
to copy from a designated place, then offering equivalent access to
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
automatically terminate your rights under this License. However, parties
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),
the recipient automatically receives a license from the original licensor
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)
that contradict the conditions of this License, they do not excuse
you from the conditions of this License. If you cannot distribute
so as to satisfy simultaneously your obligations under this License
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
public license practices. Many people have made generous contributions
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
copyright holder who places the Program under this License may add
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
ask for permission. For software which is copyrighted by the Free
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.
Copyright {\footnotesize{© (}}year) (name of author) 

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}
Gnomovision version 69, Copyright © year name of author Gnomovision
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
even be mouse-clicks or menu items -- whatever suits your program. 

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}
This General Public License does not permit incorporating your program
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
General Public License instead of this License.
\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.

\bibitem[8]{Conrad/Gurnis Seismic tomography}Conrad, C.P., and M.
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}