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 two different discretizations.

Space-time discontinuous Galerkin discretization (ST-DG)
The discretization uses the trial space \(V_h \subset \mathrm{H}^1(0,T;\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. Continue reading.Collapse content.

On a tensor product space-time mesh with a fixed mesh in space and a time series \(0=t_0\lt t_1\lt\cdots\lt 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 \[ V_{h,R} = \Big\{ \boldsymbol{v}_{h,R}\in \mathrm{L}_2(R;\mathbb{R}^m)\;\colon \boldsymbol{v}_{h,R}(t,\boldsymbol{x}) = \frac{t_n-t}{t_n-t_{n-1}} \boldsymbol{v}_h(t_{n-1},\boldsymbol{x}) + \frac{t-t_{n-1}}{t_n-t_{n-1}} \boldsymbol{w}_{h,R}(t,\boldsymbol{x})\,, \qquad\qquad\qquad\qquad\\ \quad\, \boldsymbol{v}_h \in V_h|_{[0,t_{n-1}]} \,,\ \boldsymbol{w}_{h,R} \in W_{h,R} \,,\ (t,\boldsymbol{x})\in R=I_n\times K \Big\} \,.\] The global test space is discontinuous in space and time and defined by \[W_h = \Big\{ \boldsymbol{w}_h\in \mathrm{L}_2((0,T);H)\;\colon \boldsymbol{w}_{h,R} = \boldsymbol{w}_h|_R \in W_{h,R} \Big\} \,.\] The global ansatz space \(V_h\) is only discontinuous in space but continuous in time. More details are given in [DFWZ18]. Continue reading.Collapse content.

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. Continue reading.Collapse content.

In [Ern17], [EW18b], we construct a non-conforming variant of the DPG method for acoustic waves in space-time and provide an framework to analyze it. Continue reading.Collapse content.

Numerical analysis

Using the variational setting, we can show inf-sup stability and a priori bounds for both discretizations.

The bilinear form \(b_h(\cdot,\cdot) = \big(L_h\cdot,\cdot)_{0,Q}\) is inf-sup stable in \(V_h\times W_h\) with \(\beta = C_1T^{-1}\), i.e., \[\sup_{\boldsymbol{w}_h\in W_h\setminus\{\boldsymbol{0}\}} \frac{b_h(\boldsymbol{v}_h,\boldsymbol{w}_h)}{\|\boldsymbol{w}_h\|_W} \geq \beta\, \|\boldsymbol{v}_h\|_{V_h}\,,\qquad \boldsymbol{v}_h\in V_h\,.\] A direct result is that for given \(\boldsymbol{f}\in\mathrm{L}_2(Q;\mathbb{R}^J)\) a unique solution \(\boldsymbol{u}_h\in V_h\) exists solving \[(L_h\boldsymbol{u}_h,\boldsymbol{w}_h)_{0,Q} = (\boldsymbol{f},\boldsymbol{w}_h)_{0,Q}\,,\qquad \boldsymbol{w}_h\in W_h\] and satisfying the a priori bound \(\|\boldsymbol{u}_h\|_{V_h} \leq C_2T \|M_h^{-1}\Pi_h\,\boldsymbol{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], [EW18b] is restricted to homogeneous materials in space, [EW18a] extends the arguments to the case of inhomogeneous material distributions. Continue reading.Collapse content.

Numerical experiments

Numerical experiment with ST-DG

The first 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 [DFWZ18]). 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.

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-DG, 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-DG and a classical time stepping scheme.

The next animation shows, that we don't lose any informations by truncating the mesh. We also verify the results by comparing them with a classical time stepping scheme (TS-DG).

Outlook

We also consider nonlocal material laws in time (e.g., acoustic and elastic waves with attenuation) and, as long term goal, this will be extended to the numerical homogenization of materials with (random) heterogeneous microstructure. We also consider heterogeneous media on a mesoscopic scale by constructing a local space-time multi-scale basis on finer meshes resulting into a hierarchy of multiscale models.

J. Ernesti and C. Wieners. A space-time discontinuousa 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]

D. Ziegler. A parallel and adaptive space-time discontinuous Galerkin method for visco-elastic and visco-acoustic waves. PhD thesis, Karlsruhe Institute of Technology (KIT), November2019.[bibtex]