Nonlinear Elimination Applied to Radiation Diffusion

Abstract

We apply a nonlinearly preconditioned, quasi-Newton framework to accelerate the numerical solution of the thermal radiative transfer (TRT) equations. This framework was inspired by the unpublished method that has existed for years in Teton, Lawrence Livermore National Laboratory’s deterministic TRT code. In this paper, we cast this iteration scheme within a formal nonlinear preconditioning framework and compare its performance against other iteration schemes in the framework. With proper choices of iteration controls for the various levels of the solver, we can recover the standard linearized one-step method, a full nonlinear Newton scheme, as well as the method in Teton.

In brief, the nonlinear preconditioning TRT scheme formally eliminates the material temperature equation from the nonlinear system in a nonlinear analog of a Schur complement. This nonlinear elimination step involves solving a decoupled nonlinear equation for each spatial degree of freedom and is therefore inexpensive. By applying a quasi-Newton iteration scheme on the new system, we obtain a three-level iteration scheme that is at least as efficient as commonly used TRT schemes. The new method allows full convergence to the nonlinear backward Euler time-discretized system, increasing accuracy and robustness, while using a similar number of linear iterations as the more common linearized one-step methods Eq. (4(4) $J_{e,e}=\frac{\mathrm{\partial }m(e,\varphi )}{\mathrm{\partial }e},J_{e,\varphi }=\frac{\mathrm{\partial }m(e,\varphi )}{\mathrm{\partial }\varphi }$(4) ).

Keywords:

• Nonlinear accelerated inexact Newton method
• nonlinear elimination
• multigroup thermal radiation transport

I. INTRODUCTION

There have been many studies on applying nonlinear solvers to the thermal radiative transfer (TRT) equations.^1,2 Here, we discuss a new nonlinear preconditioning framework for TRT that is inspired by our discrete ordinates solver Teton.^3,4 In this work, we apply the nonlinearly preconditioned framework to the multigroup diffusion approximation of TRT (Refs. 5, 6, and 7). However, the method more generally applies to all forms of TRT, independent of spatial and angular discretizations 9.

Given the need for an implicit time-stepping scheme, a major algorithmic challenge is how to efficiently and robustly solve the resulting nonlinear equations for the updated temperature and radiation energy density. Typically, the emission term is linearized with respect to a current temperature approximation, which results in a linear transport equation for the updated radiation energy density that involves an effective scattering term.⁸ Solving this nonlinear system requires performing multiple matrix inversions (or transport sweeps) to fully converge the effective scattering term. This has given rise to the development of efficient and robust solvers.^8–11 Even with the use of efficient multigroup linear solvers, the high expense of taking a Newton step has led to the common practice of performing a single Newton step for each time step. Unfortunately, not fully converging the Newton iteration can result in unphysical violations of the maximum principle for realistic time steps; see for example Ref. 12 and Sec. V.

An alternative (unpublished) nonlinear iteration scheme was developed by Nowak for discrete ordinates TRT. Importantly, this scheme has approximately the same cost as performing a single Newton step but fully converges the nonlinear system. We recently learned that we could cast the Teton algorithm as a specific variant within the nonlinear elimination framework developed in Ref. 13. There are entire families of more advanced nonlinear preconditioners^14,15 that will be the focus of future work. In particular, we use the formalism in Ref. 13 to derive a three-level iteration scheme that reduces to Nowak’s nonlinear elimination scheme, as well as the more traditional one-step Newton scheme described above, by an appropriate choice of iteration parameters within the iteration levels. We will explore these two variants, along with three more in Sec. V.

We will start with an introduction to nonlinear elimination in Sec. I.A, followed by a review of the thermal radiation diffusion equations and the nonlinear system formulation in Sec. II. We then proceed with a review of a method that closely follows previous methods^8,11,16 in Sec. III but cast it more explicitly in terms of a standard Newton iteration in order to apply the nonlinear elimination more easily. The main new algorithm where we apply nonlinear elimination to the multigroup radiation diffusion system is described in Sec. IV. We show some numerical results for two test problems in Sec. V, followed by some conclusions in Sec. VI.

I.A. Review of Nonlinear Elimination

We start with a general review of the Newton method and then review the nonlinear elimination method of Lanzkron et al.¹³ before introducing the TRT equations. Assume we have a coupled set of nonlinear equations we would like to solve, namely,

(1) $f(e,\varphi )=\left[\begin{array}{c}m(e,\varphi )\\ r(e,\varphi )\end{array}\right]=0,$(1)

where $e$ and $\varphi$ are the independent variables and $m(e,\varphi )$ and $r(e,\varphi )$ are two equations that describe the system.^a Within a Newton iteration scheme, a key step is to solve the block linear system,

for the updates ${\delta }_e$ and ${\delta }_{\varphi }$ and then compute the new state using

(3) $\left[\begin{array}{c}{e}^{l+1}\\ {\varphi }^{l+1}\end{array}\right]=\left[\begin{array}{c}{e}^l\\ {\varphi }^l\end{array}\right]+\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right],$(3)

where the superscript $l$ indicates quantities are evaluated at the last iteration and the Jacobian matrices are given by

(4) $J_{e,e}=\frac{\mathrm{\partial }m(e,\varphi )}{\mathrm{\partial }e},J_{e,\varphi }=\frac{\mathrm{\partial }m(e,\varphi )}{\mathrm{\partial }\varphi }$(4)

and

(5) $J_{\varphi ,e}=\frac{\mathrm{\partial }r(e,\varphi )}{\mathrm{\partial }e},J_{\varphi ,\varphi }=\frac{\mathrm{\partial }r(e,\varphi )}{\mathrm{\partial }\varphi }.$(5)

Note that the first subscript corresponds to one of the equations, and the second subscript indicates the derivative.

If one of the block matrices, $J_{e,e},$ for example, is easy to invert, we can use a Schur complement to reduce the block linear system with two variables into a smaller system with just one unknown, namely,

(6) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }=-r(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(6)

While this matrix is more complex, the reduction in dimension can significantly speed up the method. In fact, this is the basis of most grey and multigroup solvers; see, for example, Refs. 8, 10, and 11. After ${\delta }_{\varphi }$ has been computed, ${\delta }_e$ can be recovered using

(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7)

It stands to reason that if one of the blocks of the coupled system is easy to invert at the linear step, it may also be relatively cheap to solve the related nonlinear problem. Nonlinear elimination, as presented in Ref. 13, can be thought of as the nonlinear analog of a Schur complement, eliminating an unknown from the nonlinear system. It does this by posing a different nonlinear system to replace Eq. (1(1) $f(e,\varphi )=\left[\begin{array}{c}m(e,\varphi )\\ r(e,\varphi )\end{array}\right]=0,$(1) ), namely,

(8) $\hat{f}(\varphi )=r(\hat{e}(\varphi ),\varphi )=0,$(8)

where $\hat{e}(\varphi )$ returns the solution $m(\hat{e},\varphi )=0$ for a given $\varphi$. Lanzkron et al.¹³ define the formal steps necessary in order to ensure that the solution of $\hat{f}(\varphi )=0$ is the same as $f(e,\varphi )=0$, and we direct the reader there for the details. In this work, we outline the implications for our multigroup radiation diffusion implementation. In order to solve Eq. (8(8) $\hat{f}(\varphi )=r(\hat{e}(\varphi ),\varphi )=0,$(8) ) using a Newton method, we introduce an auxiliary variable ${\delta }_e$ and apply the chain rule to get

(9) $\left[\begin{array}{cc}{J}_{e,e}(\hat{e},{\varphi }^l)& J_{e,\varphi }(\hat{e},{\varphi }^l)\\ J_{\varphi ,e}(\hat{e},{\varphi }^l)& J_{\varphi ,\varphi }(\hat{e},{\varphi }^l)\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}0\\ r(\hat{e},{\varphi }^l)\end{array}\right].$(9)

This is remarkably similar to) except that the Jacobian matrix and the residual are evaluated at an advanced state, and one of the residuals is exactly zero.

Another key principle of nonlinear elimination to keep in mind is that the variable we eliminate is no longer an independent variable in the algorithm; it is strictly a function of the variable that was retained. Here, for example, anytime $\varphi$ is updated, we must solve $m(\hat{e},{\varphi }^l)=0$ for the updated $\hat{e}$. Practically, this means there is one solve for $\hat{e}$ before the outer Newton iterations begin as well as one per Newton step after $\varphi$ has been updated. This takes the place of computing the update ${\delta }_e$ using the linearized Eq. (7)(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7) . The full algorithm will be shown in more detail later within the context of our multigroup radiation solver 9.

II. THE EQUATIONS

The multigroup radiation diffusion equations⁵ coupled to a material energy equation are

(10) $\rho \frac{\mathrm{\partial }e}{\mathrm{\partial }t}=-c\sum _g{\sigma }_{a,g}{B}_g+c\sum _g{\sigma }_{a,g}{\varphi }_g+Q$(10)

and

(11) $\begin{array}{r}\frac{\mathrm{\partial }{\varphi }_g}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g+c{\sigma }_{a,g}{B}_g-c{\sigma }_{a,g}{\varphi }_g+S_g\end{array},$(11)

with a boundary condition of

(12) $c{\alpha }_g{\varphi }_g+{\beta }_g\frac{c}{3{\sigma }_{t,g}}\hat{\mathbf{n}}\cdot \mathrm{\nabla }{\varphi }_g={\gamma }_g,$(12)

where

$e$ =	=	material specific internal energy
$\rho$ =	=	density
$t$ =	=	time
$c$ =	=	speed of light
${\sigma }_{a,g}$ =	=	absorption opacity in group $g$
${\sigma }_{t,g}$ =	=	total (absorption plus scattering) group opacity
$B_g$ =	=	Planck function integrated17 over group $g$
${\varphi }_g$ =	=	radiation energy density in group $g$
$Q$ =	=	arbitrary material source
$S_g$ =	=	arbitrary radiation source for group $g$
${\alpha }_g$, ${\beta }_g$, ${\gamma }_g$ =	=	coefficients that determine the type of boundary condition.

Most of the coefficients can depend on $e$, $\rho$, and ${\varphi }_g$ in the most general case, but here we will only consider their dependence on $e$ and $\rho$ through the material temperature $T(e,\rho )$. In particular, we have these dependencies: ${\sigma }_{a,g}(\rho ,T)$, ${\sigma }_{t,g}(\rho ,T)$, and $B_g(T)$.

There are several simplifying assumptions in the system we make in this paper. First, without coupling to hydrodynamics, the density is constant, so further time dependence on the density will be ignored. Additionally, we will ignore physical scattering such as Compton scattering,⁵ which couples the radiation groups and material equation through an additional inelastic scattering process. This kind of scattering can be incorporated into this method but has been omitted here for simplicity Eq. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ).

II.A. Nonlinear Function Formulation

If we discretize Eqs. (10(10) $\rho \frac{\mathrm{\partial }e}{\mathrm{\partial }t}=-c\sum _g{\sigma }_{a,g}{B}_g+c\sum _g{\sigma }_{a,g}{\varphi }_g+Q$(10) ) and (11(11) $\begin{array}{r}\frac{\mathrm{\partial }{\varphi }_g}{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g+c{\sigma }_{a,g}{B}_g-c{\sigma }_{a,g}{\varphi }_g+S_g\end{array},$(11) ) in time using the backward Euler method, the nonlinear function we aim to solve is

(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13)

It is useful to define some notation, namely,

(14) $\varphi =\left[\begin{array}{c}{\text{pphi}}_1\cdots {\text{pphi}}_G\end{array}\right],$(14)

(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15)

(16) $r_g(\mathrm{e},{\varphi }_{\mathrm{g}})=\frac{1}{\mathrm{\Delta }t}\left({\varphi }_{\mathrm{g}}-{\varphi }_{\mathrm{g}}^{\mathrm{n}-1}\right)-\mathrm{\nabla }\cdot \frac{\mathrm{c}}{3{\sigma }_{\mathrm{t},\mathrm{g}}}\mathrm{\nabla }{\varphi }_{\mathrm{g}}-\mathrm{c}{\sigma }_{\mathrm{a},\mathrm{g}}{B}_{\mathrm{g}}+\mathrm{c}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-{\mathrm{S}}_{\mathrm{g}},$(16)

and

(17) $\mathbf{r}(e,\varphi )={\left[\begin{array}{c}{r}_1(e,{\varphi }_1)\cdots r_G(e,{\varphi }_G)\end{array}\right]}^{\mathrm{T}},$(17)

where

$\varphi$ =	=	vector of unknowns for all the group radiation energy densities
$m(e,\varphi )$ =	=	residual of material energy conservation equation
$r_g(e,{\varphi }_g)$ =	=	residual of the radiation diffusion equation for group $g$.

II.B. Block Jacobian Operators

The block Jacobian operator for Eq. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ) is

(18) $J(e,\varphi )=\left[\begin{array}{cc}\left[\begin{array}{c}{J}_{e,e}\end{array}\right]& \left[\begin{array}{ccc}{J}_{e,{\varphi }_1}& \cdots & J_{e,{\varphi }_G}\end{array}\right]\\ \left[\begin{array}{c}{J}_{{\varphi }_1,e}\\ ⋮\\ J_{{\varphi }_G,e}\end{array}\right]& \left[\begin{array}{ccc}{J}_{{\varphi }_1,{\varphi }_1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& J_{{\varphi }_G,{\varphi }_G}\end{array}\right]\end{array}\right]=\left[\begin{array}{cc}{J}_{e,e}& J_{e,\varphi }\\ J_{\varphi ,e}& J_{\varphi ,\varphi }\end{array}\right],$(18)

where

(19) $J_{e,e}=\frac{\rho }{\mathrm{\Delta }t}+\frac{c}{c_v}\sum _g{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T},J_{e,{\varphi }_g}=-c{\sigma }_{a,g},$(19)

(20) $\begin{array}{rl}& J_{{\varphi }_g,e}=-\frac{c}{c_v}{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T},\mathrm{a}\mathrm{n}\mathrm{d}\\ & J_{{\varphi }_g,{\varphi }_g}=\frac{1}{\mathrm{\Delta }t}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }+c{\sigma }_{a,g},\end{array}$(20)

and where $c_v=\frac{\mathrm{\partial }e}{\mathrm{\partial }T}.$ We have ignored opacity derivatives, but they could easily be incorporated. In fact, we will continue the use of the notation for the Jacobian operators throughout the paper, leaving them abstract, which has several advantages. First, derivatives of the opacity with respect to the radiation energy density or material energy are easy to add and propagate through the entire algorithm. The second advantage is that in our implementation, each of these terms is treated discretely, and we treat them as matrix operators. Third, it makes the notation much more compact. Note that only the Jacobian $J_{\varphi ,\varphi }$ contains spatial derivatives and is therefore more expensive to invert; the $J_{e,e}$ Jacobian matrix only has local coupling and is therefore easy to invert.

III. NEWTON ITERATION FORMULATION

Given the similarities between a standard Newton iteration step, given by), and a Newton step with nonlinear elimination applied, given by Eq. (9(9) $\left[\begin{array}{cc}{J}_{e,e}(\hat{e},{\varphi }^l)& J_{e,\varphi }(\hat{e},{\varphi }^l)\\ J_{\varphi ,e}(\hat{e},{\varphi }^l)& J_{\varphi ,\varphi }(\hat{e},{\varphi }^l)\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}0\\ r(\hat{e},{\varphi }^l)\end{array}\right].$(9) ), we will do a detailed review of the derivation for solving Eq. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ) efficiently by closely following earlier work.^8–11 Here, we cast the previous work more explicitly as operations within a standard Newton iteration so that it is easier to see the modifications that applying nonlinear elimination makes to the iteration scheme, shown in Sec. IV.

Inserting Eqs. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ) through (18(18) $J(e,\varphi )=\left[\begin{array}{cc}\left[\begin{array}{c}{J}_{e,e}\end{array}\right]& \left[\begin{array}{ccc}{J}_{e,{\varphi }_1}& \cdots & J_{e,{\varphi }_G}\end{array}\right]\\ \left[\begin{array}{c}{J}_{{\varphi }_1,e}\\ ⋮\\ J_{{\varphi }_G,e}\end{array}\right]& \left[\begin{array}{ccc}{J}_{{\varphi }_1,{\varphi }_1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& J_{{\varphi }_G,{\varphi }_G}\end{array}\right]\end{array}\right]=\left[\begin{array}{cc}{J}_{e,e}& J_{e,\varphi }\\ J_{\varphi ,e}& J_{\varphi ,\varphi }\end{array}\right],$(18) ) into, we get

(21) $\left[\begin{array}{cc}{J}_{e,e}& J_{e,\varphi }\\ J_{\varphi ,e}& J_{\varphi ,\varphi }\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}m(e^l,{\varphi }^l)\\ \mathbf{r}(e^l,{\varphi }^l)\end{array}\right].$(21)

We first eliminate the material energy variable from the linear system, which is easy to do since there are no spatial operators in the material equation. We solve for ${\delta }_e$ symbolically in Eq. (21(21) $\left[\begin{array}{cc}{J}_{e,e}& J_{e,\varphi }\\ J_{\varphi ,e}& J_{\varphi ,\varphi }\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}m(e^l,{\varphi }^l)\\ \mathbf{r}(e^l,{\varphi }^l)\end{array}\right].$(21) ) using a Schur complement to get

(22) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }{\delta }_{\varphi }^l\right],$(22)

which is the application of the generic Eq. (7(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7) ) to our multigroup TRT system. We can then insert Eq. (22(22) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }{\delta }_{\varphi }^l\right],$(22) ) into the second row of Eq. (21(21) $\left[\begin{array}{cc}{J}_{e,e}& J_{e,\varphi }\\ J_{\varphi ,e}& J_{\varphi ,\varphi }\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}m(e^l,{\varphi }^l)\\ \mathbf{r}(e^l,{\varphi }^l)\end{array}\right].$(21) ) to get an equation for ${\delta }_{\mathbf{f}}^l$, namely,

(23) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }^l=-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(23)

To connect this compact notation with previous work, we take just one row of this matrix for group $g$,

(24) $J_{{\varphi }_g,{\varphi }_g}{\delta }_{{\varphi }_g^l}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\delta }_{{\varphi }_{g^{\prime }}^l}=-r_g(e^l,{\varphi }_g^l)+J_{{\varphi }_g,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l),$(24)

and then insert the definitions of ${\delta }_{{\varphi }_g}^l={\varphi }_g^{l+1}-{\varphi }_g^l$, $m(e,\varphi )$, and $r_g(e,{\varphi }_g)$, as well as $J_{{\varphi }_g,{\varphi }_g}$ and $J_{e,{\varphi }_g}$, to get

(25) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ & =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}-Q\right).\end{array}$(25)

All properties, $T$, $B_g$, ${\sigma }_{a,g}$, and ${\sigma }_{t,g},$ are evaluated using $e^l$ unless otherwise annotated. This is the same as Eqs. (6(6) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }=-r(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(6) ) and (7(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7) ) in Ref. 8. The rectangular operator $J_{e,\mathbf{f}}$ collapses a change in the multigroup energy density space into an effect on the material energy equation, and the operator $J_{\mathbf{f},e}{J}_{e,e}^{-1}$ maps changes from the material energy equation to the radiation diffusion equation space. This factor, $J_{\mathbf{f},e}{J}_{e,e}^{-1}$, is the same as the so-called Fleck factor in Ref. 10 and the product of Eqs. (7c(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7) ) and (7d(7) ${\delta }_e^l=-J_{e,e}^{-1}\left[m(e^l,{\varphi }^l)+J_{e,\varphi }\left({\varphi }^{l+1}-{\varphi }^l\right)\right].$(7) ) in Ref. 8.

III.A. Grey Acceleration

Equation (23)(23) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }^l=-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(23) can be expensive to solve because it is large (coupling all groups together) and can be ill-conditioned for large time steps or large opacities. Without physical scattering, each of the radiation groups is only coupled via the absorption and reemission through the material; any information about the radiation spectrum is lost once it is converted into material energy. This fact is exploited in many numerical methods to develop efficient solvers.^8,10 Inserting ${\delta }_{\varphi }^l={\varphi }^{l+1}-{\varphi }^l$ into Eq. (23(23) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }^l=-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(23) ), we get

(26) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^{l+1}=\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^l-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(26)

Following Refs. 8 and 11, we split the solution into two components,

(27) ${\varphi }^{l+1}={\varphi }^{\star }+{\varphi }^{†},$(27)

where we choose ${\varphi }^{\star }$ to be the solution to a system where the group-to-group scattering term is lagged and therefore easier to solve, namely,

(28) $J_{\varphi ,\varphi }{\varphi }^{\star }=J_{\varphi ,\varphi }{\varphi }^{†}-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(28)

For one group, Eq. (28(28) $J_{\varphi ,\varphi }{\varphi }^{\star }=J_{\varphi ,\varphi }{\varphi }^{†}-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(28) ) expands and simplifies to

(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29)

If we then insert ${\varphi }^{l+1}={\varphi }^{\star }+{\varphi }^{†}$ into Eq. (26(26) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^{l+1}=\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^l-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(26) ) and subtract Eq. (28(28) $J_{\varphi ,\varphi }{\varphi }^{\star }=J_{\varphi ,\varphi }{\varphi }^{†}-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(28) ), we get

(30) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^{†}=J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\left({\varphi }^{\star }-{\varphi }^l\right),$(30)

or for just one group, we get

(31) $J_{{\varphi }_g,{\varphi }_g}{\varphi }_g^{†}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varphi }_{g^{\prime }}^{†}=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right).$(31)

This is no easier to solve than our original problem Eq. (26(26) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^{l+1}=\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\varphi }^l-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(26) ). Physically, ${\varphi }^{†}$ represents the photons that have been absorbed and then reemitted by the material within this Newton step.^b Inspecting the right side of Eq. (31(31) $J_{{\varphi }_g,{\varphi }_g}{\varphi }_g^{†}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varphi }_{g^{\prime }}^{†}=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right).$(31) ), the sum over groups eliminates the detailed information of the spectral error (groupwise error) that was made, and the reemitted photons represented by $\varphi$ have the same spectrum, or relative magnitude between the groups, for any source term. This suggests that if we could approximate the spectral shape somehow, it would be possible to average the multigroup Eq. (31(31) $J_{{\varphi }_g,{\varphi }_g}{\varphi }_g^{†}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varphi }_{g^{\prime }}^{†}=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right).$(31) ) for one unknown. In fact, this has been the study of many papers,^8,11 and we continue to closely follow those derivations Eq. (25(25) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ & =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}-Q\right).\end{array}$(25) ).

We can define a group sum of the correction term, ${\text{pphi}}^{†}={\sum }_g{\text{pphi}}_g^{†}$, and a prolongation operator that recovers the group values, namely,

(32) ${\varphi }_g^{†}={\varsigma }_g{\varphi }^{†}.$(32)

Formally, if we insert Eq. (32(32) ${\varphi }_g^{†}={\varsigma }_g{\varphi }^{†}.$(32) ) into Eq. (31(31) $J_{{\varphi }_g,{\varphi }_g}{\varphi }_g^{†}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varphi }_{g^{\prime }}^{†}=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right).$(31) ) and sum over groups, we get

(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33)

Some approximation so that we can easily estimate ${\varsigma }_g$ is needed in order to make using Eq. (33(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ) practical. Morel et al.⁸ use an infinite medium solution to Eq. (31(31) $J_{{\varphi }_g,{\varphi }_g}{\varphi }_g^{†}-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varphi }_{g^{\prime }}^{†}=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right).$(31) ) to compute

(34) ${\varsigma }_g\approx \frac{\frac{{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T}}{\frac{1}{\mathrm{\Delta }t}+c{\sigma }_{a,g}}}{\sum _{g^{\prime }}\frac{{\sigma }_{a,g^{\prime }}\frac{\mathrm{\partial }{B}_{g^{\prime }}}{\mathrm{\partial }T}}{\frac{1}{\mathrm{\Delta }t}+c{\sigma }_{a,g^{\prime }}}}.$(34)

Because we form discrete Jacobian matrices from Eqs. (19(19) $J_{e,e}=\frac{\rho }{\mathrm{\Delta }t}+\frac{c}{c_v}\sum _g{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T},J_{e,{\varphi }_g}=-c{\sigma }_{a,g},$(19) ) and (20(20) $\begin{array}{rl}& J_{{\varphi }_g,e}=-\frac{c}{c_v}{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T},\mathrm{a}\mathrm{n}\mathrm{d}\\ & J_{{\varphi }_g,{\varphi }_g}=\frac{1}{\mathrm{\Delta }t}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }+c{\sigma }_{a,g},\end{array}$(20) ) and perform all operations on the block matrices, it was difficult to interpret how to apply Eq. (34(34) ${\varsigma }_g\approx \frac{\frac{{\sigma }_{a,g}\frac{\mathrm{\partial }{B}_g}{\mathrm{\partial }T}}{\frac{1}{\mathrm{\Delta }t}+c{\sigma }_{a,g}}}{\sum _{g^{\prime }}\frac{{\sigma }_{a,g^{\prime }}\frac{\mathrm{\partial }{B}_{g^{\prime }}}{\mathrm{\partial }T}}{\frac{1}{\mathrm{\Delta }t}+c{\sigma }_{a,g^{\prime }}}}.$(34) ) in this context. Instead, we first ignore the group-to-group coupling and solve

(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35)

for ${\hat{\varphi }}_g$. We can think of ${\varphi }_g^{\star }$ as the portion of photons of the final solution that have not been absorbed and reemitted, which is often called the uncollided flux. In that light, ${\hat{\varphi }}_g$ can be thought of as the first-collided photons, whose spectrum should be the same as subsequent generations. We then compute the scaling pointwise by computing

(36) ${\varsigma }_g\approx \frac{{\hat{\varphi }}_g}{\sum _{g^{\prime }}{\hat{\varphi }}_{g^{\prime }}}.$(36)

Each time a grey solve of Eq. (33(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ) for ${\varphi }^{†}$ is needed in the algorithm, we first solve Eq. (35(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35) ), compute a spectrum using Eq. (36(36) ${\varsigma }_g\approx \frac{{\hat{\varphi }}_g}{\sum _{g^{\prime }}{\hat{\varphi }}_{g^{\prime }}}.$(36) ), and then compute a new group-collapsed matrix to use when solving Eq. (33(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ).

Once ${\varphi }^{\star }$ and ${\varphi }^{†}$ have been computed, we compute ${\varphi }^{l+1}$ using

(37) ${\varphi }_g^{l+1}={\varphi }_g^{\star }+{\varsigma }_g{\varphi }^{†}.$(37)

With our final estimate of ${\varphi }_g^{l+1}$, we can then compute the material update using

(38) $\begin{array}{rl}& e^{l+1}=e^l-J_{e,e}^{-1}\\ & \left[\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g^{l+1}-Q\right].\end{array}$(38)

III.B. The Nonlinear Multifrequency Grey Algorithm

We now assemble the details presented in this section into the full iteration scheme for the standard Nonlinear Multifrequency Grey (NLMFG) iteration, shown in Algorithm 1. This is similar to Algorithm 1 in Ref. 8, but we also include the obvious extension of adding an outer loop to take multiple Newton steps, instead of just one, in order to fully converge the nonlinear system Eq. (40(40) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l.\end{array}$(40) ).

Algorithm 1: NLMFG: Nonlinear Multifrequency-Grey iteration

$l=0$;

set ${\varphi }^l={\varphi }^{n-1}$;

Outer nonlinear iteration to converge

repeat

Multigroup linear iteration to converge

repeat

 solve for ${\varphi }_g^{\star }$ using Eq. (29(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) );

 compute ${\varsigma }_g$ by solving Eqs. (35(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35) –36)(36) ${\varsigma }_g\approx \frac{{\hat{\varphi }}_g}{\sum _{g^{\prime }}{\hat{\varphi }}_{g^{\prime }}}.$(36) ;

 solve for a grey correction ${\varphi }^{†}$ using Eq. (33)(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ;

 compute the ${\varphi }_g^{l+1}={\varphi }_g^{\star }+{\varsigma }_g{\varphi }^{†}$;

until;

solve for an updated $e^{l+1}$ by using Eq. (38)(38) $\begin{array}{rl}& e^{l+1}=e^l-J_{e,e}^{-1}\\ & \left[\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g^{l+1}-Q\right].\end{array}$(38) ;

$l=l+1$;

until;

Variants of this algorithm are possible, such as not fully converging either the outer or the inner loops. These variants will be explained and their efficacy explored in Sec. V.

IV. NONLINEAR ELIMINATION APPLIED TO MULTIGROUP RADIATION DIFFUSION

We will now follow the outline in Sec. I.A and apply it specifically to our multigroup radiation system in Eq. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ). This section is the main new work of this paper. Because of the similarity between the Newton step needed to solve the full system $f(e,\varphi )$ using and the reduced system $\hat{f}(\varphi )=f(\hat{e}(\varphi ),\varphi )$ using Eq. (9(9) $\left[\begin{array}{cc}{J}_{e,e}(\hat{e},{\varphi }^l)& J_{e,\varphi }(\hat{e},{\varphi }^l)\\ J_{\varphi ,e}(\hat{e},{\varphi }^l)& J_{\varphi ,\varphi }(\hat{e},{\varphi }^l)\end{array}\right]\left[\begin{array}{c}{\delta }_e\\ {\delta }_{\varphi }\end{array}\right]=-\left[\begin{array}{c}0\\ r(\hat{e},{\varphi }^l)\end{array}\right].$(9) ), we use the results from Sec. III but exploit the fact that $m(\hat{e},\varphi )=0$ to simplify several of the equations. Doing this also leads to some insights into how the method works.

The full linearized multigroup system in Eq. (23(23) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }^l=-\mathbf{r}(e^l,{\varphi }^l)+J_{\varphi ,e}{J}_{e,e}^{-1}m(e^l,{\varphi }^l).$(23) ) becomes

(39) $\left(J_{\varphi ,\varphi }-J_{\varphi ,e}{J}_{e,e}^{-1}{J}_{e,\varphi }\right){\delta }_{\varphi }^l=-\mathbf{r}(e^l,{\varphi }^l),$(39)

or equivalently, the expanded Eq. (25(25) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ & =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}-Q\right).\end{array}$(25) ) becomes

(40) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l.\end{array}$(40)

Note that only the absorption term is corrected for the reemission in this new method. The acceleration of the nonlinear elimination removes the need to scale the emission in the Newton step. As we converge, ${\varphi }_g^{l+1}\approx {\varphi }_g^l$, and the pseudoscattering terms vanish from Eq. (40(40) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l.\end{array}$(40) ).

The multigroup solve for ${\varphi }_g^{\star }$ in Eq. (29)(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) reduces to

(41) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g\end{array}.$(41)

Moving the work of the pseudoscattering to the grey equation in addition to performing the nonlinear elimination has transformed Eq. (29(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) ) back into a straightforward evaluation of the original Eq. (16(16) $r_g(\mathrm{e},{\varphi }_{\mathrm{g}})=\frac{1}{\mathrm{\Delta }t}\left({\varphi }_{\mathrm{g}}-{\varphi }_{\mathrm{g}}^{\mathrm{n}-1}\right)-\mathrm{\nabla }\cdot \frac{\mathrm{c}}{3{\sigma }_{\mathrm{t},\mathrm{g}}}\mathrm{\nabla }{\varphi }_{\mathrm{g}}-\mathrm{c}{\sigma }_{\mathrm{a},\mathrm{g}}{B}_{\mathrm{g}}+\mathrm{c}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-{\mathrm{S}}_{\mathrm{g}},$(16) ). The grey acceleration step Eq. (33(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ) does not change since it is only correcting for lagging the effective scattering term, which is the same in both cases.

The final Nonlinear Multifrequency Grey iteration with Nonlinear Elimination (NLMFG-NE) method is outlined in Algorithm 2 and is remarkably similar to Algorithm 1. The key difference is that there is a nonlinear solve for the material energy equation initially and at the end of each Newton step instead of the update to the material energy using Eq. (38(38) $\begin{array}{rl}& e^{l+1}=e^l-J_{e,e}^{-1}\\ & \left[\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g^{l+1}-Q\right].\end{array}$(38) ).

Algorithm 2: NLMFG-NE: Nonlinear Multifrequency Grey iteration with Nonlinear Elimination

$l=0$;

set ${\varphi }^l={\varphi }^{n-1}$;

solve for $e^l$ using $m(e^l,{\varphi }^l)=0$ using Eq. (15)(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) ;

Outer nonlinear iteration to converge

repeat

Multigroup linear iteration to converge

repeat

 solve for ${\varphi }_g^{\star }$ using Eq. (41(41) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g\end{array}.$(41) ) (or Eq. 29)(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) ;

 compute ${\varsigma }_g$ by solving Eqs. (35(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35) –36)(36) ${\varsigma }_g\approx \frac{{\hat{\varphi }}_g}{\sum _{g^{\prime }}{\hat{\varphi }}_{g^{\prime }}}.$(36) ;

 solve for a grey correction ${\varphi }^{†}$ using Eq. (33)(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ;

 compute the ${\varphi }_g^{l+1}={\varphi }_g^{\star }+{\varsigma }_g{\varphi }^{†}$;

until $|{\sum }_g{\varphi }_g^{l+1}-{\varphi }_g^l|<{\epsilon }_G|{\sum }_g{\varphi }_g^{l+1}|$;

solve Eq. (15)(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) for an updated $e^{l+1}$;

$l=l+1$;

until $|e^{l+1}-e^l|<{\epsilon }_e|e^{l+1}|$;

Algorithm 1 could almost be written in terms of Algorithm 2 with the proper choices of iteration controls when solving for $\hat{e}$ in $m(\hat{e},\varphi )=0$. In fact, Eq. (38(38) $\begin{array}{rl}& e^{l+1}=e^l-J_{e,e}^{-1}\\ & \left[\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g^{l+1}-Q\right].\end{array}$(38) ) is just one iteration of a Newton method to solve Eq. (15(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) ). The main difference that cannot be resolved by a clever setting of tolerances is that in the nonlinear elimination method, there is one nonlinear solve for $\hat{e}$ using Eq. (15(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) ) before the outer Newton iterations begin, while in Algorithm 1 there is no such application of Eq. (38(38) $\begin{array}{rl}& e^{l+1}=e^l-J_{e,e}^{-1}\\ & \left[\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g^{l+1}-Q\right].\end{array}$(38) ) before the main loop.

Also note that either Eq. (25(25) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ & =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}-Q\right).\end{array}$(25) ) or Eq. (40(40) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l.\end{array}$(40) ) along with Eq. (29(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) ) or Eq. (41(41) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g\end{array}.$(41) ) can be used in this iteration scheme. If you already had a method to solve Eq. (13(13) $\begin{array}{r}f\left(\left[\begin{array}{c}e\\ {\varphi }_1\\ ⋮\\ {\varphi }_G\end{array}\right]\right)=\left[\begin{array}{c}\frac{\rho }{\mathrm{\Delta }t}\left(e-e^{n-1}\right)+c\sum _g{\sigma }_{a,g}{B}_g-c\sum _g{\sigma }_{a,g}{\varphi }_g-Q\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_1-{\varphi }_1^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,1}}\mathrm{\nabla }{\varphi }_1-c{\sigma }_{a,1}{B}_1+c{\sigma }_{a,1}{\varphi }_1-S_1\\ ⋮\\ \frac{1}{\mathrm{\Delta }t}\left({\varphi }_G-{\varphi }_G^{n-1}\right)-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,G}}\mathrm{\nabla }{\varphi }_G-c{\sigma }_{a,G}{B}_G+c{\sigma }_{a,G}{\varphi }_G-S_G\end{array}\right]=0\end{array}.$(13) ) using a standard Newton iteration implemented with a variant of Algorithm 1, it can be accelerated easily with this method. In practice, we use Eqs. (25(25) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{l+1}-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{l+1}+c{\sigma }_{a,g}{\varphi }_g^{l+1}+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^{l+1}\\ & =\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}-Q\right).\end{array}$(25) ) and (29(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) ). In addition to simplifying our implementation, we have also noticed more robustness given that the nonlinear elimination is only applied approximately so that $m(\hat{e},\varphi )\approx 0$.

In Sec. III.A, we saw how the grey acceleration is used to accelerate the convergence of the absorption and reemission of photons within the Newton step. This absorption and reemission appear in the linear Newton step and can be applied in a linear manner, but the grey acceleration does nothing to accelerate the emission and reabsorption. The nonlinear elimination is exactly targeted at the emission and reabsorption component of the Newton iteration. But, since the emission term is fundamentally nonlinear through the Planck function, it requires a nonlinear method in order to address it. This implies that if nonlinear elimination is used, we must fully converge the nonlinear problem; a one-step Newton method is insufficient. One of the major reasons for this is that energy is only conserved between the radiation and material equations at full convergence of the nonlinear system. The emission term $c{\sigma }_{a,g}{B}_g$ is evaluated using two different states in Eq. (41(41) $\begin{array}{c}\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g\end{array}.$(41) ) and after solving Eq. (15(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) ) for the updated $\hat{e}$ in the nonlinear acceleration step; upon convergence in the outer Newton iterations, these become the same.

V. RESULTS

We will now test the different algorithm variations on two test problems. These problems are based on the radiation flow down a pipe used by Morel et al.⁸ The first variant is a zero-dimensional problem where the radiation and material start out of equilibrium. This problem allows us to study the temporal convergence of the different algorithms. The second problem is a radiation flow problem to show that the method works in multiple dimensions as well.

The methods outlined in Algorithm 1 and 2 allow for various choices of how many iterations are performed for each loop. Fully converging the outer loop gives the exact solution to the nonlinear time and space discretized system. Alternatively, just one outer iteration can be performed. If the time steps are small enough, this can be fast, but there is a stability constraint on this approximate solution to the nonlinear equations that can be violated. The linear system at each Newton iteration can be fully converged, or it can be approximately solved with just one multigroup iteration in an inexact Newton method. The final option is whether or not to apply the nonlinear elimination.

Not all combinations of these options yield desirable methods. Without nonlinear elimination, in order to conserve energy, either the outer Newton or inner linearized multigroup system must be fully converged. With nonlinear elimination, the outer Newton iterations must be converged in order to conserve energy. Table I lists the variations we studied and provides shorthand monikers as well as the solver tolerance values for each level of the algorithm. The general spirit behind the choices of tolerances is that the innermost loops are converged the tightest and that the relative tolerances get looser with each level. For the problems run here, fairly tight tolerances were chosen. Not shown in Table I is the tolerance on the inner nonlinear solves for $m(\hat{e},\varphi )=0$, which is set to ${\epsilon }_{NE}={10}^{-12}$. A tolerance of infinity means exactly one step was taken. For each diffusion matrix solve, we used the sparse direct linear solver UMFPACK (Ref. 18), so there is no tolerance control on each linear solve.

TABLE I Monikers for Algorithm Variations and Tolerances

Name	Outer Newton, ${\epsilon }_e$	Inner Multigroup, ${\epsilon }_G$	Nonlinear Elimination
LMFG: Linear Multifrequency Grey	$\mathrm{\infty }$	${10}^{-11}$	Algorithm 1
NLMFG: Nonlinear Multifrequency Grey	${10}^{-10}$	${10}^{-11}$	Algorithm 1
NLMFG-NE: NLMFG iteration with Nonlinear Elimination	${10}^{-10}$	${10}^{-11}$	Algorithm 2
INMFG: Inexact Newton Multifrequency Grey	${10}^{-10}$	$\mathrm{\infty }$	Algorithm 1
INMFG-NE: INMFG with Nonlinear Elimination	${10}^{-10}$	$\mathrm{\infty }$	Algorithm 2

For any given spatial mesh and time step size, there is some discretization error. The Linear Multifrequency Grey (LMFG) method injects additional error relative to any of the other methods since it does not attempt to fully converge the nonlinear system. We will see that this leads to the LMFG method needing smaller time steps than the fully converged methods to achieve the same solution quality.

We will measure the expense of each method in terms of the number of linear iterations it performs on the multigroup system in order to solve Eqs. (29(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) ) and (35(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35) ), and Eq. (33(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33) ) once. For our ten group problems here, this means 21 mesh-sized diffusion matrices to solve. These methods trade doing more inner linear iterations for doing more outer iterations. We assume the other work, including solving Eq. (15)(15) $m(\mathrm{e},\varphi )=\frac{\rho }{\mathrm{\Delta }\mathrm{t}}\left(\mathrm{e}-{\mathrm{e}}^{\mathrm{n}-1}\right)+\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\mathrm{B}}_{\mathrm{g}}-\mathrm{c}\sum _{\mathrm{g}}{\sigma }_{\mathrm{a},\mathrm{g}}{\varphi }_{\mathrm{g}}-\mathrm{Q},$(15) for the nonlinear elimination variants, is negligible.

Both problems use a common form for the material properties but use different coefficients. Our material model has a simple proportional relationship between the specific internal energy and the material temperature, namely,

(42) $e=c_{v}T,$(42)

where $c_v$ is known as the specific heat. The opacity has the following form:

(43) ${\sigma }_{a,g}={\sigma }_{t,g}=\rho C{T}^A({\bar{\nu }}_g{)}^B,$(43)

where ${\bar{\nu }}_g=\sqrt{{\nu }_g{\nu }_{g+1}}$ is the group-averaged photon energy, evaluated at the geometric average of the group bounds defining group $g$.

V.A. Infinite Medium Equilibration Problem

The first problem is an infinite medium problem where the radiation field starts hotter than the material, and then the two fields come into equilibrium. This mocks up the temporal behavior within a zone that is being heated by neighboring zones in multidimensional problems. The initial radiation field has a Planckian distribution at $T_{rad}=0.3$ keV, and the material starts at $T_{mat}=0.001$ keV, where 1 keV = $11604518.12$ K. The density is ρ = 0.05 g/cm³ = 50 kg/m³. The material specific heat is $c_v=0.01$ jerk/g·keV, where 1 jerk = 10⁹ J. The speed of light is c = $299.792458$ cm/shake (sh) = $299792458$ m/s, and 1 s = 10⁸ sh. For the opacity model, we set $C=10$, $A=-1$, and $B=-3$. We use ten logarithmically spaced groups with bounds (in kilo-electron-volts) of ($5\times {10}^{-3}$, $5\times {10}^{-2.6}$, $5\times {10}^{-2.2}$, $5\times {10}^{-1.8}$, $5\times {10}^{-1.4}$, $5\times {10}^{-1}$, $5\times {10}^{-0.6}$, $5\times {10}^{-0.2}$, $5\times {10}^{0.2}$, $5\times {10}^{0.6}$, and $5\times {10}^1$). The time step size $\mathrm{\Delta }t$ was held fixed for a given run, and we ran multiple simulations reducing the time step to check for temporal convergence starting at $\mathrm{\Delta }t={10}^{-2}$ sh down to $\mathrm{\Delta }t={10}^{-7}$ sh by factors of 10. The problem was run out to $t=0.02$ sh.

Figure 1 shows the temperature histories for this problem for the LMFG and fully converged nonlinear methods; all variations that fully converge the outer nonlinear loop get the same results. For large time steps, the LMFG method is not stable and shows violations of the maximum principle. This has been explored extensively in the context of Implicit Monte Carlo; see, for example, Ref. 12. The largest time step with a stable LMFG solution is $\mathrm{\Delta }t={10}^{-6}$ sh, and it does not appear converged until at least $\mathrm{\Delta }t={10}^{-7}$ sh. On the other hand, fully converging the outer nonlinear iterations is always stable, if inaccurate. The solution appears converged at $\mathrm{\Delta }t={10}^{-5}$ sh, which is larger than where LMFG is stable.

Fig. 1. Temperature with two of the iteration methods. All of the fully converged nonlinear methods get the same solutions, so only one is shown. The stability limits of LMFG are clearly violated for large time steps. Solid lines are the radiation temperature, and dashed lines are the material temperature

For the solutions at $\mathrm{\Delta }t={10}^{-5}$ sh and $\mathrm{\Delta }t={10}^{-6}$ sh, we plot the total number of multigroup linear iterations in Fig. 2. Generally, the LMFG method is the cheapest by this metric, but that ignores the solution quality shown in Fig. 1. Fully converging both the outer Newton and inner linear system (NLMFG) is clearly the most expensive, which was also observed by Brown and Woodward¹ and Lowrie.² Doing an inexact Newton step with only one inner linear multigroup solve also performs better. The nonlinear elimination variants NLMFG-NE and Inexact-Newton Multifrequency Grey with Nonlinear Elimination (INMFG-NE) performed better than the unaccelerated versions. Once the initial transient has passed, INMFG-NE has the same cost as LMFG. This is the true value of this method. It rapidly converges to the fully converged nonlinear solution when necessary, and it automatically transitions to an algorithm with the same cost as LMFG when the time steps are small enough for that method to work reasonably well.

Fig. 2. Linear iterations per time step, as measured by doing one solve for each of the diffusion equations in Eqs. (29)(29) $\begin{array}{rl}& \frac{1}{\mathrm{\Delta }t}{\varphi }_g^{\star }-\mathrm{\nabla }\cdot \frac{c}{3{\sigma }_{t,g}}\mathrm{\nabla }{\varphi }_g^{\star }+c{\sigma }_{a,g}{\varphi }_g^{\star }=-J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}c{\sigma }_{a,g^{\prime }}{\varphi }_{g^{\prime }}^l\\ & +\frac{1}{\mathrm{\Delta }t}{\varphi }_g^{n-1}+c{\sigma }_{a,g}{B}_g+S_g+J_{{\varphi }_g,e}{J}_{e,e}^{-1}\\ & \left(\frac{\rho }{\mathrm{\Delta }t}\left(e^l-e^{n-1}\right)+c\sum _{g^{\prime }}{\sigma }_{a,g^{\prime }}{B}_{g^{\prime }}(T^l)-Q\right).\end{array}$(29) , (35)(35) $J_{{\varphi }_g,{\varphi }_g}{\hat{\varphi }}_g=J_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left({\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right)$(35) , and (33)(33) $\left[\sum _g{J}_{{\varphi }_g,{\varphi }_g}{\varsigma }_g-\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}{\varsigma }_{g^{\prime }}\right]{\varphi }^{†}=\left[\sum _g{J}_{{\varphi }_g,e}{J}_{e,e}^{-1}\sum _{g^{\prime }}{J}_{e,{\varphi }_{g^{\prime }}}\left[{\varphi }_{g^{\prime }}^{\star }-{\varphi }_{g^{\prime }}^l\right]\right].$(33)

Table II shows the cumulative number of linear iterations for each method needed to get to the goal time of $t=0.02$ sh for the different time step sizes. By far the most effective way to reduce iterations and run time is to increase the time step size, as long as the solution is accurate. Also, note that INMFG-NE has the best performance for all the fully converged nonlinear solution methods, for the given discrete time and space system. In fact, even when LMFG is stable and producing reasonable results, INFMG-NE costs at most 40% more. Note that for $\mathrm{\Delta }t>{10}^{-6}$ sh, LMFG gives unphysical solutions, so these low iteration counts are not meaningful.

TABLE II Cumulative Linear Iterations to Reach $t=0.02$ sh in the Infinite Medium Problem

$\mathrm{\Delta }t$ (sh)	LMFG	NLMFG	NLMFG-NE	INMFG	INMFG-NE
${10}^{-7}$	276 242	488 653	476 482	412 338	400 183
${10}^{-6}$	40 000	71 009	60 841	50 170	40 602
${10}^{-5}$	4000	10 534	9107	6351	5469
${10}^{-4}$	400	1502	1404	845	775
${10}^{-3}$	40	233	226	128	124

V.B. Two-Dimensional Radiation Flow

We now test the methods on a multidimensional problem based on the one in Ref. 8 with one slight modification. Here, we will run in Cartesian X,Y coordinates while the original version used cylindrical R,Z coordinates.

Figure 3 shows the overall geometry together with sample solutions for the radiation and material temperatures. The problem is defined on a rectangle from $0\mathrm{c}\mathrm{m}<x<7\mathrm{c}\mathrm{m}$ and $0\mathrm{c}\mathrm{m}<y<1\mathrm{c}\mathrm{m}$. The mesh has uniform zone zones of $0.1$ cm except in the $y$ direction between $0.5\mathrm{c}\mathrm{m}<y<0.6\mathrm{c}\mathrm{m}$, where at $y=0.5$ cm, the zone size was $0.001$ cm. Each successive zone was $1.47392$ times bigger than the previous zone for a total of ten zones across this region.

Fig. 3. Geometry of the two-dimensional radiation flow problem. There is a low-density material surrounded by a high-density material, with a symmetry plane at $y=0$. The top half shows the material temperature at $t=1$ sh, and the bottom half shows the radiation temperature at the same time but reflected about the $x$-axis. See Fig. 5 for a more quantitative view of the temperatures

The material in the region $y<0.5$ cm has a density of $\rho =0.01$ g/cm³, and for $y>0.5$ cm, the density is $\rho =\mathrm{2.0.}$g/cm³. The specific heat is $c_v=0.05$ jerk/g·keV. The opacity model coefficients are $C=10$, $A=-0.5$, and $B=-3$. There are ten radiation groups with group bounds in kilo-electron-volts of (${10}^{-4}$, ${10}^{-3.5}$, ${10}^{-3}$, ${10}^{-2.5}$, ${10}^{-2}$, ${10}^{-1.5}$, ${10}^{-1}$, ${10}^{-0.5}$, ${10}^0$, ${10}^{0.5}$, and ${10}^1$).

Reflecting boundary conditions were applied on all boundaries, except on the left at $x=0$ cm and for $y<0.5$ cm, where an isotropic incoming flux injected energy into the problem. For the reflecting boundaries, the coefficients in Eq. (12(12) $c{\alpha }_g{\varphi }_g+{\beta }_g\frac{c}{3{\sigma }_{t,g}}\hat{\mathbf{n}}\cdot \mathrm{\nabla }{\varphi }_g={\gamma }_g,$(12) ) are ${\alpha }_g=0$, ${\beta }_g=1$, and ${\gamma }_g=0$. The incoming flux boundary source linearly increases from $T_{source}=0.05$ keV to $T_{source}=0.5$ keV over the range $0<t<2$ sh and remains constant after the initial ramp. For the coefficients in Eq. (12(12) $c{\alpha }_g{\varphi }_g+{\beta }_g\frac{c}{3{\sigma }_{t,g}}\hat{\mathbf{n}}\cdot \mathrm{\nabla }{\varphi }_g={\gamma }_g,$(12) ), this implies ${\alpha }_g=1/2$, ${\beta }_g=1$, and

(44) ${\gamma }_g(t)=2B_g\left(\mathrm{m}\mathrm{a}\mathrm{x}\left(0.5,(0.5-0.05)/2t\right)\right).$(44)

We tested the algorithms with three different time step sizes of $\mathrm{\Delta }t=0.001$ sh, $\mathrm{\Delta }t=0.01$ sh, and $\mathrm{\Delta }t=0.1$ sh. The goal time is $t=5$ sh, which is shorter than the goal time of $20$ sh in Ref. 8. At $\mathrm{\Delta }t=0.1$ sh, the LMFG method is unstable and showed overheating in the material much like in the infinite medium problem. This can be seen in Fig. 4. All of the methods that fully converge the nonlinear system get the same, more physically realistic, solution. To more quantitatively demonstrate this, the radiation and material temperatures along the line $y=0.2$ cm are shown at $t=1$ sh for several methods and time step sizes in Fig. 5. It is clear that the problem is not converged in time at the larger time step size, but the fully converged solutions are capturing the timescale of the Marshak wave traveling through the material reasonably well.

Fig. 4. Material temperature at $\mathrm{\Delta }t=0.1$ sh. The unstable LMFG method is on top. All of the other methods that fully converge the nonlinear system produce the solution on the bottom

Fig. 5. The material and radiation temperatures in the pipe problem at $t=1$ sh along the line $y=0.2$ cm. The solutions of LMFG and INMFG-NE are shown for $\mathrm{\Delta }t=0.1$ sh and of INMFG-NE for $\mathrm{\Delta }t=0.001$ sh

The per time step linear iteration counts are shown for the $\mathrm{\Delta }t=0.01$ sh in Fig. 6. For this problem, we see again that the INMFG and INMFG-NE methods outperform the NLMFG and NLMFG-NE methods. The nonlinear elimination only marginally improves the convergence in this problem. The cumulative total number of iterations to reach the end time of $t=5$ sh is shown in Table III. This demonstrates that for a variety of time step sizes, INMFG-NE has marginally more cost than LMFG at increased stability and accuracy. In fact, given the particular tolerances selected, it appears that the first iteration of the nonlinear methods was just above the tolerance level; slightly loosening the tolerances would have made the nonlinear methods converge in just one iteration too. Fully converging both the inner and outer iteration loops in NLMFG and NLMFG-NE is dramatically more expensive, which mirrors the previous observations in Refs. 1 and 2. Note that for $\mathrm{\Delta }t=0.1$ sh, LMFG gives unphysical solutions.

TABLE III Cumulative Linear Iterations to Reach $t=5$ sh in the Pipe Problem

$\mathrm{\Delta }t$(sh)	LMFG	NLMFG	NLMFG-NE	INMFG	INMFG-NE
0.001	49 464	88 043	85 354	56 431	56 476
0.010	20 082	38 790	38 464	23 191	23 187
0.100	5798	15 325	15 518	6389	6351

Fig. 6. The number of linear iterations per time step for each of the different methods at $\mathrm{\Delta }t=0.01$ sh. Note that the nonlinear elimination variants only marginally improve the unaccelerated methods for this problem, so the curves are nearly the same

VI. CONCLUSIONS AND FUTURE WORK

We recently discovered that we could apply the nonlinear elimination method in Ref. 13 to describe the algorithm in our discrete ordinates TRT code Teton. This has led to a mathematical framework for describing a family of accelerated nonlinear solution methods. By proper choices of tolerances at various levels of the new framework, we can recover the standard LMFG method in Refs. 8 and 11, the (previously unpublished) algorithm in Teton that we call Inexact-Newton Multifrequency-Grey with Nonlinear Elimination (INMFG-NE), as well as several other algorithms.

The INFMG-NE method’s efficacy is due to three key algorithmic components. First, a standard grey acceleration step significantly speeds up convergence of the absorption and reemission of photons within a Newton iteration. Second, the new nonlinear elimination step effectively accelerates the emission and reabsorption within a Newton step. Because the emission is fundamentally nonlinear through the Planck function, it must be handled with a nonlinear method. Third, combining the grey acceleration and nonlinear elimination steps with an inexact Newton method allows only a single multigroup linear iteration per Newton step to be robust. All three in combination makes fully converging the nonlinear system practical.

In two test problems, we find (unsurprisingly) that fully converging the nonlinear system (with a given space and time discretization error) allows for larger time steps and better accuracy than the linearized one-step method LMFG. We also find that the traditional NLMFG method, which fully converges both the multigroup linear system for each Newton step and the outer nonlinear iterations, is prohibitively expensive. By combining the three acceleration techniques that make up INMFG-NE, we find that we can compute a fully converged nonlinear solution with very nearly the same cost as converging only the inner multigroup linear system in the LMFG method, particularly given that the increased robustness and accuracy of the fully converging nonlinear system allows for larger time steps. The INMFG-NE method seamlessly transitions from doing more iterations where needed for large nonlinear excursions to an algorithm with the cost of the LMFG method when the problem could be robustly linearized.

In the future, we plan to formally apply this method to the full Boltzmann transport equation, including inelastic Compton scattering. We are also excited to investigate the nonlinear preconditioning literature^13–15 for additional methods that may be more effective for solving the thermal radiation transport equations.

Acknowledgments

Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the U.S. Department of Energy, National Nuclear Security Administration under contract DE-AC52-07NA27344. This document was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor Lawrence Livermore National Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or Lawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall not be used for advertising or product endorsement purposes. Document number: LLNL-JRNL-797841.

Notes

^a Later, these will be the material and radiation balance equations.

^b In methods that only do one Newton step, we often think of this as the absorption and reemission for the entire time step.