Log In Sign Up

Efficient and scalable data structures and algorithms for goal-oriented adaptivity of space-time FEM codes

by   Uwe Köcher, et al.
HSU Hamburg

The cost- and memory-efficient numerical simulation of coupled volume-based multi-physics problems like flow, transport, wave propagation and others remains a challenging task with finite element method (FEM) approaches. Goal-oriented space and time adaptive methods derived from the dual weighted residual (DWR) method appear to be a shiny key technology to generate optimal space-time meshes to minimise costs. Current implementations for challenging problems of numerical screening tools including the DWR technology broadly suffer in their extensibility to other problems, in high memory consumption or in missing system solver technologies. This work contributes to the efficient embedding of DWR space-time adaptive methods into numerical screening tools for challenging problems of physically relevance with a new approach of flexible data structures and algorithms on them, a modularised and complete implementation as well as illustrative examples to show the performance and efficiency.


page 1

page 2

page 3

page 4


Flexible goal-oriented space-time adaptivity for coupled Stokes flow and convection-dominated transport

In this work, a flexible, fully space-time adaptive finite element appro...

Space-Time-Parallel Poroelasticity Simulation

The accurate, reliable and efficient numerical approximation of multi-ph...

On the implementation of an adaptive multirate framework for coupled transport and flow

In this work, a multirate in time approach resolving the different time ...

GPU Acceleration of Hermite Methods for the Simulation of Wave Propagation

The Hermite methods of Goodrich, Hagstrom, and Lorenz (2006) use Hermite...

1 Motivation and significance

Figure 1: General workflow of the space-time adaptive numerical simulation tools DTM++.Project/dwr

-* using the dual weighted residuals (DWR) technique for the error estimation (left), an exemplary forward and backward time marching step

on +-dimensional space-time slabs in the DWR loop (middle) and an exemplary space-time adaptation update of the spatial triangulation and the temporal subinterval to prepare the DWR loop + (right).

1.1 Introduction

The accurate, reliable and efficient numerical approximation of flow in heterogeneous deformable porous media include several coupled multi-physics phenomena, such as multi-phase diffusion and convection-dominated transport with underlying chemical reactions, deformation and poroelastic wave propagation as well as fluid-structure interactions in a more general sense, and is of fundamental importance in environmental, civil, energy, biomedical and many other engineering fields to yield a cost-effective numerical simulation screening tool to support the still challenging and fundamental research of understanding such multi-physics phenomena. Considerable representatives are coupled problems based on the Navier-Stokes equations for incompressible flow, the Euler equations for compressible flow, the convection-diffusion-reaction equations and the Biot-Allard equations for coupled flow with poroelastic wave propagation, which are characterised by partial differential equations of nonlinear and instationary character. Space-time finite element method (FEM) approximations offer appreciable benefits over finite difference and finite volume methods such as the flexibility with which they can accommodate discontinuities in the model, material parameters and boundary conditions as well as for a priori and a posteriori error estimation to establish optimal computational space-time grids in a self-adaptive way. The dual weighted residual (DWR) method for goal-oriented adaptivity was introduced by Becker and Rannacher (

Becker1998 ; Becker2001 ) and further studied; cf. Endtmayer2017 ; Bruchhaeuser2018 ; Bruchhaeuser2018a ; Bangerth2003 and references therein. The DWR method facilitates to find optimal spatial mesh adaptations, optimal temporal mesh adaptations, as well as optimal local space-time polynomial degrees, tuning parameters and physical or numerical models, such that the overall computational costs are minimised for reaching a goal in a target quantity of interest.

A target quantity is a cost, error or energy functional and depends on an user-chosen quantity of interest. Examples for a quantity of interest are the stress derived from the primal displacement variables or the drag coefficient derived from the primal fluid velocity variables of the underlying problem. The goal of the DWR method is to satisfy a guaranteed error bound in the target quantity of interest instead of satisfying a classical error bound of the primal solution variable in a standard norm.

Unfortunately, the DWR method still has challenging drawbacks compared to standard a posteriori error based adaptive methods. It needs to solve the primal or forward problem and an auxiliary dual or adjoint problem of higher approximation quality in each loop of an optimisation problem, it needs variational space-time discretisations for the primal and dual problem, which yield problems on +-dimensional domains, and there is a lack in efficient data structures, algorithms as well as (non-)linear system solver and preconditioning technologies for sophisticated problems of physical interest, which will be at least partially resolved by this work.

Fortunately, the DWR method works in general situations, in which the problem itself and the quantity of interest are of nonlinear fashion, and it does not rely on generally unknown assumptions of the initial space-time mesh. The key idea of the DWR method is to embed the given problem of finding solutions of a partial differential equation into the framework of optimal control to estimate the functional error of , derive computable a posteriori error estimates with an approximation of the corresponding dual solution of an auxiliary problem, execute space-time mesh (and other) adaptations based on the information given by and loop until the goal in the target quantity is reached as it is outlined by Fig. 1. Consider to solve the constrained optimisation problem given by


for a (cost, error or energy) functional depending on , under the constraint of finding from the variational problem

which is set up from the partial differential equation problem in a standard way. and are corresponding variational trial and test space-time functional spaces such that primal and dual solutions exist at least locally, are unique and stable, but the latter ones can not be guaranteed for arbitrary nonlinear problems without restrictions. The calculus of variations theory associates a corresponding Lagrangian functional to (1),


with as a Lagrange multiplier or, equivalently, as the adjoint or dual solution variable, for the global optimisation problem (1). The primal solution is the first component of a stationary point of whereas the second component yields the corresponding dual solution ; cf. Becker1998 ; Becker2001 ; Bangerth2003 for details. Due to the mathematically challenging theory of even local existence of solutions in general nonlinear settings, we restrict ourselves further, without the loss of generality and the impact of our software, to a linear prototype model for instationary transport in heterogeneous porous media.

1.2 Exact Scientific Problem solved by the Software

The DTM++.Project/dwr-diffusion frontend simulation tool for the established finite element analysis library deal.II dealIIReferenceV90 yields a reference implementation of the prototype model for instationary transport in heterogeneous porous media. We consider the goal-oriented approximation of a target quantity subject to find from the diffusion equation


in the space-time domain with (open and bounded), dimension , and , , and equipped with appropriate initial and boundary conditions

with the partition of the boundary and .

A target quantity of interest derived from the primal solution of (3) might represent a local pollution concentration under diffusive transport. The coefficient functions and are the mass density and the permeability of the medium, respectively, and the right hand side function represents volume forces acting on the interior domain such as source or sink terms. In (3), the differential operators , and denote the partial derivative for the time variable , the divergence and the gradient for the space variables , respectively, and,

denotes the outward facing normal vector, in standard notation. The initial-boundary value problem is closed by the choice of an appropriate initial value function

and a boundary value function acting on the Dirichlet part of the boundary partition . Additionally, there is the choice for an appropriate inhomogeneous Neumann boundary type condition function if .

The weak variational primal problem, for brevity consider to find or put without loss of generality of the following discussion, reads as: Seek , , satisfying , , such that


with the left-hand side bilinear form ,

and the right-hand side linear form ,

Remark that , with having inhomogeneous Dirichlet boundary conditions , can be generated in a standard way as it will be explained in the following algorithm. Briefly, the weak variational dual problem results from the stationary condition of (2) given by


for and ; using the Euler-Lagrange method of constrained optimisation, cf. for details (Bangerth2003, , Sec. 6.1), Bruchhaeuser2018 . The dual problem for a general nonlinear goal functional reads as: Find , , satisfying , , such that


for all , , and with the left-hand side adjoint bilinear form ,

which is derived from integration by parts in time of of (4) and put this into the first equation of (5). Remark that the choice of trial and test functions for the dual problem here yields traces on and . To estimate the functional error of and to derive computable a posteriori error estimates for the space-time mesh adaptivity, we need to discretise the primal problem (4) and the dual problem (6).

Figure 2: Exemplary illustration of fully discrete primal and dual solutions for a fixed observation point in space in a DWR loop over time; a piecewise constant discontinuous in time primal solution of (7) (left), and, a piecewise linear continuous in time dual solution of (8) for a specific goal functional (right), and the drawn arrows represent the implemented path of the natural time integration on a (+)-slab approach.

The primal and dual problem are discretised with appropriate space-time variational methods, i.e. here a piecewise constant discontinuous Galerkin method in time for the primal problem and a piecewise linear continuous Petrov-Galerkin method in time for the dual problem (Fig. 2), combined with a piecewise polynomial continuous Ritz-Galerkin (FEM) method in space; cf. for details Koecher2015 and Bruchhaeuser2018 ; Bangerth2003 . The space-time cylinder is divided into non-overlapping space-time slabs , with the partition of as , , , , for the -th loop. On each consider a not necessarily conforming partition of into non-overlapping elements , denoted as geometrical quadrilaterals for or hexahedrons for . The fully discrete primal and dual solutions are represented by


In (7)-(8), and

denote the space-time degrees of freedom (DoF) of the primal and dual problem,

and denote the global spatial trial basis functions and and denote the global temporal trial basis functions (Fig. 3). The basis functions are defined on a reference space-time slab , e.g. as piecewise Lagrange polynomials, and mapped appropriately to ; cf. Koecher2015 for details.

Figure 3: Illustration of temporal reference trial basis functions on the reference subinterval for the primal solution (7) (left) and the dual solution (8) (right), which are mapped with , , from to .

The software implements the goal functional ,


for all , with denoting the fully discrete primal solution of the -th DWR loop and denoting the analytic solution (only possible for academic test problems) or an approximation of the exact primal solution generated on a finer space-time mesh or with higher-order approaches, which aims to reach


for an absolute or relative tolerance tol criterion, on the space-time control volume ; cf. Fig. 4.

Algorithm (dwr-diffusion). Loop until the goal (10) is reached; cf. Fig. 1. Solve the primal problem: Find the coefficient vector , , from


for , by marching forwardly in time through the slabs. The mass and stiffness matrix of the primal problem on the slab are denoted as and , respectively, the right hand side assemblies and correspond to volume forcing and inhomogeneous Neumann boundary terms and the vector

is the interpolation of the initial value function

from (3) for or the fully discrete primal solution of the previous (-) slab for on the primal finite element space of the current slab . Each system is modified such that the Dirichlet boundary conditions from (3) are applied strongly. To set up the dual problem with the goal (9), compute the contribution of the value (10) with a post-processing on each slab . Solve the dual problem: Find the coefficient vector , , from


for , by marching backwardly in time through the slabs. The mass and stiffness matrix of the dual problem on the slab are denoted as and , respectively, the right hand side assembly vector corresponds to the goal defined by (9) and the vector is the interpolation of the homogeneous initial value function (6) for the chosen goal (9) for or the fully discrete dual solution of the next (+) slab for on the dual finite element space of the current slab . The same type of boundary colorisation (either Dirichlet or Neumann type) is used for the dual problem but with homogeneous boundary value functions, even in the case of inhomogeneous primal boundary conditions. Each system is modified such that homogeneous Dirichlet boundary conditions on are applied strongly to the dual solution (8). The numerically approximated a posteriori space-time error estimate


is derived from the primal and dual solutions and , respectively, and each localised , , , is stored independently; cf. for details Bruchhaeuser2018a . The space-time mesh refinement update is implemented as follows: mark a space-time slab for refinement in time for which the corresponding value of (13) belongs to the top fraction of largest values, then, on each slab , mark a mesh cell for -dimensional isotropic refinement in space for which the corresponding value of (13) belongs to the top fraction or , for a slab that is not or is marked for time refinement, of largest values, with , then execute the spatial refinement and finally execute the temporal refinement.

Figure 4: Illustration of the implemented local (+)-dimensional space-time control volume , with and , for Eq. (9) and Sec. 3.

1.3 Contribution of the Software to Scientific Discovery

Modules of the DTM++.Project, including higher-order in time discretisations with sophisticated solver technologies and distributed-memory parallel implementations, already have contributed to the process of scientific discovery for acoustic, elastic and coupled wave propagation (xwave), mass conservative transport (meat) and coupled deformation with transport (biot) for instance; cf. Koecher2015 ; Bause2017 and references therein. The stationary predecessor of dwr-diffusion is used in Bruchhaeuser2018 and the successor is used in Bruchhaeuser2018a for stabilised instationary convection-dominated transport problems. The contributed software shares the modularity and flexibility of all DTM++.Project modules, it provides a freely available open-source framework for efficiently solving problems with goal-oriented adaptivity and will support further scientific discovery with numerical simulations of challenging problems with physical relevance.

1.4 Basic Setup and Usage of the Software

The user needs to prepare a deal.II v9.0 toolchain to compile and run the software, this can be done on any major platform by using candi; open a terminal and type in

it clone

dandi && ./

and follow the instructions. Then download, compile and run

it clone

d dwr-diffusion &&make . && make release


with the experimental setting given in Sec. 3.

1.5 Related work

The origin of the code is related to the step-14 tutorial code of deal.II dealIIReferenceV90 for the Laplace equation and other DTM++.Project modules Koecher2015 . Additionally, the code gallery of dealIIReferenceV90 provides some implementations for stationary elastoplasticity problems using DWR-based goal-oriented mesh adaptivity and for the instationary convection-diffusion-reaction equation without adaptivity for instance. All of the listed related implementations can gain advantages from our work since those are mostly specific implementations with partly lacking documentation and, more importantly, do not document their data structures and algorithms in an application-free way for allowing the reuse in other frameworks straight forwardly.

2 Software description

The software DTM++.Project/dwr-diffusion implements the introduced algorithm of Sec. 1.2 in an object-oriented way and is shipped with an extensive in-source documentation. DTM++ modules are written in the C++.17 language and make use of the new language features since C++.11 such as reference counted dynamic memory allocation, range-based loops, the auto specifier, strongly typed enum classes and others.

Figure 5: Pictorial DTM++.Project software module architecture overview.
Figure 6: Pictorial overview of the DTM++.Project data storage management. A list entry corresponds to a slab and stores a small array of shared pointer references. The large data Vector is stored independently and can be serialised.
Figure 7: Pictorial overview of the elements to handle a computational slab implemented by the DTM++::Grid_DWR class. Remark that loop based slab data structures is not stored for efficiency reasons.

2.1 Software Architecture

The general workflow to use the dwr-diffusion solver module is illustrated by Fig. 5, it highly rely on user-driven inputs by simple text-based parameter input files for a wide range of similar numerical examples. To keep the code clean for the intended user group, the implementation of the presented algorithm from Sec. 1.2 is done with a procedural-structured single central solver class template Diffusion_DWR<dim> : DTM::Problem for slightly different time discretisations of the primal and dual problem. All assemblers, such as for the matrices, right-hand side vectors and even the space-time error estimate , are using the deal.II workstream technology for thread-parallel and further MPI+X-parallel simulations. Thereby, the useful and efficient auxiliary classes of the DTM++ suite are provided independently, such that the input parameter handling and data output handling must not be included into the solver class. An overview on the general DTM++ software framework for PDE solver frameworks can be found in Koecher2015 and by the in-source doxygen documentation of the code.

2.2 Software Key Technologies

One of the key technologies of this work is the efficient handling of the space-time DoF data storage management as given in Fig. 6 which allow to efficiently iterate forwardly and backwardly to compute the primal and dual solutions as well as the error estimate. Another important key technology is a list approach of the new DTM++::Grid_DWR class to handle (+)-dimensional space-time computational slabs and grids as given in Fig. 7 and allows for inexpensive local time refinements.

3 Illustrative Examples

To demonstrate the functionalities, an analytic solution , which mimics a rotating cone with a time-dependent height,


with and , and, , , for and , , for , , , and, scalars , , is approximated on a two-dimensional L-shaped domain (Fig. 4) for with the minimisation goal from (9) on a time-dependent space-time control volume , with , and .

Figure 8: Solution profiles and dual solution contour lines for Sec. 3.

The setup uses , , . The functions , , and of (3) are derived from the analytic solution (14), the coefficient functions are set to and , and, the scalars are set to and . The space-time triangulation for is composed of slabs having 3 spatial cells each. The goal tolerance is . Solution profiles for and contour lines of are illustrated in Fig. 8. Selected convergence and error estimator results as well as a quality measurement , where is optimal, are given by Tab. 1.

1 5 3 6.070e-02 1.994e-02 0.33
2 6 12 2.634e-02 2.451e-02 0.93
17 227 441 7.304e-04 2.030e-04 0.28
18 295 603 5.744e-04 1.723e-04 0.30
Table 1: Selected convergence and error estimator results of loop , having slabs and the max. number of mesh cells on a slab , for Sec. 3.

By Fig. 9 the distribution of the time subinterval lengths of over are illustrated. The goal (10) is reached for the 18th loop of the optimisation problem (Tab. 1) with an absolute error .

The numerical results (Fig. 8-9) highlight the local space-time refinement characteristic of the solver in a self-adaptive way, while the finer space-time cells are located inside and close to the boundary of the control volume , . Remark that not any refinement is done for of this parabolic problem since the control volume is not active.

4 Impact

The software gives new and efficient insights to the implementation of goal-oriented mesh adaptivity which supports the development of cost-effective numerical screening tools based on finite element approaches for fundamental problems in science and engineering. Derivatives of the software have already been used for stabilised stationary and instationary convection-dominated problems with goal-oriented mesh adaptivity in Bruchhaeuser2018 ; Bruchhaeuser2018a as well as for coupled deformation and diffusion-driven transport in Bause2017 ; Koecher2015 . The intended user group of researchers and engineers using and developing finite element based numerical simulation screening tools for several problems is widespread and globally distributed; an example would be the extension of the work Ahmed2015 in which they are using higher order variational time discretisations for convection-dominated transport with non-goal oriented adaptive time step control, or, an efficient instationary extension to the approach presented in Endtmayer2017 . The software is currently not used in commercial settings, but the license allows for the integration into commercial tools.











Figure 9: Distribution of of the space-time slabs , , for Sec. 3.

5 Conclusions

This original software publication provides efficient and scalable data structures and algorithms for the implementation of goal-oriented mesh adaptivity in a highly modular way. The key ideas of the data structures and algorithms can be reused, even independently of the programming language, in any adaptive finite element code. The performance and applicability of the software is shown with an illustrative example and several others are shipped with the code. The work yields a major breakthrough as a freely available open-source implementation with extensive in-source documentation for the dual weighted residual method in an application-free way such that it can be easily adopted by other frameworks. Ongoing work is on the distributed-memory parallelisation, on the serialisation of the data storage hiding memory operations, and on finding optimal tuning parameters of the minimisation problem.

Required Metadata

Current code version

Ancillary data table required for subversion of the codebase. Kindly replace examples in right column with the correct information about your current code, and leave the left column as it is.

Nr. Code metadata description Please fill in this column
C1 Current code version Version 1.0.0
C2 Permanent link to code/repository used for this code version dtm-project/dwr-diffusion
C3 Legal Code License dtm-project/dwr-diffusion/License
C4 Code versioning system used git
C5 Software code languages, tools, and services used C++.17; cmake, gcc, mpi, optional: paraview, doxygen
C6 Compilation requirements, operating environments & dependencies deal.II with hdf5; Linux (Fedora, CentOS 7, RHEL 7, etc.), MacOS, Windows WSL
C7 If available Link to developer documentation/manual dtm-project/dwr-diffusion
C8 Support email for questions
Table 2: Code metadata (mandatory)