\subsection{The convection-diffusion problem}
\label{sec-diffusion-convection}%
%
% TODO:
%  - reduce code size 
%  - phi1 = convect (u, -dt, phi0)
%  - graphic output & anims
%
\pbindex{convection-diffusion}%
\cindex{method!Euler implicit scheme}%

\subsubsection*{Formulation}
% -----------------------
\label{pb-transport-unsteady}%
Let $T>0$  and $\nu>0$.
The convection-diffusion problem writes:
\begin{quote}
	$(P)$: {\em find $\phi$, defined in $\Omega \times ]0,T[$, such that}
\end{quote}
\begin{eqnarray*}
	\Frac{\partial \phi}{\partial t} + {\bf u}.\nabla \phi - \nu \Delta \phi 
		+ \sigma \phi &=& 0
		\ \mbox{ in } \Omega \times ]0,T[
		\\
	\phi(0) &=& \phi_0
		\ \mbox{ in } \Omega
		\\
	\phi(t) &=& \phi_\Gamma(t)
		\ \mbox{ on } \partial\Omega \times ]0,T[ 
\end{eqnarray*}
where ${\bf u}$, $\sigma \geq 0$, $\phi_0$ and $\phi_\Gamma$ being known.
Note the additional ${\bf u}.\nabla$ operator.

\subsubsection*{Time approximation}
% ------------------------------
\label{characteristic-method}%
\cindex{method!characteristic}%

This problem is approximated by the following 
first-order implicit Euler scheme:
\[
	\Frac{\phi^{n+1}-\phi^n \circ X^n}{\Delta t} - \nu \Delta \phi^{n+1} + \sigma \phi^{n+1} = 0
		\ \mbox{ in } \Omega
\]
where $\Delta t>0$, $\phi^n\approx\phi(n\Delta  t)$
and $\phi^{(0)}=\phi_0$.

\cindex{problem!transport equation!unsteady}%
Let $t_n=n\Delta t$, $n\geq 0$.
The term $X^n(x)$ is the position at $t_n$ of the particle that
is in $x$ at $t_{n+1}$ and is transported by ${\bf u}^n$.
Thus, $X^n(x)=X(t_n,x)$ where
$X(t,x)$ is the solution of the differential equation
\[
    \left\{
	\begin{array}{rcl}
	    \Frac{{\rm d}X}{{\rm d}t} &=&
		{\bf u}(X(t,x),t)
		\ \ \ p.p. \ \ t\in\ ]t_n,t_{n+1}[,
		\\
	    X(t_{n+1},x) &=& x.
	\end{array}
   \right.
\]
Then $X^n(x)$ is approximated by the first-order Euler approximation
\[
	X^n(x) \approx x-\Delta t \, {\bf n}^n(x).
\]
This algorithm has been introduced by O.~Pironneau~(see e.g.~\citealp{Pir-1988}),
and is known as the method of characteristic in the finite difference context
and as the Lagrange-Galerkin in the finite element one.
%
\pindexopt{library}{CGAL, {\rm computational geometry}}%
The efficient evaluation of $\phi_h\circ X^n(x)$ in an unstructured 
mesh involves a hierarchical $d$-tree (quadtree, octree) data structure for the localization
of the element $K$ of the mesh that contains $x$.
When $d=3$ requires also sophisticated geometric predicates to test whether
$x\in K$ without rounding errors, and avoid to conclude that no elements
contains a point $x$ close to $\partial K$ up to rounding errors.
This problems is addressed in \Rheolef\  based on the \code{cgal} library.

\cindex{benchmark!rotating hill}%
The following code implements the classical rotating Gaussian hill
test case (see e.g.~\citealp{rui-tabata-2001}).

% --------------------
\myexamplelicense{convect.cc}
% --------------------

\subsubsection*{Comments}
% --------------------
\findex{compose}%
\findex{integrate}%
\clindex{integrate_option}%
\clindex{characteristic}%
The \code{characteristic} variable \code{X} implements the 
localizer $X^n(x)$:
\begin{lstlisting}[numbers=none,frame=none]
  characteristic X (-delta_t*uh);
\end{lstlisting}
Combined with the \code{compose} function,
it perform the composition $\phi_h \circ X^n$.
The right-hand side is then computed
by using the \code{integrate} function:
\begin{lstlisting}[numbers=none,frame=none]
  field lh = integrate (compose(phi_h, X)*psi, iopt);
\end{lstlisting}
\cindex{quadrature formulae!Gauss}%
\cindex{quadrature formulae!Gauss-Lobatto}%
Note the additional \code{iopt} argument to the \code{integrate} function.
By default, when this argument is omitted,
a Gauss quadrature formulae is assumed, and the number of
point is computed such that it integrate exactly $2k+1$ polynomials,
where $k$ is the degree of polynomials in $X_h$.
The Gauss-Lobatto quadrature formule is recommended for Lagrange-Galerkin methods.
Recall that this choice of quadrature formulae guaranties inconditionnal
stability at any polynomial order.
Here, we specifies a Gauss-Lobatto quadrature formulae that should be
exact for $k$ order polynomials.
The bilinear form is computed via the same quadrature formulae:
\begin{lstlisting}[numbers=none,frame=none]
  form  a  = integrate (c1*phi*psi + c2*dot(grad(phi),grad(psi)), iopt);
\end{lstlisting}
A test case is described by~\citet{PirTab-2010}:
we take $\Omega=]-2,2[^d$ and $T=2\pi$.
This problem provides an example for a convection-diffusion
equation and a known analytical solution:
\[
    \phi(t,x) = \exp\left( -\lambda t -r(t) |x-x_0(t)|^2 \right)
\]
where $\lambda=4\nu/t_0 \geq 0$ with $t_0>0$ and $\nu \geq 0$,
$x_0(t)$ is the moving center of the hill
and $r(t)=1/(t_0 + 4\nu t)$.
The source term is time-dependent: $\sigma(t)=\lambda - 2d\nu r(t)$
and has been adjusted such that the right-hand side is zero.
The moving center of the hill $x_0(t)$ is associated to the velocity field
${\bf u}(t,x)$ as:
\[
  \begin{array}{|c|l|l|}\hline
	d & {\bf u}(t,x) & x_0(t) \\ \hline \hline
	1 & 1/(2\pi)     & t/(2\pi) -1/2 \\ \hline
	2 & (y,-x)       & (-\cos(t)/2,\,\sin(t)/2) \\ \hline
	3 & (y,-x,0)     & (-\cos(t)/2,\,\sin(t)/2,0) \\ \hline
  \end{array}
\]
% --------------------------
\myexamplelicense{rotating-hill.h}
% --------------------------

\cindex{function!class-function object}%
\pindexopt{library}{STL, {\rm standard template library}}%
Note the use of a class-function {\tt phi} for
the implementation of $\phi(t)$ as a function of $x$.
Such programming style has been introduced in the
{\em standard template library}~\citep{musser-saini-1997}, which is a part
of the standard {\tt C++} library.
By this way, for a given $t$, $\phi(t)$ can be interpolated as an usual function
on a mesh.

\subsubsection*{How to run the program}
% ----------------------------------
\pindex{mkgeo_grid}%
\pindex{branch}%
\pindex{convect}%

 \begin{figure}[htb]
     \begin{center}
          \includegraphics[height=8cm]{convect-fig.png}
     \end{center}
     \caption{Animation of the solution of the rotating hill problem.}
     \label{fig-convect}
  \end{figure}

We assume that the previous code is contained in
the file \reffile{convect.cc}.
Then, compile the program as usual (see page~\pageref{makefile}):
\begin{verbatim}
  make convect
\end{verbatim}
and enter the commands:
\pindexopt{mkgeo_grid}{-a}%
\pindexopt{mkgeo_grid}{-b}%
\pindex{gnuplot}%
Running the one-dimensional test case:
\begin{verbatim}
  mkgeo_grid -e 500 -a -2 -b 2 > line2.geo
  ./convect line2.geo P1 > line2.branch
  branch line2.branch -gnuplot
\end{verbatim}
Note the hill that moves from $x=-1/2$ to $x=1/2$.
Since the exact solution is known, it is possible to analyze the error:

% ---------------------------
\myexamplelicense{convect_error.cc}
% ---------------------------
\findex{interpolate}%
The numerical error $\phi_h-\pi_h(\phi)$ is computed as:
\begin{lstlisting}[numbers=none,frame=none]
  field pi_h_phi = interpolate (Xh, phi(d,nu,t));
  field eh = phih - pi_h_phi;
\end{lstlisting}
and its $L^2$ norm is printed on the standard error.
Observe the use of the \code{branch} class
as both input and output field stream.
\begin{verbatim}
  make convect_error
  ./convect_error < line2.branch > line2-cmp.branch
  branch line2-cmp.branch -gnuplot
\end{verbatim}
\clindex{branch}%
\cindex{error analysis}%
The instantaneous $L^2(\Omega)$ norm is printed at each
time step and the total error in $L^2(]0,T[;L^2(\Omega))$
is finally printed at the end of the stream.
  \begin{figure}[htb]
     %\begin{center}
       \begin{tabular}{cc}
	  $\|\phi_h-\pi_h(\phi)\|_{L^2(L^2)}$
          &
	  $\|\phi_h-\pi_h(\phi)\|_{L^\infty(L^\infty)}$
	  \\
          \includegraphics{convect_cvge_Pk_err_l2_l2.pdf}
          &
  	  \hspace{-0.0cm}
          \includegraphics{convect_cvge_Pk_err_linf_linf.pdf}
       \end{tabular}
    %\end{center}
     \caption{Diffusion-convection when $d=1$ and $\nu=10^{-2}$:
        convergence versus $h$ and $\Delta t$ for $P_1$ and $P_2$ elements:
	(left)  in $L^2(L^2)$ norm;
	(right) in $L^\infty(L^\infty)$ norm.
     }
     \label{fig-convect-cvge-Pk}
  \end{figure}
A \code{P2} approximation can be used as well:
\begin{verbatim}
  ./convect line2.geo P2 > line2.branch
  branch line2.branch -gnuplot
  ./convect_error < line2.branch > line2-cmp.branch
\end{verbatim}
On Fig.~\ref{fig-convect-cvge-Pk}.left we
observe the $L^2(L^2)$ convergence versus $h$ for the $P_1$ and $P_2$
elements when $d=1$: the errors reaches a plateau that decreases
versus $\Delta t$.
On Fig.~\ref{fig-convect-cvge-Pk}.right
the $L^\infty(L^\infty)$ norm of the error presents
a similar behavior.
Since the plateau are equispaced,
the convergence versus $\Delta t$ is of first order.

These computation was performed for a convection-diffusion
problem with $\nu=10^{-2}$.
\cindex{problem!transport equation!unsteady}%
The pure transport problem ($\nu=0$, without diffusion) computation
is obtained by:
\begin{verbatim}
  ./convect line2.geo P1 0 > line2.branch
  branch line2.branch -gnuplot
\end{verbatim}

Let us turn to the two-dimensional test case:
\pindexopt{mkgeo_grid}{-c}%
\pindexopt{mkgeo_grid}{-d}%
\pindex{paraview}%
\begin{verbatim}
  mkgeo_grid -t 80 -a -2 -b 2 -c -2 -d 2 > square2.geo
  ./convect square2.geo P1 > square2.branch
  branch square2.branch 
\end{verbatim}
The visualization and animation are similar to 
those of the head problem previously presented
in paragraph~\ref{sec-heat}.
Go to the \code{WrapByScalar} entry in pipeline brower
and adjust eventually the \code{scale factor},
e.g. to~$3$.
Then, play the animation and observe the rotating hill.
The result is shown on Fig.~\ref{fig-convect}.
The error analysis writes:
\begin{verbatim}
  ./convect_error < square2.branch > square2-cmp.branch
  branch square2-cmp.brancha-nofill -bw
\end{verbatim}
From the \code{paraview} menu, you can visualize simultaneously both
the approximate solution and the Lagrange interpolate of the exact one.
For that purpose, go first to the \code{WrapByScalar} entry in pipeline brower
and adjust the \code{scale factor}, e.g. to~$3$.
Next, go to the \code{File->Open} menu and select
\code{square2-cmp-..vtk}.
In the \code{Filter->Recent} menu, select \code{WrapByScalar}.
In the \code{Properties} panel,
go to the \code{Scalars} entry and select \code{pi_h_phi}
and adjust the \code{scale factor} to~$3$.
Next, in the same panel, in the \code{Representation} entry,
choose \code{wireframe}.
Finally,
in the \code{Coloring} entry, choose \code{solid color},
and click on \code{Edit} for selecting e.g. the red color.
You are ready to click on the video \fbox{play} button, at the top of the window.
Observe the difference between the solution and its approximation.
For serious problem, the characteristic method has been superseded
by the discontinuous Galerkin one, that will be presented
in chapter~\ref{chap-dg}, page~\pageref{chap-dg}.
You are strongly encouraged to definitively
turn to discontinuous Galerkin method for convection dominant
and pure transport problems.

\pindexopt{mkgeo_grid}{-f}%
\pindexopt{mkgeo_grid}{-g}%
\pindexopt{mkgeo_grid}{-T}%
Finally, the three-dimensional case:
\begin{verbatim}
  mkgeo_grid -T 15 -a -2 -b 2 -c -2 -d 2 -f -2 -g 2 > cube2.geo
  ./convect cube2.geo P1 > cube2.branch
\end{verbatim}
\pindexopt{branch}{-volume}%
The visualization is similar to the two-dimensional case:
\begin{verbatim}
  branch cube2.branch 
  branch cube2.branch -volume
  \end{verbatim}