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 S_N 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.

Keywords:

• Moment methods
• high-order low-order
• Second Moment Method
• preconditioned iterative solvers

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 high-order low-order methods, developed to solve high-dimensional 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 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 $\mathbf{A}\mathit{x}=\mathit{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 $({\mathrm{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 ${\mathrm{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–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]

Since 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] ${\mathrm{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 ${\mathrm{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 first- or 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)(1a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}+\mathit{q},\mathbf{x}\in \mathcal{D},$(1a) and (1b)(1b) $\mathit{\psi }(\mathbf{x},\mathbf{\Omega })=\bar{\psi }(\mathbf{x},\mathbf{\Omega }),\mathbf{x}\in \mathrm{\partial }\mathcal{D}\mathrm{and}\mathbf{\Omega }\cdot \mathbf{n}<0,$(1b) as our model problem:

(1a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }{ }^{\mathrm{\prime }}+\mathit{q},\mathbf{x}\in \mathcal{D},$(1a)

(1b) $\mathit{\psi }(\mathbf{x},\mathbf{\Omega })=\bar{\psi }(\mathbf{x},\mathbf{\Omega }),\mathbf{x}\in \mathrm{\partial }\mathcal{D}\mathrm{and}\mathbf{\Omega }\cdot \mathbf{n}<0,$(1b)

where $\mathbf{x}\in {\mathbb{R}}^{dim}$ and $\mathbf{\Omega }\in {\mathbb{S}}^2$ are the spatial and angular variables, respectively; $\mathit{\psi }(\mathbf{x},\mathbf{\Omega })$ is the angular flux; $\mathcal{D}\subset {\mathbb{R}}^{dim}$ is the spatial domain of the problem with $\mathrm{\partial }\mathcal{D}$ its boundary; ${\sigma }_t(\mathbf{x})$ and ${\sigma }_s(\mathbf{x})$ are the total and scattering macroscopic cross sections, respectively; $\mathit{q}(\mathbf{x},\mathbf{\Omega })$ is a fixed source; $\bar{\psi }(\mathbf{x},\mathbf{\Omega })$ is the inflow boundary function.

In this section, we derive the moment system and its boundary conditions. Because of the streaming term $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }$, 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:

(2a) $\mathrm{\nabla }\cdot \mathit{J}+{\sigma }_a\mathit{\phi }=Q_0,$(2a)

(2b) $\mathrm{\nabla }\cdot \mathbf{P}+{\sigma }_t\mathit{J}={\mathit{Q}}_1,$(2b)

where ${\sigma }_a(\mathbf{x})\equiv {\sigma }_t(\mathbf{x})-{\sigma }_s(\mathbf{x})$ = absorption macroscopic cross section; $Q_0=\int \mathit{q}\mathrm{d}\mathrm{\Omega }$ and ${\mathit{Q}}_1=\int \mathbf{\Omega }\mathit{q}\mathrm{d}\mathrm{\Omega }$ = zeroth and first moments of the fixed source $\mathit{q}$; $\mathit{\phi }$, $\mathit{J}$, and $\mathbf{P}$ = 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^{\pm }=\mathit{\Omega }\cdot \mathit{n}0\mathit{x\$}\mathit{\Omega }\cdot \mathit{n}\mathit{\psi }\text{d}\mathit{\Omega }$ denote the partial currents where $\mathbf{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

(3) $\begin{array}{rl}& \mathit{J}\cdot \mathbf{n}=J_n^{+}+J_n^{-}\\ & =2J_n^{-}+(J_n^{+}-J_n^{-})\\ & =2J_n^{-}+\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega }.& (3)\end{array}$(3)

Defining

(4) $\mathit{B}(\mathit{\psi })=\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega },$(4)

the boundary conditions for the moment system are

(5) $\mathit{J}\cdot \mathbf{n}=\mathit{B}+2J_{\mathrm{in}},$(5)

where $J_{\mathrm{in}}={\int }_{\mathbf{\Omega }\cdot \mathbf{n}<0}\mathbf{\Omega }\cdot \mathbf{n}\bar{\psi }\mathrm{d}\mathrm{\Omega }$ is the incoming partial current computed using the inflow boundary function $\bar{\psi }$.

Combined, the moment system is given by

(6a) $\mathrm{\nabla }\cdot \mathit{J}+{\sigma }_a\mathit{\phi }=Q_0,\mathbf{x}\in \mathcal{D},$(6a)

(6b) $\mathrm{\nabla }\cdot \mathbf{P}+{\sigma }_t\mathit{J}={\mathit{Q}}_1,\mathbf{x}\in \mathcal{D},$(6b)

(6c) $\mathit{J}\cdot \mathbf{n}=\mathit{B}+2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D}.$(6c)

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 $\mathit{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 $\mathit{\phi }$ (varphi) to denote the moment system’s scalar flux and $\mathit{\varphi }=\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }$ (phi) for the zeroth moment of the angular flux. Since the two solutions differ on the order of the discretization error, we can approximate $\mathit{\phi }/\mathit{\varphi }\approx 1$ and $\mathit{\phi }-\mathit{\varphi }\approx 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:

(7) $\mathbf{P}=\mathbf{E}\mathit{\phi },$(7)

(8) $\mathit{B}=E_b\mathit{\phi },$(8)

where

(9a) $\mathbf{E}=\frac{\int \mathbf{\Omega }\otimes \mathbf{\Omega }\psi \mathrm{d}\mathrm{\Omega }}{\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }},$(9a)

(9b) $E_b=\frac{\int |\mathbf{\Omega }\cdot \mathbf{n}|\psi \mathrm{d}\mathrm{\Omega }}{\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }},$(9b)

are the Eddington tensor and boundary factor, respectively. Since for the continuous equations, $\mathit{\phi }/\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }=1$, these closures are simply algebraic manipulations of the pressure and boundary functional. 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 (see Remark 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 flux—only its angular shape. If $\mathit{\psi }=\frac{1}{4\mathit{\pi }}(\mathit{f}+\mathbf{\Omega }\cdot \mathit{g})$ is a linearly anisotropic function for some spatially dependent functions $\mathit{f}(\mathbf{x})$ and $\mathit{g}(\mathbf{x})$, then $\mathbf{E}=1/3\mathbf{I}$ and $E_b=1/2$.

On the other hand, SMM uses additive closures of the following form:

(10a) $\mathbf{P}=\mathbf{T}+\frac{1}{3}\mathbf{I}\mathit{\phi },$(10a)

(10b) $\mathit{B}=\mathit{\beta }+\frac{1}{2}\mathit{\phi },$(10b)

where

(11a) $\mathbf{T}=\int \mathbf{\Omega }\otimes \mathbf{\Omega }\mathit{\psi }\mathrm{d}\mathrm{\Omega }-\frac{1}{3}\mathbf{I}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega },$(11a)

(11b) $\mathit{\beta }=\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega }-\frac{1}{2}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega },$(11b)

are the SMM correction tensor and boundary factor, respectively. Here, $\mathbf{I}$ denotes the $dim\times 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, $\mathit{\phi }-\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }$ 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 $\mathbf{T}$ and $\mathit{\beta }$ 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)(4) $\mathit{B}(\mathit{\psi })=\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega },$(4) , the VEF and SMM moment systems are given by

(12a) $\mathrm{\nabla }\cdot \mathit{J}+{\sigma }_a\mathit{\phi }=Q_0,\mathbf{x}\in \mathcal{D},$(12a)

(12b) $\mathrm{\nabla }\cdot \left(\mathbf{E}\mathit{\phi }\right)+{\sigma }_t\mathit{J}={\mathit{Q}}_1,\mathbf{x}\in \mathcal{D},$(12b)

(12c) $\mathit{J}\cdot \mathbf{n}=E_b\mathit{\phi }+2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(12c)

and

(13a) $\mathrm{\nabla }\cdot \mathit{J}+{\sigma }_a\mathit{\phi }=Q_0,\mathbf{x}\in \mathcal{D},$(13a)

(13b) $\frac{1}{3}\mathrm{\nabla }\mathit{\phi }+{\sigma }_t\mathit{J}={\mathit{Q}}_1-\mathrm{\nabla }\cdot \mathbf{T},\mathbf{x}\in \mathcal{D},$(13b)

(13c) $\mathit{J}\cdot \mathbf{n}=\frac{1}{2}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(13c)

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

(14) $-\mathrm{\nabla }\cdot \frac{1}{{\sigma }_t}\mathrm{\nabla }\cdot \left(\mathbf{E}\mathit{\phi }\right)+{\sigma }_a\mathit{\phi }=Q_0-\mathrm{\nabla }\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t},$(14)

while for SMM we have

(15) $-\mathrm{\nabla }\cdot \frac{1}{3{\sigma }_t}\mathrm{\nabla }\mathit{\phi }+{\sigma }_a\mathit{\phi }=Q_0-\mathrm{\nabla }\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}+\mathrm{\nabla }\cdot \frac{1}{{\sigma }_t}\mathrm{\nabla }\cdot \mathbf{T}.$(15)

Here, boundary conditions are specified by Eqs. (12c)(12c) $\mathit{J}\cdot \mathbf{n}=E_b\mathit{\phi }+2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(12c) and (13c)(13c) $\mathit{J}\cdot \mathbf{n}=\frac{1}{2}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(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 non-self-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, $\mathrm{\nabla }\cdot \frac{1}{{\sigma }_t}\mathrm{\nabla }\cdot \mathbf{T},$ shares the same mathematical structure as the differential term in the VEF operator, $\mathrm{\nabla }\cdot \frac{1}{{\sigma }_t}\mathrm{\nabla }\cdot \left(\mathbf{E}\mathit{\phi }\right)$, and thus will have analogous discretization challenges. In this paper, we leverage recently developed discretization techniques for the VEF moment system to derive discretizations for the SMM moment system.

II.C. The Moment Method Algorithm

Moment methods solve the coupled transport-moment system simultaneously. The transport equation is used to provide the VEF or SMM closures while the moment system is used to compute the moment-dependent 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 $\mathcal{M}(\mathit{\psi },\mathbf{X})=0$ denote one of the moment systems derived in Sec. II.B with $\mathbf{X}$ the moment system’s unknowns. For example, $\mathcal{M}(\mathit{\psi },\mathbf{X})$ could represent the VEF moment system in first-order form given by Eqs. (12)(8) $\mathit{B}=E_b\mathit{\phi },$(8) where $\mathbf{X}$ would include both the scalar flux and the current. In the case of the second-order form, we would set $\mathbf{X}=\mathit{\phi }$ since the scalar flux is the only unknown. For VEF, $\mathcal{M}(\mathit{\psi },\mathbf{X})$ is nonlinear in $\mathit{\psi }$ and linear in the moments $\mathbf{X}$. For SMM, $\mathcal{M}(\mathit{\psi },\mathbf{X})$ is linear in both arguments.

The moment algorithm solves the coupled system given by

(16a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q},$(16a)

(16b) $\mathcal{M}(\mathit{\psi },\mathbf{X})=0,$(16b)

where transport boundary conditions are specified in Eq. (1b)(1b) $\mathit{\psi }(\mathbf{x},\mathbf{\Omega })=\bar{\psi }(\mathbf{x},\mathbf{\Omega }),\mathbf{x}\in \mathrm{\partial }\mathcal{D}\mathrm{and}\mathbf{\Omega }\cdot \mathbf{n}<0,$(1b) . The moment system’s boundary conditions are given by Eq. (12c)(12c) $\mathit{J}\cdot \mathbf{n}=E_b\mathit{\phi }+2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(12c) for a VEF method and Eq. (13c)(13c) $\mathit{J}\cdot \mathbf{n}=\frac{1}{2}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(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)(16a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q},$(16a) 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:

(17) $\mathbf{L}\mathit{\psi }=\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi },$(17)

represent the streaming and collision operator. The coupled transport-moment system can then be rewritten as

(18a) $\mathbf{L}\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q},$(18a)

(18b) $\mathcal{M}(\mathit{\psi },\mathbf{X})=0.$(18b)

By linearly eliminating the angular flux, the coupled system is equivalent to

(19) $\mathcal{M}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right),\mathbf{X}\right)=0.$(19)

Observe that Eq. (19)(19) $\mathcal{M}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right),\mathbf{X}\right)=0.$(19) is now a function of the moment solution only. That is, we can define

(20) $\mathbf{F}(\mathbf{X})=\mathcal{M}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right),\mathbf{X}\right)$(20)

and equivalently solve $\mathbf{F}(\mathbf{X})=0$. In this reduced problem, the angular flux appears only as an auxiliary variable used to compute the residual $\mathbf{F}(\mathbf{X}),$ and we say that the angular flux is determined by the moment system. This reduced formulation $\mathbf{F}(\mathbf{X})=0$ has much lower dimension than the original coupled system given in Eq. (16)(16a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q},$(16a) but has the same solution. Because of this, advanced solvers for $\mathbf{F}(\mathbf{X})$ can be applied that would otherwise be impractical for Eq. (16)(16a) $\mathbf{\Omega }\cdot \mathrm{\nabla }\mathit{\psi }+{\sigma }_t\mathit{\psi }=\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q},$(16a) 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

(21) $\mathbf{V}(\mathit{\psi })\mathbf{X}=\mathit{f}$(21)

represent the VEF moment system such that $\mathcal{M}(\mathit{\psi },\mathbf{X})=\mathbf{V}(\mathit{\psi })\mathbf{X}-\mathit{f}$. We then have that

(22) $\mathbf{F}(\mathbf{X})=\mathbf{V}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right)\right)\mathbf{X}-\mathit{f}=0.$(22)

Operating by the inverse of the VEF moment system, the coupled transport-VEF system is equivalent to

(23) $\mathbf{X}=\mathbf{V}{\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right)\right)}^{-1}\mathit{f}.$(23)

For SMM, the moment system is of the form

(24) $\mathbf{DX}=\mathit{b}(\mathit{\psi }),$(24)

where $\mathbf{D}$ is a diffusion operator and $\mathit{b}(\mathit{\psi })$ includes the moments of the fixed-source and the transport-dependent correction sources. The root-finding problem $\mathbf{F}(\mathbf{X})=0$ is then equivalent to

(25) $\mathbf{X}={\mathbf{D}}^{-1}\mathit{b}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right)\right).$(25)

Thus, for both VEF and SMM, the solution of the coupled transport-moment system is the fixed point:

(26) $\mathbf{X}=\mathbf{G}(\mathbf{X}),$(26)

where $\mathbf{G}(\mathbf{X})$ is given by

(27) $\mathbf{G}(\mathbf{X})=\mathbf{V}{\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right)\right)}^{-1}\mathit{f}$(27)

for VEF and

(28) $\mathbf{G}(\mathbf{X})={\mathbf{D}}^{-1}\mathit{b}\left({\mathbf{L}}^{-1}\left(\frac{{\sigma }_s}{4\mathit{\pi }}\mathit{\phi }+\mathit{q}\right)\right)$(28)

for SMM. The fixed-point operator $\mathbf{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 $\mathbf{G}$ show the key differences between the VEF and SMM algorithms. VEF has a transport-dependent left-hand-side operator while the right-hand-side 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 $\mathbf{X}=\mathbf{G}(\mathbf{X})$ is fixed-point iteration:

(29) ${\mathbf{X}}^{\mathit{k}+1}=\mathbf{G}({\mathbf{X}}^k),$(29)

where ${\mathbf{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 $\mathit{k}$ previous evaluations of the fixed-point operator and chooses the next iterate by finding the optimal linear combination of the $\mathit{k}$ operator evaluations that minimizes the residual $\mathbf{F}(\mathbf{X})$. For a storage cost of $\mathit{k}$ moment system-sized vectors, Anderson acceleration increases iterative efficiency and robustness. We stress that the fixed-point problem in Eq. (29)(29) ${\mathbf{X}}^{\mathit{k}+1}=\mathbf{G}({\mathbf{X}}^k),$(29) corresponds to the reduced formulation where the angular flux is determined by the moment system, avoiding the need to store $\mathit{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]

Remark 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.

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.

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 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 $\mathit{f}$ evaluated at some function $\mathit{u}$ in the direction of $\mathit{v}$ is defined as

(30) $\mathit{D}[\mathit{f}](\mathit{u},\mathit{v})=\underset{\mathit{\omega }\to 0}{lim}\frac{\mathit{f}(\mathit{u}+\mathit{\omega v})-\mathit{f}(\mathit{u})}{\omega }.$(30)

Assuming continuity of $\mathit{f}$, the Gateaux derivative is equivalent to

(31) $\mathit{D}[\mathit{f}](\mathit{u},\mathit{v})={\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\omega }}\mathit{f}(\mathit{u}+\mathit{\omega v})\right]}_{\mathit{\omega }=0}.$(31)

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 $\mathbf{u}={\left[\begin{array}{ccc}{u}_1& \dots & u_n\end{array}\right]}^{\mathit{T}}$ and $\mathbf{v}={\left[\begin{array}{ccc}{v}_1& \dots & v_n\end{array}\right]}^{\mathit{T}}$ for some $\mathit{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 $\mathit{f}(\mathit{u})$ expanded about $u_0$ is

(32) $\mathit{f}(\mathit{u})\stackrel{\mathrm{TSE}}{⟶}\mathit{f}(u_0)+\mathit{D}[\mathit{f}](u_0,\mathit{u}-u_0).$(32)

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 ${\psi }_0$ be such a linearly anisotropic approximation to the transport problem at hand. Using Eq. (31)(31) $\mathit{D}[\mathit{f}](\mathit{u},\mathit{v})={\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\omega }}\mathit{f}(\mathit{u}+\mathit{\omega v})\right]}_{\mathit{\omega }=0}.$(31) , the Gateaux derivative of the Eddington tensor at ${\psi }_0$ in the direction $\mathit{\psi }$ is

(33) $\begin{array}{rl}& \mathit{D}[\mathbf{E}]({\psi }_0,\mathit{\psi })={\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\omega }}\mathbf{E}({\psi }_0+\mathit{\omega \psi })\right]}_{\mathit{\omega }=0}\\ & ={\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\omega }}\frac{{\mathbf{P}}_0+\mathit{\omega }\mathbf{P}}{{\varphi }_0+\mathit{\omega \varphi }}\right]}_{\mathit{\omega }=0}\\ & ={\left.\frac{\mathbf{P}({\varphi }_0+\mathit{\omega \varphi })-({\mathbf{P}}_0+\mathit{\omega }\mathbf{P})\mathit{\varphi }}{{({\varphi }_0+\mathit{\omega \varphi })}^2}\right|}_{\mathit{\omega }=0}\\ & =\frac{\mathbf{P}{\varphi }_0-{\mathbf{P}}_0\mathit{\varphi }}{{\varphi }_0^2}\\ & =\frac{1}{{\varphi }_0}\left(\mathbf{P}-{\mathbf{E}}_0\mathit{\varphi }\right),\end{array}$(33)

where we have used the quotient rule and have set ${\mathbf{P}}_0=\int \mathbf{\Omega }\otimes \mathbf{\Omega }{\psi }_0\mathrm{d}\mathrm{\Omega }$, ${\varphi }_0=\int {\psi }_0\mathrm{d}\mathrm{\Omega }$, and ${\mathbf{E}}_0={\mathbf{P}}_0/{\varphi }_0$. Since ${\psi }_0$ is linearly anisotropic, ${\mathbf{E}}_0=1/3\mathbf{I}$. Thus, the Gateaux derivative of the Eddington tensor is equivalent to

(34) $\mathit{D}[\mathbf{E}]({\psi }_0,\mathit{\psi })=\frac{1}{{\varphi }_0}\left(\int \mathbf{\Omega }\otimes \mathbf{\Omega }\mathit{\psi }\mathrm{d}\mathrm{\Omega }-\frac{1}{3}\mathbf{I}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }\right)\equiv \frac{1}{{\varphi }_0}\mathbf{T}(\mathit{\psi }),$(34)

where $\mathbf{T}(\mathit{\psi })$ is the SMM correction tensor defined in Eq. (11a)(11a) $\mathbf{T}=\int \mathbf{\Omega }\otimes \mathbf{\Omega }\mathit{\psi }\mathrm{d}\mathrm{\Omega }-\frac{1}{3}\mathbf{I}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega },$(11a) . By an analogous calculation, the Gateaux derivative of the Eddington boundary factor is

(35) $\mathit{D}[E_b]({\psi }_0,\mathit{\psi })=\frac{1}{{\varphi }_0}\left(\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega }-E_{\mathit{b}0}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega }\right),$(35)

where

(36) $E_{\mathit{b}0}=\frac{\int |\mathbf{\Omega }\cdot \mathbf{n}|{\psi }_0\mathrm{d}\mathrm{\Omega }}{\int {\psi }_0\mathrm{d}\mathrm{\Omega }}=\frac{\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathrm{d}\mathrm{\Omega }}{4\mathit{\pi }}=\frac{1}{2},$(36)

given that ${\psi }_0$ is linearly anisotropic. Thus, the Gateaux derivative of the Eddington boundary factor and the SMM boundary correction factor are related by

(37) $\mathit{D}[E_b]({\psi }_0,\mathit{\psi })=\frac{1}{{\varphi }_0}\mathit{\beta }(\mathit{\psi }),$(37)

where $\mathit{\beta }(\mathit{\psi })$ is the SMM boundary factor defined in Eq. (11b)(11b) $\mathit{\beta }=\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathit{\psi }\mathrm{d}\mathrm{\Omega }-\frac{1}{2}\int \mathit{\psi }\mathrm{d}\mathrm{\Omega },$(11b) .

We now linearize the pressure and boundary functional using the above functional derivatives of the VEF closures. Let $\mathbf{y}={\left[\begin{array}{cc}\mathit{\psi }& \mathit{\phi }\end{array}\right]}^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 ${\mathbf{y}}_0={\left[\begin{array}{cc}{\psi }_0& {\phi }_0\end{array}\right]}^{\mathit{T}}$ to be the linearly anisotropic approximation of the coupled system. Here, we consider the moment solution $\mathit{\phi }$ to be an independent variable. Using the Gateaux derivative, we can write the first-order Taylor series expansion of the pressure as

(38) $\begin{array}{rl}& \mathbf{P}(\mathbf{y})\stackrel{\mathrm{TSE}}{⟶}{\left[\mathbf{E}\phi \right]}_{\mathbf{y}={\mathbf{y}}_0}+\mathit{D}[\mathbf{E}\mathit{\phi }]({\mathbf{y}}_0,\mathbf{y}-{\mathbf{y}}_0)\\ & ={\mathbf{E}}_0{\phi }_0+\mathit{D}[\mathbf{E}{\phi }_0]({\psi }_0,\mathit{\psi }-{\psi }_0)+\mathit{D}[{\mathbf{E}}_0\mathit{\phi }]({\phi }_0,\mathit{\phi }-{\phi }_0)\\ & ={\mathbf{E}}_0{\phi }_0+\frac{{\phi }_0}{{\varphi }_0}\mathbf{T}(\mathit{\psi }-{\psi }_0)+{\mathbf{E}}_0(\mathit{\phi }-{\phi }_0)\\ & ={\mathbf{E}}_0\mathit{\phi }+\frac{{\phi }_0}{{\varphi }_0}\mathbf{T}(\mathit{\psi }-{\psi }_0)\\ & ={\mathbf{E}}_0\mathit{\phi }+\mathbf{T}(\mathit{\psi }),& (38)\end{array}$(38)

where $\mathbf{T}(\mathit{\psi }-{\psi }_0)=\mathbf{T}(\mathit{\psi })-\mathbf{T}({\psi }_0)=\mathbf{T}(\mathit{\psi })$ since the SMM correction tensor is linear and zero when the argument is linearly anisotropic, respectively. Further, we have simplified ${\phi }_0/{\varphi }_0\approx 1$ by ignoring terms on the order of the discretization error (see Remark 1). Using ${\mathbf{E}}_0=1/3\mathbf{I}$, observe that the extreme left-hand sides and right-hand sides of the linearized pressure exactly match the SMM closure given in Eq. (10a)(10a) $\mathbf{P}=\mathbf{T}+\frac{1}{3}\mathbf{I}\mathit{\phi },$(10a) .

Analogously, the first-order Taylor series expansion of the boundary functional is

(39) $\frac{\mathit{B}(\mathit{y})\text{TSE}}{�{\left[E_b\mathit{\phi }\right]}_{\mathit{y}=y_0}+\mathit{D}[E_b\mathit{\phi }](y_0,\mathit{y}-y_0)}=E_{\mathit{b}0}\mathit{\phi }+\frac{{\phi }_0}{{\varphi }_0}\mathit{\beta }(\mathit{\psi }-{\psi }_0)=E_{\mathit{b}0}\mathit{\phi }+\mathit{\beta }(\mathit{\psi }),(39)$(39)

where we have applied the same arguments to simplify $\mathit{\beta }(\mathit{\psi }-{\psi }_0)$ and ${\phi }_0/{\varphi }_0$ as was used for the pressure. Recognizing that $E_{\mathit{b}0}=1/2$, the linearized boundary functional is equivalent to the SMM closure in Eq. (10b)(10b) $\mathit{B}=\mathit{\beta }+\frac{1}{2}\mathit{\phi },$(10b) .

Remark 3: While angular quadrature rules can exactly evaluate ${\mathbf{E}}_0=1/3\mathbf{I}$, it is not possible to exactly integrate $E_{\mathit{b}0}$ with numerical quadrature due to the numerator’s integrand $|\mathbf{\Omega }\cdot \mathbf{n}|{\psi }_0$ having a discontinuous derivative. Generally, $\mathit{\beta }$, $J_{\mathit{in}}$, and $E_{\mathit{b}0}$ 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 $\mathbf{P}=\mathbf{E}\mathit{\phi }$ and $\mathit{B}=E_b\mathit{\phi }$ replaced with $\mathbf{T}+1/3\mathbf{I}\mathit{\phi }$ and $\mathit{\beta }+1/2\mathit{\phi }$, 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 $\mathcal{D}\subset {\mathbb{R}}^2$ be the domain of the problem tessellated into a collection of elements $\mathcal{T}$ such that

(40) $\mathcal{D}=\bigcup _{\mathit{K}\in \mathcal{T}}\mathit{K},$(40)

where $\mathit{K}\in \mathcal{T}$ is a quadrilateral element in the mesh obtained as the image of the reference square $\hat{K}=[0,1{]}^2$ under an invertible, polynomial mapping $\mathbf{T}:\hat{K}\to \mathit{K}$. Here, the element transformation belongs to $[{\mathcal{Q}}_m(\hat{K}){]}^2,$ where ${\mathcal{Q}}_m(\hat{K})$ is the space of polynomials of degree at most $\mathit{m}$ in each variable. We use a nodal basis to describe the element transformation. Let ${\xi }_i$ denote the $\mathit{m}+1$ Gauss-Lobatto points in the interval $[0,1]$. The $(\mathit{m}+1{)}^2$ points ${\xi }_i$ on the unit square $\hat{K}=[0,1{]}^2$ are given by the twofold Cartesian product of the one-dimensional points. Further, let ${\mathrm{\ell }}_i$ denote the Lagrange interpolating polynomial that satisfies ${\mathrm{\ell }}_i({\xi }_j)={\delta }_{\mathit{ij}},$ where ${\delta }_{\mathit{ij}}$ is the Kronecker delta. The set of functions $\{{\mathrm{\ell }}_i{\}}_{\mathit{i}=1}^{{(\mathit{m}+1)}^2}$ forms a basis for ${\mathcal{Q}}_m(\hat{K})$. For each element $\mathit{K}\in \mathcal{T}$, the mapping is of the form

(41) $\mathbf{x}(\mathit{\xi })=\mathbf{T}(\mathit{\xi })=\sum _{\mathit{i}=1}^{{(\mathit{m}+1)}^2}{\mathbf{x}}_{\mathit{K},\mathit{i}}{\mathrm{\ell }}_i(\mathit{\xi }),$(41)

where $\mathbf{x}\in \mathit{K}$, $\mathit{\xi }\in \hat{K}$, and $\{{\mathbf{x}}_{\mathit{K},\mathit{i}}{\}}_{\mathit{i}=1}^{{(\mathit{m}+1)}^2}$ = control points for element $\mathit{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.

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).

We define the set of unique faces in the mesh as $\mathrm{\Gamma }$ with ${\mathrm{\Gamma }}_0=\mathrm{\Gamma \setminus \partial }\mathcal{D}$ the set of interior faces in the mesh and ${\mathrm{\Gamma }}_b=\mathrm{\Gamma }\cap \mathrm{\partial }\mathcal{D}$ the set of boundary faces so that $\mathrm{\Gamma }={\mathrm{\Gamma }}_0\cup {\mathrm{\Gamma }}_b$. We will often work with piecewise functions defined element-wise on the mesh. On an interior face $\mathcal{F}\in {\mathrm{\Gamma }}_0$ between two arbitrary elements $K_1$ and $K_2$, we use the convention that $\mathbf{n}$ is a unit vector perpendicular to the shared face $K_1\cap K_2$ pointing from $K_1$ to $K_2$. For a piecewise polynomial $\mathit{u}$ described by $u_1$ for $\mathbf{x}\in K_1$ and $u_2$ for $\mathbf{x}\in K_2$, the average, $\{\{\cdot \}\}$, and jump, $\left[\left[\cdot \right]\right]$, of $\mathit{u}$ across the face $\mathcal{F}=K_1\cap K_2$ are defined as

(42) $\begin{array}{rl}& \{\{\mathit{u}\}\}=\frac{1}{2}(u_1+u_2),\left[\left[\mathit{u}\right]\right]=u_1-u_2,\\ & \mathrm{on}\mathcal{F}=K_1\cap K_2\in {\mathrm{\Gamma }}_0,& (42)\end{array}$(42)

with an analogous definition for vector-valued piecewise functions. On boundary faces, the jump and average are defined as

(43) $\left[\left[\mathit{u}\right]\right]=\mathit{u},\{\{\mathit{u}\}\}=\mathit{u},\mathrm{on}\mathcal{F}\in {\mathrm{\Gamma }}_b,$(43)

and likewise for vectors on the boundary. We note that a piecewise-continuous function $\mathit{u}$ satisfies $\left[\left[\mathit{u}\right]\right]=0$ and $\{\{\mathit{u}\}\}=\mathit{u}$ on each interior face.

Finally, we define the “broken” gradient, denoted by ${\mathrm{\nabla }}_h$, obtained by applying the gradient locally on each element. That is,

(44) $({\mathrm{\nabla }}_h\mathit{u}){|}_K=\mathrm{\nabla }(\mathit{u}{|}_K),\mathrm{\forall K}\in \mathcal{T}.$(44)

The broken gradient is important for piecewise-discontinuous functions where the global gradient is not well defined due to the discontinuities present at interior mesh interfaces.

III.B. Integration Transformations

The mesh transformations $\mathbf{T}$ are used to facilitate numerical integration on arbitrary elements. Let $\mathit{\xi }={\left[\begin{array}{cc}\mathit{\xi }& \mathit{\eta }\end{array}\right]}^T\in \hat{K}$ denote the reference coordinates and $\mathbf{x}={\left[\begin{array}{cc}\mathit{x}& \mathit{y}\end{array}\right]}^T=\mathbf{T}(\mathit{\xi })$ denote the physical coordinates. The Jacobian of the transformation is

(45) $\mathbf{F}=\frac{\mathrm{\partial }\mathbf{T}}{\mathrm{\partial }\mathit{\xi }}=\left[\begin{array}{cc}\frac{\mathrm{\partial x}}{\mathrm{\partial \xi }}& \frac{\mathrm{\partial x}}{\mathrm{\partial \eta }}\\ \frac{\mathrm{\partial y}}{\mathrm{\partial \xi }}& \frac{\mathrm{\partial y}}{\mathrm{\partial \eta }}\end{array}\right],$(45)

with $\mathit{J}=|\mathbf{F}|$ its determinant. The partial derivatives of the mesh transformation are computed using the basis expansion of the mesh transformations. In other words,

(46) $\mathbf{F}=\sum _{\mathit{i}=1}^{{(\mathit{m}+1)}^2}{\mathbf{x}}_{\mathit{K},\mathit{i}}\otimes \hat{\mathrm{\nabla }}{\mathrm{\ell }}_i=\sum _{\mathit{i}=1}^{{(\mathit{m}+1)}^2}\left[\begin{array}{cc}{x}_{\mathit{K},\mathit{i}}\frac{\mathrm{\partial }{\mathrm{\ell }}_i}{\mathrm{\partial \xi }}& x_{\mathit{K},\mathit{i}}\frac{\mathrm{\partial }{\mathrm{\ell }}_i}{\mathrm{\partial \eta }}\\ \mathit{\hspace{0.17em}}\\ y_{\mathit{K},\mathit{i}}\frac{\mathrm{\partial }{\mathrm{\ell }}_i}{\mathrm{\partial \xi }}& y_{\mathit{K},\mathit{i}}\frac{\mathrm{\partial }{\mathrm{\ell }}_i}{\mathrm{\partial \eta }}\end{array}\right],$(46)

where ${\mathbf{x}}_{\mathit{K},\mathit{i}}={\left[\begin{array}{cc}{\mathit{x}}_{\mathit{K},\mathit{i}}& y_{\mathit{K},\mathit{i}}\end{array}\right]}^{\mathit{T}};$ $\hat{\mathrm{\nabla }}$ = gradient with respect to the reference coordinates $\mathit{\xi }$. In this paper, integration over the domain is implicitly computed in reference space using the following sum:

(47) ${\int }_{\mathcal{D}}\left(\cdot \right)\mathrm{d}\mathbf{x}=\sum _{\mathit{K}\in \mathcal{T}}{\int }_K\left(\cdot \right)\mathrm{d}\mathbf{x}=\sum _{\mathit{K}\in \mathcal{T}}{\int }_{\hat{K}}\left(\cdot \right)\mathit{J}\mathrm{d}\mathit{\xi }.$(47)

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 $\hat{K}$.

We now discuss the transformations used to represent the integrand of Eq. (47)(47) ${\int }_{\mathcal{D}}\left(\cdot \right)\mathrm{d}\mathbf{x}=\sum _{\mathit{K}\in \mathcal{T}}{\int }_K\left(\cdot \right)\mathrm{d}\mathbf{x}=\sum _{\mathit{K}\in \mathcal{T}}{\int }_{\hat{K}}\left(\cdot \right)\mathit{J}\mathrm{d}\mathit{\xi }.$(47) in reference space. For a scalar function $\mathit{u}:\mathcal{D}\to \mathbb{R}$, denote by $\hat{u}:\hat{K}\to \mathbb{R}$ its representation in reference space. The functions $\mathit{u}$ and $\hat{u}$ are related by

(48) $\mathit{u}(\mathbf{x})=\hat{u}({\mathbf{T}}^{-1}(\mathbf{x})),$(48)

where ${\mathbf{T}}^{-1}:\mathit{K}\to \hat{K}$ denotes the inverse mesh transformation. Integration over the physical element is then equivalent to

(49) ${\int }_K\mathit{u}\mathrm{d}\mathbf{x}={\int }_{\hat{K}}\hat{u}\mathit{J}\mathrm{d}\mathit{\xi }.$(49)

Using the chain rule, the gradient transforms as

(50) $\mathrm{\nabla }\hat{u}=\left[\begin{array}{c}\frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial }\mathit{\xi }}\frac{\mathrm{\partial \xi }}{\mathrm{\partial x}}+\frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial \eta }}\frac{\mathrm{\partial \eta }}{\mathrm{\partial x}}\\ \mathit{\hspace{0.17em}}\\ \frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial \xi }}\frac{\mathrm{\partial \xi }}{\mathrm{\partial y}}+\frac{\mathrm{\partial }\hat{u}}{\mathrm{\partial \eta }}\frac{\mathrm{\partial \eta }}{\mathrm{\partial y}}\end{array}\right]={\mathbf{F}}^{-\mathit{T}}\hat{\mathrm{\nabla }}\hat{u}.$(50)

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]:

(51) $\mathit{v}=\frac{1}{J}\mathbf{F}\stackrel{\mathbf{ˆ}}{\mathit{v}}\circ {\mathbf{T}}^{-1},$(51)

where $\mathit{v}:\mathcal{D}\to {\mathbb{R}}^2$ and $\stackrel{\mathbf{ˆ}}{\mathit{v}}:\hat{K}\to {\mathbb{R}}^2$ = representations of a vector in physical and reference space, respectively. Members of the RT space have a continuous normal component, a requirement of any $\mathit{H}(\mathrm{div};\mathcal{D})$ - conforming space.^[30] The contravariant Piola transformation represents vectors on the so-called tangent space spanned by the columns of the Jacobian matrix $\mathbf{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 $\mathit{H}(\mathrm{div};\mathcal{D})$.

For a contravariant vector,

(52) ${\int }_K\mathrm{\nabla }\mathit{u}\cdot \mathit{v}\mathrm{d}\mathbf{x}={\int }_{\hat{K}}{\mathbf{F}}^{-\mathit{T}}\hat{\mathrm{\nabla }}\hat{u}\cdot \frac{1}{J}\mathbf{F}\hat{v}\mathit{J}\mathrm{d}\mathit{\xi }={\int }_{\hat{K}}\hat{\mathrm{\nabla }}\hat{u}\cdot \stackrel{\mathbf{ˆ}}{\mathit{v}}\mathrm{d}\mathit{\xi }.$(52)

The gradient transforms as

(53) $\mathrm{\nabla }\mathit{v}=\mathrm{\nabla }\left(\frac{1}{J}\mathbf{F}\hat{v}\circ {\mathbf{T}}^{-1}\right)=\frac{1}{J}\mathbf{F}\left(\hat{\mathrm{\nabla }}\stackrel{\mathbf{ˆ}}{\mathit{v}}-\stackrel{\mathbf{ˆ}}{\mathbf{B}}\right){\mathbf{F}}^{-1},$(53)

where

(54) $\stackrel{\mathbf{ˆ}}{\mathbf{B}}=\frac{1}{J}\hat{\mathrm{\nabla }}\left(\mathit{J}{\mathbf{F}}^{-1}\right)\mathbf{F}\stackrel{\mathbf{ˆ}}{\mathit{v}}$(54)

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 Ref. [10], Appendix A, where implementation details are also provided. We note that $\stackrel{\mathbf{ˆ}}{\mathbf{B}}$ is related to the Hessian of the mesh transformation and is thus zero when the mesh transformation is affine such that $\mathbf{T}(\mathit{\xi })=\mathbf{A}\mathit{\xi }+\mathit{b}$ for some $\mathbf{A}\in {\mathbb{R}}^{dim\times dim}$ and $\mathit{b}\in {\mathbb{R}}^{dim}$ independent of $\mathit{\xi }$. In addition, $\text{trace}(\widehat{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

(55) $\mathrm{\nabla }\cdot \mathit{v}=\mathrm{trace}\left(\mathrm{\nabla }\mathit{v}\right)=\frac{1}{J}\mathrm{trace}\left(\mathbf{F}\left(\hat{\mathrm{\nabla }}\stackrel{\mathbf{ˆ}}{\mathit{v}}-\stackrel{\mathbf{ˆ}}{\mathbf{B}}\right){\mathbf{F}}^{-1}\right)=\frac{1}{J}\left(\mathrm{trace}(\hat{\mathrm{\nabla }}\stackrel{\mathbf{ˆ}}{\mathit{v}})-\mathrm{trace}(\stackrel{\mathbf{ˆ}}{\mathbf{B}})\right)=\frac{1}{J}\hat{\mathrm{\nabla }}\cdot \stackrel{\mathbf{ˆ}}{\mathit{v}}.$(55)

Thus,

(56) ${\int }_K\mathit{u\nabla }\cdot \mathit{v}\mathrm{d}\mathbf{x}={\int }_{\hat{K}}\hat{u}\hat{\mathrm{\nabla }}\cdot \stackrel{\mathbf{ˆ}}{\mathit{v}}\mathrm{d}\mathit{\xi }.$(56)

Combining the results from Eqs. (52)(52) ${\int }_K\mathrm{\nabla }\mathit{u}\cdot \mathit{v}\mathrm{d}\mathbf{x}={\int }_{\hat{K}}{\mathbf{F}}^{-\mathit{T}}\hat{\mathrm{\nabla }}\hat{u}\cdot \frac{1}{J}\mathbf{F}\hat{v}\mathit{J}\mathrm{d}\mathit{\xi }={\int }_{\hat{K}}\hat{\mathrm{\nabla }}\hat{u}\cdot \stackrel{\mathbf{ˆ}}{\mathit{v}}\mathrm{d}\mathit{\xi }.$(52) and (56(56) ${\int }_K\mathit{u\nabla }\cdot \mathit{v}\mathrm{d}\mathbf{x}={\int }_{\hat{K}}\hat{u}\hat{\mathrm{\nabla }}\cdot \stackrel{\mathbf{ˆ}}{\mathit{v}}\mathrm{d}\mathit{\xi }.$(56) ) yields

(57) ${\int }_{\mathrm{\partial K}}\mathit{u}\mathit{v}\cdot \mathbf{n}\mathrm{d}\mathit{s}={\int }_{\mathrm{\partial }\hat{K}}\hat{u}\stackrel{\mathbf{ˆ}}{\mathit{v}}\cdot \stackrel{\mathbf{ˆ}}{\mathbf{n}}\mathrm{d}\hat{s},$(57)

where $\stackrel{\mathbf{ˆ}}{\mathbf{n}}$ is the normal vector in reference space corresponding to the physical space normal $\mathbf{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 one-dimensional 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 $\mathcal{T}$ or the interior skeleton of the mesh ${\mathrm{\Gamma }}_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(\mathcal{D})$, $H^1(\mathcal{D})$, and $\mathit{H}(\mathrm{div};\mathcal{D})$.

III.C.1. Discontinuous Galerkin

Let ${\mathbb{Q}}_p(\mathit{K})$ denote the space of polynomials of at most degree $\mathit{p}$ mapped from the reference element. That is,

(58) ${\mathbb{Q}}_p(\mathit{K})=\{\mathit{u}=\hat{u}\circ {\mathbf{T}}^{-1}:\hat{u}\in {\mathcal{Q}}_p(\hat{K})\},$(58)

where $\hat{u}$ indicates a function defined on the reference element. The delineation between $\mathcal{Q}$ and $\mathbb{Q}$ is required when non-affine mesh transformations are used. In such case, $\mathit{u}=\hat{u}\circ {\mathbf{T}}^{-1}\notin {\mathcal{Q}}_p(\mathit{K})$ even if $\hat{u}\in {\mathcal{Q}}_p(\hat{K})$. 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(\mathcal{D})$ space defined by

(59) $L^2(\mathcal{D})=\{\mathit{u}:\int u^2\mathrm{d}\mathbf{x}<\mathrm{\infty }\},$(59)

the space of square-integrable functions. We denote the degree-$\mathit{p}$ DG space with

(60) $Y_p=\{\mathit{u}\in L^2(\mathcal{D}):\mathit{u}{|}_K\in {\mathbb{Q}}_p(\mathit{K}),\mathrm{\forall K}\in \mathcal{T}\},$(60)

so that each function $\mathit{u}\in 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\times 3$ mesh for the space $Y_1$. The Gauss-Legendre interpolation points are used to define the nodal basis for the local polynomial space to highlight that degrees of freedom are not shared between elements.

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.

III.C.2. Continuous Finite Element

The continuous finite element (CG) space is a finite-dimensional subspace of the $H^1(\mathcal{D})$ Sobolev space defined as

(61) $H^1(\mathcal{D})=\{\mathit{u}\in L^2(\mathcal{D}):\int \mathrm{\nabla }\mathit{u}\cdot \mathrm{\nabla u}\mathrm{d}\mathbf{x}<\mathrm{\infty }\},$(61)

the space of square-integrable functions with square-integrable gradient. It can be shown that if a piecewise function $\mathit{u}\in C^0(\mathcal{D})$—the space of continuous functions with zero continuous derivatives—and satisfies $\mathit{u}{|}_K\in H^1(\mathit{K})$ for each $\mathit{K}\in \mathcal{T}$, then $\mathit{u}\in H^1(\mathcal{D})$ (see Ref. [30], Sec. 3.2.1). In other words, a finite-dimensional subspace of $H^1(\mathcal{D})$ must be continuous across interior mesh interfaces and have a locally smooth function space on each element. Thus, we take the degree-$\mathit{p}$ continuous finite element space to be

(62) $V_p=\{\mathit{u}\in C^0(\mathcal{D}):\mathit{u}{|}_K\in {\mathbb{Q}}_p(\mathit{K}),\mathrm{\forall K}\in \mathcal{T}\}.$(62)

Here, $V_p\subset H^1(\mathcal{D})$ since $\mathit{u}\in V_p$ is continuous and $\mathit{u}{|}_K\in {\mathbb{Q}}_p(\mathit{K})\subset H^1(\mathit{K})$ for each $\mathit{K}$. We use a nodal basis for ${\mathcal{Q}}_p(\hat{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.

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.

III.C.3. Raviart Thomas, Broken RT, and RT Trace Space

The RT space is a finite-dimensional subspace of the $\mathit{H}(\mathrm{div};\mathcal{D})$ space, where

(63) $\mathit{H}(\mathrm{div};\mathcal{D})=\{\mathit{v}\in [L^2(\mathcal{D}){]}^2:\mathrm{\nabla }\cdot \mathit{v}\in L^2(\mathcal{D})\}$(63)

is the space of square-integrable, vector-valued functions with square-integrable divergence.^[31,32] It can be shown that $\mathit{v}\in \mathit{H}(\mathrm{div};\mathcal{D})$ when $\mathit{v}$ satisfies the following: (1) $\mathit{v}$ is locally smooth on each element such that $\mathit{v}{|}_K\in [H^1(\mathit{K}){]}^2$ for each $\mathit{K}\in \mathcal{T}$ and (2) $\mathit{v}\cdot \mathbf{n}$ is continuous across each interior face. It is also desirable that the divergence of a vector in the degree-$\mathit{p}$ RT space belongs to the degree-$\mathit{p}$ DG space $Y_p$. Let ${\mathcal{Q}}_{\mathit{m},\mathit{n}}(\hat{K})$ denote the space of polynomials of at most $\mathit{m}$ and $\mathit{n}$ in the first and second variables, respectively, such that ${\mathcal{Q}}_{\mathit{p},\mathit{p}}(\hat{K})={\mathcal{Q}}_p(\hat{K})$. The RT space uses the local polynomial space ${\mathcal{Q}}_{\mathit{p}+1,\mathit{p}}(\hat{K})\times {\mathcal{Q}}_{\mathit{p},\mathit{p}+1}(\hat{K})$ on each element. In this way, the divergence of any vector $\mathit{v}\in {\mathcal{Q}}_{\mathit{p}+1,\mathit{p}}(\hat{K})\times {\mathcal{Q}}_{\mathit{p},\mathit{p}+1}(\hat{K})$ belongs to ${\mathcal{Q}}_p(\hat{K})$, the polynomial space used by $Y_p$ on each element. Finally, the contravariant Piola transformation, defined in Eq. (51)(51) $\mathit{v}=\frac{1}{J}\mathbf{F}\stackrel{\mathbf{ˆ}}{\mathit{v}}\circ {\mathbf{T}}^{-1},$(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-$\mathit{p}$ RT space is then

(64) $\begin{array}{r}\mathit{R}{T}_p=\{\mathit{v}\in [L^2(\mathcal{D})]:\mathit{v}{|}_K\in { }_p(\mathit{K}),\mathrm{\forall K}\in \mathcal{T}\\ \mathrm{and}\left[\mathit{v}\cdot \mathbf{n}\right]=0,\mathrm{\forall }\mathcal{F}\in {\mathrm{\Gamma }}_0\},\end{array}$(64)

where

(65) ${\mathbb{D}}_p(\mathit{K})=\{\mathit{v}=\frac{1}{J}\mathbf{F}\stackrel{\mathbf{ˆ}}{\mathit{v}}\circ {\mathbf{T}}^{-1}:\stackrel{\mathbf{ˆ}}{\mathit{v}}\in {\mathcal{Q}}_{\mathit{p}+1,\mathit{p}}(\hat{K})\times {\mathcal{Q}}_{\mathit{p},\mathit{p}+1}(\hat{K})\}$(65)

denotes the required RT element-local polynomial space composed with the contravariant Piola transformation. The constraint $\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]=0$ for each interior face holds only when $\mathit{v}\in \mathit{R}{T}_p$ has a continuous normal component. Thus, since ${\mathbb{D}}_p(\mathit{K})\subset [H^1(\mathit{K}){]}^2$ and each element of $\mathit{R}{T}_p$ has a continuous normal component, $\mathit{R}{T}_p\subset \mathit{H}(\mathrm{div};\mathcal{D})$.

Figure 5a depicts the degrees of freedom in $\mathit{R}{T}_1$ on a $3\times 3$ mesh. The black circles and red squares denote degrees of freedom representing the $\mathit{x}$ and $\mathit{y}$ components, respectively. Observe that on each element, the $\mathit{x}$ component is quadratic in the $\mathit{x}$ direction and linear in the $\mathit{y}$ direction with the opposite holding for the $\mathit{y}$ component. Furthermore, the degrees of freedom representing the normal component are shared between elements, ensuring that $\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]=0$ for each interior face and $\mathit{v}\in \mathit{R}{T}_p$.

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 $\mathit{x}$ and $\mathit{y}$ components, respectively. (b) The distribution of degrees of freedom corresponding to ${\mathrm{\Lambda }}_1$, the space defined as the interior trace of the first degree RT space, on a $3\times 3$ mesh. Observe that the interpolation for the normal component in (a) exactly matches the interpolation used in (b) on each face.

We also define two related spaces that arise in hybridizing a RT discretization. The first is the broken RT space ${\widehat{\mathit{RT}}}_{\mathit{p}}$, defined as the degree-$\mathit{p}$ RT space without the requirement of continuity in the normal component. That is,

(66) ${\widehat{\mathit{RT}}}_{\mathit{p}}=\{\mathit{v}\in [L^2(\mathcal{D})]:\mathit{v}{|}_K\in {\mathbb{D}}_p(\mathit{K}),\forall \mathit{K}\in \mathcal{T}\}.$(66)

Note that $\mathit{R}{T}_p\subset {\widehat{\mathit{RT}}}_{\mathit{p}}$ and that $\mathit{v}\in {\widehat{\mathit{RT}}}_{\mathit{p}}$ belongs to $\mathit{R}{T}_p$ if and only if $[[\mathit{v}.\mathit{n}]]\text{=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 $\mathit{R}{T}_p$ along the interior faces in the mesh. Let ${\mathcal{P}}_p(\hat{\mathcal{F}})$ be the space of univariate polynomials with degree at most $\mathit{p}$ defined over the reference line $\hat{\mathcal{F}}=[0,1]$ and

(67) $\mathbb{P}(\mathcal{F})=\{\mathit{u}=\hat{u}\circ {\mathbf{T}}^{-1}:\hat{u}\in {\mathcal{P}}_p(\hat{\mathcal{F}})\},$(67)

the associated polynomials mapped to physical space. The trace space is then

(68) ${\mathrm{\Lambda }}_p=\{\mathit{\mu }\in L^2({\mathrm{\Gamma }}_0):\mathit{\mu }{|}_{\mathcal{F}}\in {\mathbb{P}}_p(\mathcal{F}),\mathrm{\forall }\mathcal{F}\in {\mathrm{\Gamma }}_0\}.$(68)

The interior trace space for $\mathit{R}{T}_1$, ${\mathrm{\Lambda }}_1$, is depicted in Fig. 5b. Note that these degrees of freedom are exactly the degrees of freedom corresponding to the normal component of $\mathit{R}{T}_1$ on the interior faces of the mesh.

IV. TRANSPORT DISCRETIZATION

We discretize the ${\mathrm{S}}_N$ transport equations with a high-order DG spatial discretization compatible with curved meshes such as from Woods^[33] or Haut et al.^[34] In ${\mathrm{S}}_N$, the transport equation is collocated at discrete angles ${\mathbf{\Omega }}_d$, and integration over the unit sphere is approximated using a suitable angular quadrature rule $\{{\mathbf{\Omega }}_d,w_d{\}}_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}$. We use ${\psi }_d$ to denote the angular flux in discrete direction ${\mathbf{\Omega }}_d$ such that ${\psi }_d(\mathbf{x})=\mathit{\psi }(\mathbf{x},{\mathbf{\Omega }}_d)$. The DG discretization approximates each ${\psi }_d$ using the degree-$\mathit{p}$ DG finite element space $Y_p$. Through finite element interpolation, ${\psi }_d(\mathbf{x})$ can be evaluated at any point in the mesh. The SMM closures are computed using the ${\mathrm{S}}_N$ angular quadrature and finite element interpolation with

(69a) $\mathbf{T}(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\mathbf{\Omega }}_d\otimes {\mathbf{\Omega }}_d{\psi }_d(\mathbf{x})-\frac{1}{3}\mathbf{I}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}),$(69a)

(69b) $\mathit{\beta }(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|{\psi }_d(\mathbf{x})-E_{\mathit{b}0}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}).$(69b)

In light of Remark 3, we elect to use $E_{\mathit{b}0}$, as defined in Eq. (36)(36) $E_{\mathit{b}0}=\frac{\int |\mathbf{\Omega }\cdot \mathbf{n}|{\psi }_0\mathrm{d}\mathrm{\Omega }}{\int {\psi }_0\mathrm{d}\mathrm{\Omega }}=\frac{\int |\mathbf{\Omega }\cdot \mathbf{n}|\mathrm{d}\mathrm{\Omega }}{4\mathit{\pi }}=\frac{1}{2},$(36) , in place of the factor of $1/2$ that scales the zeroth moment term in the boundary correction. Here, $E_{\mathit{b}0}$ is approximated with ${\mathrm{S}}_N$ quadrature according to

(70) $E_{\mathit{b}0}=\frac{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|}{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d}=\frac{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|}{4\mathit{\pi }}.$(70)

With this discrete definition, the boundary conditions for the continuous system from Eq. (13c)(13c) $\mathit{J}\cdot \mathbf{n}=\frac{1}{2}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(13c) are now modified to

(71) $\mathit{J}\cdot \mathbf{n}=E_{\mathit{b}0}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },$(71)

where $\mathit{\beta }$ is defined using the discrete definition in Eq. (69b)(69b) $\mathit{\beta }(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|{\psi }_d(\mathbf{x})-E_{\mathit{b}0}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}).$(69b) . We also approximate the inflow partial current $J_{\mathrm{in}}$ using ${\mathrm{S}}_N$ quadrature with $J_{\mathrm{in}}(\mathbf{x})={\sum }_{{\mathbf{\Omega }}_d\cdot \mathbf{n}<0}{w}_d{\mathbf{\Omega }}_d\bar{\psi }(\mathbf{x},{\mathbf{\Omega }}_d)$.

In the limit as $N_{\mathrm{\Omega }}\to \mathrm{\infty }$, $E_{\mathit{b}0}\to 1/2,$ and the discrete $J_{\mathrm{in}}\to {\int }_{\mathbf{\Omega }\cdot \mathbf{n}<0}\mathbf{\Omega }\cdot \mathbf{n}\bar{\psi }\mathrm{\Omega }$ 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)(69a) $\mathbf{T}(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\mathbf{\Omega }}_d\otimes {\mathbf{\Omega }}_d{\psi }_d(\mathbf{x})-\frac{1}{3}\mathbf{I}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}),$(69a) and (69b)(69b) $\mathit{\beta }(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|{\psi }_d(\mathbf{x})-E_{\mathit{b}0}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}).$(69b) will be used.

Where needed, we compute the broken divergence of the SMM correction tensor using the finite element interpolation operator for $\mathbf{T}$. In other words,

(72) ${\mathrm{\nabla }}_h\cdot \mathbf{T}={\mathrm{\nabla }}_h\cdot \mathbf{P}-\frac{1}{3}{\mathrm{\nabla }}_h\mathit{\varphi },$(72)

where $\mathbf{P}$ and $\mathit{\varphi }$ 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, 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 $\mathbf{V}(\mathit{\psi },\mathit{\phi })=\mathit{f}$ represent the DG VEF moment system. Here, $\mathbf{V}$ is nonlinear in $\mathit{\psi }$, due to the nonlinear Eddington tensor and boundary factor, but linear in $\mathit{\phi }$. Further, we set ${\mathbf{y}}_0=[{\psi }_0,{\phi }_0]$ to be the linearly anisotropic approximation to the coupled transport-VEF system. The linearized operator is then

(73) $\frac{0=\mathit{V}(\mathit{\psi },\mathit{\phi })\text{ TSE}}{�V({\psi }_0,{\phi }_0)+\mathit{D}[\mathit{V}](y_0,\mathit{y}-y_0)}=\mathit{V}({\psi }_0,{\phi }_0)+\mathit{D}[\mathit{V}{|}_{\phi ={\phi }_0}]({\psi }_0,\mathit{\psi }-{\psi }_0)+\mathit{D}[\mathit{V}{|}_{\psi ={\psi }_0}]({\phi }_0,\mathit{\phi }-{\phi }_0)=\mathit{V}({\psi }_0,\mathit{\phi })+\mathit{D}[\mathit{V}{|}_{\phi ={\phi }_0}]({\psi }_0,\mathit{\psi }),(73)$(73)

where $\mathit{D}[\cdot ](\cdot ,\cdot )$ denotes the Gateaux derivative defined in Eq. (31)(31) $\mathit{D}[\mathit{f}](\mathit{u},\mathit{v})={\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\mathit{\omega }}\mathit{f}(\mathit{u}+\mathit{\omega v})\right]}_{\mathit{\omega }=0}.$(31) . We have simplified $\mathbf{V}({\psi }_0,{\phi }_0)+\mathit{D}[\mathbf{V}{|}_{\psi ={\psi }_0}]({\phi }_0,\mathit{\phi }-{\phi }_0)=\mathbf{V}({\psi }_0,\mathit{\phi })$ using the linearity of $\mathbf{V}$ in the second argument. Furthermore, the directional derivative of the Eddington tensor satisfies $\mathit{D}[\mathbf{E}]({\psi }_0,\mathit{\psi }-{\psi }_0)=\mathit{D}[\mathbf{E}]({\psi }_0,\mathit{\psi })$ with an analogous result for the directional derivative of the Eddington boundary factor. Since ${\psi }_0$ is linearly anisotropic, $\mathbf{V}({\psi }_0,\mathit{\phi })$ represents a radiation diffusion system. The correction sources are then given by $\mathit{D}[\mathbf{V}{|}_{\phi ={\phi }_0}]({\psi }_0,\mathit{\psi })$: the directional derivative of the VEF operator evaluated at ${\psi }_0$ in the direction $\mathit{\psi },$ where $\mathit{\phi }$ is fixed at the linearly anisotropic moment solution ${\phi }_0$. We proceed by evaluating the IP VEF method from Ref. [8], Sec. 5.2.1, at a linearly anisotropic transport solution ${\psi }_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)(59) $L^2(\mathcal{D})=\{\mathit{u}:\int u^2\mathrm{d}\mathbf{x}<\mathrm{\infty }\},$(59) , the IP VEF discretization is as follows: Find $\mathit{\phi }\in Y_p$ such that

(74) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_b\mathit{u\phi }\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \left(\mathbf{E}\mathit{\phi }\right)\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}-{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\frac{{\mathrm{\nabla }}_h\mathit{u}}{{\sigma }_t}\right\}\right\}\cdot \left[\left[\mathbf{E}\mathit{\phi }\mathbf{n}\right]\right]\mathrm{d}\mathit{s}\\ & +\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{1}{{\sigma }_t}\cdot \left(\mathbf{E}\mathit{\phi }\right)\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{\mathit{Q}{ }_1\cdot \mathbf{n}}{{\sigma }_t}\right\}\right\}\mathrm{d}\mathit{s}\\ & -2{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}{J}_{\mathrm{in}}\mathrm{d}\mathit{s},\mathrm{\forall u}\in Y_p,& (74)\end{array}$(74)

where $\mathit{\kappa }$ is the penalty parameter. This discretization is derived by approximating the scalar flux and each component of the current with a degree-$\mathit{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 ${\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}$ is added to stabilize the discretization. It is well known that $\mathit{\kappa }\propto {\sigma }_t^{-1}{p}^2/\mathit{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 $\mathit{h}$ and $\mathit{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

(75) $\mathbf{E}({\psi }_0)=\frac{1}{3}\mathbf{I},E_b({\psi }_0)=E_{\mathit{b}0},$(75)

where $E_{\mathit{b}0}$ is computed using the ${\mathrm{S}}_N$ quadrature rule following Eq. (70)(70) $E_{\mathit{b}0}=\frac{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|}{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d}=\frac{{\sum }_{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|}{4\mathit{\pi }}.$(70) . The use of numerical quadrature to evaluate this term is motivated by Remark 3. The diffusion operator is then

(76) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_{\mathit{b}0}\mathit{u\phi }\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{u}\cdot \mathbf{n}\right\}\right\}\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{{\mathit{Q}}_1\cdot \mathbf{n}}{{\sigma }_t}\right\}\right\}\mathrm{d}\mathit{s}\\ & -2{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}{J}_{\mathrm{in}}\mathrm{d}\mathit{s},& (76)\end{array}$(76)

where $\mathit{D}=\frac{1}{3{\sigma }_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 ${\mathit{Q}}_1$ since we have not made the assumption that the transport equation’s fixed source $\mathit{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)(74) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_b\mathit{u\phi }\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \left(\mathbf{E}\mathit{\phi }\right)\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}-{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\frac{{\mathrm{\nabla }}_h\mathit{u}}{{\sigma }_t}\right\}\right\}\cdot \left[\left[\mathbf{E}\mathit{\phi }\mathbf{n}\right]\right]\mathrm{d}\mathit{s}\\ & +\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{1}{{\sigma }_t}\cdot \left(\mathbf{E}\mathit{\phi }\right)\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{\mathit{Q}{ }_1\cdot \mathbf{n}}{{\sigma }_t}\right\}\right\}\mathrm{d}\mathit{s}\\ & -2{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}{J}_{\mathrm{in}}\mathrm{d}\mathit{s},\mathrm{\forall u}\in Y_p,& (74)\end{array}$(74) ] with $\mathit{\phi }={\phi }_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

(77) $\mathit{D}[{\int }_{{\mathrm{\Gamma }}_b}{E}_b\mathit{u}{\phi }_0\mathrm{d}\mathit{s}]({\psi }_0,\mathit{\psi })={\int }_{{\mathrm{\Gamma }}_b}\mathit{D}[E_b]({\psi }_0,\mathit{\psi })\mathit{u}{\phi }_0\mathrm{d}\mathit{s}={\int }_{{\mathrm{\Gamma }}_b}\mathit{u\beta }(\mathit{\psi })\mathrm{d}\mathit{s},$(77)

where $\mathit{\beta }(\mathit{\psi })$ is the correction factor defined in Eq. (69b)(69b) $\mathit{\beta }(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d|{\mathbf{\Omega }}_d\cdot \mathbf{n}|{\psi }_d(\mathbf{x})-E_{\mathit{b}0}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}).$(69b) that uses ${\mathrm{S}}_N$ quadrature to compute the angular moments of the discrete transport solution. We have used $\mathit{D}[E_b]({\psi }_0,\mathit{\psi })=\mathit{\beta }(\mathit{\psi })/{\varphi }_0$ and have made the approximation that ${\phi }_0/{\varphi }_0\approx 1$ (see Remark 1). Analogously, for terms with the Eddington tensor,

(78) $\mathit{D}[\mathbf{E}{\phi }_0]({\psi }_0,\mathit{\psi })=\frac{{\phi }_0}{{\varphi }_0}\mathbf{T}(\mathit{\psi })=\mathbf{T}(\mathit{\psi }),$(78)

where $\mathbf{T}(\mathit{\psi })$ is the correction tensor from Eq. (69a)(69a) $\mathbf{T}(\mathbf{x})=\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\mathbf{\Omega }}_d\otimes {\mathbf{\Omega }}_d{\psi }_d(\mathbf{x})-\frac{1}{3}\mathbf{I}\sum _{\mathit{d}=1}^{N_{\mathrm{\Omega }}}{w}_d{\psi }_d(\mathbf{x}),$(69a) and we have again made the approximation that ${\phi }_0/{\varphi }_0=1$. Thus, the correction terms can be derived by setting terms without angular flux dependence to zero and replacing

(79) $E_b\mathit{\phi }\to \mathit{\beta },\mathbf{E}\mathit{\phi }\to \mathbf{T}.$(79)

The IP SMM discretization is then: Find $\mathit{\phi }\in Y_p$ such that

(80) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_{\mathit{b}0}\mathit{u\phi }\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{u}\cdot \mathbf{n}\right\}\right\}\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{{\mathit{Q}}_1\cdot \mathbf{n}}{{\sigma }_t}\right\}\right\}\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}\left(2J_{\mathrm{in}}+\mathit{\beta }\right)\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \mathbf{T}\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}\\ & +{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\frac{{\mathrm{\nabla }}_h\mathit{u}}{{\sigma }_t}\right\}\right\}\cdot \left[\left[\mathbf{Tn}\right]\right]\mathrm{d}\mathit{s}-\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \mathbf{T}\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (80)\end{array}$(80)

where the local divergence of the correction tensor is computed with Eq. (72)(72) ${\mathrm{\nabla }}_h\cdot \mathbf{T}={\mathrm{\nabla }}_h\cdot \mathbf{P}-\frac{1}{3}{\mathrm{\nabla }}_h\mathit{\varphi },$(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)(80) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_{\mathit{b}0}\mathit{u\phi }\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\mathit{\kappa }\left[\left[\mathit{u}\right]\right]\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\mathit{D}{\mathrm{\nabla }}_h\mathit{u}\cdot \mathbf{n}\right\}\right\}\left[\left[\mathit{\phi }\right]\right]\mathrm{d}\mathit{s}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \mathit{D}{\mathrm{\nabla }}_h\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{{\mathit{Q}}_1\cdot \mathbf{n}}{{\sigma }_t}\right\}\right\}\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}\left(2J_{\mathrm{in}}+\mathit{\beta }\right)\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{u}\right]\right]\left\{\left\{\frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \mathbf{T}\cdot \mathbf{n}\right\}\right\}\mathrm{d}\mathit{s}\\ & +{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\frac{{\mathrm{\nabla }}_h\mathit{u}}{{\sigma }_t}\right\}\right\}\cdot \left[\left[\mathbf{Tn}\right]\right]\mathrm{d}\mathit{s}-\int {\mathrm{\nabla }}_h\mathit{u}\cdot \frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \mathbf{T}\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (80)\end{array}$(80) by setting the test and trial spaces to $V_p$, the degree-$\mathit{p}$ continuous finite element space. For any function $\mathit{v}\in V_p$,

(81) $\left[-0.15\mathit{em}\left[\mathit{v}\right]-0.15\mathit{em}\right]=0,\left\{\mathit{v}\right\}\text{}\}=\mathit{v},\forall \mathit{ℱ}\in {\Gamma }_0.$(81)

In addition, $\mathit{v}\in V_p$ satisfies ${\mathrm{\nabla }}_h\mathit{v}=\mathrm{\nabla 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 $\mathit{v}\in V_p$. Since DG is used for the transport discretization, the correction tensor is still generally discontinuous across interior mesh interfaces such that $\left[\left[\mathbf{Tn}\right]\right]$ is nonzero. Thus, the CG SMM moment discretization is: Find $\mathit{\phi }\in V_p$ such that

(82) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}{E}_{\mathit{b}0}\mathit{u\phi }\mathrm{d}\mathit{s}+\int \mathrm{\nabla u}\cdot \mathit{D}\mathrm{\nabla \phi }\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}\\ & =\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x}+\int \mathrm{\nabla }\mathit{u}\cdot \frac{{\mathit{Q}}_1}{{\sigma }_t}\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_b}\mathit{u}\left(2J_{\mathrm{in}}+\mathit{\beta }\right)\mathrm{d}\mathit{s}\\ & +{\int }_{{\mathrm{\Gamma }}_0}\left\{\left\{\frac{\mathrm{\nabla u}}{{\sigma }_t}\right\}\right\}\cdot \left[\left[\mathbf{Tn}\right]\right]\mathrm{d}\mathit{s}-\int \mathrm{\nabla u}\cdot \frac{1}{{\sigma }_t}{\mathrm{\nabla }}_h\cdot \mathbf{T}\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in V_p.& (82)\end{array}$(82)

The CG SMM moment discretization can also be derived directly from the CG VEF moment discretization in Ref. [8], Eq. (77)(77) $\mathit{D}[{\int }_{{\mathrm{\Gamma }}_b}{E}_b\mathit{u}{\phi }_0\mathrm{d}\mathit{s}]({\psi }_0,\mathit{\psi })={\int }_{{\mathrm{\Gamma }}_b}\mathit{D}[E_b]({\psi }_0,\mathit{\psi })\mathit{u}{\phi }_0\mathrm{d}\mathit{s}={\int }_{{\mathrm{\Gamma }}_b}\mathit{u\beta }(\mathit{\psi })\mathrm{d}\mathit{s},$(77) , via linearization or algebraically manipulating the closures.

V.B. Raviart Thomas and HRT

From Ref. [10], Eq. (61)(61) $H^1(\mathcal{D})=\{\mathit{u}\in L^2(\mathcal{D}):\int \mathrm{\nabla }\mathit{u}\cdot \mathrm{\nabla u}\mathrm{d}\mathbf{x}<\mathrm{\infty }\},$(61) , the RT VEF discretization is: Find $(\mathit{\phi },\mathit{J})\in Y_p\times \mathit{R}{T}_p$ such that

(83a) $\begin{array}{rl}& \int \mathit{u}\mathrm{\nabla }\cdot \mathit{J}\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}=\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (83\mathrm{a})\end{array}$(83a)

(83b) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \left\{\left\{\mathbf{En}\right\}\right\}\right]\right]\left\{\{\mathit{\phi }\right\}\}\mathrm{d}\mathit{s}-\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{E}\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_t\mathit{v}\cdot \mathit{J}\mathrm{d}\mathbf{x}\\ & +{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\left(\mathit{v}\cdot \mathbf{En}\right)\left(\mathit{J}\cdot \mathbf{n}\right)\mathrm{d}\mathit{s}=\int \mathit{v}\cdot {\mathit{Q}}_1\mathrm{d}\mathbf{x}\\ & +2{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\mathit{v}\cdot \mathbf{En}{J}_{\mathrm{in}}\mathrm{d}\mathit{s},\mathrm{\forall }\mathit{v}\in \mathit{R}{T}_p,& (83\mathrm{b})\end{array}$(83b)

where $\mathit{v},\mathit{J}\in \mathit{R}{T}_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 $\mathrm{\nabla }\cdot \left(\mathbf{E}\mathit{\phi }\right)$ prevents the straightforward application of standard mixed finite element techniques. To illustrate this, we multiply this divergence term by a test function $\mathit{v}\in \mathit{R}{T}_p$ and integrate over a single element $\mathit{K}$:

(84) ${\int }_K\mathit{v}\cdot \mathrm{\nabla }\cdot \left(\mathbf{E}\mathit{\phi }\right)\mathrm{d}\mathbf{x}={\int }_{\mathrm{\partial }\mathit{K}}\mathit{v}\cdot \widehat{\mathbf{E}\phi }\mathbf{n}\mathrm{d}\mathit{s}-{\int }_K\mathrm{\nabla }\mathit{v}:\mathbf{E}\mathit{\phi }\mathrm{d}\mathbf{x},$(84)

where $\widehat{\mathbf{E}\phi }\mathbf{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 element-boundary term gives rise to an interior face term of the form

(85) ${\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \widehat{\mathbf{E}\phi }\mathbf{n}\right]\right]\mathrm{d}\mathit{s},$(85)

along with a contribution on the boundary of the domain related to the boundary conditions that we omit for clarity of presentation. Note that $\mathit{v}\in \mathit{R}{T}_p$ has a continuous normal component but discontinuous tangential components. Thus, when $\mathbf{E}=\frac{1}{3}\mathbf{I}$, the interior face term vanishes since $\left[\left[\mathit{v}\cdot \mathbf{In}\right]\right]=\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]=0$. For a general Eddington tensor, however, this term remains, requiring the application of DG-like techniques to treat the discontinuity in $\mathit{v}\cdot \mathbf{En}$. We also note that the volumetric term $\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{E}\mathit{\phi }\mathrm{d}\mathbf{x}$ requires computing the gradient of an $\mathit{R}{T}_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 $\mathbf{E}=1/3\mathbf{I}$ 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:

(86) $\widehat{\mathbf{E}\phi }\mathbf{n}=\left\{\left\{\mathbf{En}\right\}\right\}\left\{\left\{\mathit{\phi }\right\}\right\}.$(86)

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 $\mathbf{P}=\mathbf{E}\mathit{\phi }$, we write the numerical flux term in terms of a numerical flux for the pressure denoted by $\hat{\mathbf{P}}\mathbf{n}$. By solving for $\mathit{\phi }$ in Eq. (12c)(12c) $\mathit{J}\cdot \mathbf{n}=E_b\mathit{\phi }+2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(12c) , the boundary scalar flux is given by

(87) $\mathit{\phi }=\frac{1}{E_b}\mathit{J}\cdot \mathbf{n}-2J_{\mathrm{in}},\mathbf{x}\in \mathrm{\partial }\mathcal{D},$(87)

so that

(88) ${\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\left(\mathit{v}\cdot \mathbf{En}\right)\left(\mathit{J}\cdot \mathbf{n}\right)\mathrm{d}\mathit{s}-2{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\mathit{v}\cdot \mathbf{En}{J}_{\mathrm{in}}\mathrm{d}\mathit{s}={\int }_{{\mathrm{\Gamma }}_b}\mathit{v}\cdot \mathbf{En}\mathit{\phi }\mathrm{d}\mathit{s}.$(88)

Thus, we can recast the RT VEF first moment [Eq. (83b)(83b) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \left\{\left\{\mathbf{En}\right\}\right\}\right]\right]\left\{\{\mathit{\phi }\right\}\}\mathrm{d}\mathit{s}-\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{E}\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_t\mathit{v}\cdot \mathit{J}\mathrm{d}\mathbf{x}\\ & +{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\left(\mathit{v}\cdot \mathbf{En}\right)\left(\mathit{J}\cdot \mathbf{n}\right)\mathrm{d}\mathit{s}=\int \mathit{v}\cdot {\mathit{Q}}_1\mathrm{d}\mathbf{x}\\ & +2{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{E_b}\mathit{v}\cdot \mathbf{En}{J}_{\mathrm{in}}\mathrm{d}\mathit{s},\mathrm{\forall }\mathit{v}\in \mathit{R}{T}_p,& (83\mathrm{b})\end{array}$(83b) ] in terms of the unclosed pressure as

(89) $\begin{array}{rl}& {\int }_{{\mathrm{\Gamma }}_b}\mathit{v}\cdot \mathbf{Pn}\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \hat{\mathbf{P}}\mathbf{n}\right]\right]\mathrm{d}\mathit{s}-\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{P}\mathrm{d}\mathbf{x}\\ & +\int {\sigma }_t\mathit{v}\cdot \mathit{J}\mathrm{d}\mathbf{x}=\int \mathit{v}\cdot {\mathit{Q}}_1\mathrm{d}\mathbf{x},\mathrm{\forall }\mathit{v}\in \mathit{R}{T}_p,& (89)\end{array}$(89)

where we have substituted the pressure in place of the product $\mathbf{E}\mathit{\phi }$. 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 $\mathbf{P}=\mathbf{T}+1/3\mathbf{I}\mathit{\phi }$. The volumetric term is then

(90) $\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{P}\mathrm{d}\mathbf{x}=\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{T}\mathrm{d}\mathbf{x}+\frac{1}{3}\int \mathrm{\nabla }\cdot \mathit{v}\mathit{\phi }\mathrm{d}\mathbf{x}.$(90)

Here, ${\mathrm{\nabla }}_h\mathit{v}:\mathbf{I}\mathit{\phi }={\mathrm{\nabla }}_h\cdot \mathit{v}=\mathrm{\nabla }\cdot \mathit{v}$ since ${\mathrm{\nabla }}_h\cdot \mathit{v}=\mathrm{\nabla }\cdot \mathit{v}$ for $\mathit{v}\in \mathit{R}{T}_p$ (Ref. [30], Proposition 3.2.2). For the numerical flux, we simply use the average:

(91) $\hat{\mathbf{P}}\mathbf{n}=\left\{\left\{\mathbf{Tn}+\frac{1}{3}\mathit{\phi }\mathbf{n}\right\}\right\}=\left\{\left\{\mathbf{Tn}\right\}\right\}+\frac{1}{3}\left\{\left\{\mathit{\phi }\right\}\right\}\mathbf{n}.$(91)

With this choice of numerical flux, the interior face bilinear form simplifies to

(92) ${\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \hat{\mathbf{P}}\mathbf{n}\right]\right]\mathrm{d}\mathit{s}={\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\right]\right]\cdot \left\{\left\{\mathbf{Tn}\right\}\right\}\mathrm{d}\mathit{s}+\frac{1}{3}{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]\left\{\left\{\mathit{\phi }\right\}\right\}\mathrm{d}\mathit{s}={\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\right]\right]\cdot \left\{\left\{\mathbf{Tn}\right\}\right\}\mathrm{d}\mathit{s},$(92)

where $\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]=0$ since the normal component is continuous for $\mathit{v}\in \mathit{R}{T}_p$. Finally, on the boundary of the domain, we solve the SMM boundary conditions in Eq. (71)(71) $\mathit{J}\cdot \mathbf{n}=E_{\mathit{b}0}\mathit{\phi }+2J_{\mathrm{in}}+\mathit{\beta },$(71) for the scalar flux and insert this relationship into $\mathbf{Pn}=\mathbf{Tn}+\mathbf{n}/3\mathit{\phi }$ to yield

(93) $\begin{array}{rl}& \mathbf{Pn}=\mathbf{Tn}+\frac{1}{3}\mathit{\phi }\mathbf{n}=\mathbf{Tn}+\frac{\mathbf{n}}{3{\mathit{E}}_{\mathit{b}0}}\left(\mathit{J}\cdot \mathbf{n}-2J_{\mathrm{in}}-\mathit{\beta }\right),\\ & \mathbf{x}\in \mathrm{\partial }\mathcal{D}.& (93)\end{array}$(93)

Thus, the RT SMM discretization is: Find $(\mathit{\phi },\mathit{J})\in Y_p\times \mathit{R}{T}_p$ such that

(94a) $\begin{array}{rl}& \int \mathit{u}\mathrm{\nabla }\cdot \mathit{J}\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}=\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (94\mathrm{a})\end{array}$(94a)

(94b) $\begin{array}{rl}& -\frac{1}{3}\int \mathrm{\nabla }\cdot \mathit{v}\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_t\mathit{v}\cdot \mathit{J}\mathrm{d}\mathbf{x}+{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{3{\mathit{E}}_{\mathit{b}0}}(\mathit{v}\cdot \mathbf{n})(\mathit{J}\cdot \mathbf{n})\mathrm{d}\mathit{s}\\ & =\int \mathit{v}\cdot {\mathit{Q}}_1\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_b}\mathit{v}\cdot \mathbf{Tn}\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{3{\mathit{E}}_{\mathit{b}0}}\mathit{v}\cdot \mathbf{n}\left(2J_{\mathrm{in}}+\mathit{\beta }\right)\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\right]\right]\cdot \left\{\left\{\mathbf{Tn}\right\}\right\}\mathrm{d}\mathit{s}+\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{T}\mathrm{d}\mathbf{x},\mathrm{\forall }\mathit{v}\in \mathit{R}{T}_p.& (94\mathrm{b})\end{array}$(94b)

Observe that the left-hand side of Eqs. (94)(94a) $\begin{array}{rl}& \int \mathit{u}\mathrm{\nabla }\cdot \mathit{J}\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}=\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (94\mathrm{a})\end{array}$(94a) is equivalent to an RT discretization of radiation diffusion with Marshak boundary conditions.

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 ${\widehat{\mathit{RT}}}_{\mathit{p}}$, and continuity of the current in the normal component is enforced with a Lagrange multiplier defined on the trace space of $\mathit{R}{T}_p$, ${\mathrm{\Lambda }}_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 element-local back substitution. A derivation of a hybridized radiation diffusion method is outlined in Ref. [10], Sec. 7.1.

The hybridized system is: Find $(\mathit{J},\mathit{\phi },\mathit{\lambda })\in {\widehat{\mathit{RT}}}_{\mathit{p}}\times Y_p\times {\mathrm{\Lambda }}_p$ such that

(95a) $\int \mathit{u}{\mathrm{\nabla }}_h\cdot \mathit{J}\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}=\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x},\mathrm{\forall }\mathit{u}\in Y_p,$(95a)

(95b) $\begin{array}{rl}& -\frac{1}{3}\int {\mathrm{\nabla }}_h\cdot \mathit{v}\mathit{\phi }\mathrm{d}\mathbf{x}+\int {\sigma }_t\mathit{v}\cdot \mathit{J}\mathrm{d}\mathbf{x}+{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{3{\mathit{E}}_{\mathit{b}0}}(\mathit{v}\cdot \mathbf{n})(\mathit{J}\cdot \mathbf{n})\mathrm{d}\mathit{s}\\ & +{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\cdot \mathbf{n}\right]\right]\mathit{\lambda }\mathrm{d}\mathit{s}=\mathcal{S},\mathrm{\forall }\mathit{v}\in {\widehat{\mathit{RT}}}_{\mathit{p}}& (95\mathrm{b})\end{array}$(95b)

(95c) ${\int }_{{\mathrm{\Gamma }}_0}\mathit{\mu }\left[\left[\mathit{J}\cdot \mathbf{n}\right]\right]\mathrm{d}\mathit{s}=0,\mathrm{\forall \mu }\in {\mathrm{\Lambda }}_p,$(95c)

where

(96) $\begin{array}{rl}& \mathcal{S}=\int \mathit{v}\cdot {\mathit{Q}}_1\mathrm{d}\mathbf{x}-{\int }_{{\mathrm{\Gamma }}_b}\mathit{v}\cdot \mathbf{Tn}\mathrm{d}\mathit{s}+{\int }_{{\mathrm{\Gamma }}_b}\frac{1}{3{\mathit{E}}_{\mathit{b}0}}\mathit{v}\cdot \mathbf{n}\left(2J_{\mathrm{in}}+\mathit{\beta }\right)\mathrm{d}\mathit{s}\\ & -{\int }_{{\mathrm{\Gamma }}_0}\left[\left[\mathit{v}\right]\right]\cdot \left\{\left\{\mathbf{Tn}\right\}\right\}\mathrm{d}\mathit{s}+\int {\mathrm{\nabla }}_h\mathit{v}:\mathbf{T}\mathrm{d}\mathbf{x}& (96)\end{array}$(96)

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)(95c) ${\int }_{{\mathrm{\Gamma }}_0}\mathit{\mu }\left[\left[\mathit{J}\cdot \mathbf{n}\right]\right]\mathrm{d}\mathit{s}=0,\mathrm{\forall \mu }\in {\mathrm{\Lambda }}_p,$(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] The 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 ${\mathrm{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

(97) $\mathit{\kappa }=\left\{\left\{\frac{{(\mathit{p}+1)}^2}{{\sigma }_t\mathit{h}}\right\}\right\},$(97)

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

(98) $\left[\begin{array}{cc}{\mathbf{M}}_t& \mathbf{G}\\ \mathbf{D}& {\mathbf{M}}_a\end{array}\right]\left[\begin{array}{c}\mathit{J}\\ \mathit{\phi }\end{array}\right]=\left[\begin{array}{c}\mathit{g}\\ \mathit{f}\end{array}\right],$(98)

where ${\mathbf{M}}_t$ = total interaction mass matrix; ${\mathbf{M}}_a$ = absorption mass matrix; $\mathbf{D}$ and $\mathbf{G}$ = discrete divergence and gradient, respectively; $\mathit{f}$ and $\mathit{g}$ = source terms. For SMM, $\mathbf{G}=-1/3{\mathbf{D}}^T$. For VEF, $\mathbf{G}$ depends on the Eddington tensor, and thus, $\mathbf{G}$ is in general not proportional to the transpose of the divergence matrix $\mathbf{D}$. We use block diagonal and lower block triangular preconditioners of the form

(99) ${\mathbf{P}}_{\mathrm{diag}}=\left[\begin{array}{cc}{\mathbf{M}}_t& \mathit{\hspace{0.17em}}\\ \mathit{\hspace{0.17em}}& \stackrel{\mathbf{ˆ}}{\mathbf{S}}\end{array}\right],{\mathbf{P}}_{\mathrm{tri}}=\left[\begin{array}{cc}{\mathbf{M}}_t& \mathit{\hspace{0.17em}}\\ \mathbf{D}& \stackrel{\mathbf{ˆ}}{\mathbf{S}}\end{array}\right],$(99)

where $\stackrel{\mathbf{ˆ}}{\mathbf{S}}={\mathbf{M}}_a-{\mathbf{D}\stackrel{\mathbf{ˆ}}{\mathbf{M}}}_{\mathit{t}}^{-1}\mathbf{G}$ is an approximate Schur complement formed using the lumped total interaction mass matrix ${\stackrel{\mathbf{ˆ}}{\mathbf{M}}}_{\mathit{t}}$ computed by summing the off-diagonal entries of ${\mathbf{M}}_t$ into the diagonal and setting the off-diagonal entries to zero. Lumping an $\mathit{R}{T}_p$ mass matrix produces a diagonal matrix, and thus, ${\stackrel{\mathbf{ˆ}}{\mathbf{M}}}_{\mathit{t}}^{-1}$ 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)(98) $\left[\begin{array}{cc}{\mathbf{M}}_t& \mathbf{G}\\ \mathbf{D}& {\mathbf{M}}_a\end{array}\right]\left[\begin{array}{c}\mathit{J}\\ \mathit{\phi }\end{array}\right]=\left[\begin{array}{c}\mathit{g}\\ \mathit{f}\end{array}\right],$(98) by negative three so that the off-diagonal blocks are symmetric and the scaled ${\mathbf{M}}_t$ is negative definite. For the block diagonal preconditioner, ${\mathbf{M}}_t$ is scaled by positive three so that ${\mathbf{P}}_{\mathrm{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.

VI.A. Method of Manufactured Solutions

The accuracy of the methods is determined with the MMS. The solution is set to

(100) $\mathit{\psi }=\frac{1}{4\mathit{\pi }}\mathit{\alpha }(\mathbf{x})+\mathbf{\Omega }\cdot \mathit{\beta }(\mathbf{x})+\mathbf{\Omega }\otimes \mathbf{\Omega }:\mathrm{\Theta }(\mathbf{x}),$(100)

where

(101a) $\mathit{\alpha }(\mathbf{x})=sin(\mathit{\pi x})sin(\mathit{\pi y})+\mathit{\delta },$(101a)

(101b) $\mathit{\beta }(\mathbf{x})=\left[\begin{array}{c}sin\frac{2\mathit{\pi }(\mathit{x}+\mathit{\omega })}{1+2\mathit{\omega }}sin\frac{2\mathit{\pi }(\mathit{y}+\mathit{\omega })}{1+2\mathit{\omega }}\\ \mathit{\hspace{0.17em}}\\ sin\frac{2\mathit{\pi }(\mathit{x}+\mathit{\omega })}{1+2\mathit{\omega }}sin\frac{2\mathit{\pi }(\mathit{y}+\mathit{\omega })}{1+2\mathit{\omega }}\end{array}\right],$(101b)

(101c) $\mathrm{\Theta }(\mathbf{x})=\left[\begin{array}{cc}\frac{1}{2}sin\frac{3\mathit{\pi }(\mathit{x}+\mathit{\zeta })}{1+2\mathit{\zeta }}sin\frac{3\mathit{\pi }(\mathit{y}+\mathit{\zeta })}{1+2\mathit{\zeta }}& sin\frac{2\mathit{\pi }(\mathit{x}+\mathit{\omega })}{1+2\mathit{\omega }}sin\frac{2\mathit{\pi }(\mathit{y}+\mathit{\omega })}{1+2\mathit{\omega }}\\ sin\frac{2\mathit{\pi }(\mathit{x}+\mathit{\omega })}{1+2\mathit{\omega }}sin\frac{2\mathit{\pi }(\mathit{y}+\mathit{\omega })}{1+2\mathit{\omega }}& \frac{1}{4}sin\frac{3\mathit{\pi }(\mathit{x}+\mathit{\zeta })}{1+2\mathit{\zeta }}sin\frac{3\mathit{\pi }(\mathit{y}+\mathit{\zeta })}{1+2\mathit{\zeta }}\end{array}\right].$(101c)

Here, $\mathit{\delta }=1.25$ is used to ensure $\mathit{\psi }>0,$ and $\mathit{\zeta }=0.1$ and $\mathit{\omega }=0.05$ are used to test spatially dependent, nonisotropic inflow boundary conditions. The domain is $\mathcal{D}=[0,1{]}^2$. With these definitions

(102a) $\mathit{\varphi }(\mathbf{x})=\mathit{\alpha }(\mathbf{x})+\frac{1}{3}\mathrm{trace}\mathrm{\Theta }(\mathbf{x}),$(102a)

(102b) $\mathit{J}(\mathbf{x})=\frac{1}{3}\mathit{\beta }(\mathbf{x}),$(102b)

(102c) $\mathbf{P}(\mathbf{x})=\frac{\mathit{\alpha }(\mathbf{x})}{3}\mathbf{I}+\frac{1}{15}\left[\begin{array}{cc}3{\mathrm{\Theta }}_{11}(\mathbf{x})+{\mathrm{\Theta }}_{22}(\mathbf{x})& {\mathrm{\Theta }}_{12}(\mathbf{x})\\ {\mathrm{\Theta }}_{21}(\mathbf{x})& {\mathrm{\Theta }}_{11}(\mathbf{x})+3{\mathrm{\Theta }}_{22}(\mathbf{x})\end{array}\right].$(102c)

This leads to an exact Eddington tensor $\mathbf{E}=\mathbf{P}/\mathit{\varphi }$ that is dense and spatially varying. The MMS $\mathit{\psi }$ and $\mathit{\varphi }$ are substituted into the transport equation to solve for the MMS source $\mathit{q}$ that forces the solution to Eq. (100)(100) $\mathit{\psi }=\frac{1}{4\mathit{\pi }}\mathit{\alpha }(\mathbf{x})+\mathbf{\Omega }\cdot \mathit{\beta }(\mathbf{x})+\mathbf{\Omega }\otimes \mathbf{\Omega }:\mathrm{\Theta }(\mathbf{x}),$(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 ${\mathit{Q}}_0$ and ${\mathit{Q}}_1$ to the moments of the MMS source. This is accomplished by projecting the MMS angular flux onto a degree-$\mathit{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_{\mathrm{in}}$. The SMM moment systems are then solved as if $\mathbf{T}$, $\mathit{\beta }$, ${\mathit{Q}}_0$, ${\mathit{Q}}_1$, and $J_{\mathrm{in}}$ are given data. Errors are calculated with the $L^2(\mathcal{D})$ norm for scalars and the $[L^2(\mathcal{D}){]}^2$ norm for vectors given by

(103) $\parallel \mathit{u}\parallel =\sqrt{\int u^2\mathbf{x}},$(103)

and

(104) $\parallel \mathit{v}\parallel =\sqrt{\int \mathit{v}\cdot \mathit{v}\mathrm{d}\mathbf{x}},$(104)

respectively. We also use the $L^2(\mathcal{D})$ projection operator ${\mathrm{\Pi }}_p:L^2(\mathcal{D})\to Y_p$ such that

(105) $\int \mathit{u}(\mathit{v}-{\mathrm{\Pi }}_p\mathit{v})\mathrm{d}\mathbf{x}=0,\mathrm{\forall }\mathit{u}\in Y_p,$(105)

for some $\mathit{v}\in L^2(\mathcal{D})$. In particular, ${\mathrm{\Pi }}_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

(106) $\mathbf{x}={\int }_0^T\mathbf{v}\mathrm{d}\mathit{t},$(106)

where the final time $\mathit{T}=0.3\mathit{\pi }$ and

(107) $\mathbf{v}=\left[\begin{array}{c}sin(\mathit{x})cos(\mathit{y})\\ -cos(\mathit{x})sin(\mathit{y})\end{array}\right]$(107)

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.

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.

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 $\mathcal{O}(h^{\mathit{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.

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.

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(\mathcal{D})$ projection operator ${\mathrm{\Pi }}_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 $\mathcal{O}(h^{\mathit{p}+2}),$ and the current should be $\mathcal{O}(h^{\mathit{p}+1})$. However, on this quadratically anisotropic transport problem, the superconvergence property is lost as the projected error measure converges at only $\mathcal{O}(h^{\mathit{p}+1})$. Furthermore, the current converges at only $\mathcal{O}(h^p)$ for odd polynomial degrees and $\mathcal{O}(h^{\mathit{p}+1/2})$ for even. This behavior was also seen in Ref. [10] for the RT 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] This may be due to the ability to exactly integrate the closure terms whereas with VEF the closures cannot be integrated exactly (see Remark 5).

TABLE I MMS Order of Accuracy and Constants

$\mathit{p}$	Value	$\parallel \mathit{\phi }-{\phi }_{\mathrm{ex}}\parallel$		$\parallel \mathit{\phi }-\mathrm{\Pi }{\phi }_{\mathrm{ex}}\parallel$		$\parallel \mathit{J}-{\mathit{J}}_{\mathrm{ex}}\parallel$
		RT	HRT	RT	HRT	RT	HRT
1	Order	2.002	2.002	2.175	2.175	0.993	0.993
	Constant	0.608	0.608	0.145	0.145	0.439	0.439
2	Order	2.989	2.989	2.964	2.964	2.521	2.521
	Constant	0.396	0.396	0.118	0.118	0.605	0.605
3	Order	4.006	4.006	4.254	4.254	2.971	2.971
	Constant	0.309	0.309	0.098	0.098	0.286	0.286

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

(108) ${\sigma }_t=1/\mathit{ϵ},{\sigma }_a=\mathit{ϵ},{\sigma }_s=1/\mathit{ϵ}-\mathit{ϵ},\mathit{q}=\mathit{ϵ},$(108)

where $\mathit{ϵ}\in (0,1]$ is a parameter that characterizes the mean free path and the thick diffusion limit corresponds to the limit $\mathit{ϵ}\to 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\times 8$ mesh with $\mathcal{D}=[0,1{]}^2$. The second is the triple point mesh shown in Fig. 8. On the triple point mesh, the angular flux 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₄ angular quadrature. The SMM and VEF algorithms are compared using $\mathit{p}=2$. The coupled transport-moment systems are solved with fixed-point iteration to a tolerance of ${10}^{-6}$.

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.

Table II shows the number of iterations to convergence on the $8\times 8$ orthogonal mesh as $\mathit{ϵ}\to 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 shown in Fig. 9 demonstrating that the methods are converging to the physically realistic diffusion solution as $\mathit{ϵ}\to 0$.

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

$\mathit{ϵ}$	IP		CG		RT		HRT
	VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
${10}^{-1}$	8	10	8	10	8	10	8	10
${10}^{-2}$	6	8	6	8	6	8	6	8
${10}^{-3}$	4	5	4	5	4	6	4	6
${10}^{-4}$	3	4	3	4	3	4	3	4

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

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)(94a) $\begin{array}{rl}& \int \mathit{u}\mathrm{\nabla }\cdot \mathit{J}\mathrm{d}\mathbf{x}+\int {\sigma }_a\mathit{u\phi }\mathrm{d}\mathbf{x}=\int \mathit{u}{Q}_0\mathrm{d}\mathbf{x},\\ & \mathrm{\forall }\mathit{u}\in Y_p,& (94\mathrm{a})\end{array}$(94a) , 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 $\mathit{\epsilon }$ for both VEF and SMM. In particular, RT and HRT SMM required five more iterations to converge compared to IP and CG for both $\mathit{ϵ}={10}^{-2}$ and $\mathit{ϵ}={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 $\mathit{ϵ}\to 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 $\mathit{ϵ}$.

TABLE III Iterations to Convergence in the Thick Diffusion Limit on the Triple Point Mesh

$\mathit{ϵ}$	IP		CG		RT		HRT
	VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
${10}^{-1}$	18	21	18	21	21	20	21	20
${10}^{-2}$	11	12	11	12	19	17	19	17
${10}^{-3}$	8	8	8	8	13	13	13	13
${10}^{-4}$	6	6	6	6	8	8	8	8

Fig. 10. Lineouts of the two-dimensional solution at $\mathit{x}=3.5$ as $\mathit{ϵ}\to 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.

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 $\mathit{c}\mathrm{\Delta t}={10}^3$ and the initial condition is ${\psi }_0={10}^{-4}$. The absorption and source are then ${\sigma }_a=1/\mathit{c}\mathrm{\Delta t}={10}^{-3}\frac{1}{\mathrm{cm}}$ and $\mathit{q}={\psi }_0/\mathit{c}\mathrm{\Delta t}={10}^{-1}\frac{1}{\mathrm{c}{\mathrm{m}}^3\cdot \mathrm{s}\cdot \mathrm{str}}$. The boundary conditions are set so that isotropic inflow of magnitude $1/2\mathit{\pi }$ enters on the left entrance of the pipe with vacuum on all other surfaces. A Level Symmetric S₁₂ 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.

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

The efficiencies of the outer fixed-point and inner preconditioned linear iterations are investigated by refining in $\mathit{h}$ and $\mathit{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 $\mathit{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, 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 Anderson-accelerated 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.

TABLE IV Anderson-Accelerated Fixed-Point Iterations on the Crooked Pipe Problem

	$N_e$	IP		CG		RT		HRT
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
$\mathit{p}=1$	112	10	9	10	10	12	12	12	12
	448	11	11	11	11	13	12	13	12
	1792	13	12	13	13	16	14	16	14
	7168	14	14	14	14	16	15	16	15
$\mathit{p}=2$	112	13	13	13	14	14	13	14	13
	448	14	14	14	14	16	15	16	15
	1792	15	14	15	14	16	15	16	15
	7168	15	14	15	14	16	16	16	16
$\mathit{p}=3$	112	14	14	14	14	16	14	16	14
	448	15	14	16	14	16	15	16	15
	1792	15	16	15	16	16	16	16	16
	7168	15	18	16	18	17	15	17	15

The average number of inner iterations per outer iteration is presented in Table V. All methods were scalable in $\mathit{h}$ and $\mathit{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 $\mathit{p}=2$ and $\mathit{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 V The Average Number of Inner, Preconditioned Linear Iterations per Outer Iteration

	$N_e$	IP		CG		RT		HRT
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
$\mathit{p}=1$	112	11.27	19.70	5.18	7.82	10.23	28.92	4.38	8.50
	448	11.75	18.75	5.33	8.92	11.36	37.23	4.64	9.83
	1792	11.64	19.62	5.36	9.71	11.24	36.67	4.76	9.43
	7168	11.93	18.73	5.13	9.60	12.35	37.12	5.06	9.73
$\mathit{p}=2$	112	10.93	18.71	6.57	10.47	15.13	50.71	6.27	12.46
	448	11.73	18.67	6.53	11.20	16.41	54.75	6.47	12.33
	1792	12.06	19.00	6.69	11.00	17.24	57.88	7.24	14.47
	7168	12.06	20.40	6.75	12.00	18.35	60.06	7.82	15.19
$\mathit{p}=3$	112	13.73	24.60	8.00	13.60	15.00	50.00	5.71	11.93
	448	14.44	25.53	8.65	14.40	17.12	51.06	5.65	11.87
	1792	15.31	26.00	9.38	14.65	18.00	52.53	5.82	11.56
	7168	15.69	25.26	9.88	14.63	18.50	54.38	6.33	13.27

TABLE VI The Average Amount of Time Spent Solving the Moment Systems per Outer Iteration

	$N_e$	IP		CG		RT		HRT
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
$\mathit{p}=1$	112	2.15	2.09	1.62	1.59	5.20	4.88	2.30	2.24
	448	7.89	7.48	6.01	5.91	22.33	24.73	8.97	9.03
	1792	30.75	29.62	23.57	23.10	89.39	101.30	36.57	36.76
	7168	125.88	114.97	93.18	93.30	398.94	402.85	149.14	146.81
$\mathit{p}=2$	112	4.41	4.20	2.79	2.64	22.57	27.52	3.57	3.49
	448	17.33	15.58	10.73	10.52	99.97	126.22	13.82	13.76
	1792	72.03	64.75	43.52	41.12	440.46	525.82	60.56	59.62
	7168	286.13	268.64	173.53	167.71	2053.12	2377.09	246.79	247.66
$\mathit{p}=3$	112	11.65	11.13	6.54	6.16	65.32	80.92	4.92	5.15
	448	48.72	46.28	26.39	25.11	300.24	340.78	21.28	21.19
	1792	202.41	181.64	112.72	102.71	1415.54	1518.80	85.99	85.48
	7168	903.56	772.37	471.34	400.64	6312.46	6829.70	364.36	361.74

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 $\mathit{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.

TABLE VII The Average Amount of Time Spent Assembling the Moment Systems per Outer Iteration

	$N_e$	IP		CG		RT		HRT
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
$\mathit{p}=1$	112	13.08	6.98	13.07	5.64	22.51	15.31	22.88	21.71
	448	49.89	25.08	49.78	17.52	70.83	50.81	72.37	63.78
	1792	195.30	96.38	198.01	59.61	242.67	182.26	256.36	207.78
	7168	784.84	364.82	788.22	220.13	932.15	678.96	963.29	730.83
$\mathit{p}=2$	112	24.91	11.32	24.81	8.46	39.69	31.76	43.39	40.83
	448	95.98	42.90	97.21	28.86	133.31	101.62	145.04	119.96
	1792	384.43	169.98	395.10	107.31	496.14	366.28	558.86	395.79
	7168	1549.05	678.84	1586.76	412.36	1947.54	1344.27	2157.27	1383.17
$\mathit{p}=3$	112	46.98	19.69	46.85	15.15	75.94	50.64	94.05	74.56
	448	189.96	74.65	185.81	52.21	274.23	178.82	353.88	223.77
	1792	747.61	289.70	765.16	193.40	1078.80	654.50	1393.50	743.02
	7168	3110.29	1131.54	3152.60	732.44	4465.49	2543.70	5585.02	2630.79

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 $\mathit{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.

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

	$N_e$	IP		CG		RT		HRT
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
$\mathit{p}=1$	112	2.08	1.82	2.07	1.96	2.64	2.49	2.55	2.55
	448	7.76	7.45	7.74	7.34	9.63	8.57	9.28	8.56
	1792	34.20	30.89	34.71	32.31	43.49	37.60	42.96	37.39
	7168	148.10	139.61	146.48	137.63	176.14	157.74	170.21	155.24
$\mathit{p}=2$	112	4.30	4.11	4.27	4.35	5.20	4.67	4.89	4.51
	448	16.60	15.88	16.53	15.65	21.27	19.54	19.65	18.43
	1792	71.46	63.11	71.05	61.38	85.62	76.62	78.42	70.60
	7168	284.71	251.91	283.51	247.63	343.61	328.58	310.90	294.44
$\mathit{p}=3$	112	9.19	8.85	9.14	8.66	12.21	10.35	11.07	9.61
	448	39.06	34.57	40.05	33.67	47.79	42.81	42.90	38.48
	1792	153.92	153.98	152.63	151.51	194.26	181.54	172.90	161.17
	7168	643.23	711.21	668.00	700.95	858.34	719.65	749.42	620.24

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 $\pm 1$ iterations indicating that SMM can converge robustly without a negative flux fixup.

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

	$N_e$	IP		CG		RT		HRT
		Fixup	No Fixup	Fixup	No Fixup	Fixup	No Fixup	Fixup	No Fixup
$\mathit{p}=1$	112	9	10	10	10	12	12	12	12
	448	11	11	11	11	12	12	12	12
	1792	12	13	13	13	14	14	14	14
	7168	14	14	14	14	15	15	15	15
$\mathit{p}=2$	112	13	13	14	12	13	13	13	13
	448	14	14	14	14	15	15	15	15
	1792	14	14	14	14	15	16	15	16
	7168	14	15	14	15	16	16	16	16
$\mathit{p}=3$	112	14	14	14	14	14	15	14	15
	448	14	14	14	14	15	15	15	15
	1792	16	15	16	15	16	16	16	16
	7168	18	17	18	17	15	16	15	16

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 ${\mathrm{S}}_N$ discretization and DSA preconditioner of Haut et al.^[21] as the mesh is refined. We use the thick diffusion limit problem from Sec. VI.B with $\mathit{ϵ}={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 ${\mathrm{S}}_N$ solutions in isolation be written as follows:

(109a) $\parallel {\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{ex}}\parallel =C_{\mathrm{SMM}}{h}^{\mathit{p}+1},$(109a)

(109b) $\parallel {\phi }_{\mathrm{SN}}-{\phi }_{\mathrm{ex}}\parallel =C_{\mathrm{SN}}{h}^{\mathit{p}+1},$(109b)

where ${\phi }_{\mathrm{ex}}$ is the exact solution; ${\phi }_{\mathrm{SMM}}$ and ${\phi }_{\mathrm{SN}}$ are the solutions generated by SMM and ${\mathrm{S}}_N$, respectively; $C_i$ are the error constants; $\mathit{h}$ is the mesh size; $\mathit{p}$ is the finite element polynomial degree. Using the triangle inequality, observe that

(110) $\begin{array}{cc}\parallel {\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{SN}}\parallel & =\parallel ({\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{ex}})+({\phi }_{\mathrm{ex}}-{\phi }_{\mathrm{SN}})\parallel \\ \mathit{\hspace{0.17em}}& \le \parallel {\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{ex}}\parallel +\parallel {\phi }_{\mathrm{SN}}-{\phi }_{\mathrm{ex}}\parallel \\ \mathit{\hspace{0.17em}}& =(C_{\mathrm{SMM}}+C_{\mathrm{SN}})h^{\mathit{p}+1}.\end{array}$(110)

Thus, we expect the difference between the solutions generated by SMM and ${\mathrm{S}}_N$ to converge with order $\mathit{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}^{\mathit{p}+1}$. Because of this, 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.

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(\mathcal{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 $\mathit{ϵ}={10}^{-2}$. Both methods used $\mathit{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 }.$

Figure 12b shows the $L^2(\mathcal{D})$ norm between the SMM and ${\mathrm{S}}_N$ solutions as a function of $h_{\max }$, the maximum characteristic mesh length in the mesh. We use Level Symmetric S₄ angular quadrature and $\mathit{p}=2$. The solutions are computed on four meshes built from 61, 81, 101, and 121 Chebyshev points in each direction. Using $\mathit{p}=2$, Eq. (110)(110) $\begin{array}{cc}\parallel {\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{SN}}\parallel & =\parallel ({\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{ex}})+({\phi }_{\mathrm{ex}}-{\phi }_{\mathrm{SN}})\parallel \\ \mathit{\hspace{0.17em}}& \le \parallel {\phi }_{\mathrm{SMM}}-{\phi }_{\mathrm{ex}}\parallel +\parallel {\phi }_{\mathrm{SN}}-{\phi }_{\mathrm{ex}}\parallel \\ \mathit{\hspace{0.17em}}& =(C_{\mathrm{SMM}}+C_{\mathrm{SN}})h^{\mathit{p}+1}.\end{array}$(110) predicts that the two solutions will converge at $\mathcal{O}(h^3)$. Using logarithmic regression, the experimentally observed order of convergence was $h_{\mathrm{m}\mathrm{a}\mathrm{x}}^{2.77}$. This result indicates that SMM does in fact converge to the ${\mathrm{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 $\mathit{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.

TABLE X Weak Scaling the Iterations Required by IP VEF and IP SMM Algorithms

Processors	Degrees of Freedom	Outer Iterations		Minimum Inner Iterations		Maximum Inner Iterations		Average Inner Iterations
		VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
1	4032	30	32	2	3	18	27	9.27	13.44
4	16128	51	55	2	3	22	32	10.45	16.13
16	64512	76	80	2	2	30	34	10.72	15.62
64	258048	135	135	1	2	35	37	10.88	15.64
256	1032192	261	277	1	1	50	37	10.59	14.78
1024	4128768	540	592	1	1	95	37	10.53	13.63

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 with 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

Processors	Degrees of Freedom	Average Sweep		Average Moment		Total
		VEF	SMM	VEF	SMM	VEF	SMM
1	4032	0.95	0.95	0.11	0.06	32.37	32.75
4	16128	1.12	1.06	0.12	0.07	64.27	62.74
16	64512	1.44	1.31	0.16	0.09	121.74	112.39
64	258048	1.62	1.48	0.20	0.11	247.57	215.81
256	1032192	1.64	1.48	0.21	0.12	484.85	445.19
1024	4128768	1.70	1.50	0.23	0.13	1048.52	972.07

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)(111) ${\epsilon }_n=\frac{\mathit{solve}\mathrm{time}\mathrm{with}one\mathrm{processor}}{\mathit{solve}\mathrm{time}\mathrm{with}n\mathrm{processors}}.$(111) denote the weak scaling efficiency where ideal weak scaling is characterized by ${\epsilon }_n=1$:

(111) ${\epsilon }_n=\frac{\mathit{solve}\mathrm{time}\mathrm{with}one\mathrm{processor}}{\mathit{solve}\mathrm{time}\mathrm{with}n\mathrm{processors}}.$(111)

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.

Fig. 13. The weak scaling efficiency for the average inner solve cost, average moment assembly cost, and total cost of the algorithm for the IP VEF and IP SMM algorithms. Scaling of the total algorithm is limited by the parallel block Jacobi transport sweep.

VI.F. Strong Scaling

Finally, we investigate strong scaling for the IP VEF and IP SMM algorithms with $\mathit{p}=2$ on the crooked pipe problem from Sec. VI.C. The problem size was fixed at 28 672 equally sized elements with S₁₂ 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 used 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.

TABLE XII Iterations to Convergence for a Strong Scaling Study of IP VEF and IP SMM

Processors	Outer Iterations		Minimum Inner Iterations		Maximum Inner Iterations		Average Inner Iterations
	VEF	SMM	VEF	SMM	VEF	SMM	VEF	SMM
1	47	53	1	1	18	31	8.26	12.60
2	57	67	1	1	25	31	9.37	13.96
4	68	77	1	2	33	33	10.96	14.99
8	76	84	1	2	44	35	11.70	15.64
16	83	86	1	2	42	36	11.07	15.98
32	112	113	1	2	40	36	11.51	16.28

The strong scaling speedup, defined equivalently to ${\epsilon }_n$ in Eq. (111)(111) ${\epsilon }_n=\frac{\mathit{solve}\mathrm{time}\mathrm{with}one\mathrm{processor}}{\mathit{solve}\mathrm{time}\mathrm{with}n\mathrm{processors}}.$(111) , is plotted in Fig. 14. Here, ideal speedup is ${\epsilon }_n=\mathit{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.

Fig. 14. Strong scaling speedup as a function of the number of processors on the crooked pipe problem with 28 672 elements, S₁₂ angular quadrature, and $\mathit{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.

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 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.^[8] and Olivier and Haut^[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 AMG-based block preconditioners for RT. The resulting discretizations and solvers are coupled to a high-order DG discretization of the ${\mathrm{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 $\mathcal{O}(h^{\mathit{p}+1})$ accuracy on refinements of a third-order 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-to-solution 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 ${\mathrm{S}}_N$ transport solution on a thick diffusion limit problem. Using a fixed angular quadrature rule, the SMM and ${\mathrm{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-to-solution 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.

For computational efficiency and software design purposes, we recommend the CG-based method since the CG SMM moment system had the fewest unknowns, was the cheapest to assemble and solve, and required the simplest discretization and preconditioning techniques. The IP- and RT-based SMMs are recommended for our intended multiphysics application in the hydrodynamics framework of Ref. [11] since both methods produce solutions in a compatible, discontinuous finite element space and are cheaper to solve and simpler to implement than the corresponding VEF-based algorithms. Furthermore, RT SMM has the distinct advantage that the left-hand-side diffusion operator is equivalent to the RT-based diffusion operator already implemented in Ref. [11], allowing the design of a transport algorithm that can reuse existing software.

We stress that suboptimal convergence of the current was observed for both RT VEF and RT SMM, and thus, it is not clear whether an RT-based moment method would truly improve physics compatibility with other mixed finite element–based multiphysics components. In the future, we plan to investigate whether the algorithmic flexibility seen here for the steady-state, linear transport problem extends to more complex and numerically demanding applications such as thermal radiative transfer and radiation hydrodynamics.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

One of the authors, S. O., was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, and the U.S. Department of Energy Computational Science Graduate Fellowship under award number DE-SC0019323 and by the U.S. Department of Energy through the Los Alamos National Laboratory. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (contract number 89233218CNA000001).