Discontinuous Galerkin Methods for a Class of Nonvariational Problems

We extend the finite element method introduced by Lakkis and Pryer (SIAM J. Sci. Comput. 33(2): 786–801, 2011) to approximate the solution of second-order elliptic problems in nonvariational form to incorporate the discontinuous Galerkin (DG) framework. This is done by viewing the “finite element Hessian” as an auxiliary variable in the formulation. Representing the finite element Hessian in a discontinuous setting yields a linear system of the same size and having the same sparsity pattern of the compact DG methods for variational elliptic problems. Furthermore, the system matrix is very easy to assemble; thus, this approach greatly reduces the computational complexity of the discretisation compared to the continuous approach. We conduct a stability and consistency analysis making use of the unified framework set out in Arnold et al. (SIAM J. Numer. Anal. 39(5): 1749–1779, 2001/2002). We also give an a posteriori analysis of the method in the case where the problem has a strong solution. The analysis applies to any consistent representation of the finite element Hessian, and thus is applicable to the previous works making use of continuous Galerkin approximations. Numerical evidence is presented showing that the method works well also in a more general setting.


Introduction
Linear, second-order, nonvariational partial differential equations (PDEs) are those which are given in the form where X∶Y = trace X ⊺ Y is the Frobenious inner product between matrices. If the matrix A is differentiable, then there is an equivalence between this problem and its variational sibling where and d is the dimension under consideration. Rewriting in this form is sometimes undesirable. For example, if the coefficient matrix A has near singular derivatives, the problem will become advection dominated and possibly unstable for conforming finite element methods. There is a wealth of material on the treatment of advection-dominated problems [c.f., 18,19]. If A is not differentiable, then the problem has no variational structure. In this case, standard finite element methods cannot be applied.
In a previous work [32], a finite element method for the approximation of the nonvariational problem (1) was introduced. This involved the introduction of a finite element Hessian represented in the same finite element space as the solution (modulo boundary conditions). The applications of the discrete representation of a Hessian of a piecewise function are becoming broader, for example, it can be used to drive anisotropic adaptive algorithms [1,40], as a notion of discrete convexity [2] and in the design of finite element methods for nonlinear fourth-order problems [37].
The algebraic formulation of the C 0 Galerkin approximation of the nonvariational problem requires the solution of large sparse (d + 1) 2 N 2 linear system [32,Lem 3.3], where N is the number of degrees of freedom. Equivalently, using a Schur complement argument, this can be reduced to an N 2 full linear system. The reason that this system is full is due to the global nature of the L 2 ( ) projection operator into a continuous finite element space. The motivation for extending the nonvariational finite element method into the discontinuous setting is the massive gain in computational efficiency over the continuous case. Indeed, due to the local representation of the projection operators in these discontinuous spaces, we are able to make massive computational savings, in that the system matrix becomes sparse and is the same size as that of a standard discontinuous Galerkin stiffness matrix for the Laplacian.
There has been a plethera of other finite element methods posed for linear nonvariational PDEs exist including [38] in which the authors pose a stable discontinuous Galerkin scheme for linear nonvariational problems using discrete analogs of the classical methods used to prove existence of strong solutions, see also [29] where similar techniques are used over curved domains. In [22][23][24] the author considers a least squares formulation, see also [31,35] for related formulations. Many of these formulations share the same disadvantage in the conditioning of the system matrix-the condition number scales like a fourth-order problem. As already mentioned, one of the desirable features of the formulation presented here is that it scales in the same fashion as a discontinuous Galerkin method for a second-order variational problem, a property shared by the method in [25] under the Cordes condition.

3
We are particularly interested in nonvariational PDEs due to their relation to general fully nonlinear PDEs which are of significant current research. In the literature finite element methods have been presented to solve this general class of problem. For example in [8] the author presents a C 1 finite element method shows stability and consistency (hence convergence) of the scheme which requires a high degree of smoothness on the exact solution. In [20,21] the authors give a method in which they approximate the general second-order fully nonlinear PDE by a sequence of fourth-order quasilinear PDEs. This is reminiscent of the vanishing viscosity method introduced for classically studying first-order fully nonlinear PDEs. Efficiency of any method used to approximate a problem such as this is key. Each of the methods is computationally costly due to their reliance on C 1 finite elements [8,21] or mixed methods [20].
In [5], a generic framework was set up to prove convergence of numerical approximations to the viscosity solutions of degenerate elliptic fully nonlinear PDEs. This involved constructing monotone sequences of approximations which are typically applied to finite difference approximations of the nonlinear problem [c.f., 36]. The assumption of consistency made in the [5] framework is incompatible with finite element methods; however, an extremely important observation made in [28] is that the consistency condition may be weakened to incorporate the finite element case using a localisation argument (in the case of isotropic diffusion).
A posteriori analysis of linear nonvariational problems is less standard than for those of variational type. Typically, results are based on "closeness" conditions of Cordes type [12] which guarantees existence of a unique smooth, H 2 solution. Under these assumptions, it is relatively straightforward to derive a posteriori bounds in the natural norm for the problem, H 2 , or a mesh-dependent equivalent. For other methods, this has been done in [24]. It has even been shown that adaptive methods for these problems converge optimally in terms of the number of degrees of freedom [30].
We show similar a posteriori bounds that make it straightforward to incorporate the method into well-developed software package for finite element methods, shown here using Dune [6,7]. In this work, we are interested in the asymptotic behaviour of the discontinuous approximation and the computational gains using this method. In a subsequent work, we will study the computational gains using the discontinuous framework presented over the continuous one given in [32], as well as exploit the powerful parallelisation capabilities of the package.
The rest of the paper is set out as follows. In Sect. 2, we formally introduce the model problem and give a brief review of known classical facts about nonvariational PDEs. In Sect. 3, we examine the discretisation of the nonvariational method in the discontinuous Galerkin framework, making use of the unified framework set out in [4] to derive a very general formulation of the finite element Hessian represented as a discontinuous object. We present some examples and examine the natural question of what happens when we try to eliminate the finite element Hessian from the formulation. In Sect. 4, we conduct an a posteriori analysis of the method and show upper and lower bounds to the error in an H 2 -like mesh-dependent norm. Finally, in Sect. 5, we detail a summary of extensive numerical experiments aimed at examining convergence and robustness of the method presented.

Problem Formulation
In this section, we formulate the model problem, fix notation and give some basic assumptions. In addition, we review the existence and uniqueness of the nonvariational problems. Let ⊂ ℝ d , d = 2, 3 , be a connected domain with polygonal boundary. The Lebesgue spaces are defined as and the Sobolev and Hilbert spaces These are equipped with the norms and derivatives D are understood in a weak sense. We pay particular attention to the cases k = 1, 2 and The model problem in strong form is: find u ∈ H 2 ( ) ∩ H 1 0 ( ) such that where the data f ∈ L 2 ( ) are prescribed and L is a general linear, second-order, uniformly elliptic partial differential operator. Let A ∈ L ∞ ( ) d×d , we then define We assume that A is uniformly positive definite, i.e., there exists a > 0 such that for all x and we call the ellipticity constant.

Definition 1 (Strong solution)
A strong solution of (1) is a function u ∈ H 2 ( ) ∩ H 1 0 ( ) , that is a twice weakly differentiable function, which satisfies the problem almost everywhere.
Theorem 1 (Existence and regularity of a strong solution of (1) [12]) Let ⊂ ℝ d be a convex polytope. Suppose A ∈ L ∞ ( ) d×d is uniformly elliptic and f ∈ L 2 ( ) . Suppose further that A satisfies the Cordes condition, that there exists 0 < ∈ L ∞ ( ) and 0 < , ∈ ℝ with + < 1 such that Then, (1) admits a unique strong solution. There also exists a constant independent of u such that Remark 1 (Less regular solutions) Note that the theory of viscosity solutions has been developed for non-classical solutions of (10), if the problem data do not satisfy the regularity assumed above, see [27].
Assumption 1 (Inf-sup condition) From hereon in, we will assume that the problem satisfies an inf-sup condition, that is, with then, for all w ∈ H 2 ( ) ∩ H 1 0 ( ) This is true under a variety of conditions. For example, those of Theorem 1. Note that for d = 2 , this criterion is satisfied for all essentially bounded, symmetric positive definite A.

Discretisation
Let T be a shape regular triangulation of , namely, T is a finite family of sets such that , 2) for any K, J ∈ T we have that K ∩ J is a full subsimplex (i.e., it is either ∅ , a vertex, an edge, a face, or the whole of K and J ) of either K and J and 3) We use the convention where h ∶ → ℝ denotes the mesh size function of T , i.e., where h K is the diameter of K. We let E be the skeleton (set of common interfaces) of the triangulation T and say e ∈ E if e is on the interior of and e ∈ if e lies on the boundary .
Let ℙ k (T) denote the space of piecewise polynomials of degree k over the triangulation T , i.e., and introduce the finite element spaces to be the usual spaces of discontinuous piecewise polynomial functions.
Remark 2 (Generalised Hessian) Assume v ∈ H 2 ( ) , let n ∶ → ℝ d be the outward pointing normal of , then the Hessian D 2 v of v, satisfies the following identity: If v ∈ H 1 ( ) , the right-hand side of (21) is still well defined in view of duality, in this case we set where the last term is understood as a duality pairing.

Definition 2 (Broken Sobolev spaces, trace spaces) We introduce the broken Sobolev space
We also make use of functions defined in these broken spaces restricted to the skeleton of the triangulation. This requires an appropriate trace space Definition 3 (Jumps, averages and tensor jumps) We define average, jump and tensor jump operators for arbitrary scalar functions v ∈ T(E) , vectors v ∈ T(E) d and matrices V ∈ T(E) d×d as Note that on the boundary of the domain , the jump and average operators are defined as We will often use the following Proposition which we state in full for clarity but whose proof is merely using the identities in Definition 3.

Construction of an Appropriate Discrete Hessian
We now use the framework set out in [4] to construct a general notion of discrete Hessian. We first give a definition using a flux formulation.
for all ∈ .
We now present the primal formulation for the generalised finite element Hessian. Proof Note that in view of Definition 3 for generic vector fields q ∈ and v ∈ , we have the following identity: Then summing (37) over K ∈ T and making use of the identity (40) we see Using the same argument for (38) Note that, again making use of (40), we have for each q ∈ H 1 (T) d and v ∈ H 1 (T) that Taking v = u in (43) and substituting into (38), we see Now choosing q = ∇ h and substituting (44) into (37), we arrive at the fully generalised finite element Hessian given by (39).

Remark 3 (Consistent representations of the gradient operator)
If one were interested in consistent representations of other derivatives, for example the gradient operator, one would need to modify the proof of Theorem 2. Examples of consistent gradient representations can be found in [4]. See also [10,11,17]. Using this methodology, it should be possible to construct an entire hierarchy of derivatives.
Example 1 An example of a DG formulation for the approximation to the Hessian, D 2 u , can be derived by taking the fluxes in the following way: The result is a discrete representation of the Hessian H[u h ] as a unique element of d×d such that We can also derive unsymmetric approximations with fluxes similar to those used for the NIPG method. This is given by taking = −1 in the following, while the symmetric version is given by = 1:

The Discontinuous Nonvariational Finite Element Method
We are now in a position to state the numerical method for the approximation of (1). We look to find u h ∈ 0 together with H[u h ] ∈ d×d such that with where the penalisation parameter > 0 is to be chosen sufficiently large. Using the L 2 projection operator P k ∶ L 2 ( ) → defined for v ∈ L 2 ( ) through it is possible to elliminate the finite element Hessian from the bilinear form for sufficiently smooth A. Proof This follows from the following identity:

Remark 4
The solution of the problem in this form is nontrivial due to the global L 2 ( ) projection appearing in the formulation. However, in the discontinuous setting, the global L 2 ( ) projection is in fact computable locally. We may actually exploit this fact to optimise our schemes efficiency. We will discuss this further in the sequel. (1) we have that A = I , then we have that and our bilinear form reduces to since P k ( A) = I. The nonvariational finite element method, thus, coincides with the classical (symmetric) interior penalty method for the Laplacian [16].

Remark 5 (Relation to standard DG methods)
It is not difficult to prove that choosing to numerical fluxes in the same way as presented in [4, Table 3.2] results in the same correlation to the DG methods summarised in the aforementioned paper for the case that A is constant. For brevity, we will not prove this here. Note that when A is not constant, we have that the nonvariational finite element method does not coincide with its standard variational finite element counterpart. Indeed, the method is able to successfully cope with classes of advection dominated problems not falling under the variational framework without special treatment [32, §4.2] which is illustrated by the result of Lemma 1.
We conclude this section with a proof of consistency of the method and then show that Galerkin orthogonality holds. with the error functional given by

A Posteriori Analysis
The discrete bilinear form (49) only makes sense over H 2 (T) × H 2 (T) . To allow for an a posteriori bound we require an extension to ensure the appropriate stability arguments can be applied. We require the bilinear form to be extended to H 2 (T) × L 2 ( ) . To do this for (u, v) ∈ H 2 (T) × L 2 ( ) , we define (64) .
Proof Beginning from Proposition 3 we note that the bound can be split into a residual component I 1 , a nonconformity component I 2 and an inconsistency component I 3 as follows: Now, we proceed to bound these term by term. In view of Cauchy-Schwartz, we have For the nonconformity term, we note that by the properties of the reconstruction given in (69). For the inconsistency term, we have that Making use of Lemma 4 and the properties of the reconstruction (69), we see that (77) (78) Again from Lemma 4 we have that and finally (69) gives that Hence, we have that Collecting (77)

Numerical Experiments
In this section, we detail numerical experiments carried out in the finite element package Dune-Fem [15] which is based on the Dune software framework [6,7]. The code makes use of the newly developed Python frontend [15] and the unified form language [3] is used to provide the problem data. The code will be made freely available within the Dune-Fem-tutorial in a future release [13]. We present some benchmark problems designed such that the exact solution is known. In each of the experiments the domain = [0, 1] 2 and we consider the coefficient matrix to be varying a(x) and b(x) . We study three different choices using x = x 1 , x 2 as follows. . (83) . (86) In this test, we take the components of A such that the differential operator can be written in variational form and is coercive, fitting into a standard analytical framework: 2) (Continuous not H 2 ) In this test, we take A such that it is comparable to [32, §4.4] 3) (Discontinuous not W 1,∞ ) In our third test, a is discontinuous and not in W 1,∞ . We take b ≡ 0 and choose with Note that the initial grid is chosen so that it aligns with the discontinuity. We study the behaviour of our method for polynomial degrees k = 1, 2, 3 , choosing the forcing such that the exact solution is given by the following two choices.
The penalty parameter is chosen as = A,max 5k(k + 1) where A,max is the maximum eigenvalue of A. Benchmarking for these tests is shown in Figs. 1, 2 and 3. For linear polynomials, the H 2 (T)-error does not converge since the piecewise Hessian of u h vanishes. So we only show the errors in the L 2 ( ) and H 1 (T) norms. For the smooth solution (left column), we clearly see that the method converges optimally in both norms. While for the less smooth solution convergence is slowed as would be expected when approximating an H 2 ( ) solution. For polynomial orders 2 and 3 convergence is optimal for all norms studied here when approximating a smooth solution and between 1 and 1.5 for the H 2 ( ) solution. The convergence rate hardly depends on the smoothness of A. The only case where some clear dependency is visible is in the L 2 ( ) errors for quadratic polynomials with the discontinuous A . In this case, the order of the L 2 ( ) error seems to equal the convergence in the H 1 (T) norm, i.e., is not optimal (see Fig. 2).
We also study the condition number of the system matrix generated when assembling (49). The results are shown in Fig. 4. In contrast to many other methods, which rewrite the Fig. 1 Convergence rates for the error measured in H 1 ( ) and L 2 ( ) for k = 1 testing the case when u is either a prescribed smooth solution or u ∈ H 2 ( )∕H 3 ( ) . We also test three different choices for diffusion coefficient A that are coercive, continuous but not H 2 and discontinuous. The rates are optimal in H 1 ( ) in all cases and in L 2 ( ) in most cases. Note we did not present the H 2 ( ) rates nor the estimate since neither will converge for k = 1 nonvariational model as a fourth-order problem, the condition number depends on the grid spacing in the same way as it does for second-order variational problem, i.e., it is O(h −2 ).
Finally, we study the behaviour of the residual error indicator and an adaptive scheme for the approximation of the solutions to these problems. The indicator is included in Figs. 2, 3 where its reliability is clearly visible. A comparison of these globally refined simulations with an adaptive simulation using an equal distribution strategy for locally refining the grid are given in Figs. 5 and 6. As to be expected, no advantage can be gained when approximating a smooth solution.  2 Convergence rates for the error measured in H 2 ( ) , H 1 ( ) and L 2 ( ) for k = 2 and the estimator given in Theorem 3. We test the case when u is either a prescribed smooth solution or u ∈ H 2 ( )∕H 3 ( ) . We also test three different choices for diffusion coefficient A that are coercive, continuous but not H 2 and discontinuous. The rates are optimal in H 2 ( ) and H 1 ( ) in all cases and in L 2 ( ) in most cases. When the solution is not smooth the rates are slowed to an expected rate in-line with the regularity of u. In all cases, the estimate is efficient and robust For the H 2 ( ) solution, the adaptive simulation approximately doubles the convergence rate of the scheme.
We conclude with a simulation for which we do not have an exact solution. We use the discontinuous A and choose a constant forcing f ≡ 1 000 . As before boundary conditions are equal to zero. Results are shown in Fig. 6 where we show the convergence of the residual indicator under global and local refinement. A visualisation of the function a, the resulting discrete solution, and the adaptive grid are given in Fig. 7.   Fig. 3 Convergence rates for the error measured in H 2 ( ) , H 1 ( ) and L 2 ( ) for k = 3 and the estimator given in Theorem 3. We test the case when u is either a prescribed smooth solution or u ∈ H 2 ( )∕H 3 ( ) . We also test three different choices for diffusion coefficient A that are coercive, continuous but not H 2 and discontinuous. The rates are optimal in H 2 ( ) and H 1 ( ) in all cases and in L 2 ( ) in most cases. When the solution is not smooth the rates are slowed to an expected rate in-line with the regularity of u. In all cases, the estimate is efficient and robust. Fig. 4 Condition number estimates of the system matrix given by (49). Notice that the condition number grows in the same asymptotic fashion as the IP-DG method for the Laplacian ( A =Identity). Also, the complexity does not change asymptotically as k is increased Fig. 5 Comparison of simulations carried out on adaptive and globally refined grids driven by the estimate from Theorem 3. We test the case when u is either a prescribed smooth solution or u ∈ H 2 ( )∕H 3 ( ) and when A is discontinuous. For the smooth solution, we find the uniform and adaptive schemes behave similarly, for the solution that is not H 3 ( ) , the adaptive scheme clearly outperforms the uniform one, although there is no gain from using k = 3 Fig. 6 Convergence of estimator under adaptive and globally refined grids. Here, we consider the case when A is discontinuous, shown in Fig. 7a, and the forcing is chosen f = 100 . No exact solution is known for this case. There is a clear indication that the adaptive scheme outperforms the uniform one Fig. 7 Visualisation of the problem coefficient, the solution and the underlying adaptive mesh refinement level. Notice the algorithm refines where the problem data are discontinuous and well as near boundary layers caused by the anisotropy of the diffusion tensor

Conclusions and Outlook
In this work, we have extended the framework from [32] for linear nonvariational problems to incorporate discontinuous approximations. We have derived a posteriori bounds for this problem and shown they are useful to drive adaptive algorithms. We would like to point out that the approach presented here can be directly applied to the case where the approximation u h is chosen in a continuous finite element space but the finite element Hessian is defined in the discontinuous fashion described here.
In the numerical experiments, we note the method is well posed and converges optimally even for A that are discontinuous. The method is well suited to solve nonlinear problems [33] and this will be the topic of ongoing research.

Compliance with Ethical Standards
Conflict of Interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.