11_aderdg-user-fluxes.tex 34.3 KB
Newer Older
1
\chapter{Generic ADER-DG and Finite Volume solvers}
2
\label{chapter:aderdg-user-defined-fluxes}
3

4

tobias's avatar
tobias committed
5
6
7
8
9
10
11
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.
12
13
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
tobias's avatar
tobias committed
14
15
kernels. 
Therefore, the present coding strategy is particulary well-suited for rapid
16
prototyping of new solvers. 
17

tobias's avatar
tobias committed
18
\begin{design}
tobias's avatar
tobias committed
19
\exahype\ realises the Hollywood principle: ``Don't call us, we call you!''
tobias's avatar
tobias committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
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.


46
\section{Establishing the data structures}\label{sec:euler-starting}
tobias's avatar
tobias committed
47
48
49


We follow the Euler equations in the form 
50
\begin{equation}\label{eq:euler}
51
52
53
54
\frac{\partial}{\partial t} \textbf{Q}
+
\boldsymbol{\nabla}\cdot\textbf{F}(\textbf{Q})
= 0
55
\quad \mbox{with} \quad
56
\textbf{Q} = \begin{pmatrix}
57
\rho\\\textbf{j}\\\ E
58
\end{pmatrix}
59
\quad \mbox{and} \quad
60
61
\textbf{F} = 
\begin{pmatrix}
62
63
64
\textbf{j}\\
\frac{1}{\rho}\textbf{j}\otimes\textbf{j} + p I \\
\frac{1}{\rho}\textbf{j}\,(E + p)
65
\end{pmatrix}
66
\end{equation}
tobias's avatar
tobias committed
67
68
69
70
71


\noindent
on a domain $\Omega $ supplemented by initial values
$\textbf{Q}(0)=\textbf{Q}_0$ and appropriate boundary conditions.
tobias's avatar
tobias committed
72
$\rho$ denotes the mass density,
73
74
75
$\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.
tobias's avatar
tobias committed
76
Introducing the adiabatic index $\gamma$, the fluid pressure is defined as
77
78
%
\begin{equation}
79
p = (\gamma-1)\,\left(E - \frac{1}{2}\textbf{j}^2/\rho\right).
80
\end{equation}
tobias's avatar
tobias committed
81
82
83
84
85
86
87
88


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

Anne Reinarz's avatar
Anne Reinarz committed
94
  peano-kernel-path const    = ./Peano/
95
  exahype-path const         = ./ExaHyPE
96
  output-directory const     = ./ApplicationExamples/EulerFlow
tobias's avatar
tobias committed
97

98
  computational-domain
99
    dimension const          = 2
tobias's avatar
tobias committed
100
    width                    = 1.0, 1.0
101
    offset                   = 0.0, 0.0
di25cox's avatar
di25cox committed
102
    end-time                 = 0.4
103
  end computational-domain
Tobias Weinzierl's avatar
Tobias Weinzierl committed
104
  
105
  solver ADER-DG MyEulerSolver
106
107
    variables const          = 5
    order const              = 3
ga96nuv's avatar
ga96nuv committed
108
109
    maximum-mesh-size        = 0.1
    time-stepping            = global
110
111
112
    type const               = nonlinear
    terms const              = flux
    optimisation const       = generic
113
    language const           = C
tobias's avatar
tobias committed
114

Anne Reinarz's avatar
Anne Reinarz committed
115
116
117
118
119
    plot vtu::Cartesian::cells::ascii EulerWriter
      variables const = 5
      time            = 0.0
      repeat          = 0.5E-1
      output          = ./variables
120
121
    end plot
  end solver
Tobias Weinzierl's avatar
Tobias Weinzierl committed
122
end exahype-project  
123
124
125
\end{code}

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

In the next lines, a solver of type \texttt{ADER-DG} is specified and assigned the name
Tobias Weinzierl's avatar
Tobias Weinzierl committed
131
\texttt{MyEulerSolver}.  The kernel type of the solver is set to
132
133
134
135
136
137
138
139
\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}.
tobias's avatar
tobias committed
140
Other options are possible (Sec. \ref{sec:solver-configuration}).
141

tobias's avatar
tobias committed
142
143
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
144
145
of the associated solver every $0.05$ time intervals.  
The first snapshot is set to be written at time $t=0$.
146
147
The above plotter statement creates a plotter file \texttt{MyPlotter} that
allows you to alter the plotter's behaviour.
148
149

Once we are satisfied with the parameters in our \exahype\ specification file, we hand it over to the \exahype\ toolkit:
150
\begin{code}
Jean-Matthieu Gallard's avatar
Jean-Matthieu Gallard committed
151
  ./Toolkit/toolkit.sh Euler2d.exahype
152
153
\end{code}

tobias's avatar
tobias committed
154

155
156
157
158
159
160
161
162
163
164
165
166
\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}

tobias's avatar
tobias committed
167

168
\section{Study the generated code}
169

tobias's avatar
tobias committed
170
171
172
173
174
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:

175
\begin{figure}[h]
tobias's avatar
tobias committed
176
177
178
\begin{center}
  \includegraphics[width=0.5\textwidth]{sketches/simplest-class-architecture.pdf}
\end{center}
179
\caption[Sketch (Inkscape): Class architecture]{Simplest class architecture for ExaHyPE solvers.}
180
\end{figure}
tobias's avatar
tobias committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204

\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. 
di25cox's avatar
di25cox committed
205
This should all be straightforward. 
tobias's avatar
tobias committed
206
207
208
209
210
211
212
213
214
215
216
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.
217
 \exahype\ also provides support for the latter. Note that symbolic identifiers 
218
 may degrade performance.
tobias's avatar
tobias committed
219
220
221
222
223
224
225
226
227
\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 [...]
228
     variables const = 5
tobias's avatar
tobias committed
229
230
231
232
233
234
235
236
237
238
     [...]
   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.
239
240
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).
241
242


tobias's avatar
tobias committed
243
244

\subsection*{Variant B: Working with symbolic identifiers}
245
246
247
\begin{warning}
Working with symbolic identifiers may degrade performance!
\end{warning}
tobias's avatar
tobias committed
248
249
250
251
252
253

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

\begin{code}
  solver [...]
254
    variables const = rho:1,j:3,E:1
tobias's avatar
tobias committed
255
256
257
258
259
260
    [...]
  end solver
\end{code}

\noindent
This specification is exactly equivalent to the previous declaration. 
261
262
263
264
265
266
267
268
269
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:
tobias's avatar
tobias committed
270

271
272
273
274
275
276
277
278
279
\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}
tobias's avatar
tobias committed
280

281
\noindent
282
The wrapper allows us to ``forget'' about the mapping of the unknown onto array
283
284
285
286
287
288
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.

tobias's avatar
tobias committed
289
We note that there are subtle syntactic differences between the plain array
290
291
292
293
294
295
296
297
298
299
300
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.
tobias's avatar
tobias committed
301
  \item All wrapper vector and matrix classes stem from the linear algebra
302
303
304
305
306
  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}
tobias's avatar
tobias committed
307
308


309
310
311
\noindent
From hereon, most of our examples rely on Variant B, i.e.~on the symbolic
names.
312
313
This makes a comparison to text books easier, but be aware that it be less 
efficient than a direct implementation with arrays.
314
315
316
317
318
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
319
wrapper and always is optional. Please note that you might also prefer to wrap
320
321
322
coordinates, e.g., with \texttt{tarch::la::Vector} classes to have Matlab-style
syntax.  
\end{design}
tobias's avatar
tobias committed
323

324
325
326
327
%\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
tobias's avatar
tobias committed
328

di25cox's avatar
di25cox committed
329

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

332
We return to the Euler flow solver:
333
To set up the experiment, we open the implementation file 
ga96nuv's avatar
ga96nuv committed
334
\texttt{MyEulerSolver.cpp} of the class \texttt{MyEulerSolver}. 
tobias's avatar
tobias committed
335
336
337
338
339
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.
di25cox's avatar
di25cox committed
340
The implementation of the initial values might look as follows\footnote{The
341
342
343
exact names of the parameters might change with newer \exahype\ versions and
additional parameters might drop in, but the principle always remains the
same.}:
344
\begin{code}
tobias's avatar
tobias committed
345
void Euler2d::MyEulerSolver::adjustPointSolution(
346
347
348
349
350
  const double *const x,
  const double w,
  const double t,
  const double dt,
  double *Q
351
) {
ga96nuv's avatar
ga96nuv committed
352
 if (tarch::la::equals( t,0.0,1e-15 )) {
353
   Variables vars(Q);
di25cox's avatar
di25cox committed
354
   const double GAMMA = 1.4;
355
356
357
358
359
360
361
362
   
   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;
di25cox's avatar
di25cox committed
363
   }
364
365
}
\end{code}
tobias's avatar
tobias committed
366
367
368



369
\noindent
di25cox's avatar
di25cox committed
370
371
The above routine enables us to specify time dependent solution values.
It further enables us to add source term contributions to the solution values.
372
373
374
375

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
376
algebra subroutines of Peano to check for ``almost equal''  up to machine precision.
377

378
379

We conclude our experimental setup with a compile run and a first test run.
380
381
382
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
383
specified any fluxes yet (a plot should be similar to fig \ref{fig:10_ic}).
384

385
\begin{figure}[h]
386
\begin{center}
di25cox's avatar
di25cox committed
387
  \includegraphics[width=0.75\textwidth]{screenshots/10_initial-condition.png}  
388
\end{center}
389
\caption[Paraview initial conditions of an EulerFlow 2D experiment]{The rest mass density $Q_0$ in the EulerFlow 2D setup. The
390
391
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.}
392
393
\label{fig:10_ic}
\end{figure}
394

395
396
397
398
399
400
401
402
403
404
405
406

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



407
\section{Realising the fluxes}
408

409
410
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
411
\texttt{eigenvalues} in file \texttt{MyEulerSolver.cpp} with code. 
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466


\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}
467
468
469
\begin{warning}
Working with symbolic identifiers may degrade performance!
\end{warning}
470
471
472

Alternatively, we can work with symbolic identifiers if we have specified the
unknowns via \texttt{rho:1,j:3,E:1}:
473
\begin{code}
474
475
void Euler2d::MyEulerSolver::flux(
  const double* const Q,
ga96nuv's avatar
ga96nuv committed
476
  double** F
477
) {
478
479
  ReadOnlyVariables vars(Q);
  Fluxes f(F);
ga96nuv's avatar
ga96nuv committed
480

481
482
483
484
  tarch::la::Matrix<3,3,double> I;
  I = 1, 0, 0,
      0, 1, 0,
      0, 0, 1;
485

486
487
488
  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() );
489

490
491
492
  f.rho ( vars.j()                                 );
  f.j   ( irho * outerDot(vars.j(),vars.j()) + p*I );
  f.E   ( irho * (vars.E() + p) * vars.j()         );
493
494
}
\end{code}
495

496

497
\begin{code}
498
499
500
501
502
void Euler2d::MyEulerSolver::eigenvalues(
  const double* const Q,
  const int normalNonZeroIndex,
  double* lambda
) {
503
504
505
  ReadOnlyVariables vars(Q);
  Variables eigs(lambda);

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

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

513
514
515
516
517
  eigs.rho()=u_n - c;
  eigs.E()  =u_n + c;
  eigs.j(u_n,u_n,u_n);
}
\end{code}
518

519

520
The  implementation of function \texttt{flux} is very straightforward. 
521
522
523
524
525
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
526
527
\texttt{Fluxes} were generated according to the variable list we have specified
in the specification file.
528
529
530
531
532
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.
533

534
535

The pointer \texttt{lambda} appearing in the signature and body of function 
536
537
\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.
538

di25cox's avatar
di25cox committed
539
The normal vectors that are involved in \exahype 's ADER-DG kernels always coincide with the (positively signed) 
540
541
542
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
Anne Reinarz's avatar
Anne Reinarz committed
543
\texttt{normalNonZeroIndex}. 
544
545
546
547


\section{Supplementing boundary conditions}

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

555
\begin{itemize}
556
557
558
559
\item outflow boundary conditions (also sometimes refered to as "null boundary
conditions" or "none boundary conditions"),
\item exact boundary conditions or 
\item reflection symmetries.
560
\end{itemize}
561
562


tobias's avatar
tobias committed
563
\noindent
564
565
566
567
568
569
570
571
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)
    {
Anne Reinarz's avatar
Anne Reinarz committed
572
      ...
573
    }
574
\end{code}
tobias's avatar
tobias committed
575

Anne Reinarz's avatar
Anne Reinarz committed
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
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}



tobias's avatar
tobias committed
592
593
594
595
\begin{figure}[htb]
\begin{center}
  \includegraphics[width=0.5\textwidth]{sketches/cube-names.png}
\end{center}
596
\caption[Rendering (Blender): Peano face indices]{The face indexes as named in the ExaHyPE code}
tobias's avatar
tobias committed
597
598
599
600
601
\label{fig:cube-names}
\end{figure}


\noindent
602
603
604
605
606
607
608
609
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
610
611
612
613
614
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}
615
616
617
618
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}
tobias's avatar
tobias committed
619
620
621
622
623
624
625
626
627
628
629
630
631
632
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.


633

tobias's avatar
tobias committed
634
\begin{figure}[htb]
635
\begin{center}
tobias's avatar
tobias committed
636
  \includegraphics[width=0.5\textwidth]{sketches/boundary-conditions.pdf}
637
\end{center}
638
  \caption[Sketch (Inkscape): 1d boundary conditions]{
tobias's avatar
tobias committed
639
640
641
642
643
644
645
    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}
646
\end{figure}
Tobias Weinzierl's avatar
Tobias Weinzierl committed
647
648


tobias's avatar
tobias committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
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.

667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
\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}.

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

734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
\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];
     }
  }
768
}
769
\end{code}
770

771
772
773
774
775
776
777
778
779
780
781
782
\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}
783

tobias's avatar
tobias committed
784
785
786
787
788
789
790
791
792

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


793
\section{Finite Volumes}\label{sec:FV-intro}
794
795
796
797
798
799
800
801
802

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:

tobias's avatar
tobias committed
803
804
805
806
807
808
\begin{code}
  solver Finite-Volumes MyFVEulerSolver
    variables         = rho:1,j:3,E:1
    patch-size        = 3
    maximum-mesh-size = 0.1
    time-stepping     = global
Tobias Weinzierl's avatar
Tobias Weinzierl committed
809
    kernel            = generic::godunov
tobias's avatar
tobias committed
810
811
812
813
    language          = C
  end solver
\end{code}

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

816
817
818
819
820
\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}
tobias's avatar
tobias committed
821
822

\noindent
823
824
825
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.
tobias's avatar
tobias committed
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
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}

tobias's avatar
tobias committed
842