Project A3 • Adaptive implicit space-time discretization for wave equations
Principal investigators
Prof. Dr. Willy Dörfler
(7/2015 - )
Prof. Dr. Christian Wieners
(7/2015 - )
Project summary
For acoustic, elastic, and electromagnetic wave equations, we developed fully adaptive discretizations in space and time with integrated error control and efficient implicit solution methods with, at least numerically validated (almost) optimal complexity. The system can have still millions of degrees of freedom. Therefore it must be parallelized on several computational cores and solved using a space-time multigrid preconditioner.
One motivation for developing space-time methods is the design of modern computer facilities with an enormous number of processor cores, where the parallel realization of conventional methods becomes inefficient. Since these machines allow a fully implicit space-time approach, new parallel solution techniques are required to solve the huge linear systems.
Linear wave equations
Acoustic waves: find pressure \(p\) and velocity \(\boldsymbol{v}\) with \[\partial_t p - \kappa \operatorname{div} \boldsymbol{v} = 0 \,, \qquad \rho\partial_t \boldsymbol{v} - \nabla p = \boldsymbol{f}\] for permeability \(\kappa\) and density \(\rho\).
Elastic waves: find stress \(\boldsymbol{\sigma}\) and velocity \(\boldsymbol{v}\) with \[\partial_t \boldsymbol{\sigma} - \boldsymbol{C} \boldsymbol{\varepsilon}(\boldsymbol{v}) = \boldsymbol{0} \,, \qquad \rho\partial_t \boldsymbol{v} - \boldsymbol{\operatorname{div}} \boldsymbol{\sigma} = \boldsymbol{f}\] with strain rate \(\boldsymbol{\varepsilon}(\boldsymbol{v}) = \operatorname{sym} (\mathrm{D} \boldsymbol{v})\) and the constitutive relation \[\boldsymbol{C}\boldsymbol{\varepsilon} = 2\mu\boldsymbol{\varepsilon} + \lambda\operatorname{trace}(\boldsymbol{\varepsilon})\textbf{1}\,.\]
We also realized space-time methods for electromagnetic waves [[Fin16]] and for viscoacoustic and viscoelastic waves.
Discretization
We consider three different discretizations.
(symmetric) Space-time discontinuous Petrov-Galerkin discretization (ST-DPG)
The Discontinuous Petrov–Galerkin Method is a methodology to construct robust discretization schemes for linear PDEs. Due to its guaranteed stability independently of the mesh, it is well suited for adaptive applications. Moreover, as a Least-Squares type method, the resulting linear system is symmetric and positive definite which can be exploited by the iterative solver.
Space-time discontinuous Petrov-Galerkin discretization (ST-PGDG)
The discretization uses the trial space $V_h \subset \mathrm{H}^1(0,1;\mathrm{L}_2(\Omega;\mathbb{R}^m))$, which is discontinuous in space but continuous in time. The test space $W_h \subset \mathrm{L}_2((0,T) \times \Omega ;\mathbb{R}^m))$ is discontinuous in space and time. On a tensor product space-time mesh with a fixed mesh in space and a time series $0 = t_0 < \dots < t_N=T$ with a local selection of polynomial degrees in space and time $p_R$ and $q_R$ in every cell, we set for the local test space $W_{h,R}=\mathbb{P}_{q_R-1}((t_{n-1},t_n); \mathbb{R}^m) \otimes \mathbb{P}_{p_R}(K;\mathbb{R}^m)$. Then, the local ansatz spaces $V_{h,R} = V_h|_R$ take the form \begin{align*} V_{h,R} = \{\mathbf{v}_{h,R} \in \mathrm{L}_2(R,\mathbb{R}^m) \colon& \mathbf{v}_{h,R}(t,\mathbf{x}) = \frac{t_n-t}{t_n-t_{n-1}} \mathbf{v}_{h,R}(t_{n-1},\mathbf{x})+ \frac{t_n-t}{t_n-t_{n-1}}\mathbf{w}_{h,R}(t,\mathbf{x}), \\ &\mathbf{v}_{h,R} \in V_h|_{[0,t_{n-1}]}, \mathbf{w}_{h,R} \in W_{h,R}, (t,\mathbf{x}) \in R = I_n \times K\} \end{align*} The global test space is discontinuous in space and time and defined by \begin{align*} W_h = \{\mathbf{w}_h \in \mathrm{L}_2((0,T);\mathrm{H}) \colon \mathbf{w}_{h,R} = \mathbf{w}_{h}|_R \} \end{align*} The global ansatz space $V_h$ is only discontinuous in space but continuous in time. More details are given in [DFWZ19].
Space-time discontinuous Galerkin discretization (ST-DGDG)
Similar to above, this discretization also uses a tensor product space-time mesh. We have each cell $R = (t_{n-1},t_n) \times K$ the polynomial degrees in space and time $(q_R, p_q)$ for its local ansatz and test space. Hence, we have the following local polynomial space $V_{h,R}=\mathbb{P}_{q_R}((t_{n-1},t_n); \mathbb{R}^m) \otimes \mathbb{P}_{p_R}(K;\mathbb R^m) \subset \mathrm{L}_2((0,T) \times \Omega ;\mathbb{R}^m))$. Globally, this yields a discontinuous space in space and time.
Numerical analysis
Using the variational setting, we can show inf-sup stability and a priori bounds for the first two discretizations.
The bilinear form $b_h(\cdot,\cdot)=(L_h \cdot,\cdot)_{0,Q}$ is inf-sup stable in $V_h \times W_h$ with $\beta=C_1T^{-1}$, i.e., \begin{align*} \operatorname{sup}_{\mathbf{w_h} \in W_h\setminus \{0\} } b_h(\mathbf{v_h},\mathbf{w_h}) \geq \beta \|\mathbf{v_h}\|_{V_h},\quad \mathbf{v_h} \in V_h. \end{align*} A direct result is that for given $\mathbf{f}\in L_2(Q;\mathbb{R}^J)$ a unique solution $\mathbf{u_h}\in V_h$ exists solving \begin{align*} (L_h\mathbf{u_h},\mathbf{w_h})_{0,Q}=(f,\mathbf{w_h})_{0,Q}, \quad \mathbf{w_h} \in W_h \end{align*} and satisfying the a priori bound $\|\mathbf{w_h}\|_{V_h} \leq C_2 T\|M_h^{-1}\Pi_h \mathbf{f}\|_W$.
An extensive numerical study in [Ern17] demonstrates that the theoretically predicted convergence rates are achieved in practice for various benchmark problems in one and two spatial dimensions. While the analysis of DPG in [Ern17] and [EW18b] is restricted to homogeneous materials in space, [EW18a] extends the arguments to the case of inhomogeneous material distributions.
Using a mesh-dependent DG-norm $\|\cdot\|_{h,DG}$, we obtain the following convergence result for a sufficiently smooth solution $\mathbf{u} \in H^s(Q,\mathbb{R})$ with $s \geq 1$ and $s \leq \min_{n,K}\{p_{n,K}, q_{n,K}\}+1$, we have \begin{align*} \|\mathbf{u} - \mathbf{u}_h\|_{h,\mathrm{DG}} \leq C h^{s-1/2}\|D^s\mathbf{u}\|_Q + CTh^{-1/2}\|M_h^{-1/2}(M_h-M) \partial_t \mathbf{u}\|_Q \end{align*} For a proof and the construction of a residual error estimator for this setting as well as numerical experiments we refer to [CDW23].
Numerical experiments
Numerical experiment with ST-DPG
Here, we also demonstrate that the space-time DPG method can be applied successfully to an application driven benchmark problem that has been inspired by the Marmousi benchmark. For this benchmark, we use the DPG method with tensor-product polynomials of order \(k\) in each space-time cell and tensor-product elements of order \(k+1\) on the space-time interfaces.
The Marmousi model is a two space dimensions seismic model devised by the Institut Français du Petrole to test seismic data inversion.
In contrast to the \(pq\)-adaptive algorithm used with ST-PGDG, we decided to reduce the computational effort by truncating the space-time mesh.
The next animation shows, that we don't lose any informations by truncating the mesh. We also verify the results by comparing them with ST-PGDG and a classical time stepping scheme.
Numerical experiment with ST-PGDG
The second example illustrates seismic tunnel exploration: An artificially generated surface wave (at \(x_\text{mid}\)) in the tunnel propagates into the solid and the reflected waves are measured in a certain region (marked red in the figure). We solve the acoustic wave equation in two space dimensions using a space-time discontinuous Galerkin discretization. In space we consider discontinuous polynomials over quadrilaterals of degree \(p\) and continuous polynomials of degree \(q\) in time (for details see [DFWZ19]). On the right we see a sketch of the computational domain \(\Omega\). The region of interest is marked in red. Starting with a finite volume discretization \((p=0, q=1)\) we perform four adaptive refinement steps using a dual weighted goal-oriented error estimator.
The animation shows the space-time solution being sliced in time. On the bottom is the polynomial degree with 0 (blue) to 4 (red). The adaptive algorithm is controlled by a goal functional with support in the region of interest.
We compared the uniform and adaptive refinement in the acoustic case observe that the adaptive algorithm saves 70% of the degrees of freedom while achieving the same accuracy compared with uniform refinement. Continue reading.Collapse content.
ref-step
#DoF (effort)
GMRES steps with MG-PC
\(\vartriangle\,E_{\text{ex}}\) in %
uniform refinement
r=1
534 528 (00%)
7
1.85
r=2
2 405 376 (00%)
13
0.22
r=3
6 414 336 (00%)
19
0.49
r=4
13 363 200 (00%)
27
0.25
adaptive refinement
r=0
133 632 (00%)
5
91.01
r=1
291 411 (55%)
7
1.27
r=2
819 279 (34%)
13
0.58
r=3
1 875 753 (29%)
20
0.56
r=4
3 594 969 (27%)
28
0.38
Acoustic wave: Uniform vs. adaptive refinement on \(44\,544 =928 \times 48\) space-time cells distributed on 64 processor cores. The error \(\vartriangle E_\text{ex}(\mathbf{u_h}) = |E(\mathbf{u_h}) -E_\text{ex}|\) of the goal functional is approximately estimated with respect to a linear extrapolation of the uniform results.
The parallel scaling behavior of the parallel multilevel preconditioner is tested for different numbers of processes. On mesh level 4 we have 2 850 816 space-time cells, a linear discretization in space and time results in 34 209 792 DoFs for the acoustic case. The computing time for solving this huge linear system system with the parallel multigrid method scales nearly optimal.
Tunnel Experiment: Strong scaling for \(\sim34\) Mio. DoFs
computed on the bwUniCluster.
We consider a model problem in order to verify our proofs given in [CDW23]. With a smooth initial value and a discontinuous material distribution we obtain the expected convergence rate in the DG-Norm.
Further, a Riemann problem was considered with a discontinuous initial value and constant material parameters. We observe convergence with a reduced rate, since the above smoothness assumption is not satisfied.
J. Ernesti and C. Wieners. A space-time discontinuous Petrov–Galerkin method for acoustic waves. In U. Langer and O. Steinbach, editors, Space-Time Methods: Applications to Partial Differential Equations, volume 25 of Radon Series on Computational and Applied Mathematics, chapter 3, pages 89–116. De Gruyter, Berlin/Boston, September2019.[preprint][bibtex]
R. S. Nejad and C. Wieners. Parallel inelastic heterogeneous multi-scale simulations. In S. Diebels and S. Rjasanow, editors, Multi-scale simulation of composite materials — Results from the project MuSiKo, pages 57–96. Springer, Berlin, February2019.[bibtex]
