High-Order Finite Element Second Moment Methods for Linear Transport

Abstract We present high-order, finite element–based Second Moment Methods (SMMs) for solving radiation transport problems in two spatial dimensions. We leverage the close connection between the Variable Eddington Factor (VEF) method and SMM to convert existing discretizations of the VEF moment system to discretizations of the SMM moment system. The moment discretizations are coupled to a high-order Discontinuous Galerkin discretization of the SN transport equations. We show that the resulting methods achieve high-order accuracy on high-order (curved) meshes, preserve the thick diffusion limit, and are effective on a challenging multimaterial problem both in outer fixed-point iterations and in inner preconditioned iterative solver iterations for the discrete moment systems. We also present parallel scaling results and provide direct comparisons to the VEF algorithms from which the SMM algorithms were derived.


I. INTRODUCTION
The Second Moment Method (SMM), introduced by Lewis and Miller Jr., [1] is an iterative, moment-based algorithm for solving the Boltzmann transport equation in applications such as nuclear reactor design, medical physics, and high energy density physics.The original concept was to iteratively couple the transport equation to a diffusion equation with linear, transport-dependent source terms that vanish when the solution is a linear function in angle. [2]The transport-dependent source terms "correct" the diffusion approximation such that upon iterative convergence, the corrected diffusion equation can reproduce the transport solution.Using the corrected diffusion solution to compute slow-to-converge physics, such as scattering, provides iterative acceleration, leading to a robust algorithm where the cost scales independent of the mean free path.
SMM can also be viewed as a member of a broader class of multiscale methods, known as moment or highorder low-order methods, developed to solve highdimensional kinetic equations in the context of multiphysics simulations (see Ref. [3] for Chacón et al.'s 60-year history review) where the kinetic equation is iteratively coupled to a small number of its moments through discrete closures.In the case of SMM, the corrected diffusion equation can be viewed as a two-moment moment system closed with additive closures.Such methods are *E-mail: solivier@lanl.govThis material is published by permission of the Los Alamos National Laboratory, for the U.S. Department of Energy under Contract No. 89233218CNA000001.The US Government retains for itself, and others acting on its behalf, a paid-up, non-exclusive, and irrevocable worldwide use in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.automatically asymptotic preserving (i.e., preserve the thick diffusion limit); are conservative even when a nonconservative discretization of the kinetic equation is used; and can be directly coupled to other physics components in the simulation, isolating the expensive, high-dimensional kinetic equation from the evolution of stiff multiphysics.Furthermore, the rapid and robust convergence of moment methods is maintained even when the kinetic and moment equations are discretized independently. [4]These so-called independent moment methods [5] allow the flexibility to choose the discretization of the moment system to meet the requirements of the overall algorithm such as computational efficiency and multiphysics compatibility.
SMM is particularly closely related to another moment method known as the Variable Eddington Factor (VEF) [6] or Quasi-Diffusion [7] method.The SMM and VEF algorithms are differentiated only by their choice of closure: SMM uses additive, linear closures while VEF uses multiplicative, nonlinear closures.VEF's multiplicative closures affect the moment operator (i.e., the left-hand side of the linear system Ax ¼ b), resulting in a non-self-adjoint moment system that has been difficult to discretize and precondition effectively (Ref.[8] and the references therein).SMM's use of additive closures results in a moment system where the closures appear as source terms only, leading to clear computational advantages over VEF.In particular, the self-adjoint structure of the diffusion approximation is preserved, allowing the use of simpler iterative schemes, such as conjugate gradient, to solve the discretized moment system.By contrast, VEF's generally nonsymmetric discrete moment system requires the use of general-purpose solvers, such as GMRES or BiCGStab, which typically require more computation and memory.SMM can also directly leverage existing discretization and preconditioning strategies developed for standard elliptic problems that would not be possible for the unusual structure of the VEF moment system.Furthermore, where the VEF closures are well defined only when the transport solution is strictly positive, the SMM closures are valid for any transport solution, including solutions that possess negativities arising in underresolved calculations.
In the process of performing a Fourier stability analysis of the VEF method, Cefus and Larsen [9] found that SMM can be viewed as a VEF method that has been linearized about a linearly anisotropic solution and that the SMM and VEF algorithms had similar iterative convergence properties.We stress that even though SMM can be viewed as a linearized VEF method, the SMM moment system is still simply an algebraic reformulation of the transport equation and is thus able to achieve the full transport physics even when the problem is not diffusive.The methods' close iterative convergence rates indicate that SMM has the potential to converge as quickly as VEF while using a moment system that is less expensive to invert at each iteration and more robust to negativities in the transport solution.
In this paper, we leverage the close connection between SMM and VEF-both through linearization and by algebraically reformulating the closures-to convert the Discontinuous Galerkin (DG) and mixed finite element independent VEF methods presented in Refs.[8]  and [10], respectively, to the corresponding independent SMMs.The SMM moment discretizations are coupled to a high-order DG discretization of the Discrete Ordinates ðS N Þ transport equations to solve problems from linear transport.
Our motivation for this work is in the context of high energy density physics simulations where the tightly coupled simulation of hydrodynamics and thermal radiative transfer is required, the latter of which typically employs the S N angular model.We are interested in designing discretizations for the SMM moment system that can be scalably solved with existing preconditioning technology and are compatible with the high-order finite element, Lagrangian hydrodynamics framework being developed at the Lawrence Livermore National Laboratory (LLNL), [11] where the thermodynamic variables are represented using discontinuous finite elements.Both the DG and mixed finite element approaches satisfy this compatibility requirement by producing the SMM scalar flux solution in a DG finite element space.The mixed finite element discretization has the added benefit that the SMM diffusion operator exactly matches the mixed finite element techniques used for radiation diffusion calculations at LLNL. [12,13] In such case, the existing diffusion package could be elevated to a transport algorithm by simply supplying a modified, transport-dependent source term to the linear and nonlinear solvers already in place for diffusion, thereby easing the implementational burden associated with multiphysics coupling and maintaining disparate code bases for diffusion and transport.
We note that mixed finite elements are also used for radiation diffusion-based reactor physics calculations [14][15][16] for which the mixed finite element SMM could also serve as a drop-in replacement.We also present a continuous finite element [Continuous Galerkin (CG)] discretization for the SMM moment system as a CG-based VEF algorithm was found to perform as well as DG-based algorithms while solving for fewer moment unknowns and requiring simpler preconditioning strategies. [8]ince the left-hand side of the SMM moment system is equivalent to radiation diffusion, consistent SMMs can be designed by using any of the diffusion discretizations developed for consistent Diffusion Synthetic Acceleration (DSA) [17] methods such as by Warsa et al., [18] Adams and Martin, [19] Wang and Ragusa, [20] and Haut et al. [21] Furthermore, many consistent DSA schemes can be scalably solved with preconditioned iterative solvers, e.g., Refs.[20] and [21].The design of efficient, consistent SMM algorithms then only requires developing consistent discretizations for the SMM correction sources.
Thus, in the case of SMM, the consistent approach is likely to be effective in achieving our research goals of multiphysics compatibility and scalable preconditioned iterative solvers for the discretized moment system.However, such a method would not enjoy the flexibilities discussed above and, in particular, would preclude the possibility of a SMM algorithm acting as a drop-in replacement for an existing diffusion package.
Previous work on discrete SMMs is limited in the literature.Stehle et al. [22] developed a domain decomposition method where SMM was used to couple transport and diffusion domains through interface conditions.The SMM moment system and correction sources were discretized to be algebraically consistent with the simple corner balance [23] S N transport method.Anistratov et al. [24] investigated a multilevel-in-energy algorithm for neutron transport problems.The discretization of the moment system was designed to be consistent with a lowest-order DG discretization of the S N transport equations.This method can be viewed as a linearization of the DG VEF method from Anistratov and Warsa. [25]Neither of these methods attain higher than second-order accuracy or are compatible with curved meshes.Furthermore, independent SMMs have not yet been investigated.
We proceed as follows.First, we introduce a general framework for the iterative moment algorithm in the context of the monoenergetic, fixed-source transport problem with isotropic scattering.We illustrate the iteration in a general context by abstracting away the choice of closure (i.e., VEF versus SMM) as well as casting the moment system in firstor second-order form.We also demonstrate that the SMM and VEF algorithms are related through linearization using the Gateaux derivative.Next, we provide background on the relevant finite element technology used to discretize the transport and SMM moment equations on curved meshes and discuss the discretization overlap that occurs in computing the SMM correction sources and the transport equation's scattering source.We derive DG, CG, Raviart Thomas (RT) mixed finite element, and hybridized Raviart Thomas (HRT) mixed finite element SMM moment system discretizations through their close connection to the corresponding discrete VEF methods derived in Refs.[8] and [10].We then present numerical results.We show that the SMMs converge the scalar flux with optimal accuracy on a manufactured transport problem and preserve the diffusion limit on an orthogonal and a curved mesh.We also compare the performance of the methods on a multimaterial problem and show parallel scaling results.Finally, we give conclusions and recommendations for future work.

II. MOMENT METHODS FOR RADIATION TRANSPORT
We take the monoenergetic, fixed-source transport problem with isotropic scattering given by Eqs.(1a) and (1b) as our model problem: where x 2 R dim and Ω 2 S 2 are the spatial and angular variables, respectively; ψðx; ΩÞ is the angular flux; D � R dim is the spatial domain of the problem with qD its boundary; σ t ðxÞ and σ s ðxÞ are the total and scattering macroscopic cross sections, respectively; qðx; ΩÞ is a fixed source; ψðx; ΩÞ is the inflow boundary function.
In this section, we derive the moment system and its boundary conditions.Because of the streaming term Ω � Ñψ, angular moments always produce more unknowns than equations.We define two types of closures that give rise to the VEF and SMM moment systems.We then describe the iterative schemes used to solve the coupled transport-moment systems.The section concludes with a discussion of the close connection between the VEF and SMM algorithms established by Cefus and Larsen. [9]

II.A. Derivation of Moment System
The moment system comprises the zeroth and first angular moments of the transport equation: where σ a ðxÞ;σ t ðxÞ À σ s ðxÞ = absorption macroscopic cross section; and first moments of the fixed source q; φ, J, and P = HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1181 zeroth, first, and second angular moments of the angular flux, respectively.We refer to the first three moments of the angular flux as the scalar flux, current, and pressure, respectively.Boundary conditions are derived by manipulating partial currents.Let J � n ¼ ò Ω�n_0 Ω � nψdΩ denote the partial currents where n is the outward unit normal to the boundary of the domain.The net current can be expressed in terms of the partial currents and is manipulated to yield Defining the boundary conditions for the moment system are where J in ¼ ð Ω�n < 0 Ω � nψdΩ is the incoming partial current computed using the inflow boundary function ψ.
Combined, the moment system is given by In three dimensions, there are 10 unknowns corresponding to the scalar flux, the three components of the current, and the six unique components of the symmetric pressure tensor but only four equations arising from the scalar zeroth moment and vector first moment equations.In addition, we have only one equation on the boundary but two unknowns corresponding to the normal component of the current and the boundary functional B. Thus, we need to provide additional symmetric tensor and scalar-valued equations on the interior and boundary of the domain, respectively, to form a solvable system of equations.These additional equations are referred to as closures.For the moment methods considered in this paper, the closures are exact such that the moment system is an equivalent reformulation of the transport equation and trivial in that the transport solution must already be known to define the closures.
Remark 1: When the discretized moment and transport equations are not algebraically consistent (e.g., in the context of an independent moment method), the moments of the angular flux and the solution of the moment system will not be equivalent; they will differ on the order of the spatial and angular discretization errors associated with the discrete transport and moment systems.To notationally separate the two scalar flux solutions, we use φ (varphi) to denote the moment system's scalar flux and ϕ ¼ òψdΩ (phi) for the zeroth moment of the angular flux.Since the two solutions differ on the order of the discretization error, we can approximate φ=ϕ � 1 and φ À ϕ � 0 without altering the overall transport algorithm's order of accuracy in space or angle.

II.B. VEF and SMM Closures
We consider two types of closures for the pressure and boundary functional that give rise to the VEF and SMM algorithms.For VEF, the pressure and boundary functional are closed by multiplying and dividing by the scalar flux: where are the Eddington tensor and boundary factor, respectively.Since for the continuous equations, φ=òψdΩ ¼ 1, these closures are simply algebraic manipulations of the

1182
OLIVIER and HAUT On the other hand, SMM uses additive closures of the following form: where are the SMM correction tensor and boundary factor, respectively.Here, I denotes the dim � dim identity tensor.In the same way that the VEF closures are derived by multiplying and dividing by the scalar flux, the SMM closures are formulated by adding and subtracting the scalar flux.Analogously to VEF, this closure procedure is simply an algebraic reformulation of the pressure and boundary functional for the continuous equations.When discretized independently, φ À òψ dΩ is on the order of the discretization error, and by the same arguments used for VEF, this numerical discrepancy can be ignored.Due to their additive formulation, the SMM closures are linear functions of the transport solution that depend on both the magnitude and angular shape of the transport solution.The SMM closures T and β vanish when the transport solution is linearly anisotropic.The factors of 1=3 and 1=2 used in the correction tensor and boundary factor, respectively, cause the SMM moment system to match the radiation diffusion equation with Marshak boundary conditions when the transport solution is linearly anisotropic.
Using these closures in Eqs.(6), the VEF and SMM moment systems are given by and respectively.By eliminating the first moment equation, the moment system can be recast as a second-order partial differential equation.For VEF, the second-order form is while for SMM we have Here, boundary conditions are specified by Eqs.(12c) and (13c) for the VEF and SMM second-order forms, respectively.In both cases, the moment system is equivalent to radiation diffusion with Marshak boundary conditions when the transport solution is linearly anisotropic.Note that the left-hand side of the VEF operator is nonself-adjoint due to the presence of the Eddington tensor inside the divergence while the SMM left-hand side is simply the self-adjoint radiation diffusion approximation.Observe that the SMM correction source, Ñ �

II.C. The Moment Method Algorithm
Moment methods solve the coupled transportmoment system simultaneously.The transport equation is used to provide the VEF or SMM closures while the moment system is used to compute the momentdependent physics.In the present discussion, the moment system's scalar flux solution is used to compute the isotropic scattering source.In this way, the coupling of the angular phase space induced by integrating over angle is avoided, allowing use of the transport sweep to efficiently invert the transport equation's streaming and collision operator.Simple iterative schemes often converge rapidly and robustly due to the closures' weak dependence on the transport solution.
We first introduce notation that abstracts away the choice of the closures and casting the moment system in first-or second-order form.Let Mðψ; XÞ ¼ 0 denote one of the moment systems derived in Sec.II.B with X the moment system's unknowns.For example, Mðψ; XÞ could represent the VEF moment system in first-order form given by Eqs.(12) where X would include both the scalar flux and the current.In the case of the second-order form, we would set X ¼ φ since the scalar flux is the only unknown.For VEF, Mðψ; XÞ is nonlinear in ψ and linear in the moments X.For SMM, Mðψ; XÞ is linear in both arguments.
The moment algorithm solves the coupled system given by where transport boundary conditions are specified in Eq. (1b).The moment system's boundary conditions are given by Eq. (12c) for a VEF method and Eq.(13c) for SMM.Here, the moment system is coupled to the transport equation through the closures, and the transport equation's scattering source is coupled to the moment system through the moment system's scalar flux.We have increased the complexity of the problem by adding the moment system's unknowns.In the case of VEF, the coupled system in Eq. ( 16) is also nonlinear due to the use of nonlinear closures.However, solving the coupled system is still advantageous due to the ability to use the transport sweep and the rapid convergence of the closures. Let: represent the streaming and collision operator.The coupled transport-moment system can then be rewritten as By linearly eliminating the angular flux, the coupled system is equivalent to Observe that Eq. ( 19) is now a function of the moment solution only.That is, we can define and equivalently solve FðXÞ ¼ 0. In this reduced problem, the angular flux appears only as an auxiliary variable used to compute the residual FðXÞ; and we say that the angular flux is determined by the moment system.This reduced formulation FðXÞ ¼ 0 has much lower dimension than the original coupled system given in Eq. ( 16) but has the same solution.Because of this, advanced solvers for FðXÞ can be applied that would otherwise be impractical for Eq. ( 16) due to the storage and computational costs associated with the high dimensionality of the angular flux.We now leverage the structure of the VEF and SMM moment systems to further simplify the above algorithm.Let

1184
OLIVIER and HAUT represent the VEF moment system such that Mðψ; XÞ ¼ VðψÞX À f .We then have that Operating by the inverse of the VEF moment system, the coupled transport-VEF system is equivalent to For SMM, the moment system is of the form where D is a diffusion operator and bðψÞ includes the moments of the fixed-source and the transport-dependent correction sources.The root-finding problem FðXÞ ¼ 0 is then equivalent to Thus, for both VEF and SMM, the solution of the coupled transport-moment system is the fixed point: where GðXÞ is given by for VEF and for SMM.The fixed-point operator G is applied in two stages.Stage 1: Apply a transport sweep to invert the streaming and collision operator on a scattering source formed from the moment system's scalar flux.Stage 2: Solve the moment system using the closures computed with the angular flux from Stage 1.The definitions of the fixed-point operator G show the key differences between the VEF and SMM algorithms.VEF has a transportdependent left-hand-side operator while the right-handside sources are fixed.On the other hand, SMM has transport-dependent sources but a fixed left-hand-side operator corresponding to radiation diffusion.The simplest algorithm to solve X ¼ GðXÞ is fixedpoint iteration: where X 0 is an initial guess.This process is repeated until the difference between successive iterates is small enough.We will also use Anderson acceleration, [26] a more advanced algorithm for solving fixed-point problems that stores a set of k previous evaluations of the fixed-point operator and chooses the next iterate by finding the optimal linear combination of the k operator evaluations that minimizes the residual FðXÞ.For a storage cost of k moment system-sized vectors, Anderson acceleration increases iterative efficiency and robustness.We stress that the fixed-point problem in Eq. ( 29) corresponds to the reduced formulation where the angular flux is determined by the moment system, avoiding the need to store k angular flux-sized vectors.The coupling between the transport and moment equations for the VEF and SMM algorithms is depicted in Fig. 1.Convergence of the fixed-point iteration is expected to be rapid since the closures are weak functions of the transport solution. [7]mark 2: The changing left-hand-side operator in the VEF algorithm means that costs associated with forming a new sparse matrix, such as finite element assembly and solver setup [e.g., lower-upper factorization or the algebraic multigrid (AMG) setup phase], must be incurred at each iteration.Since SMM's left-hand side is iteratively fixed, setup costs need be incurred only once, amortizing them over the entire iteration.Furthermore, less computational work is needed to assemble a right-hand-side vector than is needed to assemble a left-hand-side sparse matrix.Finally, SMM requires solving only the self-adjoint diffusion operator at each iteration whereas VEF requires the solution of a generally nonsymmetric operator.In this way, simpler iterative solvers can be applied to the SMM moment system that would not be convergent for the VEF moment system.For example, we will use conjugate gradient to solve three of the proposed SMM discretizations while the analogous VEF methods require the use of nonsymmetric solvers such as BiCGStab or GMRES.Thus, it is expected that SMM's fixed-point operator will be less expensive to evaluate than VEF's.

II.D. Connection via Linearization
In the process of performing a Fourier stability analysis of the VEF algorithm, Cefus and Larsen [9] showed that SMM is equivalent to the VEF algorithm linearized about a linearly anisotropic solution.Here, we use this connection between VEF and SMM to provide an HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1185 alternate method for deriving the SMM closures.We stress that the pressure and boundary functional are both linear functions of the angular flux.Thus, the following linearizations of the pressure and boundary functional do not possess linearization error and are instead simply algebraic reformulations.We first derive the functional derivatives of the Eddington tensor and boundary factor.We use the Gateaux derivative, [27] a generalization of the directional derivative that supports taking derivatives of functionals.The Gateaux derivative of a functional f evaluated at some function u in the direction of v is defined as Assuming continuity of f , the Gateaux derivative is equivalent to This second definition is preferred as it leads to simpler calculations by leveraging the familiar machinery of the partial derivative.For multivariate arguments, such as when � T for some n > 1, the Gateaux derivative is applied in a manner analogous to the partial derivative where differentiation is applied to one variable at a time with the rest fixed.Using the Gateaux derivative, the first-order Taylor series of f ðuÞ expanded about u 0 is In this way, the Gateaux derivative can be used to apply the action of the Jacobian in the linearization process.
Here, we seek to expand the pressure and boundary functional about a linearly anisotropic solution using the Gateaux derivative.Let ψ 0 be such a linearly anisotropic approximation to the transport problem at hand.Using Eq. ( 31), the Gateaux derivative of the Eddington tensor at ψ 0 in the direction ψ is where we have used the quotient rule and have set Since ψ 0 is linearly anisotropic, E 0 ¼ 1=3I.Thus, the Gateaux derivative of the Eddington tensor is equivalent to where TðψÞ is the SMM correction tensor defined in Eq. (11a).By an analogous calculation, the Gateaux derivative of the Eddington boundary factor is where βðψÞ is the SMM boundary factor defined in Eq. (11b).
We now linearize the pressure and boundary functional using the above functional derivatives of the VEF closures.Let y ¼ ψ φ ½ � T represent the solution of the coupled transport-moment system where, without loss of generality, we consider the moment system to be written in second-order form.Further, we set y 0 ¼ ψ 0 φ 0 ½ � T to be the linearly anisotropic approximation of the coupled system.Here, we consider the moment solution φ to be an independent variable.Using the Gateaux derivative, we can write the first-order Taylor series expansion of the pressure as where Tðψ À ψ 0 Þ ¼ TðψÞ À Tðψ 0 Þ ¼ TðψÞ since the SMM correction tensor is linear and zero when the argument is linearly anisotropic, respectively.Further, we have simplified φ 0 =ϕ 0 � 1 by ignoring terms on the order of the discretization error (see Remark 1).Using E 0 ¼ 1=3I, observe that the extreme left-hand sides and right-hand sides of the linearized pressure exactly match the SMM closure given in Eq. (10a).
Analogously, the first-order Taylor series expansion of the boundary functional is where we have applied the same arguments to simplify βðψ À ψ 0 Þ and φ 0 =ϕ 0 as was used for the pressure.
Recognizing that E b0 ¼ 1=2, the linearized boundary functional is equivalent to the SMM closure in Eq. (10b).
Remark 3: While angular quadrature rules can exactly evaluate E 0 ¼ 1=3I, it is not possible to exactly integrate E b0 with numerical quadrature due to the numerator's integrand jΩ � njψ 0 having a discontinuous derivative.Generally, β, J in , and E b0 must all be computed in the same manner, either analytically or by angular quadrature, to ensure the discrete boundary conditions reduce to the Marshak case when the transport solution is linearly anisotropic.
In the VEF algorithm, the Eddington tensor and boundary factor are the only sources of nonlinearity.Thus, the linearized algorithm is equivalent to the original VEF algorithm with P ¼ Eφ and B ¼ E b φ replaced with T þ 1=3Iφ and β þ 1=2φ, respectively.Such an algorithm is equivalent to SMM.Because of this, SMM can be viewed as a VEF algorithm that has been linearized about a linearly anisotropic solution.
Remark 4: Although SMM can be viewed as a VEF algorithm that has been linearized about a linearly anisotropic solution, we stress that the linearization is actually exact since the pressure and boundary functional are linear functions of the transport solution.In other words, SMM is still simply an algebraic reformulation of the transport equation and is thus accurate even in transport regimes.

III. FINITE ELEMENT PRELIMINARIES
In this section, we provide background on the finite element technology used to discretize the transport and SMM moment systems.We describe the machinery used to represent high-order meshes and the related transformations used to integrate the bilinear and linear forms that comprise the discretizations.Finally, we define the finite element spaces used to approximate the angular flux, scalar flux, and current in subsequent sections.

III.A. Description of the Mesh
Let D � R 2 be the domain of the problem tessellated into a collection of elements T such that where K 2 T is a quadrilateral element in the mesh obtained as the image of the reference square b K¼ ½0; 1� 2 under an invertible, polynomial mapping T : b K !K. Here, the element transformation belongs to ½Q m ð b KÞ� 2 ; where Q m ð b KÞ is the space of polynomials of degree at most m in each variable.We use a nodal basis to describe the element transformation.Let ξ i denote the m þ 1 Gauss-Lobatto points in the interval ½0; 1�.The ðm þ 1Þ 2 points ξ i on the unit square b K¼ ½0; 1� 2 are given by the twofold Cartesian product of the one-dimensional points.Further, let , i denote the Lagrange interpolating polynomial that satisfies forms a basis for Q m ð b KÞ.For each element K 2 T , the mapping is of the form where x 2 K, ξ 2 b K, and fx K;i g ðmþ1Þ 2 i¼1 = control points for element K.A quadratic quadrilateral mesh is depicted in Fig. 2a, where, for example, the left element is described by the control points labeled ð0; 1; 2; 5; 6; 7; 10; 11; 12Þ: Observe that the control points along the interface between the two elements are shared, ensuring that the mesh is continuous.The reference space to physical space transformation for the left element is depicted in Fig. 2b.
We define the set of unique faces in the mesh as Γ with Γ 0 ¼ ΓnqD the set of interior faces in the mesh and Γ b ¼ Γ \ qD the set of boundary faces so that Γ ¼ Γ 0 [ Γ b .We will often work with piecewise functions defined element-wise on the mesh.On an interior face F 2 Γ 0 between two arbitrary elements K 1 and K 2 , we use the convention that n is a unit vector perpendicular to the shared face For a piecewise polynomial u described by u 1 for x 2 K 1 and u 2 for x 2 K 2 , the average, ff�gg, and jump, � ½ � ½ �, of u across the face F ¼ K 1 \ K 2 are defined as with an analogous definition for vector-valued piecewise functions.On boundary faces, the jump and average are defined as and likewise for vectors on the boundary.We note that a piecewise-continuous function u satisfies u ½ � ½ � ¼ 0 and ffugg ¼ u on each interior face.
Finally, we define the "broken" gradient, denoted by Ñ h , obtained by applying the gradient locally on each element.That is, The broken gradient is important for piecewisediscontinuous functions where the global gradient is not well defined due to the discontinuities present at interior mesh interfaces.

1188
OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS

III.B. Integration Transformations
The mesh transformations T are used to facilitate numerical integration on arbitrary elements.Let ξ ¼ ξ η ½ � T 2 b K denote the reference coordinates and x ¼ x y ½ � T ¼ TðξÞ denote the physical coordinates.The Jacobian of the transformation is with J ¼ jFj its determinant.The partial derivatives of the mesh transformation are computed using the basis expansion of the mesh transformations.In other words, where x K;i ¼ x K;i y K;i ½ � T ; b Ñ = gradient with respect to the reference coordinates ξ.In this paper, integration over the domain is implicitly computed in reference space using the following sum: This provides a systematic way to integrate over arbitrary domains composed of arbitrarily shaped elements as well as the use of numerical quadrature rules defined on the reference element b K.We now discuss the transformations used to represent the integrand of Eq. (47) in reference space.For a scalar function u : D !R, denote by b u : b K !R its representation in reference space.The functions u and b u are related by where T À 1 : K !b K denotes the inverse mesh transformation.Integration over the physical element is then equivalent to Using the chain rule, the gradient transforms as In this way, the gradient in physical space can be computed using the Jacobian of the mesh transformation and the gradient in reference space.
For the vector-valued functions used in the RT space, the contravariant Piola transformation is used [28,29] : representations of a vector in physical and reference space, respectively.Members of the RT space have a continuous normal component, a requirement of any Hðdiv; DÞ -conforming space. [30]The contravariant Piola transformation represents vectors on the so-called tangent space spanned by the columns of the Jacobian matrix F. Using this basis, the components of a contravariant vector represent the normal and tangential components of the vector along a face.This facilitates the sharing of degrees of freedom that represent the normal component across interior mesh interfaces, enforcing the normal continuity required by Hðdiv; DÞ.
For a contravariant vector, The gradient transforms as where represents the portion of the transformation of the gradient arising from the derivatives of the Piola transformation.This result is derived by direct computation in HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1189 Ref. [10], Appendix A, where implementation details are also provided.We note that b B is related to the Hessian of the mesh transformation and is thus zero when the mesh transformation is affine such that TðξÞ ¼ Aξ þ b for some A 2 R dim � dim and b 2 R dim independent of ξ.In addition, traceð b BÞ ¼ 0, a result known as the Piola identity. [28]Using the Piola identity, the linearity of the trace, and the invariance of the trace under similarity transformations, the divergence of a contravariant vector transforms as Thus, Combining the results from Eqs. ( 52) and (56) yields where b n is the normal vector in reference space corresponding to the physical space normal n.In other words, the contravariant Piola transformation preserves the normal component.
In this paper, integration is implicitly computed using numerical quadrature on the reference element.Integration over surfaces is performed over the onedimensional reference element using the transformed element of length.

III.C. Finite Element Spaces
Here, we define the finite element spaces used to approximate the transport and SMM moment systems.These finite element spaces are defined on the mesh T or the interior skeleton of the mesh Γ 0 and consist of an element-local function space and a set of interelement matching conditions.The combination of a smooth element-local space and suitable matching conditions allows these finite-dimensional spaces to be subspaces of Sobolev spaces such as L 2 ðDÞ, H 1 ðDÞ, and Hðdiv; DÞ.

III.C.1. Discontinuous Galerkin
Let Q p ðKÞ denote the space of polynomials of at most degree p mapped from the reference element.That is, where b u indicates a function defined on the reference element.The delineation between Q and Q is required when non-affine mesh transformations are used.In such case, In other words, the solution can be nonpolynomial due to the composition with the inverse mesh transformation.
The DG space is a finite-dimensional subspace of the L 2 ðDÞ space defined by the space of square-integrable functions.We denote the degree-p DG space with so that each function u 2 Y p is a piecewise polynomial mapped from the reference element with no continuity requirements enforced between elements.Figure 3 shows the distribution of degrees of freedom in a 3 � 3 mesh for the space Y 1 .The Gauss-Legendre interpolation points are used to define the nodal basis for the local polynomial Fig. 3.A depiction of the distribution of degrees of freedom in the linear DG space.The nodal basis on each element is built from the Legendre nodes to illustrate that degrees of freedom are not shared between elements.

1190
OLIVIER and HAUT

III.C.2. Continuous Finite Element
The continuous finite element (CG) space is a finitedimensional subspace of the H 1 ðDÞ Sobolev space defined as the space of square-integrable functions with squareintegrable gradient.It can be shown that if a piecewise function u 2 C 0 ðDÞ-the space of continuous functions with zero continuous derivatives-and satisfies uj K 2 H 1 ðKÞ for each K 2 T , then u 2 H 1 ðDÞ (see Ref. [30], Sec.3.2.1).In other words, a finite-dimensional subspace of H 1 ðDÞ must be continuous across interior mesh interfaces and have a locally smooth function space on each element.Thus, we take the degree-p continuous finite element space to be Here, V p � H 1 ðDÞ since u 2 V p is continuous and uj K 2 Q p ðKÞ � H 1 ðKÞ for each K.We use a nodal basis for KÞ that employs points on the boundary of the element such as the Gauss-Lobatto points.This allows the sharing of degrees of freedom between elements that strongly enforces the continuity requirements of the discrete space.The distribution of unknowns for the space V 2 is shown in Fig. 4. Observe that degrees of freedom are shared between neighboring elements along the interior interfaces between elements.

III.C.3. Raviart Thomas, Broken RT, and RT Trace Space
The RT space is a finite-dimensional subspace of the Hðdiv; DÞ space, where is the space of square-integrable, vector-valued functions with square-integrable divergence. [31,32]It can be shown that v 2 Hðdiv; DÞ when v satisfies the following: (1) v is locally smooth on each element such that vj K 2 ½H 1 ðKÞ� 2 for each K 2 T and (2) v � n is continuous across each interior face.It is also desirable that the divergence of a vector in the degree-p RT space belongs to the degree-p Figure 5a depicts the degrees of freedom in RT 1 on a 3 � 3 mesh.The black circles and red squares denote We also define two related spaces that arise in hybridizing a RT discretization.The first is the broken RT space c RT p , defined as the degree-p RT space without the requirement of continuity in the normal component.That is, Note that RT p � c RT p and that v 2 c RT p belongs to RT p if and only if ½½ v:n �� ¼ 0 on all interior mesh faces.For the broken space, the degrees of freedom are distributed as in Fig. 5a; however, degrees of freedom are not shared between elements.Instead, these degrees of freedom are duplicated so that the normal component is doubly defined on interior interfaces.
Second, we define the interior trace of the RT space.This space represents the normal component of RT p along the interior faces in the mesh.Let P p ð b F Þ be the space of univariate polynomials with degree at most p defined over the reference line b F ¼ ½0; 1� and the associated polynomials mapped to physical space.The trace space is then The interior trace space for RT 1 , Λ 1 , is depicted in Fig. 5b.Note that these degrees of freedom are exactly the degrees of freedom corresponding to the normal component of RT 1 on the interior faces of the mesh.

IV. TRANSPORT DISCRETIZATION
We discretize the S N transport equations with a highorder DG spatial discretization compatible with curved meshes such as from Woods [33] or Haut et al. [34] In S N , the transport equation is collocated at discrete angles Ω d , and integration over the unit sphere is approximated using a suitable angular quadrature rule fΩ d ; w d g N Ω d¼1 .We use ψ d to denote the angular flux in discrete direction Ω d such that ψ d ðxÞ ¼ ψðx; Ω d Þ.The DG discretization approximates each ψ d using the degree-p DG finite element space Y p .Through finite element interpolation, ψ d ðxÞ can be evaluated at any point in the mesh.The SMM closures are computed using the S N angular quadrature and finite element interpolation with In light of Remark 3, we elect to use E b0 , as defined in Eq. ( 36), in place of the factor of 1=2 that scales the zeroth moment term in the boundary correction.Here, E b0 is approximated with S N quadrature according to With this discrete definition, the boundary conditions for the continuous system from Eq. ( 13c) are now modified to where β is defined using the discrete definition in Eq. (69b).
We also approximate the inflow partial current J in using S N quadrature with J in ðxÞ ¼ In the limit as N Ω ! 1, E b0 !1=2; and the discrete Ω � nψΩ so that the continuous boundary condition is recovered.This discrete prescription ensures that the boundary correction factor will vanish when the discrete transport solution is a linear function in angle and was found to be a requirement to obtain the correct moment solution and preserve the optimal convergence rate of the moment algorithm in the thick diffusion limit.Note that we use the same notation to denote the analytic and discrete SMM closures for notational simplicity.For the remainder of this paper, only the discrete definitions of the closures given by Eqs.(69a) and (69b) will be used.
Where needed, we compute the broken divergence of the SMM correction tensor using the finite element interpolation operator for T. In other words, where P and ϕ are the second and zeroth moments of the discrete angular flux, respectively, and the broken derivatives are applied to the finite element interpolation operators used to represent the angular flux on each element.
Remark 5: For VEF, the discrete closures are locally represented in space as ratios of polynomials and are thus rational functionals (see Ref. [8], Sec. 4).The bilinear and linear forms that comprise the VEF moment system discretization will then contain numerical integration error since rational functions cannot be integrated exactly with numerical quadrature.For SMM, the linear closures simply inherit the finite element representation used for the angular flux.In this way, the SMM moment system can be integrated exactly with numerical quadrature.
Remark 6: The VEF closures are normalized and thus are only well defined for strictly positive angular fluxes, requiring the use of a negative flux fixup for robustness.The SMM closures are not normalized and are thus well defined even when the transport solution contains zeros or nonphysical negativities arising from numerical error.This could be particularly beneficial in shielding or other deep-penetration problems where the solution is significantly attenuated.Here, we leverage this behavior to build an algorithm that remains robust to numerically induced negativities both with and without the use of a negative flux fixup.
The SMM scalar flux is coupled to the transport equation through the scattering source.We use a mixed-space scattering mass matrix to support the use of differing finite element spaces for the SMM scalar flux and the angular flux.In other words, the scattering mass matrix is assembled using test functions from the space used to approximate the angular flux and trial functions from the space used to approximate the SMM scalar flux.

V. DISCRETE SMMs
In this section, we leverage the close connection between VEF and SMM to systematically convert the existing VEF methods derived in Refs.[8] and [10] to discretizations of the SMM moment system.In particular, we derive analogs of the Interior Penalty (IP) DG, Continuous Finite Element, RT mixed finite element, and HRT mixed finite element methods.The IP method is derived by linearizing the IP VEF discretization about a linearly anisotropic solution.The RT mixed finite element method is derived by algebraically manipulating the RT VEF discretization.The continuous finite element and HRT methods are found by leveraging their close connection to the IP and unhybridized RT methods, HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1193 respectively.We have elected to derive the IP and RT methods via linearization and algebraic manipulation, respectively, to give examples of these two equivalent strategies.However, either method can be used to convert any discretization of the VEF moment system to a discretization of the SMM moment system.

V.A. Discontinuous Galerkin and Continuous Finite Element
Here, we linearize an existing DG VEF method about a linearly anisotropic solution to derive an analogous SMM moment system discretization.Note that since the transport equation is linear, the linearization process does not alter the transport equation, and thus, we need only to linearize the VEF moment system to derive the SMM moment system.Let Vðψ; φÞ ¼ f represent the DG VEF moment system.Here, V is nonlinear in ψ, due to the nonlinear Eddington tensor and boundary factor, but linear in φ.Further, we set y 0 ¼ ½ψ 0 ; φ 0 � to be the linearly anisotropic approximation to the coupled transport-VEF system.The linearized operator is then where D½��ð�; �Þ denotes the Gateaux derivative defined in Eq. (31).We have simplified Vðψ 0 ; φ 0 Þ þ D½Vj ψ¼ψ 0 �ðφ 0 ; φ À φ 0 Þ ¼ Vðψ 0 ; φÞ using the linearity of V in the second argument.Furthermore, the directional derivative of the Eddington tensor satisfies D½E�ðψ 0 ; ψ À ψ 0 Þ ¼ D½E�ðψ 0 ; ψÞ with an analogous result for the directional derivative of the Eddington boundary factor.Since ψ 0 is linearly anisotropic, Vðψ 0 ; φÞ represents a radiation diffusion system.The correction sources are then given by D½Vj φ¼φ 0 �ðψ 0 ; ψÞ: the directional derivative of the VEF operator evaluated at ψ 0 in the direction ψ; where φ is fixed at the linearly anisotropic moment solution φ 0 .We proceed by evaluating the IP VEF method from Ref. [8], Sec.5.2.1, at a linearly anisotropic transport solution ψ 0 to define the diffusion discretization and then derive the correction source by applying the Gateaux derivative to each term in the VEF discretization.
From Ref. [8], Eq. ( 59), the IP VEF discretization is as follows: Find φ 2 Y p such that where κ is the penalty parameter.This discretization is derived by approximating the scalar flux and each component of the current with a degree-p DG finite element space.The first moment equation is integrated-by-parts twice in order to derive a discrete elimination of the current.The numerical flux is chosen so that the current can be eliminated on each element, leading to a discretization of the second-order form of the VEF equations.The symmetric, positive semi-definite penalty bilinear form ½ �ds is added to stabilize the discretization.It is well known that κ / σ À 1 t p 2 =h is required to ensure that the penalty bilinear form dominates the negative definite terms in the discretization, resulting in an overall linear system that is positive definite and stable with respect to h and p. [35,36] We note that the penalty parameter's proportionality constant is often problem dependent.While the penalty term is required for stability, the presence of the penalty term degrades the performance of multigrid preconditioners, and thus, specialized solvers, such as the uniform subspace correction preconditioner developed by Pazner and Kolev, [37] must be used to achieve a preconditioned iterative solver that converges independent of the mesh size, polynomial degree, and penalty parameter.
Olivier et al. [8] also present discretizations of the VEF moment system that are analogs of the second method of Bassi and Rebay (BR2) [38] and the Local Discontinuous Galerkin (LDG) method. [39]These methods use alternate stabilization techniques that avoid the problem-dependent prescription of the penalty parameter.However, the BR2 and LDG stabilization terms are more expensive to compute (and more complicated to implement).Only minor differences in solution quality and iterative efficiency were seen among IP, BR2, and LDG, and CG.Thus, for clarity of presentation, we present only the IP and CG methods.
Evaluating the VEF data when the angular flux is linearly anisotropic gives

1194
OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS where E b0 is computed using the S N quadrature rule following Eq.( 70).The use of numerical quadrature to evaluate this term is motivated by Remark 3. The diffusion operator is then where D ¼ 1 3σ t is the diffusion coefficient.The above corresponds to a standard IP discretization of radiation diffusion with Marshak boundary conditions.Note that the right-hand side includes additional sources depending on Q 1 since we have not made the assumption that the transport equation's fixed source q is at most linearly anisotropic, as is customary for the diffusion approximation.Inclusion of the first moment of the source facilitates investigating the accuracy of the method with the method of manufactured solutions (MMS) on a transport problem as performed in Sec.VI.A. Next, we must determine the correction terms by computing the directional derivative of each term in the VEF discretization [Eq.( 74)] with φ ¼ φ 0 fixed.For terms without VEF data, the directional derivative is zero as they do not depend on the transport solution.For the boundary factor-dependent term, the correction term is where βðψÞ is the correction factor defined in Eq. (69b) that uses S N quadrature to compute the angular moments of the discrete transport solution.We have used D½E b �ðψ 0 ; ψÞ ¼ βðψÞ=ϕ 0 and have made the approximation that φ 0 =ϕ 0 � 1 (see Remark 1).Analogously, for terms with the Eddington tensor, where TðψÞ is the correction tensor from Eq. (69a) and we have again made the approximation that φ 0 =ϕ 0 ¼ 1.Thus, the correction terms can be derived by setting terms without angular flux dependence to zero and replacing The IP SMM discretization is then: Find φ 2 Y p such that where the local divergence of the correction tensor is computed with Eq. (72).The above represents a standard IP discretization of diffusion with Marshak boundary conditions that is corrected by transport-dependent volumetric and boundary source terms.
A continuous finite element discretization can be extracted from Eq. ( 80) by setting the test and trial spaces to V p , the degree-p continuous finite element space.For any function v 2 V p , In addition, v 2 V p satisfies Ñ h v ¼ Ñv (Ref.[30], Proposition 3.2.1).In other words, the space V p has the required continuity to allow the broken and global gradients to be equal for any v 2 V p .Since DG is used for the transport discretization, the correction tensor is still generally discontinuous across interior mesh interfaces such that Tn ½ � ½ � is nonzero.Thus, the CG SMM moment discretization is: Find φ 2 V p such that HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1195 The CG SMM moment discretization can also be derived directly from the CG VEF moment discretization in Ref.

V.B. Raviart Thomas and HRT
From Ref. [10], Eq. ( 61), the RT VEF discretization is: Find ðφ; JÞ 2 Y p � RT p such that where v; J 2 RT p are implicitly represented with the contravariant Piola transformation.Here, the presence of the Eddington tensor inside the strong form of the first moment's divergence term Ñ � Eφ ð Þ prevents the straightforward application of standard mixed finite element techniques.To illustrate this, we multiply this divergence term by a test function v 2 RT p and integrate over a single element K: where c Eφn is the numerical flux for the product of the Eddington tensor and VEF scalar flux that arises due to the discontinuous approximations used for the angular flux, and thus the Eddington tensor, and VEF scalar flux.
When assembled into a global system, the elementboundary term gives rise to an interior face term of the form along with a contribution on the boundary of the domain related to the boundary conditions that we omit for clarity of presentation.Note that v 2 RT p has a continuous normal component but discontinuous tangential components.Thus, when E ¼ 1 3 I, the interior face term vanishes since For a general Eddington tensor, however, this term remains, requiring the application of DG-like techniques to treat the discontinuity in v � En.
We also note that the volumetric term ð Ñ h v : Eφdx requires computing the gradient of an RT p vector, a procedure significantly complicated by the use of the Piola transformation (see Ref. [10], Sec.3.2 for more details).Note that Olivier and Haut [10] also present a VEF discretization where each component of the current is approximated with continuous finite elements.However, such a method could not be scalably preconditioned and solved even for a radiation diffusion problem where E ¼ 1=3I and E b ¼ 1=2.Thus, we do not consider this method here.
We first undo the choice of numerical flux and the application of the boundary conditions.We chose the following numerical flux for the product of the Eddington tensor and VEF scalar flux: This choice allows the RT VEF discretization to limit to a standard RT discretization of radiation diffusion in the thick diffusion limit.Noting that VEF closes the pressure as P ¼ Eφ, we write the numerical flux term in terms of a numerical flux for the pressure denoted by b Pn.By solving for φ in Eq. (12c), the boundary scalar flux is given by so that

1196
OLIVIER and HAUT where we have substituted the pressure in place of the product Eφ.In this form, the SMM discretization can be derived by defining the pressure on the interior of the domain, the numerical flux of the normal component of the pressure on interior faces, and the boundary conditions.On the interior of the domain, we use the SMM closure P ¼ T þ 1=3Iφ.The volumetric term is then The HRT discretization can be derived by directly applying the standard hybridization process (see Ref. [13]  and the references therein) to the RT SMM moment system.That is, the current is approximated in the broken RT space c RT p , and continuity of the current in the normal component is enforced with a Lagrange multiplier defined on the trace space of RT p , Λ p .The use of a discontinuous approximation for the scalar flux and current allows eliminating these variables locally on each element in favor of the Lagrange multiplier.This reduced system is both significantly smaller than the block system and can be effectively preconditioned with AMG.Once the Lagrange multiplier has been solved for, the scalar flux and current can be found through elementlocal back substitution.A derivation of a hybridized radiation diffusion method is outlined in Ref. [10], Sec.7.1.
The hybridized system is: Find ðJ; φ; λÞ 2 c RT p � Y p � Λ p such that

HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1197
where is the source term for the RT SMM discretization.As noted in Ref. [13], hybridization can be viewed as an algebraic process similar to static condensation.Furthermore, it was shown in Ref. [10], Sec.7.1, that the weak continuity condition in Eq. (95c) implies that the normal component of the current is continuous in a pointwise, strong sense.Thus, the unhybridized RT and HRT discretizations are equivalent.The details of an efficient implementation are discussed in Ref. [10], Sec.7.3, and Ref. [13].

VI. RESULTS
We now present numerical results demonstrating the iterative efficiency and computational performance of the discrete SMMs.The SMMs are compared to the VEF algorithms from Olivier et al. [8] and Olivier and Haut. [10]he methods were implemented using the Modular Finite Element Method (MFEM) finite element framework. [40]We use MFEM's conjugate gradient and BiCGStab implementations along with BoomerAMG from the hypre sparse linear algebra library. [41]KINSOL, from the Sundials package, [42] provided the Anderson-accelerated fixed-point solvers.As described in Ref. [42], Sec. 2, the fixed-point iteration is terminated when the maximum norm of the difference between successive iterates is below the iterative tolerance.We use the high-order DG S N transport solver from Haut et al. [34] to invert the streaming and collision operator needed in the moment algorithms.We refer to the methods as IP, CG, RT, and HRT for the interior penalty, continuous finite element, Raviart Thomas, and hybridized Raviart Thomas methods, respectively.For IP, we use following the standard prescription used for model elliptic problems [35,43] unless otherwise specified.
The IP methods are preconditioned using the subspace correction method introduced in Ref. [37] and extended to the nonsymmetric VEF system in Ref. [8], Sec. 6.This method is a two-stage preconditioner that applies AMG to the DG system assembled onto a continuous finite element space and uses a simple Jacobi smoother to attack errors present in the degrees of freedom located on mesh interfaces.The CG and HRT methods are preconditioned by AMG.Finally, the RT methods are preconditioned with block preconditioners. [44]The RT methods admit a block system of the form where M t = total interaction mass matrix; M a = absorption mass matrix; D and G = discrete divergence and gradient, respectively; f and g = source terms.For SMM, G ¼ À 1=3D T .For VEF, G depends on the Eddington tensor, and thus, G is in general not proportional to the transpose of the divergence matrix D. We use block diagonal and lower block triangular preconditioners of the form where b t G is an approximate Schur complement formed using the lumped total interaction mass matrix b M t computed by summing the off-diagonal entries of M t into the diagonal and setting the off-diagonal entries to zero.Lumping an RT p mass matrix produces a diagonal matrix, and thus, b M À 1 t can be formed efficiently without fill-in.The approximate inverse of the diagonal preconditioner is applied with approximate inverses of the diagonal blocks while the lower block triangular inverse is applied with block forward substitution.We use a single iteration of Gauss-Seidel and AMG to approximately invert the total interaction mass matrix and approximate Schur complement, respectively.
Note that for MINRES to effectively solve the RT SMM system, the system must be represented in a symmetric indefinite form and the preconditioner must be symmetric positive definite.Thus, for RT SMM, we multiply the top row of Eq. (98) by negative three so that the off-diagonal blocks are symmetric and the scaled M t is negative definite.For the block diagonal preconditioner, M t is scaled by positive three so that P diag is symmetric, positive definite.Finally, we note that the lower block triangular preconditioner is a nonsymmetric operator and thus cannot be used with MINRES.More details on the preconditioning strategies used to solve the moment systems are provided in Ref. [8], Sec. 6, and Ref. [10], Sec.6.4.

1198
OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS

VI.A. Method of Manufactured Solutions
The accuracy of the methods is determined with the MMS.The solution is set to where Here, δ ¼ 1:25 is used to ensure ψ > 0; and ζ ¼ 0:1 and ω ¼ 0:05 are used to test spatially dependent, nonisotropic inflow boundary conditions.The domain is D¼ ½0; 1� 2 .With these definitions This leads to an exact Eddington tensor E ¼ P=ϕ that is dense and spatially varying.The MMS ψ and ϕ are substituted into the transport equation to solve for the MMS source q that forces the solution to Eq. ( 100).The accuracy of the SMM discretizations is investigated in isolation by computing the SMM closures from the MMS angular flux and setting the sources Q 0 and Q 1 to the moments of the MMS source.This is accomplished by projecting the MMS angular flux onto a degree-p DG finite element space and using Level Symmetric S 4 angular quadrature to compute the SMM closures, the moments of the MMS source, and the inflow partial current J in .The SMM moment systems are then solved as if T, β, Q 0 , Q 1 , and J in are given data.Errors are calculated with the L 2 ðDÞ norm for scalars and the ½L 2 ðDÞ� 2 norm for vectors given by and HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1199 respectively.We also use the L 2 ðDÞ projection operator Π p : L 2 ðDÞ !Y p such that for some v 2 L 2 ðDÞ.In particular, Π p is used to project the exact MMS scalar flux onto a Y p finite element grid function in order to investigate a superconvergence property of mixed finite elements.
We use refinements of a third-order mesh created by distorting an orthogonal mesh according to the velocity field of the Taylor-Green vortex.This mesh distortion is generated by advecting the mesh control points with where the final time T ¼ 0:3π and is the analytic solution of the Taylor-Green vortex.Three hundred forward Euler steps were used to advect the mesh.An example mesh is shown in Fig. 6.
Figure 7 shows the error in the scalar flux on this MMS problem for a range of finite element polynomial degrees.All four methods achieve optimal Oðh pþ1 Þ convergence.The IP and CG, and RT and HRT methods have nearly identical error behavior, respectively, with the RT and HRT methods having a lower constant.This is expected since the RT and HRT methods solve for both the scalar flux and current and thus do more computational work than IP and CG.
Table I investigates a superconvergence property of mixed finite elements and the error behavior of the current.The order of accuracy and constant for the RT and HRT methods was experimentally determined by computing the error on four mesh sizes and applying logarithmic regression.We compare the scalar flux accuracy, the scalar flux accuracy when the exact solution is first projected onto Y p using the L 2 ðDÞ projection operator Π p , and the accuracy of the current.As seen in Fig. 7, the scalar flux converges optimally.For standard mixed finite element radiation diffusion, the projected error measure should superconverge at Oðh Þ; and the current should be Oðh pþ1 Þ.However, on this quadratically anisotropic transport problem, the superconvergence property is lost as the projected error measure converges at only Oðh pþ1 Þ.Furthermore, the current converges at only Oðh p Þ for odd polynomial degrees and Oðh pþ1=2 Þ for even.This behavior was also seen in Ref. [10] for the RT Fig. 6.A depiction of a third-order mesh generated by distorting an orthogonal mesh according to the Taylor-Green vortex.Refinements of this mesh are used in calculating the error with the MMS.

1200
OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS and HRT VEF methods.We note that the RT and HRT SMM moment discretizations produce solutions that are equivalent to machine precision.This differs from the behavior observed for the VEF methods where RT and HRT differed on the order of the discretization error. [10]his may be due to the ability to exactly integrate the closure terms whereas with VEF the closures cannot be integrated exactly (see Remark 5).

VI.B. Thick Diffusion Limit
The iterative convergence rates of the SMM algorithms are investigated in the thick diffusion limit.The material data are set to where � 2 ð0; 1� is a parameter that characterizes the mean free path and the thick diffusion limit corresponds to the limit � !0. We use two coarse meshes that do not resolve the mean free path to stress test the convergence of the SMM algorithm.The first is an orthogonal 8 � 8 mesh with D¼ ½0; 1� 2 .The second is the triple point mesh shown in Fig. 8. On the triple point mesh, the angular is only approximately inverted due to the presence of reentrant faces that give rise to upper block triangular entries in the streaming and collision operator.The upper block triangular portion is iteratively lagged so that the classical transport sweep can be applied.Because of this approximation, iterative convergence of the moment algorithms is expected to degrade on the triple point mesh.On both meshes, we use Level Symmetric S 4 angular quadrature.The SMM and VEF algorithms are compared using p ¼ 2. The coupled transport-moment systems are solved with fixed-point iteration to a tolerance of 10 À 6 .
Table II shows the number of iterations to convergence on the 8 � 8 orthogonal mesh as � !0. Compared to VEF, SMM converged one to two iterations slower.All of the SMMs required the same number of iterations, indicating the algorithm is insensitive to the choice of moment discretization on this single material problem.Lineouts of the solution are Fig. 8.A depiction of the triple point mesh used to stress the VEF algorithms on a severely distorted, third-order mesh.This mesh was generated with a Lagrangian hydrodynamics simulation.shown in Fig. 9 demonstrating that the methods are converging to the physically realistic diffusion solution as � !0. Iteration counts and lineouts are presented in Table III and Fig. 10, respectively, for the thick diffusion limit test on the triple point mesh.Following the prescription of the penalty parameter used for IP VEF on this problem given in Ref. [8], Eq. ( 94), the IP SMM penalty parameter was scaled by a factor of 169 to ensure stability of the moment system's discretization on the highly distorted triple point mesh.Here, the IP and CG methods and the RT and HRT methods converge equivalently, but unlike on the orthogonal mesh problem, these two sets of methods show different performance for the intermediate values of ε for both VEF and SMM.In particular, RT and HRT SMM required five more iterations to converge compared to IP and CG for both � ¼ 10 À 2 and � ¼ 10 À 3 .This behavior is also seen for the RT and HRT VEF methods.Generally, convergence was slower on the triple point mesh than on the orthogonal mesh.However, this discrepancy is reduced as � !0. The lineouts indicate that all methods attain the nontrivial, diffusion solution.Nonmonotonic oscillations are observed at the center of the solutions due to imprinting of the highly distorted mesh.Qualitatively, the RT and HRT methods show more oscillations, suggesting they are more sensitive to mesh distortion than the IP and CG methods, which might also explain their slower convergence compared to IP and CG on intermediate values of �.

VI.C. Crooked Pipe
We now show convergence in outer fixed-point iterations and inner preconditioned linear solver iterations on a more realistic, multimaterial problem.The geometry and materials are shown in Fig. 11.The problem consists of two materials, the wall and the pipe, which have an 1000× difference in total interaction cross section.Time dependence is mocked by including artificial absorption and sources that correspond to backward Euler time integration.The time step is set so that cΔt¼ 10 3 and the initial condition is ψ 0 ¼ 10 À 4 .The absorption and source are then σ a ¼ 1=cΔt¼ 10 À 3 1 cm and q ¼ ψ 0 =cΔt¼ 10 À 1 1 cm 3 � s � str .The boundary conditions are set so that isotropic inflow of magnitude 1=2π enters on the left entrance of the pipe with vacuum on all other surfaces.A Level Symmetric S 12 angular quadrature set is used.The zero and scale negative flux fixup, [45] a sweep-compatible fixup that zeros out negativities and rescales so that balance is preserved, is used inside the transport sweep to ensure positivity of the angular flux unless otherwise noted.Timing data are presented as the minimum time recorded across five repeated runs.The efficiencies of the outer fixed-point and inner preconditioned linear iterations are investigated by refining in h and p on a uniform mesh of quadrilateral elements that are aligned with the materials in the problem.The outer solver is Anderson-accelerated fixed-point iteration with a small Anderson space of size two.Anderson acceleration is not required for convergence but does provide more uniform convergence with respect to refining in h.The outer and inner iterative tolerances were 10 À 6 and 10 À 8 , respectively.The IP methods were preconditioned by the uniform subspace correction method, CG and HRT with AMG, and RT with block preconditioners.The nonsymmetric lower block triangular and symmetric block diagonal preconditioners were used for the RT VEF and RT SMM moment systems,

HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1203
respectively.The VEF moment systems were solved with BiCGStab while the IP, CG, and HRT SMM moment systems were solved with conjugate gradient.The RT SMM moment system was solved with MINRES.Note that the lower block triangular preconditioner used to precondition the RT VEF moment system is more expensive, and thus more effective, than the block diagonal preconditioner used for RT SMM.However, the lower block triangular preconditioner cannot be used with MINRES since it is not symmetric.The previous outer iteration's moment solution was used as an initial guess for the inner iteration so that the initial guess becomes progressively more accurate as the outer iteration converges.
Table IV presents the number of Andersonaccelerated fixed-point iterations until convergence for the four SMM algorithms, each compared to the corresponding VEF algorithm with an analogous discretization for the moment system.Convergence was identical for the IP and CG, and RT and HRT SMM methods, respectively.Compared to VEF, SMM converged between two iterations faster and three iterations slower.
The average number of inner iterations per outer iteration is presented in Table V.All methods were scalable in h and p. Recall that different solvers were used to solve the VEF and SMM moment systems for each discretization type.In particular, BiCGStab performs roughly twice as much work per iteration as conjugate gradient, and thus, the IP, CG, and HRT SMM moment solves need only take fewer than twice as many iterations as the corresponding VEF moment system to be cheaper in cost.This is corroborated by Table VI, where the average amount of time spent solving the moment systems is shown.Here, the IP and CG SMM moment solves are shown to be cheaper than the corresponding VEF moment solve, especially for p ¼ 2 and p ¼ 3. The VEF and SMM HRT moment solves were roughly equal in cost.The RT SMM moment solve was more expensive than the RT VEF moment solve due to the use of differing preconditioners.Block diagonal-preconditioned MINRES applied to the symmetric SMM moment operator did not outperform lower block triangular-preconditioned BiCGStab applied to the nonsymmetric VEF moment system.This indicates that it might be advantageous to use lower block triangular preconditioning and BiCGStab to solve the RT SMM moment system even though the system is symmetric.
Table VII shows the average cost per outer iteration of assembling the VEF and SMM moment systems.Recalling Remark 2, SMM is expected to have lower assembly costs due to the reduced cost associated with assembling a linear form versus a bilinear form.For the most refined problem with p ¼ 3, assembling the VEF moment system was 2.7×, 4.3×, 1.8×, and 2.2× as expensive as assembling the SMM moment system for the IP, CG, RT, and HRT discretizations, respectively.Assembly is a much larger cost than solving the resulting linear system iteratively.For example, the assembly cost for IP VEF ranges between 3.4× and 6.3× more expensive than the solve cost.For IP SMM, this ratio ranged only between 1.5× and 3.4×, suggesting that assembly occupies a relatively smaller portion of the IP SMM cost at each iteration compared to IP VEF.

1204
OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS Next, we discuss the total cost of solving the crooked pipe problem with VEF and SMM.Table VIII shows the total cost of the full VEF and SMM algorithms for each discretization type on refinements of the crooked pipe.For a given discretization type, the algorithm that required the fewest sweeps was cheapest.On the most refined problem with p ¼ 3, the discretization choice and closure choice led to a relative standard deviation of only 10% in total run time.In other words, because of the high cost of the transport sweep, the algorithmic choices investigated here lead to only a small change in total run time.
In light of Remark 6, we investigate the sensitivity of the SMMs to the use of a negative flux fixup.Table IX presents the number of Anderson-accelerated fixed-point iterations to convergence for the four SMMs with and without the negative flux fixup.Use of the fixup resulted in equivalent convergence to within � 1 iterations indicating that SMM can converge robustly without a negative flux fixup.

VI.D. Spatial Convergence to S N Solution
In this section, we compare the solutions generated by IP SMM and a reference transport solution taken to be the high-order DG S N discretization and DSA preconditioner of Haut et al. [21] as the mesh is refined.We use the   we use refinements of a nonuniform mesh built from a tensor product of the one-dimensional Chebyshev points.The clustering of the Chebyshev points at the end points allows such a mesh to capture the boundary layer in the thick diffusion limit problem.An example mesh is shown in Fig. 12a. Figure 12b shows the L 2 ðDÞ norm between the SMM and S N solutions as a function of h max , the maximum characteristic mesh length in the mesh.We use Level Symmetric S 4 angular quadrature and p ¼ 2. The solutions are computed on four meshes built from 61, 81, 101, and 121 Chebyshev points in each direction.Using p ¼ 2, Eq. ( 110) predicts that the two solutions will converge at Oðh 3 Þ.Using logarithmic regression, the experimentally observed order of convergence was h 2:77 max .This result indicates that SMM does in fact converge to the S N solution and can converge with high-order accuracy on a problem that is smooth in space and angle, corroborating the claims made in Remark 4.

VI.E. Weak Scaling
We present a weak scaling study of the IP SMM and IP VEF algorithms with p ¼ 2 on the crooked pipe problem from Sec. VI.C. Uniform refinements are used in tandem with increasing the parallel partitioning by a factor of 4 so that the number of degrees of freedom per processor remains fixed.The following results were generated on 29 nodes of the rztopaz machine at LLNL, which has two 18-core Intel Xeon E5-2695 CPUs and 128 Gbytes of memory per node.Timing data are presented as the minimum time measured across three repeated runs.
We use a parallel block Jacobi transport sweep to approximately invert the streaming and collision operator at each iteration.That is, each processor performs a local sweep on its processor-local domain using incoming angular flux information from the previous iteration.This allows each processor to sweep independently from each other at the expense of no longer exactly inverting the streaming and collision operator.In many problems of interest, a parallel block Jacobi-based algorithm can outperform a full parallel sweep-based algorithm since parallel block Jacobi converges quickly when the mean free path is small, the previous time step's solution often provides a good initial guess, and parallel block Jacobi has lower communication costs compared to the full parallel sweep.
Results are presented for SMM without a negative flux fixup applied to the angular flux since convergence was not affected by the lack of a negative flux fixup in the serial crooked pipe problem from Sec. VI.C.However, the zero and scale fixup is still used for comparisons to VEF since a fixup is required for the VEF closures to be well defined.
Table X shows the inner and outer iteration counts in the weak scaling study.Here, the outer iteration scales with the problem size due to the use of the approximate parallel block Jacobi sweep.Compared to VEF, SMM required more iterations to converge, suggesting SMM is more sensitive to the parallel block Jacobi sweep than VEF.This sensitivity could arise from SMM's use of unnormalized closures that depend on both the transport solution's magnitude and angular shape whereas VEF uses more stable, normalized closures that depend only on the angular shape of the solution.
The VEF and SMM moment systems were solved with uniform subspace correction-preconditioned BiCGStab and conjugate gradient, respectively.The previous outer iteration's solution was used as an initial guess for the inner solver.Both methods were scalable in terms of the average number of iterations to converge at each outer iteration.However, VEF shows a clear dependence on the weak scaled problem size in terms of the maximum number of inner iterations required, rising from 18 BiCGStab iterations at the smallest problem size to 95 at the largest.On the other hand, the maximum number of conjugate gradient iterations in the SMM algorithm varied only between 27 and 37.Because of the lagging of angular flux data on parallel boundaries, the closures have nonphysical discontinuities in the early stages of the iteration.For VEF, these nonphysical discontinuities affect the left-hand side of the moment system, degrading the effectiveness of the preconditioned iterative solver.The SMM algorithm avoids this behavior because the nonphysical closure terms are present in the source terms only.
Timing data are provided in Table XI for the average parallel block Jacobi sweep cost, the average cost of forming and solving the moment system, and the total cost.Observe that the average sweep cost is higher for VEF due to the additional costs associated with applying the negative flux fixup within the sweep.As expected, the moment solve cost is lower for SMM compared to VEF We stress that we have used a fixup only to ensure that the moment system's closures are computed using a positive solution as discussed in Remark 6.For moment methods in a multiphysics setting, the intent is to use the moment system's solution to couple to other physics components.Thus, applying the fixup to the angular flux within the sweep increases the cost of the sweep but does not alter the overall simulation's robustness to negativity.Because of this, SMM's ability to converge without a fixup provides a computational advantage over VEF.
Finally, we discuss the weak scaling efficiency of the algorithms.Let Eq. (111) denote the weak scaling efficiency where ideal weak scaling is characterized by ε n ¼ 1: Because of the unavoidable communication costs associated with solving elliptic equations, ideal weak scaling is not expected.In addition, the problem size per processor was chosen based on the computation and storage costs associated with the high-dimensional transport sweep not the lower-dimensional moment solve.This leads to lower parallel efficiency in the moment solve due to the small number of moment unknowns per processor inducing higher relative communication overhead.Figure 13 shows the weak scaling efficiency of the average inner solve cost, the average assembly moment system assembly cost, and the total cost of the algorithm for the IP VEF and IP SMM algorithms.The total solve weak scales poorly due to the dependence on the number of outer iterations on the parallel partitioning from the use of the parallel block Jacobi sweep.Both the VEF and SMM moment solves saturate near 20% efficient while assembly scales at 60% for both methods.

VI.F. Strong Scaling
Finally, we investigate strong scaling for the IP VEF and IP SMM algorithms with p ¼ 2 on the crooked pipe problem from Sec. VI.C.The problem size was fixed at 28 672 equally sized elements with S 12 angular quadrature.This led to 285 048 scalar flux unknowns and 21 676 032 angular flux unknowns.Fixed-point iteration without Anderson acceleration was used.The outer and inner tolerances were 10 À 6 and 10 À 8 , respectively.As with the weak scaling from the previous section, a parallel block Jacobi sweep was to invert the streaming and collision operator, and a negative flux fixup is used for VEF only.Timing data are reported as the minimum time measured across three repeated runs.
Performance is characterized on a single node of the rztopaz machine.Table XII presents the outer and inner iteration data for the strong scaling study.SMM took between 1 and 10 outer, fixed-point iterations more to converge than VEF.As with weak scaling, VEF shows an increase in the maximum number of inner iterations as the parallel partitioning increases, which is not present for the SMM moment solve.Both methods were uniform in average number of inner iterations.
The strong scaling speedup, defined equivalently to ε n in Eq. ( 111), is plotted in Fig. 14.Here, ideal speedup is ε n ¼ n.The average moment system assembly cost scales well for both VEF and SMM due to its low communication overhead with VEF and SMM achieving speedups using 32 processors of 25.2× and 24.2×, respectively.The inner solve requires more communication, and thus, the speedup on 32 processors for the average inner solve cost was reduced to 12× and 13.2× for VEF and SMM, respectively.The total cost is primarily hindered by the scaling of outer iterations with processors.The total speedup on 32 processors was 9.8× for VEF and 11.4× for SMM.

VII. CONCLUSIONS
We have presented a framework for moment methods that encompasses the VEF method and the SMM.Both methods are iterative schemes to solve the Boltzmann transport equation centered around the use of the first two angular moments of the transport equation with suitable closures to accelerate slow-to-converge physics such as scattering.The VEF and SMM algorithms are differentiated only by their choice of closure: VEF uses nonlinear, multiplicative closures while SMM uses linear, additive closures.We demonstrated the close connection between these two algorithms both by way of algebraically reformulating the closures that define each moment method and through linearization of the VEF algorithm about a linearly anisotropic solution.We stress that this linearization process is actually exact since the transport equation itself is linear, and thus, the SMM algorithm can obtain the transport solution even in transport regimes.
The close connection between the VEF and SMM algorithms was used to convert existing independent VEF methods-where the closed moment system is not discretized to be algebraically equivalent to the moments of the discrete transport equation-into corresponding  [10] into discretizations of the SMM moment system.Due to iteratively lagging the additive SMM closures in the SMM algorithm, the discrete SMM moment systems represent standard IP, CG, RT, and HRT discretizations of radiation diffusion with transport-dependent correction sources and can thus directly leverage existing preconditioned iterative solver technology.We use the subspace correction preconditioner from Pazner and Kolev [37] for the IP discretization, AMG for CG and HRT, and AMGbased block preconditioners for RT.The resulting discretizations and solvers are coupled to a high-order DG discretization of the S N transport equations to form robust and efficient independent SMMs for linear, steady-state transport problems.The SMMs were verified to converge the scalar flux with optimal Oðh pþ1 Þ accuracy on refinements of a thirdorder mesh using a manufactured, quadratically anisotropic transport solution.As with the mixed finite element VEF methods from Olivier and Haut, [10] the RT and HRT SMMs exhibited suboptimal convergence for the current.The iterative efficiency of the algorithms was tested in the thick diffusion limit on an orthogonal mesh and a severely distorted, third-order mesh generated with a Lagrangian hydrodynamics code.In both cases, the SMMs converged robustly and performed similarly to the analogous VEF methods.
The methods were also tested on a multimaterial problem designed to emulate the first time step of a thermal radiative transfer calculation.This problem had a 1000× difference in total cross section and sharp material discontinuities that were aligned with the mesh.Using solvers designed for symmetric problems, the SMM moment systems were efficiently solved uniformly with respect to the mesh size and polynomial degree.Using a small Anderson space of size two, the outer fixed-point iteration converged robustly as well.Compared to VEF, the SMMs converged nearly identically with some discretizations and problem sizes converging up to two iterations faster and others no more than three iterations slower.While the SMM moment system was cheaper to assemble and solve than the corresponding VEF moment system, the overall time-tosolution varied by only a relative standard deviation of 10% across the four discretizations and two closures.This invariance in cost is due to the high cost of the transport sweep, which dominates the cost of the operations associated with the moment system.We also demonstrated that the SMMs were able to robustly converge even without a negative flux fixup, an advantage over VEF which requires a negative flux fixup to guarantee the VEF closures are well defined.
We demonstrated that the IP SMM solution converged to the S N transport solution on a thick diffusion limit problem.Using a fixed angular quadrature rule, the SMM and S N solutions were compared as the mesh was refined and seen to converge at the expected, optimal order of convergence.This result indicates that the SMM algorithm does produce high-fidelity, transport solutions.
Finally, we conducted weak and strong scaling studies to demonstrate the scalability of the IP SMM algorithm.We used a parallel block Jacobi transport sweep to invert the streaming and collision operator in parallel.Parallel block Jacobi allows each processor to sweep its local domain independently at the expense of no longer being an exact inversion of the streaming and collision operator.The SMMs generally required more outer iterations than the corresponding VEF method.The IP SMM inner solve was more robust to the parallel partitioning, avoiding the scaling of the maximum number of inner iterations required to converge at each outer iteration exhibited by the IP VEF method.These discrepancies may be due to differences in how the VEF and SMM closures interact with the parallel block Jacobi sweep.Because of SMM's faster moment assembly and moment solve and SMM's ability to run without a negative flux fixup, SMM was able to outperform VEF in total time-tosolution by 8% on the largest weak scaling problem.On 32 processors, SMM attained a strong scaling speedup of 11.4×.
Overall, the SMM iteration was insensitive to the choice of discretization with fixed-point iteration counts roughly independent of the choice of the discretization for the SMM moment system.Furthermore, the convergence rates of the SMM and VEF algorithms were similar with only minor differences in iterations to convergence.Thus, this study indicates that the independent moment algorithm will be robust and efficient regardless of the chosen discretization technique for the moment system or the choice of closure.The discretization and closure then have the flexibility to be chosen to satisfy the requirements of the larger algorithm such as to maximize computational efficiency, improve multiphysics compatibility, or simplify software design and maintenance.In these regards, SMM is preferred since the SMM moment systems were cheaper to assemble and solve and simpler to implement than the corresponding VEF method.

Fig. 1 .
Fig. 1.A depiction of the iterative coupling used in the (a) VEF and (b) SMM algorithms.The transport equation informs the moment system through the closures while the moment system drives the transport equation by computing the scattering physics.By lagging the scattering term, the transport equation can be efficiently inverted.Rapid convergence occurs because the closures are weak functions of the solution.We omit boundary conditions for clarity of presentation.

Fig. 2 .
Fig.2.Depictions of (a) a quadratic, quadrilateral mesh where the mesh control points have been labeled and (b) the element transformation used to describe the left element of the mesh in (a).
DG space Y p .Let Q m;n ð b KÞ denote the space of polynomials of at most m and n in the first and second variables, respectively, such that Q p;p ð b KÞ ¼ Q p ð b KÞ.The RT space uses the local polynomial space Q pþ1;p ð b KÞ � Q p;pþ1 ð b KÞ on each element.In this way, the divergence of any vector v 2 Q pþ1;p ð b KÞ � Q p;pþ1 ð b KÞ belongs to Q p ð b KÞ, the polynomial space used by Y p on each element.Finally, the contravariant Piola transformation, defined in Eq. (51), is used to facilitate the sharing of the degrees of freedom that represent the normal component across an interior face in the mesh.The degree-p RT space is then where denotes the required RT element-local polynomial space composed with the contravariant Piola transformation.The constraint v � n ½ � ½ � ¼ 0 for each interior face holds only when v 2 RT p has a continuous normal component.Thus, since D p ðKÞ � ½H 1 ðKÞ� 2 and each element of RT p has a continuous normal component, RT p � Hðdiv; DÞ.

Fig. 4 .
Fig. 4. A depiction of the distribution of degrees of freedom for the quadratic continuous finite element space.Continuity of members of the finite element space is enforced by sharing degrees of freedom across neighboring elements.

Fig. 5 .
Fig. 5. (a) The distribution of degrees of freedom corresponding to the first degree RT space.Continuity of the normal component is enforced by sharing the degrees of freedom corresponding to the normal component along the interior face between neighboring elements.The circles and squares denote degrees of freedom representing the x and y components, respectively.(b) The distribution of degrees of freedom corresponding to Λ 1 , the space defined as the interior trace of the first degree RT space, on a 3 � 3 mesh.Observe that the interpolation for the normal component in (a) exactly matches the interpolation used in (b) on each face.
[30], Proposition 3.2.2).For the numerical flux, we simply use the average: With this choice of numerical flux, the interior face bilinear form simplifies to where v � n ½ � ½ � ¼ 0 since the normal component is continuous for v 2 RT p .Finally, on boundary of the domain, we solve the SMM boundary conditions in Eq. (71) for the scalar flux and insert this relationship into Pn ¼ Tn þ n=3φ to yield Thus, the RT SMM discretization is: Find ðφ; JÞ 2 Y p � RT p such that Observe that the left-hand side of Eqs.(94) is equivalent to an RT discretization of radiation diffusion with Marshak boundary conditions.

Fig. 7 .
Fig. 7. Plots of the error in the scalar flux as the mesh is refined for (a) linear, (b) quadratic, and (c) cubic finite elements.A quadratically anisotropic MMS transport problem is used.

Fig. 9 .
Fig. 9. Lineouts of the two-dimensional solution at y ¼ 1=2 as � !0 for the (a) IP, (b) CG, (c) RT, and (d) HRT methods on an orthogonal 8 � 8 mesh.The methods all converge to the asymptotic solution indicating they preserve the thick diffusion limit.

Fig.
Fig. The geometry, material data, and boundary conditions for the linearized crooked pipe problem.

Fig. 10 .
Fig. 10.Lineouts of the two-dimensional solution at x ¼ 3:5 as � !0 for the (a) IP, (b) CG, (c) RT, and (d) HRT methods on the triple point mesh.Nonmonotonic oscillations are observed due to imprinting of the highly distorted mesh onto the solution.All methods obtain a nontrivial diffusion solution.

Fig. 12 .
Fig. 12.(a) An example of a mesh built from a tensor product of Chebyshev points in the interval ½0; 1� used to resolve the steep gradients in the solution at the boundary of the domain.(b) A plot of the L 2 ðDÞ norm difference between the solutions generated by the IP SMM method and a DG S N transport method preconditioned with DSA on the thick diffusion limit problem with � ¼ 10 À 2 .Both methods used p ¼ 2. The solutions are compared on four meshes generated with a tensor product of 61, 81, 101, and 121 Chebyshev points in each direction.The errors are presented as a function of the maximum characteristic mesh length h max : When the transport and moment systems are discretized independently, the closure procedure induces errors on the order of the discretization error that can safely be ignored without altering the discrete transport algorithm's spatial or angular order of accuracy (seeRemark 1).Note that the VEF closures are bounded, nonlinear functions of the transport solution due to their normalization by the scalar flux.Furthermore, this normalization means the VEF closures do not depend on the magnitude of the angular fluxonly its angular shape.If ψ ¼ 1 4π ð f þ Ω � gÞ is a linearly anisotropic function for some spatially dependent functions f ðxÞ and gðxÞ, then E ¼ 1=3I and E b ¼ 1=2.
• HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS NUCLEAR SCIENCE AND ENGINEERING • VOLUME 198 • JUNE 2024 pressure and boundary functional.
• HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS NUCLEAR SCIENCE AND ENGINEERING• VOLUME 198 • JUNE 2024 Thus, we can recast the RT VEF first moment [Eq.(83b)] in terms of the unclosed pressure as • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS NUCLEAR SCIENCE AND ENGINEERING • VOLUME 198 • JUNE 2024

TABLE I MMS
Order of Accuracy and Constants NUCLEAR SCIENCE AND ENGINEERING • VOLUME 198 • JUNE 2024

TABLE II The
Number of Fixed-Point Iterations in the Thick Diffusion Limit

TABLE VI The
Average Amount of Time Spent Solving the Moment Systems per Outer Iteration HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1205 thick diffusion limit problem from Sec. VI.B with � ¼ 10 À 2 on an orthogonal mesh.This comparison is facilitated by the following bound on the error.Let the asymptotic spatial error for the SMM and S N solutions in isolation be written as follows:where φ ex is the exact solution; φ SMM and φ SN are the solutions generated by SMM and S N , respectively; C i are the error constants; h is the mesh size; p is the finite element polynomial degree.Using the triangle inequality, observe that Thus, we expect the difference between the solutions generated by SMM and S N to converge with order p þ 1.Note that this behavior is dependent on both solutions being in the asymptotic error regime such that their errors are well characterized by C i h pþ1 .Because of this,

TABLE VIII The
Cost of the Full Algorithm to Solve the Linearized Crooked Pipe Problem

TABLE VII The
Average Amount of Time Spent Assembling the Moment Systems per Outer Iteration NUCLEAR SCIENCE AND ENGINEERING• VOLUME 198 • JUNE 2024

TABLE IX SMM
Iterations to Convergence on the Crooked Pipe with and Without a Fixup

TABLE X Weak
Scaling the Iterations Required by IP VEF and IP SMM Algorithms HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODSwith VEF approximately 1.7× more expensive per iteration.On the largest problem, the total cost of the algorithm was 8% lower for SMM compared to VEF.Note that this is despite SMM requiring 52 more sweeps than VEF on this problem size.This discrepancy is caused by both the faster sweep times from not using a fixup and the faster SMM moment solve.

TABLE XI Timing
Data for the Weak Scaling Study of IP VEF and IP SMM HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS • OLIVIER and HAUT 1209 NUCLEAR SCIENCE AND ENGINEERING • VOLUME 198 • JUNE 2024

TABLE XII Iterations
to Convergence for a Strong Scaling Study of IP VEF and IP SMM Strong scaling speedup as a function of the number of processors on the crooked pipe problem with 28 672 elements, S 12 angular quadrature, and p ¼ 2. The speedup of the average inner solve cost, average moment system assembly cost, and total cost of the algorithm are compared for the IP VEF and IP SMM algorithms.1210 OLIVIER and HAUT • HIGH-ORDER FINITE ELEMENT SECOND MOMENT METHODS discrete SMMs.In particular, we use algebraic manipulations and linearization to systematically convert the IP, continuous finite element (CG), RT mixed finite element, and HRT mixed finite element VEF moment system discretizations from Olivier et al. and Olivier and Haut