configuring_fluidity.tex

\chapter{Configuring \fluidity}\label{chap:configuration}

\section{Overview}
A \fluidity\ simulation is configured by creating a \fluidity\ options (or flml) file using
Diamond, which is the Graphical User Interface (GUI). 
The left-hand pane of Diamond allows users to browse the
\emph{options tree}, while the right-hand pane provides brief documentation
about the option and is where users enter input data. 

This chapter aims to provide a detailed description of all the options in
the tree.  From section~\ref{sec:OptionsTree} onwards, \fluidity\ options are
described in the order in which they appear in Diamond.  Prior to this are
some important general notes about the different types of options and, in
particular, how to work with fields in \fluidity.

\section{Options syntax}
\index{options!syntax}
\fluidity\ options files, or flml files, are XML files whose grammar is defined by the
fluidity options schema \verb+fluidity_options.rng+. XML files have a
tree-like structure of elements containing other elements. This structure is
reflected in the left hand pane of Diamond, the GUI which
is used to write flml files.

The flml file can also be edited with a standard editor and the options written
in text and in code using the Spud library on which the \fluidity\ options
system is based. A location in the options tree can be written as a path,
much like the path of a file on disk. So, for example, the option
controlling the physical dimension of the problem domain is written as
\option{/geometry/dimension}. This should be read as the \option{dimension}
option which is found under \option{geometry} which is in turn at the top
level of the options tree. In Diamond, and in the flml file, the absolute
top level element is always \option{fluidity\_options} but this element is
always discarded from paths. Figure \ref{fig:geometry_dimension}\ shows a
Diamond screen shot open at the \option{/geometry/dimension} option. Further
documentation of the Spud system is available in \citet{ham2009}\ and from
\href{http://amcg.ese.ic.ac.uk/spud}{the Spud website}.

\begin{figure}[ht]
  \centering
  \fig[width=.7\textwidth]{geometry_dimension}
  \caption{A Diamond screenshot showing the \option{/geometry/dimension}
    option. Note that the option path is displayed at the bottom of the
    diamond window}
  \label{fig:geometry_dimension}
\end{figure}

\subsection{Allowed Characters}
Only certain characters are recognised by the options dictionary, which contains the flml input once it is read into fluidity. Therefore only the following letters are allowed in any input field:

\begin{verbatim}
 /_:[]1234567890qwertyuioplkjhgfdsazxcvbnmMNBVCXZASDFGHJKLPOIUYTREWQ
\end{verbatim}

Comment boxes may contain any characters. 

\subsection{Named options}
\index{options!names}
Some options in the tree, such as fields and meshes, have name
attributes. The name attribute is represented in the flml file with a double colon so that,
for example, the coordinate mesh has options path\onlypdf\linebreak
\option{/geometry/mesh::CoordinateMesh}. Note that this differs from the
convention in the Diamond interface in which name attributes are given in brackets. Figure
\ref{fig:mesh_name}

\begin{figure}[ht]
  \centering
  \fig[width=.7\textwidth]{mesh_name}
  \caption{A Diamond screenshot showing the \option{/geometry/mesh::Coordinate}
    option. Note that the name is shown in brackets in the main Diamond
    window but after double colons in the path in the bottom bar.}
  \label{fig:mesh_name}
\end{figure}

Names of objects (fields, material phases, meshes, etc.) should be camel cased (MyOwnField) and not contain spaces or underscores. Furthermore, the characters \verb+/:[]+ are prohibited as these have special meanings in the options dictionary inside \fluidity. 

\section{The options tree}\label{sec:OptionsTree}
The top level of the options tree contains the following compulsory elements:

\begin{itemize}
\item Simulation Name
\item Problem Type
\item Geometry
\item IO
\item Timestepping
\item Physical Parameters
\item Material/Phase
\end{itemize}

The first six of these are described here.

\subsection{Simulation Name}
The simulation name is the base name for all output files. For example if you set the simulation name to foo then your statistics output file will be called foo.stat

\subsection{Problem Type}
Setting problem type gives fluidity a hint as to what sort of simulation you are conducting and therefore what combinations of options are likely to be valid. If you do not know which category applies to your problem, choose "fluids". 
\subsection{Geometry}
This element contains all the options required to specify the geometry of the mesh and the accuracy of the finite element discretisation.

\subsubsection{Dimension}
\index{dimension}
The dimension of the domain of your problem. This can be 1, 2 or 3.
Be careful, once you set the dimension you can't change it again!
This is necessary to ensure that all vectors and tensors are of the correct
dimension.

\subsubsection{Meshes} \label{sec:mesh_configuration}
Meshes are the finite element spaces on which your problem is solved. Meshes
are either read from file or are derived from other meshes. Mesh options are
described in detail in section~\ref{sec:Mesh}.  There is only one required
mesh: the CoordinateMesh. Some settings or fields have
specific mesh requirements. These are discussed under the appropriate
options.

\subsubsection{Quadrature}
\index{quadrature!options}
\fluidity\ uses numerical quadrature to integrate the equations over each
element. There is a performance/accuracy trade off in quadrature: the more
quadrature points are employed, the more accurate the integrals are but
the more expensive the assembly operation is. Quadrature rules in \fluidity\ are categorised by the degree of the polynomial which they will integrate
exactly. The higher the degree of the quadrature rule, the more quadrature
points will be employed. As a general rule of thumb, the quadrature
degree employed should be at least $\max(2n_{\vec{u}+1},2n_p)$ where
$n_{\vec{u}}$ is the degree of the elements employed for velocity and $n_p$
is the degree of the elements employed for pressure. This means that degree
4 quadrature is sufficient for most of the fluidity configurations currently
in use.

The quadrature degree is specified by setting \option{/geometry/quadrature/degree}.

\subsubsection{Spherical Earth}
\label{sec:spherical_earth}
\index{spherical earth}
\index{mesh!superparametric}

Enabling \option{/geometry/spherical\_earth} informs \fluidity that your
simulation is being carried out in an Earth like geometry. This means that the mesh domain is assumed to approximate (a part of) an
N-sphere (sphere in three dimensions, circle in two dimensions). The domain is still described in
general $x,y,z$ coordinates (no spherical/cylindrical coordinates)
and the centre of the N-sphere must coincide with the origin.  Turning on this option has a number of effects.

In combination with extrusion options (see section \ref{sec:extruded_meshes}), 
specifying a \option{spherical\_earth} will extrude radially from the top surface of the
domain towards the origin.  Generally, the sphere is only extruded to a certain
depth, leaving a void around the origin. The two main applications of this
option are large-scale ocean simulations on the sphere, and geodynamic simulations of the
Earth's mantle.

When using this option the gravity direction should be defined manually in the radial
direction (pointing to the origin) using python
(\option{/physical\_parameters/gravity/vector\_field::GravityDirection/prescribed/value::WholeMesh/python}).  If used in conjunction with the
\option{buoyancy/radial\_gravity\_direction\_at\_gauss\_points} option underneath the \option{spatial\_discretisation} option for
Velocity then gravity will be automatically defined to be radial at the gauss points (rather than the nodes) in the buoyancy term
(ignoring whatever gravity direction you have specified manually in the buoyancy term only).  Specifying the gravity direction at the
gauss points in the buoyancy term helps eliminate spurious interpolation errors and non-physical imbalance in spherical simulations (see
\ref{sec:configuring_fluidity!spatial_discretisation!buoyancy}).

For ocean simulations, the option has implications for various options and terms such as wind forcing (see section \ref{sec:wind_forcing}), the calculation of buoyancy and the `direction' of absorptions (see \ref{sec:Source}).
If this option is checked, wind forcing and \eg momentum forcing from bulk formulae (see \ref{sec:bulk_formulae}) will automatically be rotated and applied in the direction tangential to the Earth's surface. It will also result in absorption terms set through the options tree being rotated and applied in the longitudinal, latitudinal and radial directions respectively. Note however that the viscosity and diffusion operators are currently not rotated automatically and thus the user must carry such rotations out themselves. An example in which the viscosity is rotated is giving in section \ref{sec:tides_in_the_med}.

Two options, \option{spherical\_earth/superparametric\_mapping} and
\option{spherical\_earth/analytical\_mapping},
are available to improve the approximation of the spherical domain
by the mesh, in particular for coarse meshes.
The general idea is that instead of using straight, linearly
interpolated tets and triangles, the elements are bent to better fit the curved
geometry (see section \ref{sec:coordinate_mesh_and_superparametric} for the
exact formula). The mesh vertices are kept in place, but all points in the interior of
an element are moved a little in the radial direction,
such that the radial coordinate varies linearly
between the vertices. For example, in the case of a triangle with all three
points on a perfect sphere, this means the entire triangle is now curved
to be exactly on that sphere. 

Using the \option{superparametric\_mapping} option
the curved elements are approximated by a higher order polynomial by using a
higher order \option{CoordinateMesh}. To use this option, you therefore need a
seperate \option{CoordinateMesh} that is derived from the usual continuous P1
mesh (typically called \option{BaseMesh}). For example, when combining this with
extrusion in 3D, one needs a 2D input mesh \option{InputMesh}, an extruded
\option{BaseMesh} with extrusion options, and a \option{CoordinateMesh} derived
from this that has set a higher order \option{mesh\_shape/polynomial\_degree}.
Other meshes (function spaces) used in the discretisation, e.g. a \PoDG or a $P2$ mesh, are
also derived from the \option{BaseMesh}.

For the \option{analytical\_mapping} option a higher order
\option{CoordinateMesh} is not necessary; the bending of the element is done
directly at the gauss points using an analytical expression. Thus the
conformance of a bended element with a perfect sphere is only limited by the
quadrature degree. Nonetheless, for calculations that are done node-wise, in
particular when setting a higher order prescribed field from a python
function, it may still make sense to combine \option{superparametric\_mapping}
and \option{analytical\_mapping}, so that the coordinates of the bended
elements are available in the nodes of the higher order mesh.

\subsubsection{Ocean Boundaries}\label{sec:ocean_boundaries}
\index{boundary conditions!ocean}
These options are required if you are running an ocean simulation with a
free surface or with various other ocean options which require the code to
know where the ocean surface and bed lie.
\option{/geometry/ocean\_boundaries/top\_surface\_ids} and
\onlypdf\linebreak\option{/geometry/ocean\_boundaries/bottom\_surface\_ids}
are lists of boundary tags from your input mesh which lie on the ocean
surface and bed respectively.

It is not usually necessary to change the settings for either of the scalar fields under this option. 

\subsection{IO}
These options control the frequency and form of model outputs.
\subsubsection{Dump format}
\index{vtk}
The file format used to output fields to disk. At this stage, vtu is the
only allowed format.

\subsubsection{Dump period}
This is the interval between the state fields being output to disk. You should usually start by setting this to a rather low value (possibly as short as your timestep) for testing and then increase it for production runs once you know how your configuration works. The value can be specified as either a constant or python function.

It is possible to swap \option{/io/dump\_period} for \option{/io/dump\_period\_in\_timesteps} to specify that you wish to have a dump every fixed number of timesteps.

\subsubsection{Output mesh}
\index{mesh!output}

All fields will be interpolated onto the same mesh for output. Usually the
CoordinateMesh is the right choice. If you have fields that are of
higher order than the selected output mesh you will lose 
accuracy. Interpolating all fields to a higher order mesh for output may give
very large dump files however. If any 
of the fields in the output is discontinuous, the mesh in the output file
will be a discontinuous version of the mesh selected here.

\subsubsection{Disable dump at start}

A dump is normally performed at the start of the simulation. This options disables that.

\subsubsection{Disable dump at end}

A dump is normally performed at the end of the simulation. This options disables that.

\subsubsection{CPU dump period}

This outputs dumps at specified CPU times. Not recommended.

\subsubsection{Wall time dump period}

Outputs at specified walltime (real time) periods. Not recommended.

\subsubsection{Max dump file count}

Limits the number of dumps by overwriting the previous dumps after the number specified here.

\subsubsection{Convergence}

You can check certain fields for convergence during nonlinear iterations.  To do this switch on \option{/timestepping/nonlinear\_iterations} and \option{/timestepping/nonlinear\_iterations/tolerance} and switch on the convergence option under the fields that you want to check for convergence. 

It is possible to enable the creation of a convergence file, giving details of the convergence of each field over the global nonlinear iteration loop. The .convergence file is in the same format as the .stat file. In order to do this, switch on 
\option{/io/convergence}  and \option{/io/convergence/convergence\_file}.  You still need the options in the above paragraph.

\subsubsection{Checkpointing}
\index{checkpointing}
\label{sec:configuring_fluidity_checkpointing}
Enables checkpointing, which saves sufficient information (including a new flml options file) to restart a simulation, i.e., to continue the simulation after it stopped. You must specify how often to checkpoint in terms of the number of dumps. There
are also options to checkpoint at the start of the simulation and at the
end. This latter is useful when running on batch systems that have a time
limit.

Up to five sets of files are created when checkpointing:
\begin{enumerate}
\item Mesh files - The from\_file meshes, in Gmsh format. Surface IDs
  and (if present) region IDs are written to the mesh files, and adaptivity
  is supported. In parallel a Gmsh mesh is written for each process.
\item Halo files (in parallel) - Halo information for each process.
\item Field files - A vtu is written for each mesh with prognostic fields in
  each state and (in parallel) for each process.
\item Checkpointed option file - A new FLML file, with the from\_file mesh
  set to read the checkpoint mesh files and the prognostic fields set
  to initialise from the checkpoint field files.
\item Checkpointed detector files - Two files related to checkpointing of
  detectors are created.
\end{enumerate}

The first checkpointed detector file has the extension .groups and contains
a header with the names of the groups of detectors in the order they were
read in the simulation. It also contains information about the number of
detectors in each group. This is to guarantee that when restarting from a
checkpoint, the detectors will be read in the same order and consequently
will be saved in that same order into the output detector file for
consistency.  The second checkpointed detector file has the extension
.positions.dat and contains the last position of the detectors at the time 
of checkpointing, in binary format.

No checkpointed detector files will be created if only static detectors are
defined in Diamond, since the position of these detectors remains always the
same. This is mainly used when Lagrangian detectors are set that are
advected by the flow and hence, their position changes as the simulation
proceeds.

At the same time as creating the two files related to checkpointing of
detectors, the detector options in the options tree or Diamond are updated
so that in the new flml file the detectors are set to initialise from the
checkpoint detectors files \\
(\option{/io/detectors/detector\_array/from\_checkpoint\_file} or \\
\option{/io/detectors/lagrangian\_detector/from\_checkpoint\_file}). For
simplicity, the static detectors are also read from the checkpoint file that
contains the position of all detectors, static and Lagrangian
(\option{/io/detectors/static\_detector/from\_checkpoint\_file}).

Checkpoint filenames all end with [dump
no.]\_checkpoint[-[process]].[extension], where the process number is added
to the mesh, halo and field files in parallel.  The checkpoint detectors
filenames contain \_det after [dump no.]\_checkpoint. 

A script is available at \lstinline[language = bash]+scripts/rename_checkpoint.py+ 
that can be used to easily rename these filenames and there contents to continue 
the naming convention of the original run. For more information see section 
\ref{sec:rename_checkpoint}.

\index{checkpointing!restarting}
To restart from a checkpoint, specify the checkpointed FLML file as input.

\subsubsection{Stat}
\index{stat file}
Contains additional options for handling stat files, for example outputting
a stat file at the start (timestep zero) and outputting stat data before and
after an adapt.

There are further options under individual fields, for example to exclude data from the stat file.

\subsubsection{Detectors}
\index{detectors!options}\label{detectors_options}

Detectors are set in Diamond with the \option{/io/detectors} option. The detectors can be set to be static detectors, \option{/io/detectors/static\_detector}, Lagrangian detectors \\
\option{/io/detectors/lagrangian\_detector} or an array of detectors, \\
\option{/io/detectors/detector\_array}.

When choosing to set detectors using an array, the total number of detectors
needs to be specified in
\option{/io/detectors/detector\_array/number\_of\_detectors} and the type of
detectors is indicated with the option
\option{/io/detectors/detector\_array/static} or \\
\option{/io/detectors/detector\_array/lagrangian}.

\index{Python!detector positions}
Examples \ref{examp:python_function_detectors} and
\ref{examp:python_function_detectors_1} illustrate the use of a Python
function to set an array of detectors, that can be static or Lagrangian.

\begin{example}
  \begin{lstlisting}[language=Python]
def val(t):
            import math

            ret=[]
            for i in range(100):
                    ret.append([-2.5,(-0.495 + i * 0.01)])

            return ret
  \end{lstlisting}
  \caption{A Python function setting 100 detectors. This
  example illustrates that it is possible to use a Python function to set an array of detectors.}
  \label{examp:python_function_detectors}
\end{example}

\begin{example}
  \begin{lstlisting}[language=Python]
def val(t):
            import math

            ret=[]
            for k in range(100,2000,100):
                for j in range(7000,25100,100):
	               for i in range(7000,25100,100):
		             ret.append([i,j,k])

	return ret
  \end{lstlisting}
  \caption{A Python function setting 622459 detectors uniformly distributed
    at intervals of 100 m in the three orthogonal directions. They cover 19 z planes, from z=100 to z=1900, with 32761 detectors in each plane, from
    x=7000 to x=25000 and y=7000 to y=25000.}
  \label{examp:python_function_detectors_1}
\end{example}

If one or more Lagrangian detectors are selected the option \option{/lagrangian\_timestepping} must be set to define detector movement. 
The user can define the order of the Runge-Kutta method to be used by defining the Butcher tableau and timestepping weights under \option{/explicit\_runge\_kutta\_guided\_search}. 
For convenience the first- and fourth-order Runge-Kutta method (\option{/forward\_euler\_guided\_search} and \option{/rk4\_guided\_search}) are available as pre-defined options. 
See section \ref{sec:lagrangian_trajectories} for more information on Lagrangian detector advection.

The output of the detectors is an ascii file called
\option{name\_of\_simulation.detectors} where the name of the simulation has
been indicated in \option{/simulation\_name}. If binary output
\option{/io/detectors/} \option{binary\_output} option is enabled then the file containing the detector data is called
\option{name\_of\_simulation.detectors.dat} and in this case \option{name\_of\_simulation.detectors} contains only the header with the
information about the detectors (name of each detector, in which column the
position of each detector is stored, etc.).

Note that the algorithm used to determine the element containing a detector (of
any kind) assumes that the detector is known to be within the simulation domain.
The option \option{/fail\_outside\_domain} will cause Fluidity to fail if a detector
is found to be outside the domain boundaries. 
This option should be the default choice when creating new simulations.
If detectors are intended to temporarily reside outside of the domain the option
\option{/write\_nan\_outside\_domain} will cause Fluidity to write "NaN" values 
to the detector output in this case.

If \option{/move\_with\_mesh} is selected with any mesh movement algorithm the 
detectors will move according to the mesh displacement. 
That is to say that a static detector will remain static with reference to the 
domain boundaries rather than its true physical coordinates.

\subsubsection{Log output}

Enables additional output to the screen or log file. Usually controlled
using the -v option when running \fluidity. However, a useful option for
logging memory diagnostics can be switched on here.

\subsection{Timestepping}
\index{time!step}
These options control the start and end time of the simulation as well as options regarding timestep size.

\subsubsection{Current time}
This is the model time at the start of the simulation. In most cases this is likely to be zero. It can be non-zero when continuing a simulation from a checkpoint.

\textbf{Time units}

If your simulation contains real data, for example when using
\option{ocean\_forcing}, \fluidity\ must know how to map simulated time onto
real time. This option allows the user to specify the ``real-world'' start
time of the simulation.  The input is a string of the form:
\begin{lstlisting}[language=bash]
seconds since 1992-10-8 15:15:42.5 -6:00 
\end{lstlisting} 

which indicates seconds since October 8th, 1992 at 3 hours, 15 minutes and
42.5 seconds in the afternoon in the time zone which is six hours to the
west of Coordinated Universal Time (i.e. Mountain Daylight Time). The time
zone specification can also be written without a colon using one or
two-digits (indicating hours) or three or four digits (indicating hours and
minutes).


\subsubsection{Timestep}
The simulation timestep. If adaptive timestepping is not used this will
define the size of the timestep used throughout the simulation.  If adaptive
timestepping is used this option defines only the size of the first
timestep.

\subsubsection{Finish time}
The model time at which the simulation should halt. Note that the simulation
may overrun slightly due to roundoff in calculating the current time or if
the timestep does not divide the simulation time exactly.

\subsubsection{Final timestep}

Rather than specify a finish time, the final timestep may be specified. This is the number of timestep after which the simulation will stop.

\subsubsection{CPU time limit}

This option will stop the simulation after the CPU time reaches this limit.
This option is useful when coupled with
\option{/io/checkpointing/checkpoint\_at\_end} enabled.

\subsubsection{Wall time limit}

This option will stop the simulation after the wall time (real time) reaches
this limit. This option is useful when coupled with
\option{/io/checkpointing/checkpoint\_at\_end} enabled.

\subsubsection{Nonlinear iterations}
Nonlinear quantities in the equations are represented by their last known
values. It may be necessary to solve the equations more than once to produce
better approximations to those last known values for reasons of accuracy or
stability. Unless there are reasons for doing this, set this value to 2.

\subsubsection{Adaptive timestep}
\label{section:config_adaptive_timestep}
This option allows the timestep, $\Delta T$, to vary throughout the run, depending on the 
Courant-–Friedrichs-–Lewy (CFL) number.  There are several sub-options here. The
\option{\ldots/requested\_cfl} is the desired upper limit of the CFL. A value of 5-10 is usual
here. \fluidity\ will increase the timestep if the CFL number is less than this value and decrease
it if the CFL is greater than this. The \option{\ldots/minimum\_timestep} and 
\option{\ldots/maximim\_timestep} options limit the timestep. The option 
\option{\ldots/increase\_tolerance} determines
the rate of growth in the timestep. A value of 1.5 indicates the timestep can grow by
at most, 50\%. There is no limit on the rate of decrease. Note that if a timestep fails
to meet the CFL limit imposed it is not re-run, but the timestep is decreased for the next
iteration. Finally, a desired $\Delta T$ can be calculated at the first time iteration by
switching on the option: \option{\ldots/at\_first\_timestep}

\subsubsection{Steady state}
It is possible to run \fluidity until it converges to a steady state; this is sometimes useful for initialising a problem. In order to do this, switch on \option{/timestepping/steady\_state} and set a tolerance. 

\subsection{Physical parameters}
These options control global physical quantities.

\subsubsection{Gravity}\label{sec:Gravity}
\index{gravity}
The importance of buoyancy is discussed in section \ref{sec:buoyancy_hydrostacy}. This
requires a gravitational field to be set and involves both its magnitude
(\eg \mss[9.8]) and a vector field specifying the direction in which gravity
points. For a 3D simulation in a flat domain with gravity pointing in the
negative z direction you would set \verb+value(WholeMesh)+ for this field to
the vector (0.0, 0.0, -1.0). For a gravitational force with spatially
varying direction, e.g. on the Earth considered in Cartesian space with
gravity pointing in the negative radial direction you could use a Python
function of the form
\begin{example}
  \begin{lstlisting}[language=Python]
def val(X, t):
   from math import sqrt
   radius=sqrt(X[0]**2+X[1]**2+X[2]**2)
   rx=X[0]/radius
   ry=X[1]/radius
   rz=X[2]/radius
   return (-rx, -ry, -rz)
  \end{lstlisting}
  \caption{A Python function returning a vector pointing in the negative radial direction.}
\end{example}

\subsubsection{Coriolis}
\index{Coriolis!options}

\fluidity\ supports the specification of the Coriolis term (section
\ref{sec:coriolis}) in a number of different ways. The following options
are available:

\begin{enumerate}
  \item \option{f\_plane} -- a single float is prescribed which corresponds to
        $f_0$ in \eqref{eq:f-plane};
  \item \option{beta\_plane} -- here two floats are prescribed, $f_0$ and
        $\beta$ in \eqref{eq:beta-plane};
  \item \option{sine\_of\_latitude} -- here the Coriolis parameter from
        \eqref{eq:coriolis_parameters} is used and $\Omega$, 
        $R_{\textrm{earth}}$ and $\textrm{latitude}_0$ are defined as floats with
        latitude calculated via 
        $\phi = y/R_{earth} + \textrm{latitude}_0$;
  \item \option{on\_sphere} -- here $\Omega$ the rotation vector pointing in the
        inertial frame $z$ direction  \eqref{eq:on_sphere_rotation} is set,
        note this is the direction pointing from the centre of mass to the North
        Pole on the Earth;
 \item \option{python\_f\_plane} -- time dependent python input
       prescribing a single float which corresponds to
       $f_0$ in \eqref{eq:f-plane} - see example \ref{ex:python_f_plane}.
\end{enumerate}

Recall that there is a factor $2$ relationship between $f$ and $\Omega$
\eqref{eq:f_omega} --- make sure you don't get caught out by this.

\begin{example}
\begin{lstlisting}[language = Python]
if t < 4000.0:
  omega = 3.0
elif t < 6000.0:
  omega = 2.5
elif t < 8000.0:
  omega = 2.0
elif t < 10000.0:
  omega = 1.5
elif t < 12000.0:
  omega = 1.0
elif t < 14000.0:
  omega = 0.5
else:
  omega = 0.0

return 2.0 * omega
\end{lstlisting}
\caption{\option{python\_f\_plane} definition, sweeping through a number of
         rotation rates. Note the factor of $2$ between $f$ and $\Omega$ (see
         equation \eqref{eq:f_omega}).}
\label{ex:python_f_plane}
\end{example}


\section{Meshes}\label{sec:Mesh}

A mesh defines the discrete
function space in which the values of one or more fields lie. For example, the mesh
defines what degree of polynomials are employed on each element, whether the
field is continuous or discontinuous between elements, and whether the
domain is periodic in any direction. 

Meshes are defined in the flml file by \option{/geometry/mesh} options. The
mesh associated with each field is referred to by name in the
\option{\ldots/mesh} option under that field.

\subsection{Reading meshes from file}
\index{mesh!input}
There must always be one mesh which is read in from a file in
Gmsh format. This is usually the \option{CoordinateMesh}. To specify the
Gmsh file from which the coordinate mesh should be read, set the
\option{file\_name} attribute of
\option{/geometry/dimension/mesh::CoordinateMesh/from\_file} to the basename
of the Gmsh file (that is, the filename without .msh)

The coordinate mesh read in from file will always have linear elements and
the Coordinate field is always continuous between elements.

\fluidity\ also has native Gmsh support, which loads in Gmsh files directly into
\fluidity, and works with binary and ASCII Gmsh formats. To enable native 
support, \fluidity\ needs to be told to expect a Gmsh file, which is achieved 
by setting the \onlypdf\option{/geometry/mesh/from\_file/format}\ option 
to \onlypdf\option{gmsh}.  \fluidity\ will now look for a file with the extension 
\lstinline[language=bash]+.msh+ when it runs.

For information on generating meshes in Gmsh format, see chapter
\ref{chap:meshes}. 

\subsection{Deriving meshes from other meshes}
\index{mesh!derived}
The alternative to reading a mesh from a file is to derive it from another
mesh. This is necessary when for instance we wish to derive a mesh with
different continuity or elements than the original mesh. 
For example, if we have a \option{CoordinateMesh} as our input mesh read 
from file, it is possible to derive a \option{VelocityMesh} from it
by adding \option{/geometry/mesh::VelocityMesh}, selecting 
\option{from\_mesh} under it and there selecting
\option{mesh::CoordinateMesh}. If nothing further is specified 
under the new mesh, the derived mesh will be exactly the same as the mesh
it is derived from.

The more interesting case occurs where we wish to derive a mesh with
different continuity or elements from the original mesh. To specify a
discontinuous mesh, under \option{\ldots/from\_mesh} enable
\option{mesh\_continuity} and select \option{discontinuous}.

Similarly, to specify a mesh with higher polynomial degree elements, enable
\option{mesh\_shape} under \option{\ldots/from\_mesh} and set
\option{polynomial\_degree}.

Meshes with any name can be added. Only the name \option{CoordinateMesh} 
is special, as it will be used to store the coordinates of the mesh. This 
is also the only required mesh. \option{VelocityMesh} and \option{PressureMesh}
are only provided as suggested names as quite often the pressure 
and velocity fields need to be on a different mesh 
than the coordinate mesh, e.g. for \PoDGPt. It is however 
not required that the velocity is defined on a mesh with the name
\option{VelocityMesh}, nor does the pressure field have to be on a mesh
with the name \option{PressureMesh}. If for instance the mesh needed 
for pressure is the same as your \option{CoordinateMesh}, e.g. 
for \Poo or \Pzero\Pone, the pressure can be defined directly on
the \option{CoordinateMesh} and no extra mesh is needed.

\subsubsection{Shape function}

This option is used to specify the degree of polynomial which should be used
for the shape functions on each element of the mesh. If not selected, the
shape functions will be the same as those on the mesh from which this mesh
is derived.

\subsubsection{Continuity}

This option can be set to discontinuous to derive a discontinuous mesh from
a continuous one. Note that it is not possible to derive a continuous mesh
from a discontinuous mesh.

\subsubsection{Periodic}\label{sec:periodic}
\index{mesh!periodic} 
\index{periodic domain} 

To specify a periodic domain in Diamond, add a new mesh under
\option{/geometry/mesh}.  Select \onlypdf\linebreak
\option{\ldots/from\_mesh/mesh::CoordinateMesh} and then turn
on \onlypdf\\
\option{\ldots/from\_mesh/periodic\_boundary\_conditions} for each dimension
that is periodic. There are three fields that need to be completed: the
surface IDs of one side, the surface IDs of the opposite side and a python
function (\option{coordinate\_map}) which contains the necessary mapping
function.

For example, suppose that the domain is the unit square shown in figure
\ref{fig:periodic}\ which is to be periodic in the $x$ direction. We
designate surface ID 1 as the \emph{physical boundary}\ and surface ID 2 as the
\emph{aliased boundary}. We therefore enter 1 in
\option{\ldots/physical\_boundary\_ids}\ and 2 in
\option{\ldots/aliased\_boundary\_ids}. The \emph{coordinate map}\ function
takes a point on the \emph{aliased}\ boundary to the corresponding point on
the \emph{physical}\ boundary. In this case, the appropriate function is:
\begin{lstlisting}[language=Python]
def val(X,t):
    result = list(X)
    result[0]=result[0]-1.0
    return result
\end{lstlisting}

\begin{figure}[ht]
  \centering
  \xfig{configuring_fluidity_images/periodic_domain}
  \caption{Periodic unit square with surface IDs 1-4 shown.}
  \label{fig:periodic}
\end{figure}

Meshes that are required to be periodic can now be derived from this
periodic mesh. Note that the periodic mesh must be directly derived from the
mesh which has been read \option{from\_file}. It is not possible to derive a
periodic mesh from a \option{from\_mesh}\ mesh.

\subsubsection{Extruded meshes}\label{sec:extruded}

It can be advantageous to have a mesh in which all the nodes line up in
vertical lines. To achieve this effect within \fluidity, it is possible to
read in a mesh in $n-1$ dimensions and extrude it along the $n$-th
dimension. An extruded mesh is specified using the
\option{\ldots/from\_mesh/extrude}\ option. 

Under this option it is necessary to set the \option{regions/bottom\_depth}. This is
a scalar value which gives the depth, a positive value. The extent of the
domain in $n$-th dimension will be $(0,-bottom\_depth)$. This value may be
set either as a constant or as a Python function. In the latter case,
function will be a function of space and time. Note in this case that the
space argument \lstinline[language=Python]+X+ will be $n-1$-dimensional. See
section \ref{sec:setting_with_python}\ for a full explanation of the use of Python functions to
prescribe field values. In this case, the depth is essentially a scalar
field over the $n-1$ dimensional parent mesh. The time argument which
will be passed to the function is the simulation start time and the function
will \emph{not}\ be re-evaluated during the simulation.

The second option which must be set is the
\option{\ldots/sizing\_function}. This specifies the mesh spacing
along the $n$-th dimension. It may once again be a constant or a Python
function. In the latter case, it will be a function of all $n$ dimensions
which facilitates the mesh spacing varying in depth as well as in the
horizontal. Once again, the function will be evaluated only at simulation
start. 

It will usually be advantageous to specify the surface ID to be associated
with the top and bottom boundaries, so that boundary conditions can be
associated with them. This is achieved using the
\option{\ldots/top\_surface\_id}\ and
\option{\ldots/bottom\_surface\_id}\ options. The lateral boundaries
of the extruded mesh will inherit the surface IDs associated with the the
boundaries of the parent (non-extruded) mesh.

It is possible to specify different options for different regions of the
mesh by adding multiple \option{\ldots/extrude/regions}\ options and changing
them to the generic rather than the \option{regions::WholeMesh}\
version. The new regions need to be named and the region IDs to which they
apply specified.

As well as extruding 2D meshes (where only the $x$ and $y$ coordinates are
specified in the Gmsh file), given a pseudo 2D mesh on a spherical
shell, \fluidity\ can perform and extrusion in the radial direction. To
perform such an extrusion simply enable the options as noted above and
additionally check the \option{/geometry/spherical\_earth option}. With this
option enabled \fluidity\ will then perform the specified extrusion towards
the centre of the sphere.

Extrusions on the sphere can be performed such that the `depth' of the extrusion
conforms to bathymetric data. To extrude according
to bathymetry \fluidity\ must be provided with a netCDF data file containing
three columns of data giving the longitude, latitude and depth
of each point respectively. The name and location of this data file must then be
entered under \option{\ldots/extrude/regions/bottom\_depth/from\_map}.
To avoid the depth at coast dropping to zero the user may also enter a minimum
depth under the \option{\ldots/from\_map/min\_depth} option. Note however, that
if a minimum depth is specified, this will be applied throughout the domain.

\subsubsection{Extruded periodic meshes}\label{sec:extrudedperiodic}

If an extruded periodic mesh is required then the periodic mesh must first
be derived from the \option{from\_file}\ mesh. The extruded mesh is then
derived from the periodic mesh. All other meshes are next derived from the
extruded mesh. A special case is the \option{CoordinateMesh}. This must be
derived from the extruded mesh by specifying
\option{periodic\_boundary\_conditions/remove\_periodicity}. At this stage it
is necessary to re-specify the \option{physical\_boundary\_ids},
\option{aliased\_boundary\_ids} and \option{coordinate\_map}. Additionally, the
\option{inverse\_coordinate\_map}\ must be given. As the name suggests, this
function must be the inverse of the original \option{coordinate\_map}.


\section{Material/Phase}
The final compulsory element in the top level of the options tree is
\option{/material\_phase}.  A \option{/material\_phase} element groups all
of the fields which pertain to one phase or one material. See
section~\ref{sec:config_multimatph} for an explanation of the distinction
between a phase and material in this context.

When configuring ocean problems (or single material/single phase fluids
problems), only one \option{/material\_phase} is required.  Multi-material
and multiphase problems will require one \option{/material\_phase} for each
phase or material in the problem.

Note that you must give each of your material phases a name.

% NOTE: if you add anything below the Sediments sections here, change the ref...
The following sections (\ref{config:spatial} to \ref{config:sediments}) describe the options below the
\option{/material\_phase} option.

\section{Fields}
\index{field}
\subsection{Types of field}

A field associates a value with every node in the domain. Examples of fields
in a fluids simulation include the velocity and pressure. Fields in \fluidity\
are distinguished by the rank of the data on the field and the way in which
that field is calculated. 

\begin{description}
\item[Scalar fields] have a scalar value at each node. Common examples
  include temperature, pressure and density.
\item[Vector fields] have a vector value, in other words a list of numbers,
  at each node. The rank of a vector field is 1 and the length of the
  vector is given by the dimension of the problem.
\item[Tensor fields] have a value given by a square matrix at each
  node. The side length of the matrix is the problem dimension and the rank
  is naturally 2. The diffusivity of a tracer is a typical example of a
  tensor-valued field.
\end{description}

Fields can also be characterised by the manner in which their value is
calculated. \fluidity\ recognises three such categories:

\begin{description}
\item[Prognostic fields] are the result of solving a partial differential
  equation. In a typical fluids simulation, the velocity and pressure are
  prognostic and are calculated by solving some variant of the Navier-Stokes
  equations. Similarly, tracers such as temperature and salinity are usually
  the result of solving an advection-diffusion equation. Prognostic fields
  typically have specified initial and boundary conditions and it will be
  necessary to specify spatial and temporal discretisation options. If an
  implicit timestepping scheme is in use (and it almost always is in
  \fluidity), it is also necessary to specify solver options. 
\item[Diagnostic fields] are calculated from other fields without solving a
  partial differential equation. A typical example is the CFL number which
  may be calculated from the timestep, the mesh spacing and the velocity
  field. 
\item[Prescribed fields] receive their values from sources external to
  \fluidity. This might be a constant or varying function specified by the
  user, or it might be interpolated from some external data set. Fields such
  as diffusivity and viscosity are often prescribed as are source and
  absorption terms.
\end{description}

An additional field type - aliased - is also available.  This links the values in one field to those in another, using no extra computational resources during the simulation (i.e. it is not an independent field).  This is useful when sharing fields between material\_phases.  For example if two material\_phases share a common velocity field then only one should contain a prognostic field while the other is aliased to the other material\_phase.

\subsection{Setting field values}\label{sec:setting_field_values}
\index{field!values}
\index{initial conditions!setting}
Field values must be specified by the user in two circumstances: the initial
value of most prognostic fields and the value throughout the simulation of
all prescribed fields. 

The initial value of prognostic fields is set with the
\option{\ldots/prognostic/initial\_condition} option while the value of
prescribed fields is set with the \option{\ldots/prescribed/value} option.


\subsubsection{Constant fields}
\index{field!constant}
Fields which are constant in space and (for prescribed fields) time may be
specified by simply providing a constant value in the \option{constant}
option under \option{\ldots/prognostic/initial\_condition},
\option{\ldots/prescribed/value}. For a scalar field, this is a single
floating point (real) value while for a vector field this is a list of reals
of length equal to the field dimension.

For tensor valued fields there are more options. It is possible to specify
an isotropic, or rotation invariant, value by choosing
\option{\ldots/value/isotropic/constant} and specifying a single real which
will be used for all the diagonal entries of the tensor field at all mesh
nodes. The off-diagonal entries of an isotropic tensor field are always
zero. A constant anisotropic field may be specified by choosing
\option{\ldots/value/anisotropic\_asymmetric/constant} and providing the
entire matrix. Finally, a constant symmetric anisotropic tensor field may be
specified by selecting \onlypdf\linebreak
\option{\ldots/value/anisotropic\_symmetric/constant}. In this case, the
user must specify all of the entries in the upper half of the matrix and
those in the lower half will be filled automatically by symmetry.

\subsubsection{Setting fields with a Python function}\label{sec:setting_with_python}
\index{field!Python function}
\index{Python!prescribed field values}
The value of a field which varies in space and (for prescribed fields) in
time may be specified by providing an appropriate function written in
Python. The Python function will be evaluated for each node in the mesh at
the start of the simulation to populate the field values. For time-varying
prescribed fields, the function will be evaluated again at the beginning of
every timestep to update the field value. If it is known that the value of
the field does not in fact vary in time, then the re-evaluation of the
Python function on each timestep can be inhibited by setting the
\onlypdf\linebreak \option{\ldots/prescribed/do\_not\_recalculate} option.

The Python function must be provided as the value of the
\option{\ldots/python} option which may be chosen as an alternative to the
\option{\ldots/constant} function. The option may contain any Python code
but it must define a function \lstinline[language=Python]+val(X,t)+ where
the sequence \lstinline[language=Python]+X+ is the coordinates of the point
at which the field is being evaluated and \lstinline[language=Python]+t+ is
the current time. 

For a scalar field, the function must return a single floating point
value. Similarly for a vector field, a sequence of values of length equal to
the field dimension must be returned. For a tensor field, there are two
cases. For an isotropic field specified with
\option{\ldots/value/isotropic/python}, the function must return a single
float which will be used for all the diagonal entries of the tensor at that
point. The off-diagonal entries will be set to zero. For the anisotropic case,
the function must return a two-dimensional array (a sequence of
sequences). It is the user's responsibility to ensure that the tensor is
symmetric in cases where it should be.

\begin{example}
  \begin{lstlisting}[language=Python]
def val(X,t):
    return (-X[1],X[0])
  \end{lstlisting}
  \caption{A Python function returning a two-dimensional solid rotating
    vector field about the origin.}
\end{example}

\subsubsection{Reading fields from a file (using the \option{from\_file} option)}
\index{field!input}
A field can be populated using saved data from a file. This is intended primarily
for picking up prescribed fields from previously run prognostic simulations
(checkpointing) and may be specified by providing the file name in
the \option{from\_file} option under \option{\ldots/prognostic/initial\_condition},
\option{\ldots/prescribed/value}.

For prescribed fields the format of the input file containing field data must be
vtu, and this will only work for those prescribed fields directly underneath
\option{\ldots/material\_phase}. For prognostic fields it is possible to select
the type of input file, under \option{\ldots/initial\_condition/from\_file/format};
the available supported formats for this include \option{vtu} and \option{NetCDF-CF 1.x}.

The file mesh must match the mesh of this field (except for piecewise constant
fields which will be remapped back from the discontinuous nodal values). In
parallel the process number is appended to the filename, e.g. if the file name
attribute is set to \option{input.vtu}, process 0 reads from \option{input-0.vtu}.

\subsubsection{Setting an initial condition from a NetCDF file}\label{sec:setting_from_netcdf}

The initial state of certain fields can be set from external data contained within a NetCDF file.
This functionality can be used by selecting the \option{\ldots/initial\_condition/from\_netcdf/format}
option.
This option will not currently work with multi-layered data files. Supported
NetCDF file conventions include the NetCDF-CF 1.x convention.

\subsubsection{Setting fields from NEMO data}\label{sec:setting_from_nemo}
\index{field!from Nemo}
Initial conditions of prognostic fields and the values of prescribed fields can also be set from an external NEMO
data file. The external data file is in the NETCDF format and data is currently available for pressure, temperature, salinity
and velocity. To set the initial condition of a prognostic field from NEMO data, choose the option 
\option{\ldots/prognostic/initial\_condition/NEMO\_data} and then under \option{format} select the required data format. For scalar fields
the formats available are `Temperature', `Salinity' and `Free-surface height', for vector fields `Velocity' is the only available format.
Setting the value of prescribed fields from NEMO data works similarly. Set the option \option{\ldots/prescribed/value/NEMO\_data} and then proceed as above.

\subsubsection{Setting an initial free surface height}\label{sec:setting_free_surface_height}
\index{inital conditions!free surface height}
\index{free surface!initial condition}
The free surface height is contained within the Pressure field.  To apply an initial condition on free surface height, 
choose the \option{\ldots/free\_surface} under the relevant Pressure initial condition option.
With this option, it is possible to set the initial free surface elevation in a tsunami simulation, for example.
The initial condition can be applied using the approaches outlined above, in this section~\ref{sec:setting_field_values}.
It is recommended that a diagnostic FreeSurface field is included if this option is used.

If set from a NetCDF file using the option \option{\ldots/initial\_condition/from\_netcdf}, and the file provides exactly the
free surface height, it is important that the child option \option{\ldots/initial\_condition/from\_netcdf/format} is set to `raw'.

If the \option{no\_normal\_stress} option is used, to impose 
boundary condition \eqref{NormalStressBC} instead of $p=\patm$, you should add a
prognostic FreeSurface (instead of a diagnostic). The initial condition is then
also set for that field.

\subsection{Region IDs}
\index{region ID}
If the input mesh defines a number of region IDs then these may be employed
to specify different field values for each region. For a prescribed field,
this is achieved by changing the \option{\ldots/value::WholeMesh} element to
the unnamed \option{\ldots/value}. The user must then specify a new name for
that value element. Next, enable the  \option{\ldots/value/region\_ids}
option and set it to a list of region ids to which this value should
apply. Any number of  \option{\ldots/value} elements may be added to allow
different values to be specified in different regions. For prognostic
fields, analogous behaviour for initial conditions may be achieved by
switching \option{\ldots/initial\_condition} from \option{WholeMesh} to a
user-specified name and specifying the region IDs appropriately.

See section \ref{sec:region_ids}\ for information on including region IDs
in meshes.

\subsection{Mathematical constraints on initial conditions}\label{sec:ICs}

For well-posedness, the initial condition of the velocity field must
satisfy both continuity (\ref{conty}) and the boundary conditions
imposed on the problem. If the normal component of velocity is
imposed on the entire boundary then the additional
compatibility constraint of global mass conservation must be
satisfied:
\begin{equation*}
\int_{\partial\Omega}\bmn\cdot\bmu=0.
\end{equation*}

\section{Advected quantities: momentum and tracers}
\label{config:spatial}

\subsection{Spatial discretisations}\label{sec:Spatial discretisations}
\index{advection-diffusion equation!discretisation options}
\index{momentum equation!discretisation options}

A number of underlying finite element schemes are available for tracer
fields and velocity. In each case there are restrictions on the mesh
continuity which must be employed. In addition, the \onlypdf\linebreak
\option{conservative\_advection} option is applicable to all discretisation
types. See chapter \ref{chap:numerical_discretisation} for details.

For each field, the spatial discretisations can be selected using
\option{\ldots/prognostic/spatial\_discretisation}.  Once selected a number
of other options will open underneath this option.

\subsubsection{Continuous Galerkin}
\label{sec:configuring_fluidity_continuous_galerkin}

Continuous Galerkin (CG) implements the CG scheme detailed in \ref{sec:balancing_diffusion}. If stabilisation methods are needed, the user can either select streamline upwind or streamline upwind Petrov-Galerkin. Other options include integration of advection by part, lumping of the mass matrix, or direct exclusion of both advection and mass terms.

\subsubsection{Control Volumes}\label{sec:CVs}

The control volume options (\option{control\_volumes}) implements the advection scheme described in section \ref{ControlVolumeAdvection}. There are 
several options to control the face value and how diffusion is implemented.

The face value (\option{face\_value}) can be set to one of:
\begin{description}
\item[FirstOrderUpwind] - see section \ref{sec:cv_fou}.  Note that first order upwinding does not require nonlinear advection iterations as the low order pivot solution uses first order upwinding itself.  However in this case it is necessary that the implicitness factor, $\theta$, is the same as the pivot implicitness factor, $\theta_p$ (see sections \ref{sec:cvtemp} and \ref{sec:configuring_fluidity_temporal_discretisation}).
\item[Trapezoidal] - see section \ref{sec:trap}, should be used with the suboptions describing a face value limiter (see section \ref{sec:cvlimiting}) active
\item[FiniteElement] - see section \ref{sec:cvfe}, should be used with the suboptions describing a face value limiter (see section \ref{sec:cvlimiting}) active
\item[FirstOrderDownwinding] - see section \ref{sec:fod}, intended for demonstration purposes only, not recommended for general use (unconditionally unstable)
\item[HyperC] - see section \ref{sec:hyperc}
\item[UltraC] - see section \ref{sec:ultrac}
\item[PotentialUltraC] - see section \ref{sec:potultrac}
\item[None] - turns off the advective terms
\end{description}

For the diffusion scheme (\option{diffusion\_scheme}) one can choose either:
\begin{description}
\item[ElementGradient] - see section \ref{sec:cvegdiff}
\item[BassiRebay] - works in two configurations equal order field and diffusivity or, for fields on a linear parent mesh, with a piecewise constant (element centred) diffusivity (staggered finite volumes), see section \ref{sec:cvbrdiff}
\end{description}

For steady state problems the mass terms may be disabled using \option{\ldots/mass\_terms/exclude\_mass\_terms}.  Note that this also requires the suitable setting of the temporal discretisation options ($\theta=1$).

\subsubsection{Coupled CV}\label{sec:CoupledCVs}

The coupled CV options (\option{coupled\_cv}) implement another control volume discretisation with face value limits enforced in such a way to give boundedness both in the field and across the sums of fields.  Section \ref{sec:coupledlimiter} contains more details on this method.

Options must be selected to describe the face value scheme, which include most of the algorithms described in section \ref{sec:CVs}.  Additionally it is necessary to prescribe the maximum and minimum bounds on the sum of this and the previous fields.  Because the coupled scheme depends on a priority ordering a priority (outside of the spatial discretisation options at \option{\ldots/scalar\_field/prognostic/priority}) must also be set with higher values having the highest priority and lower values the lowest.

Related fields to be used together during coupled limiting are grouped together based on their names from successive material\_phases.  For example, if a field called MaterialVolumeFraction has coupled\_cv options then all other fields in all other material\_phases called MaterialVolumeFraction using coupled\_cv options will be advected together in order of their priority.  Spatial discretisation options within coupled\_cv may vary between the fields but temporal discretisation options must be identical.

\subsubsection{Discontinuous Galerkin method for the
  advection-diffusion equation}

The Discontinuous Galerkin option implements the advection-diffusion
algorithm described in Section
\ref{sec:ND_discontinuous_galerkin_advection}. There are two
compulsory options to set, as well as a number of non-compulsory options.

\paragraph{Advection scheme (\option{advection\_scheme})} This
\emph{compulsory option} selects the approximation for the flux of
scalar across a face. Select one from:
\begin{itemize}
\item \option{upwind}: Use the upwind flux as described in Section
  \ref{sec:ND_discontinuous_galerkin_advection}. This is the
  \emph{recommended flux} for DG advection.
\item \option{lax\_friedrichs}: Use the Lax-Friedrichs flux as described in
  Section \ref{sec:ND_discontinuous_galerkin_advection}. This is an
  attempt to produce a bounded flux when the advecting velocity is
  discontinuous. This option is only for testing, if you have a
  discontinuous advecting velocity it is recommended to use the upwind
  flux combined with the option to project the velocity to a
  continuous space described below.
\item \option{none}: This option switches off the advection term completely.
\end{itemize}

\option{project\_velocity\_to\_continuous} \\
\option{integrate\_advection\_by\_parts} \\
\option{integrate\_conservation\_term\_by\_parts} 

\paragraph{Diffusion scheme (\option{diffusion\_scheme})} 
This \emph{compulsory option} selects the discretisation method used
for the diffusivity term. This selection is important for performance
since various different options have different stencil sizes, which
affects memory signature and hence the number of elements you can use
per processor.

Select one from:
\begin{itemize}
\item \option{bassi\_rebay}: The classical scheme of Bassi and Rebay (see
  section \ref{BassiRebay}). This scheme results in a large stencil
  for the diffusion matrix, which can reduce computational speed and
  increase memory use. If possible one should use a different option
  with a smaller stencil.
\item \option{compact\_discontinuous\_galerkin}: The compact
  discontinuous Galerkin scheme (CDG) from Peraire and Persson
  \citep{peraire2008} (see section \ref{CDG}). This scheme has the
  smallest stencil of any diffusion scheme and hence is the most
  efficient and uses the least memory resource. \emph{Recommended
    option}.

  Optionally, it is possible to set the
  \option{penalty\_parameter}.\index{penalty parameter!CDG} This
  optional option adds an extra term which penalises jumps across
  faces. You must supply a multiplicative constant which is scale
  independent (typical value is 10). This term is required to prove
  theoretical results about the CDG scheme, but we experimentally
  observe that it is not necessary, and hence, it is \emph{recommended
    not to use this option}.
\item \option{interior\_penalty}: Symmetric interior penalty (IP)
  scheme. This scheme simply integrates the diffusion operator by
  parts in each element, averages the fluxes and adds a term which
  penalises jumps. You must set the \option{penalty\_parameter}
  \index{penalty parameter!interior penalty method} which sets the
  multiplicative constant, and the \option{edge\_length\_parameter}.
  which specifies the scaling with the edge-length $h$. You must also
  select an \option{edge\_length\_option} which is either
  \option{use\_face\_integral} which computes a length scale from the
  face integral, or \option{use\_element\_centres} which uses the
  distance between centres of the two elements on either side of the
  face. Both of these options only function well for nearly isotropic
  meshes and hence CDG is the recommended diffusion choice since it 
  is compact and requires no such parameters.
\end{itemize}

\paragraph{Slope limiter (\option{slope\_limiter})}
Need to mention about subcycling.
\paragraph{Mass terms (\option{mass\_terms})}

\subsubsection{Conservative advection}

The momentum equation can be discretised conservatively by setting the
BETA factor equal to 1 (corresponding to a divergence form of the
equation). If BETA is set to zero, the discretisation is left
non-conservative. An intermediate value can alternatively be
selected. Please refer to section \ref{sec:ND_momentum_equation} for a comprehensive discussion on the influence of this parameter.

\subsubsection{Buoyancy}
\label{sec:configuring_fluidity!spatial_discretisation!buoyancy}
\index{spherical earth}

When solving for momentum on a spherical or cylindrical domain it is possible to specify that the gravity direction should be automatically
calculated to be radial (towards the origin) and evaluated directly at the gauss points using the
\option{buoyancy/radial\_gravity\_direction\_at\_gauss\_points} option (underneath Velocity \option{spatial\_discretisation}).  This
ignores the direction specified in \option{/physical\_parameters/gravity} but still uses the magnitude specified there.  Note
that the direction is only ignored in the buoyancy term of the momentum equation and the version specified under
\option{/physical\_parameters/gravity} may be used elsewhere in the model.

Evaluating the gravity direction directly at gauss points, rather than interpolating it from the values specified at the nodes,
avoids interpolation errors that introduce non-physical imbalances between buoyancy and pressure in spherical simulations.  As this
is mostly used for simulations on N-spheres it is generally combined with the \option{/geometry/spherical\_earth} option (see
\ref{sec:spherical_earth}).

\subsection{Temporal discretisations}
\label{sec:configuring_fluidity_temporal_discretisation}
Under temporal discretisation, you can set the value of theta, where $0$ is explicit, $0.5$ is Crank-Nicolson and $1$ is implicit. For scalar fields, the control volumes option may be selected if you are using control volumes or coupled cv spatial discretisation.  It contains options to set up nonlinear advection iterations, subcycling and the value of the pivot implicitness factor (see section \ref{sec:cvtemp}). The discontinuous Galerkin option can be used if you are using discontinuous galerkin spatial discretisation to set the maximum Courant number per subcycle, or the number of subcyles. 

\subsection{Source and absorption terms}\label{sec:Source}
\index{source term}
\index{absorption term}

The source and absorption terms allow for external forcing of the tracer and
momentum equations. The source is a rate of change of the tracer which is
independent of the system state while the absorption term is linear in the
tracer. The source and absorption terms modify the tracer
equation as follows (cf. \eqref{eq:general_scalar_eqn}):
\begin{equation}
  \frac{\partial c}{\partial t} = F(c,u,t) - \sigma c + F,
\end{equation}
where $F(c,u,t)$ represents the advection and diffusion terms, $\sigma$ is the
absorption and $F$ is the source. The source and absorption are usually
prescribed fields supplied by the user but in some cases it may be necessary
to provide a diagnostic field which will be set by a parameterisation
somewhere in the model. If this is the case then this will be specified in
the documentation of that parameterisation.

For tracer fields, the source and absorption are specified by the\onlypdf\\
\option{\ldots/scalar\_field/prognostic/scalar\_field::Source} and\onlypdf\\
\option{\ldots/scalar\_field/prognostic/scalar\_field::Absorption} options
respectively. For velocity, the corresponding fields are naturally
vector-valued and are set by options\onlypdf\linebreak
\option{\ldots/vector\_field::Velocity/prognostic/vector\_field::Source} and\onlypdf\\
\option{\ldots/vector\_field::Velocity/prognostic/vector\_field::Absorption}
respectively.

\subsection{Sponge regions}\label{sec:Sponge}
\index{sponge regions}
It is often useful to be able to relax momentum or a field variable to a
given state in regions of the domain, typically in regions close to
boundaries. This may be done using a combination of source and absorption
terms. If $F$ is the value of the source at a point and $\sigma$ is the
absorption, then the value of the scalar field $c$ will tend to relax to a value
of $F/\sigma$. The absorption, $\sigma$, has units $1/\mathrm{time}$ and controls the
time over which the relaxation occurs.

\section{Solving for pressure}
\index{pressure!options}
\label{sec:configuring_fluidity_pressure}

\subsection{Geostrophic pressure solvers}
\index{pressure!geostrophic balance}
\label{sec:config_geostrophic_balance}

\subsection{First guess for Poisson pressure equation} \label{sec:poisson_pressure_solution}
\fluidity's solution procedure for velocity and pressure can use a pressure Poisson guess to speed up the convergence. In order to use a pressure guess, set \option{\ldots/scalar\_field::Pressure/prognostic/scheme/poisson\_pressure\_solution} from \option{never} to \option{only\_first\_timestep}.


\subsection{Removing the null space of the pressure gradient operator} \label{Nullspaceremove}
\index{pressure!null space}

If the normal component of velocity is imposed on all boundaries then the
appropriate boundary condition for pressure 
\citep[see][]{gresho87} is obtained by taking the normal component of
(\ref{mtm}), this yielding a Neumann boundary condition for
pressure. This only serves to define the pressure field up to an
arbitrary additive constant.

There are two different and mutually exclusive options which may be used to
fix the additive constant in the pressure field. The first is that the
pressure at a single point may be set to 0. This is achieved by setting the
\option{\ldots/scalar\_field::Pressure/prognostic/reference\_node} option. The
value of the option is the number of a node at which the pressure is to be
constrained to be zero. It is an error for the node number specified here to
be greater than the number of nodes in the simulation.

The second method is to instruct the linear solver to remove the null space
of the pressure equation as a part of the solution procedure. This is
achieved by enabling the\linebreak
\option{\ldots/scalar\_field::Pressure/prognostic/solver/remove\_null\_space}
option. This approach often leads to better convergence rates than setting
the \option{reference\_node}.

If however there is a single location on the boundary where the normal
component of velocity is not specified then there is no free constant in the
pressure and neither \option{\ldots/reference\_node} nor
\option{\ldots/remove\_null\_space} should be set. An example
may be stress free outflow or the presence of a free surface.

\subsection{Continuous Galerkin pressure with control volume tested continuity}
\index{pressure!CG with CV tested continuity}
\label{sec:config_cg_pressure_cv_continuity}

As described in section \ref{sec:cg_pressure_cv_continuity} when using a continuous Galerkin 
discretisation of pressure the continuity equation can be tested with the corresponding 
dual control volume mesh. This is achieved by including the \linebreak 
\option{\ldots/scalar\_field::Pressure/prognostic/spatial\_discretisation/continuous\_galerkin/test\_continuity\_with\_cv\_dual}
option. As described in the theory section \ref{sec:cg_pressure_cv_continuity} this will imply a 
non symmetric pressure correction matrix which must be considered when selecting the pressure 
matrix solver options.

The current limitations of this method are:
\begin{enumerate}
  \item It can only be used for incompressible flow \ref{sec:ND_pressure_equation}.
  \item It cannot be used with the standard $p=0$ free surface model (but can be used with the zero normal stress free surface) \ref{sec:free_surface}.
  \item It cannot be used with the wetting and drying model \ref{sec:wetting_and_drying}.
  \item It can only be used if the pressure has a mesh associated with Lagrangian shape functions.
  \item It can only be used with control volume shape functions that are available, of which only P1CV are considered reliable. 
\end{enumerate}

\section{Solution of linear systems}\label{sec:Solve}
\index{linear solvers!options}

\subsection{Iterative Method}
\index{GMRES}
\index{conjugate gradient}
\index{multigrid}
\index{PETSc}
As described in Section \ref{ND_Linear_solvers}, for the solution of large sparse linear systems, the so called iterative methods are usually employed. These methods avoid having to explicitly construct the inverse of the matrix, which is generally dense and therefore costly to compute (both in memory and computer time). FLUIDITY is linked to PETSc: a suite of data structures and routines for the scalable (parallel) solution of scientific applications modelled by partial differential equations.  It employs the MPI standard for parallelism. FLUIDITY therefore supports any iterative method provided by the PETSc library (http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/docs/manualpages/KSP/KSPType.html --- available methods may depend on the PETSc library installed on your system). Examples include Conjugate Gradient (CG), GMRES and FGMRES (Flexible GMRES). Some options are explicitly listed under \option{solver/iterative\_method}, for example CG: \option{solver/iterative\_method::cg}, whereas 
others can be selected entering the name of the chosen scheme in \option{solver/iterative\_method}.

\subsection{Preconditioner}
\index{linear solvers!preconditioners}

The requirement for a suitable preconditioner is described in Section \ref{ND_Preconditioners}. In a manner analogous to the selection of the iterative method, some common preconditioning options are explicitly listed under \option{solver/preconditioner}, for example MG: \option{solver/preconditioner::mg}, whereas others can be selected by entering the name of the chosen scheme in \option{solver/preconditioner}. 

\subsubsection{Direct Solve}

Note that the option to solve a system exactly is available in FLUIDITY. For this, \option{solver/iterative\_method::preonly} must be selected (preonly: preconditioner only) and the preconditioner must be set to \option{solver/preconditioner::LU}. A full LU decomposition of the system is then carried out.

\subsection{Relative Error}
\index{linear solvers!convergence criteria}

The solver finishes if the preconditioned error becomes smaller than the original preconditioned error times this value. 

\subsection{Absolute Error}

The solver finishes if the preconditioned error becomes smaller than this value.

\subsection{Max Iterations}

The maximum number of iterations allowed for the linear solver before quitting.

\subsection{Start from Zero}

Switch on to start a solve with a zero vector and not a guess from a previous solve. Note that some solves always start at zero in which case this switch will have no effect (to check this, the user should refer to the log output). 

\subsection{Remove Null Space}

As documented in Section \ref{Nullspaceremove}, this option removes the null space.

\subsection{Solver Failures}
\index{errors!linear solver}

Three options are available here:

\begin{enumerate}
\item Never ignore solver failures: Solver failures are always treated as fatal errors. The model stops at the end of the time step in order to allow for the latest output to be written. 
\item Ignore non-convergence during spin-up: Allow for an initial period in which solver failures caused by non-convergence in the maximum number of iterations are ignored. 
\item Ignore all solver failures: Ignore all solver failures. This is a dangerous option that should only be used in exceptional cases. 
\end{enumerate}

It is recommended that users use the first option: Never ignore solver failures, however, on occasions (e.g. challenging initial conditions) the second might also be applicable. 

\subsection{Reordering RCM}

A bandwidth reduction algorithm --- reverse Cuthill-McKee reordering --- is used to improve cache performance.

\subsection{Solver Diagnostics}

This subsection includes a series of extra diagnostic options to help debug solver problems. 

\subsubsection{Print norms}
Print out the norm of vectors and matrices before the solve, and that of the solution vector afterwards. Norms are printed at verbosity level 2, so run \fluidity\ with -v2 or -v3.

\subsubsection{Monitors}
Options to give extra information for each iteration of the the solve. Note that some of those may really slow down your computation. 


\section{Equation of State (EoS)}\label{sec:ConfigEOS}
\index{equation of state!options}

The equation of state is a relation between state 
variables. For incompressible flows it is used to derive the density
from other variables such as temperature and salinity (cf. section \ref{sec:IncompressibleFlow}). For compressible
flows it can be a more general relation between the state variables
including density and pressure.

The following EOS are available:

\begin{description}
\item\option{\ldots/equation\_of\_state/fluids/linear} 
Is a simple linear equation of state,
where density is a function of temperature and salinity.

\item\option{\ldots/equation\_of\_state/fluids/ocean\_pade\_approximation} Is a complex EOS for ocean modelling where density is a function of temperature, salinity and pressure.
\item\option{\ldots/equation\_of\_state/compressible/miegrunneisen} Is a simple compressible material EOS.
\end{description}

\subsubsection{Linear fluid EOS}
\index{equation of state!linear}
The density is a linear function of temperature, salinity and any number of generic scalar fields:
\begin{equation}
  \rho=\rho_0 \left(1-\alpha(T-T_0)+\beta (S-S_0) - \gamma (F - F_0) \right),
\end{equation}
where $\rho_0, \alpha, T_0, \beta, S_0, \gamma$ and $F_0$ are set by the following 
options:
\begin{description}
\item \option{\ldots/linear/reference\_density} sets $\rho_0$
\item \option{\ldots/linear/temperature\_dependency/thermal\_expansion\_coefficient} sets $\alpha$
\item \option{\ldots/linear/temperature\_dependency/reference\_temperature} sets $T_0$
\item \option{\ldots/linear/salinity\_dependency/salinity\_contraction\_coefficient} sets $\beta$
\item \option{\ldots/linear/salinity\_dependency/reference\_salinity} sets $S_0$
\item \option{\ldots/linear/generic\_scalar\_field\_dependency/expansion\_coefficient} sets $\gamma$
\item \option{\ldots/linear/generic\_scalar\_field\_dependency/reference\_value} sets $T_0$
\end{description}
Note that for Boussinesq the reference density does not 
influence any of the terms in the 
momentum equation (see \eqref{boussinesq}). It may influence the outcome
of diagnostic fields depending on density.

The option \option{subtract\_out\_hydrostatic\_level} only changes 
the buoyancy term. For LinearMomentum it changes to 
$g(\rho-\rho_0)$ and does not affect the density in the $Du/Dt$ term. For
Boussinesq it changes to 
$g\rho'=g(\rho-\rho_0)/\rho_0$ (again see \eqref{boussinesq}), 
and this option should always be used.
In both cases the diagnostic
``Density'' field and all other diagnostic fields depending on density
still represent the full density.

\subsubsection{Pade ocean EOS}
\index{equation of state!Pade approximation}
This EOS is described in section~\ref{sec:PadeDescription}.
This option uses mean hydrostatic pressure based on depth to calculate the
pressure (hence why you need to provide the value of z on the top surface).
For this option, the temp field represents potential temperature \emph{not} in
situ temperature - so beware (see \citet{mcdougall2003} for a formula for converting
from in-situ to potential). The units are degrees Centigrade for potential
temperature, \PSU{} for salt, \kgmm{} for density. The reference density is
\kgmm[1000] and the momentum equation is Boussinesq using this reference density.

\subsubsection{Compressible EOS}\label{sec:Multi-material compressible EOS}
\index{equation of state!stiffened gas}
\option{\ldots/equation\_of\_state/compressible/miegrunneisen} defines a simple compressible equation of state that can be used to describe gases, liquids and solids, known as the stiffened gas EOS. %(see section~\ref{sec:StiffenedGas}).

\begin{description}
\item\option{\ldots/miegrunneisen/reference\_density} Specifies the reference density of the material in SI units.
\item\option{\ldots/miegrunneisen/ratio\_specific\_heats} Specifies the ratio of
  specific heats of the gas minus 1 $(c_p/c_V-1)$ in the perfect gas EOS and the
  Gr\"uneisen parameter in the stiffened gas equation of state. Not activating this option simplifies the compressible EOS to that of a compressible liquid.
\item\option{\ldots/miegrunneisen/bulk\_sound\_speed\_squared} Specifies the bulk sound speed squared for the material $c_B^2$. Not activating this option simplifies the compressible EOS to that of a perfect gas.
\end{description}

\section{Sub-grid Scale Parameterisations}
\label{sec:sub-grid-scale-parameterisations}
\index{turbulence models}

\fluidity\ contains a number of sub-grid scale parameterisations which model physical process below the resolution of the mesh.

\subsection{GLS}

This option enables the model described in section \ref{sec:GLS}. There are a few different 
sub-options to configure. First, you must choose which GLS method to use 
from $k-\epsilon$, $k-kl$, $k-\omega$ and $gen$. Next, the stability functions
can be chosen. CanutoA or CanutoB are recommended. If you are running a 3D model, then switching on  
\option{\ldots/calculate\_boundaries} is recommended in order for the boundary conditions to be set correctly. 
Finally, you can enable a number of optional diagnostic fields.

The user can also choose to relax the diffusivity and viscosity calculated by switching on the\linebreak
\option{\ldots/relax\_diffusivity}. The value specified must be between 0 and 1.0. A value of 0 indicates
no relaxation, 1.0 would indicate no changes to be made. If this option is activated, the diagnostic fields
\option{GLSTurbulentVerticalDiffusivity} and \option{GLSTurbulentVerticalViscosity} must also
be activated. In addition, if adaptivity is enabled, these two fields must have an interpolation
method set, e.g. \option{\ldots/GLSVerticalViscosity/diagnostic/consistent\_interpolation}.

For each field that will be effected by the subgrid scale parameterisations, 
you must enable the correct diffusivity. This
is done by specifying \option{\ldots/subgridscale\_parameterisation} in the
field to \option{GLS}. Normally, this would be the temperature, salinity and any biology fields active.

Finally, fields that are altered by the GLS model, such as the Viscosity, need to be switched
to a \option{diagnostic/algorithm::Internal}. The list of fields to switch is:
\begin{enumerate}
\item \option{GLSTurbulentKineticEnergy/Diffusivity}
\item \option{GLSTurbulentKineticEnergy/Source}
\item \option{GLSTurbulentKineticEnergy/Absorption}
\item \option{GLSGenericSecondQuantity/Diffusivity}
\item \option{GLSGenericSecondQuantity/Source}
\item \option{GLSGenericSecondQuantity/Absorption}
\item \option{Velocity/Viscosity}
\end{enumerate}

If these fields are not set correctly, a user error will occur.

\subsection{k-$\epsilon$ Turbulence Model}\label{sec:kepsilon_usage}
\index{turbulence model}
\index{k-$\epsilon$ model}

\option{../subgridscale\_parameterisations/k-epsilon} enables the turbulence model
described in \ref{sec:kepsilon}. It must not be confused with the $k-\epsilon$ option in
the GLS model (see \ref{sec:GLS}) which is only for oceans like problems.

In \fluidity\ the $k$ field is called \option{TurbulentKineticEnergy} and the $\epsilon$
field is called \option{TurbulentDissipation}.  The model works well with the following
spatial discretisation options for $k$ and $\epsilon$:
\begin{itemize}
\item
  \option{\ldots/control\_volumes/face\_value::FiniteElement/limit\_face\_value/limiter::Sweby}
\item \option{/control\_volumes/diffusion\_scheme::ElementGradient}.
\item significant speeded up is achieved by selecting
  \begin{itemize}
  \item \option{/project\_upwind\_value\_from\_point/}
  \item \option{store\_upwind\_elements\_store\_upwind\_quadrature}.
  \end{itemize}
\end{itemize}

Under temporal discretisation option:
\begin{itemize}
\item Fully implicit or Crank-Nicholson temporal discretisation is recommended.
\item \option{control\_volume} discretisation option should be selected
\item \option{number\_advection\_iterations} = 3.
\end{itemize}

The following fields are altered by the k-epsilon model and need to be switched on and set
to \option{diagnostic/algorithm::Internal}:
\begin{enumerate}
\item \option{Velocity/Viscosity}
\item \option{TurbulentKineticEnergy/Diffusivity}
\item \option{TurbulentKineticEnergy/Source}
\item \option{TurbulentKineticEnergy/Absorption}
\item \option{TurbulentDissipation/Diffusivity}
\item \option{TurbulentDissipation/Source}
\item \option{TurbulentDissipation/Absorption}
\end{enumerate}

Although it is possible to set some of these fields to \option{prescribed} or off
altogether for the purposes of debugging results and/or the code itself. If these fields
are not set correctly, a user error will occur and/or warnings will be displayed.

Boundary conditions for the $k$ and $\epsilon$ fields should be set to type
\option{k\_epsilon}. Read section \ref{sec:kepsilon} for information about how to choose
the correct settings for boundary conditions.

The molecular (or laminar) viscosity must be prescribed under
\option{\ldots/k-epsilon/tensor\_field::BackgroundViscosity}. This must be set to
\option{prescribed/value::WholeMesh/anisotropic\_symmetric/constant} with all values set
to the isotropic viscosity.

If using additional scalar fields such as temperature, salinity etc, an option is
available under \option{\ldots/subgridscale\_parameterisation::k-epsilon} to use the eddy
diffusivity scaled by a user-specified Prandtl number. A background diffusivity is
required and this can either be set under
\option{\ldots/k-epsilon/tensor\_field::BackgroundDiffusivity} or within the additional
scalar field under
\option{\ldots/subgridscale\_parameterisation::k-epsilon/background\_diffusivity}, with
the latter taking priority. Buoyancy effects can also be accounted for by setting
\option{\ldots/subgridscale\_parameterisation::k-epsilon/buoyancy\_effects} within scalar
fields.

There are options available to allow the treatment of each of the source terms in the
$k-\epsilon$ model as either source or absorption terms within Fluidity. Choosing implicit
calculation causes the terms to be calculated semi-implicitly as absorption
terms. Choosing explicit calculation calculates the term completely explicitly as a source
term.

In calculation of the source terms it is required to invert a mass matrix. For P1
discretisations this can be done using mass lumping. For other discretisations it is
necessary to solve a system of equations. When using non P1 discretisations for the $k$
and $\epsilon$ fields the option
\option{\ldots/k-epsilon/mass\_lumping\_in\_diagnostics/solve\_using\_mass\_matrix} should
be selected, along with suitable solver settings.

Within \option{\ldots/k-epsilon/debugging\_options} there are a number of options for
debugging the model. It is possible to output each term in the $k$ and $\epsilon$
equations separately, provide prescribed sources for the $k$ and $\epsilon$ equations, and
also to disable specific terms in the equations.

When using the Low Reynolds number model, which is enabled whenever a
\option{k\_epsilon/lowRe} boundary condition is specified, a \option{DistanceToWall} field
must be specified. For simple geometries the simplest method of providing this information
is to use a python function. For complex geometries where this is not possible precursor
Eikonal equation or Poisson equation simulations must be run to determine the values for
this field. These can be carried out quite easily using Fluidity. Details of this process
can be found in Tucker, P 2011: "Hybrid Hamilton/Jacobi/Poisson wall distance function
model" and Elias et al 2007: "Simple finite element-based computation of distance
functions in unstructured grids"

\subsection{Large Eddy Simulation Models}

LES models are available as options under \option{\ldots/Velocity/prognostic/spatial\_discretisation/continuous\_galerkin/les\_model}. See \ref{sec:LES} for details of the various LES models available. These models require a prescribed viscosity (\ldots/Velocity/prognostic/tensor\_field::Viscosity/prescribed), to which an eddy viscosity is added to account for subgrid-scale turbulence. The LES models are currently restricted to use in incompressible flow cases, where the discrete velocity is divergence-free and the eddy viscosity tensor is traceless.

\subsubsection{Second-Order Smagorinsky}

The modified second-order Smagorinsky model of \citet{bentham2003} is available under \option{\ldots/les\_model/second\_order}. The Smagorinsky coefficient (\option{\ldots/second\_order/smagorinsky\_coefficient}) must be set; its value should be that suggested by the literature for a particular flow type. A reasonable all-round figure is 0.1. The eddy viscosity is available as an optional diagnostic field (\option{\ldots/second\_order/tensor\_field::EddyViscosity}).

\subsubsection{Fourth-Order Smagorinsky}

The fourth-order Smagorinsky model of \citet{bentham2003} is available under \option{\ldots/les\_model/fourth\_order}. The Smagorinsky coefficient (\option{\ldots/second\_order/smagorinsky\_coefficient}) must be set; 0.1 is recommended. A fine mesh is required to get good results from this model.

\subsubsection{WALE}

The wall-adapted local eddy viscosity (WALE) model is available under \option{\ldots/les\_model/wale}. The Smagorinsky coefficient (\option{\ldots/second\_order/smagorinsky\_coefficient}) must be set; 0.1 is recommended.

\subsubsection{Dynamic LES}

The Germano dynamic LES model is available under \option{\ldots/les\_model/dynamic\_les}. The following options have to be set: first, the filter width ratio $\alpha$ (\option{\ldots/dynamic\_les/alpha}); 2 is recommended. Second, the solver options (\option{\ldots/dynamic\_les/solver}) are for solving the inverse Helmholtz equation for the test-filtered velocity; cg/SOR is recommended.

Optional options:
\begin{itemize}
\item \option{\ldots/dynamic\_les/enable\_lilly}: use the Lilly modification to the Germano model. It is recommended.
\item \option{\ldots/dynamic\_les/enable\_backscatter}: allows negative eddy viscosity, which may result in more realistic turbulent flow if the mesh resolution is fine enough.
\end{itemize}

Several diagnostic fields are available if desired:

\begin{itemize}
\item \option{\ldots/dynamic\_les/vector\_field::FilteredVelocity}: the velocity field filtered with the test filter.
\item \option{\ldots/dynamic\_les/tensor\_field::FilterWidth}: the mesh size tensor
\item \option{\ldots/dynamic\_les/tensor\_field::StrainRate}: the strain rate $\overline S_{ij}$.
\item \option{\ldots/dynamic\_les/tensor\_field::FilteredStrainRate}: the filtered strain rate $\widetilde{\overline S}_{ij}$.
\item \option{\ldots/dynamic\_les/tensor\_field::EddyViscosity}: the eddy viscosity $\nu_T$.
\end{itemize}

\section{Boundary conditions}\label{Sec:BCs_configure}

\index{boundary conditions}

The simulated system requires suitable boundary conditions for full closure.
An example could be the amount of sunlight at the ocean surface, a specified value of
temperature heating material from below, or a momentum stress in the form of wind for velocity.
It is also possible to leave boundary conditions undefined, in which case "stress-free" conditions are
applied. See section \ref{sec:BCs} for further details.

\subsection{Adding a boundary condition}\label{sec:BCs:adding}

Boundary conditions are set for each field contained in state under \option{\ldots/boundary\_conditions}. 
Multiple boundary conditions can be set for each field, such that the sides, surface and bottom can 
have different conditions. A boundary conditions is added by clicking the "+" symbol
in the appropriate field

\subsection{Selecting surfaces}\label{sec:BCs:selecting}
\index{surface ID}
To each boundary condition a set of domain surfaces is assigned on which it is applied to. The surfaces are identified by a surface ID specified during the mesh generation procedure (see section \ref{sec:surface_ids}). For example if the top and bottom of your mesh is defined as surface
1, then simply add a 1 to \option{\ldots/boundary\_conditions/surface\_ids}. Multiple surfaces 
can be added, separated by a space.

\subsection{Boundary condition types}\label{sec:BCs:types}
\index{boundary conditions!setting}
\fluidity\ supports a wide range of boundary conditions which will be introduced in the next sections.

\subsubsection{Dirichlet}
\index{boundary conditions!Dirichlet}

A Dirichlet condition sets the value of the field ($c$) at each location over the surface $\partial\Omega$:
\begin{equation*}
c(\bmx) = f(\bmx) \quad \textrm{on}\; \partial\Omega.
\end{equation*}

Dirichlet boundary conditions can also be applied weakly by selecting the \option{\ldots/apply\_weakly}
option. Unlike the strong form of the Dirichlet conditions, weak Dirichlet
conditions do not force the solution on the boundary to be pointwise equal
to the boundary condition. 

\subsubsection{Neumann}
\index{boundary conditions!Neumann}

A Neumann boundary condition sets a flux term $q$ to the normal ($\vec n$) of the surface $\partial\Omega$:
\begin{equation*}
\int_{\partial\Omega} \phi (\kaptens\nabla c)\cdot\bmn \;d\Gamma,
\end{equation*}
where $\phi$ is a test function (see section \ref{chap:numerical_discretisation}).
The Neumann condition is specified by assigning a value to the $q$, where
\begin{equation*}
q = (\kaptens\nabla c)\cdot\bmn,\quad \textrm{on}\quad \partial\Omega.
\end{equation*}

\subsubsection{Robin}
\index{boundary conditions!Robin}

A Robin boundary condition sets the diffusive flux term  $\vec{n}\cdot\tensor{\kappa}\cdot\grad c$ in weak form as:
\begin{equation*}
- \int_{\partial\Omega} \phi~\vec{n}\cdot\tensor{\kappa}\cdot\grad c =
    \int_{\partial\Omega} \phi~h(c-c_{a}),
\end{equation*}
where $\phi$ is a test function (see section \ref{chap:numerical_discretisation}).
The Robin condition is specified by assigning a value to the order zero $C_{0}$ and order one $C_{1}$ coefficient fields, where
\begin{align}
   C_{0} &= hc_{a}, \\
   C_{1} &= h.
\end{align}
Currently, the Robin boundary condition is only available for Continuous Galerkin and Control Volume 
spatial discretisations of a scalar field. 

\subsubsection{Bulk formulae}\label{sec:bulk_formulae}

These boundary conditions can be used on:
\begin{itemize}
\item Salinity
\item Temperature
\item Velocity
\item PhotosyntheticRadiation
\end{itemize}

They use meteorological data and convert it into a Neumann or Dirichlet boundary condition
as appropriate for the fields above. You do not need to have all the above fields; only 
velocity and temperature are required. More information can be found in section \ref{sec:BCs:special:oceans}.

\subsubsection{Zero flux}
\index{boundary conditions!zero flux}

For control volume discretisations only, this option prevents the field fluxing from the boundary.

\subsubsection{Flux}
\index{boundary conditions!flux}

For control volume discretisations only, this option allows a given flux $h$ of field $c$ through the boundary. In other words, we have
\begin{equation*}
   \frac{\partial c}{\partial t} = h
\end{equation*}

\subsubsection{Free surface}\label{subsec:free_surface_bc}
\index{boundary conditions!free surface}
\index{free surface!boundary condition}

The \option{\ldots/free\_surface} option allows the upper surface height to vary according to the pressure and velocity fields. This boundary condition is available on the velocity field only. When using a free surface, it is recommended that you active a diagnostic free surface (though this is optional). This option is also available at \option{\ldots/scalar\_field::FreeSurface}. 

Note that the default free surface treatment implements the physical boundary
condition p=0 (this is full pressure without subtracting the hydrostatic
component). For viscous fluids, the correct boundary condition is a no normal
stress condition, \eqref{NormalStressBC}, which includes a viscosity term. When using this option a prognostic free surface field is also required. This option is also available at \option{\ldots/scalar\_field::FreeSurface}. 

By default, the mesh geometry is not influenced by the free-surface calculation, however \fluidity\ can deform the mesh according to the free-surface elevation. 
This option is available at \option{/mesh\_adaptivity/mesh\_movement/free\_surface}.

\subsubsection{Wetting and drying}\label{subsec:wetting_drying_bc}
\index{boundary conditions!wetting and drying}

In order to use wetting and drying, first switch on the mesh deformation as described in \ref{subsec:free_surface_bc}.

Secondly, if the mesh is extruded within fluidity, the extrusion parameters have to be changed such that areas above sea level are included.
For example if a bathymetry map file is used for the extrusion, the option \option{/geometry/mesh/from\_mesh/extrude/regions/bottom\_depth/from\_map/surface\_height}
can be used to shift down the domain such that the whole bathymetry is below zero.
A non-zero initial pressure together with the relationship between pressure and free-surface elevation $p = \rho \eta$ can be used to shift the initial free-surface down accordingly as well.

Finally, wetting and drying is activated under \option{/mesh\_adaptivity/mesh\_movement/free\_surface/wetting\_and\_drying}.
The only required parameter is the wetting and drying threshold value $d_0$, which specifies the minimum layer-thickness that is retained in dry areas. 
Following equation can be used to determine the threshold value:
\begin{equation*}
\d_0 = \frac{l\Delta x}{r},
\end{equation*}
where $\Delta x$ and $l$ are the maximum horizontal element size and number of mesh layers in the dry areas, respectively and $r$ is the maximum aspect ratio. A typical value for latter is between $500-1000$.


\subsubsection{Drag}
\index{boundary conditions!drag}

This option applies a quadratic or linear drag to the Velocity field. Both the value and the type of drag need to be set. A Manning-Strickler drag can be used by activating \option{\ldots/quadratic\_drag/manning\_strickler}

\subsubsection{Wind forcing}\label{sec:wind_forcing}
\index{boundary conditions!wind stress}
\index{wind forcing}

A wind forcing can be applied to the Velocity field as either a stress or
velocity. For stress values, the physical units should match those of the
simulation, so for example, if you use the non-dimensional value of $\rho$
as 1.0, your stresses (in \unit{kgm\ensuremath{^{-1}}s\ensuremath{^{-2}}})
should be divided by the reference density.  If using wind velocity
(at 10m height) the density of the air needs to be specified in the same
units, i.e. $\rho_{\textrm{air}} = 1.3\times10^{-3}$.

Alternatively
\option{\ldots/Velocity/boundary\_conditions/wind\_forcing/wind\_stress}
sets the value of wind forcing from a NETCDF file. The NETCDF file must
contain East-West and North-South components, along with times locations
(latitude/longitude) for the values. In addition, one must set
\option{/timestepping/current\_time/time\_units} in order for the simulated
time to be matched to the NETCDF data.

\subsubsection{No normal flow}
\index{boundary conditions!no normal flow}

When using \option{\ldots/control\_volumes} under Pressure \option{\ldots/spatial\_discretisation} or when using 
\option{\ldots/integrate\_continuity\_by\_parts} with CG Pressure and Velocity this boundary condition type 
imposes a weak no normal flow boundary condition on the surfaces specified.

\subsubsection{Prescribed normal flow}
\index{boundary conditions!prescribed normal flow}

When using \option{\ldots/control\_volumes} under Pressure \option{\ldots/spatial\_discretisation}, when using 
\option{\ldots/integrate\_continuity\_by\_parts} with CG Pressure and Velocity, or DG Velocity this boundary condition type 
imposes a weak prescribed normal flow boundary condition on the surfaces specified.


\subsection{Internal boundary conditions}\label{sec:internal_bcs}
\index{boundary conditions!internal boundaries}
Fluidity has some support (with limitations) for applying boundary conditions
along internal boundaries. This is done in the normal way 
by applying physical line/surface ids to the internal boundary, and associating
boundary conditions with these in the options file. It should be
noted that this functionality is not as well tested as other boundary
conditions, so you should always first try it out on a simplified version of
your problem to see if it behaves as expected. Strong Dirichlet boundary
conditions should work correctly, other boundary condition types may not give
the right answer. Since the gmsh .msh format does not give you the
option of setting a different surface id on either side of the internal
boundary, it is not possible to apply a different boundary condition type or
value on either side.

When using mesh adaptivity, it will maintain the internal boundary and its
surface id. Mesh adaptivity in three dimensions with internal boundaries 
is currently not supported. For parallel runs you have to use flredecomp (see
section \ref{mesh!meshing tools!flredecomp}) to
decompose the problem, as fldecomp will refuse any mesh with internal boundaries.

\subsection{Special input date for boundary conditions}\label{sec:BCs:specialised}

When running free surface simulations the surface elevation at the boundary is specified by applying a pressure Dirichlet condition. Since the free surface elevation is often measured data, there are some special possibilities to specify a pressure Dirichlet condition: 

\option{\ldots/from\_file} allows the specification of a single file
containing something useful. This option is available on the Pressure (Free
Surface) field. Tidal boundary conditions can be applied by setting this option
and referencing a relevant NetCDF file containing appropriate amplitude and phase data for the desired
tidal constituent(s). The file is referenced under:
\begin{itemize}
\item \option{\ldots/tidal/file\_name},
\end{itemize}
with the amplitude and phase names (as specified in the NetCDF file) set under:
\begin{itemize}
\item \option{\ldots/tidal/variable\_name\_amplitude},
\item \option{\ldots/tidal/variable\_name\_phase}
\end{itemize}
respectively. Finally, the constituent should be selected from the list under:
\begin{itemize}
\item \option{\ldots/tidal/name}.
\end{itemize}
A separate tidal boundary condition needs to be set for each constituent.

\index{boundary conditions!NEMO data}
\option{\ldots/NEMO\_data} will set the field according to a specified NEMO
input file. This option is available for the Pressure (Free Surface) field.
In order to use this option a prescribed field containing NEMO pressure
field data must first be created. See section \ref{sec:setting_from_nemo}
for information on setting prescribed fields from NEMO data. Then, under
\option{\ldots/NEMO\_data/field\_name}, set the string to that of the prescribed
field containing the NEMO pressure data to enable this option.

\index{boundary conditions!synthetic eddy method}
\option{\ldots/synthetic\_eddy\_method} Available for velocity.
This generates statistically realistic turbulent flow at an inflow using a
statistical method (for a full explanation see \cite{jarrin_06}.
The user specifies a mean velocity (e.g. python profile),
turbulence lengthscale, Reynolds stress profile and number of samples.
This is useful for high-Reynolds-number industrial CFD flow, and/or if using an LES model.

\subsection{Special cases}\label{sec:BCs:special}

There are a few special cases of boundary conditions that are not applied
using the methods described above.  These include ocean surface forcing and
the boundary conditions on the General Length Scale (GLS) turbulence model.

\subsubsection{Ocean surface forcing}\label{sec:BCs:special:oceans}
\index{boundary conditions!ocean}

Ocean surface forcing takes parameters from ERA40 datasets, passes them
through bulk formulae and gives a boundary condition for the salinity,
temperature, photosynthetic radiation and velocity fields. The settings for
these options are in \option{/ocean\_forcing/bulk\_formulae}. However, you must also set up
\option{/timestepping/current\_time/time\_units}.

Under \option{/ocean\_forcing/bulk\_formulae} an input file must be defined. The fields on 
which bulk formulae are to be imposed should have their upper surface set to the correct
boundary condition type (\option{bulk\_formulae}).
The input file must contain the following ERA40 parameters for the
duration of the simulated time:
\begin{itemize}
 \item 10 metre U wind component (\ms)
 \item 10 metre V wind component (\ms)
 \item \m[2] temperature (\K)
 \item Surface solar radiation downward (\unit{Wm\ensuremath{^{-2}}s})
 \item Surface thermal radiation downward (\unit{Wm\ensuremath{^{-2}}s})
 \item Total precipitation (\unit{ms})
 \item Run off (\unit{ms})
 \item \m[2] dew point temperature (\K)
 \item Mean sea-level pressure (\Pa)
\end{itemize}

These variables are surface variables as defined by data files from the ERA40 website. Note that some parameters are accumulated values
and as such are required to be divided by the ERA40 temporal resolution - \fluidity { } assumes 6 hour temporal resolution. 
These parameters are used as input to the default bulk forcing formulae of \citet{large2004} included in \fluidity. Other
formulae are available: COARE 3.0 \citep{fairall2003} and those of \citet{kara2005} which are based on the COARE data.

Other options under ocean surface forcing include specifying a latitude and longitude, and using a single position for the
forcing data. These options are only really useful when simulating pseudo-1D columns (see the gls-StationPapa test for an example
of a pseudo-1D column). Enabling the \option{position} option allows the user to specify a latitude and longitude as
two real numbers (e.g. 50.0 -145.0 for 50$^\circ$ N and 145$^\circ$ W). These co-ordinates are translated into Cartesian
co-ordinates, which are then added to the positions of the surface of the mesh. This allows the use of simple mesh geometries
and co-ordinates, whilst still specifying where the forcing data should originate. Moreover, the \option{single\_location}
option forces \emph{all} surface nodes to receive the same forcing.

Finally, it is possible to output the fluxes that are imposed on the ocean surface, by enabling the
\option{output\_fluxes\_diagnostic} option. Here, the user can enable diagnostic fields for momentum, heat, salinity and
photosynthetic radiation downwards. The fluxes will then be included in the output as normal scalars or vectors, but with values
confined to the upper surface.

\subsubsection{GLS sub-grid scale parameterisation}\label{sec:BCs:special:gls}
\index{generic length scale model!boundary conditions}
\index{boundary conditions!generic length scale model}

The GLS model (see section \ref{sec:GLS}) requires that Neumann boundary
conditions are set for stability, however, the boundary conditions on the
Generic Second Quantity ($\Psi$) depend on other modelled variables. In
order for the boundary conditions to be set correctly, enable the \linebreak
\option{\ldots/subgridscale\_parameterisations/GLS/calculate\_boundaries}
option.

\subsubsection{k-epsilon sub-grid scale parameterisation} \label{sec:BCs:special:kepsilon}
\index{k-$\epsilon$ model!boundary conditions}
\index{boundary conditions!k-$\epsilon$ model}

The k-epsilon turbulence model (see section \ref{sec:kepsilon}) should apply zero
Dirichlet boundary conditions to the TurbulentKineticEnergy ($k$) field. The
TurbulentDissipation ($\epsilon$) field should use the special type of Dirichlet condition
called \option{k\_epsilon} which is calculated in the k-epsilon module. To enable
calculation of the boundary conditions on both fields, set the \linebreak
\option{\ldots/subgridscale\_parameterisations/k-epsilon/calculate\_boundaries} option.

\section{Astronomical tidal forcing}
\label{config:tides}
\index{tides}

Astronomical tidal forcing can be switched on for 11 different constituents
under:
\begin{itemize}
\item \option{/ocean\_forcing/tidal\_forcing}, 
\end{itemize}
(see \citealp{Wells2008} for descriptions of
the different constituents). These can be switched either
individually or in combination. In addition, a body tide correction can be stipulated
under:
\begin{itemize}
\item \option{\ldots/tidal\_forcing/love\_number},
\end{itemize}
for which the suggested value is 0.3 (assuming 
Love numbers of $k$=0.3 and $h$=0.61; see section \ref{astronomical}).   

Note that for many cases, specifically those involving open boundaries, it is often
desirable to combine astronomical tidal forcing with a co-oscillating boundary tide condition
(see section \ref{sec:BCs:specialised}).

\section{Ocean biology}
\index{biology}

Enabling this turns on the ocean biology model. In addition you also need to add several scalar fields in the first material phase:
\begin{itemize}
\item Phytoplankton
\item Zooplankton
\item Nutrient
\item Detritus
\item Primary production
\end{itemize}

There are several items that need configuring before biology can be used.
First a relationship between sources and sinks needs encoding. This is best
done by importing fluidity.ocean\_biology into
\option{/ocean\_biology/pznd/source\_and\_sink\_algorithm} and calling the
models from there. An example is given below.

\begin{example}
  \begin{lstlisting}[language=Python]
import fluidity.ocean_biology as biology

day=1./(3600*24)

p={}
p["alpha"]=0.015*day
p["beta"]=0.75
p["gamma"]=0.5
p["g"]=1*day
p["k_N"]=0.5
p["k"]=0.5
p["mu_P"]=0.1*day
p["mu_Z"]=0.2*day
p["mu_D"]=0.05*day
p["p_P"]=0.75
p["v"]=1.5*day

biology.pznd(state, p)  
  \end{lstlisting}
  \caption{A Python function that imports the biology module and sets the algorithm to use.}
\end{example}

The final thing to change is to add absorption coefficients in the photosynthetic radiation field for water and plankton concentration.

\section{Sediment model}
\label{config:sediments}
\index{sediments}

\fluidity\ contains a sediment model in which sediment is treated as a
tracer with a settling velocity. 
It is possible to specify multiple sediment fields to represent a distribution of sediment
characteristics.
Sediment that falls out of the domain due to settling can be recorded using a
Bedload field

Note: To use sediment, a linear equation of state must also be enabled
\option{\ldots/equation\_of\_state/fluids/linear}

\section{Population balance model}
\label{config:pbe}

The population balance model can be enabled by activating the \option{population\_balance (DQMOM)} option under the material phase.
The number of \option{abscissa}, \option{weights} and \option{weighted-abscissa} scalar fields must be chosen to be equal to the desired number of mesh nodes $N$ in the DQMOM quadrature approximation (Eq \eqref{eq:DQMOM_ndf}). 

\option{weights} and \option{weighted-abscissa} are prognostic fields and the \option{equation(AdvectionDiffusion)} option much be selected for all of these fields. The source terms in these prognostic scalars must be set to \option{diagnostic} with \option{algorithm(Internal)}. \option{abscissa} fields on the other hand are diagnostic fields. All scalars in the population balance tree must be defined on the same mesh. When the population balance model is used in conjunction with the multiphase flow equations for modelling polydispersed flows, it is suggested that the control volume discretisation should be used for the population balance scalars to ensure conservation and to ensure a realisable moment set (see the section on \emph{moment corruption} in \citet{bhutani2016}). 

Spatial diffusion coefficient must be identically specified for all prognostic scalars in the population balance model, if the spatial diffusion is present. 

\subsection{Source terms}
The following source terms are available for the PBE to be used in \fluidity\ :
\subsubsection{Growth}
\begin{itemize}
\item power law growth
\end{itemize}
\subsubsection{Internal Dispersion}
\begin{itemize}
\item constant internal dispersion coefficient
\end{itemize}
\subsubsection{Aggregation}
\begin{itemize}
\item constant aggregation
\item hydrodynamic aggregation
\item sum aggregation
\item turbulent aggregation -- \citet{laakkonen2007modelling}
\end{itemize}
\subsubsection{Breakage}
Breakage frequency:
\begin{itemize}
\item constant
\item power law breakage
\item turbulent breakage -- \citet{laakkonen2007modelling}
\end{itemize}

Breakage daughter distribution function:
\begin{itemize}
\item symmetric fragmentation
\item \citet{mccoymadras2003analytical}
\item \citet{laakkonen2007modelling}
\end{itemize}

\subsection{Initial conditions}
It is possible to specify initial condition for the prognostic scalars, i.e. weights and the weighted-abscissas (\option{use\_prognostic\_field\_initial\_conditions}), or the moments (\option{calculate\_initial\_conditions\_from\_moments}). 

\subsection{Moments and statistics}
Diagnostic scalar fields such as moments and statistics of the NDF can be calculated by enabling these options under \option{moments} and \option{statistics}, respectively. The moments must be named as \emph{Moment\_X} where \emph{X} is the moment number. The following statistics are available currently: mean, standard deviation, skew, Sauter mean diameter ($m_3$/$m_2$), and mean diameter ($m_1$/$m_0$). 

Further details of the population balance model can be found in \citet{bhutani2016polydispersed} and \citet{bhutani2016}.

\section{Configuring ocean simulations}
\label{sec:large-scale-ocean}

This section contains advice for running a large scale ocean simulation in
three dimensions. The flow domain is characterised by the large aspect ratio
between horizontal and vertical length scales. The following is a set of
recommended discretisation options for such simulations. It is based on the use
of the \PoDGPt velocity, pressure element pair \citep{cotter2009}. The
recommendations are applicable to both global or ocean-scale domains as well as
to smaller, regional or coastal domains. For the first case, see the additional
instructions in section \ref{sec:ocean_on_the_sphere}, to solve the equations
on the sphere.

\subsection{Meshes}

The mesh that you use should be unstructured in the horizontal and structured
in the vertical (this is to ensure that hydrostatic and geostrophic balances can
be maintained). It can either be constructed in gmsh and read into fluidity, or
a two-dimensional mesh can be made in gmsh and extruded within fluidity (see
section \ref{sec:extruded_meshes}). To enable the \PoDGPt element pair you will need
the two following additional meshes that are derived from the three-dimensional
CoordinateMesh:
\begin{itemize}
\item The VelocityMesh must be made discontinuous galerkin, by setting
  \option{/geometry/mesh::VelocityMesh/from\_mesh/mesh\_shape/mesh\_continuity} to discontinuous.  
\item The PressureMesh must be made quadratic, by setting
\\* \option{/geometry/mesh::PressureMesh/from\_mesh/mesh\_shape/polynomial\_degree}  to $2$.  
\end{itemize}
Also under \option{/geometry/}, the \option{ocean\_boundaries} option
should be switched on, with the surface ids of the top and bottom specified (if
you are using extrusion these should correspond to the surface ids specified
there). This is necessary to make the diagnostic FreeSurface field, and the
\option{vertical\_lumping} option for the \option{mg} preconditioner work. In
parallel, this requires a 2D input mesh and extrusion within fluidity.

\subsection{Time stepping}
\begin{itemize}
\item The timestep is typically chosen based on some knowledge of the time scales
  involved in the simulation. An adaptive timestep is less applicable here as
  this only looks at advective time scales. Other time scales that may be
  applicable: a Courant condition based on the barotropic wave speed, $c=\sqrt{gH}$ where
  $H$ is the water depth, the frequencies of the forcings, baroclinic timescales
  (e.g. the Brunt-V\"ais\"al\"a frequency), etc. Since all equations are solved in
  an implicit manner, a relatively large timestep can be chosen. This will
  however impact the temporal accuracy, so the timestep should be chosen with
  the timescale of the physical phenomena of interest in mind. Despite the
  implicit solution technique, choosing too large a timestep may still hamper
  the stability of the model.
\item \option{/timestepping/nonlinear\_iterations} is typically set to 2, to
  deal with the nonlinearities of the advection term (and coupling of the
  temperature and salinity equations and the buoyancy term in baroclinic
  simulations).
\item To ensure that the equations are indeed stable for large timesteps choose
  a value of $\theta>=0.5$ under \option{\ldots/temporal\_discretisation}
  under Velocity (and Temperature and Salinity if applicable). It is
  recommended to start with a value of 1.0 (most stable choice), and only
  reduce this after it is established that the model is sufficiently stable, in
  order to improve the temporal accuracy of the model (reduces the artificial 
  dissipation/wave damping).
\item As explained above, the timestep that is chosen can often be much larger
  than would be allowed by the Courant condition. In addition, with DG
  discretisations the maximum allowed Courant condition is much smaller than
  1.0. Therefore, we recommend switching on sub-cycling by choosing
  \option{discontinuous\_galerkin} under
  \option{\ldots/temporal\_discretisation}. A value of 0.1 is safe for the
  \option{maximum\_courant\_number\_per\_sub\_cycle}, a value of 0.5 can be a
  bit more efficient.
\item To enable the last option, you also need to add a
  diagnostic \option{scalar\_field::DG\_CourantNumber} under the
  \option{/material\_phase}.
\end{itemize}

\subsection{Velocity options}
The velocity discretisation is discontinuous galerkin. Recommended options under
\option{\ldots/vector\_field::Velocity/prognostic} are:

\begin{itemize}
  \item The \option{mesh} should of course be the VelocityMesh.
  \item \option{\ldots/equation} should be set to Boussinesq.
  \item \option{\ldots/spatial\_discretisation} must be discontinuous galerkin
  \item \option{\ldots/spatial\_discretisation/discontinuous\_galerkin/advection} is upwind
  \item
    \option{\ldots/advection\_scheme/project\_velocity\_to\_continuous} is
    recommended. The mesh that is chosen should be CoordinateMesh, unless
    using \option{super\_parametric\_mapping} on the sphere (see below).
  \item
    \option{\ldots/spatial\_discretisation/advection\_scheme/integrate\_advection\_by\_parts}
    is twice.
  \item For \option{\ldots/solver/} options, \option{gmres} for the
    \option{iterative\_method} and \option{sor} for the preconditioner are
    usually sufficient.
\end{itemize}

\subsection{Choosing Viscosity values}
The choice of viscosity in a realistic, large-scale ocean domain is a difficult
topic. The molecular (kinematic) viscosity of water, 1e-6 m$^2/$s, is too small
to have any influence at these length scales. Rather, viscosity should be thought
of as a parameterisation of the energy dissipation that happens at length scales
smaller than those that are resolved (the sub-grid scale). For smaller 
scale, CFD-type problems this viscosity can be calculated by standard turbulence
closure models (see section \ref{sec:sub-grid-scale-parameterisations}). For
larger scale problems the situation is less clear. Often, a fixed ``eddy
viscosity'' is chosen as a tunable parameter. For the vertical, one of the many
GLS configurations can be used (section
\ref{sec:sub-grid-scale-parameterisations}). At the moment, we do not have the
capability for a separate, LES-type model for the horizontal. 

In some sense, viscosity can also be thought of as a stabilisation method.
Although many numerical schemes are designed to remain stable, putting in just as
much numerical diffusion as necessary, this usually only holds for ideal
circumstances. In practice, due to under-resolving, poor coastal and bathymetry
data, and unnatural (not well-posed) boundary conditions that have to be
applied, additional viscosity is needed to keep the model stable. Thus a good
way of choosing a viscosity value is to start with a value that is not too low
(this depends on mesh resolution) to obtain a stable model. Then, if the results
are deemed too smooth with excessive mixing, the viscosity value can be lowered.

\subsection{Pressure options}
This uses the Continuous Galerkin discretisation, with a mesh of polynomial order two (already specified above).

\begin{itemize}
\item The \option{mesh} be the PressureMesh and \option{\ldots/spatial\_discretisation} must be \option{continuous\_galerkin}.
\item Under that last option, select both \option{remove\_stabilisation\_term}
  and \option{integrate\_continuity\_by\_parts}.
\item Under \option{\ldots/scheme}, \option{poisson\_pressure\_solution} should
  be set to \option{only first timestep}, and \option{use\_projection\_method}
  should be selected.
\item Under \option{\dots/solver}, select \option{cg} for the
  \option{iterative\_method}, and \option{mg} for the \option{preconditioner}.
  Under \option{\ldots/preconditioner::mg}, switch on
  \option{vertical\_lumping}.
\end{itemize}

Having a boundary condition \option{type::free\_surface} requires the option
\option{equation\_of\_state/fluids/linear} under \option{/material\_phase}, 
and the option \option{subtract\_out\_hydrostatic\_level} under it. The \option{reference\_density}
should not have an effect on a \option{equation::Boussinesq} configuration, but
it is safest to simply set it to 1.0. Also without a free surface boundary
condition, it is a good idea to set these \option{equation\_of\_state} options.

\subsection{Boundary conditions}

\begin{itemize}
  \item For the bottom boundary condition, add both a \option{type::drag} and a
    \option{type::no\_normal\_flow} boundary condition under Velocity. A
    typical value for a dimensionless bottom friction coefficient for the ocean
    is $0.0025$.
  \item For the free surface, add a \option{type::free\_surface} boundary
    condition under Velocity. To be able to read the free surface elevations as
    a field, add a diagnostic \option{scalar\_field::FreeSurface} under
    \option{/material\_phase}.
  \item For closed boundaries (coast), add a \option{type::no\_normal\_flow}
    boundary condition, and, optionally, a \option{type::drag} boundary
    condition.
  \item For open boundary conditions, you should either prescribe the velocity
    or the pressure by adding a \option{type::dirichlet} boundary condition
    under \option{vector\_field::Velocity} or \option{scalar\_field::Pressure}
    respectively. Note, that in simulations with constant density where the
    hydrostatic pressure has been subtracted, the main component of pressure is
    the barotropic pressure $p=g\eta$ (in Boussinesq we have also divided out
    the reference density $\rho_0$). Thus, neglecting non-hydrostatic effects, a
    free surface boundary condition can be applied by multiplying the desired
    free surface elevation by $g$ and setting this as a \option{type::dirichlet}
    Pressure boundary condition. For baroclinic simulations, the vertical
    pressure profile should be taken into account and added to the barotropic
    pressure.
\end{itemize}

\subsection{Baroclinic simulations}
For baroclinic simulations we need to specify how the Temperature and Salinity
fields influence the Density. This is done under
\option{/material\_phase/equation\_of\_state/fluids/linear}. 
Although not strictly required, it is a good idea
to add a diagnostic \option{scalar\_field::Density} under
\option{/material\_phase} to be able to visualise the Density field (for
Boussinesq simulations this will have been divided by the reference density).

Obviously we should also add a Temperature and/or Salinity field under
\option{/material\_phase}. For the solution of their advection-diffusion
equation we have the choice of a Continuous Galerkin, a Discontinuous Galerkin
or a Control Volume discretisation. The following are recommended options for a
Continuous Galerkin discretisation of these fields:
\begin{itemize}
  \item Choose \option{continuous\_galerkin} under
    \option{\ldots/spatial\_discretisation}. 
  \item If the fields contain discontinuities (shocks or fronts), you can
    choose \option{streamline\_upwind} under
    \option{\ldots/spatial\_discretisation/stabilisation} to stabilise the
    solution.
  \item To ensure conservation of the scalar, choose a value of 1.0 for
    \option{\dots/spatial\_discretisation/conservative\_advection}. 
  \item For \option{\ldots/solver/} options, choose \option{gmres} for the
    \option{iterative\_method} and \option{sor} as \option{preconditioner}.
\end{itemize}

As mentioned in the time stepping section, our ability to use large timesteps
relies on the fact that the equations are solved implicitly. Because we do not
solve the Temperature and Navier Stokes equations in a directly coupled way
however, the timescale associated with this coupling, the buoyancy or
Brunt-V\"ais\"al\"a frequency still poses a, very strict, limitation on the timestep. To
avoid this, add the option
\option{\ldots/vertical\_stabilization/implicit\_buoyancy} under the Velocity
field.

\subsection{Spherical domains}
\label{sec:ocean_on_the_sphere}
For ocean scale simulations, in which the curvature of the Earth is important,
the equations should be solved in a spherical domain. Follow the instructions
in section \ref{sec:spherical_earth}. Note, that when using
\option{super\_parametric\_mapping} the CoordinateMesh should be P2. This is
recommended for large-scale baroclinic simulations, in particular when
under-resolved. It is
achieved by adding another mesh, typically called \option{BaseMesh}, that has the
\option{extrude} option under it instead of the CoordinateMesh. The
CoordinateMesh is now derived from the BaseMesh with the additional option of
\option{polynomial\_order} set to 2 (in the same way as the PressureMesh is
derived). The VelocityMesh and PressureMesh should now also be derived from the
BaseMesh. When using the option \option{project\_to\_continuous} under
\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/advection\_scheme},
choose BaseMesh as the mesh to project to.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% GFD problems %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Geophysical fluid dynamics problems}
\label{sec:gfd_problems}

This section contains advice for running Geophysical Fluid Dynamics (GFD) problems, such as laboratory-scale flows e.g. the lock-exchange and the annulus or smaller-scale ocean problems e.g. a gravity current on an incline. 
  Both continuous-Galerkin (\Poo) and discontinuous-Galerkin (\PoDGPt)
discretisations may be used. The \PoDGPt discretisation is better at
maintaining geostrophic balance, and more accurate. A run with the \Poo
discretisation on the same mesh will be faster. It does however require
stabilisation, which may lead to excessive wave energy dissipation. It can be a
good choice in a highly-energy turbulent regime.

For the \PoDGPt discretisation most of the recommendations in section
\ref{sec:large-scale-ocean} should be followed. The option
\option{vertical\_lumping} under \option{\ldots/preconditioner::mg} is no longer
necessary.

For the \Poo discretisation the following differences apply:
\begin{itemize}
  \item A separate VelocityMesh and PressureMesh are not necessary. Instead you
    can put all fields on the same CoordinateMesh or, alternatively, remove the
    \option{mesh\_continuity} and \option{polynomial\_degree} options under
    these meshes.
  \item The \Poo discretisation does not provide the option of subcycling
    (\option{\ldots/temporal\_discretisation/discontinuous\_galerkin}). Instead
    you could consider using an \option{adaptive\_timestep} under
    \option{/timestepping}.
  \item Under Velocity select
    \option{\ldots/spatial\_discretisation/continuous\_galerkin}. Under it,
    you need to have \option{mass\_terms/lump\_mass\_matrix}. You may need
    to add \option{stabilisation/streamline\_upwind} to stabilise non-smooth
    flows.
  \item The options for Pressure are the same as before, except the option
    \option{vertical\_lumping} is no longer necessary under
    \option{\ldots/preconditioner::mg}.
  \item Because the P1 pressure discretisation is not very good for maintaining
    geostrophic balance, it is usually necessary to solve this balance in a
    separate solve (see section \ref{sec:balance_pressure}). The part of pressure that balances the geostrophic
    Coriolis and buoyancy terms, is solved for in a separate field called
    GeostrophicPressure that is on a higher order mesh. Add a new mesh under
    \option{/geometry} with the name \option{GeostrophicPressureMesh}, with
    option \option{from\_mesh/mesh::CoordinateMesh} and add
    the option \option{polynomial\_order} with value 2.
    Add a prognostic \option{scalar\_field::GeostrophicPressure} under
    \option{/material\_phase} using that same mesh. Under
    \option{\ldots/spatial\_discretisation}, choose \option{include\_buoyancy}
    as the \option{geostrophic\_pressure\_option} to include both Coriolis and
    buoyancy in the balance. The \option{\ldots/solver} options can be chosen
    the same as for Pressure. If your problem includes a
    \option{type::free\_surface} boundary condition under Velocity, you should
    add a zero \option{type::dirichlet} boundary condition under the
    Geostrophic field at the same free surface boundary. If not, leave all
    boundaries unspecified (this implies a Neumann condition). In that case
    however, you need the option \option{remove\_null\_space} under
    \option{\ldots/solver}.
\end{itemize}

\subsection{Control Volumes for scalar fields}
For smaller scale problems using Control Volumes (CV) may be a robust and efficient
option for solving the scalar advection diffusion equations. For larger scale problems in which maintaining stratification is
important it is not recommended. The recommended options for scalar fields 
discretised with CV are:
\begin{itemize}
  \item The mesh used should be the \option{CoordinateMesh} and the equation type \option{AdvectionDiffusion}, cf. \ref{sec:ND_advection_diffusion_discretisation}. 
  \item For the spatial discretisation a control-volumes discretisation with a finite-element face value discretisation and Sweby limiter are recommended which are selected with the options, cf. \ref{sec:CVs}:
\begin{itemize}
\item \option{\ldots/spatial\_discretisation/control\_volumes/face\_value::FiniteElement}
\item \option{\ldots/prognostic/spatial\_discretisation/control\_volumes/face\_value::FiniteElement/limit\_face\_value/limiter::Sweby}
\end{itemize}
\item To help increase speed it is possible to store upwind elements so they do
  not have to be recalculated every time step (only after adapts). To do this
  activate the option
\option{\ldots/spatial\_discretisation/control\_volumes/face\_value::FiniteElement/limit\_face\_value/limiter::Sweby/project\_upwind\_value\_from\_point/store\_upwind\_elements}.
\item For diffusion, the Element Gradient diffusion scheme is recommended, selected under \option{\ldots/spatial\_discretisation/control\_volumes/diffusion\_scheme::ElementGradient}
\item For the temporal discretisation a Crank-Nicolson scheme is recommended, with the control volume options of 3 advection iterations and limit theta, cf. \ref{sec:CVs}. These are set with the options:
\begin{itemize}
\item{\option{\ldots/temporal\_discretisation/theta} set to 0.5}
\item{\option{\ldots/temporal\_discretisation/control\_volumes/number\_advection\_iterations}}
\item{\option{\ldots/prognostic/temporal\_discretisation/control\_volumes/limit\_theta}}
\end{itemize}
\end{itemize}

\subsection{Discontinuous Galerkin for scalar fields}
When using \PoDGPt for Velocity and Pressure, using P1DG for the scalar fields
is also an option that may be more accurate but is also more expensive. If stratification is important however, there is a known
issue with maintaining stratification near the boundaries with the DG slope
limiter.
\begin{itemize}
  \item The mesh used should be the \option{VelocityMesh} and the equation type \option{AdvectionDiffusion}, cf. \ref{sec:ND_advection_diffusion_discretisation}. 
  \item For the spatial discretisation a discontinuous-Galerkin discretisation
is recommended with a Lax-Friedrichs advection scheme, velocity
projected to continuous space, advection integrated by parts once and
a compact-discontinuous-Galerkin diffusion scheme with a vertex-based
slope limiter. These options are selected with:
\begin{itemize}
\item{\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/advection\_scheme/lax\_friedrichs}, cf. \ref{sec:ND_discontinuous_galerkin_advection}}
\item{\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/advection\_scheme/project\_velocity\_to\_continuous/mesh::CoordinateMesh}}
\item{\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/advection\_scheme/integrate\_advection\_by\_parts/once}, cf. \ref{sec:ND_discontinuous_galerkin_advection}}
\item{\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/diffusion\_scheme/compact\_discontinuous\_galerkin}, cf. \ref{sec:NM_DG_diffusion}}
\item{\option{\ldots/spatial\_discretisation/discontinuous\_galerkin/slope\_limiter::Vertex\_Based}, cf. \ref{sec:ND_DG_slope_limiters}.}
\end{itemize}

\item For the temporal discretisation a Crank-Nicolson scheme with subcycling is recommended. This can be set with:
\begin{itemize}
\item{\option{\ldots/temporal\_discretisation/theta} set to 0.5}
\item{\option{\ldots/temporal\_discretisation/discontinuous\_galerkin/maximum\_courant\_number\_per\_subcycle} set to an appropriate value.}
\end{itemize}
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Shallow Water Equations           %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Shallow Water Equations} \label{sec:shallow_water_equations}
\index{shallow water equations!configuration}
To solve for the shallow water equations \eqref{eq:swe} select equation type
\option{equation::ShallowWater} under a prognostic Velocity field. You will
also need a prognostic Pressure field which solves the barotropic pressure
$g\eta$, where $\eta$ is the Free Surface elevation. 

The setup of the temporal
and spatial discretisation options is mostly the same as for 3D ocean
flow problems (see section \ref{sec:large-scale-ocean}). We recommend the \PoDGPt pressure element pair. The
main differences with a 3D setup are:
\begin{itemize}
  \item Under Velocity, the equation type should be set to
    \option{equation::ShallowWater}. Here you should specify the bottom depth
    of your domain as a prescribed \option{scalar\_field::BottomDepth}.
  \item The problem is solved on a 2D mesh. This typically means that the
    CoordinateMesh is set as the input mesh (with the \option{from\_file}
    option). As before, the VelocityMesh and PressureMesh are derived from it (using
    \option{from\_mesh}) with the options \option{mesh\_continuity} set to
    \option{discontinuous} for the VelocityMesh and 
    \option{mesh\_shape/polynomial\_degree} set to 2 for the PressureMesh.
    The option \option{/geometry/ocean\_boundaries} should not be used.
  \item Under \option{/physical\_parameters} you still need to set
    \option{gravity/magnitude}. The settings under
    \option{gravity/vector\_field::GravityDirection} will be ignored.
  \item A Coriolis force may be added under \option{/physical\_parameters} with
    either the option \option{f\_plane} or \option{beta\_plane}.
  \item Since density is not taken into account in the shallow water equations
    that are solved, equation \eqref{eq:swe}, you should not have an
    \option{equation\_of\_state} or a \option{scalar\_field::Density} under
    \option{/material\_phase}.
  \item Since the shallow water equations are solved in non-conservative form,
    the momentum advection term needs to be in non-conservative form. Set
    \option{\ldots/spatial\_discretisation/conservative\_advection} to 0.0
    under Velocity.
  \item No boundary conditions for the bottom or top free surface should be
    added under the Velocity field. Instead, to add a bottom drag select this
    option under \option{equation::ShallowWater}. Lateral boundary conditions
    are set as previously. As before, a free surface boundary condition is
    applied as a Pressure boundary condition, where the desired free surface
    elevation is multiplied by $g$.
  \item For Pressure, the option \option{\ldots/solver/iterative\_method} should
    be set to \option{iterative\_method::gmres} instead of \option{cg}.
    The \option{vertical\_lumping} under \option{preconditioner::mg} is not
    valid in 2D domains. The \option{preconditioner::mg} is still recommended.
  \item A diagnostic \option{scalar\_field::FreeSurface} (should be on the
    PressureMesh) can be added to
    visualise the free surface elevation $\eta$. This is calculated by dividing
    the Pressure value by $g$.
\end{itemize}

There is currently no capability to solve the shallow water equations on the
sphere. For ocean-scale and global models, where one is only interested in the
barotropic, depth-averaged flow, you need to set up the problem using a 3D
configuration as described in section \ref{sec:large-scale-ocean} using a
single layer.

\section{Mesh adaptivity} \label{sec:config_adapt}

The configuration on mesh adaptivity occurs in two places: under \option{mesh\_adaptivity} where 
the overall adaptive settings are configured, and on a per-field basis where both the interpolation
method is set and if that field should be considered when creating the error metric. See chapter 
\ref{chap:Adaptivity} for the background to adaptivity and more detailed information.

\subsection{Field settings}

For each field present in the simulation there are up to two options that should be set. The first
is the interpolation method that should be used to transfer the values of a field from the old to the new mesh, section \ref{sec:interpolation_algorithms}.
Second, in order to form the error metric by which the mesh is adapted, section \ref{sec:norm_choice}, the user must set which fields
should form the error metric and how the error for that field should be calculated.

\subsubsection{Interpolation method} \label{sec:config_adapt_interp}

For each prognostic field in the current state, an interpolation type, section \ref{sec:interpolation_algorithms}, must be set. These can be set by selecting an
option \option{\ldots/prognostic/<interpolation type>} where \texttt{<interpolation type>} is one of:
\begin{itemize}
\item Consistent interpolation - the default and quick interpolation method, but is non-conservative and dissipative.
\item Pseudo-consistent interpolation - not recommended at present.
\item Galerkin interpolation - Conservative and non-dissipative, but requires the construction of a supermesh \citep{farrell2009a,farrell2010a}
\item Grandy interpolation - Conservative, but highly diffusive. See \citet{grandy1999}.
\end{itemize}

For some fields, such as Pressure and Velocity other interpolation methods are available.

For diagnostic and prescribed fields an interpolation method is not required. However, if an output dump
occurs immediately following a mesh adapt, diagnostic fields may not have correct values depending
on the method by which they are calculated. In these instances, it is worth setting an interpolation type
for these fields which will ensure that the values are set correctly before an output dump occurs.

The Galerkin projection also requires some further settings depending on the mesh type. For discontinuous
meshes there are no other required settings.  For continuous meshes a solver is required in order
to perform the supermesh projection. The solver settings are configured as with any other solver, see
section \ref{sec:Solve} for more details.

With piecewise linear continuous fields additional options are available to bound the result following a Galerkin projection:\\
\option{\ldots/galerkin\_projection/continuous/bounded::Diffuse}\\
The \option{Diffuse} bounding algorithm is internal to \fluidity\ and the most frequently used.

To use the \option{Diffuse} bounding algorithm \citep{farrell2009a} the number of iterations the algorithm is allowed to take must be specified.  Additionally an optional tolerance can be specified to terminate this iteration loop early.  Furthermore if the bounds on the field are known in advance then these can be specified through:\\
\option{\ldots/bounded::Diffuse/bounds/upper\_bound}\\
and\\
\option{\ldots/bounded::Diffuse/bounds/lower\_bound}.\\
If the diffusion bounding algorithm fails to locally redistribute the unboundedness then a conservative but non-local redistribution can be activated using:\\
\option{\ldots/bounded::Diffuse/repair\_deviations}\\
again with an optional tolerance:\\
\option{\ldots/bounded::Diffuse/repair\_deviations/tolerance}.

\subsubsection{Creating an error metric}
\label{sec:configuring_fluidity_error_metric}

The second step for configuring adaptivity is to set up the fields that are to form the error metric used
to adapt the mesh. For each field that should be considered when forming the metric the option
\option{\ldots/adaptivity\_options} needs to be enabled. The type or error norm on which the metric is based (absolute or relative) is set with \option{\ldots/adaptivity\_options/absolute\_measure} or \option{\ldots/adaptivity\_options/relative\_measure}. For a $p$--norm \option{\ldots/adaptivity\_options/absolute\_measure} should be selected and the value of $p$ set with the option \option{\ldots/adaptivity\_options/absolute\_measure/p\_norm} ($p=2$ is recommended).

The InterpolationErrorBound field must be set and as with any other prescribed field can take a constant value or vary in space and time (by prescribing a python function for example). The error bound is set as separate fields within 
state, so for Temperature, for example, the acceptable error is stored an a field called 
TemperatureInterpolationErrorBound. This field is output as any other field too.

For relative interpolation error bounds a tolerance value also has to be set under \newline
\option{\ldots/adaptivity\_options/relative\_measure/tolerance}. This value prevents division by zero
and should be set to a small enough number that the field can effectively be considered zero at this value.

For discussion of the different metrics and error norms see section \ref{sec:norm_choice}.

\subsection{General adaptivity options}
\label{sec:configuring_fluidity_adaptivity_options}
\index{adaptivity options}

These are found under \option{mesh\_adaptivity}. Here the user can specify whether to use mesh movement methods, 
prescribed adaptivity (serial only) or hr adaptivity. hr adaptivity is the normal method for most applications.

Under \option{/mesh\_adaptivity/hr\_adaptivity} there are a number of mandatory options, which are:
\begin{itemize}
\item \option{period} - how often should the mesh be adapted. This can be set in number of simulation seconds, 
or in number of timesteps. It is recommended that adapt happen every 10-20 timesteps.
\item \option{maximum\_number\_of\_nodes} - sets the maximum possible number of nodes in the domain. In parallel this is the
global maximum number, but can be altered to be the number per process. If the maximum number of nodes is reached the mesh is coarsened everywhere until this is achieved.
\item \option{enable\_gradation} - is on by default and set to a value of 1.5. This constrains the jump in desired edge lengths along an edge and therefore controls how fast the mesh size may change.
\item \option{tensor\_field::MinimumEdgeLengths} - a tensor specifying the minimum edge length of an element.
\item \option{tensor\_field::MaximumEdgeLengths} - a tensor specifying the maximum edge length of an element.
\end{itemize}

In addition to these mandatory settings, there are a number of other configuration options.
\begin{itemize}
\item \option{cpu\_period} - sets the time interval for the mesh adapt in cpu time.
\item \option{minimum\_number\_of\_nodes} - sets the minimum possible number of nodes in the domain. In parallel this is the
global minimum number. The mesh is refined until this is achieved.
\item \option{adaptive\_timestep\_at\_adapt} - used in conjunction with adaptive timestep (see section \ref{section:config_adaptive_timestep}), this
option resets the timestep back to the minimum value under \option{\ldots/adaptive\_timestep/minimum\_timestep} immediately following
a mesh adapt.
\item \option{maximum\_node\_increase} - the maximum ratio by which the number of nodes is allowed to increase. A value of 1.1
indicates the number of nodes can increase by at most 10\%.
\item \option{node\_locking} - allows the locking of nodes via a python function that cannot be moved by adaptivity.
\item \option{functional\_tolerance} - specifies the minimum element functional value for which elements are considered for adaptivity, section \ref{sec:adaptive_remeshing_technology}. Default value is $0.15$.
\item \option{geometric\_constraints} - this applies geometric constraints to the metric formation which aims to prevent the metric demanding edge length that are inappropriately large in comparison to the resolution required to preserve the geometric accuracy of the boundaries. If you get `knife elements' near the boundaries try turning this option on. This only works in 3D. 
\item \option{bounding\_box\_factor} - this option bounds the edge lengths requested by the metric by bounding box of the domain, multiplied by the specified factor. The default value is $2$.
\item \option{reference\_mesh} - supply a reference mesh which supplies the minimum or maximum edge length to the metric.
\item \option{aspect\_ratio\_bound} -  maximum aspect ration of elements in the adapted mesh.
\item adapt at first timestep - perform mesh adaptivity before the first timestep occurs. This can occur a specified number of times.
\item \option{preserve\_mesh\_regions} - ensures that regions in your mesh, specified by region IDs, are preserved through adaptivity.
is adapted, then the mesh is extruded using the adaptivity metric in the 3rd dimension. You must use an extruded mesh with this option, section \ref{sec:extruded_meshes}.
\item \option{adaptivity\_library} - choose which adaptivity library to use. In 2D you are restricted to libmba2d. In 3D you can choose either libmba3D or libadaptivity (default).
\item \option{adapt\_iterations} - this options controls the number of intermediate adapt iteration during parallel adaptive simulations, section \ref{sec:parallel_adaptivity}. The default value is 3. Higher values may give you 
better meshes, especially when the number of elements per process is low.
\item \option{debug} - options for output that is useful for debugging adaptivity.
\end{itemize}

\subsubsection{Vertically structured and 2+1D adaptivity}
\label{sec:vertically_structured_adaptivity}
For some problems it can be advantageous to apply adaptivity in the horizontal
and vertical as separate steps. This means a horizontal (surface) mesh is
adapted first after which a column of nodes is created under each surface node.
The resolution in the vertical columns is either specified under the extrusion
options, or determined via a vertical adaptivity step. This functionality is switched
on using the \option{vertically\_structured\_adaptivity} option. An extruded
initial mesh is required for this option (see section \ref{sec:extruded_meshes},
and section \ref{sec:extruded} for its configuration). The horizontal adaptivity
stage is then applied to the horizontal input mesh, and the
\option{bottom\_depth} and \option{sizing\_function} extrusion options are
 reapplied for the creation of the vertical columns.
If an extruded initial mesh is
provided only \option{vertically\_structured\_adaptivity} can be used. Adapting the
mesh is all dimension will not work correctly.
Further options under
\option{vertically\_structured\_adaptivity} are:
\begin{itemize}
\item \option{inhomogeneous\_vertical\_resolution} -
This option switches on vertical adaptivity. This means it will no longer create
layers based on the \option{sizing\_function} option. Instead, the distance
between the nodes in the vertical columns is based on the vertical component of
the error metric. The vertical resolution will therefore vary over the depth and 
in each column independently. With the combination of
\option{vertically\_structured\_adaptivity} and
\option{inhomogeneous\_vertical\_resolution}, adaptivity can thus focus resolution
in both horizontal and vertical, while maintaining a columnar nodal structure.
This combination is referred to as 2+1D adaptivity.
\begin{itemize}
\item \option{adapt\_in\_vertical\_only} - With this option vertical adaptivity is
applied, but the horizontal mesh is kept fixed.
\end{itemize}
\item \option{split\_gradation} - Instead of applying gradation to the full metric
before splitting into a horizontal and vertical metric, with this option the
gradation is applied after the split. Thus in particular when specifying
anisotropic gradation, the gradation in the horizontal and vertical is applied
completely independently.
\item \option{vertically\_align\_metric} or \option{use\_full\_metric} -
The metric applied in the horizontal adaptivity stage is assembled by merging
the 3D metric in each column and then projecting to the horizontal plane. 
Typically the 3D metric for large aspect ratio problems already decomposes 
in an (almost) vertical eigenvector and 2 horizontal ones. 
However, even the slightest tilt causes vertical error bounds 
to be ``leaked'' into the merged horizontal metric, leading to unexpected small
horizontal edge lengths. Therefore for large aspect ratio problems the option
\option{vertically\_align\_metric}, which decomposes the metric \emph{before}
merging in the horizontal, is recommended.
\item \option{include\_bottom\_metric} - When constructing the horizontal metric
incorporate the components of the full metric tangential to the bottom boundary.
For example, this is useful when horizontal contours of a field intersect the
bathymetry and this information is not automatically incorporated into the
horizontal metric leading to the contact point being under-resolved.
\end{itemize}

\subsubsection{Zoltan options}
\label{sec:configuring_fluidity_zoltan_options}
\index{adaptivity options}
\index{zoltan}

There are a number of options available for controlling Zoltan's behaviour when re-partitioning the mesh during
and after adaptivity, which can be found under \option{mesh\_adaptivity/zoltan\_options}. The options are:
\begin{itemize}
\item partitioner - this is the partitioner used in the intermediate adapt iterations. It can be one of Scotch, ParMetis, Zoltan, or Zoltan Hypergraph. Default is Zoltan.
\item final partitioner - the partitioner used for the final adapt iteration where load balancing is important. Same choices as above. Default is ParMetis. 
\item element quality cutoff - at what value of element quality is an element deemed ``bad''. Default is 0.6.
\item load imbalance tolerance - a value of 1 means each processor will have exactly the same number of elements. However, smaller numbers here mean that the intermediate adapts may
not be able to move the mesh sufficiently to get a good quality mesh from adaptivity.
\item additional adapt iterations - increases the number of intermediate adapt iterations during parallel adaptive simulations.
\item zoltan debug - debugging options.
\end{itemize}

For more information on the approach to parallel adaptivity adopted in \fluidity\ see section \ref{sec:parallel_adaptivity}.
\subsubsection{Metric advection}
\label{sec:configuring_fluidity_metric_advection}
Metric advection advects the metric along with the flow, ensuring the resolution can be pushed ahead of any flow, rather than lagging behind, section \ref{section:metric_advection_general}. The advection equation is discretised with a control volume method,  section \ref{sec:ND_control_volume_advection}. For spatial discretisation a first order upwind scheme for calculation the face values (the default) and non--conservative form are generally recommended. These are selected with options
\begin{itemize}
\item \option{/mesh\_adaptivity/hr\_adaptivity/metric\_advection/spatial\_discretisation/control\_volumes/face\_value::FirstOrderUpwind}
\item \option{/mesh\_adaptivity/hr\_adaptivity/metric\_advection/spatial\_discretisation/conservative\_advection} $=0.0$
\end{itemize}
For temporal discretisation a semi--implicit discretisation in time is recommended, section \ref{sec:ND_time_disct_adv_diff}, with option
\begin{itemize}
\item \option{/mesh\_adaptivity/hr\_adaptivity/metric\_advection/temporal\_discretisation/theta} $=0.5$
\end{itemize}
The time step is controlled by the choice of CFL number, specified in \option{/mesh\_adaptivity/hr\_adaptivity/metric\_advection/temporal\_discretisation/maximum\_courant\_number\_per\_subcycle}. The metric is advected over the time period between the current and the next adapt. This time period can be scaled with the option \option{/mesh\_adaptivity/hr\_adaptivity/metric\_advection/temporal\_discretisation/scale\_advection\_time} which has a default value of $1.1$.

\section{Multiple material/phase models} \label{sec:config_multimatph}

This section contains advice on setting up simulations with multiple material\_phase options.  This enables related fields to be grouped together into related materials or phases.  For example a prognostic scalar field in one material\_phase will be advected using the Velocity field from that material\_phase, while a prognostic scalar field in another material\_phase will be advected according to the Velocity field in its material\_phase.

We refer to two typical scenarios: a \emph{multiple material} model and a \emph{multiple phase} model. A multiple phase model is one in which the Velocity field in each material\_phase is in some way independent of the velocities in the other material\_phases.  This means that scalar fields (for example phase volume fractions) in each material\_phase are advected independently.  A multiple material model is one in which the Velocity field is shared between all material\_phases so that all scalar fields are advected similarly.

\subsection{Multiple material models}

Models with a single prognostic velocity field that is shared between material\_phases (using the \option{aliased} field type) are referred to as multiple material models.  These are generally used to describe systems of nearly immiscible materials with different material properties contained within the same domain.  This section focusses on this type of multiple material\_phase simulation.

In a multiple material simulation each material\_phase requires:
\begin{itemize}
\item an equation of state, and
\item a MaterialVolumeFraction scalar field.
\end{itemize}

The equation of state provides the density properties of the material described in the current material\_phase.  For incompressible simulations a linear equation of state is used, which only requires a reference density:\\
\option{\ldots/equation\_of\_state/fluids/linear/reference\_density}\\
to be set.  For Boussinesq multimaterial simulations, where a material's density depends on temperature and/or salinity, then the same dependencies exist between the equation of state and these fields as in single material simulations.  For example a Temperature field must be present in the material\_phase where it is needed (although it may be aliased between material\_phases).  If the \option{subtract\_out\_hydrostatic\_level} option is selected, it must only be select in a single material\_phase.  Fully compressible multimaterial simulations are not supported.

The MaterialVolumeFraction field describes the location of the material, varying from $1$ in regions where the cells are entirely the current material to $0$ where none of this material is present.  As the materials are generally treated as being nearly immiscible, the prognostic MaterialVolumeFraction field should be discretised using a control volume spatial discretisation with one of the face value schemes designed for advecting step functions:
\begin{itemize}
\item \option{\ldots/spatial\_discretisation/control\_volumes/face\_value::HyperC}
\item \option{\ldots/spatial\_discretisation/control\_volumes/face\_value::UltraC}
\item \option{\ldots/spatial\_discretisation/control\_volumes/face\_value::PotentialUltraC}
\end{itemize}
as described in sections \ref{sec:hyperc}--\ref{sec:potultrac}.  These schemes are only guaranteed to be bounded for explicit advection so the implicitness factor, $\theta$:\\
\option{\ldots/temporal\_discretisation/theta}\\
and the pivot implicitness factor, $\theta_p$:\\
\option{\ldots/temporal\_discretisation/control\_volumes/pivot\_theta}\\
should be set to zero and, for high Courant number flows, advection subcycling should be used:\\
\option{\ldots/temporal\_discretisation/control\_volumes/maximum\_courant\_number\_per\_subcycle}\\
or\\
\option{\ldots/temporal\_discretisation/control\_volumes/number\_advection\_subcycles}.

For an $N$ material problem, $N$ material\_phases are required and hence $N$ MaterialVolumeFractions, $c^i, i =1, \ldots, N$, and $N$ equations of state.  However, only $N-1$ of the MaterialVolumeFraction fields, $c^i, i =1, \ldots, N-1$, need be prognostic.  The final volume fraction field, $c^N$, should always be set to \option{diagnostic}, as it can be recovered using the internal algorithm:
\begin{equation} \label{eq:config_diagvfrac}
c^N = 1 - \sum_{i=1}^{N-1}c^i\textrm{.}
\end{equation}
So, for example, in the case when $N=2$ there need only be a single prognostic MaterialVolumeFraction field and a single diagnostic MaterialVolumeFraction field.  In this case it makes no difference which material\_phase contains the prognostic volume fraction and which contains the diagnostic field.  In more complicated scenarios with $N>2$ a coupled control volume discretisation (see section \ref{sec:coupledlimiter}) becomes necessary to ensure that not only each of the $N-1$ prognostic MaterialVolumeFractions remain bounded but also that their sum, $\sum_{i=1}^{N-1}c^i$, is bounded.  This ensures, through equation \ref{eq:config_diagvfrac}, that the final diagnostic MaterialVolumeFraction is also bounded.  As discussed in section \ref{sec:coupledlimiter} this process requires a priority ordering for the fields, which must be specified at:\\
\option{\ldots/scalar\_field::MaterialVolumeFraction/prognostic/priority}.\\
The diagnostic field is always treated as the lowest priority volume fraction so in this case the choice of priority ordering and diagnostic field may affect the results if the interfaces between the materials are in the vicinity of one another.  Priority ordering and coupled limiting do not affect the advection process if the material interfaces are separated from each other.

If adaptive remeshing is being used then the bounded and minimally dissipative behaviour of the above advection must be preserved through the interpolation between successive meshes.  As discussed in section \ref{sec:config_adapt} several interpolation algorithms are available.  We discuss them again here in terms of their suitability for multiple material modelling.  Consistent interpolation on piecewise linear parent meshes guarantees boundedness of the interpolated volume fraction field and of the sum of the volume fractions.  However it tends to introduce excessive amounts of numerical diffusion and it is not conservative.  Galerkin projection guarantees conservation of the field and is not excessively dissipative.  However it does not guarantee boundedness.

To ensure minimal dissipation, conservation and boundedness it is necessary to use a bounding algorithm following the Galerkin projection.  The \option{Diffuse} bounding algorithm (see section \ref{sec:config_adapt_interp}) is generally used.  This redistributes unbounded values in the field locally, guaranteeing boundedness of each volume fraction individually.  It does not however guarantee boundedness of the sum of the volume fractions and this must be enforced by coupling each MaterialVolumeFraction field together through the interpolation with the option:\\
\option{\ldots/bounded::Diffuse/bounds/upper\_bound/coupled}\\
under all $N-1$ prognostic MaterialVolumeFractions.  As with coupled control volume advection this uses the priority numbering of the fields to determine the order in which they are bounded.  The local bounds enforced on successive fields are then modified to ensure boundedness of their sum.  This redistribution of materials during the bounding procedure introduces some relative movement between materials, which, by equation \ref{eq:config_diagvfrac}, is filled in by the diagnostic MaterialVolumeFraction.  Despite this problem bounded Galerkin projection is recommended to transfer field data during mesh adaptivity.

The advected (and interpolated) volume fractions describe volume averaged the locations of the materials.  In combination with the equation of state they can therefore be used to define global bulk values for the density using:
\begin{equation}
\rho = \sum_{i=1}^N\rho^ic^i
\end{equation}
where $\rho$ is the bulk density and $\rho^i$ are the individual material densities, given by their respective equations of state.  This bulk density can be seen in the diagnostic Density field in whichever material\_phase has the prognostic Velocity field in it.

In addition to the density the volume fractions may be used to specify a bulk Viscosity that varies between the materials according to:
\begin{equation}
\tensor{\mu} = \sum_{i=1}^N\tensor{\mu}^ic^i
\end{equation}
where $\tensor{\mu}$ is the bulk viscosity and $\tensor{\mu}^i$ is the individual material's viscosity.  To use this it is necessary to activate a MaterialViscosity field in every material\_phase with a nonzero viscosity.  Additionally the Viscosity field underneath the prognostic Velocity must be activated and set to the \option{bulk\_viscosity} diagnostic algorithm.

\subsection{Multiple phase models}
Models with one prognostic velocity field per \option{material\_phase} are referred to as multiple phase models. The use of these multiple velocity fields permits the inter-penetration and interaction between different phases.

\subsubsection{Simulation requirements}
In a multiphase simulation, each \option{material\_phase} requires:
\begin{itemize}
 \item an equation of state, and
 \item a PhaseVolumeFraction scalar field.
\end{itemize}

As per multi-material simulations, the equation of state provides the density properties of the phase described in the current \option{material\_phase}. For incompressible simulations a linear equation of state is used, which only requires a reference density:\\
\option{\ldots/equation\_of\_state/fluids/linear/reference\_density}\\
to be set.

For an $N$ phase problem, $N$ material\_phases are required and hence $N$ PhaseVolumeFractions, $\alpha_i, i = 1, \ldots, N$, and $N$ equations of state. Just as in multi-material simulations, only $N-1$ of the PhaseVolumeFraction fields, $\alpha_i, i = 1, \ldots, N-1$, need be prognostic. The final PhaseVolumeFraction field, $\alpha_N$, should always be set to \option{diagnostic}, as it can be recovered using the internal algorithm:
\begin{equation}
\alpha_N = 1 - \sum_{i=1}^{N-1}\alpha_i\textrm{.}
\end{equation}
For compressible multiphase simulations, the diagnostic PhaseVolumeFraction field should always be in the compressible (fluid) phase. This is because we do not include the Density in the advection-diffusion equation for prognostic PhaseVolumeFraction fields, so solving this equation in the compressible phase would not be correct. The particle phases on the other hand are always incompressible where the density is constant.

\subsubsection{Inter-phase momentum interaction}
Currently, \fluidity\ only supports fluid-particle drag between the continuous phase and dispersed phase(s). This can be enabled by switching on the \option{/multiphase\_interaction/fluid\_particle\_drag} option in Diamond and specifying the drag correlation: Stokes, \cite{wen_yu_1966}, \cite{ergun1952}, \cite{schiller1935drag}, \cite{lain1999experimental} or \cite{lain2002modelling}.

The drag force using the Stokes drag correlation is given by:
\begin{equation}\label{eq:stokes_drag_force}
\mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d},
\end{equation}
where $f$ and $p$ denote the fluid (i.e. continuous) and particle (i.e. dispersed) phases respectively, and $d$ is the diameter of a single particle in the dispersed phase. The drag coefficient $C_D$ is defined as:
\begin{equation}\label{eq:stokes_drag_coefficient}
C_D = \frac{24}{\mathrm{Re}},
\end{equation}
with
\begin{equation}\label{eq:particle_reynolds_number}
\mathrm{Re} = \frac{\alpha_f\ \rho_f\ d\ |\mathbf{u}_f-\mathbf{u}_p|}{\mu_f}.
\end{equation}
Note that $\mu_f$ denotes the isotropic viscosity of the fluid (i.e. continuous) phase.

With the drag correlation by \cite{wen_yu_1966}, $\mathbf{F}_D$ and $C_D$ become:
\begin{equation}\label{eq:wen_yu_drag_force}
\mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d\alpha_f^{2.7}},
\end{equation}
\begin{equation}\label{eq:wen_yu_drag_coefficient}
C_D = \frac{24}{\mathrm{Re}}\left(1.0 + 0.15\mathrm{Re}^{0.687}\right).
\end{equation}
In contrast to the Stokes drag correlation, the \cite{wen_yu_1966} drag correlation is more suitable for larger particle Reynolds number flows. 

For dense multiphase flows with $\alpha_p > 0.2$, the drag correlation by \cite{ergun1952} is often the most appropriate:
\begin{equation}\label{eq:ergun_drag_force}
\mathbf{F}_D = \left(150\frac{\alpha_p^2\mu_f}{\alpha_f d_p^2} + 1.75\frac{\alpha_p\rho_f|\mathbf{u}_f-\mathbf{u}_p|}{d_p}\right)\left(\mathbf{u}_f-\mathbf{u}_p\right).
\end{equation}

The interphase drag force $\mathbf{F}_D$, and the drag coefficient $C_D$ for the standard \cite{schiller1935drag} drag correlation is given by:
\begin{equation}\label{eq:schillernaumann_drag_force}
\mathbf{F}_D = \frac{3\ \alpha_p\ C_D\ \alpha_f\ \rho_f\ |\mathbf{u}_f-\mathbf{u}_p|\ (\mathbf{u}_f-\mathbf{u}_p)}{4\ d},
\end{equation}
\begin{equation}\label{eq:schillernaumann_drag_coefficient}
C_D=
\begin{cases}
  \dfrac{24}{Re} \left( 1+0.15 {Re}^{0.687} \right)& \text{if}\ Re<1000, \\[2ex]
  0.44                                               & \text{otherwise},
\end{cases}
\end{equation}
with a slightly different definition for the Reynolds number than above:
\begin{equation}\label{eq:particle_reynolds_number_2}
\mathrm{Re} = \frac{\rho_f\ d\ |\mathbf{u}_f-\mathbf{u}_p|}{\mu_f}.
\end{equation}
The above drag force correlation was developed for a spherical and solid dispersed phase but dispersed fluids with slight impurities can be modelled using the same correlation. \cite{schiller1935drag} drag has been used for the modelling of dilute bubble columns successfully. 

Note that within each dispersed phase, a value for $d$ must be specified in \option{../multiphase\_properties/particle\_diameter} regardless of which drag correlation is used. A scalar field can also be used instead of a constant dispersed phase diameter by choosing the \option{../multiphase\_properties/particle\_dia\_use\_scalar\_field} option instead. The latter option is useful when the multiphase model is coupled to the population balance equation \citep{bhutani2016polydispersed}.

Additionally, it is also possible to cap the dispersed phase diameter, which is useful when using a scalar field for the diameter. For instance, a very small diameter can result in a very large drag force which can result in a failing simulation.

\subsubsection{Inter-phase energy interaction}
This subsection considers compressible multiphase flow simulations where each phase has its own InternalEnergy field. Users can choose to include the inter-phase heat transfer term by \cite{gunn1978}, denoted $Q$ in Section \ref{sec:multiphase_equations}, on the RHS of each internal energy equation:

\begin{equation}\label{eq:gunn_heat_transfer_term}
Q = \frac{6 k \alpha_p \mathrm{Nu}_p}{d_p^2}\left(\frac{e_f}{C_{v,f}} - \frac{e_p}{C_{v,p}}\right),
\end{equation}
where 
\begin{equation}
\mathrm{Nu}_p = \left(7 - 10\alpha_f + 5\alpha_f^2\right)\left(1 + 0.7\mathrm{Re}_p^{0.2}\mathrm{Pr}^{\frac{1}{3}}\right) + \left(1.33 - 2.4\alpha_f + 1.2\alpha_f^2\right)\mathrm{Re}_p^{0.7}\mathrm{Pr}^{\frac{1}{3}},
\end{equation}
\begin{equation}
 \mathrm{Pr} = \frac{C_{v,f} \gamma \mu_f}{k},
\end{equation}
and
\begin{equation}
 \mathrm{Re} = \frac{\rho_f d_p |\mathbf{u}_f-\mathbf{u}_p|}{\mu_f},
\end{equation}
are the Nusselt, Prandtl and Reynolds number respectively. $C_{v,i}$ denotes the specific heat of phase i at constant volume, and $k$ denotes the effective conductivity of the fluid phase; these must be specified in \option{../multiphase\_properties/effective\_conductivity} and \option{../multiphase\_properties/specific\_heat}. Note that we have written the above in terms of internal energy (rather than temperature) by dividing $e_i$ by the specific heat at constant volume.

To include this term, switch on the \option{/multiphase\_interaction/heat\_transfer} option in Diamond.

\subsubsection{Current limitations}
\begin{itemize}
 \item Boussinesq multiphase simulations are not yet supported.
 \item The momentum equation for each \option{material\_phase} can only be discretised in non-conservative form.
 \item \option{bassi\_rebay} and \option{compact\_discontinuous\_galerkin} are currently the only DG viscosity schemes available for multiphase flow simulations.
 \item Discontinuous \option{PhaseVolumeFraction} fields are not yet supported.
 \item For compressible multiphase flow simulations, the Pressure field only supports a \option{continuous\_galerkin} discretisation.
 \item Prescribed velocity fields cannot yet be used in multiphase simulations.
 \item Fluid-particle drag can currently only support one continuous (i.e. fluid) phase.
\end{itemize}


\section{Compressible fluid model}
\index{compressible fluid model}

Enabling \option{\ldots/material\_phase/equation\_of\_state/compressible} allows the compressible equations described in sections 
\ref{sec:compressible_conservative} and \ref{sec:compressible_nonconservative} to be solved. At the moment there is one available option for the required 
compressible equation of state: Mie-Gr\"uneisen (see section \ref{sec:Multi-material compressible EOS}).  Compressible functionality is not yet fully supported and this is intended as a stub upon which further developments will be described.

This section contains advice for running a compressible simulation, by describing the necessary options to set up the problem. The options required for prognostic fields are:

\subsection{Pressure options}
\begin{itemize}
\item The value for the atmospheric pressure can be added by switching on \option{\ldots/atmospheric\_pressure}, otherwise a default of zero is used.
\item A Poisson pressure equation should not be used to calculate a first guess, therefore \\* \option{\ldots/scheme/poisson\_pressure\_solution} should be set to $never$.
\item \option{\ldots/scheme/use\_compressible\_projection\_method} should be selected, so the calculated pressure satisfies the continuity equation and the EOS.
\end{itemize}

\subsection{Density options}
\option{\ldots/prognostic/equation} and \option{\ldots/prognostic/solver} do not need to be enabled. If the equation type is not turned on, the density will make 
use of the pressure solve, so no solver options are needed either. By having an equation type turned on, the density is not only incorporated into the pressure solve 
but also an Advection-Diffusion equation is solved (and solver options need to be specified).

\subsection{Velocity options}\label{sec:VelocityOptions}
\subsubsection{Mass lumping}
For continuous velocities \option{\ldots/spatial\_discretisation/\ldots/lump\_mass\_matrix} should be turned on.
\subsubsection{Selecting the correct form of the viscous term}
In the presence of viscosity a form of the stress tensor must be chosen, as discussed in \ref{sec:ViscosityTerms}. 
\begin{itemize}
\item for compressible flow the full '\emph{stress form}' is required: 
\begin{itemize} 
\item CG discretisation: \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/stress\_form}
\item DG discretisation: not available
\end{itemize}
\item for incompressible flow with spatially varying viscosity the '\emph{partial stress form}' is required:
\begin{itemize} 
\item CG discretisation: \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/partial\_stress\_form}
\item DG discretisation: only available when using the Bassi-Rebay formulation of the viscosity term and for isotropic viscosity \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/bassi-rebay/partial\_stress\_form}
\end{itemize}
\item for incompressible flow with homogeneous viscosity, the '\emph{tensor form}' is valid:
\begin{itemize} 
\item CG discretisation: \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/tensor\_form}
\item DG discretisation: \option{\ldots/spatial\_discretisation/continuous\_galerkin/stress\_terms/bassi-rebay/tensor\_form} or any other viscosity term formulation.
\end{itemize}
\end{itemize}

\emph{Important note when using 'stress form' and 'partial stress form' with a CG velocity field:} when using isotropic viscosities, and not the 'tensor form' of the viscosity term, all components of viscosity must be set to equal to the isotropic viscosity due to the method of implementation. If using an anisotropic viscosity with the 'partial stress form' or 'stress form' consult the schema and query the code to understand how the viscosity tensor must be defined. This does not apply when using DG discretisations. 

\subsubsection{Hydrostatic balance}
The option \option{subtract\_out\_reference\_profile} splits up the Density and Pressure fields into a hydrostatic (reference) component (') and a perturbed component (''). The hydrostatic (reference) components, denoted $p^{\prime}$ and $\rho^{\prime}$, should satisfy the balance:
\begin{equation}
  \nabla p^{\prime} = \rho^{\prime}\mathbf{g}
\end{equation}
Enabling this option will subtract the hydrostatic components, specified here, from the pressure and density used in the pressure gradient and buoyancy terms in the momentum equation. This helps to maintain hydrostatic balance and prevent spurious oscillations in the pressure field when using unbalanced finite element pairs.

To use this option you will also need to create two prescribed scalar fields, called HydrostaticReferencePressure and HydrostaticReferenceDensity, which define $p^{\prime}$ and $\rho^{\prime}$. These must be on the same mesh as pressure and density, respectively, and are meant to be set up as vertical profiles (i.e. constant in the horizontal direction). 

Note that unlike all the other hydrostatic/geostrophic balance options in Fluidity (i.e. \option{subtract\_out\_hydrostatic\_level} under the linear incompressible EoS option, or with HydrostaticPressure or GeostrophicPressure fields), the hydrostatic pressure is not subtracted from the Pressure field itself. In other words, the Pressure field that gets solved for (and output in the .vtu files) is still the combined Pressure ($p = p^{\prime} + \rho^{\prime}$), and the hydrostatic pressure $p^{\prime}$ is only subtracted in the momentum equation.

\subsection{InternalEnergy options}
\begin{itemize}
\item The optional heat flux term can be included via the InternalEnergy's Diffusivity field. Users will need to specify an isotropic value for $\frac{k}{C_v}$, where $k$ is the effective conductivity and $C_v$ is the specific heat at constant volume.
\end{itemize}

\subsection{Restrictions: discretisation options and element pairs}
Either continuous Galerkin or control volumes can be used as discretisation options for pressure and density (both fields need to have the same option). When using control volumes pressure and density have to be on the same order of parent mesh.