Hybridizable discontinuous Galerkin method with mixed-order spaces for non-linear diffusion equations with internal jumps

We formulate a hybridizable discontinuous Galerkin method for parabolic equations with non-linear tensor-valued coefficients and jump conditions (Henry’s law). The analysis of the proposed scheme indicates the optimal convergence order for mildly non-linear problems. The same order is also obtained in our numerical studies for simplified settings. A series of numerical experiments investigate the effect of choosing different order approximation spaces for various unknowns.


Introduction
Mathematical models involving jump conditions at phase or material interfaces are widespread in many CFD, environmental, and biomedical applications.Examples include coupled subsurface/free-surfaces flow systems (Sochala et al. 2009), incom-pressible two-phase flows with different densities (Abels et al. 2012,) interfaces between cells or cell parts (Jäger et al. 2009), simulations of different phases within concrete (Muntean and Böhm 2009), and many more.For Darcy flows considered in the current work, the position of the interface is often fixed and known in advance, and the diffusion coefficient is a symmetric tensor with a uniformly symmetric positive definite inverse.This is precisely the setup investigated and analyzed in this paper for a numerical scheme based on a hybridizable discontinuous Galerkin (HDG) method.
The HDG methods first proposed by Cockburn and co-workers in Cockburn et al. (2008bCockburn et al. ( , 2009) ) inherit many benefits of discontinuous Galerkin (DG) methods such as natural support for non-conforming meshes and hanging nodes, hp-adaptivity, local mass conservation, etc.Moreover, due to the hybrid nature of these methods, they intrinsically support static condensation, which allows decreasing the size of the system of linear equations while keeping this linear system sparse and symmetric positive definite.Static condensation can be interpreted both in an algebraic sense or as an alternative interpretation of the numerical scheme in terms of skeleton unknowns and numerical fluxes.The latter interpretation suggests more general numerical fluxes (that can otherwise only be achieved at higher computational costs) and thus can be utilized to build HDG schemes with an increased accuracy, i.e. optimally convergent and even superconvergent methods, which can be exploited by post-processing techniques (see Nguyen et al. 2009 for an overview of the early developments of the HDG methods).
The aforementioned advantages of the HDG methodology gave rise to a large number of applications for different problems.To name just a few: incompressible Stokes (Nguyen et al. 2010) and Navier-Stokes (Nguyen et al. 2011a) equations, Maxwell equations (Nguyen et al. 2011b), linear elasticity (Di Pietro and Ern 2015), linearized shallow-water equations (Bui-Thanh 2016).Furthermore, comparisons of accuracy, robustness, and computational efficiency between HDG and DG (Woopen et al. 2014;Jaust et al. 2018), HDG and classical Finite Elements (Kirby et al. 2012), HDG and Finite Volumes (Ahnert and Bärwolff 2014) also received some attention.Notably, some recent studies (Vila-Pérez et al. 2022;Kronbichler and Wall 2018) suggest that HDG schemes are not necessarily faster than the corresponding DG discretizationsthis aspect must be rather evaluated on a case-by-case basis.However, in a decade of development, only a few studies focusing on porous media problems appeared in the HDG literature: For example, Samii et al. (2016) investigates an adaptive HDG method to deal with anisotropic diffusion in an elliptic setting, Costa-Solè et al. (2019), Fabien et al. (2018) apply hybrid methods to two-phase flows in porous media, Gatica and Sequeira (2018) analyzes the Brinkman problem, and Moon et al. (2019) investigates multiscale approaches that are related to HDG.However, discontinuities in the solution (arising from physics) are disregarded in their respective models.
Our main interest in connection with HDG schemes concerns the flow and transport in porous media, where their DG counterparts proved to be very useful with their affinity to unstructured meshes and adaptivity combined with low regularity requirements on the solution and coefficients-making the DG method well-suited for complex geometries encountered in real-world problems (Ray et al. 2017).Another feature distinguishing the HDG schemes from their DG counterparts is the convergence for approximation spaces of order zero.This capability is often desirable for porous media applications, where many numerical models utilize lowest order spaces (cf.Samii et al. 2016).
In the current work, we investigate an LDG-H method-a particular type of HDG scheme (Cockburn et al. 2009) related to the local DG (LDG) discretizations.We analyze stability and convergence when applied to instationary equations with (mildly) non-linear tensor-valued diffusion coefficients exemplified by Darcy flows.Numerically, we confirm our analytical findings for linear, optionally anisotropic diffusion tensors on polytopically bounded domains with flat interfaces.Of particular interest are approximation spaces of different orders for different unknowns such as in Cockburn et al. (2009), Rupp et al. (2018).Our formulation in principal allows for non-linear diffusion tensors and enforces jumps in the primary unknown described by Henry's law.Similar problems have already been analyzed for LDG (Rupp et al. 2018) and enriched Galerkin (Rupp and Lee 2020) schemes.The analysis relies in part on techniques employed for an LDG formulation in Aizinger et al. (2018), Rupp and Knabner (2017), Rupp et al. (2018), Reuter et al. (2019) and extends them to the LDG-H method.
The first analysis for a HDG-type method was given by Cockburn et al. (2008b) for a special LDG-H variant called the single face hybridizable (SFH) method; its distinctive feature is the penalty coefficient τ chosen as a positive constant on one face of each simplex and zero everywhere else.A more general result for elliptic problems followed in Cockburn et al. (2010) that relied on a special projection (exclusively for simplex-meshes) tailored to fit the characteristics of HDG methods.These results were generalized to non-conforming simplicial meshes in Chen and Cockburn (2012) and squares, cubes, and prisms in Cockburn et al. (2012).An analysis for a Poisson equation involving additive jumps in the primary and the flux unknowns (as opposed to multiplicative jumps considered here) was presented in Dong et al. (2017), and parabolic problems with non-linear diffusion tensors were considered in Moon et al. (2017).
The current paper (based, in part, on results presented in Rupp 2019) is structured as follows.The next section introduces the model problem followed by its semi-discrete LDG-H formulation in Sect.3. Section 4 contains the energy stability and convergence estimates for the semi-discrete formulation.Numerical convergence studies and a realistic problem with discontinuous diffusion coefficients are presented in Sect. 5.

Problem formulation
Similarly to Rupp et al. (2018), we consider a bounded Lipschitz domain Ω ⊂ R d subdivided into two open, disjoint, non-degenerated, Lipschitz polytopes Ω l , Ω g such that Ω = Ωl ∪ Ωg .We assume ∂Ω to be disjointly subdivided into Γ D and Γ N denoting the Dirichlet and Neumann boundaries, respectively.Moreover, Γ LG := ∂Ω l ∩ ∂Ω g is the boundary between Ω l and Ω g where the Henry jump condition is imposed.The problem can be formulated as follows: ) for ,d is assumed to be strongly Lipschitz and uniformly symmetric positive definite with constants L D and C D , respectively (see Remark 2.1).Here and in the following, ν α denotes the outward unit normal with respect to Ω α (α = l, g), and ν is the outward unit normal with respect to Ω.
Remark 2.1 1.In this work, D −1 (•) being Lipschitz means that there is a constant L D > 0 such that for all v, w ∈ L 2 (Ω) (2) In the setting of the Darcy equation, the flux unknown q = −D (u) ∇u can be interpreted as Darcy flux, and its "continuity" (1d) across Γ LG ensures the conservation of mass, while the primary unknown u can be interpreted as the hydraulic head.
Having defined our problem, we now specify the requirements for a weak solution to the above problem.
Problem 2.2 (Weak problem) A weak solution of (1) on domain Ω is defined as 2. The restriction of q to Ω α fulfills q α ∈ L 2 0, T ; H 1 (Ω α ) d for α = l, g. 3. System (3) holds for any smooth test functions ϕ u , ϕ q , any Lipschitz subdomain Ω e ⊂ Ω, and for a.e.t ∈ (0, T ): ) Note that this problem can be formulated under quite general conditions concerning the non-linearity.Our analysis covers the case of mildly non-linear problems.At the same time, our numerical studies in Sect. 5 are limited to linear problems and explore the effects of different approximation polynomial orders, discontinuity size, and diffusion coefficient anisotropy.

Basic definitions and some auxiliary results
The notation will be, as far as possible, kept analogous to the notation in Rupp et al. (2018) with some parts taken from Samii et al. (2016).This means that in the following, T h = {K i : i = 1, . . ., N el } denotes a d-dimensional non-overlapping polytopic partition of Ω (see Di Pietro and Ern 2012, Def.1.12) such that, for each element K i , either K i ⊂ Ω l or K i ⊂ Ω g is fulfilled.For simplicity, all partitions of Ω are assumed to be geometrically conformal (in the sense of Ern and Guermond 2004, Def 1.55).We denote by F = F (T h ) the set of faces, by F I the set of interior faces (that do not intersect with Γ LG ), by F D Dirichlet, and by F N Neumann boundary faces, and by F LG the set of faces that belong to Γ LG .Note that every face of the mesh is an element of one and only one of the following sets: The test-and ansatz-spaces are defined as the broken polynomial spaces of order at most k: For hybridization purposes, another ansatz space-the so-called skeleton space on the element faces (excluding those with the Dirichlet boundary conditions)-needs to be defined: For an element-wise defined scalar function w and an element-wise vector function v, we define the jump • on ∂K i ∩ ∂K j for neighboring mesh elements K i , K j ∈ T h , K i = K j in the following way: where ν K is the outward unit normal with respect to K. Note that a jump in a scalar variable is a vector, whereas a jump in a vector is a scalar.For F ∈ F I , an interior face shared by cells K − and K + , and for x ∈ F, we define the one-sided values of a scalar quantity w = w(x) by and for F ∈ F LG , a face on Γ LG shared by cells K l ⊂ Ω l and K g ⊂ Ω g , and for x ∈ F, we define the one-sided values of a scalar quantity w = w(x) by Definition 3.1 (Shape and contact regularity) A family of meshes T h is said to be shape and contact regular (short regular) if, for all h > 0, T h admits a geometrically conformal, matching simplicial submesh Th (Di Pietro and Ern 2012, Definition 1.37) such that 1. Th h>0 is shape-regular in the usual sense of Ciarlet (1990), meaning that there is a parameter λ 1 > 0, independent of h, such that for all K ∈ Th , where ρ K is the diameter of the largest ball that can be inscribed in K. 2. there is a constant λ 2 > 0, independent of h, such that for all K ∈ T h and for all K ∈ Th with K ⊂ K, Lemma 3.2 (Discrete trace inequality) Let (T h ) h>0 be a regular mesh family with parameters λ 1 , λ 2 , λ 3 .Then, for all h > 0 and all p ∈ P d k (T h ), the following holds with C tr only depending on λ 1 , λ 2 , d, and k: For an interior face, p L 2 (F) is assumed to contain traces from both elements Proof The first inequality follows from Lemma 1.46 in Di Pietro and Ern (2012) and the fact that h F ≤ h T for all faces F of element T .The second inequality follows from the first one and the fact that ν is a constant vector on each face.
To simplify notation, we denote in the following • L 2 (A) as • A and use explicit norm notation for all other norms.

Spatial discretization
For hybridization, a new variable on the skeleton space P k (F \ F D ) is introduced and denoted by λ h .In the following, we propose a modification of the LDG-H method described in Samii et al. (2016) capable of dealing with jump conditions across a submanifold Γ LG (i.e., Henry's Law) using instationary Darcy flow as an example application.An investigation of the same model problem using a classical LDG method was conducted in Rupp et al. (2018).To deal with Henry's law, the numerical trace ûh is defined as whereas the numerical flux qh is chosen as where τ ≥ 0 is a stabilization parameter (possibly depending on h).For more information on the stabilization parameter τ , the reader is referred to Cockburn et al. (2008a), Cockburn et al. (2008b).
The (semi-discrete) LDG-H formulation of the model problem (1) reads: Here, k, k, k denote orders of polynomial ansatz spaces for different unknowns. 123 Remark 3.4 (5c) is a globally coupled system that ensures the local mass conservation, whereas (5a), (5b) are element-local and can be solved elementwise.

Stability and error analysis
Our analysis heavily relies on techniques presented in Aizinger et al. (2018), Rupp and Knabner (2017), Rupp et al. (2018), Reuter et al. (2019), additional results used in the proofs are introduced in the corresponding sections.

Stability estimate
First, we introduce a new function H LG (x) that generalizes the Henry coefficient to the whole domain Ω (note that argument x is omitted later on): Theorem 4.1 (Energy stability of the semi-discrete problem) Assume that (u h , q h , λ h ) is a solution of Problem 3.3.For every regular, geometrically conformal family of meshes (T h ) h>0 and for all τ, h > 0 in (5), there exists a function C (s) such that for a.e.s ∈ (0, T ) Here, the second term is interpreted as We highlight that C(s) depends on h and becomes larger as h goes to 0. Nonetheless, for each fixed value of h, we obtain temporal stability and therefore globally unique existence of semi-discrete solutions.
Proof The proof uses ideas from the proof of Rupp et al. (2018, Theorem 4.9).Special care must be taken of the boundary with the specified jump condition Γ LG .We test Eq.( 5) with H LG u h , H LG q h , and − H LG λ h , integrate the second equation by parts, and sum the first two over all K ∈ T h .This yields Adding the above equations together and using algebraic simplifications gives Using Young's and Cauchy-Schwarz' inequalities, Ξ i terms can be estimated as follows: The last inequality is obtained by first invoking Lemma 3.2 and then using (2) pointwise.Substituting the above estimates into (6) and collecting similar terms gives us the following expression: The stability result follows from, first, using the discrete trace inequality given in Lemma 3.2 on boundary integral terms involving unknown solution and then applying Grönwall's inequality.Notably, this might create a factor h −1 F,min (with h F,min denoting the minimum face diameter) on the right-hand side of the stability estimate.
Remark 4.2 Although the factor h −1 F,min enters the right-hand side of the stability estimate, this estimate can still be used to ensure that for any fixed h > 0 the discrete solution is stable and therefore exists globally for any time.Our estimate should be understood in that sense.However, this result does not imply the existence of a weak solution by letting h 0. In order to show that the weak solution of (1) exists uniquely, one can, for example, use the results of Rupp and Lee (2020) and plug them into the general existence and uniqueness theory presented in Evans (2010, §7.1.2).

Convergence order estimate
The convergence analysis of the semi-discrete scheme given in Problem 3.3 relies on a special type of Projection (Π w , Π v ) introduced in Cockburn et al. (2010) and also used in Chabaud and Cockburn (2012).It is defined on each simplex K as The approximation properties of this projection are given in the following Then for k ≥ 0 and for given u ∈ H k+1 (K) and q ∈ H k+1 (K) d , system (7) is uniquely solvable for Π v q and Π w u.Furthermore, there exists a constant C Π independent of K and τ such that Proof This corollary trivially follows from Cockburn et al. (2010, Theorem 2.1) proved in the appendix of Cockburn et al. (2010).For k = 0, only the third equation of ( 7) makes sense and is sufficient to compute the projection.

Remark 4.4
The fact that this analysis technique only works for simplex-shaped elements certainly constitutes a limitation of our methodology; however, our convergence results shown in Sect. 5 indicate that our scheme also converges with the same order for rectangular elements.
For the analytical solution (u, q) and the discrete solution (u h , q h ) of our problem, we denote Here, Π denotes an L 2 -projection acting on a given ζ ∈ L 2 (F\F D ).Note that the restriction of Πζ to F ∈ F\F D is in P k (F) and satisfies: Note that H LG e g λ = e l λ on Γ LG since u l = H LG u g and Π is linear.To deal with non-linear diffusion coefficients we state Lemma 4.5 The following inequality holds for u ∈ H k+1 (Ω α ), , and e q defined in (8): for some arbitrary ε > 0.
Proof Using Hölder's, Young's inequalities and properties of D −1 (•) we obtain and the result then follows by Lemma 4.3.
Remark 4.6 A result related to Lemma 4.5 is mentioned in Rupp et al. (2018, Rem. 4.8) without proof.In our analysis, the projection defined in Lemma 4.3 is used instead of the standard L 2 -projection utilized in Rupp et al. (2018).
Next, we prove our main convergence result.
Theorem 4.7 (Convergence order estimate for the LDG-H method) Let (u, q) ∈ H 1 0, T ; LG e q dσ + F∈FI F H LG e q e λ dσ Most of the projection error terms vanish due to the orthogonality properties of the chosen projection specified in Lemma 4.3.Adding up the equations, simple algebraic manipulations give us Using Young's and Cauchy-Schwarz' inequalities as well as Lemma 4.5 we get The remainder of the proof boils down to some simple algebraic manipulations, omitting some non-negative left-hand-side terms, and using Grönwall's inequality.

Remark 4.8
The fact that the above estimates contain τ in both numerator and denominator indicates that choosing τ as a constant (per face F ∈ F \ F D ) is optimal in contrast to τ = 1/h, or τ = h, or any other h-dependent choice.

Numerical results
In Secs.5.1, 5.3, the discretization in time is performed via the implicit Euler scheme, and the implementation has been carried out in the C++ based Finite Element toolbox M++ (Wieners 2005).The numerical results in Sec 5.2 were obtained using the programming framework introduced in Rupp et al. (2022), Rupp and Kanschat (2021) and rely on the Crank-Nicolson time stepping method.

Numerical convergence study for piecewise constant diffusion coefficients
First, we verify the convergence estimates by using a smooth test problem with known solution.Our smooth test problem is Ω l = (0, 10) × (0, 5) , Ω g = (0, 10) × (5, 10) , We choose Ω as a square of side length ten and set Neumann boundary conditions at the top and bottom boundaries while letting the left and right boundaries be of Dirichlet type.The right-hand side function, initial-, and boundary conditions are chosen appropriately.The spatial domain is discretized by 2 i × 2 i squares (where i denotes the number of refinement steps).As the mesh skeleton space associated with the edges grows with each mesh refinement, the error on the skeleton space is scaled with respect to the measure of the skeleton space.sadov Choosing the penalty parameter τ = 1 and using equal order approximations for all unknowns, the LDG-H-discretization shows an optimal order of convergence of k + 1 in both variables for ansatz spaces of order k (see Tables 1, 2 and 3 uppermost block).Notably, this means an improved order of convergence in the flux variable compared to the LDG discretization of the same order.For τ = 1/h the order of convergence in the flux variable decreases by one.This is consistent with analytical results and also with numerical studies for elliptic problems given in Cockburn et al. (2008a).
In Rupp et al. (2018), different orders of approximation spaces for the scalar and the flux variable were investigated for an LDG method.There it was shown analytically that if the order was reduced for the flux variable by one to k − 1 while the order for the scalar variable was still k, both variables would converge at least with order k.In numerical experiments presented in Rupp et al. (2018), the order of convergence for the scalar variable appeared to be higher than the analytical result indicates whereas the estimate for the flux variable was sharp.Next, we consider similar scenarios for our LDG-H scheme with penalty parameters τ = 1 and τ = 1/h, respectively.Since there is yet another ansatz space to approximate the solution on element boundaries, there arises one additional possibility to vary the order of approximation.For τ = 1, the experimental order of convergence for all variables decreased as soon as we reduced the order of one of the approximation spaces.In the case of τ = 1/h though, the order of convergence decreased when we reduced the order of approximation for the primary u h or the hybrid variable λ h , and no convergence was detected if the order of one of the corresponding spaces was 0. However, if the order of the approximation space for the flux variable is reduced by one, the same effect as for the LDG method can be seen: The experimental order of convergence turns out to be the same as for the

Table 1
Errors (E) and estimated orders of (eoc) for constant approximations

Table 2
Errors (E) and estimated orders of convergence (eoc) for linear approximations

Table 3
Errors (E) and estimated orders of convergence (eoc) for quadratic approximations  2 and 3 third block).Thus one may be able to substantially reduce the total number of the degrees of freedom (the flux unknown is a d-dimensional vector!) without sacrificing the convergence order.Usually, this is not considered to be significant since classical HDG literature assumes the local solutions to be computationally cheap and the total performance to be dominated by the global solve.However, hybrid methods are often implemented in a matrix-free fashion (see e.g.Rupp et al. 2022).In this case, the base assumption that solving a local system of equations is cheap still holds, but since local systems have to be solved in every mesh cell per iteration of an iterative solver per time step the local solution time dominates the overall runtime of the code.The authors experienced this bottleneck when implementing their own code (Rupp and Kanschat 2021).In this specific situation, a significant speedup is obtained by our approach.

Numerical convergence study for piecewise continuous, spatially varying diffusion coefficients
The convergence study presented in Sect.5.1 focused on a simple example in two spatial dimensions.There we observed that one sees the convergence orders of k + 1 for all unknowns if one chooses all polynomial spaces order k and τ = 1.Next, we demonstrate that a variable diffusion coefficient or the number of spatial dimensions does not affect the convergence rates.To this end, we consider The numerical results shown in Table 4 indicate that the expected orders of convergence are obtained for the primary unknown (for conciseness, we omit the errors for the flux and skeleton unknowns).The first two columns (d = 1, d = 2) consider the aforementioned setting.The third column (H LG = 1000) describes the two-dimensional example in which the factor H LG has been increased to 1000.Finally, the last column ( D) corresponds to the anisotropic tensor-valued coefficient In all cases, we can see similar convergence rates.This showcases that Henry's constant does not influence the convergence rate but does influence the constant, since the solution's semi-norm scales with it.Moreover, the results validate our claim that our method is suitable for anisotropic diffusion.Remark 5.1 (Sharp convergence rates) None of the aforementioned numerical experiments indicate that the observed rates of convergence exceed those obtained in our error analysis.Thus, we deduce that our analytical findings are sharp with respect to the observed rate.

Nested rectangle problem
To demonstrate the robustness of our method we show results for a more involved non-smooth test problem given by  On a square of side length ten, we set homogeneous Neumann boundary conditions on all boundaries and set the right-hand-side function to zero.The initial condition for the nested rectangle problem is plotted in Fig. 1 (top left), the solutions at time T = 1 for piecewise linear and piecewise quadratic (equal order) approximation are shown in Fig. 1 (middle left) and Fig. 1 (bottom left), respectively, and the difference plot between these two is displayed in Fig. 1 (top right).In order to illustrate the effect of reducing by one the approximation order for the flux unknown q h , we also display the difference plots between equal and reduced (for q h ) order solutions for linear in Fig. 1 (middle right) and quadratic in Fig. 1 (bottom right) approximations.Especially for the piecewise quadratic approximation the difference between the solutions appears to be very small.

Conclusions
The current work introduced and investigated a hybridizable DG formulation for parabolic equations with non-linear tensor-valued coefficients and jump conditions described by Henry's law.Optimal convergence estimates were proved for mildly non-linear problems on simplex meshes, and numerical studies demonstrated the same convergence on rectangular meshes.Our numerical studies of mixed-order approximations indicate the same behavior as in the case of LDG discretizations, namely the possibility of reducing the approximation order for flux unknowns by one polynomial order without any loss of convergence or accuracy.Extending the analysis techniques to more general element shapes appears to be very desirable though challenging.
An important future step that needs to be undertaken in order to apply our method to real-world scenarios is the implementation and analysis of the presented approach for cases, in which the interface is curved and evolving in time.In such scenarios, it is particularly important to accurately track the interface and update the quadrature rules on the interface and in the elements.In these cases, hybrid versions of Cut-FEM or extended FEM might be useful.Another possible research direction concerns more involved non-itineraries and numerical treatment of highly anisotropic diffusion tensors.
the scalar D for the two-dimensional experiments.

Table 4
Errors (E) and estimated orders of convergence (eoc) Solution for the nested rectangle problem (left panel): initial condition (top), final states at T = 1 for P 1 × P d 1 × P 1 (middle) and P 2 × P d 2 × P 2 (bottom).Difference plots (right panel): u h |P 2 ×P d