Published online by De Gruyter June 22, 2022

# The DPG Method for the Convection-Reaction Problem, Revisited

Leszek Feliks Demkowicz , Nathan V. Roberts and Judit Muñoz-Matute

# Abstract

We study both conforming and non-conforming versions of the practical DPG method for the convection-reaction problem. We determine that the most common approach for DPG stability analysis – construction of a local Fortin operator – is infeasible for the convection-reaction problem. We then develop a line of argument based on a direct proof of discrete stability; we find that employing a polynomial enrichment for the test space does not suffice for this purpose, motivating the introduction of a (two-element) subgrid mesh. The argument combines mathematical analysis with numerical experiments.

MSC 2010: 35Q49; 65M60

Funding source: Sandia National Laboratories

Award Identifier / Grant number: 218322

Funding source: National Science Foundation

Award Identifier / Grant number: 1819101

Funding source: H2020 Marie Skłodowska-Curie Actions

Award Identifier / Grant number: 101017984

Funding statement: L. Demkowicz was partially supported with Sandia LDRD Project No. 218322 and NSF grant No. 1819101. N. Roberts was supported by Sandia LDRD Project No. 218322. J. Muñoz-Matute was supported by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie individual fellowship No. 101017984 (GEODPG). Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This is Sandia report number SAND2022-7111 J.

## A Boundedness Below of the Convection-Reaction Operator

Clearly, we have to make some assumptions about coefficients b i and 𝑐 to guarantee boundedness below of the operator with some reasonable lower bound on boundedness below constant 𝛼. If both 𝑏 and 𝑐 converge to zero, then α 0 as well. The stability has to come either from convection or reaction[6].

### Stability from Reaction

This can be accessed quickly by testing with 𝑢. We get

Ω b i u x i u + c u 2 = Ω f u .

Integrating the first term by parts,

Ω b i u x i u = Ω b i x i ( u 2 2 ) = - Ω div b u 2 2 + Γ + b n u 2 2 ,

we obtain

Γ + b n u 2 2 + Ω ( c - 1 2 div b ) u 2 = Ω f u .

Assuming

c - 1 2 div b β > 0 ,

we obtain

β u 2 f u β u f .

This sufficient condition clearly does not cover the pure advection case with div b = 0 .

This is a bit more difficult to analyze. We have to use a bit more sophisticated analysis based on characteristics and turning the advection-reaction problem into a family of ODEs.

### Characteristics

Solutions of the system of first-order nonlinear ODEs

d x d t = b ( x ( t ) )

are called characteristics. We shall assume that the family of characteristics can be extended to a curvilinear system of coordinates ( t , ξ ) covering the whole domain Ω; see Figure 11.

### Figure 11

Characteristics.

Each characteristic originates on inflow boundary Γ - and terminates on outflow boundary Γ + . We shall assume the parametrization

x = x ( t , ξ ) , ξ O , x ( t - ( ξ ) , t + ( ξ ) )

with the corresponding jacobian j ( x ) . If we furthermore assume that the system is orthogonal,

x ξ i x ξ j = δ i j x ξ i x t = 0 i , j = 1 , 2 ( 3D version ) ,

the jacobian equals the magnitude of the advection vector,

jac ( x ( t , ξ ) ) = | b ( x ( t , ξ ) ) ) | .

Implicit in the assumptions is that b ( x ) 0 in domain Ω. If we assume that b C 1 ( Ω ¯ ) , the Weierstrass Theorem implies that we must have a positive lower bound on | b ( x ) | .

### Convection-Reaction Problem with Homogeneous BC

We have

(A.1) { b u + c u = f in Ω , u = 0 on Γ - .

Let u ^ ( t ) = u ( x ( t ) ) , where x ( t ) is the characteristic. Then u ( x ) satisfies (A.1) if and only if u ^ ( t ) solves the initial-value ODE problem

(A.2) { d u ^ d t + c ^ u ^ = f ^ , t ( 0 , t x 0 ) , u ^ ( t - ) = 0 ,

where c ^ ( t ) := c ( x ( t ) ) , f ^ ( t ) := f ( x ( t ) ) . The linear ODE with variable coefficients (A.2) admits a closed form solution

u ^ ( t ) = t - t e - ( C ^ ( t ) - C ^ ( s ) ) f ^ ( s ) d s ,

where C ^ ( t ) is a primitive of c ^ ( t ) , i.e., d C ^ d t = c ^ . Standard reasoning and Cauchy–Schwarz inequality lead to the estimate

t - t + | u ^ ( t ) | 2 d t = t - t + | t - t e - ( C ^ ( t ) - C ^ ( s ) ) f ^ ( s ) d s | 2 d t t - t + t - t | e - ( C ^ ( t ) - C ^ ( s ) ) | 2 d s t - t | f ^ ( s ) | 2 d s d t t - t + t - t | e - ( C ^ ( t ) - C ^ ( s ) ) | 2 d s d t stability constant t - t + | f ^ ( s ) | 2 d s .

We need to work out now sufficient conditions on coefficients b ( x ) , c ( x ) to translate this estimate into a boundedness below estimate for the convection-reaction operator.

### Assumption on the Reaction Term

We shall assume that c ( x ) 0 . With this assumption,

C ^ ( t ) - C ^ ( s ) = s t c ( x ( η , ξ ) ) d η 0 e - ( C ^ ( t ) - C ^ ( s ) ) 1 ,

and the stability constant can be estimated by

t - t + t - t d s d t = ( t + - t - ) 2 2 .

Unfortunately, the 1D estimate does not translate immediately into the boundedness below condition as we have to include the jacobian in the computations,

t - t + | u ^ ( t ) | 2 jac ( x ( t , ξ ) ) = jac ^ ( t ) d t = t - t + | t - t e - ( C ^ ( t ) - C ^ ( s ) ) f ^ ( s ) d s | 2 jac ^ ( t ) d t
t - t + t - t | e - ( C ^ ( t ) - C ^ ( s ) ) | 2 d s t - t | f ^ ( s ) | 2 d s jac ^ ( t ) d t
t - t + ( t - t - ) t - t | f ^ ( s ) | 2 d s jac ^ ( t ) d t ( c 0 )
= t - t + s t + ( t - t - ) jac ^ ( t ) d t | f ^ ( s ) | 2 d s ( Fubini’s Theorem ) .
This leads to a rather technical assumption

s t + ( t - t - ) jac ^ ( t ) d t C jac ^ ( s ) , s < t < t + , t - < s < t + , C > 0 .

The assumption is easily satisfied if we assume lower and upper bounds on the jacobian,

0 < b min jac ^ = | b ^ | b max < .

This implies

jac ^ ( t ) b max b min jac ^ ( s ) ,

from which it follows that

s t + ( t - t - ) jac ^ ( t ) d t b max b min jac ^ ( s ) s t + ( t - t - ) d t b max b min ( t + - t - ) 2 2 jac ^ ( s ) .

## B Boundedness Below of the Convection-Reaction Operator. Second Approach

The disadvantage of the analysis shown in Appendix A is that it does not cover the case when the stability (boundedness below) may come from the advection in one part of the domain, and from reaction in the rest of the domain. The approach presented here rectifies this deficiency in the case when the advection field has a scalar potential, b ( x ) = V ( x ) . We shall also focus now on the adjoint operator

A v = - div ( b v ) + c v , D ( A ) = { v H A ( Ω ) : v = 0 on Ω + } .

Actually, we will develop a more general stability estimate for any v H A ( Ω ) . Following [14], we start by introducing an auxiliary unknown

w ( x ) := e V ( x ) v ( x ) , w = e V b v + e V v .

Let f := A v . We have

e V f = e V ( - b v + ( c - div b ) v ) = - div ( b w ) + ( | b | 2 + c ) w .

Multiplying both sides with 𝑤 and integrating over Ω, we obtain

- Ω div ( b w ) w + Ω ( | b | 2 + c ) w 2 = Ω e V f w .

The first term is now integrated by parts,

- Ω div ( b w ) w = Ω b ( w 2 2 ) - Γ b n w 2 = - 1 2 Γ b n w 2 - 1 2 Ω div b w 2 .

This gives

- 1 2 Γ b n w 2 + Ω ( | b | 2 + c - 1 2 div b = : a ) w 2 = Ω e V f w .

Under an additional assumption that coefficient a ( x ) is positive, we can use Young’s inequality to estimate the right-hand side,

f w a 2 w 2 + e 2 V 2 a f 2 .

This leads to the final estimate

1 2 Γ - | b n | w 2 + 1 2 Ω a w 2 Ω e 2 V 2 a f 2 + 1 2 Γ + b n w 2 .

In particular, for v = 0 on Γ + , we obtain

Ω a e 2 V v 2 Ω e 2 V a f 2 .

If 𝑎 admits a lower bound a min a and e 1 , e 2 are lower and upper bounds for e 2 V , we obtain

a min e 1 Ω v 2 Ω a e 2 V v 2 e 2 a min Ω f 2 .

This gives the final estimate for the boundedness below constant,

e 1 e 2 a min 2 Ω v 2 Ω f 2 .

Remark 11

(1) For an incompressible advection field div b = Δ V = 0 , a = | b | 2 + c . Thus both advection term | b | 2 and reaction coefficient 𝑐 may degenerate to zero, as long as the sum stays bounded away from zero.

(2) The assumption on coefficient 𝑎 to be bounded below away from zero is not as restrictive as it may appear. We can redefine w = e k V v and obtain a modified equation for 𝑤,

- div ( b w ) + ( k 2 | b | 2 + c ) w .

As long as div b is bounded, we can select a sufficiently large multiplier 𝑘 to ensure the positivity of coefficient 𝑎. Of course, a larger 𝑘 will result in a larger bound e 1 e 2 .

## C Local Stability Result

Let 𝐾 be an element with its outflow boundary denoted by K + . Let a ( x ) and w ( x ) be positive weights defined on 𝐾. Consider the constrained minimization problem

min v h P p + 1 ( K ) 1 2 { K w a ( A ( v h - v ) ) 2 + K + w b n ( v h - v ) 2 }

under the constraints

K δ u h A ( v h - v ) = 0 for all δ u h P p - 1 ( K ) , K + w b n δ w h ( v h - v ) = 0 for all δ w h P c p + 1 ( K + ) .

The first constraint is a Fortin-like condition enabling the replacement of the exact test function 𝑣 with its approximation v h in our stability analysis. The second constraint enforces weak conformity of the global discrete test function; see Section 6 for details. The minimization problem is equivalent to the mixed problem

{ v h P p + 1 ( K ) , u h P p - 1 ( K ) , w h P c p + 1 ( K + ) , K w a A v h A δ v h + K u h A δ v h + K + w b n ( v h + w h ) δ v h = K w a A v A δ v h + K + w b n v δ v h , δ v h P p + 1 ( K ) , K δ u h A v h = K δ u h A v , δ u h P p - 1 ( K ) , K + w b n δ w h v h = K + w b n δ w h v , δ w h P c p + 1 ( K + ) ,

or, equivalently,

(C.1) { v h P p + 1 ( K ) , u h P p - 1 ( K ) , w h P c p + 1 ( K + ) , K w a A v h A δ v h + K u h A δ v h + K + w b n w h δ v h = K w a A v A δ v h , δ v h P p + 1 ( K ) , K δ u h A v h = K δ u h A v , δ u h P p - 1 ( K ) , K + w b n δ w h v h = K + w b n δ w h v , δ w h P c p + 1 ( K + ) .

Let

V h := P p + 1 ( K ) , V h , 0 := { v h V h : K + w b n δ w h v h = 0 for all δ w h P c p + 1 ( K + ) } , V h , 00 := { v h V h , 0 : K δ u h A v h = 0 for all δ u h P p - 1 ( K ) } .

Note that, if the restriction of 𝑣 to K + is itself in space P c p + 1 ( K + ) , then V h , 0 = { v h V h : v h = 0 on K + } . Consider the corresponding decomposition of v h V h ,

v h = v h , 00 + v h , 0 + v h , v h , 00 V h , 00 , v h , 0 V h , 0 , v h V h ,

where

V h , 0 := { v h , 0 V h , 0 : ( v h , 0 , δ v h ) V = 0 for all δ v h V h , 00 } , V h := { v h V h : ( v h , δ v h ) V = 0 for all δ v h V h , 0 } ,

with the inner product

( v , δ v ) V = K w a A v A δ v + K + w b n v δ v .

Introduce norms for the Lagrange multiplier u h P h p - 1 ( K ) , w h P c p + 1 ( K + ) ,

u h K 2 := K u h 2 , w h K + 2 := K + w b n w 2 ,

and consider the corresponding inf-sup constants

α h := inf u h P p - 1 ( K ) sup v h V h , 0 Ω u h A v h u h K v h V , β h := inf w h P c p + 1 ( K + ) sup v h V h K + w b n w h v h w h K + v h V .

Lemma 8

The following estimate holds:

v h V 2 ( 1 + α h - 2 ) A v 2 + β h - 2 v K + 2 .

## Proof

Testing in (C.1)1 with δ v h = v h , 00 , we get

v h , 00 V = A v h , 00 L 2 ( K ) A v L 2 ( K ) .

By the Banach Closed Range Theorem,

α h = inf u h P p - 1 ( K ) sup v h V h , 0 Ω u h A v h u h K v h V = inf v h , 0 V h , 0 sup u h P p - 1 ( K ) Ω u h A v h u h K v h , 0 V .

This, along with (C.1)2, yields the estimate v h , 0 V = A v h , 0 K α h - 1 A v K . Similarly,

β h = inf w h P c p + 1 ( K + ) sup v h V h K + w b n w h v h w h K + v h V = inf v h V h sup w h P c p + 1 ( K + ) K + w b n w h v h w h K + v h V ,

along with (C.1)3, implies v h V = β h - 1 v K + . In conclusion,

v h V 2 = v h , 00 V 2 + v h , 0 V 2 + v h V 2 ( 1 + α h - 2 ) A v 2 + β h - 2 v K + 2 .

We present numerical experiments with the inf-sup constant 𝛽 for a quadrilateral element. Following the results for a triangular element, we consider a composite enriched space defined by partitions shown in Figure 12. In the case of a square element, only two outflow edges are possible (the first two cases shown in Figure 12), but in the case of a general quadrilateral, the third case with three outflow edges may occur as well. Also, the element may be subdivided in just two triangles (case not shown). We have experimented with the two elements shown in Figure 13. Depending upon the rotation angle, the square element may be refined into a quad and a triangle, two triangles or may not be refined at all. The quadrilateral element may undergo any of the three refinements shown in Figure 12. The minimum values (over all rotation angles) of the inf-sup constant β h for the two element shapes, p = 2 and different values of element size ℎ are shown in Table 6. The constant stays beautifully close to one, and it converges to one as h 0 . The same behavior is observed for elements of higher order 𝑝.

### Figure 12

Partitions of a general quad element into triangular and quadrilateral subelements defining the enriched test space.

### Table 6

Minimal (over angles) value of inf-sup constant β h for different values of element size ℎ for advection vector b = ( 1 , 0 ) , reaction coefficient c = 1.0 and element shapes shown in Figure 13.

quad / ℎ 1.0 0.1 0.01 0.001 0.0001
1 0.9949266746 0.9999887821 0.9999999876 0.9999999999 0.9999999999
2 0.9878674760 0.9999641904 0.9999999585 0.9999999999 0.9999999999

### References

[1] P. Bringmann and C. Carstensen, An adaptive least-squares FEM for the Stokes equations with optimal convergence rates, Numer. Math. 135 (2017), no. 2, 459–492. 10.1007/s00211-016-0806-1Search in Google Scholar

[2] D. Broersen, W. Dahmen and R. P. Stevenson, On the stability of DPG formulations of transport equations, Math. Comp. 87 (2018), no. 311, 1051–1082. 10.1090/mcom/3242Search in Google Scholar

[3] J. Brunken, K. Smetana and K. Urban, (Parametrized) first order transport equations: Realization of optimally stable Petrov–Galerkin methods, SIAM J. Sci. Comput. 41 (2019), no. 1, A592–A621. 10.1137/18M1176269Search in Google Scholar

[4] T. Bui-Thanh, L. Demkowicz and O. Ghattas, A unified discontinuous Petrov–Galerkin method and its analysis for Friedrichs’ systems, SIAM J. Numer. Anal. 51 (2013), no. 4, 1933–1958. 10.1137/110854369Search in Google Scholar

[5] C. Carstensen, L. Demkowicz and J. Gopalakrishnan, Breaking spaces and forms for the DPG method and applications including Maxwell equations, Comput. Math. Appl. 72 (2016), no. 3, 494–522. 10.1016/j.camwa.2016.05.004Search in Google Scholar

[6] C. Carstensen and F. Hellwig, Low-order discontinuous Petrov–Galerkin finite element methods for linear elasticity, SIAM J. Numer. Anal. 54 (2016), no. 6, 3388–3410. 10.1137/15M1032582Search in Google Scholar

[7] W. Dahmen, G. Kutyniok, W.-Q. Lim, C. Schwab and G. Welper, Adaptive anisotropic Petrov–Galerkin methods for first order transport equations, J. Comput. Appl. Math. 340 (2018), 191–220. 10.1016/j.cam.2018.02.023Search in Google Scholar

[8] W. Dahmen and R. P. Stevenson, Adaptive strategies for transport equations, Comput. Methods Appl. Math. 19 (2019), no. 3, 431–464. 10.1515/cmam-2018-0230Search in Google Scholar

[9] L. Demkowicz, Computing with h p Finite Elements. I. One- and Two-Dimensional Elliptic and Maxwell Problems, Chapman & Hall/CRC, Boca Raton, 2006. Search in Google Scholar

[10] H. De Sterck, T. A. Manteuffel, S. F. McCormick and L. Olson, Least-squares finite element methods and algebraic multigrid solvers for linear hyperbolic PDEs, SIAM J. Sci. Comput. 26 (2004), no. 1, 31–54. 10.1137/S106482750240858XSearch in Google Scholar

[11] L. Demkowicz and J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. Part I: The transport equation, Comput. Methods Appl. Mech. Engrg. 199 (2010), no. 23–24, 1558–1572. 10.1016/j.cma.2010.01.003Search in Google Scholar

[12] L. Demkowicz and J. Gopalakrishnan, A class of discontinuous Petrov–Galerkin methods. II. Optimal test functions, Numer. Methods Partial Differential Equations 27 (2011), no. 1, 70–105. 10.1002/num.20640Search in Google Scholar

[13] L. Demkowicz, J. Gopalakrishnan, S. Nagaraj and P. Sepúlveda, A spacetime DPG method for the Schrödinger equation, SIAM J. Numer. Anal. 55 (2017), no. 4, 1740–1759. 10.1137/16M1099765Search in Google Scholar

[14] L. Demkowicz and N. Heuer, Robust DPG method for convection-dominated diffusion problems, SIAM J. Numer. Anal. 51 (2013), no. 5, 2514–2537. 10.1137/120862065Search in Google Scholar

[15] L. Demkowicz and P. Zanotti, Construction of DPG Fortin operators revisited, Comput. Math. Appl. 80 (2020), no. 11, 2261–2271. 10.1016/j.camwa.2020.07.020Search in Google Scholar

[16] I. Ekeland and R. Temam, Convex Analysis and Variational Problems, Stud. Math. Appl. 1, North-Holland, Amsterdam, 1976. Search in Google Scholar

[17] A. Ern and J.-L. Guermond, Discontinuous Galerkin methods for Friedrichs’ systems. I. General theory, SIAM J. Numer. Anal. 44 (2006), no. 2, 753–778. 10.1137/050624133Search in Google Scholar

[18] A. Ern and J.-L. Guermond, A converse to Fortin’s lemma in Banach spaces, C. R. Math. Acad. Sci. Paris 354 (2016), no. 11, 1092–1095. 10.1016/j.crma.2016.09.013Search in Google Scholar

[19] J. Gopalakrishnan, P. Monk and P. Sepúlveda, A tent pitching scheme motivated by Friedrichs theory, Comput. Math. Appl. 70 (2015), no. 5, 1114–1135. 10.1016/j.camwa.2015.07.001Search in Google Scholar

[20] J. Gopalakrishnan and W. Qiu, An analysis of the practical DPG method, Math. Comp. 83 (2014), no. 286, 537–552. 10.1090/S0025-5718-2013-02721-4Search in Google Scholar

[21] M. Jensen, Discontinuous Galerkin methods for Friedrichs systems with irregular solutions, PhD thesis, Corpus Christi College, University of Oxford, 2004. Search in Google Scholar

[22] P. Joly, Some trace theorems in anisotropic Sobolev spaces, SIAM J. Math. Anal. 23 (1992), no. 3, 799–819. 10.1137/0523042Search in Google Scholar

[23] J. Muñoz Matute, D. Pardo and L. Demkowicz, A DPG-based time-marching scheme for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 373 (2021), Paper No. 113539. 10.1016/j.cma.2020.113539Search in Google Scholar

[24] J. Muñoz Matute, D. Pardo and L. Demkowicz, Equivalence between the DPG method and the exponential integrators for linear parabolic problems, J. Comput. Phys. 429 (2021), Paper No. 110016. 10.1016/j.jcp.2020.110016Search in Google Scholar

[25] S. Nagaraj, S. Petrides and L. F. Demkowicz, Construction of DPG Fortin operators for second order problems, Comput. Math. Appl. 74 (2017), no. 8, 1964–1980. 10.1016/j.camwa.2017.05.030Search in Google Scholar

[26] J. T. Oden and L. F. Demkowicz, Applied Functional Analysis, 3rd ed., CRC Press, Boca Raton, 2018. Search in Google Scholar