# Accurate solutions to time dependent transport problems with a moving mesh
and exact uncollided source treatments

William Bennett
Ryan G. McClarren

###### Abstract

To find benchmark quality solutions to time dependent $S_{n}$ transport
problems, we develop a numerical method in a Discontinuous Galerkin (DG)
framework that utilizes time dependent cell edges, called a moving mesh, and
an uncollided source treatment. The moving mesh and uncollided source
treatment is devised to circumvent discontinuities in the solution and better
realize the potential of the DG method to converge on smooth problems. The
resulting method achieves spectral convergence on smooth problems. When
applied to problems with discontinuity-inducing sources, our combined method
returns a significantly better solution than the standard DG method. On smooth
problems, we observe spectral convergence even in problems with wave fronts.
In problems where the angular flux is inherently non-smooth, we do not observe
high-order of convergence when compared with static meshes, but there is a
quantitative reduction in error of nearly three orders of magnitude.

##  1 Introduction

Implementations of the Discontinuous Galerkin method for spatial
discretization to solve the time dependent $S_{n}$ neutral particle transport
equation must reckon with the inevitable nonsmoothness due to the finite
number of wave speeds akin to contact discontinuities in fluid equations.
Either the method will fail to achieve the high-order convergence, or
computational scientists must find a way to deal with the inevitable
discontinuities that arise.

This is not a trivial undertaking. Different discontinuities are possible in
the solution and the first derivative in the angular flux and in the solution
and derivative of the scalar flux. Moreover, these discontinuities are time
dependent. Ganapol’s plane pulse benchmark solution [1], ubiquitous in the
transport literature, is a good example of the behavior of nonsmooth regions
in transport solutions. The so called uncollided flux, which refers to
particles that have not collided with other particles since being emitted from
the source, is more structured and nonsmooth. In this case, the uncollided
angular flux is a decaying delta function in each angular direction traveling
away from the initial pulse in the direction that it was emitted. The time
decay is a result of scattering, which has an overall smoothing effect on the
solution. The scalar flux solution for this problem is smooth everywhere
except for a travelling wavefront moving out from the origin at the particle
speed.

Other source configurations also exhibit the behaviors shown by the Ganapol
problem: the uncollided flux is the cause of discontinuities, scattering
smooths the solution over time, and the angular flux discontinuities can
coalesce into discontinuities in the scalar flux. Since the Ganapol problem is
a Green’s function for the time dependent transport equation, uncollided
solutions for other sources can be understood as superpositions of this
uncollided solution. This leads to another insight, if the uncollided solution
is known, then the location of any discontinuities can be determined at any
time. Also, we see from this problem that the uncollided flux in any
configuration will be the most “unsmooth” when compared to the flux from
particles that have experienced collisions. Therefore, finding the uncollided
flux has a twofold advantage. First, it gives the location of discontinuities
and second, it gives the least smooth part of the solution. Knowing the
location of discontinuities presents the opportunity to resolve them with mesh
edges and knowing the uncollided solution presents the possibility to
significantly smooth the system being solved. This is the motivation for our
method: a moving mesh to resolve discontinuities and an uncollided source
treatment to reduce nonsmoothness in the system being solved.

It is necessary to note that, as described in the plane pulse example, the
discontinuities in the angular flux are responsible for less than optimal
convergence of a DG method. Tracking and resolving each one of these
discontinuities would require a different mesh defined for each angle in the
$S_{n}$ discretization. Since this is unrealistic, especially for solutions
with large numbers of angles, we choose to only attempt to resolve
nonsmoothness in the scalar flux and abandon the hubristic quest for a
spectrally convergent method on fundamentally nonsmooth problems. The method
we propose of time dependent cell edges is similar to methods used in the
field of Computational Fluid Dynamics (CFD). For example, [2] solves an
optimization problem to find the location of discontinuities and [3] uses
shock fitting to achieve high order accuracy solutions. The principles in
those works are the same: align non-smooth features with the mesh so that
highly accurate and efficient solutions can be obtained.

Since our moving mesh is only meant to make the problem smoother inside mesh
elements, and not to resolve all discontinuities, we employ an uncollided
source treatment to make the problem smoother still. It has already been
mentioned that the solution can be split into collided and uncollided parts.
This decomposition of the solution can be extended to defining the flux based
on collisions. There is a uncollided flux (i.e., the zeroth collided flux), a
first collided flux, a second, and so on. This is called Multiple Flux
Decomposition and is used by numerical methods to distribute computational
resources more efficiently (see [4, 5, 6]). Our implementation involves using
the analytic solution for the uncollided flux as a source term to solve for
the remaining collided flux and adding the uncollided and collided solutions
at the final step. We adopt the uncollided flux solutions from [7].

In order to verify the effectiveness of the moving mesh and uncollided source
methods in solving smooth and nonsmooth problems, we chose a set of
representative sources with varying degrees of solution nonsmoothness.
Recently, we presented calculations for these sources in [7]. Also, to
quantify the effectiveness of each method on its own and in combination, we
implement four different methods. The first has both an uncollided source and
a moving mesh. One is a standard DG implementation. The final two test each
method individually: a static mesh with an uncollided source and a moving mesh
with a standard source treatment. Conveniently, the latter three methods can
be easily implemented with small modifications to the system that we derive to
find the moving mesh, uncollided solutions case.

Given that there are semi-analytic results for the problems we present here,
it may seem that computing highly accurate numerical solutions to these
problems is to carry coal to Newcastle (that is perform a pointless action).
Our ultimate goal is to develop benchmark quality solutions to nonlinear
radiative transfer without utilizing linearizations (as in the specific form
of the heat capacity widely used [8, 9, 10, 11, 12, 13] or resorting to the
equilibrium diffusion limit as the Marshak wave [14, 15, 16, 17, 18, 19, 20]).
By demonstrating that our approach is sound on problems of linear particle
transport, we build a foundation for confidence in solutions to nonlinear
problems, where there are no known semi-analytic solutions, in future work.

The uncollided solution method is presented in Section 2. The following
section contains a derivation of our DG method with moving mesh edges. The
implementation, a short description of our convergence analysis methods, and a
discussion of the results follow in Sections 4, 5, and 6 respectively.

##  2 The transport model and uncollided solution

We begin with the single-speed neutral particle transport equation in slab
geometry with isotropic scattering [1]:

$\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi=\frac{c}{2}\phi+\frac{1}{2}S(x,t,\mu).$ (1)

Here, $\psi(x,t,\mu)$ is the angular flux of particles, and
$\phi(x,t)=\int^{1}_{-1}\\!d\mu^{\prime}\,\psi(x,t,\mu^{\prime})$ is the
scalar flux given by the integral over the comment angle. The cosine of the
angle between the polar angle of the direction of travel and the $x$-axis is
represented by $\mu\in[-1,1]$. $S$ is a source and $c$ is the scattering
ratio, $c=\sigma_{\mathrm{s}}/\sigma_{\mathrm{t}}$ where $\sigma_{\mathrm{s}}$
is the scattering cross section and $\sigma_{\mathrm{t}}$ is the total
absorption plus scattering cross section. The coordinates $x$ and $t$ are
measured in units of mean free path lengths and mean free times, respectively.

We represent the angular flux as the sum of an uncollided and collided angular
flux, $\psi\equiv\psi_{\mathrm{u}}+\psi_{\mathrm{c}}$. This allows us to
obtain from Eq. (1) an expression for the the particles that have not collided
after being emitted by the source (i.e., the uncollided flux,
$\psi_{\mathrm{u}}$),

$\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi_{\mathrm{u}}=\frac{1}{2}S(x,t,\mu).$ (2)

It is possible to solve Eq. (2) analytically and recover an expression for the
uncollided flux in many cases. This is done by integrating the Green’s
solution provided by Ganapol [1] over an arbitrary source. That uncollided
Green’s function solution,
$\psi_{\mathrm{u}}=\frac{1}{2}\exp(-t)\delta\left(\mu-\frac{x}{t}\right)$, is
a piecewise constant function that has a non-zero part expanding out from the
origin at speed $\mu$. We adopt the uncollided solutions presented in [7] for
our chosen source configurations: a plane pulse, a square pulse, a Gaussian
pulse, a square source, and a Gaussian source.

With an expression for the uncollided scalar flux, the equation for the
collided flux has the same form as the original transport equation ,

$\displaystyle\left(\frac{\partial}{\partial t}+\mu\frac{\partial}{\partial
x}+1\right)\psi_{c}(x,t,\mu)$
$\displaystyle=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\left(\psi_{c}(x,t,\mu^{\prime})+\psi_{u}(x,t,\mu^{\prime})\right)$
$\displaystyle=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\psi_{c}(x,t,\mu^{\prime})+S_{u}(x,t).$
(3)

Here the uncollided source is given by

$S_{u}(x,t)=\frac{c}{2}\int^{1}_{-1}\\!d\mu^{\prime}\psi_{u}(x,t,\mu^{\prime})\\\
=\frac{c}{2}\phi_{u}(x,t).$ (4)

To find the full solution, it is necessary to sum the uncollided and collided
solutions once the collided solution has been calculated.

Since Eq. (1) and Eq. (2) have the same form, a numerical method that handles
arbitrary sources can easily be applied to both. One of the objectives of this
work is to demonstrate the efficiency improvement of analytically representing
the more discontinuous, anisotropic uncollided flux and using that to obtain
benchmark-quality solutions for the total scalar flux.

##  3 Moving mesh DG spatial discretization

The second ingredient in our benchmark solutions is the use of a moving mesh.
We develop a Discontinuous Galerkin (DG) scheme to solve equations of the form
Eq. (1). First, we use the method of discrete ordinates to approximate Eq. (1)
as a system of coupled partial differential equations,

$\left(\frac{\partial}{\partial t}+\mu_{l}\frac{\partial}{\partial
x}+1\right)\psi^{l}=\frac{c}{2}\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}\psi^{l^{\prime}}+\frac{1}{2}S(x,t,\mu)\qquad\mathrm{for}\quad
l=1\dots N,$ (5)

where the scalar flux is a weighted sum of the angular flux,
$\phi\approx\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}\psi^{l^{\prime}}$. The
number of discrete directions the angular flux is evaluated at is $N$;
$\mu_{l}$ is a discrete angle and $\psi^{l}(x,t)$ is the angular flux in the
direction $\mu_{l}$; the quadrature weights are $w_{l}$. We use standard
Gauss-Lobatto quadrature rules.

We define our solution domain on a mesh of $K$ non-overlapping cells, each
with time dependent edges $x_{L}(k,t)$ and $x_{R}(k,t)$. The left edge of the
$k^{\mathrm{th}}$ cell is always at the same position as the right edge of the
$(k-1)^{\mathrm{st}}$ cell, and so on. We define a new variable $z$ to map to
each cell to [-1,1],

$z(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.$

We define an orthonormal basis function in $z$ for each cell,

$B_{i,k}(z)=\frac{\sqrt{2i+1}}{\sqrt{x_{R}(k,t)-x_{L}(k,t)}}P_{i}(z),$ (6)

where $P_{i}$ is the $i^{th}$ Legendre Polynomial. Now we may define the weak
solution on each cell for every angle as a sum of basis functions with time
dependent coefficients,

$\psi^{l}(x,t)\approx\sum_{j=0}^{M}B_{j,k}(z)\,u^{l}_{k,j}.$ (7)

We obtain the weak formulation of Eq. (5) by multiplying by a basis function
and integrating over cell $k$,

$\underbrace{\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\frac{\partial\psi^{l}}{\partial
t}}_{\mathrm{I}}+\underbrace{\mu_{l}\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\frac{\partial\psi^{l}}{\partial
x}}_{\mathrm{II}}+\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\psi^{l}=\\\
\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\left(\frac{c}{2}\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}\psi^{l^{\prime}}+\frac{1}{2}S(x,t,\mu)\right).$
(8)

In term II, integration by parts is used to push the derivative on to the
basis function as is typical with the DG method. Since the integration domain
in term I is time dependent, we invoke a special case of the Reynolds
Transport Theorem [21] to write the total derivative as

$\frac{d}{dt}\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\psi^{l}(x,t)B_{i,k}(z)=\\\
\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\left(\frac{\partial\psi^{l}}{\partial
t}B_{i,k}(z)+\psi^{l}\frac{\partial B_{i,k}(z)}{\partial t}\right)\\\
+\frac{dx_{R}(k,t)}{dt}\psi^{l}(x_{R},t)B_{i,k}(z=1)-\frac{dx_{L}(k,t)}{dt}\psi^{l}(x_{L},t)B_{i,k}(z=-1).$
(9)

Solving Eq. (9) for term I and substituting into Eq. (8) with term II
integrated by parts, we obtain

$\frac{d}{dt}\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\psi^{l}B_{i,k}(z)+\\\
\underbrace{\frac{dx_{L}(k,t)}{dt}{\psi^{l}}(x_{L},t)B_{i,k}(z=-1)-\frac{dx_{R}(k,t)}{dt}\psi^{l}(x_{R},t)B_{i,k}(z=1)}_{III}\\\
-\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\boldsymbol{\psi_{\mathrm{c}}}^{l}\frac{dB_{i,k}(x,t)}{dt}+\underbrace{\mu_{l}{\psi^{l}}B_{i,k}(z)\bigg{|}_{x_{L}(k,t)}^{x_{R}(k,t)}}_{IV}\\\
-\mu_{l}\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\boldsymbol{\psi^{l}}^{l}\frac{dB_{i,k}(z)}{dx}+\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,\boldsymbol{\psi^{l}}^{l}B_{i,k}(z)\\\
=\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\left(\frac{c}{2}\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}\psi^{l^{\prime}}+\frac{1}{2}S(x,t,\mu)\right).$
(10)

III and IV both involve evaluating the solution at the edges of the cell and
are combined to create a numerical flux (not to be confused with the scalar
flux or the angular flux) term that governs the flow of information based on
the speed of the particles relative to the mesh
($\boldsymbol{LU}^{\mathrm{surf}}$). Then, substituting Eq. (7) into Eq. (10)
and exploiting the orthonormality of the chosen basis functions to simplify
the mass matrix and the scalar flux term, we obtain

$\frac{d\boldsymbol{U}_{l}}{dt}-\boldsymbol{\underline{\underline{G}}U}_{l}+\left(\boldsymbol{LU}_{l}\right)^{(\mathrm{surf})}-\mu_{l}\boldsymbol{\underline{\underline{L}}U}_{l}+\boldsymbol{U}_{l}=\frac{c}{2}\sum_{l^{\prime}=1}^{N}w_{l^{\prime}}\boldsymbol{U}_{l^{\prime}}+\frac{1}{2}\boldsymbol{Q},$
(11)

where the time dependent solution vector is

$\boldsymbol{U}_{l,k}=[u^{l}_{k,0},u^{l}_{k,1},...,u^{l}_{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}(z)\,\frac{dB_{i,k}(z)}{dx},$
(12) $G_{i,j}=\int^{x_{R}}_{x_{L}}\\!dx\,B_{j,k}(z)\,\frac{dB_{i}(z)}{dt},$
(13) $Q_{i}=\int^{x_{R}}_{x_{L}}\\!dx\,B_{i,k}(z)\,S(x,t,\mu),$ (14)

and

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

$\boldsymbol{\psi^{l+}}$ and $\boldsymbol{\psi^{l-}}$ are found by evaluating
Eq. (7) with an upwinding scheme relative to the mesh motion at the right and
left cell edges respectively.

The initial condition is found from

$u_{k,j}^{l}=\int^{x_{R}(k,0)}_{x_{L}(k,0)}\\!dx\,B_{i,k}(z)\psi(x,t=0,\mu_{l}).$
(16)

Equation (11) is a system of coupled ordinary differential equations for the
solution in an given cell that requires a time integration algorithm to
update. To capture particles traveling with the wavefront, a Gauss-Lobatto
quadrature scheme that includes the endpoints [-1,1] was used to calculate the
angles and the weights [22]. These we calculated with the Python package
quadpy [23]. This scheme easily handles arbitrary sources, which is useful
since the uncollided solutions we use as source terms are usually complicated
functions of space and time. To apply this scheme on a moving mesh, it is
necessary to specify the velocities of the cell edges. We detail our approach
to this in Section 4.

###  3.1 $M=1$ Example equations

To illustrate the method in a more tangible way, we show the method for $M=1$.
This choice makes Eq. (11) for an angle, $l$ and cell, $k$, have the form

$\frac{d}{dt}\left(\begin{array}[]{c}u^{l}_{k,0}\\\ \\\ u^{l}_{k,1}\\\
\end{array}\right)-\underbrace{\left(\begin{array}[]{cc}-\frac{\frac{dx_{R}}{dt}-\frac{dx_{L}}{dt}}{2(x_{R}(t)-x_{L}(t))}&0\\\
\frac{\sqrt{3}\left(\frac{dx_{L}}{dt}+\frac{dx_{R}}{dt}\right)}{x_{L}(t)-x_{R}(t)}&\frac{3\left(\frac{dx_{R}}{dt}-\frac{dx_{L}}{dt}\right)}{2(x_{L}(t)-x_{R}(t))}\\\
\end{array}\right)}_{\boldsymbol{G}}\boldsymbol{\cdot}\left(\begin{array}[]{c}u^{l}_{k,0}\\\
\\\ u^{l}_{k,1}\\\
\end{array}\right)-\mu_{l}\underbrace{\left(\begin{array}[]{cc}0&0\\\
-\frac{2\sqrt{3}}{x_{L}(t)-x_{R}(t)}&0\\\
\end{array}\right)}_{\boldsymbol{L}}\boldsymbol{\cdot}\left(\begin{array}[]{c}u^{l}_{k,0}\\\
\\\ u^{l}_{k,1}\\\ \end{array}\right)\\\ \\\
+\underbrace{\left(\begin{array}[]{c}u^{l}_{k,0}\\\ \\\ u^{l}_{k,1}\\\
\end{array}\right)}_{\boldsymbol{U}^{l}}+\underbrace{\left(\begin{array}[]{c}\left(\mu_{l}-\frac{dx_{R}}{dt}\right)\frac{1}{\sqrt{x_{R}(t)-x_{L}(t)}}\boldsymbol{\psi^{l+}}-\left(\mu_{l}-\frac{dx_{L}}{dt}\right)\frac{1}{\sqrt{x_{R}(t)-x_{L}(t)}}\boldsymbol{\psi^{l-}}\\\
\\\
\left(\mu_{l}-\frac{dx_{R}}{dt}\right)\frac{2}{\sqrt{x_{R}(t)-x_{L}(t)}}\boldsymbol{\psi^{l+}}+\left(\mu_{l}-\frac{dx_{L}}{dt}\right)\frac{2}{\sqrt{x_{R}(t)-x_{L}(t)}}\boldsymbol{\psi^{l-}}\\\
\end{array}\right)}_{\boldsymbol{LU}^{\mathrm{surf}}}\\\ \\\
=\frac{c}{2}\left(w_{0}\left(\begin{array}[]{c}u^{0}_{k,0}\\\ \\\
u^{l}_{k,1}\\\ \end{array}\right)+w_{1}\left(\begin{array}[]{c}u^{1}_{k,0}\\\
\\\ u^{1}_{k,1}\\\
\end{array}\right)+\dots+w_{N}\left(\begin{array}[]{c}u^{N}_{k,0}\\\ \\\
u^{N}_{k,1}\\\
\end{array}\right)\right)+\underbrace{\left(\begin{array}[]{c}\int^{x_{R}}_{x_{L}}\\!dx\,B_{0,k}(z)\,S(x,t,\mu)\\\
\\\ \int^{x_{R}}_{x_{L}}\\!dx\,B_{1,k}(z)\,S(x,t,\mu)\\\
\end{array}\right)}_{\boldsymbol{Q}}.$ (17)

The terms $(\mu_{l}-\frac{dx_{R}}{dt})$ and $(\mu_{l}-\frac{dx_{L}}{dt})$ give
the particle velocity relative to the right and left mesh edges respectively
and determine where the solution is evaluated in the upwinding scheme. If
$(\mu_{l}-\frac{dx_{R}}{dt})>0$,

$\boldsymbol{\psi^{l+}}=\sum_{j=0}^{M}B_{j,k}(z)\,u_{k,j}^{l}.$ (18)

If the relative velocity is negative, the right edge solution is evaluated in
the next cell,

$\boldsymbol{\psi^{l+}}=\sum_{j=0}^{M}B_{j,k+1}(z)\,u_{k+1,j}^{l}.$ (19)

For the solution at the left edge, if $(\mu_{l}-\frac{dx_{L}}{dt})>0$, the
solution from the previous cell is used.

$\boldsymbol{\psi^{l-}}=\sum_{j=0}^{M}B_{j,k-1}(z)\,u_{k-1,j}^{l}.$ (20)

For a positive relative velocity,

$\boldsymbol{\psi^{l-}}=\sum_{j=0}^{M}B_{j,k}(z)\,u_{k,j}^{l}.$ (21)

When dealing with edges at the end of the domain, $k=1$ or $k=K$, the boundary
condition is required to evaluate Eq. (19) and Eq. (20). For the infinite
medium problems we explored, the boundary is zero with the exception of the
MMS problem (Section 6.1).

##  4 Implementation

To verify the effectiveness of our moving mesh DG scheme with an uncollided
source treatment, we implement the method in a code written in Python with
Numba [24] to solve transport problems with six representative sources: a
Method of Manufactured Solutions (MMS) source, a Gaussian pulse and source, a
square pulse and source, and a plane pulse. For error quantification, which is
addressed in the next Section, the full solution for the MMS source is known
from the problem setup, and we adopt the full solutions from [7] to calculate
the accuracy of our method in the remaining five source configurations. All of
our solutions are available on
Github111www.github.com/wbennett39/moving_mesh_radiative_transfer.

Since we intend to compare the uncollided source treatment and the moving mesh
individually and in unison, it is necessary to develop a code that readily
switches between four different methods: (1) a moving mesh using an uncollided
source, (2) a moving mesh with out an uncollided source, a static mesh with
(3) and without (4) the uncollided source. Each of these cases requires
solving the system of ODE’s given by Eq. (11) with different functions for the
source term and the time dependent mesh edges.

In this section, the methods we used to describe the moving and static meshes
and the source treatments are presented in a general way. The subsequent
results section includes the particular implementations for each source. Also
in this section is a discussion of the time integration method.

###  4.1 Moving mesh

While our method admits mesh motion that is any function of time that does not
result in mesh edges crossing or zero width cells, we have chosen to restrict
our investigation to one simple method for moving the cells. This simple
method requires the mesh to be subdivided into an even number of cells. The
initial mesh spans a finite width of $2x_{0}$, centered on zero. The edges
move with a constant velocity away from the origin that is dependent on their
initial location. If the edges are defined as a vector of location values
where $x_{0}$ is the leftmost cell edge, $x_{1}$ is the right edge of that
cell or the left edge of the adjacent cell, etc., then

$\boldsymbol{X}(t)=\left[x_{0}(t),x_{1}(t),...,x_{K}(t)\right],$ (22)

and the location of the edges can be found with,

$x_{k}(t)=x_{k}(0)+vt\frac{x_{k}(0)}{x_{K}(0)},$ (23)

where $x_{k}(0)$ is the initial location of a given edge and $v$, the particle
velocity, is unity.

If the initial widths are chosen to span a finite source, Eq. (23) moves the
outermost edges at speed one, matching the solution wavefront. Since a cell
interface will be initialized at $x=0$, that interface will never move. For
our static mesh calculations, we simply span a chosen width with evenly-spaced
cells and set $v=0$.

###  4.2 Source treatment

The pulsed sources, i.e., the Gaussian pulse, the square pulse, and the plane
pulse, are equivalent to initial conditions for Eq. (1). For these cases in
the standard source treatments that do not employ the method of uncollided
solutions, the source $S$ in Eq. (14) is set to zero and the initial condition
is found by letting $\psi(x,t=0)=S(x,t=0)/2$ and inserting into Eq. (16). For
the uncollided source treatment of the pulsed sources, the initial condition
is set to zero and the uncollided solution is used in Eq. (14) to find the
source in the weak formulation, $\boldsymbol{Q}$. Similarly for the source
cases (Gaussian source, square source), $\boldsymbol{Q}$ is found by
integrating either the uncollided solution or the source term, depending on
the method used.

###  4.3 Time integration

To solve the system of coupled ordinary differential equations from Eq. (11),
we employed an $8^{\mathrm{th}}$ order, explicit Runge-Kudda algorithm, DOP853
[25], as implemented in scipy [26]. The relative tolerance parameter was set
to $5\times 10^{-13}$ and the absolute tolerance to $10^{-12}$.

##  5 Error characterization

We judge the effectiveness of our scheme by characterizing the convergence of
the solution on a test problem. We use the root mean squared error (RMSE) of
the computed scalar flux as our error metric,

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

where $\phi_{i}$ is the calculated scalar flux at a given node,
$\hat{\phi_{i}}$ is the corresponding benchmark solution, and $N$ is the total
number of nodes in the computational solution. To characterize the
convergence, we solve a benchmark problem and increase the degrees of freedom.
For a problem that is algebraically convergent, holding the number of basis
functions constant and increasing the cell divisions leads to an error
behavior that limits, as $K\rightarrow\infty$, to the form

$\mathrm{RMSE}=C\,K^{-A},$ (25)

where $C$ is the $y$ intercept, and the constant $A$ is the rate of
convergence. If $A$ is $2$, the method is said to converge at second order.
The curve in Eq. (25) is a straight line on a graph where both axes have a
logarithmic scale.

In characterizing problems where the method shows algebraic convergence, the
intercept is significant. Two separate methods may have the same algebraic
convergence order, but wildly different errors when tried on a problem if one
method has a smaller intercept value. Therefore, we use our data for the
errors as a function of the number of cells, $K$ to estimate the values of $A$
and $C$ based on Eq. (25).

For problems that demonstrate geometric spectral convergence, the error can be
modeled as

$\mathrm{RMSE}=C\,\exp(-c_{1}M),$ (26)

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 spectral problems,
the coefficient $C$ is less consequential as our results show that the error
is exceedingly small for modest values of $M$.

##  6 Results

Table 1: Parameters for each test case in section 6

.  Section Problem Initial Condition Source Uncollided source $c$ $x_{0}$
$t_{0}$ $\sigma$ 6.1 MMS
$\psi(x,t=0,\mu)=\frac{e^{-x^{2}/2}}{2}\Theta(t-|x|+x_{0})$
$S(x,t,\mu)=-\frac{e^{-\frac{x^{2}}{2}}(\mu,(t+1)x+1)}{(t+1)^{2}}\Theta(x-|t|+x_{0})$
- 1.0 $0.1$ - - 6.2 Gaussian pulse
$\psi(x,t=0,\mu)=\frac{1}{2}\exp\left(\frac{-x^{2}}{\sigma^{2}}\right)$ - Eq.
(31) 1.0 - - 0.5 6.3 Gaussian source $\psi(x,t=0,\mu)=0$
$S(x,t)=\exp\left({\frac{-x^{2}}{\sigma^{2}}}\right)\,\Theta(t_{0}-t)$ Eq.
(33) 1.0 - $5.0$ 0.5 6.4 Plane pulse
$\psi(x,t=0,\mu)=\frac{1}{2}\delta(x)\delta(t)$ - Eq. (35) 1.0 0.5 - - 6.5
Square pulse $\psi(x,t=0,\mu)=\frac{1}{2}\Theta\left(x_{0}-|x|\right)$ - Eq.
(39) 1.0 0.5 - - 6.6 Square Source $\psi(x,t=0,\mu)=0$
$S(x,t)=\Theta(x_{0}-|x|)\Theta(t_{0}-t)$ Eq. (43) 1.0 0.5 $5.0$ - 6.7 Square
pulse $\psi(x,t=0,\mu)=\frac{1}{2}\Theta\left(x_{0}-|x|\right)$ - Eq. (39)
0.8, 1.2 $0.625$, $0.417$ - - 6.7 Gaussian pulse
$\psi(x,t=0,\mu)=\frac{1}{2}\exp\left(\frac{-x^{2}}{\sigma^{2}}\right)$ - Eq.
(31) 0.8, 1.2 - - $0.625$, $0.417$

###  6.1 MMS

In the Method of Manufactured solutions, a solution is chosen and inserted
into the governing equations to solve for a source term [27, 28, 29]. This
source is then used in a numerical implementation to converge to the already
known solution. Here, we specify a solution that will mimic the behavior of a
plane pulse source, where a discontinuous wavefront smooths into a solution
with Gaussian characteristics. The solution is:

$\psi_{\mathrm{MMS}}(x,\mu,t)=\frac{e^{-x^{2}/2}}{2(1+t)}\Theta(t-|x|+x_{0}).$
(27)

where $\Theta$ is a step function. Notice that $\psi_{\mathrm{MMS}}$ does not
depend on $\mu$ so that the quadrature order will not affect the numerical
solution. For three times, $t=1$, $t=5$, and $t=10$ with $x_{0}=0.1$, the
solution is plotted in Figure 1. The source term that yields this manufactured
solution is

$S_{\mathrm{MMS}}=-\frac{e^{-\frac{x^{2}}{2}}(\mu(t+1)x+1)}{(t+1)^{2}}\Theta(x-|t|+x_{0}).$
(28)

The source for Eq. (11) is found by inserting Eq. (28) into Eq. (14).

It is important to note that the wavefront in this case is different than the
wavefronts that appear in the solutions for the finite width sources (plane
pulse, square pulse and source). In this MMS problem, the wavefront is not a
feature of the finite wavespeed in the governing equation, but is instead
imposed by the step function. Therefore, it is necessary in this problem, and
this problem only, to impose a boundary condition at the wavefront, $x=\pm
t\pm x_{0}$. The value of the solution vector at the edges is found by
integrating

$u_{k,j}^{l}=\int^{x_{R}(k,t)}_{x_{L}(k,t)}\\!dx\,B_{i,k}(z)\psi_{\mathrm{MMS}}(x,t,\mu_{l}),$
(29)

at the wavefront. The initial condition is found from Eq. (16).

The moving mesh is ideal for enforcing this boundary condition since a cell
edge can always be matched to the travelling wavefront. While it is possible
to enforce the boundary condition, Eq. (29), with a static mesh, we implement
this MMS problem not to compare the different schemes but to test the rate of
convergence of the moving mesh method on a problem with a known smooth
solution. We also do not attempt to find the uncollided solution in this
configuration. The moving mesh cell edges are governed by Eq. (23).

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 1: MMS analytic solution, $\phi$, with $x_{0}=0.1$ at times $t=1,5,$
and $10$.

The MMS problem is the only problem we consider that we are assured the
angular flux solution is a smooth function. The Gaussian source and pulse
considered later have smooth source functions, but there is no analytic
expression for the full solution. Figure 2 shows linear convergence on a log-
linear scale which, as explained in Section 5, indicates spectral convergence
(Eq. (26)).

Figure 2: MMS problem convergence results on a logarithmic-linear scale with
increasing number of basis functions for $t=1$ and $4$ cells in the moving
mesh. The uncollided source treatment is not used in this problem. Using Eq.
(26), the decay rate in $M$ is $c_{1}\approx 1.3$.

If we fix $M$, we expect the error to converge at an algebraic rate equal to
the number of basis functions $(M+1)$. In Figure 3 we observe that this is
indeed the case for the MMS problem. Also, this figure shows the phenomena
which will be repeated in subsequent tests of a lowering of the $y$ intercept
with increasing number of basis functions. This intercept is $C$ in Eq. (25).
Although the MMS source is anisotropic, this solution required only $S_{32}$
to converge to machine precision. The MMS solution results do not necessarily
prove the merit of the moving mesh since the wavefront is created in this
problem with a boundary condition and says nothing about uncollided source
treatments. However, the MMS test verifies that this method can converge at
optimal rates when applied to a smooth problem.

Figure 3: MMS problem convergence results for a moving mesh, standard source
treatment

on a logarithmic scale at $t=1$ for $M=2$ (circles), $M=4$ (triangles), and
$M=6$ (squares).

###  6.2 Gaussian pulse

For the next problem we consider a Gaussian pulse of the form,

$S_{\mathrm{gp}}(x,t)=\exp\left({\frac{-x^{2}}{\sigma^{2}}}\right)\,\delta(t),$
(30)

where $\sigma$, not to be confused with a cross section, is the standard
deviation. The solution for $\sigma=0.5$ is plotted in Figure 4. The
uncollided solution is significant at early times, but has decayed to
approximately zero by $t=5$. Unlike the finite width sources examined later,
the uncollided solution for this source is a smooth function.

The solution for the uncollided scalar flux in this configuration is [7],

$\phi_{u}^{\mathrm{gs}}(x,t)=\sigma\,\sqrt{\pi}\,e^{-t}\,\frac{\text{erf}\left(\frac{t-x}{\sigma}\right)+\text{erf}\left(\frac{t+x}{\sigma}\right)}{4t}.$
(31)

For the uncollided case, Eq. (31) is treated as $S$ in Eq. (14) and the
initial condition is zero everywhere. For the methods that do not use an
uncollided source, the initial condition is found by inserting
$\psi(x,t=0)=\frac{1}{2}\exp\left({\frac{-x^{2}}{\sigma^{2}}}\right)$ into Eq.
(16) and setting $S(x,t)=0$.

While the initial condition does have an infinite width, setting the initial
mesh width ($x_{0}$) to be the location where the initial condition is below a
specified small tolerance in Eq. (23) is sufficient to achieve accurate
results. For example, the initial condition is less than $1\times 10^{-16}$ at
$x=3.1$ with the standard deviation, $\sigma=0.5$. We use the same $x_{0}$ for
static mesh calculations and initialize the mesh to span the space
$[-t_{\mathrm{final}}-x_{0},t_{\mathrm{final}}+x_{0}]$.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 4: Gaussian pulse semi-analytic solution, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed) with $\sigma=0.5$ and $c=1$. The uncollided
solution is not shown for times where it is negligible.

The Gaussian pulse and Gaussian source are a more realistic test than the MMS
problem and still have the guarantee of a smooth source and smooth uncollided
solution. All combinations of moving or static mesh and uncollided source or
standard source that we applied to this problem achieved spectral convergence,
shown in Figure 5. The uncollided, moving mesh method performed the best,
consistently achieving lower error than the other methods.

For the Gaussian pulse, only $S_{256}$ is required to achieve accuracies of
$\mathrm{RMSE}\approx 10^{-6}$. After that point, more angles are required to
ensure that the angular discretization error is lower than the errors in
discretizing space or time.

Figure 6(b) demonstrates that each method can achieve sixth-order convergence
by holding $M=6$ and increasing the number of mesh subdivisions, $K$. Since
our methods are spectrally convergent, with seven basis functions ($M=6$) one
would expect the convergence to approach seventh-order as
$K\rightarrow\infty$. This convergence test is limited by the accuracy of the
benchmark solution before $K$ is sufficiently large to achieve seventh order
convergence. However, Figure 6(a) achieves fourth-order convergence with four
basis functions ($M=3$) before reaching the accuracy limit of the benchmark.

Figure 5: Gaussian pulse convergence with logarithmic- linear scaling at $t=1$
increasing number of basis functions with $4$ cells. Blue lines indicate the
uncollided solution is used, red that no uncollided source is used. Dashed
lines are for a static mesh and solid lines are for the moving mesh. Using Eq.
(26), the estimated decay rate in $M$ is $c_{1}\approx 1.0$.

(a) $M=3$

(b) $M=6$

Figure 6: Gaussian pulse convergence results on a logarithmic scale with
$c=1$, $\sigma=0.5$ at $t=1$. Blue lines indicate the uncollided solution is
used, red that no uncollided source is used. Dashed lines are for a static
mesh and solid lines are for the moving mesh.

###  6.3 Gaussian source

A Gaussian source turned on at $t=0$ and left on until $t=t_{0}$ is a
superposition of Gaussian pulses,

$S_{\mathrm{gs}}(x,t)=\exp\left({\frac{-x^{2}}{\sigma^{2}}}\right)\,\Theta(t_{0}-t).$
(32)

Figure 7 shows this solution with $\sigma=0.5$ and $t_{0}=5$. Notice that the
uncollided solution is still significant at $t=5$. Like the Gaussian pulse,
the uncollided angular flux and the collided angular flux are smooth
functions.

Reference [7] gives the uncollided scalar flux for this source as a
convolution of the Gaussian pulse uncollided flux,

$\phi_{u}^{\mathrm{gs}}(x,t)=\int_{0}^{\mathrm{min}(t,t_{0})}\\!d\tau\,\sigma\,\sqrt{\pi}\,e^{-(t-\tau)}\,\frac{\text{erf}\left(\frac{t-\tau-x}{\sigma}\right)+\text{erf}\left(\frac{t-\tau+x}{\sigma}\right)}{4(t-\tau)}.$
(33)

For the two solution methods that use an uncollided source treatment, Eq. (33)
is integrated in Eq. (14). For the two cases that do not employ the uncollided
source, $S$ in Eq. (14) is Eq. (32). The angular flux is initialized to be
zero everywhere for both source treatments.

The moving mesh and static mesh are treated in the exact same way as Section
6.2, where the initial width is set to be the initial width of a Gaussian
pulse with the same standard deviation.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 7: Gaussian source semi-analytic solution, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed) with $\sigma=0.5$, $t_{0}=5$, and $c=1$. The
uncollided flux is not shown for times where it is negligible.

Our solutions for the Gaussian source achieve similar levels of accuracy as in
the Gaussian pulse. Figure 8 shows spectral convergence for all four methods.
Compared with the Gaussian pulse (Figure 5), these results show a more drastic
difference in the intercepts between the methods that use the uncollided
solution and those that do not. This could be because, as shown in Figure 7,
the uncollided solution is a significant portion of the solution at later
times.

Figure 9(b) shows these methods achieving almost sixth order convergence with
$M=6$ and improving to seventh order for $K=16$ before the error approaches
the number of accurate digits for the benchmark solution. Figure 9(a) better
illustrates how each method can achieve an order of convergence equal to the
number of basis functions ($M+1$) as $K\rightarrow\infty$ with fourth-order
convergence with $M=3$. These solutions have an angular error requiring
$S_{512}$ to achieve $\mathrm{RMSE}\approx 10^{-7}$.

Figure 8: Gaussian source convergence on a logarithmic-linear scale at $t=1$
increasing number of basis functions with $4$ cells. Blue lines indicate the
uncollided solution is used, red that no uncollided source is used. Dashed
lines are for a static mesh and solid lines are for the moving mesh. Using Eq.
(26), the decay rate in $M$ is estimate to be $c_{1}\approx 0.95$

(a) $M=3$

(b) $M=6$

Figure 9: Gaussian source convergence results on a logarithmic scale with
$c=1$, $\sigma=0.5$ at $t=1$. Blue lines indicate the uncollided solution is
used, red that no uncollided source is used. Dashed lines are for a static
mesh and solid lines are for the moving mesh.

###  6.4 Plane pulse

Ganapol’s frequently used benchmark solution for an infinite plane source is a
Green’s function for the source,

$S_{\mathrm{pl}}=\delta(x)\delta(t).$ (34)

The full solution is plotted for three times in Figure 10. This solution is
shows how a discontinuous uncollided scalar flux can cause discontinuities in
the first derivative of the full solution. In this case, the discontinuity is
manifested as traveling wavefront at early time. The solution smooths to be
redolent of a Gaussian at later times.

The uncollided solution for this configuration, also given by [1], is

$\phi_{\mathrm{u}}^{\mathrm{pl}}(x,t)=\frac{\exp\left(-t\right)}{2t}\Theta\left(1-\left|\frac{x}{t}\right|\right).$
(35)

Equation (35) is substituted for $S$ in Eq. (14) for the uncollided source
treatment. For the two cases where an uncollided source is not used, it is
necessary to approximate the delta function initial condition. For the static
mesh case, an initial condition of,

$\psi^{\mathrm{pl}}(x,t=0)=\frac{1}{2}\frac{\Theta(x_{0}-|x|)}{2x_{0}},$ (36)

where $x_{0}$ is a very small number, is used to approximate the initial
condition in Eq. (16). The moving mesh, no uncollided source case did not
converge due to how we represent the initial condition and is not included
here.

The moving mesh is governed by Eq. (23) with $x_{0}$ set to approximately
zero, which ensures that the step function in Eq. (35) is always one. The
static mesh for this source spans $[-t_{\mathrm{final}},t_{\mathrm{final}}]$
with evenly spaced cells.

For $j>0$, our $B_{j}$ basis are orthogonal to functions that are constant in
the range $x_{L}$ to $x_{R}$. The uncollided solution always satisfies this
criterion since the wavefronts are never inside any of the cells. Therefore,
it is simple to integrate the uncollided solution for the moving mesh over a
cell in the source term. The source becomes,

$Q_{0}(x,t)=\int^{x_{R}}_{x_{L}}\\!dx\,B_{0,k}(z)\,\frac{\exp\left(-t\right)}{2t}=\sqrt{x_{R}-x_{L}}\,\frac{\exp\left(-t\right)}{2t},$
(37)

This simplification is not possible in the static mesh case due to the step
function, which requires integration at the moving wavefront.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 10: Plane pulse semi-analytic solution, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed) with $c=1$. The uncollided flux is not shown for
times where it is negligible.

For the finite width sources, the plane pulse and the square source and pulse,
spectral convergence is unrealistic with our methods. The discontinuities in
the scalar flux can be easily matched with mesh edges, resulting in a smooth
problem, but the angular flux poses a more difficult problem. In the angular
flux of the plane pulse for example, the angular collided flux has a
discontinuity in the first derivative travelling outward from the origin at
speed $\mu t$ that is caused by the uncollided angular flux. This means that
there is as many discontinuities in the angular fluxes as there are discrete
angles.

Though the highly nonsmooth nature of the solution to the plane pulse at early
time restricts all methods in Figure 11 to first order convergence, there is a
substantial difference in the accuracy of each method. Table 2 shows that in
the $M=6$ case, there is an almost $900$ times reduction in the intercept of
the uncollided source, moving mesh case compared to the standard DG
implementation. These results show that using an uncollided source is
effective in reducing the error, in this case $36$ times, but using the
uncollided source and the moving mesh is a significantly more effective
method.

The Gaussian nature of the plane pulse solution at later times allows for
higher order convergence, as shown in Figure 12. The no-uncollided solution
method in this case is restricted by the inherent error in approximating a
delta function and the uncollided source, static mesh method is restricted by
the necessity to resolve the discontinuous wavefront at early times with a
polynomial. The uncollided, moving mesh method shows nearly optimal, sixth-
order convergence with a significantly lower error level.

The plane pulse solution requires more angular resolution than all of the
other source configurations tested, with $S_{4096}$ required for errors less
than $10^{-4}$. This angular error is reduced as the solution smooths at later
times.

(a) $M=4$

(b) $M=6$

Figure 11: Plane pulse convergence results on a logarithmic scale with $c=1$ at $t=1$. Blue lines indicate the uncollided solution is used, red that no uncollided source is used. Dashed lines are for a static mesh and solid lines are for the moving mesh. Table 2: Plane pulse intercepts from Figure (11) where the intercept, is found from the curve fit, $\mathrm{RMSE}=C\,K^{-A}$, where $K$ is the number of mesh subdivisions and $A$ is the order of convergence. The improvement from the baseline is found by dividing the intercept from the no uncollided solution, static mesh case by the intercept from the given case. | M = 4 | M = 6
---|---|---
| intercept | improvement | intercept | improvement
no uncollided $+$ static mesh | 0.41393 | – | 0.4140 | –
uncollided $+$ static mesh | 0.01848 | 22 | 0.01145 | 36
uncollided $+$ moving mesh | 0.0008613 | 481 | 0.0004731 | 875
Figure 12: Plane pulse convergence results on a logarithmic scale with $c=1$
at $t=10$ with $M=6$ . Blue lines indicate the uncollided solution is used,
dashed that the mesh is static, and solid that the mesh is moving.

###  6.5 Square pulse

We next consider a finite width square pulse source of the form,

$S_{sp}=\Theta(x_{0}-|x|)\delta(t).$ (38)

[7] gives the piecewise linear uncollided solution for this problem,

$\phi^{\mathrm{sp}}_{u}(x,t)=\begin{cases}0&|x|-t>x_{0}\\\
x_{0}\frac{e^{-t}}{t}&t>x_{0}\>\&\>x_{0}-t\leq x\leq t-x_{0}\\\ e^{-t}&t\leq
x_{0}\>\&\>t-x_{0}\leq x\leq x_{0}-t\\\
\frac{e^{-t}(t+x+x_{0})}{2t}&-t-x_{0}<x<t+x_{0}\>\&\>x_{0}+x\leq t\leq
x_{0}-x\\\ \frac{e^{-t}(t-x+x_{0})}{2t}&-t-x_{0}<x<t+x_{0}\>\&\>x_{0}-x\leq
t\leq x_{0}+x.\end{cases}$ (39)

The full solution and the uncollided piecewise solution for an initial width
$x_{0}=0.5$ at $t=1$ is plotted in Figure 13. The uncollided flux is
discontinuous in the first derivative at $x=\pm x_{0}$. These discontinuities
travel towards the origin at the wavespeed, meet at $x=0$, then travel
outwards again. Tracking this discontinuity involves moving edges towards the
origin, which preformed poorly in test implementation, so a hybrid static-
moving mesh was adopted for the moving mesh implementation. This method
requires the number of cells to multiples of two greater than or equal to
four. Half of the zones span the source and never move, and half begin at the
edges with width zero and move outwards. Therefore, edges on the left edges of
the source obey

$x_{k}(t)=x_{ko}-\frac{x_{k}}{x_{0}}\times
vt\qquad\mathrm{for}\>k<\frac{K}{4}.$ (40)

Mesh edges from $k=\frac{K}{4}$ to $k=\frac{3K}{4}$ stay at their initial
values. For the rest,

$x_{k}(t)=x_{ko}+\frac{x_{k}}{x_{K}}\times
vt\qquad\mathrm{for}\>k>\frac{3K}{4}.$ (41)

While this mesh method does not track the interior discontinuity, it ensures
better resolution of the solution inside the source where the uncollided
solution presents difficulty by allocating two-thirds of the cells to the
source region. The static mesh in this case spans
$[-t_{\mathrm{final}}-x_{0},t_{\mathrm{final}}+x_{0}]$ with an edge at the
initial source width.

For our uncollided source treatment, Eq. (39) is substituted into Eq. (14) and
the initial condition is set to zero. For the standard source treatment, the
initial condition is found by substituting
$\psi(x,t=0)=\frac{1}{2}\Theta\left(x_{0}-|x|\right)$ into Eq. (16) and
setting $S$ in Eq. (14) to zero.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 13: Square pulse semi-analytic solution, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed) with $x_{0}=0.5$ and $c=1$. The uncollided flux
is not shown for times where it is negligible.

The square pulse source is significantly smoother than the plane pulse. The
former has an uncollided angular flux that is piecewise smooth, where the
latter has an uncollided angular flux made of of travelling delta functions.
Figure 14 shows that the uncollided methods are able to achieve second order
convergence at early time, with the standard source treatments doing slightly
worse. The square pulse requires similar angular error resolution as do the
smooth Gaussian problems, with $S_{512}$ required for RMSE less than
$10^{-6}$.

As with the plane pulse, the intercept value is important. For six basis
functions, the uncollided, moving mesh and the uncollided static mesh have the
same order of convergence, but the former has a 10 times intercept reduction
from the baseline and the latter is reduced 5 times (Table 3). Table 3 also
shows a feature that is less discernible in Figure 14, the intercept for $M=6$
is approximately 3 times smaller than for $M=4$.

(a) $M=4$

(b) $M=6$

Figure 14: Square pulse convergence results on a logarithmic scale with $c=1$ at $t=1$. Blue lines indicate the uncollided solution is used, red that no uncollided source is used. Dashed lines are for a static mesh and solid lines are for the moving mesh.  Table 3: Square pulse intercepts from Figure 11 where the intercept, $C$ is found from the curve fit, $\mathrm{RMSE}=C\,\Delta x^{-A},$ where $\Delta x$ is the number of mesh subdivisions and $A$ is the order of convergence. The improvement from the baseline is found by dividing the intercept from the no uncollided solutions, static mesh case by the intercept from the given case | M = 4 | M = 6
---|---|---
| intercept | improvement | intercept | improvement
no uncollided $+$ static mesh | 0.14594 | - | 0.03220 | -
uncollided $+$ static mesh | 0.02664 | 6 | 0.006501 | 5
no ucollided $+$ moving mesh | 0.0078136 | 19 | 0.00414 | 8
uncollided $+$ moving mesh | 0.006379 | 23 | 0.002266 | 10

###  6.6 Square source

We define a square source which is a superposition of square pulses while
$t<t_{0}$,

$S_{ss}=\Theta(x_{0}-|x|)\Theta(t_{0}-t).$ (42)

The solution for this configuration is shown in Figure (15) where the
uncollided solution is continuous, but not everywhere smooth and is still
significant by the time the source is turned off, $t=5$. The uncollided
solution for this source, given by [7], is

$\phi_{u}^{\mathrm{sp}}(x,t)=\left[-x_{0}\mathrm{E_{i}}(\tau-t)\right]\bigg{|}_{0}^{b}+\frac{1}{2}\left[(|x|-x_{0})\text{Ei}(\tau-t)+e^{\tau-t}\right]\bigg{|}_{b}^{c}+\left[e^{-(t-\tau)}\right]\bigg{|}_{c}^{d}.$
(43)

Where the intervals are defined by

$b=\left[\mathrm{min}\left(d,t-|x|-x_{0}\right)\right]_{+},$ (44)
$c=\left[\mathrm{min}\left(d,t+|x|-x_{0}\right)\right]_{+},$ (45)
$d=\left[\mathrm{min}\left(t_{0},t,t-|x|+x_{0}\right)\right]_{+}.$ (46)

$[\cdot]_{+}$ returns the positive part of its argument and $\mathrm{E_{i}}$
is the exponential integral.

To find the source term in the weak formulation, Eq. (43) is inserted into Eq.
(14). For the standard source evaluation, Eq. (42) is used as the source in
Eq. (14).

The square source uncollided flux creates a discontinuity in the first
derivative at $x=\pm x_{0}$ while the source is on. The moving mesh defined in
Section 6.5 is ideal for handling this discontinuity since a mesh edge is
always at $\pm x_{0}$. Also, a two-thirds of the cells available cells are
inside of the source where the solution is more difficult to resolve. The
static mesh in this case spans
$[-t_{\mathrm{final}}-x_{0},t_{\mathrm{final}}+x_{0}]$ with evenly spaced
cells. An edge is always located at the source edge.

(a) $t=1$

(b) $t=5$

(c) $t=10$

Figure 15: Square source semi-analytic solution, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed) with $x_{0}=0.5$, $t_{0}=5$, and $c=1$. The
uncollided flux is not shown for times where it is negligible.

The uncollided source methods for the square source achieves second order
convergence at early times and the standard source treatments do slightly
worse, shown in Figure 16. Like the square pulse, discontinuities in the
uncollided angular flux restrict the order of convergence of the uncollided
source methods. The angular error dependence in this case is similar to the
square pulse, with $S_{512}$ required for errors less than $10^{-6}$.

Since the uncollided solution in this case lasts longer than the square pulse
case before decaying to zero, using the uncollided source, moving mesh
treatment creates a more significant intercept reduction from the standard
method. Table 4 shows that for the $M=6$ case, the intercept is reduced 70
times from the standard method. It is interesting to note that, according to
Eq. (25), for $M=6$, the standard DG method would require 1300 hundred spatial
subdivisions to achieve an error that the uncollided, moving mesh method
returns with 32 cells.

The $t=1$ convergence results for every case except the MMS problem for the
uncollided source, moving mesh and the standard source, moving mesh are shown
side by side in Figure 21. This figure is a good illustration of the
difference in convergence between sources within the same method and the
difference in accuracy between the uncollided, moving mesh and the standard
moving mesh methods.

Table 4: Square source intercepts from Figure (16) where the intercept, $C$ is found from the curve fit $\mathrm{RMSE}=C\,\Delta x^{-A}$, where $\Delta x$ is the number of mesh subdivisions and $A$ is the order of convergence. The improvement from the baseline is found by dividing the intercept from the no uncollided solutions, static mesh case by the intercept from the given case | M = 4 | M = 6
---|---|---
| intercept | improvement | intercept | improvement
no uncollided $+$ static mesh | 0.1713 | - | 0.04705 | -
uncollided $+$ static mesh | 0.03235 | 5 | 0.009397 | 5
no ucollided $+$ moving mesh | 0.05331 | 3 | 0.003389 | 14
uncollided $+$ moving mesh | 0.01674 | 10 | 0.0006706 | 70

(a) $M=4$

(b) $M=6$

Figure 16: Square source convergence results at $t=1$. Blue lines indicate the
uncollided solution is used, red that no uncollided source is used. Dashed
lines are for a static mesh and solid lines are for the moving mesh.

###  6.7 Results for non-purely scattering problems

All of the test problems considered so far have been pure scattering problems
($c=1$). In this section, we briefly show that for the Gaussian pulse and
square pulse problems, our method performs well with partially absorbing and
multiplying scattering ratios. For these problems, we are able to avoid
recalculating a benchmark by rescaling already calculated results. In an
underappreciated footnote, Case and Zweifel ([30, p. 175]) presented a clever
scaling to relate the solution for $c=1$ ($\psi^{1}$) to solutions for any
other nonzero scattering ratio ($\psi^{c}$),

$\psi^{c}(x,\mu,t)=c\exp\left(-(1-c)t\right)\psi^{1}(cx,\mu,ct).$ (47)

This scaling was derived for initial value problems with no source term, hence
the choice of the pulse-type problems in this section.

It is also necessary to scale the initial condition of each source in the code
implementation so that our code is running the same problem as the benchmark.
This is done by scaling the parameters, $x_{0}$ and $\sigma$ for the square
pulse and Gaussian pulse respectively, in the initial condition and the
uncollided source,

$x_{0}^{\prime}=\frac{x_{0}}{c}\qquad\sigma^{\prime}=\frac{\sigma}{c},$ (48)

where $x_{0}^{\prime}$ and $\sigma^{\prime}$ are the new parameters for our
source term. Also, the scaled evaluation time becomes,

$t_{\mathrm{final}}=\frac{t}{c},$ (49)

where $t$ is the evaluation time of the benchmark problem.

For $c=0.8$ and $c=1.2$, the benchmarks for the square pulse and Gaussian
pulse are plotted in Figures 17 and 18 respectively. We choose the evaluation
times $t_{\mathrm{final}}=1.25$ and $t_{\mathrm{final}}=0.83$ so we may use
the $t=1$ benchmark.

Convergence tests showed that our moving mesh, uncollided solutions method
handles these partially absorbing or multiplying problems just as well as
purely scattering problems. The results for the plane pulse behave similarly
to $c=1$ results, with second order convergence at early times and a
significant difference in magnitude between the uncollided, moving mesh case
and the standard DG implementation. These results are shown in Figure 20.
Changing the scattering ratio also did not impact the convergence behavior of
the method on the Gaussian pulse problem, with the method achieving fifth
order convergence with $M=4$ for the tested values of $c$ (Figure 19). In both
tests, the multiplying scattering ratio showed similar convergence
characteristics to the absorbing problem, but with a higher relative error.
This is not necessarily an effect of the larger scattering ratio, but more
likely of earlier evaluation time, when the solution is more challenging.

(a) $c=0.8$, $t=1.25$

(b) $c=1.2$, $t=0.83$

Figure 17: Square pulse semi-analytic solutions, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed), at two times for different values of $c$, scaled
from $t=1$, $c=1$ solution (13(a)).

(a) $c=0.8$, $t=1.25$

(b) $c=1.2$, $t=0.83$

Figure 18: Gaussian pulse semi-analytic solutions, $\phi$ (solid) and
$\phi_{\mathrm{u}}$ (dashed), at two times for different values of $c$, scaled
from $t=1$, $c=1$ solution (4(a))

(a) $M=4$, $c=0.8$, $t=1.25$

(b) $M=4$, $c=1.2$, $t=0.83$

Figure 19: Gaussian pulse convergence results on a logarithmic scale with
$c=0.8$ evaluated at $t=1.25$ and $c=1.2$ evaluated at $t=0.83$. The standard
deviation of the initial condition, $\sigma$ is $0.625$ in the first case and
$0.417$ in the second. Blue lines indicate the uncollided solution is used,
red that no uncollided source is used. Dashed lines are for a static mesh and
solid lines are for the moving mesh.

(a) $M=4$, $c=0.8$, $t=1.25$

(b) $M=4$, $c=1.2$, $t=0.83$

Figure 20: Square pulse convergence results on a logarithmic scale with
$c=0.8$ evaluated at $t=1.25$ and $c=1.2$ evaluated at $t=0.83$. The standard
deviation of the initial condition, $\sigma$ is $0.625$ in the first case and
$0.417$ in the second. Blue lines indicate the uncollided solution is used,
red that no uncollided source is used. Dashed lines are for a static mesh and
solid lines are for the moving mesh.

(a) Moving mesh

(b) Static mesh

Figure 21: Convergence results for an increasing number of cell divisions,
$K$, for $t=1$ for $M=6$ for the uncollided source moving mesh (21(a)) and
static mesh (21(b)) cases for every test problem except the MMS problem.

###  6.8 Computational efficiency

For the square source problem of Section 6.6, benchmarks were created storing
the average computational time over $5$ runs holding $M=6$ constant and
increasing number of cell divisions and holding the number of cells constant
($4$) while increasing the number of basis functions, Figure 22. The
calculations were performed on a 2020 MacBook Pro Laptop with an 8 core M1 CPU
and 16GB of RAM. For both cases, the uncollided source, moving mesh method
returned more accurate solutions for less calculation time. The results also
show that the most efficient way to obtain accurate solutions with this method
is to use less mesh subdivisions and a higher order polynomial. The most
accurate solution in Figure 22(a) takes about three times less time than the
most accurate solution in Figure (22(b)). This can also be understood in light
of the trend apparent in Figures 11, 14, and 16 where increasing the number of
basis functions from 5 to 7 returns the same convergence but starting from a
smaller intercept value.

(a) $4$ cells, $M=2,4,8$ and $16$

(b) $M=6$, $2,4,8$ and $16$ cells

Figure 22: Calculation time vs RMSE for five runs of the square source
problem. Blue lines indicate the uncollided solution is used, dashed that the
mesh is static, and solid that the mesh is moving. Panel (22(a)) contains data
from calculations holding the number of mesh subdivisions constant and
increasing the number of basis functions and Panel (22(b)) contains data from
calculations holding $M$ constant and increasing the number of mesh
subdivisions.

##  7 Conclusions

For smooth problems, like the MMS problem, the Gaussian pulse, the Gaussian
source, and the nonsmooth sources at later times, the method of uncollided
solutions with a moving mesh achieves spectral convergence and produces highly
accurate solutions for relatively few degrees of freedom. In the Gaussian
source case, there is a definite benefit to using the uncollided solution, but
it does not seem that using a moving mesh gives a large added benefit since
the solution does not contain discontinuous wavefronts. The moving mesh,
uncollided source treatment was most useful in problems with finite width,
nonsmooth sources. For the plane pulse problem, the highly nonsmooth angular
flux prohibited high order convergence, but also amplified the difference in
efficiency between the uncollided, moving mesh case and the standard DG
implementation. The comparatively smooth square pulse and square source
achieved second order convergence at early times, also with the uncollided and
moving mesh performing the best by a significant margin. These cases also
showed an improvement of the order of convergence from using the uncollided
solutions. For the square source case, we also achieved second order
convergence with significant difference in solution accuracy achieved with a
moving mesh and the uncollided source. Also for this case, we showed that
using a higher order polynomial to interpolate the solution and a lower number
of mesh subdivisions returned accurate solutions with low computational
expense compared to a method with a lower order polynomial and more mesh
subdivisions. Additionally, we showed for that the method preforms well for
scattering ratios not equal to one ($c\neq 1$) with tests of the square and
Gaussian pulse source problems.

Given the encouraging results from this study, we plan for future work to
utilize the moving mesh and uncollided flux treatments to solve problems of
nonlinear radiative transfer to create novel verification-quality benchmark
solutions. Since there are no inherent challenges to applying the
Discontinuous Galerkin method (and spectral methods in general) to nonlinear
problems, we do not anticipate any special difficulties other than those
introduced by the allowance of negative solutions in this undertaking. Systems
with anisotropic scattering, finite medium problems, and problems in spherical
and cylindrical geometry are also of interest. We also suggest that these
approaches could be used for other, transport problems not amenable to
(semi-)analytic solution technique such as neutronics problems with thermal
feedback or rarefied gas dynamics.

## Acknowledgement

This work was supported by a Los Alamos National Laboratory under contract
#599008, “Verification for Radiation Hydrodynamics Simulations”.

## References

  * [1] B. D. Ganapol, R. S. Baker, J. A. Dahl, R. E. Alcouffe, Homogeneous infinite media time-dependent analytical benchmarks, Tech. Rep. LA-UR-01-1854, Los Alamos National Laboratory (2001).
URL https://www.osti.gov/biblio/975281

  * [2] D. Amsallem, M. J. Zahr, C. Farhat, Nonlinear model order reduction based on local reduced-order bases, International Journal for Numerical Methods in Engineering 92 (10) (2012) 891–916.
  * [3] A. K. Henrick, T. D. Aslam, J. M. Powers, Simulations of pulsating one-dimensional detonations with true fifth order accuracy, Journal of Computational Physics 213 (1) (2006) 311–329.
  * [4] R. E. Alcouffe, R. D. O’Dell, F. W. Brinkley, A first-collision source method that satisfies Discrete Sn Transport balance, Nuclear Science and Engineering 105 (2) (1990) 198–203.  doi:10.13182/nse90-a23749.
  * [5] C. D. Hauck, R. G. McClarren, A collision-based hybrid method for time-dependent, linear, kinetic transport equations, Multiscale Modeling & Simulation 11 (4) (2013) 1197–1227.  doi:10.1137/110846610.
  * [6] W. J. Walters, A. Haghighat, The adaptive collision source method for discrete ordinates radiation transport, Annals of Nuclear Energy 105 (2017) 45–58.  doi:10.1016/j.anucene.2017.02.013.
  * [7] W. Bennett, R. G. McClarren, Benchmarks for infinite medium, time dependent transport problems with isotropic scattering (2022).  doi:10.48550/ARXIV.2205.15783.
URL https://arxiv.org/abs/2205.15783

  * [8] B. Su, G. L. Olson, An analytical benchmark for non-equilibrium radiative transfer in an isotropically scattering medium, Annals of Nuclear Energy 24 (13) (1997) 1035–1055.  doi:https://doi.org/10.1016/S0306-4549(96)00100-4.
URL https://www.sciencedirect.com/science/article/pii/S0306454996001004

  * [9] S. Bingjing, G. L. Olson, Benchmark results for the non-equilibrium marshak diffusion problem, Journal of Quantitative Spectroscopy and Radiative Transfer 56 (3) (1996) 337–351.  doi:https://doi.org/10.1016/0022-4073(96)84524-9.
URL https://www.sciencedirect.com/science/article/pii/0022407396845249

  * [10] R. McClarren, J. Holloway, T. Brunner, Analytic $\mathrm{P}_{1}$ solutions for time-dependent, thermal radiative transfer in several geometries, Journal of Quantitative Spectroscopy & Radiative Transfer - J QUANT SPECTROSC RADIAT 109 (2008) 389–403.  doi:10.1016/j.jqsrt.2007.08.006.
  * [11] R. G. McClarren, J. G. Wöhlbier, Solutions for ion–electron–radiation coupling with radiation and electron diffusion, Journal of Quantitative Spectroscopy and Radiative Transfer 112 (1) (2011) 119–130.  doi:https://doi.org/10.1016/j.jqsrt.2010.08.015.
URL https://www.sciencedirect.com/science/article/pii/S002240731000333X

  * [12] G. Pomraning, An extension of the eddington approximation, Journal of Quantitative Spectroscopy and Radiative Transfer 9 (3) (1969) 407–422.
  * [13] B. Ganapol, G. Pomraning, The non-equilibrium marshak wave problem: A transport theory solution, Journal of Quantitative Spectroscopy and Radiative Transfer 29 (4) (1983) 311–320.  doi:https://doi.org/10.1016/0022-4073(83)90049-3.
URL https://www.sciencedirect.com/science/article/pii/0022407383900493

  * [14] A. G. Petschek, R. E. Williamson, J. K. Wooten Jr., The penetration of radiation with constant driving temperature, Tech. rep., Los Alamos Scientific Lab. (1960).
  * [15] R. E. Marshak, Effect of radiation on shock wave behavior, The Physics of Fluids 1 (1) (1958).
  * [16] W. Bennett, R. 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) (2021) 20210119.
  * [17] T. K. Lane, R. G. McClarren, New self-similar radiation-hydrodynamics solutions in the high-energy density, equilibrium diffusion limit, New Journal of Physics 15 (9) (2013) 095013.  doi:10.1088/1367-2630/15/9/095013.
URL https://doi.org/10.1088/1367-2630/15/9/095013

  * [18] B. G. Verma, J. P. Vishwakarma, An exact similarity solution for a spherical shock wave in a magneto-radiative gas, Astrophysics and Space Science 58 (1978) 139–147.
  * [19] J. H. Hammer, M. D. Rosen, A consistent approach to solving the radiation diffusion equation, Physics of Plasmas 10 (1829) (2003).  doi:10.1063/1.1564599.
  * [20] S. I. Heizler, Asymptotic telegrapher’s equation (p1) approximation for the transport equation, Nuclear Science and Engineering 166 (1) (2010) 17–35.  doi:10.13182/nse09-77.
  * [21] J. Marsden, A. Tromba, Vector Calculus, no. pp. 1-93 in Vector Calculus, W. H. Freeman, 2003.
URL https://books.google.com/books?id=LiRLJf2m_dwC

  * [22] S. Chandrasekhar, Radiative transfer, Dover Publications, 1960.
  * [23] N. Schlömer, N. Papior, D. Arnold, J. Blechta, R. Zetter, nschloe/quadpy: None (2021).  doi:10.5281/ZENODO.5541216.
URL https://zenodo.org/record/5541216

  * [24] S. K. Lam, A. Pitrou, S. Seibert, Numba, in: Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC - LLVM ’15, ACM Press, 2015, pp. 1–6.  doi:10.1145/2833157.2833162.
  * [25] E. Hairer, S. P. Norsett, G. Wanner, Solving Ordinary Differential Equations I, Springer Berlin Heidelberg, 1993.  doi:10.1007/978-3-540-78862-1.
  * [26] P. Virtanen, R. Gommers, et al., SciPy 1.0: fundamental algorithms for scientific computing in Python, Nature Methods 17 (3) (2020) 261–272.  doi:10.1038/s41592-019-0686-2.
  * [27] C. Lingus, Analytical test case’s [sic] for neutron and radiation transport codes, in: Second conference on transport theory, Los Alamos, New Mexico, 1971, p. 655.
  * [28] P. Knupp, K. Salari, Verification of computer codes in computational science and engineering, Chapman and Hall/CRC, 2002.
  * [29] R. G. McClarren, R. B. Lowrie, Manufactured solutions for the p1 radiation-hydrodynamics equations, Journal of Quantitative Spectroscopy and Radiative Transfer 109 (15) (2008) 2590–2602.
  * [30] K. Case, P. Zweifel, Linear Transport Theory, Addison-Wesley series in nuclear engineering, Addison-Wesley Publishing Company, 1967.
URL https://books.google.com/books?id=uQtRAAAAMAAJ