<commit_before>\chapter{Examples}
\label{chap:examples}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------WRITING A TEST CASE--------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Writing a test case} 

\subsection{Options file validation}

Since the grammar of \fluidity\ options files is defined in terms of a RelaxNG
schema, it is possible to check whether a given options file is valid. The
medium test ``diamond\_validation'' validates all test problem options files
with the .flml file extension. Hence if the options schema is updated in such
a way that an options file becomes invalid, the diamond\_validation test will
flag an error. Test problem XML files are also tested for validation against a
RelaxNG schema defining their format.

Some validation failures, such as the addition of new compulsory options, can
be repaired automatically using Diamond. The \lstinline[language = Bash]+update_options+
script can be used to attempt an automatic update of all test options files.

\subsection{The method of manufactured solutions}\label{sect:mms}
\index{method of manufactured solutions}
Problems with analytical solutions provide a powerful mechanism for testing the implementation
of numerical algorithms (code verification). The problem can be run on a series of meshes
with increasing resolution (\eg $\Delta x \rightarrow \Delta x/2$), or decreasing time step,
and the exact error between numerical and analytical solution computed in any norm. The theoretical
order of convergence of the algorithm can then be confirmed for the implentation code. Any significant
deviation from the theoretical order of convergence can be identified and is often due to coding bugs
or because the treatment of boundary conditions is reducing accuracy.

However, analytical solutions are not generally available, especially for complex equations sets
and for problems with physical realism. These are of course the problems that we are primarily 
interested in tackling numerically.
The method of manufactured solutions (MMS) takes the opoosite approach to deriving a solution to 
a given equation set: it starts from a given solution and derives the equations that admit this
solution. This achieved by substituting the solution into the eqution set and evaluating the residual
symbolically. This residual is then taken as a prescribed source term in the equation set being considered
and the identical procedure used for verification with an analytical solution used. The flexibility of the MMS
approach means that individual terms and discretisations can be rigorously tested through a suite of MMS problems. 

This procedure is especially useful in \fluidity\ as it is very easy to prescibe source terms in the equations
being solved. When substituting in a relatively simple analytical form for a solution into the governing
equations one will typically obtain a very complex expression for the residual or source term. While this
can result in a problem for many codes, in \fluidity\ one is able to prescribe anaytical source terms easily
within Diamond's Python function interface. The complex source term can be obtained using a Python symbolic algebra
library such as SAGE and then simply cut-and-paste into the Diamond file.

Some background material on the method of manufactured solutions, and guidelines for selecting good solutions,
is given in \cite{salari2000}. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------THE LOCK-EXCHANGE---------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The lock-exchange}
\label{sect:lock_exchange}
\subsection{Overview}
\label{sect:lock_exchange_overview}

The lock-exchange is a classic laboratory-scale fluid dynamics problem, \citep{fannelop_94, huppert_06, simpson_87}. A flat-bottomed tank is separated into two portions by a vertical barrier. One portion, the `lock', is filled with the source fluid. This is of different density to the ambient fluid which fills the other portion. As the barrier is removed, the denser fluid collapses under the lighter. Two gravity currents form and propagate in opposite directions, one above the other, along the tank. After an initial acceleration, the gravity current fronts travel at a constant speed until the end walls exert an influence or viscous forces begin to dominate, \citep{cantero_07, hartel_99, huppert_80}. At the current front a bulbous head may develop and become taller than the trailing fluid. A shear instability can manifest at the density-interface (hereafter interface) between the two fluids, \citep{turner_73}, and this leads to the formation of Kelvin-Helmholtz billows that enhance mixing.

The lock-exchange has been the subject of many theoretical, experimental and numerical studies, and the front speed (or Froude number) is commonly calculated, making it an excellent diagnostic for verification of \fluidity\, e.g. \citep{benjamin_68, klemp_94, hartel_00}. Furthermore the same physical processes are that are encountered in gravity currents over a range of scales are incorporated. Therefore, the lock-exchange presents a tractable means of studying the processes involved and contributes to our understanding of real-world flows, such as sediment-laden density currents and oceanic overflows, and their impact.

In this example, \fluidity\ is used to simulate a lock-exchange and the following functionality is demonstrated:

\begin{itemize}
\item 2D flow
\item Non-hydrostatic flow
\item Incompressible flow
\item Boussinesq flow
\item Mesh adaptivity

\end{itemize}

\subsection{Configuration}
\label{sect:lock_exchange_configuration}

The domain and physical parameters are set up after \cite{fringer_06} and \cite{ozgokmen_07}, table \ref{tab:le_physical_parameters}. The domain is a 2D rectangular box, $0\leq x \leq 0.8\,$m, $0 \leq z \leq 0.1\,$m. Initially, dense, cold water of $T = -0.5\,^\circ$C fills one half of the domain, $x<0.4\,$m, and light, warm water of $T = 0.5\,^\circ$C fills the other half, $x\geq0.4\,$m, figure \ref{fig:lock_exchange}. At $t=0\,$s, ${\bf u} = {\bf 0}$ everywhere. A no-slip boundary condition is applied along the bottom of the domain, ${\bf u = 0}$ at $z=0$, and a free-slip boundary condition is applied to the top of the domain and the side walls, ${\bf u\cdot\bf n} = 0$ at $z = 0.1\,$m and $x = 0.0,\, 0.8\,$m.

The basic choices for the numerical set-up are outlined in table \ref{tab:le_numerical_configuration}. They comprise a set of standard options for a buoyancy-driven flow such as the lock-exchange. The mesh is adapted to both the velocity and the temperature fields.  

\begin{table}[th]
\centering
\begin{tabular}[h]{l   c  c }  \hline
gravitational acceleration (ms$^{-2}$)			& $g$ 								& $10$  \\
kinematic viscosity (m$^2$s$^{-1}$)			& $\nu$								& $10^{-6}$  \\
thermal diffusivity (m$^2$s$^{-1}$)					& $\kappa$ 					& 0  \\
thermal expansion coefficient ($^\circ$C$^{-1}$)	& $\alpha$ 							& $10^{-3}$ \\ 
domain height (m)					& $H = 2h$							& $0.1$ \\ %\hline \hline
reduced gravity	(ms$^{-2}$)				& $g' = g\frac{\rho_1 - \rho_2}{\rho_0} = -g\alpha (T_1 - T_2)$	& $10^{-2}$  \\
buoyancy velocity					& $u_b = \sqrt{g'H}$ 						& $\sqrt{10^{-3}}$ \\
Grashof number						& $Gr = \left( \frac{h \sqrt{g'h}}{\nu} \right)^2$ 		& $1.25 \times 10^{6}$\\ \hline
\end{tabular}
\caption{Physical parameters for the lock-exchange set-up.}
\label{tab:le_physical_parameters}
\end{table}

\begin{table}[th]
\centering
\begin{tabular}[h]{lll}  \hline
Numerical component                           & Configuration                   & Section\\ \hline
Geometry                                      & from triangle file, via gmsh    & \ref{sect:triangle_format}, \ref{sect:meshing_tools_non_fluidity_gmsh}\\
Time step                                     & 0.025\,s                        & \\
Time loop                                     & 2 non-linear Picard iterations  & \ref{sect:ND_time_loop} \\
Equation of state                             & linear                          & \ref{sect:equation_of_state} \\
Momentum and Pressure spatial discretisation  & \Poo                            & \ref{Sect:ND_momentum_equation}, \ref{Sect:ND_pressure_equation}, \ref{sect:configuring_fluidity_continuous_galerkin}, \ref{sect:configuring_fluidity_pressure} \\
Momentum temporal discretisation              & $\theta = 0.5$ (Crank-Nicholson)  & \ref{sect:ND_time_theta_scheme}, \ref{sect:configuring_fluidity_temporal_discretisation} \\
                                              & non-linear relaxation term $=0.5$ & \ref{sect:relax} \\
Temperature advection                         & advection-diffusion equation    & \ref{Sect:MP-AdvdifEqn} \\
                                              & control volumes                 & \ref{Sect:CVs} \\
Temperature temporal discretisation           & $\theta = 0.5$ (Crank-Nicholson)& \ref{sect:ND_time_theta_scheme}, \ref{sect:configuring_fluidity_temporal_discretisation} \\
                                              & advection iterations, limit theta & ??? \\\hline
\end{tabular}
\caption{Numerical configuration for the lock-exchange}
\label{tab:le_numerical_configuration}
\end{table}

\subsection{Results} 
\label{sect:lock_exchange_results}

The expected dynamics of a lock-exchange flow are observed, figure \ref{fig:lock_exchange}: two gravity currents propagate in opposite directions with the foremost point of the no-slip front raised above the lower boundary. Kelvin-Helmholtz billows form at the interface leading to mixing of the two fluids which would not be observed with a hydrostatic formulation. The mesh adapts well, increasing the resolution around the interface and Kelvin-Helmholtz billows with anisotropic elements. Similar images can be generated by visualising the vtu files using Paraview or Mayavi2, see the \href{http://amcg-www.ese.ic.ac.uk/}{AMCG website}\ for more information.

Two diagnostics are considered. First the Froude number, $Fr = U/u_b$, is the ratio of speed, $U$, to the the buoyancy velocity, table \ref{tab:le_physical_parameters} is calculated. The speed with which the no-slip and free-slip fronts propagate along the domain, $U_{ns}$ and $U_{fs}$, are calculated from the model output and are used to give the corresponding no-slip and free-slip Froude numbers, $Fr_{ns}$ and $Fr_{fs}$, figure \ref{fig:lock_exchange_Fr}. For this basic set up the values are comparable to but generally smaller than previously published values, table \ref{tab:lock_exchange_Fr_comparison}. Spatially varying the horizontal velocity adaptivity settings as suggested in the Exercises, section \ref{sect:le_exercises}, can increase the Froude number, showing good agreement between the \fluidity\ values and previously published values. The second diagnostic, the ``domain fraction'', gives an illustration of the mixing taking place, \ref{fig:lock_exchange_mixing}. The temperature ($T$) range is divided into three classes: $T<-0.4$, $-0.4<T<0.4$ and $0.4<T$. For each class the fraction of fluid in the domain that has temperature within the given range gives the domain fraction. The class with range $-0.4<T<0.4$ is considered representative of the mixed fluid. As the simulation progresses, the amount of mixed fluid increases, more rapidly at first as Kelvin-Helmholtz billows form and just after the fronts hit the end wall. Later, as the dynamics become less turbulent and the fluid oscillates back and forth in the tank the mixing is less vigorous. The increase in mixed fluid is compensated for by a decrease in non-mixed fluid.

\begin{figure}[ht]
  \centering
  % For some reason tex4ht doesn't like these images.
  \onlypdf{
  \subfigure[$t = 0\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_0_T}}
  \subfigure[$t = 0\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_0_mesh_nice}} \\
  \subfigure[$t = 12.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_10_T}}
  \subfigure[$t = 12.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_10_mesh}} \\
  \subfigure[$t = 37.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_30_T}}
  \subfigure[$t =
  37.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_30_mesh}}
  }
  \caption{Lock-exchange temperature distribution (colour) with meshes, over time ($t$)}
  \label{fig:lock_exchange}
\end{figure}

\begin{figure}[ht]
  \centering
  \fig[width=1.0\textwidth, clip=True, trim = 32.5mm 20mm 32.5mm 15mm]{./examples_images/lock_exchange/front_speed}
  \caption{Distance along the domain ($X$) and Froude number ($Fr$) for the no-slip and free-slip fronts in the lock-exchange.}
  \label{fig:lock_exchange_Fr}
\end{figure}

\begin{table}[th]
\centering
\begin{tabular}{|c|l|c|c|c|} \hline
& Reference 					& $Gr$ 			& no-slip $Fr$ 	& free-slip $Fr$	\\ [0.1ex] \hline
& Adaptive 1 (\fluidity)			& $1.25 \times 10^{6}$	& 0.381	        & 0.457			\\ [0.1ex]
& Adaptive 2 (\fluidity)			& $1.25 \times 10^{6}$	& 0.417	        & 0.482			\\ [0.1ex]\hline
\multirow{4}{*}{\begin{sideways}\parbox{18mm}{numerical}\end{sideways}} &\cite{hartel_00}, 2D	& $1.25 \times 10^{6}$	& 0.406		& 0.477			\\ [0.1ex]
& \cite{fringer_06}, 2D		                & $1.25 \times 10^{6}$ 	& 0.397 	& 0.428 		\\ [0.1ex]
& \cite{cantero_07}, 3D	                        & $1.50 \times 10^{6}$	& 0.407		& 			\\ [0.1ex]
& \cite{elias_08}, 3D		                & $1.50 \times 10^{6}$	& 0.409		&			\\ [0.1ex] \hline
\multirow{4}{*}{\begin{sideways}\parbox{19.5mm}{experiment}\end{sideways}}& \cite{keulegan_57} & $7.16 \times 10^{7}$ & 0.440&	\\ [0.1ex]
& \cite{britter_78} 	                        & $6.3 \times 10^{7}$	& 	 	& 0.518		 	\\ [0.1ex]
& \cite{simpson_79}                             & $4.8 \times 10^{6}$   & 0.432	        &	 		\\ [0.1ex]
& \cite{shin_04}			        & $>5 \times 10^5$	& 0.464		&			\\ [0.1ex] \hline
\multirow{2}{*}{\begin{sideways}\parbox{11.5mm}{theory}\end{sideways}} & \cite{benjamin_68} &	inviscid	&  & 0.5	\\[0.1ex]
& \cite{klemp_94} 		                &	inviscid	& 		& 0.527			\\ [1.0ex]\hline
\end{tabular}
\caption{Comparison of \fluidity\ Froude numbers with previously published results. \newline Adaptive 1: basic set up as described \newline Adaptive 2: setup with spatially varying interpolation error for the horizontal velocity as suggested in the Exercises, section \ref{sect:le_exercises}.}
\label{tab:lock_exchange_Fr_comparison}
\end{table} 

\begin{figure}[ht]
  \centering
  \fig[width=0.6\textwidth, clip=True, trim = 32.5mm 20mm 32.5mm 20mm]{./examples_images/lock_exchange/mixing}
  \caption{Time evolution of fraction of domain that contains fluid in three temperature classes.}
  \label{fig:lock_exchange_mixing}
\end{figure}

\subsection{Exercises}
\label{sect:le_exercises}
To explore the functionality of \fluidity, the following variations on this example would be constructive learning exercises:
\begin{itemize}
\item Run from the checkpoint to look at the mixing later in the simulation, section \ref{sect:configuring_fluidity_checkpointing}
\item In the adaptivity settings include a spatially varying interpolation error for the horizontal velocity and compare Froude numbers (python can be found in the comment section for the velocity adaptivity settings), section \ref{sect:configuring_fluidity_error_metric}
\item Change the frequency of the adapt, section \ref{sect:configuring_fluidity_adaptivity_options}
\item Turn on metric advection, section \ref{sect:configuring_fluidity_metric_advection}
\item Change the diffusivity and viscosity values
\item Run with a fixed mesh (note this will require making a new input mesh), sections \ref{sect:triangle_format}, \ref{sect:meshing_tools_non_fluidity_gmsh}
\item Try adding some detectors to visualise the particle trajectories, section \ref{sect:diagnostics_detectors} or \newline cf. \option{\ldots/fluidity-longtests/lock\_exchange\_2d\_}
\end{itemize}
Note to compare Froude numbers and mixing between different runs the remember to copy the images otherwise they will be overwritten.

\subsection{Lagrangian trajectories illustrated in the lock-exchange}
\label{sect:lock_exchange_Lagrangian_trajectories}
Lagrangian trajectories from moving or Lagrangian detectors have been calculated in a lock-exchange test case. The test case can be found in \option{\ldots/fluidity-longtests/lock\_exchange\_2d\_} \\ \option{Lagrangian\_paths}. The domain of the 2D rectangular box is in this case $-5\leq x \leq 5\,$m and $-0.5 \leq z \leq 0.5\,$m. In this case a free-slip boundary condition is applied to the top and bottom of the domain ($w = 0$) and also to the end walls ($u = 0$).
Also, the denser (cooler) fluid is now on the right ($0<x \leq 5\,$m) as opposed to the previous example. This can be observed in Figure \ref{fig:lock_exchange_Lag_traject_Temp_mesh}, where the Temperature distribution and the adapted mesh at two different times are displayed.

Figure \ref{fig:lock_exchange_Lagrangian_trajectories} illustrates the Lagrangian trajectories from three particles initially at three different positions in the interface between the fluids.

\begin{figure}[ht]
  \centering
\begin{minipage}{1\linewidth}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Temp_at_1000vtu}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Mesh_adaptivity_at_1000vtu}
\end{minipage}\\
\begin{minipage}{1\linewidth}
\vspace{0.5cm}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Temp_at_1500vtu}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Mesh_adaptivity_at_1500vtu}
\end{minipage}
  \caption{Lock-exchange temperature distribution and meshes at two different instants in time}
  \label{fig:lock_exchange_Lag_traject_Temp_mesh}
\end{figure}

%\vspace*{-1.8cm}

\begin{figure}[ht]
  \centering
\vspace*{-2.6cm}
  \fig[width=0.6\textwidth]{./examples_images/lock_exchange/Lag_paths_lock_exchange_labels}
\vspace{-2.8cm}
  \caption{An example of Lagrangian trajectories in a lock-exchange. Note that in this test case the denser fluid is on the right as opposed to the lock-exchange test case from Figure \ref{fig:lock_exchange} where the denser fluid is on the left. Here, the denser fluid flows below the lighter from right to left. The particles presented at initially three different positions in the interface between the fluids, follow the vortical structures resulting from the shear instability in the interface.}
  \label{fig:lock_exchange_Lagrangian_trajectories}
\end{figure}

Note that the results indicated in Figures \ref{fig:lock_exchange_Lag_traject_Temp_mesh} and \ref{fig:lock_exchange_Lagrangian_trajectories} have been calculated with a time step $\Delta t$=0.1 s. This is a large time step for this simulation but has been chosen for this test case since it runs faster and illustrates the Lagrangian trajectories. However, when accuracy is sought, it is recommended to use a smaller time step, i.e., $\Delta t$=0.01 s.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------ADVECTION--------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Advection example}
\label{sect:advection}

In this problem the advection of a passive tracer by a prescribed velocity field is considered.
The problem is taken from \cite{hanert2004}, see also \cite{hecht2000}. The problem domain
is an $L\times L$ box with $L=\km[1000]$. The prescribed velocity field is obtained from the
Stommel streamfunction
\begin{equation*}
\Psi(x,y) = \frac{FL}{\pi \gamma \rho H}\sin\left(\frac{\pi y}{L}\right)\left(p e^{Ax} + q e^{Bx}-1\right),
\end{equation*}
\begin{equation*}
p = 1-q = \frac{1-e^{LB}}{e^{LA}-e^{LB}},\quad
A=-\frac{\alpha}{2}+\sqrt{\frac{\alpha^2}{4}+\left(\frac{\pi}{L}\right)^2},\quad
B=-\frac{\alpha}{2}-\sqrt{\frac{\alpha^2}{4}+\left(\frac{\pi}{L}\right)^2}.
\end{equation*}
where the magnitude of the wind stress is $F = \Nmm[0.1]$, the frictional
coefficient is $\gamma = \invs[10^{-6}]$, the density is $\rho = \kgmmm[1000]$, the
fluid depth is $H=\m[200]$ and $\alpha = \beta/\gamma$ is chosen such that $\alpha^{-1}=\km[100]$.

The initial condition for the passive tracer is given by the Gaussian function:
\begin{equation}
s_0(x,y) = 10\exp\left({-\frac{(x-2L/3)^2 + (y-L/3)^2}{(L/50)^2}}\right).
\label{eqn:guassian}
\end{equation}

The reference solution at time $t$ used to compare numerical results against
is constructed as follows.  The ordinary differential equation (ODE) system
$\dot{\pmb{x}}=\pmb{u}$ is integrated repeatedly from time $t$ to $0$, with
initial condition given by the node positions on the mesh at time $t$.  This
yields the departure point, at time $0$, for each node on the current mesh.
The reference solution is found by evaluating the function
(\ref{eqn:guassian}) at these departure points.  Integration of the system
of ODEs for departure points is computed using the SciPy (Scientific tools
for Python) odeint function, which itself is a wrapper around the LSODE ODE
solver \cite{hindmarsh1983}.

In this example equation (\ref{PDE}) with $\kappa = 0$ is solved for the evolution of the tracer field until a
time of $\s[10^7]$ . The time step is taken to be $\s[1000]$ and $\theta=0.5$.
For the adaptive simulations the mesh is adapted every 20 time steps.
Linear interpolation is again used to transfer the solution field between adapted meshes.

\subsubsection{Results}

Figure \ref{fig:stommel_fix} shows results at four time levels for three
different fixed uniform isotropic mesh resolutions. Figures \ref{fig:stommel3}
and \ref{fig:stommel4} show the corresponding results computed with anisotropic
adapted meshes. Figure \ref{fig:stommel4} shows the results using a lower interpolation
error in the construction of the metric and hence has higher resolution. While there
are some qualitative differences between the results from the highest resolution fixed
mesh results in figure \ref{fig:stommel_fix} with those shown in figure \ref{fig:stommel3},
the differences are minimal compared with figure \ref{fig:stommel4}.

For quantitative comparison between the fixed and adapted mesh results figures \ref{fig:Stommel2Errors}
and \ref{fig:Stommel2_L2_Errors} show the time evolution of the maximum value of the advected field,
for which the analytical solution preserves a value of ten,
and the relative \Ltwo norm \cite{hanert2004} of the difference between the numerical and the reference
solution whose construction is described above. The results from figures \ref{fig:stommel3} and \ref{fig:stommel4}
correspond to the lowest and middle resolution adaptive results in figures \ref{fig:Stommel2Errors} and
\ref{fig:Stommel2_L2_Errors}.


%\begin{figure}
%\centering
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_2}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_10}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_12}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_20}}
%\caption{Adapted mesh (top) and isolines of the tracer field
%in the advection through a Stommel gyre simulation (section 3(b)) at
%(left to right) time levels $1,5,6,10\times 10^6$s. This simulation uses a tighter tolerance in the definition
%of the metric compared to those results presented in figure \ref{fig:stommel3}.
%At each time level there are 2397, 10433, 9893, 4303 nodes respectively which were obtained from a metric (\ref{eqn:ME}) with $\epsilon_u=0.05$.
%The maximum value for the tracer field at each time level is 9.76, 9.50, 9.46, 9.35 respectively.}
%\label{fig:stommel4}
%\end{figure}

The fixed mesh results can be seen to diffuse more than the adaptive mesh
when the tracer is close to the western boundary current, even for the highest resolution mesh
where over $150\,000$ nodes are used. The adaptive mesh results lose
some of the initial magnitude in the tracer field at the start of the simulation due to the resolution of the
initial mesh used before the first adapt, i.e. in the first 20 time steps.
Improved results would be obtained if a higher resolution initial mesh was used,
or if the mesh was adapted at the start of the simulation to better represent the initial condition.
At the end of the simulated time the magnitude of the tracer field for the highest resolution mesh ($157\, 045$ nodes)
can be seen to have dropped to below that of the middle adaptive resolution simulation, which corresponds to
figure \ref{fig:stommel4}, i.e. ($4303$ nodes).

%\begin{figure}
%\centering
%\fig[width=12.5cm,clip]{images/Stommel2Errors}
%\caption{Left: Maximum value of the advected field in the fixed isotropic mesh simulations for problem 3(b).
%The upper, middle and lower lines correspond to the fine
%($157\,045$ nodes), medium ($39\,490$ nodes) and coarse ($10\,056$ nodes)
%scale meshes respectively. The maximum value for the medium and coarse resolution at the end of the simulated time
%are 7.94 and 5.67 respectively. Middle: Maximum value of the advected field in the adaptive mesh simulations;
%the lines (top-bottom) correspond to the five interpolation errors used, the
%variation in time of the number of nodes in these simulations is shown in the panel on the right
%with the upper line corresponding to the lowest interpolation error and the upper line in the middle panel,
%and so on. The maximum number of nodes used in the highest resolution adaptive simulation is
%$41\,778$ and occurs at time $5.3\times 10^6$s.}
%\label{fig:Stommel2Errors}
%\end{figure}

In the adaptive mesh results the number of nodes used increases as the tracer field
enters the western boundary current, this is due to the gradients (and curvatures)
increasing as the tracer field is stretched, the number of nodes then decreases as
the tracer field leaves the western boundary current region and the simulation becomes
less computationally demanding. This demonstrates an important capability of adaptive methods,
the ability to remove resolution where it is no longer required.

The loss in magnitude of the tracer can be seen to be most significant while the
tracer field is in the western boundary current region, this is due to the resolution
being comparably lower as the scale of the advected field reduces, hence numerical
diffusion increases. The effect is lower in the adaptive simulations as the mesh resolution is increased
more significantly in this area at this time.

These results are consistent with the error plots in figure \ref{fig:Stommel2_L2_Errors}.
The errors grow more rapidly when the tracer is being advected through the western boundary current.
This is the case for both fixed and adaptive mesh results, where again the highest fixed resolution
mesh results match most closely those obtained from the middle adaptive mesh simulation (figure \ref{fig:stommel4}),
i.e. using over an order of magnitude fewer nodes.

With both the \Ltwo error and the maximum value of the advected field, the fixed mesh results
actually tend to perform better before and after the tracer is passing through the western boundary current.
For example, as can be seen from the growth rates of the plotted lines in figures \ref{fig:Stommel2Errors} and \ref{fig:Stommel2_L2_Errors}
 at the start and end of the simulation. The adaptive mesh runs perform better, in terms of
growth of errors, when the problem becomes more challenging numerically and the model is able to automatically use
increased anisotropic resolution close to the western boundary.

%\begin{figure}
%\centering
%\fig[width=8cm,clip]{images/Stommel2_L2_Errors}
%\caption{Relative \Ltwo norm of error between the numerical and reference solutions (for problem 3(b))
%plotted against time using left: fixed uniform mesh and right: adapted meshes.
%Each line, with decreasing error, corresponds to a mesh with increased number of nodes. The details
%of the meshes are the same as those described in the caption to figure \ref{fig:Stommel2Errors}.}
%\label{fig:Stommel2_L2_Errors}
%\end{figure}

The computational cost of the metric construction and the mesh optimisation procedure are approximately
equal for this problem. When combined they correspond approximately to the time taken for the model
to complete 3 time steps for this problem. Hence, the computational overhead of adaptivity can be estimated
at 15 percent of total run time. The middle adaptive mesh simulation (corresponding to figure \ref{fig:stommel4})
uses a maximum number of nodes which corresponds closely to the number used by the coarsest fixed mesh simulation.
However, the adaptive simulation runs approximately 60 percent faster than the fixed mesh simulation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------LID DRIVEN CAVITY---------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Lid-driven cavity}
\label{sect:lid_driven_cavity}

\subsection{Overview}
The lid-driven cavity is a problem that is often used as part of the verification procedure for
CFD codes. The geometry and boundary conditions are simple to prescribe and in
two dimensions there are a number of highly accurate numerical benchmark solutions
available for a wide range of Reynolds numbers \citep{botella1998,bruneau2006,erturk2005}. Here the two-dimensional
problem at a Reynolds number of 1000 is given as an example.

\subsection{Configuration}
The unsteady momentum equations with nonlinear advection and viscosity
terms are solved in a unit square in the $x$ and $y$
directions along with the continuity equations expressing incompressibility. 
No-slip velocity boundary conditions are imposed on boundaries $x=0,1$ and $y=0$, and the prescribed velocity $u=1$, $v=0$ are
set on the boundary $y=1$ (the `lid'). Here the problem is initialised with a zero velocity field and the solution
allowed to converge to steady state via time-stepping.
A subset of the benchmark data available from the literature is then used to
test for numerical convergence. Here this involves the calculation of the kinetic energy
\begin{equation}
\int_{\Omega} (u^2+v^2) \; d\Omega,
\end{equation}
which is compared against the value 0.044503 taken from \cite{bruneau2006}. In addition, the x-component of velocity and
pressure are evaluated at 17 points along the line $x=0.5$ and compared against the data from \cite{botella1998}.



%\begin{table}
%\centering
%\caption{Results from Botella along the $x=0.5$ line ... need to double check these values}
%\smallskip
%\renewcommand{\tabcolsep}{1cm}
%\renewcommand{\arraystretch}{1.0}
%\begin{tabular}{@{\hspace{0.3cm}}c@{\hspace{0.3cm}} @{\hspace{0.3cm}}c@{\hspace{0.3cm}} @{\hspace{0.3cm}}c@{\hspace{0.3cm}}}
%\hline
%$y$     \ &     $u$          \  &  $p$  \\
%
%\hline
%
%0.0000 \   &  0.0000000 \   &  0.110591   \\
%0.0547 \   & -0.1812881 \   &  0.109689   \\
%0.0625 \   & -0.2023300 \   &  0.109200   \\
%0.0703 \   & -0.2228955 \   &  0.108566   \\
%0.1016 \   & -0.3004561 \   &  0.104187   \\
%0.1719 \   & -0.3885691 \   &  0.081925   \\
%0.2813 \   & -0.2803696 \   &  0.040377   \\
%0.4531 \   & -0.1081999 \   &  0.004434   \\
%0.5000 \   & -0.0620561 \   &  0.000000   \\
%0.6172 \   &  0.0570178 \   & -0.000827   \\
%0.7344 \   &  0.1886747 \   &  0.012122   \\
%0.8516 \   &  0.3372212 \   &  0.034910   \\
%0.9531 \   &  0.4723329 \   &  0.050329   \\
%0.9609 \   &  0.5169277 \   &  0.050949   \\
%0.9688 \   &  0.5808359 \   &  0.051514   \\
%0.9766 \   &  0.6644227 \   &  0.052009   \\
%1.0000 \   &  1.0000000 \   &  0.052987   \\
%\hline
%\end{tabular}
%\label{Tab:drivcav}
%\end{table}

\subsection{Results}
Plots of the solutions and benchmark data are given in figure \ref{fig:drivcav1}. Also shown is a plot
of the error convergence with mesh spacing. Here a regular triangular mesh is used with progressive uniform
refinement in the $x,y$ plane. Second order spatial convergence can clearly be seen for the three quantities compared.
Adaptive refinement is not particularly advantageous for the problem at this reasonably low Reynolds number, but
yields significant improvements in efficiency at higher Reynolds number where boundary layers and
recirculating eddies are more dynamic, anisotropic and smaller in size compared to the entire domain.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------BACKWARD FACING STEP------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Backward facing step}
\label{sect:backward_facing_step}


\subsection{Overview}
The backward-facing step is a classical CFD example and one of the most frequently selected
problems for simulating the separation and reattachment of turbulent flows.
It is also often used as a test problem for validating and benchmarking numerical codes.
At high Reynolds numbers and in three spatial dimensions the problem has substantial
computing requirements, this makes it an ideal HPC benchmark problem for use here.
In the context of ocean modelling flow separation is important
in large scale western boundary currents such as the Gulf Stream.
<commit_msg>Update on the text for the driven cavity problem. Could do with a proof read from someone else.<commit_after>\chapter{Examples}
\label{chap:examples}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%--------------------WRITING A TEST CASE--------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Writing a test case} 

\subsection{Options file validation}

Since the grammar of \fluidity\ options files is defined in terms of a RelaxNG
schema, it is possible to check whether a given options file is valid. The
medium test ``diamond\_validation'' validates all test problem options files
with the .flml file extension. Hence if the options schema is updated in such
a way that an options file becomes invalid, the diamond\_validation test will
flag an error. Test problem XML files are also tested for validation against a
RelaxNG schema defining their format.

Some validation failures, such as the addition of new compulsory options, can
be repaired automatically using Diamond. The \lstinline[language = Bash]+update_options+
script can be used to attempt an automatic update of all test options files.

\subsection{The method of manufactured solutions}\label{sect:mms}
\index{method of manufactured solutions}
Problems with analytical solutions provide a powerful mechanism for testing the implementation
of numerical algorithms (code verification). The problem can be run on a series of meshes
with increasing resolution (\eg $\Delta x \rightarrow \Delta x/2$), or decreasing time step,
and the exact error between numerical and analytical solution computed in any norm. The theoretical
order of convergence of the algorithm can then be confirmed for the implentation code. Any significant
deviation from the theoretical order of convergence can be identified and is often due to coding bugs
or because the treatment of boundary conditions is reducing accuracy.

However, analytical solutions are not generally available, especially for complex equations sets
and for problems with physical realism. These are of course the problems that we are primarily 
interested in tackling numerically.
The method of manufactured solutions (MMS) takes the opoosite approach to deriving a solution to 
a given equation set: it starts from a given solution and derives the equations that admit this
solution. This achieved by substituting the solution into the eqution set and evaluating the residual
symbolically. This residual is then taken as a prescribed source term in the equation set being considered
and the identical procedure used for verification with an analytical solution used. The flexibility of the MMS
approach means that individual terms and discretisations can be rigorously tested through a suite of MMS problems. 

This procedure is especially useful in \fluidity\ as it is very easy to prescibe source terms in the equations
being solved. When substituting in a relatively simple analytical form for a solution into the governing
equations one will typically obtain a very complex expression for the residual or source term. While this
can result in a problem for many codes, in \fluidity\ one is able to prescribe anaytical source terms easily
within Diamond's Python function interface. The complex source term can be obtained using a Python symbolic algebra
library such as SAGE and then simply cut-and-paste into the Diamond file.

Some background material on the method of manufactured solutions, and guidelines for selecting good solutions,
is given in \cite{salari2000}. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------THE LOCK-EXCHANGE---------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The lock-exchange}
\label{sect:lock_exchange}
\subsection{Overview}
\label{sect:lock_exchange_overview}

The lock-exchange is a classic laboratory-scale fluid dynamics problem, \citep{fannelop_94, huppert_06, simpson_87}. A flat-bottomed tank is separated into two portions by a vertical barrier. One portion, the `lock', is filled with the source fluid. This is of different density to the ambient fluid which fills the other portion. As the barrier is removed, the denser fluid collapses under the lighter. Two gravity currents form and propagate in opposite directions, one above the other, along the tank. After an initial acceleration, the gravity current fronts travel at a constant speed until the end walls exert an influence or viscous forces begin to dominate, \citep{cantero_07, hartel_99, huppert_80}. At the current front a bulbous head may develop and become taller than the trailing fluid. A shear instability can manifest at the density-interface (hereafter interface) between the two fluids, \citep{turner_73}, and this leads to the formation of Kelvin-Helmholtz billows that enhance mixing.

The lock-exchange has been the subject of many theoretical, experimental and numerical studies, and the front speed (or Froude number) is commonly calculated, making it an excellent diagnostic for verification of \fluidity\, e.g. \citep{benjamin_68, klemp_94, hartel_00}. Furthermore the same physical processes are that are encountered in gravity currents over a range of scales are incorporated. Therefore, the lock-exchange presents a tractable means of studying the processes involved and contributes to our understanding of real-world flows, such as sediment-laden density currents and oceanic overflows, and their impact.

In this example, \fluidity\ is used to simulate a lock-exchange and the following functionality is demonstrated:

\begin{itemize}
\item 2D flow
\item Non-hydrostatic flow
\item Incompressible flow
\item Boussinesq flow
\item Mesh adaptivity

\end{itemize}

\subsection{Configuration}
\label{sect:lock_exchange_configuration}

The domain and physical parameters are set up after \cite{fringer_06} and \cite{ozgokmen_07}, table \ref{tab:le_physical_parameters}. The domain is a 2D rectangular box, $0\leq x \leq 0.8\,$m, $0 \leq z \leq 0.1\,$m. Initially, dense, cold water of $T = -0.5\,^\circ$C fills one half of the domain, $x<0.4\,$m, and light, warm water of $T = 0.5\,^\circ$C fills the other half, $x\geq0.4\,$m, figure \ref{fig:lock_exchange}. At $t=0\,$s, ${\bf u} = {\bf 0}$ everywhere. A no-slip boundary condition is applied along the bottom of the domain, ${\bf u = 0}$ at $z=0$, and a free-slip boundary condition is applied to the top of the domain and the side walls, ${\bf u\cdot\bf n} = 0$ at $z = 0.1\,$m and $x = 0.0,\, 0.8\,$m.

The basic choices for the numerical set-up are outlined in table \ref{tab:le_numerical_configuration}. They comprise a set of standard options for a buoyancy-driven flow such as the lock-exchange. The mesh is adapted to both the velocity and the temperature fields.  

\begin{table}[th]
\centering
\begin{tabular}[h]{l   c  c }  \hline
gravitational acceleration (ms$^{-2}$)			& $g$ 								& $10$  \\
kinematic viscosity (m$^2$s$^{-1}$)			& $\nu$								& $10^{-6}$  \\
thermal diffusivity (m$^2$s$^{-1}$)					& $\kappa$ 					& 0  \\
thermal expansion coefficient ($^\circ$C$^{-1}$)	& $\alpha$ 							& $10^{-3}$ \\ 
domain height (m)					& $H = 2h$							& $0.1$ \\ %\hline \hline
reduced gravity	(ms$^{-2}$)				& $g' = g\frac{\rho_1 - \rho_2}{\rho_0} = -g\alpha (T_1 - T_2)$	& $10^{-2}$  \\
buoyancy velocity					& $u_b = \sqrt{g'H}$ 						& $\sqrt{10^{-3}}$ \\
Grashof number						& $Gr = \left( \frac{h \sqrt{g'h}}{\nu} \right)^2$ 		& $1.25 \times 10^{6}$\\ \hline
\end{tabular}
\caption{Physical parameters for the lock-exchange set-up.}
\label{tab:le_physical_parameters}
\end{table}

\begin{table}[th]
\centering
\begin{tabular}[h]{lll}  \hline
Numerical component                           & Configuration                   & Section\\ \hline
Geometry                                      & from triangle file, via gmsh    & \ref{sect:triangle_format}, \ref{sect:meshing_tools_non_fluidity_gmsh}\\
Time step                                     & 0.025\,s                        & \\
Time loop                                     & 2 non-linear Picard iterations  & \ref{sect:ND_time_loop} \\
Equation of state                             & linear                          & \ref{sect:equation_of_state} \\
Momentum and Pressure spatial discretisation  & \Poo                            & \ref{Sect:ND_momentum_equation}, \ref{Sect:ND_pressure_equation}, \ref{sect:configuring_fluidity_continuous_galerkin}, \ref{sect:configuring_fluidity_pressure} \\
Momentum temporal discretisation              & $\theta = 0.5$ (Crank-Nicholson)  & \ref{sect:ND_time_theta_scheme}, \ref{sect:configuring_fluidity_temporal_discretisation} \\
                                              & non-linear relaxation term $=0.5$ & \ref{sect:relax} \\
Temperature advection                         & advection-diffusion equation    & \ref{Sect:MP-AdvdifEqn} \\
                                              & control volumes                 & \ref{Sect:CVs} \\
Temperature temporal discretisation           & $\theta = 0.5$ (Crank-Nicholson)& \ref{sect:ND_time_theta_scheme}, \ref{sect:configuring_fluidity_temporal_discretisation} \\
                                              & advection iterations, limit theta & ??? \\\hline
\end{tabular}
\caption{Numerical configuration for the lock-exchange}
\label{tab:le_numerical_configuration}
\end{table}

\subsection{Results} 
\label{sect:lock_exchange_results}

The expected dynamics of a lock-exchange flow are observed, figure \ref{fig:lock_exchange}: two gravity currents propagate in opposite directions with the foremost point of the no-slip front raised above the lower boundary. Kelvin-Helmholtz billows form at the interface leading to mixing of the two fluids which would not be observed with a hydrostatic formulation. The mesh adapts well, increasing the resolution around the interface and Kelvin-Helmholtz billows with anisotropic elements. Similar images can be generated by visualising the vtu files using Paraview or Mayavi2, see the \href{http://amcg-www.ese.ic.ac.uk/}{AMCG website}\ for more information.

Two diagnostics are considered. First the Froude number, $Fr = U/u_b$, is the ratio of speed, $U$, to the the buoyancy velocity, table \ref{tab:le_physical_parameters} is calculated. The speed with which the no-slip and free-slip fronts propagate along the domain, $U_{ns}$ and $U_{fs}$, are calculated from the model output and are used to give the corresponding no-slip and free-slip Froude numbers, $Fr_{ns}$ and $Fr_{fs}$, figure \ref{fig:lock_exchange_Fr}. For this basic set up the values are comparable to but generally smaller than previously published values, table \ref{tab:lock_exchange_Fr_comparison}. Spatially varying the horizontal velocity adaptivity settings as suggested in the Exercises, section \ref{sect:le_exercises}, can increase the Froude number, showing good agreement between the \fluidity\ values and previously published values. The second diagnostic, the ``domain fraction'', gives an illustration of the mixing taking place, \ref{fig:lock_exchange_mixing}. The temperature ($T$) range is divided into three classes: $T<-0.4$, $-0.4<T<0.4$ and $0.4<T$. For each class the fraction of fluid in the domain that has temperature within the given range gives the domain fraction. The class with range $-0.4<T<0.4$ is considered representative of the mixed fluid. As the simulation progresses, the amount of mixed fluid increases, more rapidly at first as Kelvin-Helmholtz billows form and just after the fronts hit the end wall. Later, as the dynamics become less turbulent and the fluid oscillates back and forth in the tank the mixing is less vigorous. The increase in mixed fluid is compensated for by a decrease in non-mixed fluid.

\begin{figure}[ht]
  \centering
  % For some reason tex4ht doesn't like these images.
  \onlypdf{
  \subfigure[$t = 0\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_0_T}}
  \subfigure[$t = 0\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_0_mesh_nice}} \\
  \subfigure[$t = 12.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_10_T}}
  \subfigure[$t = 12.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_10_mesh}} \\
  \subfigure[$t = 37.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_30_T}}
  \subfigure[$t =
  37.475\,$s]{\fig[width=0.45\textwidth]{./examples_images/lock_exchange/le_basic_30_mesh}}
  }
  \caption{Lock-exchange temperature distribution (colour) with meshes, over time ($t$)}
  \label{fig:lock_exchange}
\end{figure}

\begin{figure}[ht]
  \centering
  \fig[width=1.0\textwidth, clip=True, trim = 32.5mm 20mm 32.5mm 15mm]{./examples_images/lock_exchange/front_speed}
  \caption{Distance along the domain ($X$) and Froude number ($Fr$) for the no-slip and free-slip fronts in the lock-exchange.}
  \label{fig:lock_exchange_Fr}
\end{figure}

\begin{table}[th]
\centering
\begin{tabular}{|c|l|c|c|c|} \hline
& Reference 					& $Gr$ 			& no-slip $Fr$ 	& free-slip $Fr$	\\ [0.1ex] \hline
& Adaptive 1 (\fluidity)			& $1.25 \times 10^{6}$	& 0.381	        & 0.457			\\ [0.1ex]
& Adaptive 2 (\fluidity)			& $1.25 \times 10^{6}$	& 0.417	        & 0.482			\\ [0.1ex]\hline
\multirow{4}{*}{\begin{sideways}\parbox{18mm}{numerical}\end{sideways}} &\cite{hartel_00}, 2D	& $1.25 \times 10^{6}$	& 0.406		& 0.477			\\ [0.1ex]
& \cite{fringer_06}, 2D		                & $1.25 \times 10^{6}$ 	& 0.397 	& 0.428 		\\ [0.1ex]
& \cite{cantero_07}, 3D	                        & $1.50 \times 10^{6}$	& 0.407		& 			\\ [0.1ex]
& \cite{elias_08}, 3D		                & $1.50 \times 10^{6}$	& 0.409		&			\\ [0.1ex] \hline
\multirow{4}{*}{\begin{sideways}\parbox{19.5mm}{experiment}\end{sideways}}& \cite{keulegan_57} & $7.16 \times 10^{7}$ & 0.440&	\\ [0.1ex]
& \cite{britter_78} 	                        & $6.3 \times 10^{7}$	& 	 	& 0.518		 	\\ [0.1ex]
& \cite{simpson_79}                             & $4.8 \times 10^{6}$   & 0.432	        &	 		\\ [0.1ex]
& \cite{shin_04}			        & $>5 \times 10^5$	& 0.464		&			\\ [0.1ex] \hline
\multirow{2}{*}{\begin{sideways}\parbox{11.5mm}{theory}\end{sideways}} & \cite{benjamin_68} &	inviscid	&  & 0.5	\\[0.1ex]
& \cite{klemp_94} 		                &	inviscid	& 		& 0.527			\\ [1.0ex]\hline
\end{tabular}
\caption{Comparison of \fluidity\ Froude numbers with previously published results. \newline Adaptive 1: basic set up as described \newline Adaptive 2: setup with spatially varying interpolation error for the horizontal velocity as suggested in the Exercises, section \ref{sect:le_exercises}.}
\label{tab:lock_exchange_Fr_comparison}
\end{table} 

\begin{figure}[ht]
  \centering
  \fig[width=0.6\textwidth, clip=True, trim = 32.5mm 20mm 32.5mm 20mm]{./examples_images/lock_exchange/mixing}
  \caption{Time evolution of fraction of domain that contains fluid in three temperature classes.}
  \label{fig:lock_exchange_mixing}
\end{figure}

\subsection{Exercises}
\label{sect:le_exercises}
To explore the functionality of \fluidity, the following variations on this example would be constructive learning exercises:
\begin{itemize}
\item Run from the checkpoint to look at the mixing later in the simulation, section \ref{sect:configuring_fluidity_checkpointing}
\item In the adaptivity settings include a spatially varying interpolation error for the horizontal velocity and compare Froude numbers (python can be found in the comment section for the velocity adaptivity settings), section \ref{sect:configuring_fluidity_error_metric}
\item Change the frequency of the adapt, section \ref{sect:configuring_fluidity_adaptivity_options}
\item Turn on metric advection, section \ref{sect:configuring_fluidity_metric_advection}
\item Change the diffusivity and viscosity values
\item Run with a fixed mesh (note this will require making a new input mesh), sections \ref{sect:triangle_format}, \ref{sect:meshing_tools_non_fluidity_gmsh}
\item Try adding some detectors to visualise the particle trajectories, section \ref{sect:diagnostics_detectors} or \newline cf. \option{\ldots/fluidity-longtests/lock\_exchange\_2d\_}
\end{itemize}
Note to compare Froude numbers and mixing between different runs the remember to copy the images otherwise they will be overwritten.

\subsection{Lagrangian trajectories illustrated in the lock-exchange}
\label{sect:lock_exchange_Lagrangian_trajectories}
Lagrangian trajectories from moving or Lagrangian detectors have been calculated in a lock-exchange test case. The test case can be found in \option{\ldots/fluidity-longtests/lock\_exchange\_2d\_} \\ \option{Lagrangian\_paths}. The domain of the 2D rectangular box is in this case $-5\leq x \leq 5\,$m and $-0.5 \leq z \leq 0.5\,$m. In this case a free-slip boundary condition is applied to the top and bottom of the domain ($w = 0$) and also to the end walls ($u = 0$).
Also, the denser (cooler) fluid is now on the right ($0<x \leq 5\,$m) as opposed to the previous example. This can be observed in Figure \ref{fig:lock_exchange_Lag_traject_Temp_mesh}, where the Temperature distribution and the adapted mesh at two different times are displayed.

Figure \ref{fig:lock_exchange_Lagrangian_trajectories} illustrates the Lagrangian trajectories from three particles initially at three different positions in the interface between the fluids.

\begin{figure}[ht]
  \centering
\begin{minipage}{1\linewidth}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Temp_at_1000vtu}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Mesh_adaptivity_at_1000vtu}
\end{minipage}\\
\begin{minipage}{1\linewidth}
\vspace{0.5cm}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Temp_at_1500vtu}
  \fig[width=0.48\textwidth]{./examples_images/lock_exchange/Mesh_adaptivity_at_1500vtu}
\end{minipage}
  \caption{Lock-exchange temperature distribution and meshes at two different instants in time}
  \label{fig:lock_exchange_Lag_traject_Temp_mesh}
\end{figure}

%\vspace*{-1.8cm}

\begin{figure}[ht]
  \centering
\vspace*{-2.6cm}
  \fig[width=0.6\textwidth]{./examples_images/lock_exchange/Lag_paths_lock_exchange_labels}
\vspace{-2.8cm}
  \caption{An example of Lagrangian trajectories in a lock-exchange. Note that in this test case the denser fluid is on the right as opposed to the lock-exchange test case from Figure \ref{fig:lock_exchange} where the denser fluid is on the left. Here, the denser fluid flows below the lighter from right to left. The particles presented at initially three different positions in the interface between the fluids, follow the vortical structures resulting from the shear instability in the interface.}
  \label{fig:lock_exchange_Lagrangian_trajectories}
\end{figure}

Note that the results indicated in Figures \ref{fig:lock_exchange_Lag_traject_Temp_mesh} and \ref{fig:lock_exchange_Lagrangian_trajectories} have been calculated with a time step $\Delta t$=0.1 s. This is a large time step for this simulation but has been chosen for this test case since it runs faster and illustrates the Lagrangian trajectories. However, when accuracy is sought, it is recommended to use a smaller time step, i.e., $\Delta t$=0.01 s.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------ADVECTION--------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Advection example}
\label{sect:advection}

In this problem the advection of a passive tracer by a prescribed velocity field is considered.
The problem is taken from \cite{hanert2004}, see also \cite{hecht2000}. The problem domain
is an $L\times L$ box with $L=\km[1000]$. The prescribed velocity field is obtained from the
Stommel streamfunction
\begin{equation*}
\Psi(x,y) = \frac{FL}{\pi \gamma \rho H}\sin\left(\frac{\pi y}{L}\right)\left(p e^{Ax} + q e^{Bx}-1\right),
\end{equation*}
\begin{equation*}
p = 1-q = \frac{1-e^{LB}}{e^{LA}-e^{LB}},\quad
A=-\frac{\alpha}{2}+\sqrt{\frac{\alpha^2}{4}+\left(\frac{\pi}{L}\right)^2},\quad
B=-\frac{\alpha}{2}-\sqrt{\frac{\alpha^2}{4}+\left(\frac{\pi}{L}\right)^2}.
\end{equation*}
where the magnitude of the wind stress is $F = \Nmm[0.1]$, the frictional
coefficient is $\gamma = \invs[10^{-6}]$, the density is $\rho = \kgmmm[1000]$, the
fluid depth is $H=\m[200]$ and $\alpha = \beta/\gamma$ is chosen such that $\alpha^{-1}=\km[100]$.

The initial condition for the passive tracer is given by the Gaussian function:
\begin{equation}
s_0(x,y) = 10\exp\left({-\frac{(x-2L/3)^2 + (y-L/3)^2}{(L/50)^2}}\right).
\label{eqn:guassian}
\end{equation}

The reference solution at time $t$ used to compare numerical results against
is constructed as follows.  The ordinary differential equation (ODE) system
$\dot{\pmb{x}}=\pmb{u}$ is integrated repeatedly from time $t$ to $0$, with
initial condition given by the node positions on the mesh at time $t$.  This
yields the departure point, at time $0$, for each node on the current mesh.
The reference solution is found by evaluating the function
(\ref{eqn:guassian}) at these departure points.  Integration of the system
of ODEs for departure points is computed using the SciPy (Scientific tools
for Python) odeint function, which itself is a wrapper around the LSODE ODE
solver \cite{hindmarsh1983}.

In this example equation (\ref{PDE}) with $\kappa = 0$ is solved for the evolution of the tracer field until a
time of $\s[10^7]$ . The time step is taken to be $\s[1000]$ and $\theta=0.5$.
For the adaptive simulations the mesh is adapted every 20 time steps.
Linear interpolation is again used to transfer the solution field between adapted meshes.

\subsubsection{Results}

Figure \ref{fig:stommel_fix} shows results at four time levels for three
different fixed uniform isotropic mesh resolutions. Figures \ref{fig:stommel3}
and \ref{fig:stommel4} show the corresponding results computed with anisotropic
adapted meshes. Figure \ref{fig:stommel4} shows the results using a lower interpolation
error in the construction of the metric and hence has higher resolution. While there
are some qualitative differences between the results from the highest resolution fixed
mesh results in figure \ref{fig:stommel_fix} with those shown in figure \ref{fig:stommel3},
the differences are minimal compared with figure \ref{fig:stommel4}.

For quantitative comparison between the fixed and adapted mesh results figures \ref{fig:Stommel2Errors}
and \ref{fig:Stommel2_L2_Errors} show the time evolution of the maximum value of the advected field,
for which the analytical solution preserves a value of ten,
and the relative \Ltwo norm \cite{hanert2004} of the difference between the numerical and the reference
solution whose construction is described above. The results from figures \ref{fig:stommel3} and \ref{fig:stommel4}
correspond to the lowest and middle resolution adaptive results in figures \ref{fig:Stommel2Errors} and
\ref{fig:Stommel2_L2_Errors}.


%\begin{figure}
%\centering
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_2}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_10}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_12}}
%\subfigure{
%\fig[width=2.95cm,clip]{images/Adapt_C_20}}
%\caption{Adapted mesh (top) and isolines of the tracer field
%in the advection through a Stommel gyre simulation (section 3(b)) at
%(left to right) time levels $1,5,6,10\times 10^6$s. This simulation uses a tighter tolerance in the definition
%of the metric compared to those results presented in figure \ref{fig:stommel3}.
%At each time level there are 2397, 10433, 9893, 4303 nodes respectively which were obtained from a metric (\ref{eqn:ME}) with $\epsilon_u=0.05$.
%The maximum value for the tracer field at each time level is 9.76, 9.50, 9.46, 9.35 respectively.}
%\label{fig:stommel4}
%\end{figure}

The fixed mesh results can be seen to diffuse more than the adaptive mesh
when the tracer is close to the western boundary current, even for the highest resolution mesh
where over $150\,000$ nodes are used. The adaptive mesh results lose
some of the initial magnitude in the tracer field at the start of the simulation due to the resolution of the
initial mesh used before the first adapt, i.e. in the first 20 time steps.
Improved results would be obtained if a higher resolution initial mesh was used,
or if the mesh was adapted at the start of the simulation to better represent the initial condition.
At the end of the simulated time the magnitude of the tracer field for the highest resolution mesh ($157\, 045$ nodes)
can be seen to have dropped to below that of the middle adaptive resolution simulation, which corresponds to
figure \ref{fig:stommel4}, i.e. ($4303$ nodes).

%\begin{figure}
%\centering
%\fig[width=12.5cm,clip]{images/Stommel2Errors}
%\caption{Left: Maximum value of the advected field in the fixed isotropic mesh simulations for problem 3(b).
%The upper, middle and lower lines correspond to the fine
%($157\,045$ nodes), medium ($39\,490$ nodes) and coarse ($10\,056$ nodes)
%scale meshes respectively. The maximum value for the medium and coarse resolution at the end of the simulated time
%are 7.94 and 5.67 respectively. Middle: Maximum value of the advected field in the adaptive mesh simulations;
%the lines (top-bottom) correspond to the five interpolation errors used, the
%variation in time of the number of nodes in these simulations is shown in the panel on the right
%with the upper line corresponding to the lowest interpolation error and the upper line in the middle panel,
%and so on. The maximum number of nodes used in the highest resolution adaptive simulation is
%$41\,778$ and occurs at time $5.3\times 10^6$s.}
%\label{fig:Stommel2Errors}
%\end{figure}

In the adaptive mesh results the number of nodes used increases as the tracer field
enters the western boundary current, this is due to the gradients (and curvatures)
increasing as the tracer field is stretched, the number of nodes then decreases as
the tracer field leaves the western boundary current region and the simulation becomes
less computationally demanding. This demonstrates an important capability of adaptive methods,
the ability to remove resolution where it is no longer required.

The loss in magnitude of the tracer can be seen to be most significant while the
tracer field is in the western boundary current region, this is due to the resolution
being comparably lower as the scale of the advected field reduces, hence numerical
diffusion increases. The effect is lower in the adaptive simulations as the mesh resolution is increased
more significantly in this area at this time.

These results are consistent with the error plots in figure \ref{fig:Stommel2_L2_Errors}.
The errors grow more rapidly when the tracer is being advected through the western boundary current.
This is the case for both fixed and adaptive mesh results, where again the highest fixed resolution
mesh results match most closely those obtained from the middle adaptive mesh simulation (figure \ref{fig:stommel4}),
i.e. using over an order of magnitude fewer nodes.

With both the \Ltwo error and the maximum value of the advected field, the fixed mesh results
actually tend to perform better before and after the tracer is passing through the western boundary current.
For example, as can be seen from the growth rates of the plotted lines in figures \ref{fig:Stommel2Errors} and \ref{fig:Stommel2_L2_Errors}
 at the start and end of the simulation. The adaptive mesh runs perform better, in terms of
growth of errors, when the problem becomes more challenging numerically and the model is able to automatically use
increased anisotropic resolution close to the western boundary.

%\begin{figure}
%\centering
%\fig[width=8cm,clip]{images/Stommel2_L2_Errors}
%\caption{Relative \Ltwo norm of error between the numerical and reference solutions (for problem 3(b))
%plotted against time using left: fixed uniform mesh and right: adapted meshes.
%Each line, with decreasing error, corresponds to a mesh with increased number of nodes. The details
%of the meshes are the same as those described in the caption to figure \ref{fig:Stommel2Errors}.}
%\label{fig:Stommel2_L2_Errors}
%\end{figure}

The computational cost of the metric construction and the mesh optimisation procedure are approximately
equal for this problem. When combined they correspond approximately to the time taken for the model
to complete 3 time steps for this problem. Hence, the computational overhead of adaptivity can be estimated
at 15 percent of total run time. The middle adaptive mesh simulation (corresponding to figure \ref{fig:stommel4})
uses a maximum number of nodes which corresponds closely to the number used by the coarsest fixed mesh simulation.
However, the adaptive simulation runs approximately 60 percent faster than the fixed mesh simulation.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------LID DRIVEN CAVITY---------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Lid-driven cavity}
\label{sect:lid_driven_cavity}

\subsection{Overview}
The lid-driven cavity is a problem that is often used as part of the verification procedure for
CFD codes. The geometry and boundary conditions are simple to prescribe and in
two dimensions there are a number of highly accurate numerical benchmark solutions
available for a wide range of Reynolds numbers \citep{botella1998,erturk2005,bruneau2006}. 
Here the two-dimensional problem at a Reynolds number of 1000 is given as an example.

\subsection{Configuration}
The problem domain is $(x,y) \in [0,1]^2$ and this is represented by an unstructured triangluar mesh
of approximately uniform resolution. To enable convergence analysis a sequence of meshes are considered
with spacing of $h=1/2^n$, $n=4,5,6,7$. The two-dimensional incompressible Navier-Stokes equations are
considered with unit density.

No-slip velocity boundary conditions are imposed on the boundaries $x=0,1$ and $y=0$, 
and the prescribed velocity $u=1$, $v=0$ are set on the boundary $y=1$ (the `lid'). 

Here the problem is initialised with a zero velocity field and the solution
allowed to converge to steady state via time-stepping. The steady state convergence criterion is
that in the infinity norm the velocity varies by less that $10^{-6}$ between time levels.

The time step is constant between all runs, $\Delta t = 0.05$, and data is dumped to disk every
1.0 time units (20 time steps). The Crank-Nicolson $(\theta=0.5)$ discretisation is used in time. 
A \Poo Galerkin method is used for the discretisation of velocity and pressure, with no stabilisation
of velocity. Based on a domain size of $L=1.0$ and a velocity scale taken from the lid $(U=1.0)$, to achive
a Reynolds number of 1000 viscosity is taken to be isotropic with a value of $\nu=0.001$.

 

\begin{figure}
\centering
\subfigure{
\includegraphics[width=7cm,clip]{examples_images/driven_cavity_2d/driven_cavity_streamfunction.png}}
\subfigure{
\includegraphics[width=7cm,clip]{examples_images/driven_cavity_2d/driven_cavity_vorticity.png}}
\caption{Diagnostic fields from the lid-driven cavity problem at steady state at $1/128$ resolution.
Left: the streamfunction. Right: the vorticity. The contour levels are taken from those given by \cite{botella1998} in their tables
7 and 8.}
\label{fig:driven_cavity1}
\end{figure}


\subsection{Results}
Plots of the streamfunction and vorticity from the $h=1/128$ simulation are shown in figure \ref{fig:driven_cavity1}.
Observe the good qualititative agreement with \citep{botella1998,erturk2005,bruneau2006}. 
For a quantitative comaprison we take advantage of the tabulated benchmark data available in these papers.
Only a subset of these are computed: the $u$ velocity at a series of points along the line $x=0.5$ and the
$v$ velocity at a series of points along the line $y=0.5$ from \citep{erturk2005}, the same quantities at
different points along the same lines as well as the pressure along both the $x=0.5$ and $y=0.5$ lines from 
\citep{botella1998}, and the kinetic energy and the minimum streamfunction value from \citep{bruneau2006}.
The RMS difference in the case of the first six sets of benchmark data are taken with values extracted from 
the numerical solution, and the absolute difference taken in the final two. These eight error values are
defined as error1, \ldots error8 respectively and and plotted for the four mesh resolutions in figure 
\ref{fig:driven_cavity2}.
Second order spatial convergence can clearly be seen.
Adaptive refinement is not particularly advantageous for the problem at this reasonably low Reynolds number, but
yields significant improvements in efficiency at higher Reynolds number where boundary layers and
recirculating eddies are more dynamic, anisotropic and smaller in size compared to the entire domain.

\begin{figure}
\centering
\includegraphics[width=10cm,clip]{examples_images/driven_cavity_2d/driven_cavity_error_plot.png}
\caption{Convergence of the eight error metrics computed for the lid-driven cavity problem with mesh spacing. The eight
metrics are described in the text.}
\label{fig:driven_cavity2}
\end{figure}

\subsection{Exercises}
\begin{enumerate}
\item Examine the way that the $u=1$ lid condition is applied in the flml file. Why has it not been set as simply a constant?
Try changing it to a constant and see what happens to the errors that you achieve. [Hint: some people consider to "regularised" 
lid-driven cavity problem, try finding some papers that discuss this, try updating the boundary condition so it matches the 
regularised problem and compare to benchmark data if you can find it].
\item The references given above include data from other Reynolds numbers. Try updating the problem set-up and the python which computes
the errors for a higher Reynolds number.
\item Try switching on mesh adaptivity. Test adapting based on different metrics, e.g. try weighting $u$, $v$ and $p$
differently and see what meshes you get. Try varying these weights as well as the maximum and minimum allowed element
sizes to see how they affect each other and the mesh that results. Can you get a metric that results in a lower
error for the same number of nodes compared to the fixed mesh [Hint: it may be easier to achieve this at higher Reynolds numbers].
\end{enumerate}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------BACKWARD FACING STEP------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Backward facing step}
\label{sect:backward_facing_step}


\subsection{Overview}
The backward-facing step is a classical CFD example and one of the most frequently selected
problems for simulating the separation and reattachment of turbulent flows.
It is also often used as a test problem for validating and benchmarking numerical codes.
At high Reynolds numbers and in three spatial dimensions the problem has substantial
computing requirements, this makes it an ideal HPC benchmark problem for use here.
In the context of ocean modelling flow separation is important
in large scale western boundary currents such as the Gulf Stream.

{'fraction_non_alphanumeric': 0.11256991806989616, 'fraction_numerical': 0.02669665415565259, 'mean_word_length': 4.1933257489571485, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1414591', 'n_tokens_mistral': 20767, 'n_tokens_neox': 18754, 'n_words': 8119}