# Benchmark solutions for radiative transfer with a moving mesh and exact
uncollided source treatments

William Bennett
Ryan G. McClarren

(August 28, 2024)

###### Abstract

The set of benchmark solutions used in the thermal radiative transfer
community suffer some coverage gaps, in particular nonlinear, nonequilibrium
problems. Also, there are no nonequilibrium, optically thick benchmarks. These
shortcomings motivated the origination of a numerical method unfettered by
preconditions of linearity and easily able to converge on smooth optically
thick problems, a moving mesh Discontinuous Galerkin (DG) framework that
utilizes an uncollided source treatment. Having already proven this method on
time dependent scattering transport problems, we present here solutions to
non-equilibrium thermal radiative transfer problems for familiar linearized
systems and more physical nonlinear systems in both optically thin and thick
regimes, including both the full transport and the $S_{2}$/$P_{1}$ solution.
Geometric convergence is observed for smooth sources at all times and some
nonsmooth sources at late times when there is local equilibrium and accurate
solutions are achieved for step sources when the solution is not smooth.

##  1 Introduction

The Stefan-Boltzmann law [1], which describes the relationship between
radiation emitted from a material and its temperature as proportional to
temperature to the fourth power, is to blame for the obdurate nonlinearity in
high-energy density radiative heat transfer models. For this reason, the
extant analytic benchmarks with space and time dependence in this field are
predicated on assumptions of linearity or equilibrium. There are solutions
that assume the radiation energy and internal energy of the material instantly
equilibrate, inducing a Marshak Wave [2, 3, 4], and solutions for non-
equilibrium problems that linearize the system in $T^{4}$ by invoking a form
for the the material heat capacity that is proportional to temperature cubed,
an innovation of Pomraning [5]. This technique of defining a heat capacity to
linearize the system has been used to produce an abundance of solutions,
including transport treatments for the $P_{1}$ equations [6], full transport
solutions with one speed [7, 8], and non-grey problems [9]. In the diffusion
limit, benchmarks solutions have been provided for one temperature [10], three
temperature [11], and a non-grey treatment [12].

While these solutions are invaluable to code developers for verification, it
is necessary to point out that there are certain drawbacks to using linear
problems to verify codes whose purpose it is to solve nonlinear systems. While
ideally the numerical code in question would solve the fully nonlinear
equations and implement a special equation of state when running these
verification problems, there is nothing to prohibit the curators of these
codes from simply solving a linearized system when the benchmark is being run.
The result is the same in both cases. Also, solutions to linear systems can be
scaled to match benchmarks, unlike the unforgiving solutions to nonlinear
systems, and the solution to the linearized equations equilibrates more
quickly than a nonlinear problem as the temperature increases as a result of
the special equation of state. Solving a linear problem does not completely
verify the functionality of a radiative transfer code.

Although nonlinearity is an impediment for analytic methods, it is not
necessarily a source of difficulty for spectral methods. This was the impetus
for our development of a moving mesh, uncollided source treatment
Discontinuous Galerkin (DG) method for solving transport problems [13]. The
time dependent cell edges, which we call a moving mesh, and the uncollided
source treatment were added on to the DG implementation because the transport
equation with its finite wavespeeds admits discontinuities that inhibit DG
methods from attaining their higher order convergence potential. The moving
mesh and uncollided source can present a smoother problem for the method to
solve: the moving mesh by matching edges to moving wavefront discontinuities
and the uncollided source by analytically resolving the most structured part
of the solution.

As documented in [13], we have already conducted extensive tests with this
method on time-dependent transport problems, which allowed for a detailed
analysis of the efficacy of the moving mesh and uncollided mesh for different
source types. For example, for finite width, nonsmooth sources that induce a
nonsmooth solution that is smoothed over time, the method proved the most
beneficial when compared to a standard DG implementation, but displayed only
algebraic error convergence, not the optimal geometric convergence that DG
methods are capable of. For smooth, Gaussian sources however, we we able to
achieve spectral convergence, but the importance of the moving mesh was
diminished.

With an understanding of the effectiveness of this method on linear systems,
we apply it to nonlinear radiative transfer problems and obtain results with
accuracy comparable to an analytic solution, which is the stated intent of
this work. Before attempting fully nonlinear problems however, we first apply
our method the existing linear radiative transfer problems. This will allow us
to test our method on problems with known solutions and uncover deficiencies
in a more forgiving arena. For nonlinear problems, we can still gauge the
precision of our solution by inspecting magnitude of the expansion
coefficients and the accuracy by checking against existing numerical $S_{n}$
solvers.

We selected the Su and Olson transport benchmark [8] as an ideal verification
solution for our method. Unfortunately, the results are not given to enough
digits to fully demonstrate the effectiveness of our scheme and recalculation
of these results is non-trivial. Therefore, we rely on integration of a
$P_{1}$ version of this benchmark [6] to create solutions which are not
necessarily as physically accurate as the full transport solution, but can be
calculated out to more digits.

There are no existing radiative transfer transport verification solutions for
optically thick problems outside of the equilibrium diffusion limit. By
optically thick, we mean that the source width or the support of an initial
condition is orders of magnitude larger than a mean free path. Conversely,
optically thin problems have source widths comparable to a mean free path. For
a transport code to have sufficient coverage of verification problems,
converging to a diffusion benchmark of a thick problem while not resolving a
mean free path is a good test. If the code resolves a mean free path, however,
it will converge to the diffusion problem plus a transport correction. It is
for the purpose of verifying this transport correction that we include
transport solutions and the $S_{2}$/$P_{1}$ solutions for optically thick
problems.

The remaining sections of the paper are organized as follows. Section 2
contains an introduction to our model equations, nondimensionalization, and
derivation of the uncollided source. Our DG implementation is laid out briefly
in Section 3, but a more detailed derivation is left to [13]. Section 4 is
devoted to the calculation of $S_{2}$ benchmarks and the corresponding
uncollided solutions used in verifying our method. Following this is a
description of how the convergence of the error in the results section is
calculated (Section 5. then our results (Sections 6 and 7). The results
sections also contain specific details of the methods used in each problem and
discussion of the solution characteristics.

##  2 Equations

We study non-equilbrium time dependent radiative heat transfer in an infinite,
purely absorbing, stationary medium with an internal radiation source. The
radiation transport and material balance equations for this system are,

$\left(\frac{1}{v}\frac{\partial}{\partial\tau}+\mu\frac{\partial}{\partial
z}+\sigma_{a}\right)\psi(z,\tau,\mu)=\sigma_{a}\left(\frac{1}{2}avT(z,\tau)^{4}\right)+\frac{1}{2}S(z,\tau),$
(1)
$\frac{\partial}{\partial\tau}e(z,\tau)=\sigma_{a}\left(\phi(z,\tau)-avT(z,\tau)^{4}\right),$
(2)

where the general form of the equation of state is,

$e=\int_{0}^{T}\\!dT^{\prime}\,C_{\mathrm{v}}(T^{\prime}).$ (3)

The variables in these equations are, $\psi$, the angular flux or intensity,
$\phi=\int_{-1}^{1}\\!d\mu^{\prime}\psi(x,t,\mu^{\prime})$ the scalar flux,
$T$, the temperature, and $e$, the material energy density. $\psi$ and $\phi$
have units of energy per area per time ([GJ$\cdot$cm-2ns-1]) and $e$ has units
of energy density ([GJ$\cdot$cm-3]). $S$ is a source term with units of energy
density per time. $\mu\in[-1,1]$ is the cosine of the particle direction with
respect to the $z$ axis. $v$ is the particle velocity, which is the speed of
light in a vacuum for our application, $v=29.998$ cm ns-1. The radiation
constant is $a=4\sigma_{\mathrm{SB}}/v=0.0137225$ GJ cm-3 keV-4, where
$\sigma_{\mathrm{SB}}$ is the Stefan-Boltzman constant. The absorption cross
section, $\sigma_{a}$, is in units of inverse length.

We seek a non-dimensionalization for these equations that is compatible with
the non-dimensionalization given the in Su-Olson benchmark [9] and that may be
used in optically thick problems without enlarging the non-dimensionalized
length to accommodate the larger opacity,

$x=l\sigma_{a}z\qquad t=lv\sigma_{a}\tau.$ (4)

$l$ is a dimensionless scaling variable that is set to one for thin problems
and a small number to offset the greater $\sigma_{a}$ in optically thick
problems.

Each equation is transformed into the new variables and divided by
$av\sigma_{a}T_{H}^{4}$, where $T_{H}$ is the reference temperature, called
the hohlraum temperature in previous work,

$\left(l\frac{\partial}{\partial t}+\mu l\frac{\partial}{\partial
x}+1\right)\overline{\psi}(x,t,\mu)=c_{a}\left(\frac{1}{2}\overline{T}(x,t)^{4}\right)+\frac{1}{2}Q(x,t),$
(5) $l\frac{\partial}{\partial
t}\overline{e}(x,t)=c_{a}\left(\overline{\phi}(x,t)-\overline{T}(x,t)^{4}\right).$
(6)

Our non-dimensional dependent variables are now,

$\overline{\psi}(x,t)=\frac{\psi(x,t)}{avT_{H}^{4}}\qquad\overline{\phi}(x,t)=\frac{\phi(x,t)}{avT_{H}^{4}}\qquad\overline{T}(x,t)=\frac{T(x,t)}{T_{H}}\qquad\overline{e}(x,t)=\frac{e(x,t)}{aT_{H}^{4}},$
(7)

the non-dimensional source is,

$Q(x,t)=\frac{S(x,t)}{\sigma_{a}avT_{H}^{4}}$ (8)

and the absorption ratio is defined,

$c_{a}=1.$ (9)

In this work, we consider two functional forms for $C_{\mathrm{v}}$. To solve
the Su-Olson benchmark problem, we use the familiar form,

$C_{\mathrm{v}}=\alpha T^{3}.$ (10)

which renders Eq. (6) linear in $T^{4}$. With the conventional choice of
$\alpha=4a$, now

$\overline{e}_{\mathrm{SU}}=\overline{T}^{4},$ (11)

where the subscript “SU” indicates that this is the equation of state for the
linear Su-Olson problem.

While it is important for our investigation to solve these linear problems the
novel aspect of this paper is results for nonlinear problems. For these, we
choose a more physical, constant specific heat,
$C_{\mathrm{v}}=\mathrm{C}_{\mathrm{v0}},$ with units of of energy density per
temperature. This choice renders $e=\mathrm{C}_{\mathrm{v0}}T.$ To find the
relationship between the nondimensional variables with this equation of state,
we define
$\overline{\mathrm{C}_{\mathrm{v0}}}=\frac{\mathrm{C}_{\mathrm{v0}}}{aT_{H}^{3}}.$
Now we can write

$\overline{e}_{\mathrm{N}}=\overline{\mathrm{C}}_{\mathrm{v0}}\overline{T},$
(12)

where the subscript “N” indicates that this is our equation of state for the
nonlinear problems.

Other forms of $C_{\mathrm{v}}$ are certainly allowable. For instance, we
could implement a specific heat that is dependent on a non-constant opacity.
However, this will be left for future work.

###  2.1 Uncollided solutions

In time dependent transport trials, we found that the deployment of an an
uncollided source treatment, where using the solution to the equation,

$\left(l\frac{\partial}{\partial t}+\mu l\frac{\partial}{\partial
x}+1\right)\overline{\psi}_{u}(x,t,\mu)=\frac{1}{2}Q(x,t),$ (13)

as a source term to solve for the collided flux is a significant boon for
accuracy when the solution is not smooth. The Green’s solution to Eq. (13)
with $l=1$ was provided by [14]. This solution is integrated for different
source configurations in [15], including a square and a Gaussian source. Using
these solutions, we can say that $\overline{\psi}_{u}$ is known and can be
integrated analytically to find $\overline{\phi}_{u}$. For problems where
$l\neq 1$, a simple scaling is required. For optically thick problems when
$l\ll 1$, the uncollided solution is not as useful since it has decayed to
zero by the pertinent evaluation times.

To solve for the remaining collided portion of the flux, we have the system,

$\left(l\frac{\partial}{\partial t}+\mu l\frac{\partial}{\partial
x}+1\right)\overline{\psi}_{c}(x,t,\mu)=c_{a}\left(\frac{1}{2}\overline{T}(x,t)^{4}\right)$
(14) $l\frac{\partial}{\partial
t}\overline{e}(x,t)=c_{a}\left(\overline{\phi}_{c}(x,t)+\overline{\phi}_{u}(x,t)-\overline{T}(x,t)^{4}\right).$
(15)

In linear transport applications, it is possible to decompose the flux
infinitely, not just into uncolllided and collided flux, but uncollided, first
collided, second collided, etc. Even though the radiative transfer equations
are nonlinear, we are able to use this linear solution technique since the
uncollided flux has no interaction with the material. However, we cannot not
further decompose the flux as we could in a linear transport problem.

Answers obtained with an “uncollided source” treatment refer to solutions to
Eqs. (14) and (15). A “standard source” treatment refers to Eqs. (5) and (6).
The most useful source treatment used in a specific problem is determined by
the behavior of the uncollided flux during the solution time window. Tests run
in [13] showed that integrating the uncollided source could require more
computation time than a standard source. This is because the uncollided source
is a complex function of space and more difficult to integrate with quadrature
than the standard source. As a rule, problems at times where the uncollided
flux has not decayed enough to be a negligible portion of the flux are good
candidates for an uncollided source treatment. In these problems, [13] showed
increase in accuracy and rate of convergence. At times where the uncollided
solution has decayed, the uncollided source treatment is not as helpful.

While the problems investigated in this paper are in purely absorbing media,
the coupling between the material energy density and the radiation energy
density acts as a scatterer in that it can smooth discontinuities over time.
For this reason, we expect that the insights derived from solving purely
scattering transport problems with uncollided source treatments will extend to
these purely absorbing radiative transfer problems.

##  3 Moving Mesh DG spatial discretization

Similar to the procedure in [13], we define a DG spatial discretization with a
moving mesh to solve equations of the form (5) and (6). We leave some of the
details of the derivation to [13]. To solve for the integral over $\mu$ to
find the scalar flux, we discretize in angle via the method of discrete
ordinates, where the solid angle $\mu\in[-1,1]$ is discretized by choosing the
points with a Gauss-Lobatto quadrature rule [16] for our full transport
solution or, in the case of the $S_{2}$ solution, a Gauss-Legendre rule. With
the corresponding weights from our chosen quadrature, we can define the scalar
flux as a weighted sum,

$\overline{\phi}\approx\sum^{N}_{n^{\prime}=1}w_{n^{\prime}}\overline{\psi}^{n^{\prime}},$
(16)

where $w_{n}$ are the weights and $\psi^{n}$ is the scalar flux evaluated at a
given angle. This choice makes Eq. (5) and (6),

$\left(l\frac{\partial}{\partial t}+\mu_{n}l\frac{\partial}{\partial
x}+1\right)\overline{\psi}^{n}(x,t)=c_{a}\left(\frac{1}{2}\overline{T}(x,t)^{4}\right)+\frac{1}{2}Q(x,t)\>\>\>\mathrm{for}\>n=1\dots
N,$ (17) $\frac{\partial}{\partial
t}\overline{e}(x,t)=c_{a}\left(\sum^{N}_{n^{\prime}=1}w_{n^{\prime}}\psi^{n^{\prime}}-\overline{T}(x,t)^{4}\right).$
(18)

To discretize the spatial domain, we define $K$ non-overlapping cells with
time dependent edges $x_{L}(k,t)$ and $x_{R}(k,t)$. To allow simplifications
to the coming weak form of the equations, we define a mapping variable
$x^{\prime}(k,t)$ that maps $x$ to [-1,1] inside a cell,

$x^{\prime}(k,t)\equiv\frac{x_{L}(k,t)+x_{R}(k,t)-2x}{x_{L}(k,t)-x_{R}(k,t)},\qquad
k=1\dots K.$ (19)

Now we define an orthonormalized Legendre polynomial basis function in
$x^{\prime}$ for each cell,

$B_{i,k}(x^{\prime})=\frac{\sqrt{2i+1}}{\sqrt{x_{R}(k,t)-x_{L}(k,t)}}P_{i}(x^{\prime}).$
(20)

Therefore, the weak solution of the angular flux in a cell for a given angle
is

$\overline{\psi}^{n}(x,t)\approx\sum_{j=0}^{M}B_{j,k}(x^{\prime})\,u^{n}_{k,j}.$
(21)

where $u$ is an entry in our three dimensional solution matrix. Likewise, the
solution for the energy density in a given cell is,

$\overline{e}(x,t)\approx\sum_{j=0}^{M}B_{j,k}(x^{\prime})\,u^{N+1}_{k,j},$
(22)

The standard DG procedure for finding the weak form of the equations involves
multiplying each equation by a basis function, integrating over a cell,
invoking integration by parts to shift the spatial derivative onto the basis
function, and taking advantage of orthogonality to simplify the system. Our
moving mesh method is similar to this, but with the added step of invoking the
Reynolds Transport Theorem [17] since our integration domain is time
dependent. Leaving the general outline of this procedure to [13], we arrive
at,

$\frac{d}{dt}\boldsymbol{U_{n}}-\boldsymbol{\underline{\underline{G}}U_{n}}+\left(\boldsymbol{\underline{\underline{L}}U_{n}}\right)^{(\mathrm{surf})}-\mu_{n}\boldsymbol{LU_{n}}+\frac{1}{l}\boldsymbol{U_{n}}=\frac{c_{a}}{2l}\boldsymbol{H}+\frac{1}{2l}\boldsymbol{Q}\>\mathrm{for}\>n=1\dots
N,$ (23)
$\frac{d}{dt}\boldsymbol{U_{N+1}}+\boldsymbol{RU}^{\mathrm{surf}}-\boldsymbol{\underline{\underline{G}}U}_{N+1}=\frac{c_{a}}{l}\left(\sum_{n^{\prime}=1}^{N}w_{n}^{\prime}\boldsymbol{U}_{n^{\prime}}-\boldsymbol{H}\right),$
(24)

where the time dependent solution vector is

$\boldsymbol{U}_{n,k}=[u^{n}_{k,0},u^{n}_{k,1},...,u^{n}_{k,M}]^{T},$

where $M+1$ is the number of basis functions. We also define

$L_{i,j}=\int_{x_{L}}^{x_{R}}\\!\,dx\,B_{j,k}(x^{\prime})\,\frac{dB_{i,k}(x^{\prime})}{dx},$
(25)
$G_{i,j}=\int^{x_{R}}_{x_{L}}\\!dx\,B_{j,k}(x^{\prime})\,\frac{dB_{i}(x^{\prime})}{dt},$
(26) $Q_{i}=\int^{x_{R}}_{x_{L}}\\!dx\,B_{i,k}(x^{\prime})\,Q(x,t),$ (27)
$H_{i}=\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i}(x^{\prime})\,\overline{T}^{4}(x,t).$
(28)

The numerical flux terms, which calculate the direction of flow of the
solution with an upwinding scheme based on the relative velocity of a particle
with the cell edges

$\left(LU\right)^{\mathrm{surf}}_{i}=\left(\mu_{n}-\frac{dx_{R}}{dt}\right)B_{i,k}(x^{\prime}=1)\boldsymbol{\overline{\psi}^{n+}}-\left(\mu_{l}-\frac{dx_{L}}{dt}\right)B_{i,k}(x^{\prime}=-1)\boldsymbol{\overline{\psi}^{n-}}.$
(29)

$\left(RU\right)^{\mathrm{surf}}_{i}=\left(-\frac{dx_{R}}{dt}\right)B_{i,k}(x^{\prime}=1)\boldsymbol{\overline{e}^{+}}-\left(-\frac{dx_{L}}{dt}\right)B_{i,k}(x^{\prime}=-1)\boldsymbol{\overline{e}^{-}}.$
(30)

$\boldsymbol{\overline{\psi}^{l+}}$ and $\boldsymbol{\overline{\psi}^{l-}}$
are found by evaluating Eq. (21) and $\boldsymbol{\overline{e}^{+}}$ and
$\boldsymbol{\overline{e}^{-}}$ are found by evaluating Eq. (22).

If we choose to employ an uncollided source treatment, Eqs. (23) and (24)
change slightly in that the source term $\boldsymbol{Q}$ disappears from the
RHS of Eq. (23) and $\frac{c_{a}}{l}\boldsymbol{\overline{\phi}}$ is added to
the RHS of Eq. (24) where,

$\overline{\boldsymbol{\phi_{u}}}=\int_{x_{L}}^{x_{R}}\\!\,dx\,B_{j,k}(x^{\prime})\,\overline{\phi_{u}}(x,t).$
(31)

In this case, the numerical solution is for the collided flux, so it is
necessary to add the uncollided flux at the final step to obtain the full
solution.

The general solution procedure is as follows. First, parameters such as the
number of basis functions, the number of spaces, and the $S_{n}$ order are
set. A source is specified and depending on the source treatment, the
uncollided solution or the standard source is integrated at each timestep with
a standard Gaussian integrator with points equal to $2M+1$. The edges of the
mesh are governed by a function designed to optimize the solution for the
specific source. This function also returns velocities of the mesh edges in
order to calculate the numerical flux. The temperature balance terms are found
by the equation of state and integrated in the same way as the source term.
The solver returns the coefficient arrays and the scalar flux and material
energy are reconstructed via the expansions defined in Eqs. (21) and (22) and,
depending on the source treatment, the uncollided flux is added onto the
scalar flux.

To obtain solutions from our equations (23) and (24), we calculate the
quadrature weights with the python package quadpy [18] and integrate the ODES
in with a built in integrator from scipy [19]. Our python implementation can
be found on Github111www.github.com/wbennett39/moving_mesh_radiative_transfer.

##  4 Benchmarks and uncollided solutions for the $S_{2}$ radiative transfer
equations

In order to show more decimal places of accuracy in our linear problems than
are given in [8], it was expedient to calculate our own analytic benchmarks.
Since we already include results to each problem calculated with $S_{2}$, we
choose to verify our solver by using a $S_{2}/P_{1}$ benchmark given by [6].
This benchmark gives the analytic expression for the scalar flux and energy
density solutions to Eq. (17) and (18) with $N=2$, angles and Gauss-Legendre
weighting, ($[\mu^{1},\mu^{2}]=[\frac{-1}{\sqrt{3}},\frac{1}{\sqrt{3}}]$, and
$[w_{1},w_{2}]=[1,1]$) when the source is a delta function in space and time.

The Green’s function given by [6] for
$\overline{\phi}=\overline{\psi}_{1}+\overline{\psi}_{2}$ for a delta function
source at position $s$ is,

$G(x,s,t)=\frac{v}{2\sqrt{3}}e^{-t}\left(\frac{tI_{1}\left(\sqrt{t^{2}-3(x-s)^{2}}\right)}{t^{2}-3(x-s)^{2}}\Theta\left(t-\sqrt{3}|x-s|\right)+I_{0}\left(\sqrt{t^{2}-3(x-s)^{2}}\right)\delta\left(t-\sqrt{3}|x-s|\right)\right),$
(32)

where $\Theta$ is a step function, $\delta$ is a Dirac delta function and
$I_{0}$ and $I_{1}$ are modified Bessel functions of the first kind. The
Green’s function for the material energy density is,

$G_{U}(x,s,t)=\frac{\sqrt{3}}{2}e^{-t}\left(I_{0}\left(\sqrt{t^{2}-3(x-s)^{2}}\right)\Theta\left(t-\sqrt{3}|x-s|\right)\right).$
(33)

We choose to find solutions for a square source and a Gaussian source, to test
our method on both smooth and nonsmooth problems. Therefore, the solution to
the integral,

$\overline{\phi}_{ss,gs}=\int^{\infty}_{-\infty}\\!ds\int^{\infty}_{0}\\!dt^{\prime}\,S_{ss,gs}(s,t^{\prime})\,G(x,s,t-t^{\prime})$
(34)

gives the scalar flux. The energy density is likewise obtained by

$\overline{e}_{ss,gs}=\int^{\infty}_{-\infty}\\!ds\int^{\infty}_{0}\\!dt^{\prime}\,S_{ss,gs}(s,t^{\prime})\,G_{u}(x,s,t-t^{\prime}).$
(35)

Where the subscript in the solution and the source is either “ss” for square
source or “gs” for Gaussian. The source term is,

$S^{ss}(x,t)=\Theta(x_{0}-x)\Theta(t_{0}-t),$ (36)

for the square source, or

$S^{gs}(x,t)=\exp\left(\frac{x^{2}}{x_{0}^{2}}\right)\Theta(t_{0}-t),$ (37)

for the Gaussian source.

Since finding a benchmark for an optically thick problem requires evaluating
these integrals at extremely late times where the integrand is not well
behaved, we only calculate $S_{2}$ benchmarks for our thin problems.

###  4.1 $S_{2}$ uncollided solutions

The uncollided solutions that we have utilized so far for the uncollided
source treatment have been full transport solutions from [15]. We cannot use
these to solve the $S_{2}$ transport equation, since the two uncollided fluxes
are not equal. The full transport solutions are based on the assumption that
the $S_{n}$ order of the ODE’s sufficiently resolves the angular error and
that the collided flux calculated with quadrature is a good approximation of
the analytic integral over $\mu$, i.e.

$\sum^{N}_{n^{\prime}=1}w_{n^{\prime}}\overline{\psi}^{n^{\prime}}\approx\int^{1}_{-1}\\!d\mu^{\prime}\,\psi(x,t,\mu^{\prime}),$
(38)

so that it is acceptable to employ the uncollided scalar flux found by
integrating analytically the solution for the uncollided angular flux. In the
$S_{2}$ equations, the assumption of Eq. (38) does not hold and the uncollided
scalar flux must be found by numerical quadrature of the angular flux.

Therefore, the process for finding the uncollided scalar flux to use as a
source in the $S_{2}$ solutions to our radiative transfer problems is to find
the Green’s solution for the angular flux, integrate that solution with
quadrature, and then integrate again over the given source. The uncollided
solution to Eq. 13 with a delta function source $(\delta(x)\delta(t))$ is
[14],

$\overline{\psi}_{u}(x,t)=\frac{e^{-t}}{2t}\delta\left(\mu-\frac{x}{t}\right).$
(39)

To find the $S_{2}$ uncollided scalar flux, the integral is done by Gauss-
Legendre quadrature with $N=2$ to give the uncollided scalar flux,

$\overline{\phi}_{u}^{pl}(x,t)=\frac{e^{-t}}{2t}\left(\delta\left(-\frac{1}{\sqrt{3}}-\frac{x}{t}\right)+\delta\left(\frac{1}{\sqrt{3}}-\frac{x}{t}\right)\right).$
(40)

To finally find the uncollided scalar flux that corresponds to the benchmark
solutions calculated with Eqs. (34), we integrate

$\overline{\phi}^{ss,gs}_{u}(x,t)=\int^{\infty}_{-\infty}\\!ds\int^{\infty}_{0}\\!dt^{\prime}\,\overline{\phi}_{u}^{pl}(x-s,t-t^{\prime})\,S^{ss,gs}(s,t^{\prime}),$
(41)

where $S^{ss,gs}$ is given by Eq. (36) or Eq. (37). Solutions to Eq. (41) are
given in Appendix A.

##  5 Error estimation methods

In the problems presented, two methods are used to estimate the solution
accuracy. For problems with a benchmark solution, we use the root mean square
error (RMSE) as our error metric. This is calculated by,

$\mathrm{RMSE}=\sqrt{\sum^{N}_{i}\frac{|y_{i}-\hat{y_{i}}|^{2}}{N}},$ (42)

where $y_{i}$ is either the calculated scalar flux or the calculated material
energy density at a given node, $\hat{y_{i}}$ is the corresponding benchmark
solution, and $N$ is the total number of nodes in the computational solution.

For problems that demonstrate geometric spectral convergence, as
$M\rightarrow\infty$, the error can be modeled as

$\mathrm{ERROR}=C\,\exp(-c_{1}M),$ (43)

where $M$ is the highest polynomial order of the basis and $C$ and $c_{1}$ are
constants that could depend on the number of cells used in the problem. This
curve is a straight line on a logarithmic-linear scale.

For all of the problems in the following section, we plot the average of the
absolute value of the coefficients in the solution expansion to characterize
the solution convergence. We define the average value of the $j^{\mathrm{th}}$
coefficient in the solution expansion,

$|c_{j}|=\frac{\sum_{k=1}^{K}|a_{k}|}{K},$ (44)

where $j$ corresponds to the order of the Legendre polynomial in the basis and
$K$ is the number of cells.

When characterizing the error of $\overline{\phi}$, since we are interested in
the residual error of scalar flux, $a_{k}$ is the weighted average using the
weights from Eq. (16),

$a_{k}=\frac{\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}u_{(l^{\prime},k,j)}}{\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}}.$
(45)

For the material energy density, $a_{k}$ is,

$a_{k}=u_{(N+1,k,j)}.$ (46)

##  6 Optically thin results

The results in this section are for problems where the source width is equal
to a mean free path, meeting our definition of an optically thin problem.
These problems are characterized by solutions where the uncollided solution is
a significant portion of the flux and travelling wavefronts. Therefore, the
problems in this section all use an uncollided source and the square sources
whhich have travelling discontinuities employ a moving mesh.

###  6.1 Su-Olson problem with a square source

We first replicate the the Su-Olson problem using the same square source
originally presented in [8] with $\sigma_{a}=1$ cm-1, the source width,
$x_{0}=0.5$ and the source duration $t_{0}=10$. The uncollided solution for
this source has already been presented in [15]. For the $S_{2}$ treatment of
this problem, the uncollided source is given by Eq. (58). The temperature is
calculated by Eq. (11).

Some modifications were made to the original mesh function invented to solve
the square source transport problem in [13]. In that mesh, the mesh edges
inside the source never moved while the edges outside travelled outwards with
the wavespeed. This was done to resolve the static discontinuities at the
source edge and the travelling discontinuities at the wavefront. In the
original Su-Olson results, the source turns off at $t_{0}=10$ and solutions
are required long afterwards ($t=31.6228,100$). With our previous square
source mesh, the edges would remain clustered around the source region long
after the source has ceased to introduce nonsmoothness. This is not the
optimal distribution of computational zones.

Therefore, the mesh function used in this problem is as follows. If the mesh
edges are defined as the vector,

$\boldsymbol{X}(t)=\left[x^{0}(t),x^{1}(t),...,x^{K}(t)\right],$ (47)

and initialized to be,

$\displaystyle\mathrm{if}\>\frac{K}{4}\leq
k\leq\frac{3K}{4},\>\>x^{k}_{o}=\frac{1}{2}y_{j},$ (48)
$\displaystyle\mathrm{if}\>k<\frac{K}{4},\>\>x^{k}_{o}=\frac{s_{k}(\delta
x)+2x_{0}+\delta x}{2}$ (49)
$\displaystyle\mathrm{if}\>\>k>\frac{3K}{4},\>\>x^{k}_{o}=\frac{s_{l}(\delta
x)-2x_{0}-\delta x}{2}$ (50)

where $y_{j}$ are the Gauss-Lobatto evaluation points with $N$, the number of
points, equal to $\frac{K}{2}+1$ numbered from $0$ and $s_{m}$ are the Gauss-
Lobatto evaluation points for $N=\frac{K}{4}+1$. The indices $j$ and $m$ are
equal to $k-\frac{K}{4}$ and $k-\frac{K}{2}+1$ respectively. $\delta x$ is a
small initial width, and $K$ is always an even number. This initialization
assigns one third of the edges to the source and the other two thirds to cover
the rest of the solution domain. Each subdomain is spanned by edges with
Gauss-Lobatto spacing, which has the effect of concentrating cells near the
source edges and the outgoing wavefronts, where discontinuities are most
likely.

As time progresses and the outside edges move outwards with the solution,
their position is defined as,

$\mathrm{if}\>\>t\leq
t_{0},\>\>x^{k}(t)=x^{k}_{o}+\frac{x^{k}}{x^{k}_{o}}\times vt,$ (51)

where $v$ is the wavespeed, one for the transport problems and
$\frac{1}{\sqrt{3}}$ for the $S_{2}$ problems. The edge velocity is defined,

$\mathrm{if}\>\>t\leq
t_{0},\>\>\frac{dx^{k}}{dt}=\frac{x^{k}}{x^{k}_{o}}\times vt,$ (52)

Defining the edge positions and velocities this way preserves the relative
spacing of the initialized edges, meaning that the edges are clustered at the
source edges and the leading wavefronts.

At later times when the source is off, the solution to a square source in an
optically thin problem will behave more like the solution for a Gaussian
source since the solution will become smoother without the source emitting
uncollided particles. Information flow is no longer dominated by the
wavespeed. The solution will be practically zero some distance from the origin
that is much less than $vt$. For instance, in the Su-Olson problem at $t=100$
the solution is practically zero past $x=\pm 30$.

For these reasons, when the source turns off, a constant acceleration will
divert the trajectory of each edge so that at the final time, they are evenly
spaced over a specified width. This width is an estimate to how far the
solution will have traveled by the evaluation time. We chose a constant
acceleration instead of a instantaneous velocity change because the latter
induced numerical errors which resulted in failure to converge to the
benchmark solutions. The acceleration for each edge is found by,

$c^{k}=2\frac{\left(\frac{dx^{k}}{dt}\bigg{|}_{t_{0}}(t_{0}-t_{\mathrm{final}})-x^{k}\bigg{|}_{t_{0}}+x^{k}\bigg{|}_{t_{\mathrm{final}}}\right)}{\left(t_{0}-t_{\mathrm{final}}\right)^{2}}$
(53)

where we find $x^{k}\bigg{|}_{t_{\mathrm{final}}}$ by specifying that the
final positions vector, $\boldsymbol{X}(t_{\mathrm{final}})$, evenly span
$[-x_{f}/2,x_{f}/2]$ where $x_{f}$ is our estimate for the width of the
solution domain at the evaluation time.

With the acceleration defined, we calculate positions of the edges after the
source has turned off with,

$\mathrm{if}\>\>t>t_{0},\>\>x^{k}(t)=\frac{1}{2}c^{k}\left(t-t_{0}\right)^{2}+\frac{dx^{k}}{dt}\bigg{|}_{t_{0}}\left(t-t_{0}\right)+x^{k}\bigg{|}_{t_{0}},$
(54)

and the velocities,

$\mathrm{if}\>\>t>t_{0},\frac{dx^{k}}{dt}=c^{k}\left(t-t_{0}\right)+\frac{dx^{k}}{dt}\bigg{|}_{t_{0}}.$
(55)

We refer to this method for governing the mesh edges and velocities as the
“thin square source mesh”. An example $x$ vs $t$ diagram of the mesh edges is
given in Figure 1 for clarification.

Figure 1: Edge position for each edge in the thin square source mesh run out
to $t_{\mathrm{final}}=31.6228$ with $8$ spaces, a wavespeed,
$v=\frac{1}{\sqrt{3}}$, and a final domain width, $x_{f}=30$.

After completing the necessary steps of defining a source, choosing a
functional form for the temperature, and defining a mesh, we may present our
results for this problem. Tables 1 and 2 give our solutions for the same
points and evaluation times as in Table 1 and 2 of [8]. The convergence
results for the coefficient expansions of these results are plotted in Figures
2 and 3 . The solutions at a few selected times are plotted in Figure 4 For
each case, a moving mesh and the uncollided source was used except in the
$S_{2}$ solution at times greater than $t_{0}$ since the $S_{2}$ uncollided
solution becomes sharp and difficult to resolve via quadrature. In these
cases, the standard square source was integrated.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 2: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
Su-Olson square source problem where $x_{0}=0.5$, $t_{0}=10$. The quadrature
order for all results is $S_{256}$. All results were calculated with a moving
mesh and uncollided source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 3: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44))) for the optically thin ($\sigma_{a}=1$
cm-1) $S_{2}$ Su-Olson square source problem where $x_{0}=0.5$, $t_{0}=10$.
All results were calculated with a moving mesh and uncollided source treatment
except for the $t=31.6228$ and $t=100$ cases where a standard source treatment
was used.

(a) $t=0.1$

(b) $t=0.31623$

(c) $t=1$

(d) $t=3.16228$

(e) $t=10$

(f) $t=100$

Figure 4: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$)
solutions for the optically thin Su-Olson square source problem with
$x_{0}=0.5$, $t_{0}=10$. Solid lines are scalar flux, $\overline{\phi}$, and
dashed are material energy density, $\overline{e}$.

The convergence results (Figures 2 and 3) show that for this problem, $S_{2}$
is considerably smoother at early times. Twice as many spatial divisions were
required at early times in the full transport solution $256$ to achieve
similar levels of convergence to the $S_{2}$ solutions. Both cases exhibited
similar behavior over time. At early times the significantly nonsmooth
uncollided flux induced discontinuities in the material energy density and
required far more spatial divisions to resolve, stiffening the problem and
limiting the number of basis functions that could reasonably be used in the
solution. After the source turned off, the solution smoothed and equilibrated
locally (complete equilibrium is impossible with an infinite material). The
solution became smoother and could easily be resolved with fewer spaces and
more basis functions. We also note that at later times after the source has
turned off, the full transport and $S_{2}$ solutions become more similar.

Table 1: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$ for the thin square source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. RMSE values included in the last row of the $S_{2}$ table are calculated from the $S_{2}$ benchmark from [6]. Convergence results for these answers are plotted in Figures 2 and 3. Bolded digits are those that agree with the published Su-Olson solution in [8]. x/t | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 |  0.095317 |  0.275294 |  0.643151 |  1.20069 |  2.235815 |  0.690187 |  0.357195
0.1 |  0.095317 |  0.275294 |  0.635943 |  1.188724 |  2.219553 |  0.689743 |  0.357137
0.17783 |  0.095317 |  0.275294 |  0.619626 |  1.162044 |  2.183558 |  0.688773 |  0.357011
0.31623 |  0.095317 |  0.262715 |  0.561896 |  1.071861 |  2.064534 |  0.685719 |  0.356612
0.45 |  0.08824 |  0.203128 |  0.447114 |  0.909526 |  1.860758 |  0.681168 |  0.356016
0.5 |  0.047658 |  0.137647 |  0.358083 |  0.799027 |  1.731816 |  0.679072 | 0.35574
0.56234 |  0.003762 |  0.062776 |  0.253722 |  0.666804 |  1.574955 |  0.67616 |  0.355355
0.75 | - |  0.002793 |  0.114315 |  0.446752 |  1.273984 |  0.665459 |  0.353929
1.0 | - | - |  0.036471 |  0.275396 |  0.987815 |  0.646922 |  0.351409
1.33352 | - | - |  0.002894 |  0.145309 |  0.708221 |  0.615381 |  0.346972
1.77828 | - | - | - |  0.059674 |  0.450163 |  0.563509 |  0.339223
3.16228 | - | - | - |  0.001155 |  0.096453 |  0.369659 |  0.303466
5.62341 | - | - | - | - |  0.003632 |  0.108305 |  0.213818
10.0 | - | - | - | - | - |  0.003914 |  0.072059
17.78279 | - | - | - | - | - | - |  0.002721
$\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 |  0.095317 |  0.275294 |  0.661668 |  1.0451 | 1.918396 |  0.659852 |  0.352728
0.1 |  0.095317 |  0.275294 |  0.645629 |  1.034885 | 1.906562 |  0.659502 |  0.352673
0.17783 |  0.095317 |  0.275294 |  0.604866 |  1.012559 | 1.88071 |  0.658738 |  0.352555
0.31623 |  0.095317 |  0.275294 |  0.515344 |  0.941841 |  1.798924 |  0.656328 |  0.352181
0.45 |  0.089213 |  0.179984 |  0.405572 |  0.835497 |  1.67622 |  0.652732 |  0.351621
0.5 |  0.047658 |  0.137647 |  0.358083 |  0.786079 |  1.619313 |  0.651071 |  0.351362
0.56234 |  0.000000 |  0.085419 |  0.299433 |  0.722866 |  1.545756 |  0.648764 |  0.351001
0.75 | - |  0.0000000 |  0.155119 |  0.552606 |  1.338003 |  0.640252 |  0.349662
1.0 | - | - |  0.027203 |  0.369788 | 1.092239 |  0.625397 |  0.347295
1.33352 | - | - |  0.000000 |  0.19493 |  0.816475 |  0.599795 |  0.343124
1.77828 | - | - | - |  0.056328 |  0.532488 |  0.556761 |  0.335829
3.16228 | - | - | - | - |  0.09838 |  0.384521 |  0.30198
5.62341 | - | - | - | - |  0.00028 |  0.116245 |  0.215656
10.0 | - | - | - | - | - |  0.002009 |  0.073828
17.78279 | - | - | - | - | - | - |  0.002307
RMSE | 2.656e-07 | 1.747e-07 | 1.642e-06 | 3.589e-07 | 2.647e-07 | 9.157e-08 | 6.128e-09
Table 2: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$ for the thin square source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. RMSE values included in the last row of the $S_{2}$ table are calculated from the $S_{2}$ benchmark from [6]. Convergence results for these answers are plotted in Figures 2 and 3. Bolded digits are those that agree with the published Su-Olson solution in [8]. x/t | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 |  0.004682 |  0.040935 |  0.271307 |  0.94687 |  2.111923 |  0.704991 |  0.359136
0.1 |  0.004682 |  0.040935 |  0.268692 |  0.937154 |  2.095970 |  0.704514 |  0.359078
0.17783 |  0.004682 |  0.040935 |  0.26264 |  0.915402 |  2.060646 |  0.703474 |  0.358949
0.31623 |  0.004682 |  0.04034 |  0.239814 |  0.840926 |  1.943709 |  0.700198 |  0.358544
0.45 |  0.004552 |  0.033142 |  0.188264 |  0.702883 |  1.742967 | 0.69532 |  0.357938
0.5 |  0.002342 |  0.020469 |  0.141918 |  0.604935 |  1.615402 |  0.693073 |  0.357657
0.56234 | 0.00005 | 0.00635 | 0.08838 |  0.48846 |  1.460394 |  0.689954 |  0.357266
0.75 | - |  0.000063 |  0.030141 |  0.306558 |  1.165912 |  0.678498 |  0.355816
1.0 | - | - | 0.00625 |  0.175192 |  0.889908 |  0.658685 |  0.353254
1.33352 | - | - |  0.000162 | 0.08352 |  0.625213 |  0.625066 |  0.348744
1.77828 | - | - | - |  0.029349 |  0.386884 |  0.570027 |  0.34087
3.16228 | - | - | - |  0.000183 |  0.076146 |  0.367269 |  0.304561
5.62341 | - | - | - | - |  0.002412 |  0.103114 |  0.213768
10.0 | - | - | - | - | - |  0.003426 |  0.071226
17.78279 | - | - | - | - | - | - |  0.002609
x/t | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 |  0.004682 |  0.040935 |  0.280241 |  0.847357 | 1.808991 |  0.672725 |  0.354597
0.1 |  0.004682 |  0.040935 |  0.273782 |  0.837875 | 1.797308 |  0.672354 |  0.354541
0.17783 |  0.004682 |  0.040935 |  0.261727 |  0.817144 | 1.771785 |  0.671544 |  0.354421
0.31623 |  0.004682 |  0.040936 |  0.225609 |  0.751409 |  1.691033 |  0.668989 |  0.354041
0.45 |  0.004642 |  0.030469 |  0.169249 |  0.652361 |  1.569863 |  0.665177 |  0.353472
0.5 |  0.002341 |  0.020467 |  0.141916 |  0.606256 |  1.513661 |  0.663418 |  0.353209
0.56234 |  0.000000 |  0.008542 |  0.10839 |  0.547588 |  1.441078 |  0.660972 |  0.352842
0.75 | - |  0.000000 |  0.038396 |  0.393514 |  1.23688 |  0.651954 |  0.351482
1.0 | - | - |  0.001762 |  0.236801 |  0.997162 |  0.636229 |  0.349077
1.33352 | - | - |  0.000000 |  0.101119 |  0.731352 |  0.609161 |  0.344841
1.77828 | - | - | - |  0.019412 |  0.462775 |  0.563768 |  0.337432
3.16228 | - | - | - |  0.000000 |  0.074037 |  0.383609 |  0.303078
5.62341 | - | - | - | - |  0.000086 |  0.110507 |  0.215664
10.0 | - | - | - | - | - |  0.001626 |  0.072987
17.78279 | - | - | - | - | - | - |  0.002196
RMSE | 1.233e-08 | 2.372e-08 | 9.65e-08 | 1.366e-07 | 4.052e-08 | 9.937e-08 | 5.936e-09

Since we claim to present benchmark quality results, we are obliged to discuss
the accuracy of Tables 1 and 2. While [8] claims to converge their solution to
four digits, the observant reader will see that in some cases only 3 digits
match. Given that our solutions are converged beyond this point, we believe
that our reported digits are correct.

###  6.2 Constant $C_{\mathrm{v}}$ problem with a square source

This problem uses the same source as the problem of the last section but with
a different functional form of the heat capacity. Using Eq. (12) our system
becomes nonlinear. We choose
$C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$.
This value was chosen to see an appreciable change in temperature during the
selected time window. Now that we no longer have the convenient condition that
$\overline{e}=\overline{T}^{4}$, the local equilibrium condition is not
$\overline{\phi}=\overline{e}$ as in the Su-Olson problem but
$\overline{\phi}^{1/4}=\overline{T}$. For this reason, the solution plots for
this problem and all subsequent constant $C_{\mathrm{v}}$ problems do not show
scalar flux and material energy density but rather radiation temperature and
material temperature.

Though we can no longer rely on benchmark solutions for this problem,, we can
be confident that our solution is converged by plotting the magnitude of the
coefficients and check for systematic errors with a $S_{n}$ solver. Also,
since the mesh method employed here is the same method described in Section
6.1, we can be confident that the mesh is not introducing error.

Solutions to this problem are plotted in Figure 7. We note that the problem is
not everywhere at equilibrium by $t=100$ as the Su-Olson problem is. Also, we
note that the solution does not travel as far. In the Su-Olson problem, the
specific heat is very small when temperature is small and increases with the
cube of the temperature. This has the effect of attracting the solution to
equilibrium. This effect is not present in a constant $C_{\mathrm{v}}$ case
and there is less incentive for the solution to fall into local equilibrium.
It is also noteworthy that at very early times ($t<1$) the scalar flux has not
interacted with the material as much as in the Su-Olson problem and is mostly
made up of the uncollided flux. This is apparent in Figures 5(a) and 6(a).

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 5: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
constant $C_{\mathrm{v}}$ square source problem where $x_{0}=0.5$, $t_{0}=10$.
The quadrature order for all results is $S_{256}$. All results were calculated
with a moving mesh and uncollided source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 6: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
$S_{2}$ constant $C_{\mathrm{v}}$ square source problem where $x_{0}=0.5$,
$t_{0}=10$. All results were calculated with a moving mesh and uncollided
source treatment except for the $t=31.6228$ and $t=100$ cases where a standard
source treatment was used.

(a) $t=0.1$

(b) $t=0.31623$

(c) $t=1$

(d) $t=3.16228$

(e) $t=10$

(f) $t=100$

Figure 7: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$)
solutions for the optically thin constant $C_{\mathrm{v}}$ square source
problem with $x_{0}=0.5$, $t_{0}=10$. Solid lines are radiation temperature
$\overline{\phi}^{1/4}$, and dashed are temperature, $\overline{T}$.

Since the solution has not fully equilibrated at later times, the solutions
are less smooth compared to the Su-Olson problem. The repercussions of this
can be observed by comparing the convergence results at late times for this
problem in Figures 5 and 6 to the convergence results of the Su-Olson problem
in Figure 2. Also, the convergence results show that the material energy
density is generally more nonsmooth than the scalar flux. Nevertheless, we are
satisfied with the convergence of these results and present them in Tables 3
and 4.

The difference between the full transport solution and our $S_{2}$ result is
also of interest, as it provides insight into the physical characteristics of
the system. We note that the two solutions only begin to look similar at later
times as the solution equilibrates. This tells us that the solution becomes
less angularly dependent and better approximated by only two angles.

Table 3: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thin square source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.5$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 5 and 6 $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.095162 | 0.271108 | 0.563683 | 0.765084 | 1.96832 | 0.267247 | 0.085108
0.1 | 0.095162 | 0.271108 | 0.557609 | 0.756116 | 1.950367 | 0.266877 | 0.085054
0.17783 | 0.095162 | 0.271108 | 0.543861 | 0.736106 | 1.910675 | 0.266071 | 0.084937
0.31623 | 0.095162 | 0.258592 | 0.495115 | 0.668231 | 1.779896 | 0.263527 | 0.084565
0.45 | 0.08809 | 0.199962 | 0.396442 | 0.543721 | 1.558248 | 0.259729 | 0.084008
0.5 | 0.047581 | 0.135554 | 0.316071 | 0.453151 | 1.420865 | 0.257976 | 0.08375
0.56234 | 0.00376 | 0.061935 | 0.222261 | 0.349209 | 1.252213 | 0.255538 | 0.083392
0.75 | - | 0.002788 | 0.102348 | 0.21078 | 0.908755 | 0.246543 | 0.082061
1.0 | - | - | 0.034228 | 0.124305 | 0.562958 | 0.230831 | 0.079715
1.33352 | - | - | 0.002864 | 0.067319 | 0.27752 | 0.203718 | 0.075591
1.77828 | - | - | - | 0.031357 | 0.120054 | 0.158039 | 0.068419
3.16228 | - | - | - | 0.001057 | 0.013737 | 0.022075 | 0.036021
5.62341 | - | - | - | - | 0.000413 | 0.000814 | 0.001068
10.0 | - | - | - | - | - | 5e-06 | 5e-06
17.78279 | - | - | - | - | - | - | -
$\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.095162 | 0.271108 | 0.579404 | 0.649946 | 1.57957 | 0.243368 | 0.081882
0.1 | 0.095162 | 0.271108 | 0.56606 | 0.641494 | 1.567037 | 0.243074 | 0.081833
0.17783 | 0.095162 | 0.271108 | 0.529958 | 0.623092 | 1.539665 | 0.242432 | 0.081725
0.31623 | 0.095162 | 0.271108 | 0.45241 | 0.565417 | 1.453129 | 0.240406 | 0.081383
0.45 | 0.08906 | 0.177033 | 0.357556 | 0.480232 | 1.323479 | 0.237381 | 0.080871
0.5 | 0.047581 | 0.135554 | 0.31607 | 0.441134 | 1.263425 | 0.235984 | 0.080634
0.56234 | - | 0.084378 | 0.264889 | 0.392645 | 1.185561 | 0.234042 | 0.080305
0.75 | - | - | 0.140336 | 0.276699 | 0.962617 | 0.226873 | 0.079082
1.0 | - | - | 0.02637 | 0.175485 | 0.692975 | 0.214346 | 0.076927
1.33352 | - | - | - | 0.097064 | 0.392628 | 0.192725 | 0.073141
1.77828 | - | - | - | 0.033465 | 0.155709 | 0.156395 | 0.066564
3.16228 | - | - | - | - | 0.007509 | 0.03175 | 0.03718
5.62341 | - | - | - | - | 4.7e-05 | 0.000413 | 0.001075
10.0 | - | - | - | - | - | - | -
17.78279 | - | - | - | - | - | - | -
Table 4: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thin square source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.5$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 5 and 6 $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.004837 | 0.045121 | 0.354022 | 1.613529 | 2.57461 | 1.592549 | 1.190296
0.1 | 0.004837 | 0.045121 | 0.350958 | 1.601467 | 2.568476 | 1.591998 | 1.190108
0.17783 | 0.004837 | 0.045121 | 0.343803 | 1.573757 | 2.554747 | 1.590795 | 1.189698
0.31623 | 0.004837 | 0.044507 | 0.316063 | 1.47078 | 2.507772 | 1.586979 | 1.188398
0.45 | 0.004705 | 0.036765 | 0.249325 | 1.238666 | 2.421019 | 1.581228 | 1.186445
0.5 | 0.002419 | 0.022562 | 0.183937 | 1.025219 | 2.361647 | 1.578549 | 1.185538
0.56234 | 5.1e-05 | 0.006779 | 0.108887 | 0.759317 | 2.280932 | 1.5748 | 1.184271
0.75 | - | 6.4e-05 | 0.034842 | 0.416175 | 2.069946 | 1.56071 | 1.179537
1.0 | - | - | 0.006872 | 0.214491 | 1.68516 | 1.535052 | 1.171036
1.33352 | - | - | 0.000168 | 0.094966 | 1.028758 | 1.487096 | 1.155611
1.77828 | - | - | - | 0.032116 | 0.471906 | 1.391456 | 1.127131
3.16228 | - | - | - | 0.000196 | 0.049604 | 0.471468 | 0.954827
5.62341 | - | - | - | - | 0.001163 | 0.019493 | 0.082189
10.0 | - | - | - | - | - | 0.000113 | 0.000487
17.78279 | - | - | - | - | - | - | -
$\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.004837 | 0.045121 | 0.364106 | 1.512533 | 2.431781 | 1.555683 | 1.178842
0.1 | 0.004837 | 0.045121 | 0.357116 | 1.497067 | 2.426675 | 1.555215 | 1.178665
0.17783 | 0.004837 | 0.045121 | 0.343553 | 1.462337 | 2.415402 | 1.554193 | 1.178277
0.31623 | 0.004837 | 0.045121 | 0.299658 | 1.343527 | 2.378606 | 1.550952 | 1.177048
0.45 | 0.004796 | 0.033907 | 0.223636 | 1.138096 | 2.319758 | 1.546073 | 1.175203
0.5 | 0.002418 | 0.02256 | 0.183934 | 1.031348 | 2.290735 | 1.543804 | 1.174346
0.56234 | - | 0.00909 | 0.135476 | 0.89324 | 2.251138 | 1.54063 | 1.173149
0.75 | - | - | 0.043524 | 0.56287 | 2.120765 | 1.528729 | 1.16868
1.0 | - | - | 0.001805 | 0.289458 | 1.885442 | 1.507188 | 1.160663
1.33352 | - | - | - | 0.106017 | 1.298617 | 1.467469 | 1.146151
1.77828 | - | - | - | 0.018336 | 0.530154 | 1.391381 | 1.119483
3.16228 | - | - | - | - | 0.023197 | 0.536229 | 0.964849
5.62341 | - | - | - | - | 3.1e-05 | 0.005921 | 0.056159
10.0 | - | - | - | - | - | - | 2.4e-05
17.78279 | - | - | - | - | - | - | -

###  6.3 Su-Olson problem with a Gaussian source

Returning to the linearized Su-Olson problem, we consider a Gaussian source.
Here the source is defined by Eq. 37 where the uncollided solution is taken
from [15] for the full transport solution or Eq.(56) for the $S_{2}$ solution.
We set $x_{0}=0.5$ and the source duration is still $t_{0}=10$. In [13],
smooth Gaussian sources allowed for geometric convergence of the solution at
all times. We expect the same result in this application.

Since there are no discontinuities induced by nonsmoothness in the source, we
are able to employ a far simpler mesh function than what was used for the thin
square source problems. We only guess the edge of the problem domain and span
the given space evenly with stationary edges. The moving mesh was not used in
this case because earlier tests in [13] revealed the mesh to be non-useful in
smooth problems. The uncollided solution however, was employed. We include
Gaussian sources though they do not prove to be challenging enough problems to
require the full application of our method because we can achieve very
accurate solutions.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 8: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
Su-Olson Gaussian source problem where $x_{0}=0.5$, $t_{0}=10$. The quadrature
order for all results is $S_{16}$. All results were calculated with a static
mesh and uncollided source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 9: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
$S_{2}$ Su-Olson Gaussian source problem where $x_{0}=0.5$, $t_{0}=10$. All
results were calculated with a moving mesh and uncollided source treatment
except for the $t=31.6228$ and $t=100$ cases where a standard source treatment
was used. The dashed lines represent solutions found with a moving mesh.

(a) $t=0.1$

(b) $t=0.31623$

(c) $t=1$

(d) $t=3.16228$

(e) $t=10$

(f) $t=100$

Figure 10: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$)
solutions for the optically thin Su-Olson Gaussian source problem with
$x_{0}=0.5$, $t_{0}=10$. Solid lines are scalar flux, $\overline{\phi}$, and
dashed are material energy density, $\overline{e}$.

While $S_{256}$ was used for the full transport solutions for the thin square
sources, we only use $S_{64}$ on the Gaussian sources. This choice is informed
by tests run in [13] that showed that far fewer quadrature points are required
to resolve the angular error.

With the temperature defined by Eq (11) we present solutions in Tables 5 and 6
with convergence results shown in Figures 8 and 18. Solutions are shown in
Figure 10.

Table 5: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thin Gaussian source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. RMSE values included in the last row of the $S_{2}$ table are calculated from the $S_{2}$ benchmark from [6]. Convergence results for these answers are plotted in Figures 8 and 18. $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.094869 | 0.264712 | 0.571441 | 1.053107 | 1.956705 | 0.610025 | 0.316341
0.1 | 0.091217 | 0.255181 | 0.556247 | 1.033176 | 1.933017 | 0.609634 | 0.31629
0.17783 | 0.083721 | 0.235541 | 0.524582 | 0.991429 | 1.883176 | 0.608783 | 0.316178
0.31623 | 0.063836 | 0.182853 | 0.436813 | 0.873994 | 1.741091 | 0.6061 | 0.315826
0.45 | 0.042514 | 0.125118 | 0.334025 | 0.732182 | 1.564765 | 0.602104 | 0.315298
0.5 | 0.035215 | 0.104949 | 0.29568 | 0.67763 | 1.49511 | 0.600263 | 0.315054
0.56234 | 0.027081 | 0.082142 | 0.250049 | 0.611096 | 1.408369 | 0.597705 | 0.314714
0.75 | 0.010198 | 0.033042 | 0.137011 | 0.433926 | 1.163525 | 0.588305 | 0.313452
1.0 | 0.001799 | 0.006564 | 0.050023 | 0.267122 | 0.898484 | 0.572018 | 0.311224
1.33352 | 8.2e-05 | 0.000371 | 0.008965 | 0.139587 | 0.64098 | 0.544297 | 0.3073
1.77828 | - | 2e-06 | 0.000387 | 0.05756 | 0.40731 | 0.498676 | 0.300446
3.16228 | - | - | - | 0.00137 | 0.087854 | 0.327865 | 0.268818
5.62341 | - | - | - | - | 0.003363 | 0.096556 | 0.189497
10.0 | - | - | - | - | - | 0.003523 | 0.063953
17.78279 | - | - | - | - | - | - | 0.002424
$\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.00467 | 0.040073 | 0.240699 | 0.727949 | 1.572225 | 0.594815 | 0.314048
0.1 | 0.00449 | 0.03858 | 0.23398 | 0.717893 | 1.560265 | 0.594488 | 0.313999
0.17783 | 0.004119 | 0.035509 | 0.219946 | 0.696571 | 1.534799 | 0.593775 | 0.313893
0.31623 | 0.003137 | 0.027317 | 0.180837 | 0.6345 | 1.45974 | 0.591526 | 0.313557
0.45 | 0.002085 | 0.018437 | 0.134764 | 0.554569 | 1.360669 | 0.58817 | 0.313054
0.5 | 0.001726 | 0.015365 | 0.117566 | 0.522032 | 1.319396 | 0.586621 | 0.312821
0.56234 | 0.001326 | 0.011918 | 0.09716 | 0.480715 | 1.266036 | 0.584468 | 0.312496
0.75 | 0.000497 | 0.004632 | 0.047669 | 0.359391 | 1.101318 | 0.576529 | 0.311293
1.0 | 8.7e-05 | 0.000863 | 0.013277 | 0.223358 | 0.893342 | 0.562683 | 0.309166
1.33352 | 3e-06 | 4.3e-05 | 0.001308 | 0.100634 | 0.657392 | 0.538848 | 0.305419
1.77828 | - | - | 1.8e-05 | 0.022908 | 0.417809 | 0.498866 | 0.298865
3.16228 | - | - | - | - | 0.068504 | 0.34004 | 0.268475
5.62341 | - | - | - | - | 0.000121 | 0.098551 | 0.191126
10.0 | - | - | - | - | - | 0.001487 | 0.064775
17.78279 | - | - | - | - | - | - | 0.001958
$\boldsymbol{RMSE}$ | 9.105e-10 | 8.561e-10 | 2.604e-09 | 7.839e-09 | 1.321e-08 | 6.808e-09 | 6.899e-05
Table 6: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thin Gaussian source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. RMSE values included in the last row of the $S_{2}$ table are calculated from the $S_{2}$ benchmark from [6]. Convergence results for these answers are plotted in Figures 8 and 18. $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.00467 | 0.040089 | 0.245869 | 0.833299 | 1.847977 | 0.623026 | 0.318058
0.1 | 0.00449 | 0.038593 | 0.238404 | 0.815919 | 1.824555 | 0.622608 | 0.318006
0.17783 | 0.004119 | 0.035517 | 0.222913 | 0.779578 | 1.775289 | 0.621695 | 0.317892
0.31623 | 0.003137 | 0.027313 | 0.18051 | 0.677838 | 1.63499 | 0.61882 | 0.317534
0.45 | 0.002085 | 0.018426 | 0.132108 | 0.556161 | 1.461245 | 0.614538 | 0.316998
0.5 | 0.001726 | 0.015355 | 0.114504 | 0.509791 | 1.392751 | 0.612566 | 0.316749
0.56234 | 0.001326 | 0.011908 | 0.093967 | 0.453643 | 1.307592 | 0.609828 | 0.316404
0.75 | 0.000497 | 0.004629 | 0.04582 | 0.307085 | 1.068295 | 0.59977 | 0.315121
1.0 | 8.7e-05 | 0.000865 | 0.013544 | 0.175481 | 0.812011 | 0.58237 | 0.312856
1.33352 | 3e-06 | 4.4e-05 | 0.001743 | 0.082634 | 0.567439 | 0.552834 | 0.308867
1.77828 | - | - | 5e-05 | 0.029302 | 0.35102 | 0.504445 | 0.301903
3.16228 | - | - | - | 0.000297 | 0.069582 | 0.325807 | 0.269788
5.62341 | - | - | - | - | 0.002245 | 0.091971 | 0.189454
10.0 | - | - | - | - | - | 0.003086 | 0.063216
17.78279 | - | - | - | - | - | - | 0.002324
$\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.00467 | 0.040073 | 0.240699 | 0.727949 | 1.572225 | 0.594815 | 0.314048
0.1 | 0.00449 | 0.03858 | 0.23398 | 0.717893 | 1.560265 | 0.594488 | 0.313999
0.17783 | 0.004119 | 0.035509 | 0.219946 | 0.696571 | 1.534799 | 0.593775 | 0.313893
0.31623 | 0.003137 | 0.027317 | 0.180837 | 0.6345 | 1.45974 | 0.591526 | 0.313557
0.45 | 0.002085 | 0.018437 | 0.134764 | 0.554569 | 1.360669 | 0.58817 | 0.313054
0.5 | 0.001726 | 0.015365 | 0.117566 | 0.522032 | 1.319396 | 0.586621 | 0.312821
0.56234 | 0.001326 | 0.011918 | 0.09716 | 0.480715 | 1.266036 | 0.584468 | 0.312496
0.75 | 0.000497 | 0.004632 | 0.047669 | 0.359391 | 1.101318 | 0.576529 | 0.311293
1.0 | 8.7e-05 | 0.000863 | 0.013277 | 0.223358 | 0.893342 | 0.562683 | 0.309166
1.33352 | 3e-06 | 4.3e-05 | 0.001308 | 0.100634 | 0.657392 | 0.538848 | 0.305419
1.77828 | - | - | 1.8e-05 | 0.022908 | 0.417809 | 0.498866 | 0.298865
3.16228 | - | - | - | - | 0.068504 | 0.34004 | 0.268475
5.62341 | - | - | - | - | 0.000121 | 0.098551 | 0.191126
10.0 | - | - | - | - | - | 0.001487 | 0.064775
17.78279 | - | - | - | - | - | - | 0.001958
$\boldsymbol{RMSE}$ | 9.105e-10 | 8.561e-10 | 2.604e-09 | 7.839e-09 | 1.321e-08 | 6.808e-09 | 6.899e-05

As illustrated in the aforementioned convergence plots, the problem converges
geometrically even with a standard static mesh. We quickly note that for
$t=31.6228$ and $t=100$ in the $S_{2}$ results in Figure 18, a moving mesh was
employed. In this case, the mesh moved with a constant speed from the initial
width to the specified final width. This was not done out of necessity, but
rather to ascertain weather the moving mesh was useful in these problems.

On the difference between the full and $S_{2}$ solutions, we point out that it
can be seen in Figure 10 that the $S_{2}$ solution is not as accurate in
resolving the source region after the source has been on for some time then
becomes more accurate after the source is off. This is not surprising since
the source emits uncollided particles in every direction while it is on. At
later times, the two solutions are visually indistinguishable.

###  6.4 Constant $C_{\mathrm{v}}$ problem with a Gaussian source

We include a problem with the same source and parameters as the last section
but with a constant specific heat so that the temperature to material energy
density conversion is given by Eq. (12). We specify the the dimensional
specific heat to be
$C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$.
While we have less certainty in forecasting the behavior of these nonlinear
results, we still expect geometric convergence since there are no sources of
nonsmoothness.

Like the Gaussian source in the linearized system, we specified a static mesh
that evenly spans some estimated width of the actual solution domain. The
convergence results in Figures 11 and 12 show that this was sufficient to
achieve geometric convergence at all chosen times. We also see the phenomenon
first observed in Section 6.2 of extremely fast convergence for the scalar
flux at early times. This is once again the result of the scalar flux being
most uncollided at these times, making solving for the collided portion
simpler.

During the time the source is on there is a large discrepancy between the 2
and transport solutions; this discrepancy recedes after the source is turned
off. Similar to the nonlinear square source, the solution is not in
equilibrium by $t=100$. This leads us to draw the conclusion that equilibrium
is not as impacted by the nonsmoothness of the source as it is by the
functional form of the specific heat, since the Su-Olson problem with a
nonsmooth source goes more quickly into equilibrium than this constant
$C_{\mathrm{v}}$ smooth source problem. We present the solutions in Tables 7
and 8.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 11: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
constant $C_{\mathrm{v}}$ Gaussian source problem where $x_{0}=0.5$,
$t_{0}=10$. The quadrature order for all results is $S_{16}$. All results were
calculated with a static mesh and uncollided source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 12: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thin ($\sigma_{a}=1$ cm-1)
$S_{2}$ constant $C_{\mathrm{v}}$ Gaussian source problem where $x_{0}=0.5$,
$t_{0}=10$. All results were calculated with a static mesh and uncollided
source treatment.

(a) $t=0.1$

(b) $t=0.31623$

(c) $t=1$

(d) $t=3.16228$

(e) $t=10$

(f) $t=100$

Figure 13: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$) solutions for the optically thin constant $C_{\mathrm{v}}$ Gaussian source problem with $x_{0}=0.5$, $t_{0}=10$. Solid lines are radiation temperature $\overline{\phi}^{1/4}$, and dashed are temperature, $\overline{T}$. Table 7: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thin Gaussian source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.5$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 11 and 12. $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.094715 | 0.260632 | 0.501868 | 0.653649 | 1.639873 | 0.218959 | 0.069516
0.1 | 0.091069 | 0.251252 | 0.488575 | 0.636913 | 1.614763 | 0.21862 | 0.069466
0.17783 | 0.083585 | 0.231922 | 0.460879 | 0.602409 | 1.561824 | 0.217881 | 0.069356
0.31623 | 0.063732 | 0.180063 | 0.384158 | 0.509408 | 1.41001 | 0.215547 | 0.069011
0.45 | 0.042445 | 0.123229 | 0.294366 | 0.405067 | 1.219261 | 0.212063 | 0.068495
0.5 | 0.035158 | 0.103372 | 0.260874 | 0.367198 | 1.142975 | 0.210455 | 0.068256
0.56234 | 0.027037 | 0.080916 | 0.221011 | 0.322683 | 1.047037 | 0.208219 | 0.067923
0.75 | 0.010181 | 0.032561 | 0.122095 | 0.213179 | 0.768814 | 0.19997 | 0.066689
1.0 | 0.001796 | 0.006473 | 0.045352 | 0.123297 | 0.461235 | 0.185567 | 0.064514
1.33352 | 8.2e-05 | 0.000367 | 0.008352 | 0.064418 | 0.223392 | 0.160718 | 0.060695
1.77828 | - | 2e-06 | 0.000371 | 0.029558 | 0.097211 | 0.118847 | 0.054065
3.16228 | - | - | - | 0.001165 | 0.011091 | 0.014263 | 0.024308
5.62341 | - | - | - | - | 0.00034 | 0.00057 | 0.000634
10.0 | - | - | - | - | - | 4e-06 | 3e-06
17.78279 | - | - | - | - | - | - | -
x/t | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.094714 | 0.260323 | 0.477228 | 0.539628 | 1.295862 | 0.197266 | 0.066556
0.1 | 0.091067 | 0.250998 | 0.467159 | 0.53082 | 1.283232 | 0.197001 | 0.066511
0.17783 | 0.083584 | 0.231775 | 0.445829 | 0.512305 | 1.256308 | 0.196423 | 0.066411
0.31623 | 0.063733 | 0.180133 | 0.384035 | 0.459709 | 1.176709 | 0.194598 | 0.066097
0.45 | 0.042446 | 0.123414 | 0.305857 | 0.394987 | 1.071024 | 0.191874 | 0.065626
0.5 | 0.035159 | 0.103564 | 0.274806 | 0.369646 | 1.026763 | 0.190617 | 0.065409
0.56234 | 0.027038 | 0.081091 | 0.236322 | 0.338318 | 0.969319 | 0.188869 | 0.065106
0.75 | 0.010181 | 0.032613 | 0.132763 | 0.251849 | 0.790231 | 0.18242 | 0.063983
1.0 | 0.001796 | 0.006438 | 0.045179 | 0.164176 | 0.560306 | 0.171164 | 0.062005
1.33352 | 8.1e-05 | 0.000354 | 0.005814 | 0.08995 | 0.308784 | 0.151783 | 0.058534
1.77828 | - | 1e-06 | 0.000114 | 0.03286 | 0.12204 | 0.119388 | 0.052522
3.16228 | - | - | -0.0 | 2e-06 | 0.006114 | 0.019949 | 0.02605
5.62341 | - | - | - | - | 4.6e-05 | 0.000264 | 0.000564
10.0 | - | - | - | - | - | - | -
17.78279 | - | - | - | - | - | - | -
Table 8: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thin Gaussian source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.5$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 11 and 12 $\boldsymbol{x/t}$ | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
---|---|---|---|---|---|---|---
0.01 | 0.004825 | 0.044224 | 0.323377 | 1.4684 | 2.455993 | 1.518526 | 1.133064
0.1 | 0.004638 | 0.042571 | 0.313295 | 1.438115 | 2.446023 | 1.51794 | 1.132861
0.17783 | 0.004255 | 0.039172 | 0.292399 | 1.373348 | 2.424567 | 1.516657 | 1.132417
0.31623 | 0.003241 | 0.030112 | 0.235379 | 1.183116 | 2.359381 | 1.512586 | 1.131013
0.45 | 0.002154 | 0.020302 | 0.170708 | 0.944125 | 2.267966 | 1.506441 | 1.128901
0.5 | 0.001783 | 0.016913 | 0.147332 | 0.851718 | 2.227558 | 1.503576 | 1.12792
0.56234 | 0.001369 | 0.013111 | 0.120189 | 0.740141 | 2.172608 | 1.499564 | 1.126549
0.75 | 0.000513 | 0.005089 | 0.057328 | 0.458276 | 1.96946 | 1.484442 | 1.121422
1.0 | 8.9e-05 | 0.000949 | 0.016339 | 0.231272 | 1.546952 | 1.456735 | 1.112196
1.33352 | 4e-06 | 4.8e-05 | 0.002002 | 0.097311 | 0.889883 | 1.4042 | 1.095386
1.77828 | - | - | 5.5e-05 | 0.032461 | 0.400128 | 1.294119 | 1.064084
3.16228 | - | - | - | 0.000322 | 0.042204 | 0.349552 | 0.857807
5.62341 | - | - | - | - | 0.001039 | 0.014856 | 0.054666
10.0 | - | - | - | - | - | 8.8e-05 | 0.000336
17.78279 | - | - | - | - | - | - | -
x/t | 0.1 | 0.31623 | 1.0 | 3.16228 | 10.0 | 31.6228 | 100.0
0.01 | 0.004825 | 0.044206 | 0.317475 | 1.332131 | 2.309445 | 1.479441 | 1.120796
0.1 | 0.004638 | 0.042557 | 0.308253 | 1.31107 | 2.303385 | 1.478948 | 1.120606
0.17783 | 0.004255 | 0.039164 | 0.289032 | 1.265902 | 2.29029 | 1.47787 | 1.120189
0.31623 | 0.003241 | 0.030116 | 0.235788 | 1.131312 | 2.250065 | 1.474452 | 1.118869
0.45 | 0.002154 | 0.020312 | 0.173764 | 0.954656 | 2.192602 | 1.469301 | 1.116886
0.5 | 0.001783 | 0.016923 | 0.150839 | 0.882882 | 2.166882 | 1.466903 | 1.115965
0.56234 | 0.001369 | 0.013121 | 0.123828 | 0.79271 | 2.13172 | 1.463548 | 1.114678
0.75 | 0.000513 | 0.005092 | 0.059398 | 0.540414 | 2.00314 | 1.45095 | 1.109868
1.0 | 8.9e-05 | 0.000947 | 0.016015 | 0.293197 | 1.73921 | 1.428067 | 1.101228
1.33352 | 4e-06 | 4.7e-05 | 0.001515 | 0.113254 | 1.125017 | 1.385564 | 1.085536
1.77828 | - | - | 2e-05 | 0.022889 | 0.453017 | 1.302486 | 1.05652
3.16228 | - | - | - | - | 0.021263 | 0.375028 | 0.880385
5.62341 | - | - | - | - | 3.7e-05 | 0.00417 | 0.03231
10.0 | - | - | - | - | - | 1e-06 | 1.4e-05
17.78279 | - | - | - | - | - | - | -

##  7 Optically thick results

Here we include results for problems that we consider optically thick. By
this, we mean that the source width is far greater than a mean free path. We
accomplish this by specifying, $\sigma_{a}=800\>\>\mathrm{cm}^{-1}$ and
$z_{0}<1$ cm ($z_{0}$ is the dimensional $x_{0}$). In this section, we include
results for linearized Gaussian and square sources as well as a constant
$C_{\mathrm{v}}$ Gaussian source problem. We do not include a constant
$C_{\mathrm{v}}$ square source, since our method could not resolve the
nonlinear, nonequilibrium, and very sharp wave that the square source induces.

Since in an optically thick problem results are of interest only after many
mean free times, we give results for $\tau\approx 0.01,0.1,$ and $1$ ns. By
these times, the uncollided source is negligible and any discontinuous
wavefronts have decayed. For these reasons, the problems in this section do
not employ an uncollided source or a moving mesh.

###  7.1 Su-Olson problem with a square source

To keep the dimensional spatial domain manageable, we set $l=\frac{1}{800}$ in
Eqs. (5) and (6). This makes the dimensional and nondimensional domains the
same, but stiffens the system. We are again using a square source (Eq. (36))
with $x_{0}=0.5$ and $t_{0}=10$. We forgo the the use if an uncollided source
since the evaluation times of interest are long after the uncollided solution
has decayed to zero.

For the mesh in this problem, only a static mesh was necessary for
satisfactory convergence. We use a initialization outlined in Section 6.1 but
that the initial width $\delta x$ is not set to a small number, but to a guess
of the solution width at the evaluation time and the edges never move.
Essentially, the mesh is the same as the initial mesh for the thin square
source but covering a the whole domain. The Gauss-Legendre spacing of the
edges has the effect of concentrating static edges around the source edge
which makes it more likely that the region where the wavefront will be is
resolved.

The initial guess for the solution width was important and was refined with
each run increasing the number of spatial divisions. Since negative solutions
are possible in our DG formulation and more likely to occur when there is a
sharp wavefront, the temperature was calculated with
$\overline{T}=\mathrm{sign}(\overline{e})|\overline{e}|^{1/4}$.

The solution plots for this problem (Figure 16) show that for the chosen
times, the solution is in local equilibrium. Unlike the selected thin square
source solutions where a discontinuous wave travelling at the wavespeed
determines the speed the solution travels, here a wave resembling a nonlinear
heat waves moves outwards while the source is on. Geometric convergence, shown
in Figures 14 and 15, is only possible since the nonsmooth portions of the
scalar flux have decayed to zero and the solution is in equilibrium, which has
the effect of smoothing the leading edge of the wavefront.

We also take not of the similarity between the transport and $S_{2}$
solutions, apparent in the solution plots and in Tables 9 and 10. This is
expected from the optically thick results, since we saw the two solutions
converge at long times where there was equilibrium.

Table 9: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thick square source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. Convergence results for these answers are plotted in Figures 14 and 15. $\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
---|---|---|---
- | 4.999998 | 4.999998 | 4.999957
0.0579 | 4.999998 | 4.999998 | 4.999802
0.1158 | 4.999998 | 4.999998 | 4.998522
0.1737 | 4.999998 | 4.999998 | 4.991207
0.2316 | 4.999998 | 4.999998 | 4.959099
0.2895 | 4.999998 | 4.999998 | 4.850729
0.3474 | 4.999998 | 4.999957 | 4.569407
0.4053 | 4.999998 | 4.981594 | 4.007694
0.4632 | 4.997551 | 4.256673 | 3.144975
0.5211 | 0.141706 | 1.375262 | 2.125719
0.5789 | - | 0.063798 | 1.200785
0.6368 | - | 0.000274 | 0.552639
0.6947 | - | - | 0.203932
0.7526 | - | - | 0.05963
0.8105 | - | - | 0.013701
0.8684 | - | - | 0.002459
0.9263 | - | - | 0.000343
0.9842 | - | - | 3.7e-05
1.0421 | - | - | 3e-06
1.1 | - | - | -
$\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
- | 4.999999 | 5.0 | 4.999961
0.0579 | 4.999999 | 5.0 | 4.999807
0.1158 | 4.999999 | 5.0 | 4.998528
0.1737 | 4.999999 | 5.0 | 4.991218
0.2316 | 4.999999 | 5.0 | 4.959116
0.2895 | 4.999999 | 4.999999 | 4.850739
0.3474 | 4.999999 | 4.999961 | 4.56939
0.4053 | 4.999999 | 4.981698 | 4.007652
0.4632 | 4.997908 | 4.256289 | 3.144948
0.5211 | 0.141029 | 1.375703 | 2.125739
0.5789 | - | 0.063676 | 1.200832
0.6368 | - | 0.000265 | 0.552668
0.6947 | - | - | 0.203932
0.7526 | - | - | 0.059617
0.8105 | - | - | 0.013692
0.8684 | - | - | 0.002455
0.9263 | - | - | 0.000342
0.9842 | - | - | 3.7e-05
1.0421 | - | - | 3e-06
1.1 | - | - | -
Table 10: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thick square source Su-Olson problem with $x_{0}=0.5$, $t_{0}=10$. Convergence results for these answers are plotted in Figures 14 and 15. - | 0.3 | 3.0 | 30.0
---|---|---|---
- | 4.999998 | 4.999998 | 4.999957
0.0579 | 4.999998 | 4.999998 | 4.999802
0.1158 | 4.999998 | 4.999998 | 4.998522
0.1737 | 4.999998 | 4.999998 | 4.991208
0.2316 | 4.999998 | 4.999998 | 4.959105
0.2895 | 4.999998 | 4.999998 | 4.850742
0.3474 | 4.999998 | 4.999957 | 4.56943
0.4053 | 4.999998 | 4.981624 | 4.007719
0.4632 | 4.997609 | 4.256925 | 3.144988
0.5211 | 0.140405 | 1.375054 | 2.125712
0.5789 | - | 0.063722 | 1.200762
0.6368 | - | 0.000273 | 0.552615
0.6947 | - | - | 0.203916
0.7526 | - | - | 0.059622
0.8105 | - | - | 0.013699
0.8684 | - | - | 0.002458
0.9263 | - | - | 0.000343
0.9842 | - | - | 3.7e-05
1.0421 | - | - | 3e-06
1.1 | - | - | -
- | 0.3 | 3.0 | 30.0
- | 4.999999 | 5.0 | 4.999961
0.0579 | 4.999999 | 5.0 | 4.999807
0.1158 | 4.999999 | 5.0 | 4.998529
0.1737 | 4.999999 | 5.0 | 4.99122
0.2316 | 4.999999 | 5.0 | 4.959121
0.2895 | 4.999999 | 4.999999 | 4.850752
0.3474 | 4.999999 | 4.999961 | 4.569413
0.4053 | 4.999999 | 4.981729 | 4.007677
0.4632 | 4.997962 | 4.256541 | 3.144961
0.5211 | 0.139713 | 1.375495 | 2.125732
0.5789 | - | 0.063599 | 1.200809
0.6368 | - | 0.000264 | 0.552644
0.6947 | - | - | 0.203916
0.7526 | - | - | 0.05961
0.8105 | - | - | 0.013689
0.8684 | - | - | 0.002455
0.9263 | - | - | 0.000342
0.9842 | - | - | 3.6e-05
1.0421 | - | - | 3e-06
1.1 | - | - | -

.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 14: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) Su-Olson square source problem where $x_{0}=0.5$, $t_{0}=10$. The
quadrature order for all results is $S_{16}$. All results were calculated with
a static mesh and standard source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 15: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) $S_{2}$ Su-Olson square source problem where $x_{0}=0.5$, $t_{0}=10$.
All results were calculated with a static mesh and standard source treatment.

(a) $t=0.3$

(b) $t=3$

(c) $t=30$

Figure 16: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$)
solutions for the optically thick Su-Olson square source problem with
$x_{0}=0.5$, $t_{0}=10$. Solid lines are scalar flux, $\overline{\phi}$, and
dashed are material energy density, $\overline{e}$. On the scale of this
figure, the solid and dashed lines are coincident.

###  7.2 Su-Olson problem with a Gaussian source

In order to provide consistent examples across optically thin and thick
regimes and as a trial run for the constant $C_{\mathrm{v}}$ thick Gaussian of
the next section, we include here a optically thick Gaussian source. For this
and the nonlinear Gaussian, we specify the length parameter $x_{0}=0.375$. The
source duration, $t_{0}$, is still 10 and $l=\frac{1}{800}$. Once again, the
source is given by Eq. (37) and we do not use an uncollided source. Like the
thin Gaussian sources, a moving mesh is not necessary.

Figure 19 show that like the optically thick square source, the solutions are
in equilibrium during the selected time window. There is however, no wavefront
but the solution maintains Gaussian characteristics. Like the square source,
the transport and $S_{2}$ solutions are very similar. Geometric convergence of
our standard DG method is shown in Figures 17 and 18. We note that in this
problem and the Su-Olson, square source $128$ spatial cells were required to
achieve the desired rate of convergence, though the solution was smooth. This
was due to the small length scales.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 17: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) Su-Olson Gaussian source problem where $x_{0}=0.375$, $t_{0}=10$. The
quadrature order for all results is $S_{16}$. All results were calculated with
a static mesh and standard source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 18: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) $S_{2}$ Su-Olson square Gaussian problem where $x_{0}=0.375$,
$t_{0}=10$. All results were calculated with a static mesh and standard source
treatment.

(a) $t=0.3$

(b) $t=3$

(c) $t=30$

Figure 19: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$) solutions for the optically thick Su-Olson Gaussian source problem with $x_{0}=0.375$, $t_{0}=10$. Solid lines are scalar flux, $\overline{\phi}$, and dashed are material energy density, $\overline{e}$. On the scale of this figure, the solid and dashed lines are coincident. Table 11: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thick Gaussian source Su-Olson problem with $x_{0}=0.375$, $t_{0}=10$. Convergence results for these answers are plotted in Figures 17 and 9. $\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
---|---|---|---
- | 4.99565 | 4.956221 | 4.607285
0.0842 | 4.750452 | 4.71669 | 4.414224
0.1684 | 4.084736 | 4.065342 | 3.882236
0.2526 | 3.175988 | 3.173439 | 3.13421
0.3368 | 2.232954 | 2.243553 | 2.322698
0.4211 | 1.418754 | 1.435689 | 1.579266
0.5053 | 0.815508 | 0.832457 | 0.986082
0.5895 | 0.423872 | 0.437156 | 0.565183
0.6737 | 0.199217 | 0.207914 | 0.297361
0.7579 | 0.084665 | 0.089558 | 0.143615
0.8421 | 0.032536 | 0.034938 | 0.063669
0.9263 | 0.011306 | 0.012344 | 0.025911
1.0105 | 0.003552 | 0.00395 | 0.009679
1.0947 | 0.001009 | 0.001144 | 0.003319
1.1789 | 0.000259 | 0.0003 | 0.001044
1.2632 | 6e-05 | 7.1e-05 | 0.000301
1.3474 | 1.2e-05 | 1.5e-05 | 7.9e-05
1.4316 | 2e-06 | 2e-06 | 1.9e-05
1.5158 | - | - | 4e-06
1.6 | - | - | -
$\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
- | 4.995653 | 4.956229 | 4.607284
0.0842 | 4.750456 | 4.716697 | 4.414223
0.1684 | 4.084739 | 4.065349 | 3.882238
0.2526 | 3.175991 | 3.173444 | 3.134214
0.3368 | 2.232956 | 2.243557 | 2.322703
0.4211 | 1.418755 | 1.435691 | 1.57927
0.5053 | 0.815508 | 0.832459 | 0.986085
0.5895 | 0.423872 | 0.437157 | 0.565185
0.6737 | 0.199218 | 0.207914 | 0.297362
0.7579 | 0.084665 | 0.089558 | 0.143614
0.8421 | 0.032536 | 0.034938 | 0.063669
0.9263 | 0.011306 | 0.012344 | 0.02591
1.0105 | 0.003552 | 0.00395 | 0.009679
1.0947 | 0.001009 | 0.001144 | 0.003319
1.1789 | 0.000259 | 0.0003 | 0.001044
1.2632 | 6e-05 | 7.1e-05 | 0.000301
1.3474 | 1.2e-05 | 1.5e-05 | 7.9e-05
1.4316 | 2e-06 | 2e-06 | 1.9e-05
1.5158 | - | - | 4e-06
1.6 | - | - | -
Table 12: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thick Gaussian source Su-Olson problem with $x_{0}=0.375$, $t_{0}=10$. Convergence results for these answers are plotted in Figures 17 and 9. $\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
---|---|---|---
- | 4.995668 | 4.956239 | 4.6073
0.0842 | 4.750468 | 4.716705 | 4.414236
0.1684 | 4.084745 | 4.065351 | 3.882244
0.2526 | 3.17599 | 3.17344 | 3.134212
0.3368 | 2.232949 | 2.243548 | 2.322695
0.4211 | 1.418746 | 1.435681 | 1.57926
0.5053 | 0.8155 | 0.832449 | 0.986075
0.5895 | 0.423866 | 0.43715 | 0.565178
0.6737 | 0.199213 | 0.20791 | 0.297357
0.7579 | 0.084663 | 0.089556 | 0.143612
0.8421 | 0.032535 | 0.034937 | 0.063668
0.9263 | 0.011305 | 0.012343 | 0.02591
1.0105 | 0.003552 | 0.003949 | 0.009679
1.0947 | 0.001009 | 0.001144 | 0.003319
1.1789 | 0.000259 | 0.0003 | 0.001044
1.2632 | 6e-05 | 7.1e-05 | 0.000301
1.3474 | 1.2e-05 | 1.5e-05 | 7.9e-05
1.4316 | 2e-06 | 2e-06 | 1.9e-05
1.5158 | - | - | 4e-06
1.6 | - | - | -
$\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
- | 4.995672 | 4.956247 | 4.607298
0.0842 | 4.750472 | 4.716712 | 4.414236
0.1684 | 4.084748 | 4.065358 | 3.882246
0.2526 | 3.175992 | 3.173446 | 3.134216
0.3368 | 2.232951 | 2.243552 | 2.3227
0.4211 | 1.418747 | 1.435684 | 1.579264
0.5053 | 0.8155 | 0.832451 | 0.986078
0.5895 | 0.423866 | 0.437151 | 0.565179
0.6737 | 0.199214 | 0.20791 | 0.297357
0.7579 | 0.084663 | 0.089556 | 0.143611
0.8421 | 0.032535 | 0.034937 | 0.063667
0.9263 | 0.011305 | 0.012343 | 0.025909
1.0105 | 0.003552 | 0.003949 | 0.009679
1.0947 | 0.001009 | 0.001144 | 0.003319
1.1789 | 0.000259 | 0.0003 | 0.001044
1.2632 | 6e-05 | 7.1e-05 | 0.000301
1.3474 | 1.2e-05 | 1.5e-05 | 7.9e-05
1.4316 | 2e-06 | 2e-06 | 1.9e-05
1.5158 | - | - | 4e-06
1.6 | - | - | -

###  7.3 Constant $C_{\mathrm{v}}$ Gaussian problem

Finally, we provide results for the constant $C_{\mathrm{v}}$ optically thick
problem with a Gaussian source. Like the linear version of this problem,
$x_{0}=0.375$, $t_{0}=10$, and $l=\frac{1}{800}$. We choose our constant
opacity to be the same as we used for the optically thin case,
$C_{\mathrm{v0}}=0.03$. The source is again given by Eq/ (37) and the
uncollided solution is not used. Like the linear thick Gaussian, a static mesh
is employed.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 20: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) constant $C_{\mathrm{v}}$ Gaussian source problem where $x_{0}=0.375$,
$t_{0}=10$. The quadrature order for all results is $S_{16}$. All results were
calculated with a static mesh and standard source treatment.

(a) Radiation energy density, $\overline{\phi}$

(b) Material energy density, $\overline{e}$

Figure 21: Log-linear scaled average value of the solution expansion
coefficients (found by Eqs. (44)) for the optically thick ($\sigma_{a}=800$
cm-1) $S_{2}$ constant $C_{\mathrm{v}}$ Gaussian problem where $x_{0}=0.375$,
$t_{0}=10$. All results were calculated with a static mesh and standard source
treatment.

(a) $t=0.3$

(b) $t=3$

(c) $t=30$

Figure 22: $S_{2}$ (left of $x=0$) and full transport (right of $x=0$) solutions for the optically thick constant $C_{\mathrm{v}}$ Gaussian source problem with $x_{0}=0.375$, $t_{0}=10$. Solid lines are radiation temperature $\overline{\phi}^{1/4}$, and dashed are temperature, $\overline{T}$. On the scale of this figure, the solid and dashed lines are coincident. Table 13: Transport (top) and $S_{2}$ (bottom) results for the scalar flux, $\overline{\phi}$, for the thick Gaussian source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.375$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 20 and 21. $\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
---|---|---|---
- | 6.494763 | 6.374101 | 5.426524
0.0789 | 6.11596 | 6.012864 | 5.186732
0.1579 | 5.085009 | 5.026915 | 4.518059
0.2368 | 3.682628 | 3.677866 | 3.56146
0.3158 | 2.248727 | 2.285336 | 2.499432
0.3947 | 1.081952 | 1.134228 | 1.513531
0.4737 | 0.352326 | 0.389958 | 0.733651
0.5526 | 0.060476 | 0.069031 | 0.230211
0.6316 | 0.00502 | 0.005255 | 0.01702
0.7105 | 0.000253 | 0.000255 | 0.000274
0.7895 | 8e-06 | 8e-06 | 8e-06
0.8684 | - | - | -
$\boldsymbol{x/t}$ | 0.3 | 3.0 | 30.0
- | 6.494762 | 6.374098 | 5.42651
0.0789 | 6.115958 | 6.012861 | 5.18672
0.1579 | 5.085008 | 5.026914 | 4.518052
0.2368 | 3.682627 | 3.677866 | 3.561461
0.3158 | 2.248727 | 2.285337 | 2.499437
0.3947 | 1.081952 | 1.134229 | 1.513539
0.4737 | 0.352326 | 0.389959 | 0.733658
0.5526 | 0.060476 | 0.069032 | 0.230216
0.6316 | 0.00502 | 0.005255 | 0.017021
0.7105 | 0.000253 | 0.000255 | 0.000274
0.7895 | 8e-06 | 8e-06 | 8e-06
0.8684 | - | - | -
Table 14: Transport (top) and $S_{2}$ (bottom) results for the material energy density, $\overline{e}$, for the thick Gaussian source constant $C_{\mathrm{v}}$ problem with $x_{0}=0.375$, $t_{0}=10$, and $C_{\mathrm{v0}}=0.03\>\mathrm{GJ}\cdot\mathrm{cm}^{-3}\cdot\mathrm{keV}^{-1}$. Convergence results for these answers are plotted in Figures 20 and 21. x/t | 0.3 | 3.0 | 30.0
---|---|---|---
- | 3.490028 | 3.473704 | 3.33671
0.0737 | 3.444608 | 3.429805 | 3.303993
0.1474 | 3.309364 | 3.299069 | 3.20642
0.2211 | 3.086882 | 3.083965 | 3.04554
0.2947 | 2.780372 | 2.787696 | 2.82358
0.3684 | 2.389727 | 2.410536 | 2.540796
0.4421 | 1.913502 | 1.950669 | 2.194284
0.5158 | 1.363511 | 1.408883 | 1.769126
0.5895 | 0.82657 | 0.846776 | 1.213573
0.6632 | 0.436854 | 0.439279 | 0.485835
0.7368 | 0.210529 | 0.210694 | 0.212429
0.8105 | 0.093598 | 0.093606 | 0.093685
0.8842 | 0.038508 | 0.038508 | 0.038511
0.9579 | 0.014665 | 0.014665 | 0.014665
1.0316 | 0.005169 | 0.005169 | 0.005169
1.1053 | 0.001686 | 0.001686 | 0.001686
1.1789 | 0.00051 | 0.00051 | 0.00051
1.2526 | 0.000142 | 0.000142 | 0.000142
1.3263 | 3.6e-05 | 3.6e-05 | 3.6e-05
1.4 | 8e-06 | 8e-06 | 8e-06
x/t | 0.3 | 3.0 | 30.0
- | 3.490028 | 3.473704 | 3.336708
0.0737 | 3.444608 | 3.429805 | 3.303991
0.1474 | 3.309364 | 3.299069 | 3.206419
0.2211 | 3.086882 | 3.083965 | 3.04554
0.2947 | 2.780372 | 2.787696 | 2.823581
0.3684 | 2.389727 | 2.410537 | 2.540799
0.4421 | 1.913502 | 1.95067 | 2.194288
0.5158 | 1.363511 | 1.408884 | 1.769132
0.5895 | 0.82657 | 0.846774 | 1.213586
0.6632 | 0.436853 | 0.439278 | 0.485758
0.7368 | 0.210529 | 0.210693 | 0.212427
0.8105 | 0.093598 | 0.093606 | 0.093684
0.8842 | 0.038508 | 0.038508 | 0.038511
0.9579 | 0.014665 | 0.014665 | 0.014665
1.0316 | 0.005169 | 0.005169 | 0.005169
1.1053 | 0.001686 | 0.001686 | 0.001686
1.1789 | 0.00051 | 0.00051 | 0.00051
1.2526 | 0.000142 | 0.000142 | 0.000142
1.3263 | 3.6e-05 | 3.6e-05 | 3.6e-05
1.4 | 8e-06 | 8e-06 | 8e-06

Although the Gaussian is slightly misshapen in the solution plots (Figure 22)
when compared to the linearized Gaussian, the solution is smooth like the
linearized problem and spectral convergence is observed in Figures 20 and 21.
The scalar flux and material energy density are given in Tables 13 and 14.

##  8 Conclusions

We have presented benchmark solutions to time dependent radiative transfer
problems with two functional forms for the specific heat in optically thin and
thick media. These solutions will be useful to researchers who seek a
nonlinear radiative transfer benchmark for verification purposes, desire more
digits of accuracy for the Su-Olson type problem, or who intend to resolve a
mean free path in the optically thick limit with a transport code. Although
discontinuous sources are inconvenient for DG methods and required a complex
mesh treatment to find accurate results we have provided solutions for square
sources since they are simple to simulate in numerical codes and are already
implemented in codes that run the Su-Olson problem. Researchers who implement
the DG friendly Gaussian source that we have defined will be able to converge
to even more accurate results than for the square source.

To support our claim of benchmark quality results, we presented the
convergence of the magnitude of coefficients in the solution expansion for
each result. While this does not necessarily guarantee the correctness of the
system being solved, it does provide confidence that the solution is converged
to a particular value. We further increased confidence in our solutions by
running $S_{2}$ benchmarks for the linearized problems to high accuracy and
comparing to the original published results in [8]. Furthermore, the nonlinear
problems were checked for systematic errors with a $S_{n}$ code.

##  Appendix A Uncollided solutions to the $S_{2}$ transport equation

This section contains solutions to Eq. (41) for a square and a Gaussian
source, which are used for the uncollided source treatment in $S_{2}$
problems.

###  A.1 Gaussian source

With Eq. (37) as $S$ in Eq. (41), the integral evaluates to,

$\overline{\phi}_{u}^{gs}(x,t)=\frac{1}{4}\sqrt{3\pi}\sigma
e^{\frac{3\sigma^{2}}{4}-\sqrt{3}x}\left(\text{erf}\left(\frac{-3\sigma^{2}-2t+2\text{t0}+2\sqrt{3}x}{2\sqrt{3}\sigma}\right)+\right.\\\
\left.e^{2\sqrt{3}x}\left(\text{erf}\left(\frac{\sqrt{3}\sigma}{2}+\frac{t}{\sqrt{3}\sigma}+\frac{x}{\sigma}\right)-\text{erf}\left(\frac{3\sigma^{2}+2t-2\text{t0}+2\sqrt{3}x}{2\sqrt{3}\sigma}\right)\right)+\text{erf}\left(\frac{3\sqrt{3}\sigma^{2}+2\sqrt{3}t-6x}{6\sigma}\right)\right)$
(56)

Where $\sigma$ is the standard deviation.

###  A.2 Square source

To evaluate the integral, Eq. (41) for a square source, Eq. (36), each
possible case of integration limits allowed by the step functions in the
source must be considered. This will result in a piecewise function for the
uncollided scalar flux. If we define,

$F(t,\tau)=-\frac{1}{2}\exp(-t+\tau),$ (57)

then, if $t\leq t_{0}$,

$\phi^{ss}_{u}(x,t)=\\\
\begin{cases}F\bigg{|}_{\mathrm{max}(0,\tau_{a})}^{\mathrm{max}(0,\tau_{b})}&(|x|>x_{0})\\\
\\\
\begin{cases}2F\bigg{|}_{0}^{\mathrm{min}(t,t_{0})}&(t+\sqrt{3}|x|\leq\sqrt{3}x_{0})\\\
\\\
F\bigg{|}_{\tau_{a}}^{\tau_{c}}+2F\bigg{|}_{\tau_{c}}^{\mathrm{min}(t,t_{0})}&\left(t+\sqrt{3}\,|x|\geq\sqrt{3}\,x_{0}\right)\>\&\>\left(t-\sqrt{3}(|x|+x_{0})>0\right)\\\
\\\
F\bigg{|}_{0}^{\mathrm{max}(0,\tau_{c})}+2F\bigg{|}_{\mathrm{max}(0,\tau_{c})}^{\mathrm{min}(t,t_{0})}&\left(t+\sqrt{3}\,|x|\geq\sqrt{3}\,x_{0}\right)\>\&\>\left(t-\sqrt{3}(|x|+x_{0})\leq
0\right)\end{cases}&(|x|\leq x_{0})\end{cases}$ (58)

or, if $t>t_{0}$,

$\phi^{ss}_{u}(x,t)=\begin{cases}F\bigg{|}_{\mathrm{min}\left(\tau_{e},t_{0}\right)}^{\mathrm{min}\left(\tau_{d},t_{0}\right)}&x_{0}-\frac{\sqrt{3}}{3}(t-t_{0})\leq
x\leq x_{0}+\frac{\sqrt{3}}{3}(t-t_{0})\\\ \\\
F\bigg{|}_{\mathrm{min}\left(\tau_{e},t_{0}\right)}^{\mathrm{min}\left(\tau_{d},t_{0}\right)}&x>x_{0}+\frac{\sqrt{3}}{3}(t-t_{0})\\\
\\\ \begin{cases}2F\bigg{|}_{0}^{t_{0}}&t-\frac{3}{\sqrt{3}}(x_{0}-x)\leq 0\\\
\\\
2F\bigg{|}^{t_{0}}_{\tau_{d}}+F\bigg{|}_{0}^{\tau_{d}}&t-\frac{3}{\sqrt{3}}(x_{0}-x)>0\>\mathrm{and}\>t-\frac{3}{\sqrt{3}}(x_{0}-x)<t_{0}\end{cases}&x<x_{0}-\frac{\sqrt{3}}{3}(t-t_{0})\end{cases}$
(59)

Where

$\tau_{a}=t-\sqrt{3}(|x|+x_{0}),$ (60) $\tau_{b}=t-\sqrt{3}(|x|-x_{0}),$ (61)
$\tau_{c}=t-\sqrt{3}(x_{0}-|x|),$ (62)
$\tau_{d}=\left|t-\frac{3}{\sqrt{3}}(x-x_{0})\right|_{+},$ (63)

and

$\tau_{e}=\left|t-\frac{3}{\sqrt{3}}(x+x_{0})\right|_{+}.$ (64)

$|\cdot|$ returns the positive part of its argument.

## References

  * [1] J Stefan.  Uber die beziehung zwischen der warmestrahlung und der temperatur, sitzungsberichte der mathematisch-naturwissenschaftlichen classe der kaiserlichen.  Akademie der Wissenschaften, 79:S–391, 1879.
  * [2] R. E. Marshak.  Effect of radiation on shock wave behavior.  The Physics of Fluids, 1(1), 1958.
  * [3] Albert G Petschek, Ralph E Williamson, and John K Wooten Jr.  The penetration of radiation with constant driving temperature.  Technical report, Los Alamos Scientific Lab., N. Mex., 1960.
  * [4] William Bennett and Ryan G McClarren.  Self-similar solutions for high-energy density radiative transfer with separate ion and electron temperatures.  Proceedings of the Royal Society A, 477(2249):20210119, 2021.
  * [5] GC Pomraning.  The non-equilibrium marshak wave problem.  Journal of Quantitative Spectroscopy and Radiative Transfer, 21(3):249–261, 1979.
  * [6] Ryan G McClarren, James Paul Holloway, and Thomas A Brunner.  Analytic p1 solutions for time-dependent, thermal radiative transfer in several geometries.  Journal of Quantitative Spectroscopy and Radiative Transfer, 109(3):389–403, 2008.
  * [7] Barry D Ganapol and GC Pomraning.  The non-equilibrium marshak wave problem: a transport theory solution.  Journal of Quantitative Spectroscopy and Radiative Transfer, 29(4):311–320, 1983.
  * [8] Bingjing Su and Gordon L. Olson.  An analytical benchmark for non-equilibrium radiative transfer in an isotropically scattering medium.  Annals of Nuclear Energy, 24(13):1035–1055, 1997.
  * [9] Bingjing Su and Gordon L. Olson.  Non-grey benchmark results for two temperature non-equilibrium radiative transfer.  Journal of Quantitative Spectroscopy and Radiative Transfer, 62(3):279–302, 1999.
  * [10] Su Bingjing and Gordon L Olson.  Benchmark results for the non-equilibrium marshak diffusion problem.  Journal of Quantitative Spectroscopy and Radiative Transfer, 56(3):337–351, 1996.
  * [11] Ryan G McClarren and John G Wöhlbier.  Solutions for ion–electron–radiation coupling with radiation and electron diffusion.  Journal of Quantitative Spectroscopy and Radiative Transfer, 112(1):119–130, 2011.
  * [12] Ryan G. McClarren.  Two-group radiative transfer benchmarks for the non-equilibrium diffusion model.  Journal of Computational and Theoretical Transport, 50(6-7):583–597, 2022.
  * [13] William Bennett and Ryan G. McClarren.  Accurate solutions to time dependent transport problems with a moving mesh and exact uncollided source treatment, 2022.
  * [14] B. D. Ganapol, R. S. Baker, J. A. Dahl, and Raymond E. Alcouffe.  Homogeneous infinite media time-dependent analytical benchmarks.  Technical Report LA-UR-01-1854, Los Alamos National Laboratory, 2001.
  * [15] William Bennett and Ryan G McClarren.  Benchmarks for infinite medium, time dependent transport problems with isotropic scattering.  arXiv preprint arXiv:2205.15783, 2022.
  * [16] Subrahmanyan Chandrasekhar.  Radiative transfer.  Dover Publications, 1960.
  * [17] J.E. Marsden and A. Tromba.  Vector Calculus.  Number pp. 1-93 in Vector Calculus. W. H. Freeman, 2003.
  * [18] Nico Schlömer, Nick Papior, Darius Arnold, Jan Blechta, and Rasmus Zetter.  nschloe/quadpy: None, 2021.
  * [19] Pauli Virtanen, Ralf Gommers, et al.  SciPy 1.0: fundamental algorithms for scientific computing in Python.  Nature Methods, 17(3):261–272, 2020.