Variational time discretization for mixed finite element approximations of nonstationary diffusion problems

Journal of Computational and Applied Mathematics 289 (2015) 208–224

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Variational time discretization for mixed finite element approximations of nonstationary diffusion problems Markus Bause ∗ , Uwe Köcher Helmut Schmidt University, University of the Federal Armed Forces Hamburg, Department of Mechanical Engineering, Holstenhofweg 85, 22043 Hamburg, Germany

article

info

Article history: Received 12 September 2014 Received in revised form 16 January 2015 Keywords: Variational time discretization Continuous and discontinuous Galerkin methods Mixed finite element method Space–time Galerkin methods

abstract We develop and study numerically two families of variational time discretization schemes for mixed finite element approximations applied to nonstationary diffusion problems. Continuous and discontinuous approximations of the time variable are encountered. The solution of the arising algebraic block system of equations by a Schur complement technique is described and an efficient preconditioner for the iterative solution process is constructed. The expected higher order rates of convergence are demonstrated in numerical experiments. Moreover, superconvergence properties are verified. Further, the efficiency and stability of the approaches are illustrated for a more sophisticated three-dimensional application of practical interest with discontinuous and anisotropic material properties. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Numerical simulations of time dependent single and multiphase phase flow and multicomponent transport processes in porous media are desirable in several fields of natural sciences and in a large number of branches of technology; cf. e.g. [1]. The categories are far from clear-cut, but include environmental engineering (groundwater and soil contamination), hydrology, geophysics, agricultural, sanitary and chemical engineering, metal casting, ceramic engineering as well as petroleum engineering, geothermal energy recovery, carbon (dioxide) capture and storage in geological structures, gas injection as an efficient method of enhanced oil recovery and the technologies of wood, paper, cement, lightweight composite materials (aerospace engineering) and pharmaceutics. Further, temperature changes may accompany the physical and chemical processes. Temperature variation then influences the multiphase flow and multicomponent transport by altering physical and chemical properties of the phases such as densities, viscosities and reaction rates. The accurate numerical approximation of such flow and transport processes in heterogeneous and anisotropic porous media continues to be a challenging task. As a prototype model system we may consider the nonlinear set of coupled partial differential equations

∂t (Θ cl ) + ∇ · (Dl ∇ cl − Q cl ) = Rl (c1 , . . . , cL ) in Ω × I , on ΓD × I ,

(1.2)

−(Dl ∇ cl − Q cl ) · n = gF ,i on ΓF × I ,

(1.3)

cl (·, 0) =

(1.4)

cl = gD,l

∗

cl0

in Ω × {0}

Corresponding author. E-mail address: [email protected] (M. Bause). URL: http://www.hsu-hh.de/mb-mathe (M. Bause).

http://dx.doi.org/10.1016/j.cam.2015.02.015 0377-0427/© 2015 Elsevier B.V. All rights reserved.

(1.1)

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

209

for l = 1, . . . , L, with L ∈ N. Eq. (1.1) includes the effects of advective, diffusive and reactive transport in a multicomponent system with the unknown vector c = (c1 , . . . , cL )⊤ of concentrations of mobile chemical species. In (1.1)–(1.4), we denote by I = (0, T ) a finite time interval and Ω ⊂ Rd , with d = 2 or d = 3, is a domain with Lipschitz-continuous boundary Γ which is split into a flux boundary ΓF and a Dirichlet boundary ΓD . The outer unit normal vector to Ω is denoted by n. The rates Ri (c1 , . . . , cL ) model chemical reaction rates and sources and sinks. Further, Θ denotes the volumetric water content, Q the Darcy flux and Dl ∈ Rd,d , with d = 2 or d = 3, is the diffusion–dispersion tensor of the l-th chemical species. If reactive transport in porous media is considered, the water flux Q and the volumetric water content Θ are typically computed by solving Richards’ equation, equipped with appropriate initial and boundary conditions, which (in its mixed pressure formulation) reads as

∂t (Θ (ψ)) + ∇ · Q = 0,

Q = −K (Θ (ψ))∇(ψ + z ).

(1.5)

In (1.5), ψ is the hydraulic pressure head, Θ (·) is the given parametrization of the water content, K (·) is the given hydraulic conductivity and z denotes the height against the gravity direction. The applicability and value of the mixed finite element method (MFEM) and its hybrid variant (MHFEM) for solving the systems (1.1)–(1.5), respectively, have been demonstrated in a wide class of works; cf. [2–10] and the references therein. Appreciable advantages of mixed finite element methods are their local mass conservation property (cf. [11,12]) and the fact that they provide an explicit flux approximation as part of the formulation itself which is of particular importance for coupled flow and transport processes as, for instance, given by the set of Eqs. (1.1)–(1.5), where the Darcy flow field Q of the transport problems (1.1) is calculated from Richards’ equation (1.5). While the discretization in space involves a significant set of complexities and challenges, temporal approximations for (highly) transient processes in porous media have received relatively little interest (cf., e.g., [3,7,8,13]) and have been limited to traditional non-adaptive first and second order methods. Rigorous studies of higher order time discretizations are still missing in this field of application even though multiscale phenomena in the time evolution and strongly time dependent processes are often present, for instance, if complex chemical reaction mechanisms arise in Eq. (1.1), such that from the point of view of accuracy of numerical predictions the application of higher order time discretization schemes is recommendable and useful. The Galerkin method is a known approach to solve time dependent problems; cf. e.g., [14–18]. So far, this variational approach has rarely been used in practice despite of its significant advantages like a uniform variational approach for stability and error analyses, the natural construction of higher order methods, the applicability of goal-oriented error control [19] based on the dual weighted residual approach and the applicability of adaptive finite element techniques for changing the polynomial degree as well as the length of the time intervals. One reason for this might be the higher algorithmic complexity of solving the resulting algebraic systems of equations. Recently, higher order Galerkin time stepping schemes of variational type have been developed and studied for the heat equation [17,20], the nonstationary Stokes and Navier–Stokes equations [15,21,22] and wave problems [16,23,24]. In simulations for benchmarking configurations the efficiency and accuracy of continuous Galerkin–Petrov and discontinuous Galerkin time discretizations were demonstrated [15,16,20–22, 24–26]. Here, we introduce two families of continuous and discontinuous variational time discretization schemes that are combined with mixed finite element approximations of the spatial variables. The class of continuous approximations of the time variable is characterized by a Galerkin–Petrov type approach with piecewise polynomials of order r and, therefore, is referred to as the cGP(r) method. These schemes belong to the classes of A-stable methods [17]. This approach has already been studied in [14] where error estimates have been proved and superconvergence results have been observed. Further error estimates have been provided in [17,26]. The discontinuous variational approximation of the time variable results in a discontinuous Galerkin method with piecewise polynomials of order r. We refer to this approach as the dG(r) method. The resulting numerical schemes are known to belong to the class of strongly A-stable methods [26,27]. Error estimates for the semi-discretization in time of evolution equations of parabolic type are provided in [18]. For the mixed approximation in space we use the Raviart–Thomas family of inf–sup stable pairs of finite elements on quadrilateral and hexahedral elements. The fully discrete approximation schemes are developed by means of a Rothe type approach by discretizing the partial differential equation first in time and then in space. The schemes are studied numerically for a prototype diffusion model only. Their application to more general transport problems, for instance to convection–diffusion–reaction models and systems of equations, is from the algorithmic point of view straightforward. In a forthcoming work [28], an error analysis is presented for the classes of discretization schemes that are considered here. This numerical analysis is currently limited to problems with self-adjoint differential operators. An extension to more general problems is still an ongoing work. In this work, the solution process of the arising linear systems of equations is described carefully, since this step has shown to be an essential issue in the application of our variational time discretization schemes. This is done for the cGP(2) and dG(1) methods, that are used in this work for the numerical convergence studies and further numerical experiments. The extension of the solution strategy and the construction of efficient preconditioners for still higher order temporal approximations is still an ongoing work. In contrast to former works [3,4,29,30], the mixed finite element approach is applied here in a non-hybrid framework and the linear systems of equations are solved by Schur complement techniques along with nested preconditioned conjugate gradient iterations. This is done due to the experience made by one of the authors with applying hybrid mixed finite element methods to nonlinear problems and coupled systems of equations and the resulting high algorithmic complexity. In this work, to solve the linear systems of equations, the degrees of freedom of the vectorvalued variable are eliminated first. In the resulting reduced system of the equations for the cGP(2) and dG(1) approaches,

210

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

respectively, the coefficients with respect to the representation in the spatial finite element basis of one coefficient function of the temporal expansion in terms of basis functions are then further eliminated by computing the Schur complement technique. The resulting systems of equations are then solved. As we shall show below, this solution strategy requires the application of highly adapted and efficient preconditioners that are derived from preconditioning techniques for fourth order problems. This work is organized as follows. In Section 2 we introduce our family of continuous approximations of the time variable with piecewise polynomials of order r. This semi-discretization in time is then combined with a mixed finite element approximation of the spatial variables. Finally, the iterative solution of the resulting linear system of equations is described. In Section 3 our family of discontinuous variational time discretization schemes is briefly described. In Section 4 numerical experiments illustrating the experimental order of convergence of the considered schemes are presented. Moreover, the space–time Galerkin approach is applied to a sophisticated three-dimensional diffusion problem with anisotropic and heterogeneous material properties in order to demonstrate its stability, robustness and potential for solving problems of practical interest. 2. Continuous variational time discretization and mixed finite element approximation In this section we develop a family of continuous variational time discretization schemes that are combined with mixed finite element discretizations [12] of the spatial variables. In order to reduce the technical complexity we consider the prototype diffusion model

∂t u − ∇ · (D∇ u) = f in Ω × I , (2.1) on ∂ Ω × I , (2.2) u(·, 0) = u0 in Ω , (2.3) where I = (0, T ) with final time T > 0 and Ω ⊂ Rd , d = 2 or d = 3, is a polygonal or polyhedral bounded domain, respectively. We assume that the symmetric diffusion matrix D = D(x) = (dij (x))di,j=1 satisfies dij ∈ L∞ (Ω ) and dij = dji for i, j = 1, . . . , d and is (uniformly) elliptic, i.e. ξ ⊤ D(x)ξ ≥ θ|ξ|2 for almost every x ∈ Ω and all ξ ∈ Rd and some constant θ > 0. Further we let f ∈ L2 (I , L2 (Ω )) and u0 ∈ H01 (Ω ) such that u=0

the existence of a unique weak solution u ∈ H 1 (I ; L2 (Ω )) ∩ L2 (I ; H01 (Ω )) ∩ C (I ; L2 (Ω ))

(2.4)

of problem (2.1)–(2.3) is ensured; cf. [31]. Here, the initial boundary value problem (2.1)–(2.3) is studied as a prototype model for the more sophisticated problems (1.1)–(1.4) or (1.5), respectively, as well as for multiphase systems [1,6]. More general boundary conditions for (1.1) can be incorporated in the model problem (2.1)–(2.3) and in its mixed finite element discretization in the usual way. Throughout this paper, standard notations and conventions will be used. We denote by H m (Ω ) the Sobolev space of L2 functions with derivatives up to order m in L2 (Ω ) and by ⟨·, ·⟩ the inner product in L2 (Ω ). Further, for the solution spaces of the mixed problem formulation we use the following notation, V = H (div; Ω ),

W = L2 (Ω ),

(2.5)

where H (div; Ω ) = {q ∈ L (Ω ) | ∇ · q ∈ L (Ω )}. Further, we consider the following Bochner spaces of functions with values in the Banach space X , 2

2

C (I ; X ) = {w : [0, T ] → X | w is continuous},

 



 ∥w(t )∥2X dt < ∞ , I   H 1 (I ; X ) = w ∈ L2 (I ; X ) | ∂t w ∈ L2 (I ; X ∗ ) ,

L2 (I ; X ) = w : (0, T ) → X 

that are equipped with their natural norms and where the time derivative ∂t is understood in the sense of distributions on (0, T ). In order to derive our families of discretization schemes, we first define the auxiliary flux variable q = −D(x)∇ u. Since ∂t u ∈ L2 (I ; L2 (Ω )) is satisfied by means of (2.4), it holds that q ∈ L2 (I ; V ). The pair {u, q} ∈ H 1 (I ; W ) ∩ C (I ; W ) × L2 (I ; V ) then satisfies the set of variational equations T



⟨∂t u, w ⟩ dt + 0

⟨D−1 q, v ⟩ dt − 0

⟨∇ · q, w⟩ dt = 0

T



T



T



⟨f , w⟩ dt ,

(2.6)

0 T



⟨u, ∇ · v ⟩ dt = 0

(2.7)

0

for all w ∈ L2 (I ; W ) and v ∈ L2 (I ; V ) and the initial condition u(0) = u0 . The global in time variational problem (2.6), (2.7) is now the starting point for our semidiscretization in time.

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

211

2.1. Variational discretization in time For the discretization in time we decompose the time interval (0, T ] into N subintervals In = (tn−1 , tn ], where n ∈ {1, . . . , N } and 0 = t0 < t1 < · · · < tn−1 < tn = T . Further τ denotes the discretization parameter in time and is defined as the maximum time step size τ = max1≤n≤N τn , where τn = tn − tn−1 . We introduce the subspaces of piecewise polynomials of order r in time,

Xr (X ) := {u ∈ C (¯I ; X )|u| I n ∈ Pr (I n ; X ), ∀n} ⊂ {u ∈ L2 (I ; X ) | dt u ∈ L2 (I ; X ∗ )}, Yr (X ) := {w ∈ L2 (I ; X )|w | In ∈ Pr (In ; X ), ∀n}, where



  r   j j j  Pr (J ; X ) = w : J → X w(t ) = ξn t , ξn ∈ X , j = 0, . . . , r . j =0

For the family of continuous variational time discretizations the spaces Xr (X ) of continuous functions act as solution spaces whereas the spaces Yr −1 (X ) consisting of piecewise polynomials that are discontinuous at the end points of the time intervals are used as test spaces. Since in this case the solution and test spaces differ, a discretization of Galerkin–Petrov type is thus obtained. A semi-discrete variational approximation in time of the mixed form of problem (2.1)–(2.3) is now defined by solving the variational equations (2.6), (2.7) in discrete subspaces: Find {uτ , qτ } ∈ Xr (W ) × Xr (V ) such that T



⟨∂t uτ , wτ ⟩ dt + 0

⟨∇ · qτ , wτ ⟩ dt = 0

T



T



⟨D−1 qτ , vτ ⟩ dt −

⟨f , wτ ⟩ dt ,

(2.8)

0



0

T



T

⟨uτ , ∇ · vτ ⟩ dt = 0

(2.9)

0

for all wτ ∈ Yr −1 (W ) and v ∈ Yr −1 (V ) with the initial condition uτ (0) = u0 . We refer to the solution of Eqs. (2.8), (2.9) as the continuous Galerkin–Petrov method with piecewise polynomials of order r and use the notation cGP(r). To ensure the existence of solutions to (2.8), (2.9), it is sufficient to use the test spaces Yr −1 (W ) and Yr −1 (V ) with piecewise polynomials of order r − 1, since the continuity constraints at the discrete time points tn , n = 0, . . . , N, that are implied by the definition of the solution spaces Xr (W ) and Xr (V ) yield an additional constraint. Along with an error analysis the existence and uniqueness of solutions to Eqs. (2.8), (2.9) is shown in a forthcoming work [28]. By using discontinuous test basis functions wτ (t ) = wψn,i (t ) and vτ = v ψn,i (t ) for i = 1, . . . , r with arbitrary time independent functions w ∈ W and v ∈ V , respectively, and piecewise polynomial functions ψn,i : I → R that are of order r − 1 on In and vanish outside I \ I n , we can recast the variational equations (2.8), (2.9) as a time marching scheme: For n = 1, . . . , N find uτ |In ∈ Pr (In ; W ) and qτ |In ∈ Pr (In ; V ) such that



⟨∂t uτ , w⟩ ψn,i (t ) dt +



In



⟨∇ · qτ , w⟩ ψn,i (t ) dt = In

⟨D−1 qτ , v ⟩ ψn,i (t ) dt − In



T

⟨f , w⟩ ψn,i (t ) dt ,

(2.10)

0



⟨uτ , ∇ · v ⟩ ψn,i (t ) dt = 0

(2.11)

In

for all w ∈ W and v ∈ V and i = 1, . . . , r with the continuity constraints uτ |In (tn−1 ) = uτ |In−1 (tn−1 ) for n ≥ 2 and the initial condition uτ |In (tn−1 ) = u0 for n = 1. To determine uτ |In and qτ |In , respectively, we represent these functions in terms of basis functions with respect to the time variable of Xr (W ) and Xr (V ), respectively, such that uτ |In (t ) =

r 

Unj ϕn,j (t ) and qτ |In (t ) =

j =0

r 

Qnj ϕn,j (t )

(2.12)

j=0 j

j

with coefficient functions Un ∈ W and Qn ∈ V for j = 0, . . . , r and polynomial basis functions ϕn,j ∈ Pr (In ; R) that are Lagrange functions with respect to r + 1 nodal points tn,j ∈ In satisfying the conditions

ϕn,j (tn,i ) = δi,j ,

for i, j = 0, . . . , r ,

(2.13)

with the Kronecker symbol δi,j . For the treatment of the continuity constraints and the initial condition we choose tn,0 = tn−1 . This implies that Un0 = uτ |In−1 (tn−1 ) if n ≥ 2

and Un0 = u0

if n = 1.

(2.14)

212

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

The other points tn,1 , . . . , tn,r are chosen as the quadrature points of the r-point Gaussian quadrature formula on In which is exact if the function to be integrated is a polynomial of degree less than or equal to 2r − 1. The basis functions ϕn,j ∈ Pr (In ; R)

of (2.12), for j = 0, . . . , r, are defined via the affine reference transformation onto ˆI = [0, 1], tˆ :=

1

τn

(t − tn−1 ) ∈ ˆI ,

for t ∈ In ,

by means of

ϕn,j (t ) = ϕˆ j (tˆ) with ϕˆ j ∈ Pr (ˆI ; R) and ϕˆ j (tˆi ) = δi,j ,

for i, j = 0, . . . , r ,

(2.15)

where tˆ0 = 0 and tˆi , i = 1, . . . , r, are the standard r point Gaussian quadrature points on the interval ˆI = [0, 1]. For further technical details of this transformation we refer to [16,17,20,21]. Similarly, the test basis functions ψn,i ∈ Pr −1 (In ; R) are

ˆ i ∈ Pr −1 (ˆI ; R) through defined by suitable reference basis functions ψ ψn,i (t ) = ψˆ i (tˆ) with ψˆ i (tˆl ) = δi,l for i, l ∈ {1, . . . , r },

(2.16)

where tˆl , l = 1, . . . , r, again denote the quadrature points of the r point Gaussian formula on ˆI. Now we transform all the time integrals in (2.10), (2.11) to the reference interval ˆI and obtain the following system of variational problems for j j the coefficient functions Un ∈ W and Qn ∈ V of the representation (2.12): For n = 1, . . . , N find coefficient functions j j {Un , Qn } ∈ W × V , for j = 1, . . . , r, such that r 

αˆ ij ⟨Unj , w⟩ + τn βˆ i,i ⟨∇ · Qni , w⟩ = τn βˆ i,i ⟨f (tn,i ), w⟩

(2.17)

j =0

⟨D−1 Qni , v ⟩ − ⟨Uni , ∇ · v ⟩ = 0

(2.18)

for i = 1, . . . , r and all {w, v } ∈ W × V , where ∈ W is defined by means of the continuity constraints (2.14). For the derivation of (2.17), (2.18) from (2.10), (2.11) we note the following issues. Firstly, the time integration in (2.10) of the product of the source term f with the test functions is done in practice by numerical quadrature. Therefore, we tacitly replaced the integrand f in the integral on the right-hand side of Eq. (2.10) by its Lagrange interpolate Πr f ∈ Pr (In ; L2 (Ω )) that is given by Un0

Πr f (t ) =

r 

f (tn,j )ϕn,j (t ),

for t ∈ In .

j=0

This yields the term on the right-hand side of Eq. (2.17). Further, for deriving (2.17), (2.18) we also used that



d In

dt

ϕn,j (t ) · ψn,i (t ) dt =



d

ˆI

dtˆ d

= ωˆ i · 

ϕn,j (t ) · ψn,i (t ) dt = τn In

 ˆI

ϕˆ n,j (tˆ) · ψˆ n,i (tˆ) dtˆ =

dtˆ

r  m=1

ϕˆ j (tˆi ) =: αˆ i,j ,

d dtˆ

ϕˆ j (tˆm ) · ψˆ i (tˆm )

for i = 1, . . . , r , j = 0, . . . , r ,

ϕˆ n,j (tˆ) · ψˆ n,i (tˆ) dtˆ = τn ·

= τn · ωˆ i · δi,j := τn · βˆ i,j ,

ωˆ m ·

r 

(2.19)

ωˆ m · ϕˆ j (tˆm ) · ψˆ i (tˆm )

m=1

for i = 1, . . . , r , j = 0, . . . , r ,

(2.20)

by means of definitions (2.13), (2.15) and (2.16) of the trial and test basis functions, respectively. In a forthcoming work [28] the scheme (2.17), (2.18) with numerical integration of the right-hand side terms is also studied with respect to the existence and uniqueness of solutions and its approximation properties. Remark 2.1. The advantage of using the Gaussian quadrature rule for evaluating the time integrals in (2.10), (2.11) is that in the resulting scheme (2.17), (2.18) the coefficient function Qn0 does not arise such that no flux approximation at the initial time point tn−1 of the subinterval In is involved in the scheme. Therefore, the initial flux at time t = 0, that we do j j not have in practice, is not involved and needed for calculating the coefficient functions Un and Qn for j = 1, . . . , r at the intermediate Gaussian quadrature points tn,l ∈ (tn−1 , tn ) for l = 1, . . . , r on all subintervals In . For evaluating the continuity constraint (2.14) for n ≥ 2, the value uτ |In−1 (tn−1 ) is further needed. The quantity uτ |In−1 (tn−1 ) has to be computed since tn−1 does not belong to the set of the Gaussian quadrature points tn−1,l ∈ (tn−2 , tn−1 ), for l = 1, . . . , r, on the time interval In−1 . In our implementation we calculate uτ |In−1 (tn−1 ) explicitly by evaluating the fully discrete counterpart of the first of the representations in (2.12) on the subinterval In−1 . For the evaluation of the first of the representations in (2.12) on the subinterval I1 = (t0 , t1 ] we recall that the initial value for the scalar variable u ∈ C (I ; L2 (Ω )) is well-defined such that the second of the equations in (2.14) is reasonable. For sufficiently regular solutions of the initial–boundary value problem (2.1)– (2.3) such that an initial flux q(0) exists, an approximation of the flux variable at the endpoints tn of the time intervals In can

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

213

be obtained in an analogous way by using the second of the representations in (2.12) and continuity constraints similarly to (2.14). Alternatively, in [21] an order preserving post processing technique based on interpolation is suggested and used for evaluating the velocity and pressure of Stokes approximations at the endpoints of the time intervals. This might be applied here as well and also used in the case that the solution is not sufficiently regular such that an initial flux is not available. 2.2. The fully discrete cGP(r)-MFEM approach Now, we present the fully discrete approximation scheme that is obtained by discretizing the variational equations (2.17), (2.18) with respect to their spatial variables. For this we choose a pair of finite element spaces Wh ⊂ W and Vh ⊂ V = H (div; Ω ) satisfying the inf–sup stability condition; cf. [11,12]. Here, we denote by Th = {K } a finite element decomposition of mesh size h of the polyhedral domain Ω into closed subsets K , quadrilaterals in two dimensions and hexahedrals in three dimensions. Since the software library deal.ii [32] that we use for our implementation of the schemes and the calculations (cf. Section 4) allows only quadrilateral and hexahedral elements, we restrict ourselves to these types of elements in the following. Triangular and tetrahedral elements can be treated in an analogous way. The decompositions are assumed to be face to face. We use the notation ‘‘face’’ in the two- and three-dimensional case. If d = 2, then the faces are the edges of the quadrilaterals. For brevity, we restrict ourselves to the class of Raviart–Thomas elements for the twodimensional case and to the class of Raviart–Thomas–Nédélec elements in three space dimensions. An extension to other inf–sup stable pairs of mixed finite element spaces is straightforward. The construction of the discrete function spaces on quadrilateral and hexahedral finite elements is done by a transformation TK : Kˆ → K from the reference element Kˆ = [0, 1]d , with d = 2 or d = 3 respectively, through a diffeomorphism TK for all K ∈ Th . For this, let



  p1  p2   j pi,j xi1 x2 , pi,j ∈ R , Qˆ p1 ,p2 := pˆ : [0, 1]2 → Rpˆ (ˆx) = i=0 j=0    p3 p1  p2    3 i j k p1 ,p2 ,p3  ˆ pi,j,k x1 x2 x3 , pi,j,k ∈ R . Q := pˆ : [0, 1] → Rpˆ (ˆx) = i=0 j=0 k=0

p

p

We then define the discrete subspaces Wh ⊂ W and Vh ⊂ V by

    −1 p,p   ˆ  w ∈ W w ◦ T ∈ Q , for K ∈ T for d = 2, h  K  K p   Wh = Wh :=      w ∈ W wK ◦ TK−1 ∈ Qˆ p,p,p , for K ∈ Th for d = 3,

(2.21)

    −1 p+1,p p,p+1   ˆ ˆ  ∈ Q × Q , for K ∈ T for d = 2, v ∈ V v ◦ T h K  K  p   Vh = Vh :=      v ∈ V vK ◦ TK−1 ∈ Qˆ p+1,p,p × Qˆ p,p+1,p × Qˆ p,p,p+1 , for K ∈ Th

(2.22)

and

for d = 3

with space dimensions dim Wh = NWh := |Th | · (p + 1)d

 and

dim Vh = NVh :=

|Th | −

|FhI | 4



· [d(p + 1)d−1 (p + 2)],

(2.23)

where |Th | denotes the cardinality of Th and |FhI | is the number of interior faces. Then, we obtain our fully discrete cGP(r)MFEM approximation scheme, by solving the variational problem (2.17), (2.18) in the pair of discrete subspaces Wh ⊂ W j j and Vh ⊂ V : For n = 1, . . . , N find coefficient functions {Un,h , Qn,h } ∈ Wh × Vh , for j = 1, . . . , r, such that r 

αˆ ij ⟨Unj ,h , wh ⟩ + τn βˆ i,i ⟨∇ · Qni,h , wh ⟩ = τn βˆ i,i ⟨f (tn,i ), w⟩,

(2.24)

j=0

⟨D−1 Qni,h , vh ⟩ − ⟨Uni ,h , ∇ · vh ⟩ = 0

(2.25)

for i = 1, . . . , r and all {wh , vh } ∈ Wh × Vh , where constraints (2.14), i.e. Un0,h =

r 

j

Un−1,h ϕn−1,j (tn−1 ) if n ≥ 2

Un0,h

∈ Wh is defined by means of the discrete counterpart of the continuity

and Un0,h = Ph u0

j =0

In (2.26), Ph : L2 (Ω ) → Wh denotes the L2 projection onto Wh .

if n = 1.

(2.26)

214

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

In terms of globally defined basis functions {wh,l }l=1,...,NW of Wh and {vh,m }m=1,...,NV of Vh the fully discrete coefficient h

j

h

j

functions Un,h ∈ Wh , for j = 0, . . . , r, and Qn,h ∈ Vh , for j = 1, . . . , r, admit the representation j Un,h

(x) =

NW h 

j un , l

wh,l (x) and

j Qn,h

(x) =

l =1

NV h 

qjn,m vh,m (x).

(2.27)

m=1

Next, in (2.24), (2.25) we choose the basis functions of Wh and Vh , respectively, as test functions. Further, we define the N ,N N ,N N ,N matrices M ∈ R Wh Wh , B ∈ R Vh Wh and A ∈ R Vh Vh by ms,l = ⟨wh,l , wh,s ⟩, l, s = 1, . . . , NWh ,

M = (ms,l )s,l , B = (bt ,l )t ,l ,

bt ,l = ⟨wh,l , ∇ · vh,t ⟩, t = 1, . . . , NVh , l = 1, . . . , NWh ,

A = (at ,m )t ,m ,

at ,m = ⟨D−1 vh,m , vh,t ⟩, m, t = 1, . . . , NVh , j

as well as the vectors Fn ∈ R Fnj

NW

h

by

(j)

= (Fs )s = ⟨f (tn,j ), wh,s ⟩,

s = 1, . . . , NWh ,

for j = 0, . . . , r ,. These matrices and vectors are assembled, as usual, elementwise. For the implementational details we refer to [24]. For the coefficient vectors of the expansions in (2.27) we use the notation j ujn = (un,1 , . . . , un,NW )⊤ h

j and qjn = (qn,1 , . . . , qn,NV )⊤ .

(2.28)

h

Then, we recast the fully discrete variational problem (2.24), (2.25) in its algebraic form: For n = 1, . . . , N find coefficient j j N N vectors {un , qn } ∈ R Wh × R Vh , for j = 1, . . . , r, of the expansion (2.27), (2.28) such that r 

αˆ ij Mujn + τn βˆ i,i B⊤ qin = τn βˆ i,i Fni ,

(2.29)

j =0

Aqin − Buin = 0

(2.30)

for i = 1, . . . , r, where u0n ∈ R h is defined by means of the discrete continuity constraints (2.26) along with the representation (2.27) of Un0,h in terms of finite element basis functions and the definition (2.28), i.e. NW

r  j =0

j

Un−1,h ϕn−1,j (tn−1 ) = Un0,h =

NW h 

u0n,l wh,l (x) if n ≥ 2

l=1

and Ph u0 = Un0,h =

NW h 

u0n,l wh,l (x)

if n = 1.

(2.31)

l =1

Eqs. (2.29), (2.30) represent the algebraic form of our fully discrete scheme combining a continuous variational time discretization with a mixed finite element approximation in space. For solving Eqs. (2.29), (2.30) we do not apply an additional hybridization technique as it is done, for instance, in [3,4,29]. In fact, we solve the algebraic system directly by a Schur complement technique. This is shown now for the case of an approximation with piecewise polynomials of second order in time, corresponding to the parameter choice r = 2, for that convergence of third order in the L2 -norm as the natural norm of the discrete scheme can be expected (cf. Section 4). The solution process requires a highly adapted preconditioner for the iterative solver. The case r = 1 is not in the scope of our interest in this work, since it is very close to the wellknown Crank–Nicolson scheme (cf. [20,23] for details), and we focus here on higher order temporal approximations. So far, we have not studied yet solution strategies for the algebraic systems of equations that result from discretizations in time with still higher order polynomial functions such that r ≥ 3. This remains a work for the future. Nevertheless, we made the experience that choosing r = 2 is sufficient from the point of view of accuracy and a reasonable compromise between accuracy, computational costs and algorithmic complexity for applications of practical interest in the context of mixed finite element methods. In particular, in Section 4 we will demonstrate numerically that superconvergence of fourth order is obtained in the end points tn of the subintervals In by a piecewise quadratic approximation in time with r = 2. Nevertheless, applications with high temporal fluctuations and multiple time scales, which for instance arise if fast chemical reactions or high temperature variations are present, might occur. In these cases, higher order time discretizations with r ≥ 3 are reasonable, recommendable and might even be necessary. 2.3. Solution of the algebraic system for the cGP(2)-MFEM approach We shall now describe the solution of the system (2.29), (2.30) for r = 2 corresponding to a piecewise quadratic approximation in time. Putting r = 2 we recast problem (2.29), (2.30) as the following block system of equations: For

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224 j

j

× RNVh , for j = 1, 2, of the expansion (2.27), (2.28) such that   1   qn αˆ 1,2 M τn βˆ 1,1 Fn1 − αˆ 1,0 Mu0n  q2    αˆ 2,2 M  τn βˆ 2,2 Fn2 − αˆ 2,0 Mu0n     n1  =     un  0 0 −B 0 u2

n = 1, . . . , N find coefficient vectors {un , qn } ∈ R

   

τn βˆ 1,1 B⊤

αˆ 1,1 M αˆ 2,1 M −B

0

τn βˆ 2,2 B⊤

0 A 0

0 A

215

0

NW

h

(2.32)

n

N

with u0n ∈ R Wh being defined by means of the discrete continuity constraints (2.31) along with definition (2.28). The coefficients αˆ i,j and βi,j in (2.32) are given by Eqs. (2.19) and (2.20) as

√ 3

αˆ 1,0 = − √ ,

αˆ 1,1 =

3

3 2

,

αˆ 1,2 =

3 2

√

3−1

√

3+3

,

3

αˆ 2,0 = √ ,

3

αˆ 2,1 =

3

2

3−1

√

3−3

,

αˆ 2,2 =

3 2

and

βˆ 1,0 = 0,

βˆ 1,1 =

1 2

,

βˆ 1,2 = 0,

βˆ 2,0 = 0,

βˆ 2,1 = 0,

βˆ 2,2 =

1 2

.

The 4 × 4 block matrix in (2.32) is of a complex structure with respect to the iterative solution of the system and the construction of an efficient preconditioner. Moreover, the corresponding condition number depends on the mesh size such that the iterative solution of these linear systems of equations by preconditioned Krylov-space methods is supposed to become expensive for finer meshes. For a more refined study of the characteristics of the block matrix in (2.32) and for a comparative analysis of different solution strategies and iterative solvers for (2.32) we also refer to [24]. In [20], a geometrical multigrid solver is presented for variational time discretization schemes applied to the heat equation that are combined with continuous Galerkin approximations of the spatial variables. For these schemes it is shown that the multigrid solver is much more efficient than a preconditioned BiCGStab solver, since it shows a robust convergence behaviour that is nearly independent of the space mesh size and the time step such that large time steps get feasible also with respect to efficient solution methods. In this work we use a different approach. In two steps we eliminate internal degrees of freedom and, thereby, reduce the system size. For the resulting symmetric system of equations in the essential unknown u2n we construct an efficient preconditioner. This enables us to use a preconditioned conjugate gradient method for solving the condensed linear system derived from (2.32). For the solution of (2.32) we firstly eliminate the degrees of freedom q1n and q2n of the flux approximation. For this we multiply the third and fourth of the equations in (2.32) by τ2n B⊤ A−1 and subtract the resulting equations from the first and second one of the equations in (2.32). This yields the reduced system of equations



αˆ 1,1 M +

 

τn 2

B⊤ A−1 B

αˆ 2,1 M

  τ n 1  1 αˆ 1,2 M Fn − αˆ 1,0 Mu0n  un  2 , τn ⊤ −1  u2 =  τn 2 0 n αˆ 2,2 M + B A B Fn − αˆ 2,0 Mun 2

(2.33)

2

whereas the flux unknowns can be calculated in a postprocessing step from Aq1n = Bu1n

Aq2n = Bu2n .

and

(2.34)

The block system (2.33) is now solved by calculating the Schur complement operator of the lower left block in the 2 × 2 block matrix. From the second block of Eqs. (2.33) we get that u1n = −

1

αˆ 2,1

M −1



αˆ 2,2 M +

τn 2



B⊤ A−1 B u2n −

τ

n

2

Fn2 − αˆ 1,0 Mu0n



.

Substituting this identity in the first block of Eqs. (2.33) then yields the linear systems of equations Su2n = (−α2,1 )

τ

n

2





Fn1 − αˆ 1,0 Mu0n + αˆ 1,1 M +

τn 2



B⊤ A−1 B M −1

τ

n

2



Fn2 − αˆ 1,0 Mu0n ,

(2.35)

where the Schur complement operator of the lower left block in the 2 × 2 block matrix is given by



S = αˆ 1,2 (−αˆ 2,1 )M + αˆ 1,1 M +

τn 2





B⊤ A−1 B M −1 αˆ 2,2 M +

τn 2



B⊤ A−1 B .

(2.36)

The coefficient vectors u0n in (2.33) and (2.35), respectively, are calculated by evaluating the continuity constraints in (2.31) for r = 2. Summarizing, the solution of the algebraic form (2.32) of the cGP(2)-MFEM approach is determined by solving

216

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

the following subproblems for n = 1, . . . , N:









⃗ u⃗0n + αˆ 1,1 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ M ⃗ −1 ⃗2n = (−αˆ 2,1 ) τ2n F⃗n1 − αˆ 1,0 M S⃗u 2 ⃗ −1 ⃗1n = − αˆ 1 M u 2,1



⃗ + αˆ 2,2 M

τ n ⃗ ⊤ ⃗ −1 ⃗ 2



⃗2n − B u

B A



τn 2

⃗ u⃗0n F⃗n2 − αˆ 2,0 M



τn 2





⃗ u⃗0n , F⃗n2 − αˆ 2,0 M

,

⃗ q⃗1n = B⃗u⃗1n , A ⃗ q⃗2n = B⃗u⃗2n , A

(2.37)









⃗ + αˆ 1,1 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ M ⃗ −1 αˆ 2,2 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ with S⃗ = αˆ 1,2 (−αˆ 2,1 )M 2 2 and

u0n

⃗ =

2 

j un − 1

⃗

ϕn−1,j (tn−1 ) if n ≥ 2

and

j =0

P h u0 =

NW h 

u0n,l wh,l (⃗ x)

if n = 1.

l=1

The first linear system of equations in (2.37) for the essential unknown vector u2n is solved by preconditioned conjugate gradient (CG) iterations. The Schur complement S is obviously symmetric. In the CG iterations the arising factors B⊤ A−1 B have not to be determined explicitly, since the application of the preconditioned CG method can be based on evaluating matrix–vector products and inner products between vectors as well as the application of the preconditioner. The technical details of our preconditioned CG iterations are given below and summarized in Table 1. The mass matrix M of the scalar variable is a block diagonal matrix with block sizes equal to the number of degrees of freedom per element such that its inverse of M −1 in the second system of Eqs. (2.37) and in the Schur complement S is available at low computational costs. For the efficient solution of the systems (2.37) the preconditioning technique for the first of Eqs. (2.37) is vitally important. The Schur complement S is ill-conditioned since it approximates in terms of the involved products of the matrix B a fourth order differential operator. The construction of our preconditioner is based on ideas developed in [33] where preconditioned solvers for a class of fourth order problems are in the scope of interest. The approach in [33] makes use of an inexact factorization of the Schur complement for which uniform bounds can be proved. In [34] the technique is adapted to continuous and discontinuous variational time discretizations that are combined with continuous Galerkin approximations in space and applied to parabolic problems. Here we apply the technique of [34] for preconditioning the operator S in Eq. (2.35). To solve Eq. (2.35), we use a conjugate gradient (CG for brevity) method with left preconditioning [35]. The preconditioned form of the system (2.35) is then given by Pµ−∗1 Su2n = Pµ−∗1 b1n , with a vector b1n being equal to the right-hand side of Eq. (2.35). We recall that left preconditioning amounts to applying the conjugate gradient algorithm with respect to the weighted inner product (cf. [35])

⟨x, y ⟩Pµ∗ := ⟨Pµ∗ x, y ⟩RNWh . Following the lines of [34], we define our left-sided preconditioner Pµ∗ by the matrix representation Pµ∗ := Kµ∗ M −1 Kµ∗ ,

with Kµ := µM +

τn 2

B⊤ A−1 B,

(2.38)

and some parameter µ > 0. Using the arguments given in [34], we can easily determine an optimal parameter µ∗ such that the condition number κ(S∗ ) = ∥S∗ ∥2 ∥S∗−1 ∥2 of the preconditioned operator S∗ = Pµ−∗1 S is minimized. Here we skip the technical calculations and refer to [24] for further details. By the approach we find that

µ∗ =

√

√

3

with κ(S∗ ) ≤ 8 − 4 3.

(2.39)

Our CG iterations are summarized in the pseudocode of Table 1. For this we also refer to [34]. Our implementation of the algorithm is discussed in [24]. The preconditioned CG iterations in Table 1 require the evaluation of the operator P −1 r = Kµ−∗1 M Kµ−∗1 r







(2.40)

which is done in the sequence of operations indicated by the brackets in Eq. (2.40). Therefore, the numerical effort of the algorithm lies in an efficient solution of a linear system with the matrix Kµ∗ = µM + τ2n B⊤ A−1 B. This is done by nested CG iterations using the preconditioned algorithm of Table 1 again. For these nested CG iterations we choose a left-sided preconditioner that is given in its matrix representation by

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

217

Table 1 Preconditioned conjugate gradient algorithm for Eq. (2.35).

 τn  ⊤  Kµ∗ = µ∗ M + B diag(A)−1 B 2

(2.41)

with the block-diagonal mass matrix M of the scalar variable and the approximate inverse diag(A)−1 of the mass matrix A of the vector-valued variable. In (2.41), diag(A) denotes the diagonal matrix that is built from the diagonal entries of A. The application of the preconditioner  Kµ∗ in the inner preconditioned CG loop is done by solving the corresponding linear systems of equations  Kµ∗ p = r. This nested approach is further used for the evaluation of the matrix–vector products B⊤ A−1 Bx and is described now in this context.  The quantities B⊤ A−1 (Bx) are computed from the inner to the outer brackets whenever they occur. Then the main effort comes through the implementation of the application of the inverse matrix A−1 . For this we can use iterative system solvers or parallel direct system solvers within our software DTM++ based on the deal.II library [32]. For the numerical calculations presented in Section 4 we used the parallel direct SuperLU, for that we got the best numerical performance properties. For more detailed studies regarding the application of iterative and parallel direct solvers and their performance properties we also refer to [24]. This type of a linear system solver is then also used to calculate the flux quantities q1n and q2n in the system (2.37) when the unknowns u1n and u2n are available. 3. Discontinuous variational time discretization and mixed finite element approximation 3.1. The fully discrete dG(r)-MFEM approach In this section we develop a family of discontinuous variational time discretization schemes that are combined with a mixed finite element discretization of the spatial variables. Basically, we follow the lines of Section 2. Therefore, we keep the derivation short and restrict ourselves to presenting only the differences to the construction of the family of continuous variational time discretization schemes in Section 2. The basic concept in the construction of discontinuous Galerkin approximations of time dependent problems can be found in [18]. For further technical details in the application of this time discretization methods we also refer to [16,20,26]. With the notation introduced in Section 2.1, the semi-discrete discontinuous Galerkin approximation of the variational problem (2.6), (2.7) with piecewise polynomials of order r, referred to as the dG(r) method, reads as follows: Find {uτ , qτ }

218

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

∈ Yr (W ) × Yr (V ) such that    N  N   ⟨∂t uτ , wτ ⟩ dt + ⟨∇ · qτ , wτ ⟩ dt + ⟨[uτ ]n−1 , wτ (tn+−1 )⟩ = ⟨f , wτ ⟩ dt , In

n =1

In

T



⟨D−1 qτ , vτ ⟩ dt −

n =1

(3.1)

In

T



⟨uτ , ∇ · vτ ⟩ dt = 0

(3.2)

0

0

for all wτ ∈ Yr (W ) and v ∈ Yr (V ) with the initial value uτ (0) := u0 . In (3.1), (3.2) we use the notation uτ (tn− ) =

lim uτ (t ),

uτ (tn+ ) =

t → tn − 0

lim uτ (t ),

[uτ ]n = uτ (tn+ ) − uτ (tn− ).

t → tn + 0

In the derivation of Eq. (3.1) integration by parts is applied twice, first in the time-continuous case, T



⟨∂t u, w⟩ dt = −

T



⟨u, ∂t w⟩ dt − ⟨u(0), w(0)⟩

(3.3)

0

0

for sufficiently smooth functions u, w with w(T ) = 0, and then in the time-discrete counterpart of (3.3),

−

N   n =1

⟨uτ , ∂t wτ ⟩ dt = In

N  

⟨∂t uτ , wτ ⟩ dt + In

n =1

N  ⟨[uτ ]n−1 , wτ (tn+−1 )⟩ + ⟨uτ (t0+ ), wτ (t0+ )⟩, n =2

for functions uτ , wτ ∈ Y (W ), which finally leads us to the inner products in (3.1) of the jumps [uτ ]n−1 at the discrete time points tn−1 with the test functions wτ (tn+−1 ). The dG(r) approach belongs to the class of Galerkin methods since the spaces of trial and test functions coincide. Now, we follow the derivation of Section 2.1. As in (2.16), we choose trial and test basis functions of Lagrange type ϕn,j , ψn,i ∈ Pr (In , R) on In with respect to r + 1 Gaussian quadrature nodes tn,j ∈ (tn−1 , tn ) for j = 1, . . . , r + 1. In contrast to Section 2.1 the (global) trial basis functions are also discontinuous now. The semi-discrete approximations {uτ , qτ } ∈ Yr (W ) × Yr (V ) then admit the representation r

uτ |In (t ) =

r +1 

Unj ϕn,j (t ) and qτ |In (t ) =

j =1

r +1 

Qnj ϕn,j (t ).

(3.4)

j=1

This enables us to recast the variational problem (3.1), (3.2) in the following form: For n = 1, . . . , N find coefficient functions {Unj , Qnj } ∈ W × V , for j = 1, . . . , r + 1, such that r +1 

αˆ ij ⟨Unj , w⟩ + τn βˆ i,i ⟨∇ · Qni , w⟩ = τn βˆ i,i ⟨f (tn,i ), w⟩ + γˆi ⟨uτ (tn−−1 ), w⟩

(3.5)

j =1

⟨D−1 Qni , v ⟩ − ⟨Uni , ∇ · v ⟩ = 0

(3.6)

for i = 1, . . . , r + 1 and all {w, v } ∈ W × V , where uτ (tn−−1 ) ∈ W is defined by means of u− τ ,n−1 = uτ |In−1 (tn−1 ) if n ≥ 2

and

if n = 1.

u− τ ,n−1 = u0

To find Eq. (3.5), we again replaced the right-hand side function f by its Lagrange interpolate Πr f ∈ Pr (In ; L2 (Ω )) that is built with respect to the r + 1 Lagrange basis functions ϕn,j . Further, we used that



d In

dt

ϕn,j (t ) · ψn,i (t ) dt + ϕn,j (tn−1 ) · ϕn,i (tn−1 ) = =



d

ˆI dtˆ r +1 

ϕˆ n,j (tˆ) · ψˆ n,i (tˆ) dtˆ + ϕˆ n,j (0) · ψˆ n,i (0)

ωˆ m ·

m=1

= ωˆ i · 

ϕn,j (t ) · ψn,i (t ) dt = τn In

 ˆI

ϕˆ n,j (tˆ) · ψˆ n,i (tˆ) dtˆ = τn ·

= τn · ωˆ i · δi,j := τn · βˆ i,j , ϕn,i (tn+−1 )

= ϕˆ i (0) =: γˆi

for i, j = 1, . . . , r + 1.

d dtˆ

r +1 

d dtˆ

ϕˆ j (tˆm ) · ψˆ i (tˆm ) + ϕˆ n,j (0) · ψˆ n,i (0)

ϕˆ j (tˆi ) + ϕˆ n,j (0) · ψˆ n,i (0) =: αˆ i,j ,

(3.7)

ωˆ m · ϕˆ j (tˆm ) · ψˆ i (tˆm )

m=1

(3.8) (3.9)

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

219

The discretization in space of the variational problem (3.5), (3.6) is now done in the same way as for the cGP(r) family of continuous Galerkin approximations in Section 2.2. By means of the representations in (2.27) of the coefficient functions j j Un ∈ W and Qn ∈ V with respect to the spatial finite element basis functions we end up with the following algebraic form of j j N N the fully discrete dG(r)-MFEM approach: For n = 1, . . . , N find coefficient vectors {un , qn } ∈ R Wh × R Vh , for j = 1, . . . , r + 1, of the expansion (2.27), (2.28) such that r +1 

αˆ ij Mujn + τn βˆ i,i B⊤ qin = τn βˆ i,i Fni + γi Mu− n−1 ,

(3.10)

j=1

Aqin − Buin = 0

(3.11)

for i = 1, . . . , r + 1, where u− n −1 ∈ R

NW

r +1 

u− n−1 =

j

un−1 ϕn−1 (tn−1 )

h

is defined by means of

if n ≥ 2

and

NW h 

Ph u0 =

u− n−1,l wh,l (x) if n = 1,

(3.12)

l =1

j =1

where Ph again denotes the L2 projection onto Wh . Thus, for n = 1 the vector u− n−1 is the representation of Ph u0 with respect to the finite element basis of Wh . 3.2. Solution of the algebraic system for the dG(1)-MFEM approach We shall now consider the system of Eqs. (3.10), (3.11) for the parameter r = 1 where the discrete solution is constructed by means of discontinuous piecewise polynomial functions in time of degree 1. Then, we have to compute two unknown coefficient vectors in each time interval In for the scalar- and the vector-valued variable, respectively. Putting j j N N r = 1 in (3.10), (3.11) yields: For n = 1, . . . , N find coefficient vectors {un , qn } ∈ R Wh × R Vh , for j = 1, 2, such that



τn βˆ i,i B⊤

0

αˆ 1,1 M

   

0

τn βˆ i,i B⊤

A

0

αˆ 2,1 M −B

0

A

0

with u0n ∈ R

NW

  1

 τn βˆ i,i Fn1 + γˆ1 Mu− n −1    2  αˆ 2,2 M  qn  τn βˆ i,i Fn2 + γˆ2 Mu− n −1     1 =    0  un  0

αˆ 1,2 M

−B

qn



u2n

(3.13)

0

being defined by (3.12).

h

The coefficients αˆ i,j , βˆ i,j and γˆi in (3.13) are given by Eqs. (3.7)–(3.9) as

√ αˆ 1,1 = 1,

αˆ 1,2 =

3−1 2

,

αˆ 2,1 =

√ − 3−1 2

,

αˆ 2,2 = 1

and

√ βˆ 1,1 =

1 2

,

βˆ 1,2 = 0,

βˆ 2,1 = 0,

βˆ 2,2 =

1 2

,

γˆ1 =

1+ 2

√ 3

,

γ2 =

1− 2

3

.

The 4 × 4 block systems are solved in the same way as described in Section 2.3 for the cGP(2)-MFEM approach. Elimination of the unknown coefficient vectors of the vector-valued variable yields the reduced 2 × 2 block system

  

αˆ 1,1 M +

τn 2

B⊤ A−1 B

αˆ 2,1 M

 τ    n 1 0 1 αˆ 1,2 M F + γ ˆ Mu 1 n n  un 2  , τn ⊤ −1  u2 =  τn 2 0 αˆ 2,2 M + B A B Fn + γˆ2 Mun n 2

(3.14)

2

whereas the flux unknowns can be calculated in a postprocessing step from Aq1n = Bu1n

and

Aq2n = Bu2n .

(3.15)

Finally, we solve the algebraic system (3.14) by the Schur complement approach described in Section 2.3. The solution of the algebraic form (3.13) of the dG(1)-MFEM approach is then determined by solving for n = 1, . . . , N the sequence of subproblems:

220

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224













⃗ u⃗0n + αˆ 1,1 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ M ⃗ −1 τn F⃗n2 + γˆ2 M ⃗ u⃗0n , ⃗2n = (−αˆ 2,1 ) τ2n F⃗n1 + γˆ1 M S⃗u 2 2 ⃗ −1 ⃗1n = − αˆ 1 M u 2,1



⃗ + αˆ 2,2 M

τ n ⃗ ⊤ ⃗ −1 ⃗ 2

B A



⃗2n − B u



τn 2

⃗ u⃗0n F⃗n2 + γˆ2 M



,

⃗ q⃗1n = B⃗u⃗1n , A ⃗ q⃗2n = B⃗u⃗2n , A

(3.16)









⃗ −1 αˆ 2,2 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ ⃗ + αˆ 1,1 M ⃗ + τn B⃗⊤ A⃗ −1 B⃗ M with S⃗ = αˆ 1,2 (−αˆ 2,1 )M 2 2 and

u0n

⃗ =

2 

j un − 1

⃗

ϕn−1,j (tn−1 ) if n ≥ 2 and Ph u0 =

NW h 

u0n,l wh,l (⃗ x) if n = 1.

l =1

j =0

The difference in the construction of the preconditioner (2.38), Kµ∗ = µ∗ M +

τn 2

B⊤ A−1 B,

for the Schur complement S then comes through the choice of the parameter µ∗ . For the case of the dG(1)-MFEM approximation with S being defined in (3.16) we find in the same way as sketched in Section 2.3 that

√ µ∗ =

6 2

√

with κ(S∗ ) ≤ 6 − 2 6

for the left-sided preconditioned operator S∗ = Pµ−∗1 S; cf. [24,34]. The rest of the solution strategy by nested iterations and the application of the preconditioner remains unchanged compared to the cGP(2)-MFEM approach presented in Section 2.3. 4. Numerical results In this section we study numerically the proposed variational time discretization schemes combined with mixed finite element methods in space. In Section 4.1 the experimental order of convergence is analysed. Since we focus in this work on the discretization of the time variable, we restrict ourselves on analysing the discretization error in time. For further convergence studies including the spatial convergence behaviour we also refer to [16,24]. In Section 4.2 we illustrate the robustness of our approach by simulating three-dimensional diffusive transport in a strongly anisotropic and heterogeneous medium. By this application we mimic flow and transport processes in composite soil formations that are studied in hydrological and environmental engineering. For further numerical tests we refer to [24]. Our schemes have been implemented in the open-source finite element library deal.ii [32]. The code can be run in parallel on distributed processors. For grid partitioning and parallel linear algebra we use the p4est [36] and Trilinos [37] libraries. For further details of our implementation we refer to [24]. 4.1. Convergence studies As test problem we consider our model problem (2.1)–(2.3) with D = I and Ω = (0, 1)2 . The prescribed solution is u(x, y, t ) = e2π t · x(1 − x)y(1 − y).

(4.1)

The corresponding flux variable q = −D∇ u and the source term f are then derived from the exact solution. The initial data is u0 (x, y) = u(x, y, 0). Since we focus here on higher order temporal discretizations, we restrict our numerical experiments to considering the cGP(2)-MFEM approach (2.37) and the dG(1)-MFEM method (3.16). In both schemes two coefficient vectors have to be calculated in each time of the intervals In such that the methods are of the same computational complexity. However, they have different stability properties. Whereas the cGP(2) approach is known to be A-stable [17], the dG(1) method is strongly A-stable cf. [26,27] and has, thereby, better damping properties with respect to high frequency error components. This is noted here, since the damping properties might be of importance for applications that are of high complexity and sensitivity in the temporal evolution [16,24]. For convergence studies of the cGP(1) approach we refer to [16] where this scheme is applied to wave problems and carefully studied. Since the dG(0) method corresponds to the well known implicit Euler method of first order accuracy, it is not in our scope of interest here. Moreover, combined implicit Euler and mixed finite element discretizations have often been studied in the literature yet. For the discretization in space we choose the pair {Wh , Vh } of mixed finite element spaces defined in (2.21) and (2.22), respectively, with parameter p = 2, i.e. Wh = Wh2 and Vh = Vh2 . Therefore, the exact solution {u(t ), q(t )} is an element of the mixed finite element spaces Wh2 × Vh2 such that the full discretization error is, up to possible quadrature errors, equal to the time discretization error u(t ) − uτ (t ). For this we note that {u(t ), q(t )} ∈ Wh × Vh implies that the semi-discretization

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

221

Table 2 Error norms and experimental order of convergence (EOC) for the cGP(2)-MFEM approximation of problem (2.1)–(2.3) with prescribed solution (4.1).

τ

∥q − qτ ∥L2 (L 2 )

EOC

∥u − uτ ∥L2 (L2 )

EOC

∥q − qτ ∥l∞ (L 2 )

EOC

∥u − uτ ∥l∞ (L2 )

EOC

1.000e−01 5.000e−02 2.500e−02 1.250e−02 6.250e−03 3.125e−03 1.563e−03 7.813e−04 3.906e−04 1.953e−04 9.766e−05

2.5986e−02 3.4005e−03 4.2756e−04 5.3064e−05 6.6017e−06 8.2381e−07 1.0293e−07 1.2865e−08 1.6081e−09 2.0104e−10 2.6019e−11

– 2.93 2.99 3.01 3.01 3.00 3.00 3.00 3.00 3.00 2.95

5.6667e−03 7.4019e−04 9.3559e−05 1.1727e−05 1.4668e−06 1.8338e−07 2.2923e−08 2.8654e−09 3.5817e−10 4.4774e−11 5.6359e−12

– 2.94 2.98 3.00 3.00 3.00 3.00 3.00 3.00 3.00 2.99

8.8605e−02 8.1198e−03 7.5003e−04 5.9011e−05 4.0566e−06 2.5592e−07 1.6090e−08 1.0069e−09 6.1861e−11 9.7130e−12 2.0072e−11

– 3.45 3.44 3.67 3.86 3.99 3.99 4.00 4.02 2.67 −1.05

1.6441e−02 1.1326e−03 7.3454e−05 4.6686e−06 2.9367e−07 1.8382e−08 1.1495e−09 7.1690e−11 4.4128e−12 1.5324e−12 1.9739e−12

– 3.86 3.95 3.98 3.99 4.00 4.00 4.00 4.02 1.53 −0.37

Table 3 Error norms and experimental order of convergence (EOC) for the dG(1)-MFEM approximation of problem (2.1)–(2.3) with prescribed solution (4.1).

τ

∥q − qτ ∥L2 (L 2 )

EOC

∥u − uτ ∥L2 (L2 )

EOC

∥q − qτ ∥l∞ (L 2 )

EOC

∥u − uτ ∥l∞ (L2 )

EOC

1.000e−01 5.000e−02 2.500e−02 1.250e−02 6.250e−03 3.125e−03 1.563e−03 7.813e−04 3.906e−04 1.953e−04 9.766e−05

4.4915e−01 1.2375e−01 3.2421e−02 8.2954e−03 2.0987e−03 5.2802e−04 1.3247e−04 3.3179e−05 8.3033e−06 2.0769e−06 5.1938e−07

– 1.86 1.93 1.97 1.98 1.99 1.99 2.00 2.00 2.00 2.00

1.0035e−01 2.7651e−02 7.2393e−03 1.8503e−03 4.6759e−04 1.1752e−04 2.9457e−05 7.3738e−06 1.8447e−06 4.6132e−07 1.1535e−07

– 1.86 1.93 1.97 1.98 1.99 2.00 2.00 2.00 2.00 2.00

5.5554e−01 8.7719e−02 1.3552e−02 2.1334e−03 3.4149e−04 5.3954e−05 8.1517e−06 1.1627e−06 1.5762e−07 2.0625e−08 2.6431e−09

– 2.66 2.69 2.67 2.64 2.66 2.73 2.81 2.88 2.93 2.96

1.1445e−01 1.6559e−02 2.2249e−03 2.8926e−04 3.7046e−05 4.7127e−06 5.9732e−07 7.5452e−08 9.4980e−09 1.1931e−09 1.5013e−10

– 2.79 2.90 2.94 2.96 2.97 2.98 2.98 2.99 2.99 2.99

in space {uh (t ), qh (t )} is equal to the exact solution {u(t ), q(t )}. If the time discretization is then applied to {uh (t ), qh (t )}, it follows that {uh,τ (t ), qh,τ (t )} = {uτ (t ), qτ (t )}. Since space and time discretizations are of Galerkin type, they commute and we obtain that {uτ ,h (t ), qτ ,h (t )} = {uh,τ (t ), qh,τ (t )} = {uτ (t ), qτ (t )}. We apply the time discretization schemes cGP(2) and dG(1) with an equidistant time step size τ = T /N. To measure the error, the following time discrete L∞ -norm of a function v : I → L2 (Ω ) is used

∥v∥l∞ (L2 ) = max ∥v − (tn )∥L2 (Ω ) , 1≤n≤N

v − (tn ) = lim v(t ), t → tn − 0

tn = nτ .

The behaviour of the standard L2 -norm ∥ · ∥L2 (L2 ) = ∥ · ∥L2 (I ;L2 (Ω )) and the time discrete L∞ -norm ∥ · ∥l∞ (L2 ) of the time discretization errors u(t )−uτ (t ) and q(t )−qτ (t ) of the scalar- and the vector-valued variable over the time interval I = [0, 1] are summarized in Tables 2 and 3, respectively. Moreover, the experimental order of convergence (EOC) is presented. For the cGP(2) method with piecewise quadratic polynomials in time third order convergence with respect to the standard L2 -norm is expected. For the dG(1) approach based on a piecewise linear approximation in time we expect second order of convergence. These rates of convergence are nicely recovered in Table 2 and Table 3, respectively, for the scalar- and the vector-valued variable. Thus, from the point of view of convergence rates the cGP(2) method is superior to the dG(1) approach since convergence of order 3 instead of order 2 is obtained with comparable computational costs. In the time discrete L∞ -norm a superconvergence behaviour is observed for either schemes and both variables, the solution u and the flux field q, with superconvergence of order 4 for the cGP(2) method and of order 3 for the dG(1) approach, as it is shown in Tables 2 and 3. This superconvergence behaviour of the cGP(r + 1) and dG(r ), r ≥ 1, families of variational time discretization schemes was already depicted in former works; cf. e.g. [16,17,20,21,26]. 4.2. Diffusion in heterogeneous and anisotropic media Finally, we shall illustrate the application of our families of approximation schemes to a more sophisticated test problem that is of practical interest. Thereby, we shall demonstrate the potential of the approaches and show that we reached the point that the methods can be put into use. For this, we choose a test configuration that is typical for the simulation of flow and transport processes in layered porous media, for instance in the subsurface. Such processes are a prominent class of problems for the application of mixed finite element methods. Here (local) mass conservation properties of the numerical schemes are of high importance. In particular, for such type of applications mixed finite element methods have shown their superiority and value over non-mixed continuous approximations. We study the model problem (2.1)–(2.3). The only modification now is that we use a time-dependent Dirichlet boundary condition that is prescribed below in Eq. (4.4). The computational domain Ω = (0, 1)3 and the definition of the anisotropic and heterogeneous diffusion tensor D are sketched in Fig. 4.1. In particular, a cross-section plot along the z-axis is adumbrated, i.e. the cross section plane with z = const. is

222

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

Fig. 4.1. Cross section plot in z-direction of the domain Ω = (0, 1)3 with layers’ dimensions and specification of the heterogeneous diffusion parameter and angles of the rotation matrices (4.3) in the diffusion matrix (4.2).

shown. The domain is built with three layered blocks with smaller diffusion properties that are given by the shaded regions K2 . In the two regions K1 and K2 (cf. Fig. 4.1) the heterogeneous diffusion tensor D is anisotropic and defined by

 D(x) = Qx Qy Qz ⊤

⊤

⊤

kx (x) 0 0

0 ky (x) 0

0 0 Qz Qy Qx kz (x)



(4.2)

with rotation matrices

 Qx =

1 0 0

0 cos(θ ) − sin(θ )

0 sin(θ ) , cos(θ )



 Qy =

cos(ψ) 0 − sin(ψ)

0 1 0

sin(ψ) 0 cos(ψ)



cos(φ) − sin(φ) 0

 Qz =

sin(φ) cos(φ) 0

0 0 . 1



(4.3)

In Eq. (4.2), we rotate the diffusion tensor such that its principle directions of diffusive transport are different from the coordinate lines and are not aligned with the faces of the finite element decomposition of the domain Ω into hexahedral elements. This choice of the diffusion matrix increases the numerical complexity of the problem and might be an additional source of numerical inaccuracies and diffusion. A steady two-dimensional counterpart of our test setting was introduced in [38]. The final time of the simulation is T = 3 and the initial value is u0 (x) = 0. We prescribe the inhomogeneous Dirichlet boundary condition g (x, t ) = min{t , 1.0}((1 − x)yz − x(1 − y)(1 − z )),

(4.4)

such that die scalar variable increases from zero to one in the corner P1 (cf. Fig. 4.1) within the time interval [0, 1] and is kept equal to one for t > 1, whereas the scalar variable is kept equal to zero in P2 . Thereby we model a transport process that is driven, for instance, by a temperature or pressure difference on the domain’s boundary. Exemplarily, we use the continuous cGP(2) method (2.32) of the proposed classes of time discretization schemes. The mixed finite element approach with p = 1 in (2.21) and (2.22), respectively, is applied for the approximation in space. The discretization parameters in time and space are chosen as τ = 10−2 and h = 10−1 . The simulation was done on a parallel machine with 2 octo-core CPUs by using 16 processes. For the scalar variable u and the vectorial variable q the number of degrees or the dimension of the unknown coefficient vectors, respectively, (cf. Eqs. (2.28) and (2.23)) amounts to NWh = 8000 and NVh = 25 200, thus in total 33 200 degrees of freedom per unknown coefficient of the temporal expansion (2.27). For the inversion of the mass matrix A in our solution algorithm (2.37) we used the parallel direct solver SuperLU [32,36]. The wall-clock time of the simulation was 57 min (3472.01 s). For further studies of the CPU time and the parallel scalability we refer to [24]. In Fig. 4.2 we visualize the calculated profiles of the scalar-valued solution u of the initial–boundary value problem and of the modulus of the vector-valued flux variable q. Further, the approximation of the flux field q is shown by glyphs and arrows. In the visualized profiles we recognize that the layers are nicely resolved, which in particular becomes obvious in the flux field. No smearing or artificial diffusion of the solution is observed which is a characteristic feature of the mixed finite element approach and underlines its superiority for simulating transport in anisotropic and heterogeneous media. This was already pointed out in former works, cf. e.g. [38]. Here, we have demonstrated the applicability of the developed families of variational time discretization schemes to problems of practical interest. The superiority and potential of the numerical approaches become in particular evident for solutions with high temporal fluctuations and for the case of a dominating temporal discretization error. Summary In this work we developed two classes of variational time discretization schemes for mixed finite element approximations of diffusive transport. In numerical studies the expected theoretical convergence rates could be verified for both variables. Moreover, a superconvergence behaviour with respect to the temporal approximation of the suggested methods could be shown. In a three-dimensional numerical simulation the applicability of the scheme to problems of practical interest was demonstrated. Moreover, the solution of the arising algebraic systems of equations was carefully described which is an important issue in the application of the proposed variational time discretization schemes. Our strategy is based on a Schur

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

223

Fig. 4.2. Higher order variational time discretization and mass conservative mixed finite element approximation of diffusive transport: Solution (left), modulus of flux (centre) and flux field (right) for t = 3.

complement technique along with an elimination of internal degrees of freedom in order to reduce the size of the algebraic system of equations and to overcome the indefiniteness of the system matrix. Nested preconditioned conjugate gradient iterations then become feasible. A hybridization of the mixed finite element approach was not studied in this work. This was mainly done due to the algorithmic complexity of this approach and former experiences of one of the authors with hybrid mixed finite element techniques [3,4,29]. Summarizing, the suggested methods offer large potential for further developments and for the highly accurate numerical approximation of time dependent transport problems, in particular if strong temporal variations and fluctuations are present. Acknowledgement This work was supported by the German Academic Exchange Service (DAAD) under the grant ID 56435737. References [1] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface: A Contribution to the Modeling of Hydrosystems, Springer, Berlin, 1997. [2] T. Arbogast, M.F. Wheeler, N.Y. Zhang, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal. 33 (1996) 1669–1687. [3] M. Bause, P. Knabner, Computation of variably saturated subsurface flow by adaptive mixed hybrid finite element methods, Adv. Water Resour. 27 (2004) 565–581. [4] M. Bause, Higher and lowest order mixed finite element methods for subsurface flow problems with solutions of weak regularity, Adv. Water Resour. 31 (2008) 370–382. [5] M.A. Celia, E.T. Bouloutas, R.L. Zarba, A general mass-conservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26 (7) (1990) 1483–1496. [6] Z. Chen, G. Huan, Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, Philadelphia, 2005. [7] F.A. Radu, I.S. Pop, P. Knabner, Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards’ equation, SIAM J. Numer. Anal. 42 (2004) 1452–1478. [8] F.A. Radu, I.S. Pop, S. Attinger, Analysis of an Euler implicit-mixed finite element scheme for reactive solute transport in porous media, Numer. Methods Partial Differential Equations 26 (2010) 320–344. [9] C.S. Woodward, C.N. Dawson, Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, SIAM J. Numer. Anal. 37 (2000) 701–724. [10] M. Bause, J. Hoffmann, P. Knabner, First-order convergence of multi point flux approximation on triangular grids and comparison with mixed finite element methods, Numer. Math. 116 (2010) 1–29. [11] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, New York, 1991. [12] Z. Chen, Finite Element Methods and their Applications, Springer, Berlin, 2010. [13] M.W. Farthing, C.E. Kees, C.T. Miller, Mixed finite element methods and higher-order temporal approximations, Adv. Water Resour. 25 (2002) 85–101. [14] A.K. Aziz, P. Monk, Continuous finite elements in space and time for the heat equation, Math. Comp. 52 (1989) 255–274. [15] S. Hussain, F. Schieweck, S. Turek, Higher order Galerkin time discretization for nonstationary incompressible flow, in: A. Cangiani, et al. (Eds.), Numer. Math. and Adv. Appl 2011, Springer, Berlin, 2013, pp. 509–517. [16] U. Köcher, M. Bause, Variational space–time methods for the wave equation, J. Sci. Comput. 61 (2) (2014) 424–453. [17] F. Schieweck, A-stable discontinuous Galerkin–Petrov time discretization of higher order, J. Numer. Math. 18 (1) (2010) 25–57. [18] V. Thomeé, Galerkin Finite Element Methods for Parabolic Problems, Springer, Berlin, 2006. [19] W. Bangerth, R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Birkhäuser, Basel, 2003. [20] S. Hussain, F. Schieweck, S. Turek, Higher order Galerkin time discretizations and fast multigrid solvers for the heat equation, J. Numer. Math. 19 (1) (2011) 41–61. [21] S. Hussain, F. Schieweck, S. Turek, A note on accurate and efficient higher order Galerkin time stepping schemes for nonstationary Stokes equations, Open Numer. Methods J. 4 (2012) 35–45. [22] N. Ahmed, G. Matthies, Numerical studies of variational-type time-discretization techniquesfor transient Oseen problem, in: Proceedings of Contributed Papers and Posters, Slovak University of Technology, Faculty of Civil Engineering, Department of Mathematics and Descriptive Geometry, Bratislava, 2012, pp. 404–415. ISBN: 978-80-227-3742-5/pbk. [23] W. Bangerth, M. Geiger, R. Rannacher, Adaptive Galerkin finite element methods for the wave equation, Comput. Methods Appl. Math. 10 (1) (2010) 3–48. [24] U. Köcher, Variational space–time methods for the elastic wave equation and the diffusion equation (Ph.D. thesis), Helmut-Schmidt-University, Hamburg, 2015, submitted.

224

M. Bause, U. Köcher / Journal of Computational and Applied Mathematics 289 (2015) 208–224

[25] N. Ahmed, G. Matthies, Numerical studies of Galerkin-type time-discretizations applied to transient convection–diffusion–reaction equations, World Acad. Sci. Eng. Technol. 66 (2012) 586–593. [26] G. Matthies, F. Schieweck, Higher order variational time discretizations for nonlinear systems of ordinary differential equations. Preprint No. 23/2011 Otto von Guericke Universität Magdeburg, 2011, pp. 1–30. [27] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, Stiff and Differential-Algebraic Problems, Springer, Berlin, 1996. [28] M. Bause, F. Radu, U. Köcher, Error analysis for discretizations of parabolic problems using continuous finite elements in time and mixed finite elements in space, 2014 (in preparation). [29] M. Bause, F. Brunner, P. Knabner, F. Radu, An improved optimal order mixed finite element method for semilinear transport problems, in: A. Cangiani, Andrea, et al. (Eds.), Numerical Mathematics and Advanced Applications 2011: Proceedings of ENUMATH 2011, Springer, Berlin, 2013, pp. 247–256. [30] F. Brunner, F. Radu, M. Bause, P. Knabner, Optimal order convergence of a modified BDM1 mixed finite element scheme for reactive transport in porous media, Adv. Water Resour. 35 (2012) 163–171. [31] L.C. Evans, Partial Differential Equations, American Mathematical Society, Providence, Rhode Island, 2010. [32] W. Bangerth, T. Heister, G. Kanschat, deal.II differential equations analysis library, Technical reference, 2014. http://www.dealii.org. [33] E. Bänsch, P. Morin, R.H. Nochetto, Preconditioning a class of fourth order problems by operator splitting, Numer. Math. 118 (2) (2011) 197–228. [34] Stephan Weller, Steffen Basting, Efficient preconditioning of variational time discretization methods for parabolic Partial Differential Equations, ESAIM: Math. Model. Numer. Anal. 49 (2) (2015) 331–347; cf. http://dx.doi.org/10.1051/m2an/2014036. [35] Y. Saad, Iterative Methods for Sparse Linear Systems, second ed., SIAM, Philadelphia, 2003. [36] M. Heroux, An overview of Trilinos, Sandia National Laboratories, SAND2003–2927, 2003. [37] C. Burstedde, L.C. Wilcox, O. Ghattas, p4est: scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM J. Sci. Comput. 33 (3) (2011) 1103–1133. [38] R. Klausen, On locally conservative numerical methods for ellipticm problems; application to reservoir simulation (Ph.D. thesis, Dissertation 297), Unipub AS, University of Oslo, 2003.