1 Motivation and significance
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 functionand 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).
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
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.
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.
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 functionfrom (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.
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 https://ithub.com/dealii/candi
dandi && ./candi.sh
and follow the instructions. Then download, compile and run
it clone https://ithub.org/dtm-project/dwr-diffusion
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.
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 .
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.
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.
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.
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.
- (1) R. Becker and R. Rannacher, Weighted a posteriori error control in FE methods, In Proc. of the 2nd EC on Numer. Math. and Adv. Appl. ENUMATH 1997, H.G. Bock et.al. (eds), World Scientific, Singapore, 1998.
- (2) R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer. 10:1-102, doi:10.1017/S0962492901000010, 2001.
- (3) B. Endtmayer and T. Wick, A partition-of-unity dual-weighted residual approach for multi-objective error estimation applied to elliptic problems, Comput. Meth. Appl. Math. 17(4):1-25, doi:10.1515/cmam-2017-0001, 2017.
- (4) M.P. Bruchhäuser, K. Schwegler and M. Bause, Numerical study of goal-oriented error control for stabilized finite element methods, In Adv. Finite Element Meth. with Appl., T. Apel et.al. (eds), Lecture Notes in Comput. Sci. and Engrg., Springer, Berlin, p.1-19, accepted, arXiv:1803.10643, 2018.
- (5) M.P. Bruchhäuser, K. Schwegler and M. Bause, Dual weighted residual based error control for nonstationary convection-dominated equations: potential or ballast?, submitted, p.1-13, arXiv:1812.06810, 2018.
- (6) W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhäuser, Basel, 2003.
- (7) G. Alzetta, D. Arndt, W. Bangerth, V. Boddu, B. Brands, D. Davydov, R. Gassmoeller, T. Heister, L. Heltai, K. Kormann, M. Kronbichler, M. Maier, J.-P. Pelteret, B. Turcksin and D. Wells, The deal.II library, Version 9.0, J. Numer. Math. 26(4):173-183, doi:10.1515/jnma-2018-0054, 2018.
- (8) U. Köcher, Variational space-time methods for the elastic wave equation and the diffusion equation, Ph.D. thesis, Mech. Engrg. Helmut-Schmidt-University Hamburg, p. 1-188, urn:nbn:de:gbv:705-opus-31129, 2015.
- (9) M. Bause, F.A. Radu and U. Köcher, Space-time finite element approximation of the Biot poroelasticity system with iterative coupling, Comput. Meth. Appl. Mech. Engrg. 320:745-768, doi:10.1016/j.cma.2017.03.017, 2017.
- (10) N. Ahmed and V. John, Adaptive time step control for higher order variational time discretizations applied to convection-diffusion equations, Comput. Meth. Appl. Mech. Engrg. 285:83-101, doi:10.1016/j.cma.2014.10.054, 2015.
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||https://github.com/ dtm-project/dwr-diffusion|
|C3||Legal Code License||https://github.com/ 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||https://github.com/ dtm-project/dwr-diffusion|
|C8||Support email for email@example.com|