On conservation and dual
consistency for
summation-by-parts based
approximations of parabolic
problems
Fatemeh Ghasemi and Jan Nordström
On conservation and dual consistency for
summation-by-parts based approximations of parabolic
problems
Fatemeh Ghasemi & Jan Nordstr¨om
Computational Mathematics, Department of Mathematics, Link¨oping University, SE-581 83 Link¨oping, Sweden
fatemeh.ghasemi@liu.se, jan.nordstrom@liu.se
1. Introduction
We consider the coupling of parabolic problems discretized using differ-ence operators on summation-by-parts (SBP) form with interface conditions imposed weakly. In [1, 2], it was shown that conservation and dual consis-tency are equivalent concepts for linear conservation laws. Here, we show that these concepts are equivalent also for symmetric or symmetrizable parabolic problems, exemplified by the heat equation. We rewrite the heat equation as first order system as in the local discontinuous Galerkin method and show the equivalence of dual consistency and conservation for both the first and second order forms [3, 4].
2. The continuous problem
Consider the coupling of two heat equations
ut= (ux)x, −1 ≤ x ≤ 0, t > 0, vt= (vx)x, 0 ≤ x ≤ 1, t > 0, u(0, t) = v(0, t), x = 0, t > 0, ux(0, t) = vx(0, t), x = 0, t > 0,
(1)
augmented with initial conditions. In (1), u(x, t) and v(x, t) are the solutions in the left and right domain, respectively. Furthermore, (x, t) > 0 is the thermal diffusivity which is continuous across the interface x = 0. For clarity and ease of presentation, the boundary conditions at x = ±1 are ignored.
2 THE CONTINUOUS PROBLEM 2
The problem (1) is reduced to first order form by introducing the new variables ˜u = ux, ˜v = vx ˜ IUt+ AU + BUx =0, −1 ≤ x ≤ 0, t > 0, ˜ IVt+ AV + BVx =0, 0 ≤ x ≤ 1, t > 0, U (0, t) =V (0, t), x = 0, t > 0, (2)
augmented with homogeneous initial conditions. In (2), U = [u, ˜u]T, V = [v, ˜v]T, ˜I = diag(1, 0), A = diag(0, 1/(x, t)) and B = − 0 1
1 0
. 2.1. Conservation
Let ϕ(x, t) be a smooth function with compact support in −1 ≤ x ≤ 1. Multiplying (1) by ϕ(x, t) and performing integration by parts in space twice, leads to Z 0 −1 ϕutdx + Z 1 0 ϕvtdx = Z 0 −1 (ϕx)xudx + Z 1 0 (ϕx)xvdx +(0, t) ϕ(0, t) ux(0, t) − vx(0, t) − ϕx(0, t) u(0, t) − v(0, t) . (3)
The interface conditions in (1) remove the interface terms in (3) and conser-vation follows.
Multiplying (2) with the transpose of the compactly supported and smooth vector function Φ(x, t), and integrating in space leads to
Z 0 −1 ΦTIU˜ tdx + Z 1 0 ΦTIV˜ tdx + Z 0 −1 ΦTAU dx + Z 1 0 ΦTAV dx = − ΦT(0, t)B U (0, t) − V (0, t) + Z 0 −1 ΦTxBU dx + Z 1 0 ΦTxBV dx. (4)
The interface conditions in (2) remove the interface terms in (4), which (dis-regarding the zero order terms) implies that (2) is on conservative form. Remark 1. In this paper, we use a broad definition of conservation where the equations are multiplied by a smooth test function. By performing integration-by-parts we transfer all derivatives to the test function. We say that the prob-lem is conservative if all interface terms vanish. This definition is motivated by the comparison with the problem on first order form.
3 THE SEMI-DISCRETE APPROXIMATION 3
2.2. The dual problem
The dual problem [5] associated to (1) is
φτ = (φx)x, −1 ≤ x ≤ 0, τ > 0, ψτ = (ψx)x, 0 ≤ x ≤ 1, τ > 0, φ(0, τ ) = ψ(0, τ ), x = 0, τ > 0, φx(0, τ ) = ψx(0, τ ), x = 0, τ > 0,
(5)
augmented with homogeneous initial conditions at τ = 0. In (5), the dual time variable is τ = T − t, where T is the final time for the primal problem. We ignore the boundary conditions at x = ±1 as in the primal problem. In the same way as for the primal problem, (5) can be shown to be conservative. The dual problem associated to (2) is
˜ IΦτ + AΦ − BΦx =0, −1 ≤ x ≤ 0, τ > 0, ˜ IΨτ + AΨ − BΨx =0, 0 ≤ x ≤ 1, τ > 0, Φ(0, t) =Ψ(0, t), x = 0, τ > 0, (6)
augmented with homogeneous initial conditions at τ = 0. In (6), Φ = [φ, ˜φ]T and Ψ = [ψ, ˜ψ]T. Conservation of (6) is shown in the same way as for (2).
3. The semi-discrete approximation
To define a semi-discretization of (1), we discretize the left and right domains by N + 1 and M + 1 grid points, respectively. Let the vectors u, v, U and V contain approximations of u, v, U and V , respectively. Then, the first derivative is approximated using SBP operators
ux ≈DLu = PL−1QLu, vx ≈DRv = PR−1QRv, Ux ≈(DL⊗ I2)U, Vx ≈(DR⊗ I2)V,
(7) where ⊗ denotes the Kronecker product and the subscripts L, R refer to the left and right spatial intervals, respectively. For simplicity, we use the first derivative twice as our second derivative operator. In (7), PL,R are symmetric positive definite matrices and QL,R are almost skew-symmetric matrices satisfying QL,R+ QTL,R= diag(−1, 0, ..., 0, 1). I2 is the 2 × 2 identity matrix.
3 THE SEMI-DISCRETE APPROXIMATION 4
Now, we discretize (1) in space as in [6] and obtain ut− DL(LDLu) =PL−1σ 1 L(uN − v0) + σL2 (DLu)N − (DRv)0eN +PL−1DTLσ3 L(uN − v0) + σL4 (DLu)N − (DRv)0eN, vt− DR(RDRv) =PR−1σ 1 R(v0 − uN) + σR2 (DRv)0− (DLu)Ne0 +PR−1DTRσ3R(v0 − uN) + σR4 (DRv)0− (DLu)Ne0, (8)
where the diagonal matrices L and R approximate (x, t) point wise on the left and right domains, respectively. In (8), the vectors eN = [0, · · · , 0, 1]T and e0 = [1, 0, · · · , 0]T have the length of the left and right mesh, respectively. The scalars uN, (DLu)N, v0 and (DRv)0 are approximations of u, ux, v and vx at x = 0. Furthermore, the penalty terms on the right-hand side (RHS) impose the interface conditions weakly.
Next, we discretize (2) as
(IN ⊗ ˜I)Ut+ [AL+ DL⊗ B]U = (PL−1eN ⊗ Σl)(UN − V0), (IM ⊗ ˜I)Vt+ [AR+ DR⊗ B]V = (PR−1e0⊗ Σr)(V0− UN),
(9)
where IN and IM are the identity matrices of size N + 1 and M + 1 and AL and AR approximate A = A(x, t) point wise on the left and right domains, respectively.
We can rearrange the terms in (8) such that we get a first order scheme, which is equivalent to (9) if Σl= " σ1 L+ αLσ3L− αL β (σ 2 L+ αLσL4)(σL3 + σR3) −1 β (σ 2 L+ αLσ4L) −σ3 L I + αL βIσ 4 L(σL3 + σR3) 1 βIσ 4 L # , Σr= " σ1 R− αRσ3R− αR β (σ 2 R− αRσ4R)(σ3L+ σ3R) −1 β (σ 2 R− αRσR4) −σ3 R I + αR βIσ 4 R(σ 3 L+ σ 3 R) 1 βIσ 4 R # , (10) where I = (0, t), αL = eTNP −1 L eN, αR = eT0P −1 R e0 and β = −I + αLσL4 + αRσ4R.
3 THE SEMI-DISCRETE APPROXIMATION 5 is equivalent to (8) with Σl= σ1 l σ2l σ3 l σ4l and Σr = σ1 r σr2 σ3 r σr4 if σL1 =σl1+ αLIσl3− αLI ˜ β (σ 2 l + αLIσl4)(σ 3 l + σ 3 r), σ 2 L= − I ˜ β(σ 2 l + αLIσ4l) σL3 = − Iσl3+ αL2I ˜ β σ 4 l(σ 3 l + σ 3 r), σ 4 L= 2 I ˜ βσ 4 l, σR1 =σr1− αRIσr3− αRI ˜ β (σ 2 r − αRIσ4r)(σ 3 l + σ 3 r), σ 2 R= − I ˜ β(σ 2 r − αRIσr4) σR3 = − Iσr3+ αR2I ˜ β σ 4 r(σ 3 l + σ 3 r), σ 4 R= 2I ˜ βσ 4 r, (11) where ˜β = −1 + αLIσ4l + αRIσl4. 3.1. Conservation
Let ϕϕϕL,R denote smooth grid functions at the left and right domains. Multiplying (8) by ϕϕϕT
LPL and ϕϕϕTRPR and using ϕϕϕNL = ϕϕϕ0R:= ϕϕϕI leads to ϕ ϕϕTLPLut+ ϕϕϕTRPRvt = [DLL(DLϕϕϕL)]TPLu + [DRR(DRϕϕϕR)]TPRv + IT, (12) where IT = ϕϕϕI(σL1 − σ 1 R)[uN − v0] + ϕϕϕI(σL2 − σ 2 R+ I)[(DLu)N − (DRv)0] +(DRϕϕϕR)0(σ3L− I− σR3)[uN − v0] + (DRϕϕϕR)0(σ4L− σ 4 R)[(DLu)N − (DRv)0] +(DLϕϕϕL)N − (DRϕϕϕR)0σ3L(uN − v0) − IuN + σL4 (DLu)N − (DRv)0.
The last term in IT is proportional to the order of the scheme. Following [7], we say that (8) is conservative to the order of the scheme if IT = 0, which holds if (DLϕϕϕL)N = (DRϕϕϕR)0 and σL1 = σR1, σL2 = σR2 − I, σL3 = σ 3 R+ I, σ4L= σ 4 R. (13)
On the other hand, if (8) is conservative to the order of the scheme, then IT = 0 for all possible solutions and test functions which implies that (13) must hold.
We summarize the result in
Proposition 1. The interface procedure in (8) is conservative if and only if (13) holds.
3 THE SEMI-DISCRETE APPROXIMATION 6
Remark 2. It can be shown that Proposition 1 holds for narrow stencil SBP operators that can accommodate variable coefficients as defined in [8].
Next, we multiply (9) with ΦΦΦT
L(PL⊗ I2) and ΦΦΦTR(PR⊗ I2) to get
ΦΦΦTL(PL⊗ ˜I)Ut+ ΦΦΦTR(PR⊗ ˜I)Vt+ ΦΦΦTL(PL⊗ I2)ALU + ΦΦΦTR(PR⊗ I2)ARV = ΦΦΦTI(Σl− Σr− B)[UN − V0] + ΦΦΦLT(DTLPL⊗ B)U + ΦΦΦTR(DTRPR⊗ B)V,
(14) which mimics the continuous formulation in (4). In (14), ΦΦΦI := (ΦΦΦL)N = (ΦΦΦR)0, and a conservative scheme requires
Σl= Σr+ B. (15)
On the other hand, the interface term on the RHS of (14) must vanish for all possible solutions and test functions if (9) is a conservative scheme, which implies that (15) holds.
We summarize the result in
Proposition 2. The interface procedure in (9) is conservative if and only if (15) holds.
Remark 3. If the scheme (8) is conservative, then the first order scheme (9) with penalty matrices given in (10), is conservative. Similarly, if the scheme (9) is conservative, then the second order scheme (8) with penalty coefficients given in (11), is conservative.
Remark 4. Since σ4
L = σ4R 6= 0 may increase the truncation error of the scheme (8), see [6], an obvious choice is σ4
L = σR4 = 0, which is used in the rest of the paper.
3.2. Dual Consistency
To determine dual consistency, we rewrite (8) in the compact form wt+ Lw = 0, L = −(D2+ P−1Σ), (16)
where w = [u, v]T, D2 = diag(DLLDL, DRRDR), P−1 = diag(PL−1, P −1 R ) and Σ = F (Σ1+ DΣ2) + DTF (Σ3+ DΣ4), F = +EN −EI −ET I +E0 . (17)
3 THE SEMI-DISCRETE APPROXIMATION 7
In (17), Σi = diag(σLi, σRi), where σL,Ri given in (8), D = diag(DL, DR), while EN = eNeTN and E0 = e0eT0. The (N + 1) × (M + 1) matrix EI = eNeT0 picks out the parts of the vectors residing at the interface such that EIv = [0, · · · , 0, v0]T and EITu = [uN, 0, · · · , 0]T.
The discrete dual problem corresponding to (16) is (see [5]) ϑ
ϑϑτ + P−1LTP ϑϑϑ = 0, (18) where ϑϑϑ = [φφφ, ψψψ]T. By using the SBP properties of the operators, we expand (18) and write it in component form as
φφφt−DL(LDLφφφ) = PL−1σ 1 LφφφN − σ1Rψψψ0+ (σ3L− I)(DLφφφ)N − σR3(DRψψψ)0eN +PL−1DTL(I+ σ2L)φφφN − σR2 ψψψ0eN, ψ ψ ψt−DR(RDRψψψ) = PR−1σ 1 Rψψψ0− σL1φφφN + (I + σR3)(DRψψψ)0− σL3(DLφφφ)Ne0 +PR−1DTR(σ2R− I)ψψψ0− σL2 φφφNe0. (19) Dual consistency requires that (19) is a consistent approximation of (5), which leads to
Proposition 3. The scheme (19) is dual consistent if only if (13) holds. Proof. The left-hand side (LHS) of (19) approximate the equations in (5). If (13) is satisfied, then the RHS of (19) imposes the dual interface conditions in (5). On the other hand, the interface terms on the RHS of (19) should vanish if ϕϕϕ and ψψψ satisfies the numerical interface conditions exactly. In other words, if (13) is not satisfied, the numerical scheme (19) would impose interface conditions at x = 0 that do not exist in the adjoint problem.
To determine dual consistency, we rewrite (9) in compact form ˜IW
t+ LW = 0, L = D2− P−1S, ˜I = diag(IN ⊗ ˜I, IM ⊗ ˜I), (20) and W = [U, V]. In (20), D2 = diag(AL+ DL⊗ B, AR+ DR⊗ B), and
S = +EN ⊗ Σl −EI ⊗ Σl −ET I ⊗ Σr +E0⊗ Σr , P−1 = diag(PL−1⊗ I2, PR−1⊗ I2).
The discrete dual problem corresponding to (20) is ˜
4 EQUIVALENCE OF CONSERVATION AND DUAL CONSISTENCY 8
where θθθ = [ΦΦΦ, ΨΨΨ]T. By expanding (21) in component form, we find (IN ⊗ ˜I)ΦΦΦt+ [AL−DL⊗ B]ΦΦΦ = (PL−1eN ⊗ (ΣTl − B T))ΦΦΦ N−(PL−1eN ⊗ ΣTr)ΨΨΨ0, (IM ⊗ ˜I)ΨΨΨt+ [AR−DR⊗ B]ΨΨΨ = (PR−1e0⊗ (ΣTr + B T))ΨΨΨ 0− (PR−1e0⊗ ΣTl )ΦΦΦN. (22) Proposition 4. The scheme (22) is dual consistent if only if (15) holds. Proof. The terms on the LHS in (22) approximate the two systems of equa-tions in (6). If (15) is satisfied, the RHS of (22) imposes the dual interface conditions in (6). On the other hand, the interface terms on the RHS of (22) should vanish if ΦΦΦ and ΨΨΨ satisfies the numerical interface conditions exactly. In other words, if (15) is not satisfied, the numerical scheme (22) would impose interface conditions at x = 0 that do not exist in the adjoint problem.
4. Equivalence of conservation and dual consistency We are now ready to state the main results.
Proposition 5. The approximation (16) is conservative if and only if it is dual consistent.
Proof. It follows immediately from Propositions (1) and (3).
Proposition 6. The approximation (20) is conservative if and only if it is dual consistent.
Proof. It follows immediately from Propositions (2) and (4).
Remark 5. Proposition 5 and 6 imply that the conservation relations in Remark 3 also lead to dual consistency.
5. References
[1] J. Nordstr¨om, F. Ghasemi, On the relation between conservation and dual consistency for summation-by-parts schemes, Journal of Computational Physics 344 (2017) 437–439.
[2] J. Nordstr¨om, F. Ghasemi, Corrigendum to “On the relation between con-servation and dual consistency for summation-by-parts schemes”, Journal of Computational Physics 360 (2018) 247.
5 REFERENCES 9
[3] B. Cockburn, C.-W. Shu, The local discontinuous Galerkin method for time-dependent convection-diffusion systems, SIAM Journal on Numeri-cal Analysis 35 (1998) 2440–2463.
[4] D. N. Arnold, F. Brezzi, B. Cockburn, L. D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis 39 (2002) 1749–1779.
[5] J. Berg, J. Nordstr¨om, Superconvergent functional output for time-dependent problems using finite differences on summation-by-parts form, Journal of Computational Physics 231 (2012) 6846–6860.
[6] M. H. Carpenter, J. Nordstr¨om, D. Gottlieb, Revisiting and extending in-terface penalties for multidomain summation-by-parts operators, Journal of Scientific Computing 45 (2010) 118–150.
[7] K. Mattsson, M. H. Carpenter, Stable and accurate interpolation oper-ators for high-order multiblock finite difference methods, SIAM Journal on Scientific Computing 32 (2010) 2298–2320.
[8] K. Mattsson, Summation by parts operators for finite difference ap-proximations of second-derivatives with variable coefficients, Journal of Scientific Computing 51 (2012) 650–682.
