Chapter 5: Generic ADER-DG and FInite Volume solvers (finished) - improved:...

Chapter 5: Generic ADER-DG and FInite Volume solvers (finished) - improved: subsection 5.6.1 Integration of Exact Boundary Conditions, time-averaged space-time predictor and quadrature, implementation considerations

11_aderdg-user-fluxes.tex

@@ -664,7 +664,70 @@ 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}
+\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. 
 
@@ -705,10 +768,6 @@ void DemoADERDG::DemoSolver::boundaryValues(
 }
 \end{code}
 
-This code snippet assumes you have a function \texttt{setYourExactData}
-which give exact boundary conditions for a point \texttt{x}
-at time \texttt{ti}.
-
 \subsection{Outflow Boundary Conditions}
 As another example, outflow boundary conditions are archieved
 by just copying the fluxes and states.