11_aderdg-user-fluxes.tex

\chapter{Generic ADER-DG and Finite Volume solvers}
\label{chapter:aderdg-user-defined-fluxes}


In our prevous chapter, the simulation run neither computes anything nor
does it yield any output.  
In this chapter, we introduce \exahype\ generic kernels: they allow the
user to specify flux and eigenvalue functions for the kernels,
but delegate the actual solve completely to \exahype. 
Furthermore, we quickly run through some visualisation and the handling of 
material parameters.
The resulting code can run in parallel and scale, but its single node
performance might remain poor, as we do not rely on \exahype's high performance
kernels. 
Therefore, the present coding strategy is particulary well-suited for rapid
prototyping of new solvers. 

\begin{design}
\exahype\ realises the Hollywood principle: ``Don't call us, we call you!''
The user does {\em not} write code that runs through some data structures, she
does {\em not} write code that determines how operations are invoked in parallel
and she does {\em not} write down an algorithmic outline in which order
operations are performed.

\exahype\ is the commander in chief: It organises all the data run throughs,
determines which core and node at which time invokes which operation and how the
operations are internally enumerated. Only for the application-specific
routines, it calls back user code. Application-specific means 
 \begin{itemize}
  \item how do we compute the flux function,
  \item how are the eigenvalues to be estimated,
  \item how do we convert the unknowns into records that can be plotted,
  \item how and where do we have to set non-homogeneous right-hand sides,
  \item \ldots
 \end{itemize}
\end{design}


\noindent
We have adopted this design pattern/paradigm from the Peano software that 
equipts \exahype\ with adaptive mesh refinement.
Our guideline illustrates all of these steps at hands of a solver for the Euler
equations.


\section{Establishing the data structures}\label{sec:euler-starting}


We follow the Euler equations in the form 
\begin{equation}\label{eq:euler}
\frac{\partial}{\partial t} \textbf{Q}
+
\boldsymbol{\nabla}\cdot\textbf{F}(\textbf{Q})
= 0
\quad \mbox{with} \quad
\textbf{Q} = \begin{pmatrix}
\rho\\\textbf{j}\\\ E
\end{pmatrix}
\quad \mbox{and} \quad
\textbf{F} = 
\begin{pmatrix}
\textbf{j}\\
\frac{1}{\rho}\textbf{j}\otimes\textbf{j} + p I \\
\frac{1}{\rho}\textbf{j}\,(E + p)
\end{pmatrix}
\end{equation}


\noindent
on a domain $\Omega $ supplemented by initial values
$\textbf{Q}(0)=\textbf{Q}_0$ and appropriate boundary conditions.
$\rho$ denotes the mass density,
$\textbf{j}\in\mathbb{R}^d$ denotes the momentum density, $E$ denotes the energy density, and $p$ denotes the fluid
pressure.
For our 2d setup, the velocity in z-direction $v_z$ is set to zero.
Introducing the adiabatic index $\gamma$, the fluid pressure is defined as
%
\begin{equation}
p = (\gamma-1)\,\left(E - \frac{1}{2}\textbf{j}^2/\rho\right).
\end{equation}


\noindent
Our example restricts to 
\[
 \Omega = (0,1)^d
\]
as complicated domains are subject of discussion in a separate Chapter \ref{chapter:complicated-domains}.
A corresponding specification file \texttt{euler-2d.exahype} for this
setup is
\begin{code}
exahype-project  Euler2d

  peano-kernel-path const    = ./Peano/
  exahype-path const         = ./ExaHyPE
  output-directory const     = ./ApplicationExamples/EulerFlow

  computational-domain
    dimension const          = 2
    width                    = 1.0, 1.0
    offset                   = 0.0, 0.0
    end-time                 = 0.4
  end computational-domain
  
  solver ADER-DG MyEulerSolver
    variables const          = 5
    order const              = 3
    maximum-mesh-size        = 0.1
    time-stepping            = global
    type const               = nonlinear
    terms const              = flux
    optimisation const       = generic
    language const           = C

    plot vtu::Cartesian::cells::ascii EulerWriter
      variables const = 5
      time            = 0.0
      repeat          = 0.5E-1
      output          = ./variables
    end plot
  end solver
end exahype-project  
\end{code}

\noindent
The spec. file sets some paths in the preamble before it specifies the computational
domain and the simulation time frame in the environment
\texttt{computational-domain}.

In the next lines, a solver of type \texttt{ADER-DG} is specified and assigned the name
\texttt{MyEulerSolver}.  The kernel type of the solver is set to
\texttt{nonlinear}. 
We will not employ any \texttt{optimisation}. Instead we use
some generic ADER-DG kernels that can be applied to virtually
any hyperbolic problem. In this case, 
the user is only required to provide the \exahype\ engine with
problem specific flux (and eigenvalues) definitions.
Here, we tell the toolkit that we want to specify a
conservative flux via parameter \texttt{terms}.
Other options are possible (Sec. \ref{sec:solver-configuration}).

Within the solver environment, there is also a plotter specified and configured.
This plotter is further configured to write out a snapshot of the solution
of the associated solver every $0.05$ time intervals.  
The first snapshot is set to be written at time $t=0$.
The above plotter statement creates a plotter file \texttt{MyPlotter} that
allows you to alter the plotter's behaviour.

Once we are satisfied with the parameters in our \exahype\ specification file, we hand it over to the \exahype\ toolkit:
\begin{code}
  ./Toolkit/toolkit.sh Euler2d.exahype
\end{code}


\begin{design}
The formulation (\ref{eq:euler}) lacks a dependency on the spatial position $x$.
As such, it seems that \exahype\ does not support spatially varying
fluxes/equations.
This impression is wrong.  
Our design philosophy is that spatial variations actually are material
parameters, i.e.,~we recommend that you introduce material parameters for your
spatial positions $x$ or quantities derived from there.
The Section \ref{section:advanced-features:material-parameters} details the
usage of material parameters.
\end{design}


\section{Study the generated code}

We obtain a new directory equalling \texttt{output-directory} from the
specification file if such a directory hasn't existed before.
Lets change to this path. 
The directory contains a makefile, and it contains a bunch of generated files:

\begin{figure}[h]
\begin{center}
  \includegraphics[width=0.5\textwidth]{sketches/simplest-class-architecture.pdf}
\end{center}
\caption[Sketch (Inkscape): Class architecture]{Simplest class architecture for ExaHyPE solvers.}
\end{figure}

\begin{itemize}
  \item \texttt{MyAbstractEulerSolver} is a class holding solely glue code,
  i.e.~this class is not to be altered. It invokes for example the generic
  non-linear kernels. Once you rerun \exahype's toolkit, it will overwrite this
  class. The class implements some interfaces/abstract classes from \exahype's
  kernel. There is usually no need to study those at all.
  \item \texttt{MyEulerSolver} is the actual solver. This is where users
  implement the solver's behaviour, i.e.~this class is to be befilled with
  actual code. Some methods are pregenerated (empty) to make your life easier.
  Have a close look into the header file---by default, the toolkit places all
  documentation in the headers (cmp.~Section \ref{chapter:code-docu})---which
  contains hint what behaviour can be added through a particular method.
  \item \texttt{MyPlotter} is a class that allows us to tailor \exahype's
  output. It implements an interface and can quite freely be adopted. For the
  time being, the default variant dropped by the toolkit does the job. This
  class implements an interface from the \exahype\ kernel. There's usually no
  need to study this interface---again, the toolkit generates quite some
  comments to the generated headers that you may redefine and implement.
\end{itemize}

\noindent
Before we continue our tour de \exahype, it might be reasonable to compile the
overall code once and to doublecheck whether we can visualise the resulting outputs. 
This should all be straightforward. 
Please verify that any output you visualise holds only garbage
since no initial values are set so far.


\section{Working with the arrays}


\begin{design}
 \exahype\ relies on plain arrays to encode all unknowns and fluxes. We
 however leave it to the user to decide whether she prefers to work with these
 plain double arrays (similar to Fortran) or with some symbolic identifiers.
 \exahype\ also provides support for the latter. Note that symbolic identifiers 
 may degrade performance.
\end{design}

\subsection*{Variant A: Sticking to the arrays}

If you prefer to work with plain double arrays always, it is sufficient to tell
the solver how many unknowns your solution value contains:

\begin{code}
   solver [...]
     variables const = 5
     [...]
   end solver
\end{code}

\noindent
It is up to you to keep track which entry in any double array has which
semantics. 
In the Euler equations, you have for example to keep in mind that the fifth
entry always is the energy $E$.
All routines of relevance are plain double pointers to arrays.
So an operation altering the energy typically reads \texttt{Q[4]} (we work in the C/C++ language and thus start to count with 0).
All fluxes are two-dimensional double arrays, i.e.~are represented as matrices (that is, they are not guaranteed to be continous storage).



\subsection*{Variant B: Working with symbolic identifiers}
\begin{warning}
Working with symbolic identifiers may degrade performance!
\end{warning}

As alternative to plain arrays, you may instruct \exahype\ about the
unknowns you work with:

\begin{code}
  solver [...]
    variables const = rho:1,j:3,E:1
    [...]
  end solver
\end{code}

\noindent
This specification is exactly equivalent to the previous declaration. 
It tells \exahype\ that there are five unkonwns held in total.
Two of them are plain scalars, the middle one is a vector with three entries. 
Also, all operation signatures remain exactly the same. 
Yet, the toolkit now creates a couple of additional classes that can be
wrapped around the array.
These classes do not copy any stuff around or themselves actually alter the
array.
They provide however routines to access the array entries through their unknown
names instead of plain array indices:

\begin{code}
void Euler2d::MyEulerSolver::anyRoutine(..., double *Q) {
  Variables vars(Q); // we wrap the array
   
  vars.rho() = 1.0;  // accesses Q[0]
  vars.j(2)  = 1.0;  // accesses Q[3]
  vars.E()   = 1.0;  // accesses Q[4]
}
\end{code}

\noindent
The wrapper allows us to ``forget'' about the mapping of the unknown onto array
indices.
We may use variable names instead.
Besides the \texttt{Variables} wrapper, the toolkit also generates wrappers
around the matrices as well as read only variants of the wrappers that are to be
used if you obtain a \texttt{const double const*} argument.

We note that there are subtle syntactic differences between the plain array
access style and the wrappers.
The latter typically rely on \texttt{()} brackets similar to Matlab.
Without going into details, we want to link to two aspects w.r.t.~symbolic
accessors:
\begin{enumerate}
  \item The wrappers are plain wrappers around arrays. You may intermix both
  access variants---though you have to decide whether you consider this to be
  good style. Furthermore, all vectors offered through the wrapper provide a
  method \texttt{data()} that again returns a pointer to the respective
  subarray. The latter is useful if you prefer to work with symbolic identifiers
  but still have to connect to array-based subroutines.
  \item All wrapper vector and matrix classes stem from the linear algebra
  subpackage of Peano's technical architecture (\texttt{tarch::la}). We refer to
  Peano's documentation for details, but there are plenty of frequently used
  operations (norms, scalar products, dense matrix inversion, \ldots) shipped
  along with Peano that are all available to \exahype.
\end{enumerate}


\noindent
From hereon, most of our examples rely on Variant B, i.e.~on the symbolic
names.
This makes a comparison to text books easier, but be aware that it be less 
efficient than a direct implementation with arrays.
The cooking down to a plain array-based implementation is 
straightforward.

\begin{design}
All of our signatures rely on plain arrays. Symbolic access is offered via a
wrapper and always is optional. Please note that you might also prefer to wrap
coordinates, e.g., with \texttt{tarch::la::Vector} classes to have Matlab-style
syntax.  
\end{design}

%\noindent
%We do not apply such a sophisticated syntactic connvention here to keep the
%guidebook as simple to understand as possible.
% ^-- I don't understand this sentence


\section{Setting up the experiment}\label{section:setup}

We return to the Euler flow solver:
To set up the experiment, we open the implementation file 
\texttt{MyEulerSolver.cpp} of the class \texttt{MyEulerSolver}. 
There is a routine \texttt{adjustPointSolution}
that allow us to setup the initial grid.
Alternatively, we can also use \texttt{adjustSolution}.
One routine works point-wisely, the other one
hands over a whole patch.
The implementation of the initial values might look as follows\footnote{The
exact names of the parameters might change with newer \exahype\ versions and
additional parameters might drop in, but the principle always remains the
same.}:
\begin{code}
void Euler2d::MyEulerSolver::adjustPointSolution(
  const double *const x,
  const double w,
  const double t,
  const double dt,
  double *Q
) {
 if (tarch::la::equals( t,0.0,1e-15 )) {
   Variables vars(Q);
   const double GAMMA = 1.4;
   
   vars.j( 0.0, 0.0, 0.0 );
   vars.rho() = 1.0;    
   vars.E() =
        1. / (GAMMA - 1) +
        std::exp(-((x[0] - 0.5) * (x[0] - 0.5) + (x[1] - 0.5) * (x[1] - 0.5)) /
                 (0.05 * 0.05)) *
            1.0e-3;
   }
}
\end{code}



\noindent
The above routine enables us to specify time dependent solution values.
It further enables us to add source term contributions to the solution values.

As we only want to impose initial conditions so we check if the time
\texttt{t} is zero. As these values are floating point values, standard bitwise C comparison
would be inappropriate. We rely here on a routine coming along with the linear
algebra subroutines of Peano to check for ``almost equal''  up to machine precision.


We conclude our experimental setup with a compile run and a first test run.
This time we get a meaningful image and definitely not a program stop for an
assertion, as we have set all initial values properly.
However, nothing happens so far; which is not surprising given that we haven't
specified any fluxes yet (a plot should be similar to fig \ref{fig:10_ic}).

\begin{figure}[h]
\begin{center}
  \includegraphics[width=0.75\textwidth]{screenshots/10_initial-condition.png}  
\end{center}
\caption[Paraview initial conditions of an EulerFlow 2D experiment]{The rest mass density $Q_0$ in the EulerFlow 2D setup. The
scalar value is encoded both in height and color. So far, the PDE
has not been specified so it stays constant during the time evolution.}
\label{fig:10_ic}
\end{figure}


\begin{design}
The adoption of the solution (to initial values) as well as source terms are
protected by an additional query \texttt{hasToAdjustSolution}. This allows
\exahype\ to optimise the code aggressively: The user's routines are invoked
only for regions where \texttt{hasToAdjustSolution}'s result holds. For all the
remaining computational domain, \exahype\ uses numerical subroutines that are
optimised to work without solution modifications or source terms.
\end{design}



\section{Realising the fluxes}

To actually implement the Euler equations, we have to realise the fluxes
into the code. We do so by filling the functions \texttt{flux} and
\texttt{eigenvalues} in file \texttt{MyEulerSolver.cpp} with code. 


\subsection*{Variant A}

In this first part, we present the realisation based upon plain arrays:
\begin{code}
void Euler::FirstEulerSolver::flux(const double* const Q, double** F) {
  // Dimensions             = 2
  // Number of variables    = 5 (#unknowns + #parameters)
  const double GAMMA = 1.4;

  const double irho = 1.0 / Q[0];
  const double p =
      (GAMMA - 1) * (Q[4] - 0.5 * (Q[1] * Q[1] + Q[2] * Q[2]) * irho);

  double* f = F[0];
  double* g = F[1];

  // f is the flux on faces with a normal along x direction
  f[0] = Q[1];
  f[1] = irho * Q[1] * Q[1] + p;
  f[2] = irho * Q[1] * Q[2];
  f[3] = irho * Q[1] * Q[3];
  f[4] = irho * Q[1] * (Q[4] + p);
  // g is the flux on faces with a normal along y direction
  g[0] = Q[2];
  g[1] = irho * Q[2] * Q[1];
  g[2] = irho * Q[2] * Q[2] + p;
  g[3] = irho * Q[2] * Q[3];
  g[4] = irho * Q[2] * (Q[4] + p);
}


void Euler::FirstEulerSolver::eigenvalues(const double* const Q, 
  const int normalNonZeroIndex, double* lambda) {
  // Dimensions             = 2
  // Number of variables    = 5 (#unknowns + #parameters)
  const double GAMMA = 1.4;

  double irho = 1.0 / Q[0];
  double p = (GAMMA - 1) * (Q[4] - 0.5 * (Q[1] * Q[1] + Q[2] * Q[2]) * irho);

  double u_n = Q[normalNonZeroIndex + 1] * irho;
  double c = std::sqrt(GAMMA * p * irho);

  lambda[0] = u_n - c;
  lambda[1] = u_n;
  lambda[2] = u_n;
  lambda[3] = u_n;
  lambda[4] = u_n + c;
}

\end{code}

\subsection*{Variant B}
\begin{warning}
Working with symbolic identifiers may degrade performance!
\end{warning}

Alternatively, we can work with symbolic identifiers if we have specified the
unknowns via \texttt{rho:1,j:3,E:1}:
\begin{code}
void Euler2d::MyEulerSolver::flux(
  const double* const Q,
  double** F
) {
  ReadOnlyVariables vars(Q);
  Fluxes f(F);

  tarch::la::Matrix<3,3,double> I;
  I = 1, 0, 0,
      0, 1, 0,
      0, 0, 1;

  const double GAMMA = 1.4;
  const double irho = 1./vars.rho();
  const double p = (GAMMA-1) * (vars.E() - 0.5 * irho * vars.j()*vars.j() );

  f.rho ( vars.j()                                 );
  f.j   ( irho * outerDot(vars.j(),vars.j()) + p*I );
  f.E   ( irho * (vars.E() + p) * vars.j()         );
}
\end{code}


\begin{code}
void Euler2d::MyEulerSolver::eigenvalues(
  const double* const Q,
  const int normalNonZeroIndex,
  double* lambda
) {
  ReadOnlyVariables vars(Q);
  Variables eigs(lambda);

  const double GAMMA = 1.4;
  const double irho = 1./vars.rho();
  const double p = (GAMMA-1) * (vars.E() - 0.5 * irho * vars.j()*vars.j() );

  double u_n = vars.j(normalNonZeroIndex) * irho;
  double c   = std::sqrt(GAMMA * p * irho);

  eigs.rho()=u_n - c;
  eigs.E()  =u_n + c;
  eigs.j(u_n,u_n,u_n);
}
\end{code}


The  implementation of function \texttt{flux} is very straightforward. 
Again, the details are only subtle:
We wrap up the arrays \texttt{Q} and \texttt{F} in
wrappers of type \linebreak \texttt{ReadOnlyVariables} and \texttt{Fluxes}.
Similar to \texttt{Variables}, the definitions of \linebreak
\texttt{ReadOnlyVariables} and
\texttt{Fluxes} were generated according to the variable list we have specified
in the specification file.
While \texttt{Fluxes} indeed is a 1:1 equivalent to \texttt{Variables}, we have
to use a read-only variant of \texttt{Variables} here as the input array
\texttt{Q} is protected. 
The read-only symbolic wrapper equals exactly its standard counterpart but lacks
all setters.


The pointer \texttt{lambda} appearing in the signature and body of function 
\texttt{eigenvalues} has the size as the vector of conserved variables.
It thus makes sense to wrap it in an object of type \texttt{Variables}, too.

The normal vectors that are involved in \exahype 's ADER-DG kernels always coincide with the (positively signed) 
Cartesian unit vectors. Thus, the \texttt{eigenvalues} function is only supplied with the index of the single component 
that is non-zero within these normal vectors.
In function \texttt{eigenvalues}, the aforementioned index corresponds to the parameter
\texttt{normalNonZeroIndex}. 


\section{Supplementing boundary conditions}

\exahype\ offers a generic API to implement any kind of \emph{local}
boundary conditions.
Local here refers to element-wise boundary conditions where the solver in one
cell does not have to communicate with other cells.
That is, you can straightforwardly\footnote{In this guidebook, this implies
that we stick to those guys only.} implement

\begin{itemize}
\item outflow boundary conditions (also sometimes refered to as "null boundary
conditions" or "none boundary conditions"),
\item exact boundary conditions or 
\item reflection symmetries.
\end{itemize}


\noindent
The signature to implement your boundary conditions reads
\begin{code}
void Euler::MyEulerSolver::boundaryValues(
    const double* const x, const double t,  const double dt,
    const int faceIndex, const int normalNonZero,
    const double* const fluxIn,const double* const stateIn,
    double* fluxOut, double* stateOut)
    {
      ...
    }
\end{code}

For now we simply set stateOut = stateIn and fluxOut = fluxIn.
A rebuild and rerun should yield results similar to
figure \ref{fig:10_tevolution}.

\begin{figure}[h]
\begin{center}
  \includegraphics[width=0.45\textwidth]{screenshots/10_early-stage.png}  
  \includegraphics[width=0.45\textwidth]{screenshots/10_later-stage.png}  
\end{center}
\caption[Paraview screenshots of the EulerFlow time evolution.]{Time evolution of $Q_0$, now with the Hydrodynamics PDEs implemented.
The left figure shows an early time while the right figure shows a later time.}
\label{fig:10_tevolution}
\end{figure}



\begin{figure}[htb]
\begin{center}
  \includegraphics[width=0.5\textwidth]{sketches/cube-names.png}
\end{center}
\caption[Rendering (Blender): Peano face indices]{The face indexes as named in the ExaHyPE code}
\label{fig:cube-names}
\end{figure}


\noindent
The input arguments (marked with the C \texttt{const} modifier) offer
the position of a cell's boundary and time as well as the local timestep of the
cell, and the cell's state and flux variables.
Cell hereby always is inside the computational domain, i.e.~\exahype\ queries
the boundary values from the inner cell's point of view.


The two variables \texttt{faceIndex} and \texttt{normalNonzero} to answer
the questions at which boundary side we currently are. The face index
decodes as
\begin{code}
0-left, 1-right, 2-front, 3-back, 4-bottom, 5-top
\end{code}
or in other terms
\begin{code}
0 x=xmin   1 x=xmax, 2 y=ymin    3  y=ymax  4 z=zmin     5 z=zmax
\end{code}
This is also encoded in Figure \ref{fig:cube-names}.


There is a big difference between Finite Volume boundary conditions as
typically found in literature and our ADER-DG boundary conditions.
For a DG method in general it is not sufficient
to just prescribe the \texttt{stateOut}, i.e.~the values, along the boundary
that are just outside of the domain.
We further need to prescribe the normal fluxes (\texttt{fluxOut}; Figure
\ref{fig:boundary-conditions}); unless if use lowest order DG, i.e.~Finite Volumes.
For Finite Volumes, the fluxes are not required as they directly result from
the values ``outside'' of the domain.



\begin{figure}[htb]
\begin{center}
  \includegraphics[width=0.5\textwidth]{sketches/boundary-conditions.pdf}
\end{center}
  \caption[Sketch (Inkscape): 1d boundary conditions]{
    1-dimensional sketch of boundary conditions in \exahype's ADER-DG.
    The state inside the domain and the fluxes from this side are passed to the
    boundary treatment by the solver. They are therefore const. A boundary
    condition is imposed through prescribing the state outside of the domain
    (thus describing the solution's jump) and the associated flux.
  }
\label{fig:boundary-conditions}
\end{figure}


Furthermore, ADER-DG methods require us to prescribe the state just outside of
the domain (\texttt{stateOut}) and the corresponding fluxes (\texttt{fluxOut})
at all the temporal interpolation points.
Within each cell, ADER-DG relies on a local space-time formulation: it derives
all operators from integrals $\int_{t}^{t+\Delta t} \ldots dt$.
If you use a type of global time stepping where all cells march with exactly the
same time step size, then the integral formulation translates into an arithmetic
exercise.
At the end of the day, the integrals over all state and flux variables along the
face over the time span $(t,t+\Delta t)$ are required to set proper boundary
conditions.


The story is slightly different if you use some local time stepping. In this
case, \exahype\ has to know the boundary conditions over time. 
It has to know the polynomial description of the state and flux solution over
the whole time span.

\subsection{Integration of Exact Boundary Conditions}\label{sub:exact-BC}

In classical Runge-Kutta DG schemes, only a weak form in space of the PDE is obtained, and the time is kept continuous, thus reducing the problem to a nonlinear system of ODEs, which is subsequently integrated in time using classical ODE solvers. In the ADER-DG framework, a completely different paradigm is used.  Integration in time is solved by a special weak form of the PDE, which is built on a set of test functions defined on prisms in both space and time. The ADER-DG scheme can be regarded as a predictor-corrector method. First a local time evolution of each cell is calculated in a predictor step and then the influence of the neighbours is added in the single corrector step.
Let us consider the following general PDE of the form:
\begin{equation}\label{eq:general-equation}
  \frac{\partial}{\partial t} \textbf{Q}
  +
  \boldsymbol{\nabla}\cdot
  \underbrace{\textbf{F}(\textbf{Q})}
  _{
    \texttt{flux}}
  +
  \underbrace{ \textbf{B}(\textbf{Q}) \cdot \boldsymbol{\nabla}\textbf{Q} }
  _{\texttt{ncp}}
  = 
  \underbrace{
   \textbf{S}(\textbf{Q})
  }
  _{
    \texttt{source}
  }
\end{equation}
Since we adopt discontinuous Galerkin finite element method, the numerical solution $\textbf{u}_h$ restricted to a control volume $\Omega_i$ is written at time $\textit{t}^n$ in terms of some nodal basis functions $\Phi_l(\textbf{x})$ and some unknown degrees of freedom $\hat{\textbf{u}}_{i,l}^n$:
\begin{equation}\label{eq:numerical-sol}
\textbf{u}_h(\textbf{x},\textit{t}^n)|_{\Omega_i} = \sum_{i} \textbf{u}_{i,l}\Phi_l(\textbf{x}) := \hat{\textbf{u}}_{i,l}^n\Phi_l(\textbf{x})
\end{equation}
where $l = (l_1,l_2,l_3)$ is a multi-index and the spatial basis functions $\Phi_l(\textbf{x}) = \varphi_{l_1}(\xi) \varphi_{l_2}(\eta)\varphi_{l_3}(\zeta)$ are generated via tensor products of one-dimensional basis functions on the reference element $\mathbf{\xi}=(\xi,\eta,\zeta) \in [0,1]^d$.
In order to derive the ADER-DG method, we first multiply the governing PDE system with a test function $\mathbf{\Phi}_k$ and integrate over space-time control volume $\Omega_i \times [t^n;t^{n+1}]$ to get:
\begin{equation}\label{eq:integral-eq}
 \int_{t_n}^{t^{n+1}} \int_{\Omega_i} \mathbf{\Phi}_k \frac{\partial}{\partial t} \textbf{Q} \textrm{ } d\textbf{x} \textrm{ }dt
  +
  \int_{t_n}^{t^{n+1}} \int_{\Omega_i} \mathbf{\Phi}_k \left(\boldsymbol{\nabla}\cdot\textbf{F}(\textbf{Q}) +
 \textbf{B}(\textbf{Q}) \cdot \boldsymbol{\nabla}\textbf{Q}\right) \textrm{ } d\textbf{x} \textrm{ }dt
  = 
     \int_{t_n}^{t^{n+1}} \int_{\Omega_i} \mathbf{\Phi}_k \textbf{S}(\textbf{Q}) \textrm{ } d\textbf{x} \textrm{ }dt
\end{equation}
with $d\textit{x}=dx \textrm{ } dy \textrm{ } dz$. Integrating the first term by parts in time and using Eq. \ref{eq:numerical-sol} as well as integrating the flux divergence term by parts in space and introducing the element-local space-time predictor $\textbf{q}_h(\textbf{x},t)$ instead of $\textbf{Q}$, the weak formulation in Eq. \ref{eq:integral-eq} can be written as:
\begin{equation}\label{eq:weak-form}
\begin{split}
\left( \int_{\Omega_i} \mathbf{\Phi}_k \mathbf{\Phi}_l \textrm{ } d\textbf{x} \right) \left( \hat{\textbf{u}}_{i,l}^{n+1} - \hat{\textbf{u}}_{i,l}^n \right)
 \textrm{ } +\textrm{ }
  \int_{t_n}^{t^{n+1}} \int_{^\partial\Omega_i} \mathbf{\Phi}_k \mathcal{G}(\textbf{q}_h^{-},\textbf{q}_h^{+}) \cdot \textbf{n} \textrm{ }dS \textrm{ }dt 
\textrm{ }-\textrm{ }
   \int_{t_n}^{t^{n+1}} \int_{\Omega_i^{o}} \nabla\mathbf{\Phi}_k \cdot \textbf{F}(\textbf{q}_h)\textrm{ }  d\textbf{x}\textrm{ }dt \\
\textrm{ }+\textrm{ }
   \int_{t_n}^{t^{n+1}} \int_{\Omega_i^{o}} \mathbf{\Phi}_k \left(\textbf{B}( \textbf{q}_h) \cdot \nabla \textbf{q}_h \right) \textrm{ } d\textbf{x}\textrm{ }dt
 \textrm{ } = \textrm{ }
     \int_{t_n}^{t^{n+1}} \int_{\Omega_i} \mathbf{\Phi}_k \textbf{S}(\textbf{q}_h) \textrm{ } d\textbf{x} \textrm{ }dt
     \end{split}
\end{equation}
where  $\textbf{q}_h^{-}$ and $\textbf{q}_h^{+}$ are boundary extrapolated values of space-time predictor and the boundary integral contains the approximate Riemann solver $\mathcal{G}$ which accounts for the jumps across element interfaces, also in the presence of non-conservative products.\footnote{Developments presented in subsection \ref{sub:exact-BC} up to this point come from the following source: Dumbser, M.; Fambri, F.; Tavelli, M.; Bader, M.; Weinzierl, T. Efficient Implementation of ADER Discontinuous Galerkin Schemes for a Scalable Hyperbolic PDE Engine. Axioms 2018, 7, 63.}

In the ExaHyPE implementation, Riemann solver accepts input (stateOut and fluxOut) based on time-averaged space-time predictor $ \bar{\textbf{q}}_h(\textbf{x})$. Since we operate locally and in scaled reference element of size $(0,1)^d \times [0,1]$, averaging of the predictor over time, can be achieved via the following quadrature formula:
\begin{equation}\label{eq:avg}
\begin{split}
\int_{t^{n}}^{t^{n+1}} \textbf{q}_h(\textbf{x},t) dt = \Delta t \int_0^1 \hat{\textbf{q}}_h(\textbf{x},\hat{t}) d\hat{t}
\quad \quad \Rightarrow \bar{\textbf{q}}_h(\textbf{x}) := \frac{\int_{t^{n}}^{t^{n+1}} \textbf{q}_h(\textbf{x},t) dt}{\Delta t}
\quad \quad  \Rightarrow \bar{\textbf{q}}_h(\textbf{x}) = \sum_{i=0}^{N} w_i \hat{\textbf{q}}_h(\textbf{x},\hat{t}_i)
\end{split}
\end{equation}
This is exactly the form of an integral that is being calculated in spaceTimePredictor() function of ExaHyPE as the last step.

To be compliant with this logic, an input based on time-averaged solution needs to be passed to the Riemann solver also in case when the domain boundary is encountered. This can be done using exact boundary conditions. To this end the value of a sought function at a given space locations, in particular at the boundary, needs to be known at all times of interest. Actual simulation time for the reference cell is recovered from the given Gauss-Legendre quadrature points by an appropriate mapping, i.e. ti = t + xi*dt. User needs to provide a function setYourExactData(), which calculates solution values at the boundary point x at the given simulation time point ti. Subsequently, a flux() function needs to be called on the boundary. As explained both stateOut and fluxOut need to be averaged in time, which is done in the last loop of boundaryValues() function presented below, and it is done similarily to what has been presented for the averaging of the predictor in Eq. \ref{eq:avg}.

We supply here a small example how to correctly integrate
imposed Boundary Conditions in the ADER-DG scheme. 

\begin{code}
#include "kernels/GaussLegendreQuadrature.h"

void DemoADERDG::DemoSolver::boundaryValues(
   const double* const x,const double t, const double dt,
   const int faceIndex, const int normalNonZero,
   const double * const fluxIn, const double* const stateIn,
   double *fluxOut, double* stateOut) {

  // Defining your constants
  const int nVar = DemoADERDG::AbstractSolver::NumberOfVariables;
  const int order = DemoADERDG::AbstractSolver::Order;
  const int basisSize = order + 1;
  const int nDim = DIMENSIONS;

  double Qgp[nVar];
  std::memset(stateOut, 0, nVar * sizeof(double));
  std::memset(fluxOut, 0, nVar * sizeof(double));

  double F[nDim][nVar];

  // Integrate solution in gauss points (Qgp) in time
  for(int i=0; i < basisSize; i++)  { // i == time
     const double weight = kernels::gaussLegendreWeights[order][i];
     const double xi = kernels::gaussLegendreNodes[order][i];
     double ti = t + xi * dt;

     setYourExactData(x, Qgp, &ti);
     flux(Qgp, F); // calling your Solver's flux function
     for(int m=0; m < nVar; m++) {
        stateOut[m] += weight * Qgp[m];
        fluxOut[m] += weight * F[normalNonZero][m];
     }
  }
}
\end{code}

\subsection{Outflow Boundary Conditions}
As another example, outflow boundary conditions are archieved
by just copying the fluxes and states.

\begin{code}
void DemoADERDG::DemoSolver::boundaryValues(...) {
  for(int i=0; i < DemoADERDG::AbstractSolver::NumberOfVariables; i++) {
      fluxOut[i] = fluxIn[i];
      stateOut[i] = stateIn[i];
  }
}
\end{code}


\begin{design}
  We stick to $\Omega = (0,1)^d$ without MPI here. Please read carefully through
  Chapter \ref{chapter:complicated-domains} if you need more sophisticated
  boundary conditions or if precise boundary conditions form an essential part
  of your simulation.
\end{design}


\section{Finite Volumes}\label{sec:FV-intro}

If you prefer your code to run with Finite Volumes,  \exahype\ sticks to all
paradigms introduced so far. 
The user has to provide basically fluxes, eigenvalues and boundary/initial
conditions.
All such information is already available here. 
Consequently, switching from ADER-DG to Finite Volumes is a slight modification
of the configuration file and a rerun of the toolkit:

\begin{code}
  solver Finite-Volumes MyFVEulerSolver
    variables         = rho:1,j:3,E:1
    patch-size        = 3
    maximum-mesh-size = 0.1
    time-stepping     = global
    kernel            = generic::godunov
    language          = C
  end solver
\end{code}

\todo[inline]{This needs an update.}

\begin{design}
 With ADER-DG, \exahype\ embeds a higher order polynomial into each grid cell.
 With Finite Volumes, \exahype\ embeds a small patch (a regular Cartesian grid)
 into each grid cell.
\end{design}

\noindent
We switch the solver type to \texttt{Finite-Volumes} and fix a patch
resolution, before we recompile the \exahype\ application and run a Finite
Volume solver instead of the ADER-DG variant.
There are only subtle differences to take into account:

\begin{itemize}
  \item If you (as discussed later in this document) fuse the ADER-DG solver
    with a Finite Volumes solver,
    \texttt{patch-size} typically should be chosen as $2 \cdot$ \texttt{order}
    $+1$. This ensures that the admissible time step sizes of the DG scheme
    matches the time step sizes of the Finite Volumes formulation.
  \item Boundary treatment of Finite Volume solvers simpler than for Finite
    Volumes. It only requires us to prescribe the state variables (unknowns)
    outside of the domain. All fluxes then are determined from hereon. As a
    consequence, the \texttt{boundaryValues} routine is a briefer signature.
    While you can share flux/eigenvalue computations, e.g., between an ADER-DG
    and FV solver, the boundary treatment has to be realised independently.
\end{itemize}