Dominic Etienne Charrier committed Sep 01, 2017 1 \chapter{Generic ADER-DG and Finite Volume solvers}  Dominic Etienne Charrier committed Oct 05, 2016 2 \label{chapter:aderdg-user-defined-fluxes}  Tobias Weinzierl committed Jan 29, 2016 3   tobias committed Feb 04, 2017 4   tobias committed Feb 06, 2017 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.  Tobias Weinzierl committed Jan 29, 2016 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 committed Feb 06, 2017 14 15 kernels. Therefore, the present coding strategy is particulary well-suited for rapid  Dominic Etienne Charrier committed Feb 02, 2017 16 prototyping of new solvers.  Tobias Weinzierl committed Jan 29, 2016 17   tobias committed Feb 06, 2017 18 \begin{design}  tobias committed Feb 06, 2017 19 \exahype\ realises the Hollywood principle: Don't call us, we call you!''  tobias committed Feb 06, 2017 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.  SvenK committed May 24, 2017 46 \section{Establishing the data structures}\label{sec:euler-starting}  tobias committed Feb 06, 2017 47 48 49  We follow the Euler equations in the form  SvenK committed Feb 22, 2017 50 \label{eq:euler}  Dominic Etienne Charrier committed Jan 29, 2016 51 52 53 54 \frac{\partial}{\partial t} \textbf{Q} + \boldsymbol{\nabla}\cdot\textbf{F}(\textbf{Q}) = 0  Tobias Weinzierl committed Jan 29, 2016 55 \quad \mbox{with} \quad  Dominic Etienne Charrier committed Jan 29, 2016 56 \textbf{Q} = \begin{pmatrix}  Dominic Etienne Charrier committed Feb 02, 2017 57 \rho\\\textbf{j}\\\ E  Dominic Etienne Charrier committed Jan 29, 2016 58 \end{pmatrix}  Tobias Weinzierl committed Jan 29, 2016 59 \quad \mbox{and} \quad  Dominic Etienne Charrier committed Jan 29, 2016 60 61 \textbf{F} = \begin{pmatrix}  Dominic Etienne Charrier committed Feb 02, 2017 62 63 64 \textbf{j}\\ \frac{1}{\rho}\textbf{j}\otimes\textbf{j} + p I \\ \frac{1}{\rho}\textbf{j}\,(E + p)  Tobias Weinzierl committed Jan 29, 2016 65 \end{pmatrix}  Sven committed Feb 10, 2017 66   tobias committed Feb 27, 2017 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 committed Feb 06, 2017 72 $\rho$ denotes the mass density,  Dominic Etienne Charrier committed Feb 02, 2017 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 committed Feb 06, 2017 76 Introducing the adiabatic index $\gamma$, the fluid pressure is defined as  Sven committed Feb 10, 2017 77 78 %  Dominic Etienne Charrier committed Feb 02, 2017 79 p = (\gamma-1)\,\left(E - \frac{1}{2}\textbf{j}^2/\rho\right).  Sven committed Feb 10, 2017 80   tobias committed Feb 27, 2017 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}.  di25cox committed Feb 02, 2016 89 A corresponding specification file \texttt{euler-2d.exahype} for this  Tobias Weinzierl committed Jan 29, 2016 90 91 92 setup is \begin{code} exahype-project Euler2d  Tobias Weinzierl committed Jan 29, 2016 93   Anne Reinarz committed Aug 09, 2018 94  peano-kernel-path const = ./Peano/  95  exahype-path const = ./ExaHyPE  96  output-directory const = ./ApplicationExamples/EulerFlow  tobias committed Jul 07, 2016 97   Tobias Weinzierl committed Jan 29, 2016 98  computational-domain  99  dimension const = 2  tobias committed Jul 07, 2016 100  width = 1.0, 1.0  Tobias Weinzierl committed Jan 29, 2016 101  offset = 0.0, 0.0  di25cox committed Feb 05, 2016 102  end-time = 0.4  Tobias Weinzierl committed Jan 29, 2016 103  end computational-domain  Tobias Weinzierl committed Jan 29, 2016 104   Tobias Weinzierl committed Jan 29, 2016 105  solver ADER-DG MyEulerSolver  106 107  variables const = 5 order const = 3  ga96nuv committed Aug 02, 2016 108 109  maximum-mesh-size = 0.1 time-stepping = global  Dominic Etienne Charrier committed Sep 21, 2017 110 111 112  type const = nonlinear terms const = flux optimisation const = generic  113  language const = C  tobias committed Jul 07, 2016 114   Anne Reinarz committed Aug 09, 2018 115 116 117 118 119  plot vtu::Cartesian::cells::ascii EulerWriter variables const = 5 time = 0.0 repeat = 0.5E-1 output = ./variables  Tobias Weinzierl committed Jan 29, 2016 120 121  end plot end solver  Tobias Weinzierl committed Jan 29, 2016 122 end exahype-project  Tobias Weinzierl committed Jan 29, 2016 123 124 125 \end{code} \noindent  SvenK committed Feb 22, 2017 126 The spec. file sets some paths in the preamble before it specifies the computational  Tobias Weinzierl committed Jan 29, 2016 127 128 domain and the simulation time frame in the environment \texttt{computational-domain}.  di25cox committed Feb 05, 2016 129 130  In the next lines, a solver of type \texttt{ADER-DG} is specified and assigned the name  Tobias Weinzierl committed Mar 05, 2016 131 \texttt{MyEulerSolver}. The kernel type of the solver is set to  Dominic Etienne Charrier committed Sep 21, 2017 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 committed Sep 22, 2017 140 Other options are possible (Sec. \ref{sec:solver-configuration}).  di25cox committed Feb 05, 2016 141   tobias committed Feb 06, 2017 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  di25cox committed Feb 05, 2016 144 145 of the associated solver every $0.05$ time intervals. The first snapshot is set to be written at time $t=0$.  tobias committed Jan 09, 2017 146 147 The above plotter statement creates a plotter file \texttt{MyPlotter} that allows you to alter the plotter's behaviour.  di25cox committed Feb 05, 2016 148 149  Once we are satisfied with the parameters in our \exahype\ specification file, we hand it over to the \exahype\ toolkit:  Tobias Weinzierl committed Jan 29, 2016 150 \begin{code}  Jean-Matthieu Gallard committed Sep 17, 2018 151  ./Toolkit/toolkit.sh Euler2d.exahype  Tobias Weinzierl committed Jan 29, 2016 152 153 \end{code}  tobias committed Feb 06, 2017 154   tobias committed Feb 27, 2017 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 committed Feb 06, 2017 167   Tobias Weinzierl committed Jan 29, 2016 168 \section{Study the generated code}  Tobias Weinzierl committed Jan 29, 2016 169   tobias committed Feb 06, 2017 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:  Sven committed Feb 10, 2017 175 \begin{figure}[h]  tobias committed Feb 06, 2017 176 177 178 \begin{center} \includegraphics[width=0.5\textwidth]{sketches/simplest-class-architecture.pdf} \end{center}  Sven K committed Nov 07, 2017 179 \caption[Sketch (Inkscape): Class architecture]{Simplest class architecture for ExaHyPE solvers.}  Sven committed Feb 10, 2017 180 \end{figure}  tobias committed Feb 06, 2017 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 committed Feb 05, 2016 205 This should all be straightforward.  tobias committed Feb 06, 2017 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.  Jean-Matthieu Gallard committed Jun 20, 2017 217  \exahype\ also provides support for the latter. Note that symbolic identifiers  Dominic Etienne Charrier committed Sep 01, 2017 218  may degrade performance.  tobias committed Feb 06, 2017 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 committed Feb 06, 2017 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.  SvenK committed Feb 22, 2017 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).  tobias committed Feb 06, 2017 241 242   tobias committed Feb 06, 2017 243 244  \subsection*{Variant B: Working with symbolic identifiers}  Dominic Etienne Charrier committed Sep 01, 2017 245 246 247 \begin{warning} Working with symbolic identifiers may degrade performance! \end{warning}  tobias committed Feb 06, 2017 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 committed Feb 06, 2017 255 256 257 258 259 260  [...] end solver \end{code} \noindent This specification is exactly equivalent to the previous declaration.  tobias committed Feb 06, 2017 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 committed Feb 06, 2017 270   tobias committed Feb 06, 2017 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 committed Feb 06, 2017 280   tobias committed Feb 06, 2017 281 \noindent  SvenK committed Feb 22, 2017 282 The wrapper allows us to forget'' about the mapping of the unknown onto array  tobias committed Feb 06, 2017 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 committed Feb 28, 2017 289 We note that there are subtle syntactic differences between the plain array  tobias committed Feb 06, 2017 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 committed Feb 28, 2017 301  \item All wrapper vector and matrix classes stem from the linear algebra  tobias committed Feb 06, 2017 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 committed Feb 06, 2017 307 308   tobias committed Feb 06, 2017 309 310 311 \noindent From hereon, most of our examples rely on Variant B, i.e.~on the symbolic names.  Jean-Matthieu Gallard committed Jun 20, 2017 312 313 This makes a comparison to text books easier, but be aware that it be less efficient than a direct implementation with arrays.  tobias committed Feb 06, 2017 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  SvenK committed Feb 22, 2017 319 wrapper and always is optional. Please note that you might also prefer to wrap  tobias committed Feb 06, 2017 320 321 322 coordinates, e.g., with \texttt{tarch::la::Vector} classes to have Matlab-style syntax. \end{design}  tobias committed Feb 06, 2017 323   SvenK committed Feb 22, 2017 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 committed Feb 06, 2017 328   di25cox committed Feb 05, 2016 329   SvenK committed Feb 22, 2017 330 \section{Setting up the experiment}\label{section:setup}  Tobias Weinzierl committed Jan 29, 2016 331   tobias committed Feb 06, 2017 332 We return to the Euler flow solver:  Tobias Weinzierl committed Jan 29, 2016 333 To set up the experiment, we open the implementation file  ga96nuv committed Aug 02, 2016 334 \texttt{MyEulerSolver.cpp} of the class \texttt{MyEulerSolver}.  tobias committed Jun 15, 2018 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 committed Feb 05, 2016 340 The implementation of the initial values might look as follows\footnote{The  Tobias Weinzierl committed Jan 29, 2016 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.}:  Tobias Weinzierl committed Jan 29, 2016 344 \begin{code}  tobias committed Jun 15, 2018 345 void Euler2d::MyEulerSolver::adjustPointSolution(  Angelika Schwarz committed Apr 27, 2016 346 347 348 349 350  const double *const x, const double w, const double t, const double dt, double *Q  di25cox committed Feb 02, 2016 351 ) {  ga96nuv committed Aug 02, 2016 352  if (tarch::la::equals( t,0.0,1e-15 )) {  Dominic Etienne Charrier committed Feb 02, 2017 353  Variables vars(Q);  di25cox committed Mar 07, 2016 354  const double GAMMA = 1.4;  Dominic Etienne Charrier committed Feb 02, 2017 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 committed Mar 07, 2016 363  }  Tobias Weinzierl committed Jan 29, 2016 364 365 } \end{code}  tobias committed Feb 06, 2017 366 367 368   Tobias Weinzierl committed Jan 29, 2016 369 \noindent  di25cox committed Mar 07, 2016 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.  tobias committed Feb 06, 2017 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  SvenK committed Feb 22, 2017 376 algebra subroutines of Peano to check for almost equal'' up to machine precision.  tobias committed Feb 06, 2017 377   Tobias Weinzierl committed Jan 29, 2016 378 379  We conclude our experimental setup with a compile run and a first test run.  Tobias Weinzierl committed Jan 29, 2016 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  Sven committed Feb 10, 2017 383 specified any fluxes yet (a plot should be similar to fig \ref{fig:10_ic}).  Tobias Weinzierl committed Jan 29, 2016 384   Sven committed Feb 10, 2017 385 \begin{figure}[h]  Tobias Weinzierl committed Jan 29, 2016 386 \begin{center}  di25cox committed Feb 05, 2016 387  \includegraphics[width=0.75\textwidth]{screenshots/10_initial-condition.png}  Tobias Weinzierl committed Jan 29, 2016 388 \end{center}  Sven K committed Nov 03, 2017 389 \caption[Paraview initial conditions of an EulerFlow 2D experiment]{The rest mass density $Q_0$ in the EulerFlow 2D setup. The  SvenK committed Feb 22, 2017 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.}  Sven committed Feb 10, 2017 392 393 \label{fig:10_ic} \end{figure}  Tobias Weinzierl committed Jan 29, 2016 394   tobias committed Feb 06, 2017 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}  Tobias Weinzierl committed Jan 29, 2016 407 \section{Realising the fluxes}  Tobias Weinzierl committed Jan 29, 2016 408   SvenK committed Feb 22, 2017 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  di25cox committed Feb 05, 2016 411 \texttt{eigenvalues} in file \texttt{MyEulerSolver.cpp} with code.  tobias committed Feb 06, 2017 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}  Dominic Etienne Charrier committed Sep 01, 2017 467 468 469 \begin{warning} Working with symbolic identifiers may degrade performance! \end{warning}  tobias committed Feb 06, 2017 470 471 472  Alternatively, we can work with symbolic identifiers if we have specified the unknowns via \texttt{rho:1,j:3,E:1}:  Tobias Weinzierl committed Jan 29, 2016 473 \begin{code}  di25cox committed Feb 02, 2016 474 475 void Euler2d::MyEulerSolver::flux( const double* const Q,  ga96nuv committed Aug 02, 2016 476  double** F  di25cox committed Feb 02, 2016 477 ) {  Dominic Etienne Charrier committed Feb 02, 2017 478 479  ReadOnlyVariables vars(Q); Fluxes f(F);  ga96nuv committed Aug 02, 2016 480   Dominic Etienne Charrier committed Feb 02, 2017 481 482 483 484  tarch::la::Matrix<3,3,double> I; I = 1, 0, 0, 0, 1, 0, 0, 0, 1;  Tobias Weinzierl committed Jan 29, 2016 485   Dominic Etienne Charrier committed Feb 02, 2017 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() );  Tobias Weinzierl committed Jan 29, 2016 489   Dominic Etienne Charrier committed Feb 02, 2017 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() );  Tobias Weinzierl committed Jan 29, 2016 493 494 } \end{code}  Tobias Weinzierl committed Jan 29, 2016 495   di25cox committed Feb 02, 2016 496   Tobias Weinzierl committed Jan 29, 2016 497 \begin{code}  di25cox committed Feb 02, 2016 498 499 500 501 502 void Euler2d::MyEulerSolver::eigenvalues( const double* const Q, const int normalNonZeroIndex, double* lambda ) {  Dominic Etienne Charrier committed Feb 02, 2017 503 504 505  ReadOnlyVariables vars(Q); Variables eigs(lambda);  di25cox committed Feb 02, 2016 506  const double GAMMA = 1.4;  Dominic Etienne Charrier committed Feb 02, 2017 507 508  const double irho = 1./vars.rho(); const double p = (GAMMA-1) * (vars.E() - 0.5 * irho * vars.j()*vars.j() );  Tobias Weinzierl committed Jan 29, 2016 509   Dominic Etienne Charrier committed Feb 02, 2017 510 511  double u_n = vars.j(normalNonZeroIndex) * irho; double c = std::sqrt(GAMMA * p * irho);  Tobias Weinzierl committed Jan 29, 2016 512   Dominic Etienne Charrier committed Feb 02, 2017 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}  Tobias Weinzierl committed Jan 29, 2016 518   di25cox committed Feb 02, 2016 519   Dominic Etienne Charrier committed Feb 02, 2017 520 The implementation of function \texttt{flux} is very straightforward.  tobias committed Feb 06, 2017 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  Dominic Etienne Charrier committed Feb 02, 2017 526 527 \texttt{Fluxes} were generated according to the variable list we have specified in the specification file.  tobias committed Feb 06, 2017 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.  Tobias Weinzierl committed Jan 29, 2016 533   Dominic Etienne Charrier committed Feb 02, 2017 534 535  The pointer \texttt{lambda} appearing in the signature and body of function  tobias committed Feb 06, 2017 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.  di25cox committed Feb 05, 2016 538   di25cox committed Feb 07, 2016 539 The normal vectors that are involved in \exahype 's ADER-DG kernels always coincide with the (positively signed)  di25cox committed Feb 05, 2016 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 committed Aug 09, 2018 543 \texttt{normalNonZeroIndex}.  tobias committed Feb 06, 2017 544 545 546 547  \section{Supplementing boundary conditions}  SvenK committed Feb 22, 2017 548 \exahype\ offers a generic API to implement any kind of \emph{local}  tobias committed Feb 28, 2017 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.  tobias committed Feb 06, 2017 552 553 554 That is, you can straightforwardly\footnote{In this guidebook, this implies that we stick to those guys only.} implement  SvenK committed Feb 01, 2017 555 \begin{itemize}  tobias committed Feb 06, 2017 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.  SvenK committed Feb 01, 2017 560 \end{itemize}  tobias committed Feb 06, 2017 561 562   tobias committed Feb 15, 2017 563 \noindent  SvenK committed Feb 01, 2017 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 committed Aug 09, 2018 572  ...  tobias committed Feb 06, 2017 573  }  SvenK committed Feb 01, 2017 574 \end{code}  tobias committed Feb 15, 2017 575   Anne Reinarz committed Aug 09, 2018 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 committed Feb 15, 2017 592 593 594 595 \begin{figure}[htb] \begin{center} \includegraphics[width=0.5\textwidth]{sketches/cube-names.png} \end{center}  Sven K committed Nov 07, 2017 596 \caption[Rendering (Blender): Peano face indices]{The face indexes as named in the ExaHyPE code}  tobias committed Feb 15, 2017 597 598 599 600 601 \label{fig:cube-names} \end{figure} \noindent  tobias committed Feb 06, 2017 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  SvenK committed Feb 01, 2017 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}  SvenK committed Feb 09, 2017 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 committed Feb 15, 2017 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.  SvenK committed Feb 01, 2017 633   tobias committed Feb 15, 2017 634 \begin{figure}[htb]  tobias committed Feb 06, 2017 635 \begin{center}  tobias committed Feb 15, 2017 636  \includegraphics[width=0.5\textwidth]{sketches/boundary-conditions.pdf}  tobias committed Feb 06, 2017 637 \end{center}  Sven K committed Nov 07, 2017 638  \caption[Sketch (Inkscape): 1d boundary conditions]{  tobias committed Feb 15, 2017 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}  Sven committed Feb 10, 2017 646 \end{figure}  Tobias Weinzierl committed Aug 08, 2016 647 648   tobias committed Feb 15, 2017 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.  Olga Glogowska committed Oct 13, 2019 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: \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} } 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$: \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}) 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: \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 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: \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} 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: \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} 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}.  Sven Köppel committed Feb 24, 2017 731 732 We supply here a small example how to correctly integrate imposed Boundary Conditions in the ADER-DG scheme.  SvenK committed Feb 22, 2017 733   Sven Köppel committed Feb 24, 2017 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]; } }  tobias committed Feb 27, 2017 768 }  Sven Köppel committed Feb 24, 2017 769 \end{code}  SvenK committed Feb 22, 2017 770   Sven Köppel committed Feb 24, 2017 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}  SvenK committed Feb 22, 2017 783   tobias committed Feb 27, 2017 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}  SvenK committed May 24, 2017 793 \section{Finite Volumes}\label{sec:FV-intro}  tobias committed Feb 06, 2017 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 committed Feb 06, 2017 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 committed Jun 23, 2017 809  kernel = generic::godunov  tobias committed Feb 06, 2017 810 811 812 813  language = C end solver \end{code}  Sven K committed Nov 03, 2017 814 \todo[inline]{This needs an update.}  tobias committed Sep 22, 2017 815   tobias committed Feb 06, 2017 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 committed Feb 06, 2017 821 822  \noindent  tobias committed Feb 06, 2017 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 committed Feb 15, 2017 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 committed Feb 27, 2017 842