A Family of Independent Variable Eddington Factor Methods with Efficient Linear Solvers

We present a family of discretizations for the Variable Eddington Factor (VEF) equations that have high-order accuracy on curved meshes and efficient preconditioned iterative solvers. The VEF discretizations are combined with a high-order Discontinuous Galerkin transport discretization to form an effective high-order, linear transport method. The VEF discretizations are derived by extending the unified analysis of Discontinuous Galerkin methods for elliptic problems to the VEF equations. This framework is used to define analogs of the interior penalty, second method of Bassi and Rebay, minimal dissipation local Discontinuous Galerkin, and continuous finite element methods. The analysis of subspace correction preconditioners, which use a continuous operator to iteratively precondition the discontinuous discretization, is extended to the case of the non-symmetric VEF system. Numerical results demonstrate that the VEF discretizations have arbitrary-order accuracy on curved meshes, preserve the thick diffusion limit, and are effective on a proxy problem from thermal radiative transfer in both outer transport iterations and inner preconditioned linear solver iterations. In addition, a parallel weak scaling study of the interior penalty VEF discretization demonstrates the scalability of the method out to 1152 processors.



There are no comments yet.


page 25


Efficient low-order refined preconditioners for high-order matrix-free continuous and discontinuous Galerkin methods

In this paper, we design preconditioners for the matrix-free solution of...

An Efficient High-order Numerical Solver for Diffusion Equations with Strong Anisotropy

In this paper, we present an interior penalty discontinuous Galerkin fin...

Weak Scaling of DSA Preconditioning of Transport Sweeps using HYPRE

This report summarizes the weak scaling performance of the diffusion-syn...

Towards a Scalable Hierarchical High-order CFD Solver

Development of highly scalable and robust algorithms for large-scale CFD...

Fast Tensor Product Schwarz Smoothers for High-Order Discontinuous Galerkin Methods

In this article, we discuss the efficient implementation of powerful dom...

Matrix-free multigrid block-preconditioners for higher order Discontinuous Galerkin discretisations

Efficient and suitably preconditioned iterative solvers for elliptic par...

A scalable DG solver for the electroneutral Nernst-Planck equations

The robust, scalable simulation of flowing electrochemical systems is in...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The Variable Eddington Factor (VEF) method mihalas, also known as Quasidiffusion (QD) goldin

, is a rapidly-converging nonlinear scheme for solving the Boltzmann transport equation, a crucial component of high energy density physics (HEDP) simulations, nuclear reactor analysis, and medical physics. VEF has been applied to a wide range of transport and multiphysics problems including (but not limited to) nuclear reactor eigenvalue problems

airstova_eigenvalue, nuclear reactor kinetics doi:10.13182/NSE13-42, and thermal radiative transfer (TRT) anistratov1996nonlinear. It performs well in problems having both optically thick and thin regions and treats anisotropic scattering equally well anistratov_fvm; ARISTOVA2000139

. Robust convergence is achieved by iteratively coupling the transport equation to the VEF equations, a moment-based equivalent reformulation of transport. The exact closures used to form the VEF equations are weak functions of the solution meaning even simple iterative schemes, such as fixed-point iteration, converge rapidly.

VEF offers significant algorithmic flexibility in that any valid discretization of the VEF equations will yield a rapidly-converging algorithm. This is in stark contrast to Diffusion Synthetic Acceleration (DSA) which places severe restrictions on the discretization of the moment equations in order to guarantee stability A. In the case where the VEF and transport discretizations are not algebraically consistent, referred to as a VEF method with an “independent” discretization doi:10.1080/00411459308203810; two-level-independent-warsa, the discrete solutions of the transport and VEF equations will differ on the order of the discretization error and will be equivalent only in the limit as the mesh is refined. However, even in an under-resolved problem, VEF still produces a “transport solution” in that the solution of the VEF method is a discrete solution of an equivalent reformulation of the transport equation. Furthermore, VEF methods generally preserve the thick diffusion limit diflim and have conservation even if the transport discretization in isolation does not. These properties are particularly useful in multiphysics calculations since the lower-dimensional VEF equations can be directly coupled to the other physics components in place of the high-dimensional transport equation. In addition, discretizations for the transport and VEF equations can be designed independently so that they are in some sense optimal for their intended uses.

This flexibility has been exploited to improve efficiency in relation to all seven dimensions of the transport equation. GHASSEMI2020109315 showed that different order temporal discretizations can be applied to the transport and VEF equations. Ongoing work suggests that time-stepping stability and accuracy can be maintained when just one transport inversion is performed per time step yee_mc21. anistratov2021implicit used data compression techniques to reduce storage costs in time-dependent calculations. In astrophysics, VEF is used to simplify the implementation of coupling TRT to hydrodynamics and to avoid the memory cost of solving the time-dependent transport equation Jiang_2012; GNEDIN2001437; GEHMEYR1994320. Davis_2012 used a short characteristics discretization of the transport equation. me and LOU2019258 designed a spatial discretization of the VEF equations to increase multiphysics compatibility. YEE2020109696 showed that robust convergence is maintained even when positivity-preserving methods are used inside the iteration. ANISTRATOV2019186 solved the multigroup TRT equations by using a VEF method with multiple levels in frequency. It is also well-known that the multigroup eigenvalue problem can be solved with only the need for eigenvalue iterations on the one-group VEF equations AL.

The above techniques rely on the efficient solution of the discretized VEF equations. VEF methods reduce the overall cost of the simulation by trading inversions of the high-dimensional transport equation with inversions of the lower-dimensional VEF equations. In all of VEF’s applications, the solution of the discretized VEF equations is buried under multiple nested loops corresponding to time integration, Newton iterations, eigenvalue iterations, multi-group iterations, and/or fixed-point iterations. The efficient iterative solution of the VEF equations is then crucial to the efficiency of the overall algorithm and is a prerequisite for the practicality of any VEF method.

The unusual structure of the VEF equations and their lack of self-adjointness make the development of discretizations and their corresponding preconditioned iterative solvers difficult. While considerable effort has been placed into the discretization of the VEF equations, to our knowledge, existing methods either rely on expensive and unscalable preconditioners such as block incomplete LU (BILU) factorization, cannot be solved with interation counts independent of the mesh size, or do not mention solvers entirely. Previous work on discretizing the VEF equations includes finite volume anistratov_fvm; doi:10.1080/00411459308203810; QDBC; Jiang_2012; Jones2019TheQM, finite difference WIESELQUIST2014343, mixed finite element vallette; me; olivier_mandc; LOU2019258, continuous finite element wieselquist; two-level-independent-warsa, and discontinuous finite element dima_dfem techniques. Most VEF methods are designed to be algebraically consistent with their application’s discretized transport equation which typically requires discretizing the first-order form of the VEF equations. Such discretizations solve for both the zeroth and first moment of the solution and thus have significantly more unknowns than discretizations of the second-order form. In addition, block preconditioners benzi_golub_liesen_2005 are required to efficiently solve discretizations of the first-order form. Such solvers can require nested iteration for robustness (see warsa_mfem for a radiation diffusion example).


showed that VEF methods with and without algebraic consistency converge equivalently as long as the transport data is properly represented. In particular, computing the Eddington tensor and boundary factor using finite element interpolation and Discrete Ordinates () angular quadrature enables rapid convergence for any valid discretization of the VEF equations. An independent discretization of the second-order form of the VEF equations then has the potential to provide the rapid convergence of a consistent VEF method while solving for fewer unknowns and avoiding the need for block preconditioners. Such a method also has the flexibility to discretize the VEF equations in a manner that can leverage existing linear solver technology.

Our motivation for this research is in the context of HEDP experiments where the tightly-coupled simulation of hydrodynamics and TRT is required, the latter of which typically includes the transport equation. For hydrodynamics, it has been shown that, compared to low-order methods, high-order methods on curved meshes have improved robustness, symmetry preservation, and strong scaling on emerging high performance computer architectures blast; blast2; blast3. Transport methods compatible with this multiphysics framework are desired. graph_sweeps showed that adequately approximating realistic meshes generated from a high-order hydrodynamics code as straight-edged required a significant number of mesh refinements leading to an impractical increase in transport unknowns. It is also possible that high-order accurate transport methods could be beneficial in terms of multiphysics compatibility with high-order hydrodynamics. High-order transport methods compatible with curved meshes have been developed recently in woods_thesis; graph_sweeps with corresponding consistent DSA discretizations in ldrd_dsa; doi:10.1080/00295639.2020.1799603. However, high-order discretizations of the VEF equations compatible with curved meshes have not yet been developed.

In this paper, we design a family of independent VEF discretizations for the linear, steady-state transport problem that can be efficiently and scalably solved with high-order accuracy, in multiple dimensions, and on curved meshes. Our approach is to begin with discretization techniques known to have effective preconditioners on the simpler case of radiation diffusion (i.e. the model Poisson problem) and adapt them to the VEF equations. By using the Eddington tensor and boundary factor interpolation procedure established in two-level-independent-warsa, these methods achieve both rapid convergence in outer fixed-point iterations and in inner linear solver iterations when paired with a high-order Discontinuous Galerkin (DG) discretization of transport.

In particular, we extend the unified analysis of DG methods developed for elliptic problems presented by Arnold2002 to the VEF equations to derive analogues of the interior penalty (IP), second method of Bassi and Rebay (BR2), Minimal Dissipation Local Discontinuous Galerkin (MDLDG), and continuous finite element (CG) techniques. We show that the IP and BR2 VEF methods are effectively preconditioned by the subspace correction method from Pazner2021 and that Algebraic Multigrid (AMG) is effective for the CG and MDLDG discretizations. dima_dfem also applied DG techniques to the VEF equations but they discretize the first-order form of the VEF equations and only consider lowest-order elements in one dimension. We note that our CG operator is equivalent to extending the continuous finite element VEF discretization in two-level-independent-warsa to multiple dimensions, arbitrary-order, and curved meshes.

The paper proceeds as follows. First, we describe the VEF method analytically and discuss iterative schemes to solve the coupled transport-VEF system. Then, we provide background on representing high-order meshes and finite element solutions and present the mathematical notation that will be used in the remainder of the paper. We derive the extension of the unified framework for DG to the VEF equations. The IP, BR2, MDLDG, and CG VEF discretizations are derived from this framework. Section 6 discusses the design and analysis of the subspace correction preconditioners, and extends their analysis to the case of non-symmetric linear systems. We next give computational results. We show that all the methods presented achieve high-order accuracy on curved meshes through the method of manufactured solutions, preserve the thick diffusion limit both on orthogonal and a severely distorted curved mesh, and are effective on the linearized, steady-state crooked pipe problem, a challenging proxy problem from TRT, in both outer fixed-point iterations and inner linear solver iterations. We also present a parallel weak scaling study for the IP discretization which demonstrates the scalability of the algorithm out to 1152 processors and 10 million VEF scalar flux unknowns. Finally, we give conclusions and recommendations for future work.

2 The VEF Method

The steady-state, mono-energetic, fixed-source transport problem with isotropic scattering and inflow boundary conditions is:


where is the angular flux, the direction of particle flow, the spatial domain of the problem with its boundary, and the total and scattering macroscopic cross sections, respectively, the fixed-source, and the inflow boundary function. The VEF equations are given by


where is the absorption macroscopic cross section, and the zeroth and first angular moments of the angular flux, and


is the Eddington tensor. We refer to as the scalar flux and as the current. In addition, are the angular moments of the fixed-source, . The VEF equations are derived by taking the zeroth and first angular moments of the transport equation and closing the second moment of the angular flux, , with


By eliminating the current, the VEF equations can be cast as a drift-diffusion equation:


In both the first-order form (Eq. 2) and second-order form (Eq. 5), the presence of the Eddington tensor inside the divergence leads to diffusion, advection, and reaction-like terms that make applying existing discretization techniques difficult.

The Miften-Larsen transport-consistent boundary conditions QDBC are




is the incoming partial current computed from the transport boundary inflow function and


is the Eddington boundary factor. This boundary condition is derived by manipulating partial currents and using an analogous nonlinear closure. In equations, with the partial currents defined as ,


where in Eq. 6 plays the role of using the transport equation’s inflow boundary condition.

If the Eddington tensor and boundary factor are known, the VEF equations define the zeroth and first moments of the angular flux. In other words, the VEF equations with Miften-Larsen boundary conditions are an equivalent reformulation of the transport equation. However, this is a trivial closure in that the solution to the transport equation must already be known to define the VEF data. VEF methods rely on the fact that the VEF data are weak functions of the angular flux and thus simple iterative schemes can converge rapidly.

Note that when an independent discretization is used for the VEF equations, the discretized VEF scalar flux and VEF current will not be equivalent to the zeroth and first angular moments of the discrete angular flux; the two solutions will differ on the order of the spatial discretization error. To notationally separate the two scalar flux solutions, we use (varphi) to denote the VEF scalar flux and (phi) as the zeroth moment of the angular flux.

VEF methods seek the solution of the nonlinearly coupled system of equations:


where the drift-diffusion form of VEF is used for brevity. Boundary conditions are specified by Eqs. 1b and 6 for the transport and VEF drift-diffusion equation, respectively. Here, the transport equation’s scattering source is now coupled to the VEF drift-diffusion equation and the data for the VEF drift-diffusion equation are nonlinearly coupled to the transport equation. We have increased the complexity of the problem by adding the VEF scalar flux as an additional unknown and by casting the linear transport problem as nonlinear. However, properties of the VEF data allow this nonlinear, coupled system to be solved more efficiently.



be the streaming and collision operator and VEF drift-diffusion operator, respectively, where indicates a nonlinear dependence on the argument. By linearly eliminating the angular flux, the transport-VEF system is equivalent to


Applying the inverse of the drift-diffusion operator, we see that the solution of the coupled transport-VEF system is the fixed-point:




The fixed-point operator is applied in two stages: 1) solve the transport equation using a scattering source formed from the VEF scalar flux and 2) solve the VEF drift-diffusion equation using the VEF data computed from the angular flux from stage 1).

The simplest algorithm to solve Eq. 14 is fixed-point iteration:


where is an initial guess for the solution. This process is repeated until the difference between successive iterates is small enough. Since the Eddington tensor and boundary factor are weak functions of the angular flux even this simple iteration strategy will converge rapidly.

Iterative efficiency can be improved with the use of Anderson acceleration. Anderson acceleration defines the next iterate as the linear combination of the previous iterates that minimizes the residual . For the storage cost of previous iterates, Anderson acceleration increases the convergence rate and improves robustness. While it is not practical to store multiple copies of the angular flux, it is reasonable to expect that a small set of scalar

flux-sized vectors can be stored. The process of linearly eliminating the transport equation, codified in Eq. 

13, allows the Anderson space to be built from the much smaller scalar flux-sized vectors only. In the case where a subset of the angular flux unknowns are not eliminated, such as when a parallel block Jacobi sweep is used to avoid communication costs or when mesh cycles or reentrant faces are present, the solution vector can be augmented with these un-eliminated unknowns so that they are included in the Anderson space. This is the nonlinear analog to the ideas used for Krylov-accelerated source iteration doi:10.13182/NSE02-14.

In addition, defining the nonlinear residual as


root-finding methods such as Jacobian-free Newton Krylov (JFNK) can be used. JFNK builds a new Krylov space to approximate the gradient of at each iteration meaning information across iterations is not kept. JFNK typically required significantly more evaluations of than Anderson-accelerated fixed-point iteration. Thus, we present results using fixed-point iteration and Anderson-accelerated fixed-point iteration only.

The following sections present the discretizations and solvers needed to efficiently evaluate numerically.

3 Mesh and Finite Element Preliminaries

3.1 Description of the Mesh

Let with be the domain of the problem. Consider the tessellation

with the element in the mesh . Each coordinate of the mesh is represented by a piecewise continuous polynomial. In other words, the mesh itself is a member of an finite element space. This allows representation of curved surfaces and enforces continuity of the mesh coordinates along the interfaces between elements. Figure 0(a) depicts a mesh of two quadratic, quadrilateral elements where the mesh control points labeled 2, 7, and 12 are shared between the two elements to enforce continuity of the shared interior interface between them.

The mesh element is given as the image of the reference element under an invertible, polynomial mapping where for simplicial elements (triangles and tetrahedra) or for tensor product elements (quadrilaterals and hexahedra). Here, is the space of all polynomials of total degree at most in all variables and the space of all polynomials of degree at most in each variable. For example, in two dimensions,




We do not consider the use of on tensor-product elements for either the mesh or the solution.

The reference element is the unit -simplex for simplicial elements (i.e. a triangle with coordinates (0,0), (1,0), and (0,1)) or the unit -cube for tensor product elements. Figure 0(b) depicts a mesh transformation for a non-affine, linear, quadrilateral element. In the remainder of this document, we assume the use of tensor product elements however the derivations apply analogously to simplicial elements.

Figure 1: Depictions of (a) the continuity of an interior face in a high-order curved mesh and (b) the reference transformation for a non-affine, linear, quadrilateral element.

Let denote the reference coordinate. The Jacobian matrix of the mapping is


Further define as the determinant of the Jacobian matrix. As an example, the transformation, Jacobian matrix, and determinant for the transformation depicted in Fig. 0(b) are


The mesh transformations are used to perform integration in reference space using:


For integrands involving gradients, the chain rule implies that


Integration over surfaces is performed over the dimensional reference element using the transformed element of surface area. In this document, integration over the domain is implicitly performed using numerical quadrature and the relations in Eqs. 22 and 23. Finally, the characteristic mesh length, , is computed with


with .

3.2 Finite Element Spaces

On each element, we will seek solutions to the transport and VEF drift-diffusion equations in the space of polynomials mapped from the reference element defined by


where indicates a function defined on the reference element. The dilineation between and is required when non-affine111Examples of non-affine transformations include mapping the reference square to a trapezoid or any high-order, curved element. mesh transformations are used. In such a case, even if . That is, the solution can be non-polynomial due to the composition with the inverse of the element transformation. For example, the inverse of the transformation in Fig. 0(b) is


which is non-polynomial in the first coordinate.

The degree- DG finite element space is:


so that each function is a piecewise polynomial mapped from the reference element with no continuity enforced between elements. Its vector-valued analog is


which simply uses the scalar DG space for each component of the vector. We will also need the discrete , or continuous finite element space, defined as:


Here, is a piecewise continuous mapped polynomial.

A nodal basis for the element-local polynomial space is used. For a degree- element, let denote the Gauss-Lobatto or Gauss-Legendre points in the interval . The points on the unit cube are given by the -fold Cartesian product of the one-dimensional points. Let denote the Lagrange interpolating polynomial satisfying where is the Kronecker delta. The set of functions form a basis for the space . The DG and finite element spaces are built element-by-element from this local basis.

Note that the Gauss-Lobatto points include the interval end points while the Gauss-Legendre points do not. Thus, using Gauss-Lobatto points yields both points on the interior and the boundary of the element while using Gauss-Legendre leads to points on the interior of the element only. These are referred to as closed and open bases, respectively. In the case of DG, no continuity between elements is enforced so it is acceptable to use either an open or closed basis. Both Gauss-Lobatto and Gauss-Legendre have the required properties to be accurate in the limit so the choice of Gauss-Lobatto versus Gauss-Legendre is typically dictated by other aspects of the overall algorithm such as preconditioners. The basis formed from the Gauss-Legendre points has the beneficial property of diagonal mass matrices on affine meshes, while the basis formed from Gauss-Lobatto points typically leads to sparser global systems since closed bases couple fewer unknowns on interior faces. A closed basis is required for finite element spaces to enable the strong enforcement of continuity between elements.

3.3 Mathematical Notation

It is helpful to define the “broken” gradient, denoted , obtained by applying the gradient locally on each element. That is,


This distinction is important since for , is not well-defined since may be discontinuous across element interfaces. However, is well-defined since is locally differentiable on each element.

Figure 2: A depiction of a discontinuous, piecewise quadratic solution across two quadrilateral elements. The normal vector, , is defined as pointing from to along the face between and .

We will use the following notation to describe the jump and average of a discontinuous function along an interior mesh face. Let be the set of all unique faces in the mesh and the set of unique interior faces. Additionally, define as the set of faces on the boundary so that . We define as the outward unit normal to element . On an interior face between elements and , we use the convention that is the the unit vector perpendicular to the shared face pointing from to (see Fig. 2). On such an interior face, the jump, , and average, , are defined as


where with analogous definitions for vectors.

Note that in contrast to the notation of Arnold2002, our jump operator does not change the rank of its argument: the jump of a scalar is a scalar and the jump of a vector is a vector. Consequently, our notation is not invariant under element renumbering, since flipping the ordering of the elements negates the value of the jump. However, the bilinear and linear forms presented in this paper always pair the jump with another normal-dependent term. The negation of the jump induced by swapping the element ordering is then balanced by flipping the orientation of the normal vector, and so the discretizations under consideration are in fact invariant with respect to the element ordering.

On the boundary of the mesh, we set the jump and average to


and likewise for vector-valued functions on the boundary. A straightforward computation shows that


We refer to this as the “jumps and averages identity”. The restriction of the integration to interior faces for the second term in the last equality is consistent with the notation of Arnold2002 and is used so that only one term contributes on the boundary of the mesh.

Finally, we refer to a function as “single-valued” on an interior face if its values obtained from approaching from each side of the face are identical so that


Note in particular that the jump and average operators are single-valued.

4 Transport Discretizations

In this work, we assume the transport equation is discretized with the Discrete Ordinates () angular model and an arbitrary-order Discontinuous Galerkin (DG) spatial discretization compatible with curved meshes (e.g. woods_thesis; graph_sweeps). In , the transport equation is collocated at discrete angles, , and integration is numerically approximated using a suitable angular quadrature rule on the unit sphere. The VEF data are then


where is the discrete angular flux in direction . With degree- DG in space, the angular flux in each discrete direction is a member of . Through the standard finite element interpolation procedure, the Eddington tensor and boundary factor in Eq. 35 can be evaluated at any location in the mesh. Note that it is important to interpolate the numerator and denominator of the VEF data independently. That is, the boundary factor and each component of the Eddington tensor are represented as degree- improper rational polynomials on each element.



as the discrete second moment of the angular flux and using the quotient rule, the local divergence of the Eddington tensor


is well-defined assuming . Here, the divergence of a second-order tensor is the vector formed by taking the divergence of each of the columns of the tensor.

We restrict our attention to problems where inside the domain, for some . This assumption is reasonable for our applications but may be violated in shielding or deep penetration problems. Application of a positivity-preserving negative flux fixup then ensures that is bounded away from zero, so that , , and are all bounded. Thus, through angular quadrature and finite element interpolation, the Eddington tensor, boundary factor, and the local divergence of the Eddington tensor can be evaluated at any point in any element of the mesh. This completes the definition of the connection between the discrete transport equation and the VEF drift-diffusion equation. Note that since the angular flux is generally discontinuous across interior mesh interfaces, the Eddington tensor and its divergence also will be. Thus, we will carefully design the discretization of the VEF drift-diffusion equation to accommodate discontinuous data.

The VEF scalar flux connects with the transport equation in the scattering source. To support generality, we assume that the finite element space for the VEF scalar flux and the finite element space for the angular fluxes are different. The scattering source is then constructed using a mixed-space mass matrix that has test functions in the space for the angular flux and trial functions in the space for the VEF scalar flux.

5 Derivation of DG VEF

In this section, we adapt the derivation of a unified framework for DG methods designed for the Poisson equation from Arnold2002 to the VEF equations. This enables the use of any of the DG methods described there. Arnold2002 derive a family of DG methods for:


with Dirichlet boundary conditions. The present goal is to adapt their derivation to the VEF equations:


with the Robin style boundary conditions given in Eq. 6. We will see significant differences in the final equation since the Eddington tensor is inside the divergence. Additionally, the presence of a right-hand side in the first moment equation as well as non-unit coefficients introduce further complications. We will then derive analogues of the interior penalty (IP), second method of Bassi and Rebay (BR2), and Minimial Dissipation Local Discontinuous Galerkin (MDLDG) variants. Finally, we will show how to extract a continuous finite element method from this framework.

5.1 Adaption of the Unified Framework to VEF

We seek the VEF scalar flux in the degree- DG finite element space and the current in the degree- vector-valued DG finite element space . The weak form is then: find such that for all :


where the numerical fluxes and are approximations of and on the boundaries of the elements in the mesh. We group the product as the numerical flux to mimic the integration by parts of a tensor times a vector. Here, the gradient of a vector is




is the scalar contraction of two tensors. Note that if then


and the symmetric weak form for radiation diffusion is recovered.

Summing the zeroth moment over all elements:


where the jumps and averages identity (Eq. 33) was used. We will now use the discrete first moment to determine a functional form for . Integrating by parts locally over element , we have that


The first moment’s weak form on each element becomes:


Summing over all elements and using the jumps and averages identity, the weak form for the first moment is:


where is evaluated as , and the term is computed using Eq. 37.

We now wish to write all terms as volumetric integrals so that a functional form for the current can be found. To that end, define lifting operators and such that


where and are vector functions that are singled-valued on . Note that the lifting operators are finite element grid functions just as the current is and that the left hand sides are simply the total interaction mass matrix. Since is piecewise discontinuous, the mass matrix is block-diagonal by element and thus the systems of equations corresponding to Eqs. 48a and 48b are amenable to efficient direct factorization (see Appendix A).

Setting and and using the definitions of the lifting operators, Eq. 47 can be written entirely in terms of volumetric integrals as:


for all . Subtracting the right hand side and setting the integrand to zero implies that


Observe that the above represents the element-local strong form of the current, found by analytically eliminating the current, with additional terms that capture the effect of the numerical fluxes. In other words, we have derived the discrete elimination of the current.

Using this discrete form for the current and the definitions of the lifting operators to convert from volumetric integrals back to surface integrals, the zeroth moment becomes:


On boundary faces, we apply the Miften-Larsen boundary conditions by setting


All the methods we consider use so-called conservative numerical fluxes such that


Using the boundary conditions and the assumption of conservative numerical fluxes, Eq. 51 becomes:


Equation 54 defines a family of DG methods. That is, through the specification of the numerical flux for the current on interior faces, analogues of all the methods listed in Arnold2002 can be derived.

5.2 Specification of Numerical Fluxes

All the methods we consider use numerical fluxes of the form


where and are single-valued functions whose purpose is to ensure a stable discretization. The IP, BR2, and LDG methods differ only in the choice of and . With these numerical fluxes, Eq. 54 becomes:


Recall that this form has already applied boundary conditions according to Eq. 52. In other words, the above corresponds to a DG scheme with the following numerical fluxes:


5.2.1 Interior Penalty

An interior penalty (IP)-like method uses


where is the penalty parameter. IP methods require that in order to guarantee stability. The full IP weak form is then: find such that


5.2.2 Br2

The second method of Bassi and Rebay (BR2) uses an alternative penalty term. Let be a face-local lifting operator defined by


Here, is a scalar function that is single-valued on the interior face . Note that the integration on the left hand side is over the entire domain while the right hand side is localized to a single interior face. This means the right hand side, and thus , will be non-zero only for DOFs in elements that share the face .

A BR2-like discretization sets


so that the relevant term is


This BR2 numerical flux avoids the need to tune the penalty parameter while still allowing element-by-element assembly (see Appendix A).

The BR2 DG VEF discretization is then: find such that


5.2.3 Local Discontinuous Galerkin

Finally, we consider the Local Discontinuous Galerkin (LDG) method. In general, LDG uses the following numerical fluxes: