# AIR multigrid with GMRES polynomials (AIRG) and additive preconditioners for
Boltzmann transport111UK Ministry of Defence © Crown owned copyright 2023/AWE

S. Dargaville<EMAIL_ADDRESS>R.P. Smedley-Stevenson  P.N. Smith
C.C. Pain  Applied Modelling and Computation Group, Imperial College London,
SW7 2AZ, UK  AWE, Aldermaston, Reading, RG7 4PR, UK  ANSWERS Software Service,
Jacobs, Kimmeridge House, Dorset Green Technology Park, Dorchester, DT2 8ZB,
UK

###### Abstract

We develop a reduction multigrid based on approximate ideal restriction (AIR)
for use with asymmetric linear systems. We use fixed-order GMRES polynomials
to approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$ and we use these polynomials
to build grid transfer operators and perform F-point smoothing. We can also
apply a fixed sparsity to these polynomials to prevent fill-in.

When applied in the streaming limit of the Boltzmann Transport Equation (BTE),
with a P0 angular discretisation and a low-memory spatial discretisation, this
“AIRG” multigrid used as a preconditioner to an outer GMRES iteration
outperforms the lAIR implementation in hypre, with two to three times less
work. AIRG is very close to scalable; we find either fixed work in the solve
with slight growth in the setup, or slight growth in the solve with fixed work
in the setup when using fixed sparsity. Using fixed sparsity we see less than
20% growth in the work of the solve with either 6 levels of spatial refinement
or 3 levels of angular refinement. In problems with scattering AIRG performs
as well as lAIR, but using the full matrix with scattering is not scalable.

We then present an iterative method designed for use with scattering which
uses the additive combination of two fixed-sparsity preconditioners applied to
the angular flux; a single AIRG V-cycle on the streaming/removal operator and
a DSA method with a CG FEM. We find with space or angle refinement our
iterative method is very close to scalable with fixed memory use.

###### keywords:

Asymmetric multigrid , Advection , Radiation transport , Boltzmann , AIR ,
GMRES polynomials

##  1 Introduction

The Boltzmann transport equation (BTE) describes the distribution of particles
moving through an interacting medium and is used to model radiation transport,
along with spectral-waves and fluid problems through kinetic and lattice-
Boltzmann methods. The mono-energetic steady-state form of the Boltzmann
Transport Equation (BTE), with linear scattering and straight-line propagation
is written in (1) as

$\bm{\Omega}\cdot\nabla_{\bm{r}}\bm{\psi}(\bm{r},\bm{\Omega})+\sigma_{\textrm{t}}\bm{\psi}(\bm{r},\bm{\Omega})-\int_{\bm{\Omega}^{\prime}}\sigma_{\textrm{s}}(\bm{r},\bm{\Omega}^{\prime}\rightarrow\bm{\Omega})\bm{\psi}(\bm{r},\bm{\Omega}^{\prime})\textrm{d}\bm{\Omega}^{\prime}=S_{\textrm{e}}(\bm{r},\bm{\Omega}).$
(1)

Equation (1) is a 5-dimensional linear PDE, with three spatial dimensions and
two angular dimensions; we neglect the energy and time dimensions. The angular
flux, $\bm{\psi}(\bm{r},\bm{\Omega})$, describes the number of particles
moving in direction $\bm{\Omega}$, at spatial position $\bm{r}$. The
macroscopic total cross section of the material that the particles are moving
through is given by $\sigma_{\textrm{t}}$, which describes particles removed
either through absorption or scattering in the material. The source of
particles coming from scattering in many radiative processes is given by an
integral term, where $\sigma_{\textrm{s}}$ is the macroscopic scatter cross-
sections for this process that describes how particles scatter from direction
$\bm{\Omega}^{\prime}$ into direction $\bm{\Omega}$. Finally any external
sources of particles are given by $S_{\textrm{e}}$.

One of the main challenges in solving (1) is that when the scattering cross-
sections are large (i.e., the particles are interacting with the material),
the BTE tends to a diffusion equation, whereas when the scattering and total
cross-sections are zero (i.e., particles are moving through a vacuum), the BTE
is purely hyperbolic and hence stable discretisations must be used and the
resulting linear systems are asymmetric and non-normal.

In (1), the scattering cross-section, $\sigma_{\textrm{s}}$, for any given
angle-to-angle scattering event is often described by a Legendre expansion
whose coefficients we denote as $\bm{\sigma}_{\textrm{s}}$. If we discretise
in space/angle (with a discretisation like Sn, or FEM), we introduce,
$\bm{\varphi}$, which is the angular flux, $\bm{\psi}$, in Legendre space. We
then write (1) as a $2\times 2$ block system, namely

$\begin{bmatrix}\textrm{{I}}&\textrm{{D}}_{\textrm{m}}\\\
-\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}&\textrm{{L}}\end{bmatrix}\begin{bmatrix}\bm{\varphi}\\\
\bm{\psi}\end{bmatrix}=\begin{bmatrix}\bm{0}\\\
\hat{\textrm{{S}}}_{\textrm{e}}\end{bmatrix},$ (2)

where L is the streaming/removal operator, $\bm{\Sigma}_{\textrm{s}}$ is a
matrix with the scattering cross-sections for each spatial node,
$\textrm{{D}}_{\textrm{m}}$ and $\textrm{{M}}_{\textrm{m}}$ are formed from
the tensor product of the spatial mass matrices and the mapping operators
which map between our angular discretisation and the moments of the Legendre
space and $\hat{\textrm{{S}}_{\textrm{e}}}$ is the discretised source term
from (1). When discretised in space with an upwind discretisation and an
appropriate ordering of $\bm{\psi}$ is applied, L is a block-diagonal matrix,
where each of the blocks is lower triangular and corresponds to the spatial
coupling for each direction (i.e., each direction in L is not coupled to the
others; L represents a set of advection equations for each direction).

Equation (2) is typically solved by forming the Schur complement of block L,
namely

$(\textrm{{I}}+\textrm{{D}}_{\textrm{m}}\textrm{{L}}^{-1}\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}})\bm{\varphi}=-\textrm{{D}}_{\textrm{m}}\textrm{{L}}^{-1}\hat{\textrm{{S}}}_{\textrm{e}},$
(3)

and then a preconditioned Richardson iteration (known as a source iteration in
the transport community), with preconditioner $\textrm{{M}}^{-1}$ is applied
to recover

$\bm{\varphi}^{n+1}=\bm{\varphi}^{n}-\textrm{{M}}^{-1}(\textrm{{D}}_{\textrm{m}}\textrm{{L}}^{-1}(\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\bm{\varphi}^{n}+\hat{\textrm{{S}}}_{\textrm{e}})+\bm{\varphi}^{n}).$
(4)

Computing the solution to (2) therefore only requires the Legendre
representation of the angular flux (the angular flux can easily be formed,
either at a single spatial point or across the domain if needed), as each of
the components of (4) are block-diagonal. This allows for a very low memory
iterative method. If there is scattering in the problem, then typically an
additive preconditioner like
$\textrm{{M}}^{-1}=\textrm{{I}}+\textrm{{D}}_{\textrm{diff}}^{-1}$ is used;
this is known as diffusion-synthetic acceleration (DSA), where
$\textrm{{D}}_{\textrm{diff}}$ is a diffusion operator, with
diffusion/sink/boundary coefficients taken from an asymptotic analysis of the
BTE in the diffusion limit.

The spatial discretisation applied to the diffusion operator in DSA can govern
its effectiveness in accelerating convergence; research into discretisations
of the diffusion operator that are effective and/or “consistent” with
transport are extensive. The use of a Krylov method to solve (3) instead of a
Richardson method can reduce the dependence on this consistency [1] and hence
“inconsistent” discretisations can be applied to the diffusion operator as
part of a DSA resulting in SPD systems that can be solved efficiently.
Transport synthetic acceleration (TSA) [2] has also been used to accelerate
convergence in the scattering limit, where lower order transport solutions are
used instead of diffusion; see [3] for a review of both DSA and TSA method in
transport.

The solution of (3) relies on the exact inversion of L; if as mentioned above
the spatial and angular discretisations result in L with lower-triangular
structure then L can be inverted exactly with a single (matrix-free) Gauss-
Seidel (GS) iteration, also known as a sweep in the transport community. If
unstructured spatial grids are used, inverting L can become more difficult,
particularly in parallel; finding an ideal sweep ordering is then NP-complete.
Furthermore, the introduction of different spatial or angular discretisations,
time dependence or additional physics that cannot be mapped to Legendre space
compounds this problem. Similarly, if we wish to use angular adaptivity where
the angular resolution differs throughout the spatial grid, L is no longer
block-diagonal.

Our goal in the AMCG has been to investigate different discretisations,
adaptive and iterative methods for solving the BTE, that have the potential to
overcome these difficulties while remaining scalable (i.e., perform a fixed
amount of work with a fixed memory consumption, even with space/angle
refinement). We instead form the Schur complement of block I from (2) and
solve the system formed with the angular flux, namely

$(\textrm{{L}}+\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\textrm{{D}}_{\textrm{m}})\bm{\psi}=\hat{\textrm{{S}}}_{\textrm{e}}.$
(5)

There are several disadvantages to solving (5) instead of (3), the most
important of which is the considerable increase in memory required. The
angular flux, $\bm{\psi}$, is much bigger than $\bm{\varphi}$ and the matrix
$\textrm{{L}}+\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\textrm{{D}}_{\textrm{m}}$
is not block diagonal; instead the scattering operator
$\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\textrm{{D}}_{\textrm{m}}$
couples different angles together and results in dense angle-angle blocks
where the nnzs scale like $\mathcal{O}(n^{2})$ with angle size. These
disadvantages would typically preclude the development of a practical
transport algorithm. Previously we have tackled those problems through the
combination of: using a stable spatial discretisation based on static
condensation that has the stencil of a CG discretisation (reducing the size of
$\bm{\psi}$ but at the cost of making the blocks in L no longer lower
triangular), using angular adaptivity to only focus angular resolution where
required (reducing the size of $\bm{\psi}$), and using a matrix-free multigrid
to solve (5) that does not rely on the explicit construction of
$\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\textrm{{D}}_{\textrm{m}}$
or on the lower triangular structure in L.

We showed previously [4, 5] that these techniques perform well on many
transport problems, allowing the practical use of both angular adaptivity with
high levels of refinement and unstructured spatial grids. These methods do not
use GS/sweep smoothers however and do not scale well in the streaming limit
where $\bm{\sigma}_{\textrm{s}}$ tends to zero. Similarly, most multigrid
methods in the literature that achieve good performance for the BTE have
relied on either block-based smoothers, often on a cell/element, which do not
scale with increasing angle size but which perform well in the scattering
limit [6, 7, 8, 9], or GS/sweeps as smoothers which perform well in the
streaming limit [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25] but as mentioned can be difficult to parallelise on unstructured grids.
Developing multigrid methods which perform well with local smoothers is
difficult, particularly as hyperbolic problems have long proved difficult to
solve with multigrid due to the lack of developed theory for non-SPD matrices.
Recently multigrid methods based on approximate ideal restrictors (AIR) have
been developed [26, 27, 28, 29] that show good convergence in asymmetric
problems, and in particular with the BTE when using point-based smoothers
[29].

The aim of this work is hence to build an iterative method that solves (5)
with the following features:

  1. 1.

Compatible with the space/angle discretisations we have developed previously
[4, 5] and hence does not rely on L having block diagonal and/or lower
triangular structure

  2. 2.

Does not require the explicit construction of
$\textrm{{M}}_{\textrm{m}}\bm{\Sigma}_{\textrm{s}}\textrm{{D}}_{\textrm{m}}$

  3. 3.

Does not require the use of GS/sweeps

  4. 4.

Scalable in both the streaming and scattering limits with space/angle
refinement

  5. 5.

Compatible with angular adaptivity

  6. 6.

Good strong and weak scaling in parallel with unstructured grids

We examine the first four of these points in this paper, and leave
investigation of adaptivity and the parallel performance of our method to
future work. Here we present two contributions: the first is an algaebraic
multigrid method based on combining AIR with low-order GMRES polynomials,
which we call AIRG; the second is an iterative method where we apply additive
preconditioners to an outer GMRES iteration on the angular flux, based on the
streaming/removal operator and a DSA diffusion operator. Both these
contributions can be used independently of the other; for example we could use
AIRG to invert L and/or the DSA operator as part of a typical DG FEM/source
iteration.

We then compare our iterative method to the hypre implementation of lAIR and
find performance advantages. We therefore have an iterative method on
unstructured grids which never requires the assembly of the full matrix in
(5), that has both good performance and fixed memory consumption across all
parameter regimes with space/angle refinement for Boltzmann transport
problems.

##  2 Discretisations

We begin by presenting the spatial and angular discretisations used in this
work; they are based on those presented in [30, 31, 32, 33, 4, 5] and hence we
only discuss their key features.

###  2.1 Angular discretisation

We use a P0 DG FEM in angle (or equivalently a cell-centred FVM) with constant
area azi/polar elements that we normalise so that the angular mass matrix is
the identity. The first level of our angular discretisation is denoted level
1, with one constant basis function per octant. Each subsequent level of
refinement comes from splitting an angular element into four at the midpoint
of the azi/cosine polar bounds; this is structured, nested refinement. This is
equivalent to the Haar wavelet discretisations discussed in [4]; indeed as
part of the matrix-free iterative method in [4] an $\mathcal{O}(n)$ mapping
to/from this P0 space to Haar space was performed during every matvec. We
instead choose to solve in P0 space as we can form an assembled copy of the
streaming/removal matrix that has fixed sparsity with angular refinement; this
property was not needed for the matrix-free methods used in our previous work,
all we required was a scalable mapping between P0 and Haar spaces. We can also
adapt in this P0 space in the same manner as our wavelet space [4], given the
equivalence. We investigate adapting in P0 space in future work.

###  2.2 Spatial discretisation

Our spatial discretisation is a sub-grid scale FEM, which represents the
angular flux as $\bm{\psi}=\bm{\phi}+\bm{\theta}$, where $\bm{\phi}$ is the
solution on a “coarse” scale and $\bm{\theta}$ is the solution on a “fine”
scale. The finite element expansions for both the fine and coarse scales can
be written as

$\phi(\bm{r},\bm{\Omega})\approx\sum_{i=1}^{\eta_{N}}N_{i}(\bm{r})\tilde{\phi}_{i}(\bm{\Omega});\qquad\theta(\bm{r},\bm{\Omega})\approx\sum_{i=1}^{\eta_{Q}}Q_{i}(\bm{r})\tilde{\theta}_{i}(\bm{\Omega}),$
(6)

with $\eta_{N}$ continuous basis functions, $N_{i}$, and $\eta_{Q}$
discontinuous basis functions, $Q_{i}$, with $\tilde{\phi}_{i}$ and
$\tilde{\theta}_{i}$ the expansion coefficients, respectively. In this work we
use linear basis functions for both the continuous and discontinuous spatial
expansions.

As described in Section 2.1, we use a P0 discretisation, with (constant) basis
functions $G_{j}(\bm{\Omega})$, with a spatially varying number of angular
elements $\eta_{A}^{i}$ and $\eta_{D}^{i}$ on the coarse and fine scales
respectively (we enforce that DG nodes have the same angular expansion as
their CG counterparts). The expansion coefficients $\tilde{\phi}_{i}$ and
$\tilde{\theta}_{i}$ in (6) can then be written as space/angle expansion
coefficients $\tilde{\phi}_{i,j}$ and $\tilde{\theta}_{i,j}$ via the expansion

$\tilde{\phi}_{i}(\bm{\Omega})\approx\sum_{j=1}^{\eta_{A}^{i}}G_{j}(\bm{\Omega})\tilde{\phi}_{i,j};\qquad\tilde{\theta}_{i}(\bm{\Omega})\approx\sum_{j=1}^{\eta_{D}^{i}}G_{j}(\bm{\Omega})\tilde{\theta}_{i,j}.$
(7)

Following standard FEM theory the discretised form of (1) can then be written
as

$\begin{bmatrix}\textrm{{A}}&\textrm{{B}}\\\ \textrm{{C}}&\textrm{{D}}\\\
\end{bmatrix}\begin{bmatrix}\bm{\Phi}\\\ \bm{\Theta}\\\
\end{bmatrix}=\begin{bmatrix}\textrm{{S}}_{\bm{\Phi}}\\\
\textrm{{S}}_{\bm{\Theta}}\\\ \end{bmatrix},$ (8)

where ${\bm{\Phi}}$ and ${\bm{\Theta}}$ are vectors containing the coarse and
fine scale expansion coefficients, we denote the number of unknowns in
${\bm{\Phi}}$ as NCDOFs and in ${\bm{\Theta}}$ as NDDOFs. The discretised
source terms for both scales are $\textrm{{S}}_{\bm{\Phi}}$ and
$\textrm{{S}}_{\bm{\Theta}}$. We should note the matrices A and D are the
standard CG and DG FEM matrices that result from discretising (1).

We can then form a Schur complement of block D and recover

$(\textrm{{A}}-\textrm{{B}}\textrm{{D}}^{-1}\textrm{{C}}){\bm{\Phi}}=\textrm{{S}}_{\bm{\Phi}}-\textrm{{B}}\textrm{{D}}^{-1}\textrm{{S}}_{\bm{\Theta}}.$
(9)

The fine solution $\bm{\Theta}$ can then be computed

$\bm{\Theta}=\textrm{{D}}^{-1}(\textrm{{S}}_{\bm{\Theta}}-\textrm{{C}}\bm{\Phi}),$
(10)

and our discrete solution is the addition of both the coarse and fine
solutions, namely $\bm{\Psi}=\bm{\Phi}+\bm{\Theta}$ (where the coarse solution
$\bm{\Phi}$ has been projected onto the fine space). In order to solve (9) and
(10) efficiently/scalably, we must sparsify D and this sparsification is
dependent on the angular discretisation used, see [33, 34, 35, 36, 37, 38, 4,
39] for examples. In this work we replace $\textrm{{D}}^{-1}$ in (9) and (10)
with $\hat{\textrm{{D}}}^{-1}$, which is the streaming operator with removal
and self-scatter only, and vacuum conditions applied on each DG element (as
this removes the jump terms that couple the DG elements, resulting in element
blocks). We can then invert this matrix element-wise and store the result,
with a constant nnzs with space/angle refinement, as it has the same stencil
as the streaming operator. Now if we consider the streaming/removal (denoted
with a subscript $\Omega$) and scattering contributions (denoted with a
subscript S) in (9) separately, along with our sparsified
$\hat{\textrm{{D}}}^{-1}$ we can rearrange (9) and write

$\left(\textrm{{A}}_{\Omega}-\textrm{{B}}_{\Omega}\hat{\textrm{{D}}}^{-1}\textrm{{C}}_{\Omega}\right){\bm{\Phi}}+\left((\textrm{{A}}_{\textrm{S}}+\textrm{{B}}_{\textrm{S}}(y+\hat{\textrm{{D}}}^{-1}\textrm{{C}}_{\Omega})+\textrm{{B}}_{\Omega}y\right){\bm{\Phi}}=\textrm{{S}}_{\bm{\Phi}}-(\textrm{{B}}_{\Omega}+\textrm{{B}}_{\textrm{S}})\hat{\textrm{{D}}}^{-1}\textrm{{S}}_{\bm{\Theta}}.$
(11)

where $y=\hat{\textrm{{D}}}^{-1}\textrm{{C}}_{\textrm{S}}$ and our fine
component is
$\bm{\Theta}=\hat{\textrm{{D}}}^{-1}(\textrm{{S}}_{\bm{\Theta}}-(\textrm{{C}}_{\Omega}+\textrm{{C}}_{\textrm{S}})\bm{\Phi})$.
We can see that (11) is now written similarly to (5), where the left term is
(very close to) the sub-grid scale streaming/removal operator, with the
remainder being the contribution from scattering. It is not exactly the
streaming/removal operator as $\hat{\textrm{{D}}}^{-1}$ contains self-scatter,
but it has the same stencil as the streaming/removal operator and we call it
such below for simplicity; practically we can modify our stabilisation such
that $\hat{\textrm{{D}}}^{-1}$ does not contain self-scatter without any
substantial differences to either the stability of our system or the
preconditioning described in Section 3. As mentioned in Section 1, we cannot
explicitly form the scattering contribution as it has dense blocks, but given
that $\hat{\textrm{{D}}}^{-1}$ has fixed sparsity and the scatter
contributions from $\textrm{{A}}_{\textrm{S}},\textrm{{B}}_{\textrm{S}}$ and
$\textrm{{C}}_{\textrm{S}}$ can be formed in Legendre space and mapped back,
we can scalably compute a matrix-free matvec with angular refinement in P0
space (with a fixed scatter order).

Practically, there are many ways to implement a matvec for (11) (these involve
further rearrangement of (11)), depending on the amount of memory available.
Some key points include that away from boundaries, the element matrices in
$\textrm{{A}}_{\Omega}$, $\textrm{{B}}_{\Omega}$, $\textrm{{C}}_{\Omega}$ and
$\hat{\textrm{{D}}}^{-1}$ all share the same (block) sparsity; in fact given
matching angular resolution and the same order of basis functions on the
continuous and discontinuous spatial meshes (we enforce both conditions) two
of the element matrices are identical, with $\textrm{{B}}_{\Omega}$ and
$\textrm{{C}}_{\Omega}$ related by the transpose of the spatial table. A
matvec with these components can therefore consist of performing several
matvecs on small (max $3\times 3$ or $4\times 4$) angular blocks which can be
kept in cache (similar to the work performed during a DG sweep, but without
the Gauss-Seidel dependency on ordering). We can also combine a number of the
maps to/from Legendre space. We describe some of the choices we make in
Section 7, but note that a FLOP count shows our sub-grid scale matvec with
scattering can be computed with 1.8x the FLOPs of an equivalent DG matvec.

With our sub-grid scale discretisation, we are therefore choosing to increase
the number of (local) FLOPs in order to significantly decrease our memory
consumption. We can build iterative methods that only depend on the the
coarse-scale solution, ${\bm{\Phi}}$, which is on a CG stencil, which in 3D
has approximately $20\times$ fewer unknowns than a DG method; e.g., a
GMRES(20) space built on $\bm{\Phi}$ would take approximately the same space
as a single copy of $\bm{\Psi}$ in 3D. One additional benefit is as noted in
[4], is that with either a wavelet or non-wavelet angular discretisation, our
sub-grid scale discretisation does not require interpolation between areas of
different angular resolution (e.g., across faces), due to the lack of DG jump
terms.

##  3 Additively preconditioned iterative method

This section details the iterative method we use to solve (11). For simplicity
we begin by writing the unseparated form (9) with our sparsified
$\hat{\textrm{{D}}}^{-1}$ and introduce a preconditioner $\textrm{{M}}^{-1}$
on the right to give

$(\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}})\textrm{{M}}^{-1}\textrm{{u}}=\textrm{{S}}_{\bm{\Phi}}-\textrm{{B}}\textrm{{D}}^{-1}\textrm{{S}}_{\bm{\Theta}},\quad\textrm{{u}}=\textrm{{M}}\bm{\psi}.$
(12)

We solve (12) with GMRES and use a matrix-free matvec as described above to
compute the action of
$(\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}})$ (and to
compute the source and $\bm{\Theta}$). The preconditioner we apply is based on
the additive combination of a two-level method in angle and the
streaming/removal operator, which we denote as

$\textrm{{M}}^{-1}=\textrm{{M}}_{\textrm{angle}}^{-1}+\textrm{{M}}_{\Omega}^{-1}.$
(13)

The first of these uses a DSA type operator; both DSA and TSA can be thought
of as additive angular multigrid preconditioners, with DSA forming a two-level
multigrid with the coarse level represented by a diffusion equation, whereas
TSA can form multiple coarse grids if desired through lower angular resolution
transport discretisations (e.g., see [6, 14, 16, 17, 18, 20, 40, 41, 24, 25]).
We found both were very effective, but in this work we use a DSA method which
we write as

$\textrm{{M}}_{\textrm{angle}}^{-1}=\textrm{{R}}_{\textrm{angle}}\textrm{{D}}_{\textrm{diff}}^{-1}\textrm{{P}}_{\textrm{angle}},$
(14)

where the angular restrict/prolong are simply the mappings between the
constant moment (i.e., the 0th order scatter term) and our P0 angular space.
$\textrm{{D}}_{\textrm{diff}}$ is a standard DSA diffusion operator with
diffusion coefficient of $\kappa=1/3\sigma_{s}$ and Robin conditions on vacuum
boundaries that we discretise with a CG FEM; this makes our DSA “inconsistent”
but we find this performs well with a Krylov method as the outer iteration, as
in [1]. For further discussions we appeal to the wealth of literature on
different acceleration methods.

The second component of our additive preconditioner is based on the sub-grid
scale streaming/removal operator in (11), namely

$\textrm{{M}}_{\Omega}^{-1}=\left(\textrm{{A}}_{\Omega}-\textrm{{B}}_{\Omega}\hat{\textrm{{D}}}^{-1}\textrm{{C}}_{\Omega}\right)^{-1}$
(15)

The additive combination of these two preconditioners is designed to achieve
good performance in both the streaming and scattering limits. For problems
with both streaming/scattering regions we would apply each of the additive
preconditioners only in regions which require it [42]. It is trivial to form
the grid-transfer operators $\textrm{{R}}_{\textrm{angle}}$ and
$\textrm{{P}}_{\textrm{angle}}$ regardless of angular refinement, but for a
scalable iterative method we must be able to apply the inverses of both the
diffusion matrix with a fixed amount of work given spatial refinement, and the
streaming/removal operator with fixed work given space/angle refinement. We
can do so inexactly given these are applied as preconditioners. This is in
contrast to a source iteration like (4), where [29] notes that inexact
application of $\textrm{{L}}^{-1}$ changes the discrete solution and they show
that increasingly accurate solves are required with grid refinement.

The next section describes a novel reduction multigrid that we can use to
apply the inverses of our streaming/removal operator,
$\textrm{{A}}_{\Omega}-\textrm{{B}}_{\Omega}\hat{\textrm{{D}}}^{-1}\textrm{{C}}_{\Omega}$
(or the streaming operator in the limit of zero total cross-section), and a CG
diffusion operator, $\textrm{{D}}_{\textrm{diff}}$ if desired. For comparison
purposes we also test AIRG on the full matrix,
$\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}$, in the
scattering limit; we cannot form this matrix scalably but it helps demonstrate
that AIRG is applicable in both advective and diffuse problems.

##  4 AIRG multigrid

We begin this section with a summary of a reduction multigrid [43, 44, 45, 28,
26, 46]; we should note that reduction multigrids, block LDU
factorisations/preconditioners and multi-level ILU methods all use similar
building blocks (we discuss this further below). If we consider a general
linear system $\textrm{{A}}\textrm{{x}}=\textrm{{b}}$ we can form a block-
system due to a coarse/fine (CF) splitting as

$\begin{bmatrix}\textrm{{A}}_{\textrm{ff}}&\textrm{{A}}_{\textrm{fc}}\\\
\textrm{{A}}_{\textrm{cf}}&\textrm{{A}}_{\textrm{cc}}\end{bmatrix}\begin{bmatrix}\bm{x_{\textrm{f}}}\\\
\bm{x_{\textrm{c}}}\end{bmatrix}=\begin{bmatrix}\bm{b_{\textrm{f}}}\\\
\bm{b_{\textrm{c}}}\end{bmatrix}.$ (16)

If we write the prolongator and restrictor as

$\textrm{{P}}=\begin{bmatrix}\textrm{{W}}\\\
\textrm{{I}}\end{bmatrix},\quad\textrm{{R}}=\begin{bmatrix}\textrm{{Z}}&\textrm{{I}}\end{bmatrix},$
(17)

then we can form a coarse grid matrix as
$\textrm{{A}}_{\textrm{coarse}}=\textrm{{R}}\textrm{{A}}\textrm{{P}}$ and a
multgrid hierarchy can be built by applying the same technique to
$\textrm{{A}}_{\textrm{coarse}}$. To see how the operators R and P are
constructed we consider an exact two-grid method, where down smoothing (“pre”)
occurs before restriction and up smoothing occurs after prolongation (“post”).
We can write the error at the $i$th step as
$\textrm{{e}}^{i}=\bar{\textrm{{x}}}-\textrm{{x}}^{i}$, where
$\bar{\textrm{{x}}}$ is the exact solution to our linear system. Following
[26], if we have error on the top grid after the down smooth, we can form our
coarse grid residual by computing $\textrm{{R}}\textrm{{A}}\textrm{{e}}^{i}$
and hence the error after a coarse grid solve is
$\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{R}}\textrm{{A}}\textrm{{e}}^{i}$.
The error on the top grid after coarse grid correction, denoted as
$\textrm{{e}}^{i+1}$ is hence

$\textrm{{e}}^{i+1}=(\textrm{{I}}-\textrm{{P}}\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{R}}\textrm{{A}})\textrm{{e}}^{i}.$
(18)

We can consider the error, $\textrm{{e}}^{i}$, to be made up of two
components, one of which is in the range of interpolation and a remainder,
namely

$\textrm{{e}}^{i}=\begin{bmatrix}\textrm{{e}}_{\textrm{f}}\\\
\textrm{{e}}_{\textrm{c}}\end{bmatrix}=\begin{bmatrix}\textrm{{W}}\textrm{{e}}_{\textrm{c}}\\\
\textrm{{e}}_{\textrm{c}}\end{bmatrix}+\begin{bmatrix}\bm{\delta}\textrm{{e}}_{\textrm{f}}^{i}\\\
\textrm{{0}}\end{bmatrix}.$ (19)

We can then write (18) as

$\textrm{{e}}^{i+1}=(\textrm{{I}}-\textrm{{P}}\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{R}}\textrm{{A}})\begin{bmatrix}\bm{\delta}\textrm{{e}}_{\textrm{f}}^{i}\\\
\textrm{{0}}\end{bmatrix}.$ (20)

If $\bm{\delta}\textrm{{e}}_{\textrm{f}}^{i}=\textrm{{0}}$ then F-point error
is in the range of interpolation and the coarse-grid correction is exact, as
$\textrm{{e}}^{i+1}=\textrm{{0}}$. Methods based on ideal restriction however
choose Z to ensure that

$\textrm{{R}}\textrm{{A}}\begin{bmatrix}\bm{\delta}\textrm{{e}}_{\textrm{f}}^{i}\\\
\textrm{{0}}\end{bmatrix}=\textrm{{0}}.$ (21)

This enforces that any error on the top grid that is not in the range of
interpolation does not make it down to the coarse grid. A local form of (21)
is used explicitly to form the restrictor in lAIR. Expanding (21) gives the
condition

$(\textrm{{Z}}\textrm{{A}}_{\textrm{ff}}+\textrm{{A}}_{\textrm{cf}})\bm{\delta}\textrm{{e}}_{\textrm{f}}^{i}=\textrm{{0}},$

and hence
$\textrm{{Z}}=-\textrm{{A}}_{\textrm{cf}}\textrm{{A}}_{\textrm{ff}}^{-1}$ is
known as the “ideal” restrictor that satisfies this condition. The error after
coarse-grid correction (20) then becomes

$\textrm{{e}}^{i+1}=\begin{bmatrix}\textrm{{e}}_{\textrm{f}}^{i}-\textrm{{W}}\textrm{{e}}_{\textrm{c}}^{i}\\\
\textrm{{0}}\end{bmatrix}.$ (22)

Post F-point smoothing to convergence then gives an exact two-grid method.
This is one of the defining characteristics of a reduction multigrid, that
ideal restriction ensures that after coarse-grid correction, the C-point error
is zero and hence F-point smoothing is appropriate; a similar statement can be
used to construct an “ideal” prolongator given by
$\textrm{{W}}=-\textrm{{A}}_{\textrm{ff}}^{-1}\textrm{{A}}_{\textrm{fc}}$.
Both [26, 28] build approximations to the ideal restrictor and then build a
classical one-point prolongator that interpolates each F-point from it’s
strongest C-point connection with injection; by definition this preserves the
constant. [26] show that this is sufficient to ensure good convergence in
advection problems, and [27] prove more general criteria on the required
operators.

With the ideal operators the near-nullspace vectors are naturally preserved
and we don’t have to explicitly provide (or determine through adaptive methods
[47, 48]) the near-nullspace vectors; in fact any error modes (near-nullspace
or otherwise) that are not in the range of interpolation stay on the top grid
and are smoothed, with the rest being accurately restricted to the coarse grid
by construction; if we consider near-nullspace vectors this can be easily
demonstrated by considering that if
$\textrm{{A}}[\textrm{{n}}_{\textrm{f}},\textrm{{n}}_{\textrm{c}}]^{\textrm{T}}\approx\textrm{{0}}$,
where $[\textrm{{n}}_{\textrm{f}},\textrm{{n}}_{\textrm{c}}]$ is a
(partitioned) near-nullspace vector, then with the ideal operators
$\textrm{{A}}_{\textrm{coarse}}\textrm{{n}}_{\textrm{c}}\approx\textrm{{0}}$.
This is in contrast to traditional multigrid methods, which smooth high-
frequency modes preferentially on the top grid, with interpolation
specifically designed to transfer low-frequency modes (i.e., near-nullspace
vectors) to the coarse grid, where they become high-frequency and hence can be
smoothed easily.

[26] show that the error propagation matrix, $\bm{\epsilon}$ of a reduction
multigrid with F point up smooths (and no down smooths) is given by

$\bm{\epsilon}=\textrm{{I}}-\textrm{{M}}_{\textrm{LDU}}^{-1}\textrm{{A}},$
(23)

where $\textrm{{M}}_{\textrm{LDU}}$ is a block LDU factorisation of A given by

$\textrm{{M}}_{\textrm{LDU}}=\begin{bmatrix}\textrm{{I}}&\textrm{{0}}\\\
\textrm{{A}}_{\textrm{cf}}\textrm{{A}}_{\textrm{ff}}^{-1}&\textrm{{I}}\end{bmatrix}\begin{bmatrix}\textrm{{A}}_{\textrm{ff}}&\textrm{{0}}\\\
\textrm{{0}}&\textrm{{S}}\end{bmatrix}\begin{bmatrix}\textrm{{I}}&\textrm{{A}}_{\textrm{ff}}^{-1}\textrm{{A}}_{\textrm{fc}}\\\
\textrm{{0}}&\textrm{{I}}\end{bmatrix},$ (24)

where
$\textrm{{S}}=\textrm{{A}}_{\textrm{cc}}-\textrm{{A}}_{\textrm{cf}}\textrm{{A}}_{\textrm{ff}}^{-1}\textrm{{A}}_{\textrm{fc}}$
is the Schur complement of block $\textrm{{A}}_{\textrm{ff}}$. The inverse of
$\textrm{{M}}_{\textrm{LDU}}$ is given by

$\textrm{{M}}_{\textrm{LDU}}^{-1}=\begin{bmatrix}\textrm{{I}}&\textrm{{W}}\\\
\textrm{{0}}&\textrm{{I}}\end{bmatrix}\begin{bmatrix}\textrm{{A}}_{\textrm{ff}}^{-1}&\textrm{{0}}\\\
\textrm{{0}}&\textrm{{S}}^{-1}\end{bmatrix}\begin{bmatrix}\textrm{{I}}&\textrm{{0}}\\\
\textrm{{Z}}&\textrm{{I}}\end{bmatrix}.$ (25)

The ideal operators give that $\textrm{{A}}_{\textrm{coarse}}=\textrm{{S}}$.
Block LDU factorisations formed with approximate ideal operators have a long
history of being used as preconditioners. The same factorisation can be
applied to S to recover a multilevel LDU method. The benefit to using a block
LDU method, instead of a reduction multigrid is that up F-point smoothing and
coarse grid correction occur additively (as written in (25)), rather than
F-point smoothing occuring after coarse grid correction with a reduction
multigrid (as written in (22)). This is an attractive property in parallel.

If we wish to form a reduction multigrid (or an LDU method) we must
approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$; we denote the approximation as
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}\approx\textrm{{A}}_{\textrm{ff}}^{-1}$
and hence our approximate ideal restrictor and prolongator are given by

$\textrm{{P}}=\begin{bmatrix}-\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}\textrm{{A}}_{\textrm{fc}}\\\
\textrm{{I}}\end{bmatrix},\quad\textrm{{R}}=\begin{bmatrix}-\textrm{{A}}_{\textrm{cf}}\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}&\textrm{{I}}\end{bmatrix}.$
(26)

The strength of our approximate ideal operators allows the use of zero “down”
smoothing iterations, with F-point smoothing only on the “up” cycle after
coarse-grid correction. The authors in [26] however do perform some C-point
smoothing on their “up” cycle (with F-F-C Jacobi) for robustness, but we did
not find that necessary in this work.

Equation (22) suggests a simple choice for the prolongator would be
$\textrm{{W}}=\textrm{{0}}$, but in a multi-level setting where we are
approximating $\textrm{{A}}_{\textrm{ff}}^{-1}$ [28, 26] show this would
require an increasingly accurate approximation with grid refinement. Rather
than build a classical prolongator like [26, 28], as we have an approximation
of $\textrm{{A}}_{\textrm{ff}}^{-1}$ it only costs one extra matmatmult per
level to compute the ideal prolongator, P. We then keep only the largest (in
magnitude) entry, and hence we form a one-point ideal prolongator. This is
like using AIR in conjunction with AIP which [27] suggest could form a
scalable method for non-symmetric systems. This is different to lAIR [28],
where Z is computed directly; computing the ideal prolongator would therefore
require an equivalent calculation for W. We find this is a robust choice for
our prolongator (as it satisfies the approximation properties required by
[27]) while also being very simple as it doesn’t require knowledge of the
near-nullspace.

To perform up F-point smoothing on each level, we use a Richardson iteration
to apply $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$. On each level if we are
smoothing $\textrm{{A}}\textrm{{e}}=\textrm{{r}}$, where r is the residual
computed after the down smooth then

$\textrm{{e}}_{\textrm{f}}^{n+1}=\textrm{{e}}_{\textrm{f}}^{n}+\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}(\textrm{{r}}_{\textrm{f}}-\textrm{{A}}_{\textrm{fc}}\textrm{{e}}_{\textrm{c}}^{n}-\textrm{{A}}_{\textrm{ff}}\textrm{{e}}_{\textrm{f}}^{n}).$
(27)

The coarse-grid error does not change during this process, so we can cache the
result of $\textrm{{A}}_{\textrm{fc}}\textrm{{e}}_{\textrm{c}}^{n}$ during
multiple F-point smooths. The next section details how we construct our
approximation $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$.

###  4.1 GMRES polynomials

Forming good, but sparse approximations to $\textrm{{A}}_{\textrm{ff}}^{-1}$
is the key to a reduction multigrid (and LDU methods as mentioned); this is
achieved in-part by ensuring a “good” CF splitting that results in a well-
conditioned $\textrm{{A}}_{\textrm{ff}}$. In particular
$\textrm{{A}}_{\textrm{ff}}$ is often better conditioned (or more diagonally
dominant) than A. In this section we assume a suitable CF splitting has been
performed; see Section 5 for more details.

The original AMGr work [44] approximated $\textrm{{A}}_{\textrm{ff}}^{-1}$
using the inverse diagonal of $\textrm{{A}}_{\textrm{ff}}$. The nAIR method
presented by [26] on advection-diffusion equations used a matrix-polynomial
approximation of $\textrm{{A}}_{\textrm{ff}}^{-1}$ generated from a truncated
Neumann series; unfortunately this does not converge well in the diffuse limit
(i.e., when $\textrm{{A}}_{\textrm{ff}}$ is not lower-triangular). The work in
[28] however solved dense, local linear systems, which enforce that
$\textrm{{R}}\textrm{{A}}=\textrm{{0}}$ within a certain F-point sparsity
pattern to compute Z directly; this was denoted as lAIR. For advection-
diffusion equations, this showed good performance in both the advection and
diffuse limit. The work of [29] specifically examined the performance of lAIR
in parallel for the BTE, but used on L in (3). [46] showed the use of sparse
approximate inverses (SAIs) [49, 50] on diffusion equations to approximate
$\textrm{{A}}_{\textrm{ff}}^{-1}$ (they also approximated both W and
$\textrm{{A}}_{\textrm{cc}}^{-1}$ for C-point smoothing). The variant of SAI
method used in [46] solved the minimisation problem,
$\textrm{argmin}_{\textrm{{M}}}||\textrm{{I}}-\textrm{{A}}\textrm{{M}}||_{\textrm{F}}$
(which can be written as a separate least-squares problem for each column), to
generate an approximate inverse M, while enforcing the fixed sparsity of
$\textrm{{A}}_{\textrm{ff}}$ on M. For LDU methods, [51] outlined many of the
techniques used; these include using diagonal approximations, or incomplete LU
(often referred to as multi-level ILU [52]) and Cholesky factorisations.

Rather than explicitly approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$, many
multigrid methods such as smoothed aggregation [53, 54, 55] or root-node [56]
methods build their prolongation operators by minimising the energy (in an
appropriate norm) of their prolongator. This is equivalent to solving
$\textrm{{A}}\textrm{{P}}=\textrm{{0}}$ (either column-wise or in a global
sense) and given the equivalent construction of (21) for the prolongator,
these methods can therefore converge to the ideal prolongator. One of the
benefits of forming the ideal operators through a minimisation process is that
constraints on the sparsity of the resulting operator can be applied at all
points through this process. In particular, if this sparsity affects the
ability to preserve near-nullspace vectors, their preservation can be enforced
at each step as the minimisation converges (e.g., see both [57] for an
overview and [58] for a general application to asymmetric systems).

In this work we compute $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ with GMRES
polynomials [59, 60, 61, 62, 63]. Polynomial methods have been used in
multigrid for many years. Examples include, as mentioned, in nAIR [26] to
approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$, in AMLI to approximate the
inverse of the coarse-grid matrix [64, 65] and commonly with Chebychev
polynomials used as smoothers for SPD matrices, among others. The aim of using
GMRES polynomials in this work is to ensure good approximations to
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ that don’t depend on particular
orderings of the unknowns, diagonal, lower-triangular or block structure,
while still allowing good performance in parallel. As such AIRG should be
applicable to many common discretisations of advection-diffusion problems, in
both limits, or any operators where a suitable splitting can produce a well-
conditioned $\textrm{{A}}_{\textrm{ff}}$ that a GMRES polynomial can
approximate. A key difference between this work and previous work with AIR is
that we also reuse these polynomials as our F-point smoothers. This allows us
to build a very simple and strong multigrid method for asymmetric problems.

We begin by summarising the GMRES method. For a general linear system
$\textrm{{A}}\textrm{{x}}=\textrm{{b}}$, where A is $n\times n$, we specify an
initial guess $\textrm{{x}}^{0}=\textrm{{0}}$ and hence an initial residual
$\textrm{{r}}^{0}=\textrm{{b}}$. The Krylov subspace of dimension $m$:
$\textrm{span}\\{\textrm{{b}},\textrm{{A}}\textrm{{b}},\textrm{{A}}^{2}\textrm{{b}},\ldots,\textrm{{A}}^{m-1}\textrm{{b}}\\}$
is used in GMRES to build an approximate solution. This solution at step $m$
can be written as $\textrm{{x}}^{m}=q_{m-1}(\textrm{{A}})\textrm{{b}}$, where
$q_{m-1}(\textrm{{A}})$ is a matrix polynomial of degree $m-1$ known as the
GMRES polynomial. This is the polynomial that minimises the residual
$||\textrm{{r}}^{m}||=||p(\textrm{{A}})\textrm{{r}}^{0}||_{2}$ subject to
$p(0)=1$, where $p(\textrm{{A}})=1-\textrm{{A}}q_{m-1}(\textrm{{A}})$ is known
as the residual polynomial.

The coefficients of $q_{m-1}$ correspond to the required linear combinations
of the Krylov vectors and hence a typical GMRES algorithm can be modified to
generate them. [61, 62] discuss the generation and application of these
polynomials in detail and care must be taken if high order polynomials are
desired. Thankfully we are only concerned with low-order polynomials for use
within our multigrid hierarchy and as such consider two different bases.

As part of a typical GMRES algorithm, we consider a set of orthonormal vectors
which form a basis for our subspace, stored in the columns of the matrix
$\textrm{{V}}_{m}$ (i.e., the Arnoldi basis). Similarly the Krylov vectors
$\textrm{{b}},\textrm{{A}}\textrm{{b}},\ldots$ make up the columns of the
matrix $\textrm{{K}}_{m}$ (i.e., the Krylov basis). Given the vectors in
$\textrm{{V}}_{m}$ and $\textrm{{K}}_{m}$ span the same space we can write
them as linear combinations of each other and hence
$\textrm{{V}}_{m}=\textrm{{K}}_{m}\textrm{{C}}_{m}$, where $\textrm{{C}}_{m}$
is of size $m\times m$. Typically GMRES doesn’t store $\textrm{{C}}_{m}$ but a
small modification can be made to store these values at each GMRES step. The
approximate solution produced by GMRES is given by
$\textrm{{x}}^{m}=\textrm{{x}}^{0}+\textrm{{V}}_{m}\textrm{{y}}_{m}$, where
$\textrm{{y}}_{m}$ comes from the solution of the least-squares problem. The
coefficients for the polynomial $q_{m-1}$ can then be computed through
$(\alpha_{0},\ldots,\alpha_{m-1})^{\textrm{T}}=\textrm{{C}}_{m}\textrm{{y}}_{m}$
(as in [60]).

Equivalently (in exact arithmetic) a QR factorisation of
$\textrm{{K}}_{m+1}=\textrm{{Q}}\textrm{{R}}$ can be computed. If we form the
submatrix $\tilde{\textrm{{R}}}$ from R but without the first column, and note
that $\beta=\textrm{{R}}_{1,1}$ then the polynomial coefficients come from the
solution of the least-squares problem
$(\alpha_{0},\ldots,\alpha_{m-1})^{\textrm{T}}=\textrm{argmin}_{y_{m}}||\beta\textrm{{e}}_{1}-\tilde{\textrm{{R}}}\textrm{{y}}_{m}||_{2}$,
where $\textrm{{e}}_{1}$ is the first column of the $m+1$ identity. GMRES
methods normally don’t use the Krylov basis directly as it is poorly
conditioned when $m\rightarrow\infty$, although this is typically not a
concern at low-order; for example [61] found using the Krylov basis was stable
up to 10th order. In either case our GMRES polynomial of degree $m-1$ is given
by

$q_{m-1}(\textrm{{A}})=\alpha_{0}+\alpha_{1}\textrm{{A}}+\alpha_{2}\textrm{{A}}^{2}+\ldots+\alpha_{m-1}\textrm{{A}}^{m-1}.$
(28)

At each step of a GMRES method, the subspace size, $m$, grows and hence the
GMRES polynomial changes. Polynomial preconditioning methods typically freeze
the size of the subspace and perform a precompute step to generate a GMRES
polynomial of a fixed order. This polynomial is then used as a stationary
preconditioner, often for GMRES itself. In this work we want to approximate
$\textrm{{A}}_{\textrm{ff}}^{-1}$ on each multigrid level; we use GMRES
polynomials with a fixed polynomial order.

To generate the coefficients for our GMRES polynomials, during our multigrid
setup, on each level of our multigrid hierarchy we set $m$ to a fixed value,
assign an initial guess of zero and use a random rhs (see [66, 67, 62]). We
can then chose to use either the Arnoldi or Krylov basis with
$\textrm{{A}}_{\textrm{ff}}$ to form our coefficients. The Arnoldi basis
requires $m$ steps of the modified GMRES described above; this costs $m$
matvecs along with a number of dot products and norms (i.e., reductions) on
each level.

Instead we could use the Krylov basis, $\textrm{{K}}_{m+1}$, which still
requires $m$ matvecs (the setup of a matrix-power kernel may not be worth it
given we only need to compute our low-order polynomial coefficients once), but
in parallel a tall-skinny QR (TSQR) factorisation can be used to decrease the
number of reductions. This relies on QR factorisations of small local blocks
along with a single all-reduce to generate R; we don’t require Q to compute
our polynomial coefficients and hence the small local Q blocks can be
discarded. This is equivalent to modifying a communication-avoiding GMRES with
$s=1$ [68, 67] to generate the coefficients. We tested using both the Arnoldi
and Krylov basis in this work and they generate the same polynomial
coefficients to near round-off error which shows that stability is not a
concern at such low orders. Hence we can use the Krylov basis in parallel and
the generation of the polynomial coefficients can be considered as
communication-avoiding.

Rather than use a random vector, we could form a GMRES polynomial by using a
block method and solving
$\textrm{{Z}}\textrm{{A}}_{\textrm{ff}}=-\textrm{{A}}_{\textrm{cf}}$. That
polynomial would be tailored to computing Z which is not desirable given we
also use our polynomial to perform F-point smoothing and hence want to probe
all the modes of $\textrm{{A}}_{\textrm{ff}}$, as discussed above.

If we wish to use low-order GMRES polynomials to approximate
$\textrm{{A}}_{\textrm{ff}}^{-1}$, we risk not having converged our
approximate ideal operators sufficiently. Similarly, the introduction of fixed
sparsity, drop tolerances or equivalent in
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$, P, R or
$\textrm{{A}}_{\textrm{coarse}}$ may exacerbate this. We examine the impact of
both varying $m$ and sparsifying our operators below in Section 4.2.

####  4.1.1 Fixed sparsity polynomials

In the limit $m\rightarrow n$, the GMRES polynomials converge to
$\textrm{{A}}_{\textrm{ff}}^{-1}$ exactly. Practically this is impossible to
achieve given the storage requirements. Even at low-order however, we would
like to avoid the fill-in that comes from explicitly forming our polynomials
with $m>2$. To do this, we can construct our polynomials by enforcing a fixed
sparsity on each of the matrix powers; for simplicity we chose the sparsity of
$\textrm{{A}}_{\textrm{ff}}$. This is particularly suited to convection
operators given that the fill-in should be small. For example, if we consider
a third-order polynomial (i.e., $m=4$), and denoting the sparsity pattern of
$\textrm{{A}}_{\textrm{ff}}$ as
$S\subset\\{(i,j)\,|\,(\textrm{{A}}_{\textrm{ff}})_{i,j}\neq 0\\}$, we enforce
that

$(\tilde{\textrm{{A}}_{\textrm{ff}}^{2}})_{i,j}=(\textrm{{A}}_{\textrm{ff}}\textrm{{A}}_{\textrm{ff}})_{i,j},\quad(\tilde{\textrm{{A}}_{\textrm{ff}}^{3}})_{i,j}=(\tilde{\textrm{{A}}_{\textrm{ff}}^{2}}\textrm{{A}}_{\textrm{ff}})_{i,j}\quad(i,j)\in
S,$ (29)

and for $(i,j)$ not in $S$, the entries are zero. This is simply computing
$\textrm{{A}}_{\textrm{ff}}^{2}$ with no fill-in, using this approximation
when computing subsequent matrix-powers and again enforcing no fill-in on the
result. Our fixed-sparsity approximation to $\textrm{{A}}_{\textrm{ff}}^{-1}$
with $m=4$ would then be given by

$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}=\alpha_{0}+\alpha_{1}\textrm{{A}}_{\textrm{ff}}+\alpha_{2}\tilde{\textrm{{A}}_{\textrm{ff}}^{2}}+\alpha_{3}\tilde{\textrm{{A}}_{\textrm{ff}}^{3}}\approx
q_{3}(\textrm{{A}}_{\textrm{ff}}).$ (30)

Fixing the sparsity of the matrix powers reduces the memory consumption of our
hierarchy and also allows us to optimise the construction of our polynomial.
For example, with $m=4$ it costs two matrix-matrix additions and two
matmatmults to explicitly construct our polynomial, but given the shared
sparsity, the additions can be performed quickly and the matmatmults with
$m>2$ can share the same row data required when computing
$\textrm{{A}}_{\textrm{ff}}^{2}$. In parallel this means we only require the
communication of off-processor row data once with $m>2$. Computing low-order
GMRES polynomials with fixed sparsity with the Krylov basis therefore only
requires the communication associated with $m$ matvecs, a single all-reduce
and if $m>2$ the matmatmult which produces $\textrm{{A}}_{\textrm{ff}}^{2}$,
regardless of the polynomial order.

This has the potential to scale well, in contrast to some of the other methods
which approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$ described above. We cannot
use nAIR with truncated Neumann series due to the lack of lower triangular
structure in our spatial discretisation, while with lAIR we found we require
greater than distance 2 neighbours for scalability (we examine this in Section
7), but the communication required for this becomes prohibitive. Using ILU
factorisations make parallelisation difficult given the sequential nature of
the underlying Gaussian elimination; approximate ILU factorisations have more
parallelism [69, 70] but they still require triangle solves (which could then
also be approximated with truncated Neumann series for better performance in
parallel, for example). The SAIs used by [46] should scale well, given that if
the fixed sparsity of $\textrm{{A}}_{\textrm{ff}}$ is used, then the formation
of an approximate inverse only requires the same communication as computing
$\textrm{{A}}_{\textrm{ff}}^{2}$; we must however also consider the local cost
of computing a SAI and the effectiveness of its approximation; we examine this
in Section 7.

Typically with AIR the approximations $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$
on each level are thrown away once the grid-transfer operators have been built
and different F-point smoothers are used; as mentioned we save them to use as
smoothers instead. The GMRES polynomials are very strong smoothers (which
don’t require the calculation of any extra dampening parameters) but storing
them explicitly takes extra memory. We examine the total memory consumption of
the AIRG hierarchy in Section 7, but note that given the fixed sparsity of
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$, our experiments with pure streaming
show it has approximately $65\%$ as many non-zeros as R. We could instead
throw away $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ after our grid transfer
operators are computed on each level and perform F-point smoothing by applying
$q_{m-1}(\textrm{{A}}_{\textrm{ff}})$ matrix-free. With $m\leq 2$ the number
of matvecs required is the same and hence
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ could be discarded. With $m>2$
however, applying $q_{m-1}(\textrm{{A}}_{\textrm{ff}})$ matrix-free would
require more matvecs in order to apply the matrix-powers. As such we believe
the (constant sized) extra memory required is easily justified. We also
investigated using Chebyshev polynomials as smoothers. We precomputed the
required eigenvalue estimates with GMRES in order to build the bounding
circle/ellipse given our asymmetric linear systems. We found that smoothing
with these polynomials was far less efficient/robust (and practically more
difficult) than using the GMRES polynomials; indeed using the GMRES
polynomials directly is one of the key messages of [60].

####  4.1.2 Drop tolerance on $\textrm{{A}}_{\textrm{ff}}$

If our original linear system $\textrm{{A}}\textrm{{x}}=\textrm{{b}}$ is not
sparse, neglecting the fill-in with the fixed sparsity polynomials in Section
4.1.1 may not be sufficient by itself to ensure a practical multigrid method.
We can also introduce a drop tolerance to $\textrm{{A}}_{\textrm{ff}}$ that is
applied prior to constructing our polynomial approximations
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$. In general this is not necessary in
the streaming limit, but with scattering we find it helps keep the complexity
low. This is similar to the strong R threshold used in the hypre
implementation of lAIR, which determines strong neighbours prior to the
construction of Z, and in Section 7 we denote it as such. We can either apply
the dropping after the polynomial coefficients are computed (in a similar
fashion to how the fixed sparsity in Section 4.1.1 is applied), or before such
that we are forming a polynomial approximation to a sparsified
$\textrm{{A}}_{\textrm{ff}}$. With scattering, we did not find much of a
difference in convergence so we chose to apply the dropping before, as this
helps reduce the cost of the matvec used to compute the polynomial
coefficients.

###  4.2 Approximations of $\textrm{{A}}_{\textrm{ff}}^{-1}$

(a)

(b)

(c)

Figure 1: Eigenvalue distribution of
$\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{S}}$ on the second level of a
pure streaming problem, where S is the exact Schur complement formed from the
ideal restrictor and prolongator (i.e., the exact coarse-grid matrix) and
$\textrm{{A}}_{\textrm{coarse}}$ is the approximate coarse grid matrix formed
from our approximate ideal restrictor and prolongator. The colours correspond
to different GMRES polynomial orders for approximating
$\textrm{{A}}_{\textrm{ff}}^{-1}$, with $m=1$, $m=2$, $m=3$ and $m=4$.

As mentioned in Section 4.1, the low-order polynomials we use and the fixed
sparsity described in Section 4.1.1 may impact the operators in our multigrid
hierarchy. Given this we examine the impact of our approximations to
$\textrm{{A}}_{\textrm{ff}}^{-1}$ and the resulting operators by considering
the spectrum of $\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{S}}$ and
$\textrm{{A}}_{\textrm{coarse}}$, where $\textrm{{A}}_{\textrm{coarse}}$ is
the coarse matrix formed from our approximate operators and S is the exact
coarse grid matrix.

As an example, we use AIRG on a pure streaming problem with an unstructured
grid, with two levels of uniform refinement in angle (giving 16 angles). Fig.
1 plots the the spectrum of $\textrm{{A}}_{\textrm{coarse}}^{-1}\textrm{{S}}$
and ideally it should approach one. Fig. 1a shows the result of constructing
our coarse grid matrix with increasing orders of our GMRES polynomial, but
without fixing the sparsity of the matrix powers, so we are using
$q_{0}(\textrm{{A}}_{\textrm{ff}})$, $q_{1}(\textrm{{A}}_{\textrm{ff}})$,
$q_{2}(\textrm{{A}}_{\textrm{ff}})$ and $q_{3}(\textrm{{A}}_{\textrm{ff}})$
exactly. We can see increasing the order increases the accuracy of our
resulting coarse grid matrix, as would be expected, with the eigenvalues
converging to one. The radius of a circle that bounds the eigenvalues reduces
from 0.1760 to 0.0072 with $m=1$ to $m=4$, respectively.

Fig. 1b shows the result of introducing both the fixed sparsity matrix powers
discussed in Section 4.1.1 and introducing relative drop tolerances to the
resulting R and P operators. For the restrictor we drop any entry in a row
that is less than 0.1 times the maximum absolute row entry, and for the
prolongator we keep only the largest entry in each row. We can see these
approximations are detrimental to $\textrm{{A}}_{\textrm{coarse}}$ compared
with Fig. 1a, with the eigenvalues further from one. The bounding circle at
all orders has greater radius, with 0.2533 for $m=1$ and 0.0913 for $m=4$.

Finally, Fig. 1c shows the results from using the same fixed sparsity matrix
powers and drop tolerances on R and P while also introducing a relative drop
tolerance on the resulting $\textrm{{A}}_{\textrm{coarse}}$, where we drop any
entry in a row that is less than 0.1 times the maximum absolute row entry. We
again see this further degrades our approximate coarse grid matrix, with the
effect of increasing the order of our GMRES polynomials diminished; the
bounding circle goes from 0.3539 to 0.3532 with $m=1$ to $m=4$, respectively.

It is clear that introducing additional sparsity into our GMRES polynomials
and resulting operators degrades our coarse matrix. We examine this further by
plotting the smallest eigenvalues of $\textrm{{A}}_{\textrm{coarse}}$ and S in
Fig. 2. In the limit of ideal operators we know the near-nullspace vectors are
preserved, but we would like to verify that with an approximate ideal
restrictor and approximate ideal prolongator this is still the case. We can
see that in Fig. 2a that the GMRES polynomials with $m=4$ and fixed sparsity
does an excellent job capturing the smallest eigenvalues. Furthermore
introducing both fixed sparsity to the GMRES polynomial and drop tolerances on
R and P results in reasonable approximations. We can see in Fig. 2b that
introducing the drop tolerances on the resulting
$\textrm{{A}}_{\textrm{coarse}}$ results in small eigenvalues that do not
match the exact coarse grid matrix.

These results indicate that our low-order GMRES polynomials with fixed
sparsity and the introduction of drop tolerances to our approximate ideal R
and P results in an excellent approximation for the coarse grid streaming
operator. It is well known that introducing additional sparsity to multigrid
operators can harm the resulting operators, and it is clear from these results
that introducing a drop tolerance to our coarse grid matrix has the biggest
impact. In a multilevel setting however, we find that doing this often result
in the best performance, with the slight increase in iterations balanced by
the reduced complexity; care must be taken to not make the drop tolerance too
high.

(a)

(b)

Figure 2: Eigenvalue distribution for the smallest eigenvalues of
$\textrm{{A}}_{\textrm{coarse}}$ and S on the second level of a pure streaming
problem, where S is the exact Schur complement formed from the ideal
restrictor and prolongator (i.e., the exact coarse-grid matrix) whose
eigenvalues are denoted with the black “o”. The red symbols are the
eigenvalues of $\textrm{{A}}_{\textrm{coarse}}$, which is the approximate
coarse grid matrix formed from our approximate ideal restrictor and
prolongator given a GMRES polynomial approximation of
$\textrm{{A}}_{\textrm{ff}}^{-1}$ with $m=4$ and different sparsification
applied to operators.

##  5 CF splitting

All multigrid/multilevel methods require the formation of a hierarchy of
“grids”; LDU methods and reduction multigrids like in this work require the
selection of a subset of DOFs defined as “fine” and “coarse”. For asymmetric
linear systems, CF splitting algorithms often result in coarse grids with
directionality (i.e., they result in a semi-coarsening), typically through
heuristic methods that identify strong connections in matrix entries (e.g, see
[71]) with algorithms like CLJP, PMIS, HMIS, etc [72]) or through compatible
relaxation [73, 74, 75]. We would like the CF splitting to produce a well-
conditioned $\textrm{{A}}_{\textrm{ff}}$ on each level without giving a large
grid or operator complexity across the hierarchy. The effectiveness of some of
the approximations used in the literature (described in Section 4) for
$\textrm{{A}}_{\textrm{ff}}^{-1}$ also clearly depend on the sparsity of
$\textrm{{A}}_{\textrm{ff}}$ produced by a CF splitting.

Previous works have used various CF splittings, including those that produce a
maximally-independent set, giving a diagonal $\textrm{{A}}_{\textrm{ff}}$ that
is easily inverted [76]; or if a block-independent set is generated then
$\textrm{{A}}_{\textrm{ff}}$ is block-diagonal and the blocks can be inverted
directly [77, 78]. With a more general CF splitting [52] used ILU
factorisations to approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$. [79] produced
CF splittings specifically for reduction multigrids and LDU methods that are
targeted at producing a diagonally dominant $\textrm{{A}}_{\textrm{ff}}$.

In this work we use traditional CF splitting algorithms (like in [28]) as we
find they perform well enough and parallel implementations are readily
available. Section 7 presents the results from using the lAIR implementation
in hypre and we found that using the Falgout-CLJP algorithm in hypre resulted
in good CF splittings. In order to make fair comparisons, we show results from
using AIRG with the same algorithm.

##  6 Work estimates

One of the key metrics we use to quantify the performance of the iterative
methods tested is the number of Work Units (WUs) required to solve our linear
systems; this is a FLOP count scaled by the number of FLOPs required to
compute a matvec. We present several different WU calculations, each of which
is scaled by a different matvec FLOP count. This is in an attempt to show fair
comparisons against other multigrid methods, along with source iteration.

To begin, we must first establish a FLOP count for all the different
components of our iterative methods. We begin with our definition of the Cycle
Complexity (CC) of AIRG. The CC is the amount of work performed during a
single V-cycle, scaled by the number of nnzs in the top-grid matrix. Our
calculation of the CC includes the work performed during smoothing and grid-
transfer operators; we use our definition of CC and WUs in all the results
below. We define the work required to compute a matvec with our matrices on
each level as $\\{.\\}^{l}$. For an assembled matrix we set this as the nnzs.
This assumes fused-multiply-add (FMA) instructions are available and hence the
cost of multiplying by $\textrm{{A}}_{\textrm{ff}}$ for example is nnzs rather
than $2\times$nnzs (this cost scales out of the CC anyway). If we consider a
general linear system, $\textrm{{A}}\textrm{{x}}=\textrm{{b}}$, the FLOP count
for performing a single V-cycle with $l_{\textrm{max}}$ levels of AIRG is
given by

$\textrm{FLOP}_{\textrm{V}}^{\textrm{AIRG}}=\\{\hat{\textrm{{A}}}^{-1}\\}^{l_{\textrm{max}}}+\sum_{l=1}^{l=l_{\textrm{max}}-1}v_{\textrm{up}}\\{\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}\\}^{l}+v_{\textrm{up}}\\{\textrm{{A}}_{\textrm{ff}}\\}^{l}+\\{\textrm{{A}}_{\textrm{fc}}\\}^{l}+\\{\textrm{{R}}\\}^{l}+\\{\textrm{{P}}\\}^{l},$
(31)

where $v_{\textrm{up}}=2$ is the number of up F-point smooths and we perform
one application of a GMRES polynomial approximation of
$\hat{\textrm{{A}}}^{-1}$ as a coarse grid solve on $l=l_{\textrm{max}}$.

In Section 7, we also show the results from using lAIR in hypre with FCF-
Jacobi smoothing; by default the CC output by hypre doesn’t include all work
associated with smoothing, residual calculation, etc, and hence we recompute
it. Due to the use of FCF-Jacobi, the result of
$\textrm{{A}}_{\textrm{fc}}\textrm{{x}}_{\textrm{c}}$ during the F smooths
cannot be cached (similarly for the C-point smooths). We therefore compute

$\textrm{FLOP}_{\textrm{V}}^{\textrm{hypre}}=\\{\hat{\textrm{{A}}}^{-1}\\}^{l_{\textrm{max}}}+\sum_{l=1}^{l=l_{\textrm{max}}-1}v_{\textrm{up}}(2n_{\textrm{F}}^{l}+n_{\textrm{C}}^{l})+v_{\textrm{up}}\left(2\left(\\{\textrm{{A}}_{\textrm{ff}}\\}^{l}+\\{\textrm{{A}}_{\textrm{fc}}\\}^{l}\right)+\\{\textrm{{A}}_{\textrm{cf}}\\}^{l}+\\{\textrm{{A}}_{\textrm{cc}}\\}^{l}\right)+\\{\textrm{{R}}\\}^{l}+\\{\textrm{{P}}\\}^{l},$
(32)

where $n_{\textrm{F}}^{l}$ and $n_{\textrm{C}}^{l}$ are the number of F and
C-points on level $l$, respectively and given we only use one FCF-Jacobi as
our smoother on each level $v_{\textrm{up}}=1$. The cycle complexity (for
either hypre or AIRG) is then given by

$\textrm{CC}=\frac{\textrm{FLOP}_{\textrm{V}}}{\\{\textrm{{A}}\\}^{1}}.$ (33)

Any matvec that involves scatter should be computed matrix-free and we denote
that with an “mf” subscript. Given our sub-grid scale discretisation, we need
to account for the cost of computing the source on the rhs of (11), the fine-
scale solution $\bm{\Theta}$ and the addition of the coarse and fine scale
solutions to form $\bm{\psi}$. These are given by

$\textrm{FLOP}_{\textrm{source}}=\\{\textrm{{B}}\\}_{\textrm{mf}}+\\{\hat{\textrm{{D}}}^{-1}\\},\quad\textrm{FLOP}_{\textrm{SGS}}=\\{\textrm{{C}}\\}_{\textrm{mf}}+\\{\hat{\textrm{{D}}}^{-1}\\},\quad\textrm{FLOP}_{\bm{\psi}}=\textrm{NDDOFs}$
(34)

The FLOP count of one iteration of our angular preconditioner (14) is given by

$\textrm{FLOP}_{\textrm{angle}}=2\times\textrm{NCDOFs}+4.5\times\\{\textrm{{D}}_{\textrm{diff}}\\}.$
(35)

We investigated using AIRG to invert the diffusion operator, but found it
difficult to beat the default boomerAMG implementation in hypre (i.e., not
lAIR), which is unsurprising given hypre has been heavily optimised for such
elliptic operators. The factor of 4.5 comes from the cycle complexity of
running boomerAMG on a heavily refined spatial grid and as might be expected
we see the cycle complexity plateaus to around this value (hence we have an
upper bound on work on less refined grids).

We must now quantify the cost of our matrix-free matvecs. As mentioned in
Section 2.2 there are numerous ways we could compute such a matvec, depending
on how much memory we have available; Section 3 discussed that we store the
streaming/removal operator, $\textrm{{M}}_{\Omega}$ to precondition with and
hence we use that in our matvec. The additional cost therefore comes with
scattering (which we assume is isotropic in this work) and hence we have

$\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}_{\textrm{mf}}=\\{\textrm{{M}}_{\Omega}\\}+\sum_{i\in\textrm{cg
nodes}}2\times\delta(i)\times\textrm{NCDOFs}(i)+\\\ \sum_{i\in\textrm{dg
nodes}}\delta(i)\times(2\times\textrm{NDDOFs}(i)+2\times
n_{\textrm{nodes}}+1)+\sum_{e\in\textrm{eles}}3\times\delta(i)\times\\{\hat{\textrm{{D}}}^{-1}_{e}\\}$
(36)

where $n_{\textrm{nodes}}$ is the number of spatial nodes on our DG elements
(3 in 2D or 4 in 3D with tri/tets), $\\{\hat{\textrm{{D}}}^{-1}_{e}\\}$ is the
number of non-zeros in our stored block approximation on a given element (the
factor of 3 comes from one application of $\hat{\textrm{{D}}}^{-1}$ and one of
$\textrm{{B}}_{\Omega}$ and $\textrm{{C}}_{\Omega}$ which have the same
sparsity); this sparsity depends on the angles present on each DG node, but is
at a maximum when uniform angle is used and for each angle present on all
nodes of an element we have a $n_{\textrm{nodes}}\times n_{\textrm{nodes}}$
block. $\textrm{N*DOFS}(i)$ is the number of DOFs on an individual (CG or DG)
spatial node $i$ and $\delta(i)$ is 1 or 0 depending on the presence of
scatter; for a CG node this is zero if every element connected to node $i$ has
a zero scatter cross-section, and one otherwise, for a DG node this is zero if
the element containing node $i$ has a zero cross-section, one otherwise. The
calculation in (36) includes the work required to map to/from Legendre space
on both our coarse (CG) and fine (DG) spatial meshes in order to apply the
scatter component in $\textrm{{A}},\textrm{{B}}$ and C. Similar expressions
are used for the individual $\\{\textrm{{B}}\\}_{\textrm{mf}}$ and
$\\{\textrm{{C}}\\}_{\textrm{mf}}$ in (34).

We now have all the components required to calculate our WUs. We begin with a
simple definition, where we use AIRG as a preconditioner on the assembled
matrix, $\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}$, that
could include scatter and is hence non-scalable. If $n_{\textrm{its}}$ is the
number of outer GMRES iterations performed then the total FLOPs are

$\textrm{FLOP}_{\textrm{full}}=n_{\textrm{its}}\left(\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}+\textrm{FLOP}_{\textrm{V}}\right)+\textrm{FLOP}_{\textrm{source}}+\textrm{FLOP}_{\textrm{SGS}}+\textrm{FLOP}_{\bm{\psi}},$
(37)

and hence the WUs

$\textrm{WUs}^{\textrm{full}}=\frac{\textrm{FLOP}_{\textrm{full}}}{\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}}.$
(38)

If instead we use the additively preconditioned iterative method defined in
Section 3 along with our matrix-free matvec, our total FLOPs are

$\textrm{FLOP}_{\textrm{mf}}=n_{\textrm{its}}\left(\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}_{\textrm{mf}}+\textrm{FLOP}_{\textrm{V}}+\textrm{FLOP}_{\textrm{angle}}\right)+\textrm{FLOP}_{\textrm{source}}+\textrm{FLOP}_{\textrm{SGS}}+\textrm{FLOP}_{\bm{\psi}},$
(39)

We then scale these FLOP counts in several different ways. Firstly there is
the WUs required to compute a matrix-free matvec of our sub-grid scale
discretisation, namely

$\textrm{WUs}^{\textrm{mf}}=\frac{\textrm{FLOP}_{\textrm{mf}}}{\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}_{\textrm{mf}}}.$
(40)

In order to make rough comparisons with a traditional DG FEM source iteration
method, we can scale by the work required to compute a matrix-free matvec with
a DG FEM. In order to not unfairly disadvantage a DG discretisation, we assume
the DG streaming operator is stored in memory, and hence the FLOPs required to
compute a single DG matvec is

$\textrm{FLOP}_{\textrm{DG}}=\frac{5}{3}\\{\hat{\textrm{{D}}}^{-1}\\}+\sum_{i\in\textrm{dg
nodes}}\delta(i)\times(2\times\textrm{NDDOFs}(i)+n_{\textrm{nodes}}).$ (41)

The factor of 5/3 comes from the jump terms (in 2D) that are not included in
the nnzs of our sparsified DG matrix $\hat{\textrm{{D}}}^{-1}$. The work units
scaled by this quantity are therefore

$\textrm{WUs}^{\textrm{DG}}=\frac{\textrm{FLOP}_{\textrm{full}}}{\textrm{FLOP}_{\textrm{DG}}}\quad\textrm{or}\quad\frac{\textrm{FLOP}_{\textrm{mf}}}{\textrm{FLOP}_{\textrm{DG}}}.$
(42)

Again we note that with scattering we have
$\\{\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}\\}_{\textrm{mf}}\approx
1.8\times\textrm{FLOP}_{\textrm{DG}}$.

##  7 Results

Outlined below are several examples problems in both the streaming and
scattering limits, designed to test the performance of AIRG and our additively
preconditioned iterative method. We solve our linear systems with GMRES(30) to
a relative tolerance of 1$\times$ 10${}^{\text{-10}}$, with an absolute
tolerance of 1$\times$ 10${}^{\text{-50}}$ and use an initial guess of zero
unless otherwise stated. We should note that we use AIRG (and lAIR) on our
matrices without relying on any (potential) block structure, for example in
the streaming limit with uniform angle we could use our multigrid on each of
the angle blocks separately. We also don’t scale our matrices; for example we
could view a diagonal scaling as preconditioning the outer GMRES iteration
and/or the GMRES polynomials in our multigrid, but we did not find it
necessary.

When using AIRG we perform zero down smooths and two up F-point smooths (our
C-points remain unchanged). On the bottom level of our multigrid we use one
Richardson iteration to apply a GMRES polynomial approximation of the coarse
matrix. We use a 1-point prolongator which computes W and then keeps the
biggest absolute entry per row. Unless otherwise noted we use 3rd order
($m=4$) GMRES polynomials with fixed sparsity as described in Section 4.1.1.
We only use isotropic scatter in this work. For both AIRG and lAIR, we use the
row-wise infinity norm to define any drop tolerances. When using the lAIR
implementation in hypre, we use zero down smooths and one iteration of FCF-
Jacobi for up smooths, while on the bottom level we use a direct-solve, as we
found these options resulted in the lowest cost/best scaling. We searched the
parameter space to try and find the best values for drop tolerances, strength
of connections, etc when using lAIR and AIRG; these searches were not
exhaustive however, and it is possible there are more optimal values.

All timing results are taken from compiling our code, PETSc 3.15 and hypre
with “-O3” optimisation flag. We compare timing results between our PETSc
implementation of AIRG and the hypre implementation of lAIR. As such we try
and limit the impact of different implementation details. Given this, when
calculating the setup time we exclude the CF splitting time, the time required
to drop entries from matrices (as the PETSc interfaces require us to take
copies of matrices), and to extract the submatrices
$\textrm{{A}}_{\textrm{ff}}$, $\textrm{{A}}_{\textrm{fc}}$, etc. These should
all be relatively low cost parts of the setup, are shared by both AIRG and
lAIR and should scale with the nnzs. The setup time we compare is therefore
that required to form the restrictors, prolongators and coarse matrices. Also
we should note that with AIRG our setup time is an upper bound, as we have not
built an optimised matmatmult for building our fixed sparsity GMRES
polynomials (we know the sparsity of the matrix powers is the same as the two
input matrices). As such we compute a standard matmatmult in PETSc and then
drop entries (although we do provide a flop count for a matmatmult with fixed
sparsity). When timing lAIR, we use the PETSc interface to hypre and hence we
run two solves (and set the initial condition to zero in both). The hypre
setup occurs on the 0th iteration of the first solve, so a second solve allows
us to correctly measure just the solve time. We only show solve times for
problems with pure streaming, as our implementation of the P0 matrix-free
matvec with scattering is not well optimised.

All tabulations of memory used are scaled by the total NDOFs in $\bm{\psi}$ in
(8), i.e., NDOFs=NCDOFs + NDDOFs. Included in this figure is the memory
required to store the GMRES space, both additive preconditioners (if required)
and hence the AIRG hierarchy, $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$,
separate copies of
$\textrm{{A}}_{\textrm{ff}},\textrm{{A}}_{\textrm{fc}},\textrm{{A}}_{\textrm{cf}}$
and $\textrm{{A}}_{\textrm{cc}}$, the diffusion operator and temporary
storage. We do not report the memory use of hypre, although given the operator
complexities it is similar to AIRG.

###  7.1 AIRG on the matrix
$\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}$

In this section, we build the matrix
$\textrm{{A}}-\textrm{{B}}\hat{\textrm{{D}}}^{-1}\textrm{{C}}$ and use this
matrix as a preconditioner, applied with 1 V-cycle of either AIRG or lAIR. As
such we do not use the iterative method described in Section 3, nor do we use
the matrix-free matvec described in Section 2.2. As mentioned using an
iterative method that relies on the full operator is not practical with
scattering given the nonlinear increase in nnzs with angular refinement, but
we wish to show our methods are still convergent in the scattering (diffuse)
limit. Instead, Section 7.2 shows the results from using the iterative method
in Section 3.

Our test problem is a $3\times 3$ box with a source of strength 1 and size
$0.2\times 0.2$ in the centre of the domain. We apply vacuum conditions on the
boundaries and discretise this problem with unstructured triangles and ensure
that our grids are not semi-structured (e.g., we don’t refine coarse grids by
splitting elements). We use uniform level 1 refinement in angle, with 1 angle
per octant (similar to S2).

####  7.1.1 Pure streaming problem

For the pure streaming problem we set the total and scatter cross-sections to
zero. To begin, we examine the performance of AIRG and lAIR with spatial
refinement and a fixed uniform level 1 angular discretisation (i.e, with 4
angles in 2D). Table 1 shows that using distance 1 lAIR in this problem, with
Falgout-CLJP CF splitting results in growth in both the iteration count and
work. We could not find a combination of parameters that results in
scalability with lAIR; increasing the number of FCF smooths to 3 results in an
iteration count with similar growth, namely 15, 14, 14, 17, 19 and 22, but
with a cycle complexity at the finest spatial refinement of 11.2 and hence 271
WUs. Even using distance 2 lAIR didn’t result in scalability, as shown in
Table 2 where we can see a slightly decreased iteration count, with higher
cycle and operator complexities, resulting in a similar number of WUs. Using
both distance 2 lAIR and 3 FCF smooths still results in growth, with iteration
counts of 15, 13, 14, 16, 17 and 19, with a cycle complexity of 15.9 and hence
323 WUs at the finest spatial refinement.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 26 | 2.96 | 1.5 | 106 | 26.3 | -
591 | 1.6$\times$ 10${}^{\text{4}}$ | 25 | 3.5 | 1.77 | 116 | 27.8 | -
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 28 | 3.8 | 1.9 | 138 | 32.6 | -
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 31 | 4.1 | 2.0 | 159 | 37.4 | -
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 36 | 4.2 | 2.1 | 189 | 44.5 | -
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 42 | 4.3 | 2.16 | 225 | 52.6 | -
Table 1: Results from using distance 1 lAIR in hypre on a pure streaming problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP with a strong threshold of 0.2, drop tolerance on A of 0.0075 and R of 0.025 and a strong R threshold of 0.25. CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 26 | 3.2 | 1.58 | 112 | 27.9 | -
591 | 1.6$\times$ 10${}^{\text{4}}$ | 24 | 3.7 | 1.85 | 116 | 28 | -
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 27 | 4.1 | 2.04 | 141 | 33.5 | -
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 30 | 4.4 | 2.16 | 164 | 38.6 | -
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 34 | 4.6 | 2.26 | 192 | 44.9 | -
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 38 | 4.7 | 2.32 | 219 | 51.1 | -
Table 2: Results from using distance 2 lAIR in hypre on a pure streaming
problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP
with a strong threshold of 0.2, drop tolerance on A of 0.0075 and R of 0.025
and a strong R threshold 0.25.

We also tried decreasing the strong R threshold to 1 $\times$
10${}^{\text{-7}}$ in case some neighbours were being excluded, but this
resulted in very little change in the iteration count, while the number of
nnzs in the restrictor (and hence the setup time) grew considerably.
Preliminary investigation suggests we would need to include neighbours at
greater distance than two (along with only F-point smoothing, rather than
FCF), but this increases the setup cost considerably. We also note that using
nAIR (at several different orders) in this case results in a similar iteration
count (even with diagonal scaling of our operators); this is likely due to the
lack of lower triangular structure in our discretisation.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 11 | 5.58 | 1.96 | 83 | 20.4 | 11.9
591 | 1.6$\times$ 10${}^{\text{4}}$ | 10 | 5.85 | 2.48 | 79 | 18.9 | 11.7
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 8 | 6.4 | 2.88 | 70 | 16.6 | 12.2
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 8 | 6.7 | 3.16 | 73 | 17.1 | 12.4
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 9 | 6.9 | 3.36 | 82 | 19.3 | 12.5
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 9 | 7.07 | 3.48 | 84 | 19.6 | 12.7
Table 3: Results from using AIRG with $m=4$ and without fixed sparsity on a
pure streaming problem in 2D with CF splitting by the hypre implementation of
Falgout-CLJP with a strong threshold of 0.2, drop tolerance on A of 0.0075 and
R of 0.025.

Tables 3 & 4 however shows that using AIRG with a third-order GMRES polynomial
and Falgout-CLJP CF splitting results in less work with smaller growth. Using
GMRES polynomials without fixed sparsity requires 84 WUs at the highest level
of refinement, compared to 75 with fixed sparsity. The work in Table 3 has
plateaued, however the work with fixed sparsity is growing slightly with
spatial refinement. Compared to lAIR, we see that AIRG with sparsity control
needs three times less work to solve our pure streaming problem. Rather than
use our one-point ideal prolongator, we also investigated using W with the
same drop tolerances as applied to Z. With fixed-sparsity this resulted in 11,
10, 9, 9, 9 and 10 iterations, with 67, 68, 66, 68, 70, and 78 WUs. The slight
increase in cycle complexity is compensated by using one fewer iterations at
the finest level, but overall we see similar work. In an attempt to ascertain
the effect of using our one-point approximation to the ideal prolongator, we
also constructed a classical one-point prolongator as in [28, 26] where each
F-point is injected from its strongest C-point neighbour. With fixed-sparsity
this results in the same number of iterations and WUs as in Table 4, which
confirms that the classical prolongator is not responsible for the poor
performance of lAIR. The benefit of using a classical one-point prolongator in
this fashion is we require one fewer matmatmults on each level of our setup.
For the rest of this paper we use the one-point approximation to the ideal
operator, but it is clear there is a range of practical operators (with
different sparsification strategies) we could use with AIRG.

Table 4 shows we can solve our streaming problem with the equivalent of
approximately 18 DG matvecs. We can also see that the plateau in the operator
complexity results in almost constant memory use, at 10 copies of the angular
flux. We should also note that we can use AIRG as a solver, rather than as a
preconditioner. We see the same iteration count in pure streaming problems and
the lack of the GMRES space means we only need memory equivalent to approx. 5
copies of the angular flux; this is the amount of memory that would be
required to store just a DG streaming operator in 2D. This shows that our
discretisation and iterative method are low-memory in streaming problems.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 12 | 3.7 | 1.97 | 67 | 16.5 | 10
591 | 1.6$\times$ 10${}^{\text{4}}$ | 10 | 4.1 | 2.5 | 62 | 14.7 | 10
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 9 | 4.4 | 2.87 | 60 | 14.1 | 10.2
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 10 | 4.6 | 3.17 | 67 | 15.8 | 10.3
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 10 | 4.74 | 3.36 | 69 | 16 | 10.4
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 11 | 4.84 | 3.5 | 75 | 17.6 | 10.4
Table 4: Results from using AIRG with $m=4$ and fixed sparsity on a pure
streaming problem in 2D with CF splitting by the hypre implementation of
Falgout-CLJP with a strong threshold of 0.2, drop tolerance on A of 0.0075 and
R of 0.025.

The cost of a solve must also be balanced by the cost of the setup of our
multigrid. Our setup involves computing our GMRES polynomial approximations,
followed by standard AMG operations, namely computing matrix-matrix products
to form our restrictors/prolongators and coarse matrices on each level. Fig. 3
begins by showing the relative amount of work required to compute Z for AIRG
and distance 1 lAIR. For AIRG this is the sum of FLOPs required to compute
$\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ and those to compute
$-\textrm{{A}}_{\textrm{cf}}\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$. For lAIR,
it is more difficult to calculate a FLOP count for the small dense solves
which locally enforce $\textrm{{R}}\textrm{{A}}=\bm{0}$, as it is dependent on
the BLAS implementation; instead we take the size of each of the $n\times n$
dense systems and compute the sum of $n^{3}$. As such we do not try and
compare an absolute measure of the work required by both (we have timing
results below). Fig. 3 instead shows that the growth in work for both AIRG
with fixed sparsity and lAIR begins to plateau with grid refinement; AIRG
without fixed sparsity however grows with refinement, showing the necessity of
the sparsity control on the matrix powers.

Fig. 4a shows that the cost of computing our GMRES polynomial approximations
to $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ with fixed sparsity is constant, at
around 8 FLOPs per DOF. The plateauing growth seen in Fig. 3 therefore comes
from the matmatmult required to compute
$-\textrm{{A}}_{\textrm{cf}}\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$. The FLOP
count with fixed sparsity relies on a custom matmatmult algorithm that assumes
the same sparsity in each component of the matmatmult, as in (29). Computing
the matrix powers with a standard matmatmult and then dropping any fill-in
results in an almost constant number of FLOPs per DOF, at around 9.6, which is
convenient as this makes it less necessary to build a custom matmatmult
implementation. If we do not fix the sparsity of our polynomial approximations
we can see growth in the cost, due to the increasingly dense matrix powers, as
might be expected.

Fig. 4b shows that the cost of our setup with fixed sparsity plateaus, and the
total cost of our setup plus solve grows, due to the increase in work during
the solve shown in Table 4, with 29% growth from lowest to highest refinement.
Without fixed sparsity we see less growth in Table 3 from the solve, but
higher growth in the setup. This results in similar growth in the total work,
but the fixed sparsity case uses approximately 30% fewer FLOPs.

(a)  Figure 3: Sum of FLOPs required across all levels to compute Z, scaled to
the NDOFs, and then relative to the work required on the least refined spatial
grid, for AIRG with $m=4$ (see Table 4) and distance 1 lAIR (see Table 1) with
Falgout-CLJP in a 2D pure streaming problem. The relative scaling is done
separately for each line; they do not all cost the same at the coarsest
resolution, we provide timings in Section 7.1.1.

(a)

(b)

Figure 4: Sum of FLOPs required across all levels during the setup and solve,
scaled to the NDOFs, for AIRG with $m=4$ and Falgout-CLJP in a 2D pure
streaming problem (see Tables 4 and 3).

Fig. 5 shows results for the solve, Z, setup and total time taken to solve our
system with AIRG and lAIR. We can see in Fig. 5a that the results match those
in Tables 1-5, with AIRG with Falgout-CLJP coarsening and fixed sparsity
giving the lowest solve times. We found a difference in the efficiency of our
PETSc implementation vs the hypre implementations however. If we modify the
parameters used with AIRG and increase the work required to solve to roughly
match that distance 1 lAIR in this problem, we found that there was a factor
of almost 2 difference in the solve times.

Fig. 5b shows the total cost of computing the approximate ideal restrictor, Z.
For AIRG we also show the time for computing the GMRES polynomial
approximation to $\textrm{{A}}_{\textrm{ff}}^{-1}$ separately. We can see for
AIRG with fixed sparsity there is slight growth in the time to compute our
polynomials, with much greater growth without fixed sparsity. Given the FLOP
count in Fig. 4a we would expect the time with fixed sparsity to be constant.
The further growth we see in the time to compute Z comes from the matmatmult
to compute $-\textrm{{A}}_{\textrm{cf}}\textrm{{A}}_{\textrm{ff}}^{-1}$, which
Fig. 4b suggests should plateau. We also see growth in the time taken to
compute Z with lAIR, with higher growth for distance 2. Again Fig. 3 shows the
work estimate plateauing with spatial refinement. At the higher spatial
refinements, the Z with AIRG, Falgout-CLJP coarsening and fixed sparsity costs
more to compute than distance 1 lAIR, but less than distance 2. An AIRG
implementation that takes advantage of the shared sparsity of the matrices
when forming $\hat{\textrm{{A}}}_{\textrm{ff}}^{-1}$ could reduce this setup
further, with a reduction in the FLOPs required (by approx. 20% as shown in
Fig. 4a), and also in the cost of any symbolic computation. Given the
substantially improved convergence shown in Table 4 when compared to distance
2 lAIR in Table 2, this shows AIRG is an effective and relatively cheap way to
compute approximate ideal operators in convection problems.

Fig. 5c shows the setup times (which include the times to compute Z from Fig.
5b) and we can see that lAIR has the cheapest overall setup, with fixed
sparsity AIRG coming in the middle, and AIRG without fixed sparsity giving the
highest setup times. We expect AIRG without fixed sparsity to be increasing
given the increased cost of computing the GMRES polynomial with spatial
refinement, as shown in Fig. 4a, but we see increased setup time with AIRG
with fixed sparsity and distance 1 lAIR, whereas the work estimates in Figures
3 and 4b again suggests the setup times should plateau. We found that the time
taken to compute the matmatmults used to build the transfer operators and
coarse matrices is increasing more than the FLOP count would imply. Of course
a FLOP count is not necessarily perfectly indicative of the time required in a
matmatmult given the symbolic compute and memory accesses. It is possible more
efficient matmatmult implementations are available.

Fig. 5d shows the total time taken to setup and solve with our multigrid. Two
of our AIRG results, namely using GMRES polynomials without and without fixed
sparsity manage to beat lAIR. This helps show the promise of GMRES polynomials
as part of an AIR-style multigrid; even without fixed sparsity they can be
competitive. The fixed sparsity AIRG however results in a considerably
decrease in total time, taking approx. 3$\times$ less time than distance 1
lAIR. Even if we try and remove the effect of implementation differences
between lAIR and AIRG discussed above, by equating the solve time of lAIR with
an AIRG result that requires similar WUs, fixed sparsity AIRG still has a
total time of 0.65 $\times$ that of distance 1 lAIR.

(a)

(b)

(c)

(d)

Figure 5: Timings per DOF for AIRG with $m=4$ and lAIR in a 2D pure streaming
problem. The $\times$ is AIRG with fixed sparsity with Falgout-CLJP, the
$\times$ is AIRG without fixed sparsity and Falgout-CLJP, the $\otimes$ is
distance 1 lAIR with Falgout-CLJP, the $\otimes$ is distance 2 lAIR with
Falgout-CLJP.

We also examined using SAIs to build an approximation of
$\textrm{{A}}_{\textrm{ff}}^{-1}$ with the fixed sparsity of
$\textrm{{A}}_{\textrm{ff}}$, like in [46]. We found approximations of
$\textrm{{A}}_{\textrm{ff}}^{-1}$ computed with SAI comparable to our fixed
sparsity GMRES polynomials with $m=4$, with identical convergence behavior
compared to the fixed sparsity AIRG shown in Table 4 (except with the first
spatial refinement). With spatial refinement using SAI had 11, 10, 9, 10, 10
and 11 iterations and hence the same work units and solve times as AIRG given
the matching fixed sparsity. However we found it was more expensive to setup
SAIs on each level when compared to our GMRES polynomials, by approximately
2$\times$; we used the ParaSails implementation in hypre for comparisons. In
particular the fixed sparsity GMRES polynomials require fewer FLOPs as the
nnzs in $\textrm{{A}}_{\textrm{ff}}$ grow compared to fixed sparsity SAI. If
$n_{\textrm{F}}$ is the number of F-points, we require $m$ matvecs which scale
linearly with the nnzs and a single QR factorisation of size
$n_{\textrm{F}}\times(m+1)$, which doesn’t depend on the nnzs. If $m>2$, we
must also compute $m-2$ fixed sparsity matrix powers at cost
$(m-2)s_{\textrm{r}}s_{\textrm{c}}n_{\textrm{F}}$, where $s_{\textrm{r}}$ and
$s_{\textrm{c}}$ are the average number of nnzs per row and column of
$\textrm{{A}}_{\textrm{ff}}$, respectively. If we assume
$s=s_{\textrm{r}}\approx s_{\textrm{c}}$ then the cost of computing the
matrix-powers can be written as
$s(m-2)\times\textrm{nnzs}(\textrm{{A}}_{\textrm{ff}})$. In comparison, the
SAI algorithm with the fixed sparsity of $\textrm{{A}}_{\textrm{ff}}$ requires
solving $n_{\textrm{F}}$ (local) least-squares problems. If we just consider
the required $n_{\textrm{F}}$ QR factorisations of size $s\times s$, this
requires a FLOP count of
$s^{2}\times\textrm{nnzs}(\textrm{{A}}_{\textrm{ff}})$ (we have dropped the
constants). For many problems a low polynomial order is acceptable and hence
$(m-2)\ll s$; this is particularly true on lower multigrid levels where the
average row or column sparsity can grow. There may be a scale, however, at
which the extra local cost of setting up SAIs is balanced by the extra
communication required by our fixed sparsity GMRES polynomials. With $m=4$ we
require the communication associated with four matvecs, a single all-reduce
and computing $\textrm{{A}}_{\textrm{ff}}^{2}$, whereas SAI only requires that
of $\textrm{{A}}_{\textrm{ff}}^{2}$. One of the benefits of using our GMRES
polynomials however is that we can decrease the amount of communication
required by decreasing the polynomial order below 2.

Given this, we examine the role of changing the GMRES polynomial order with
AIRG and fixed sparsity, from zero to four ($m=1$ to $m=5$; the zeroth and
first order polynomials implicitly have fixed sparsity) in Fig. 6. We can see
the zeroth order polynomial is very cheap to construct, but results in the
highest total time. This is because the iteration count is higher. The
increasingly higher order polynomials take longer to setup, although the
difference between successive orders decreases, due to the shared sparsity.
With increasing polynomial order, on the most refined spatial grid we have 40,
16, 12, 11 and 11 iterations to solve, with cycle complexities of 3.25, 4.67,
4.81, 4.84 and 4.85, respectively. We can see that all the polynomials between
first and third order result in similar total times; the decreased cost of
setup for the lower orders is balanced by the increase in iterations. We see
however very similar growth in iteration count with spatial refinement with
the first through fourth order polynomials; for example with first order GMRES
polynomials ($m=2$) in Fig. 6 we have 15, 14, 12, 14, 14 and 16 iterations.
This gives a higher overall amount of work than with $m=4$ (up to
approximately 100 WUs at the highest refinement), but given the similar growth
in iteration count and reduced communication in parallel during the setup,
requiring only two matvecs and a single all-reduce, this may be a good choice
in parallel.

(a)

(b)

Figure 6: Timings per DOF for AIRG with fixed sparsity, Falgout-CLJP and with
varying GMRES polynomial order in a 2D pure streaming problem. The $\times$ is
AIRG with $m=1$, o is $m=2$, $\otimes$ is $m=3$, $\square$ is $m=4$,
$\diamond$ is $m=5$.

Given the results with spatial refinement above, we chose to examine the role
of angular refinement with AIRG with fixed sparsity and $m=4$. Table 5 shows
that with 3 levels of angular refinement on the third refined spatial grid,
the iteration count increases slightly from 9 to 11, but the cycle and
operator complexity are almost identical and hence we have fixed memory
consumption. We do not show the timings as they scale as would be expected.
Table 6 shows that distance 1 lAIR requires a fixed amount of work with
angular refinement, but this is roughly twice that of AIRG.

CG nodes | Angle lvl. | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---|---
2313 | 1 | 6.3$\times$ 10${}^{\text{4}}$ | 9 | 4.4 | 2.87 | 60 | 14.1 | 10.2
2313 | 2 | 2.5$\times$ 10${}^{\text{5}}$ | 10 | 4.4 | 2.88 | 65 | 15.4 | 10.4
2313 | 3 | 1$\times$ 10${}^{\text{6}}$ | 11 | 4.4 | 2.87 | 70 | 16.7 | 10.4
Table 5: Results from using AIRG with $m=4$ and fixed sparsity on a pure streaming problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP with a strong threshold of 0.2, drop tolerance on A of 0.0075 and R of 0.025 with different levels of angular refinement. CG nodes | Angle lvl. | NDOFs | $n_{\textrm{its}}$ | CC | Op Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---|---
2313 | 1 | 6.3$\times$ 10${}^{\text{4}}$ | 28 | 3.8 | 1.9 | 138 | 32.6 | -
2313 | 2 | 2.5$\times$ 10${}^{\text{5}}$ | 27 | 3.7 | 2.36 | 133 | 31.7 | -
2313 | 3 | 1$\times$ 10${}^{\text{6}}$ | 28 | 3.8 | 2.35 | 133 | 31.7 | -
Table 6: Results from using distance 1 lAIR in hypre on a pure streaming
problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP
with a strong threshold of 0.2, drop tolerance on A of 0.0075 and R of 0.025
and a strong R threshold of 0.25 with different levels of angular refinement.

####  7.1.2 Scattering problem

To test the performance of AIRG with diffusion, we set the total and
scattering cross-section to 10.0. Tables 7 and 8 show that both distance 1 and
distance 2 lAIR, respectively, perform similarly, with approximately 3$\times$
growth in WUs from the least to most refined spatial grid. Both the cycle and
operator complexities have plateaued though.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 23 | 2.2 | 1.3 | 77 | 49 | -
591 | 1.6$\times$ 10${}^{\text{4}}$ | 27 | 3.0 | 1.9 | 112 | 69 | -
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 30 | 3.3 | 2.2 | 134 | 82 | -
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 34 | 3.4 | 2.2 | 153 | 93 | -
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 41 | 3.4 | 2.2 | 183 | 110 | -
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 54 | 3.3 | 2.2 | 238 | 144 | -
Table 7: Results from using distance 1 lAIR in hypre on a pure scattering problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP with a strong threshold of 0.9, drop tolerance on A of 1$\times$ 10${}^{\text{-4}}$, R of 1$\times$ 10${}^{\text{-2}}$ and strong R threshold of 0.4. CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 22 | 2.5 | 1.5 | 80 | 51 | -
591 | 1.6$\times$ 10${}^{\text{4}}$ | 25 | 3.5 | 2.2 | 117 | 72 | -
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 27 | 3.8 | 2.4 | 133 | 81 | -
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 31 | 3.7 | 2.4 | 150 | 91 | -
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 37 | 3.7 | 2.5 | 178 | 108 | -
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 47 | 3.7 | 2.5 | 225 | 136 | -
Table 8: Results from using distance 2 lAIR in hypre on a pure scattering
problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP
with a strong threshold of 0.9, drop tolerance on A of 1$\times$
10${}^{\text{-4}}$, R of 1$\times$ 10${}^{\text{-2}}$ and strong R threshold
of 0.4.

AIRG performs simliarly to lAIR in this problem, as shown in Tables 9 without
fixed sparsity and 10 with fixed sparsity, with around 3$\times$ growth and
similar number of WUs. AIRG with fixed sparsity results in slightly lower
operator complexities, but both methods result in memory use of around 20
copies of the angular flux; this is higher than that in the streaming limit as
the full matrix with scattering and a uniform angular discretisation at one
level of refinement has 4$\times$ the nnzs as that in the streaming limit.
There is not a great deal of difference in the cycle complexity between AIRG
with and without fixed sparsity; this is because the strong R threshold of 0.4
results in a very sparse version of $\textrm{{A}}_{\textrm{ff}}$ being used to
construct the GMRES polynomials (and hence very little fill-in relative to the
top grid matrix). We can decrease the strong R tolerance to decrease the
iteration count (in both AIRG and lAIR), but the nnzs in (or equivalently the
number of neighbours used to construct) Z grows considerably, as might be
expected with scattering.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 22 | 1.9 | 1.3 | 68 | 43 | 17.3
591 | 1.6$\times$ 10${}^{\text{4}}$ | 26 | 2.2 | 2.1 | 88 | 54 | 18.0
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 34 | 2.5 | 2.4 | 122 | 74 | 18.7
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 41 | 2.7 | 2.7 | 155 | 94 | 19.5
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 45 | 2.8 | 2.8 | 175 | 106 | 19.9
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 61 | 2.7 | 2.9 | 232 | 140 | 19.5
Table 9: Results from using AIRG with $m=4$ and without fixed sparsity on a pure scattering problem in 2D with CF splitting by the hypre implementation of Falgout-CLJP with a strong threshold of 0.9, drop tolerance on A of 1$\times$ 10${}^{\text{-4}}$, R of 1$\times$ 10${}^{\text{-2}}$ and strong R threshold of 0.4. CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{full}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 22 | 1.9 | 1.3 | 67 | 43 | 17.3
591 | 1.6$\times$ 10${}^{\text{4}}$ | 26 | 2.2 | 2.0 | 88 | 54 | 18.0
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 34 | 2.5 | 2.4 | 122 | 74 | 18.7
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 38 | 2.7 | 2.7 | 144 | 87 | 19.4
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 51 | 2.8 | 2.8 | 196 | 118 | 19.7
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 60 | 2.7 | 2.8 | 224 | 135 | 19.3
Table 10: Results from using AIRG with $m=4$ and fixed sparsity on a pure
scattering problem in 2D with CF splitting by the hypre implementation of
Falgout-CLJP with a strong threshold of 0.9, drop tolerance on A of 1$\times$
10${}^{\text{-4}}$, R of 1$\times$ 10${}^{\text{-2}}$ and strong R threshold
of 0.4.

Fig. 7 shows the timing results from AIRG and lAIR in this problem. Again we
see a curious implementation difference in solve times in Fig. 7a, as both
lAIR and AIRG require simliar number of WUs, but the hypre implementation of
lAIR requires roughly twice the time to solve. Fig. 7b shows that the time to
compute the GMRES polynomial for AIRG is largely constant, with the time to
compute Z again between distance 1 and distance 2 lAIR; this is also true for
the total setup time in Fig. 7c. Fig. 7d shows that our AIRG implementation is
the cheapest method overall, taking around 0.4$\times$ the amount of time as
lAIR to solve this problem. If we again equate the solve time between lAIR and
AIRG given the similar amount of work required, lAIR and AIRG perform
similarly, each requiring approximately 2.9 $\mu$s total time per DOF.

(a)

(b)

(c)

(d)

Figure 7: Timings per DOF for AIRG with $m=4$ and lAIR in a 2D pure scattering
problem. The $\times$ is AIRG with fixed sparsity with Falgout-CLJP, the
$\times$ is AIRG without fixed sparsity and Falgout-CLJP, the $\otimes$ is
distance 1 lAIR with Falgout-CLJP, the $\otimes$ is distance 2 lAIR with
Falgout-CLJP.

Fig. 8 shows the results from changing the GMRES polynomial order and similar
trends to that in the streaming limit can be seen, namely the 0th order
polynomial is very cheap to setup, but results in the highest total time.
Indeed the 0th order polynomial did not converge at the two highest spatial
refinements. The first through fourth order polynomials all result in similar
total times, with 62, 61, 60 and 64 iterations, respectively. This result
indicates that the higher polynomial order does not necessarily help decrease
the iteration count with scattering. This is because the strong R threshold is
so high; decreasing this makes the effect of the polynomial order (and the
fixed sparsity) much more pronounced, but we found the lowest overall total
times by allowing heavy dropping. Given using the full matrix is not scalable
with angular refinement, we don’t show convergence results in that case;
instead the next section uses the iterative method defined in Section 3 on
this scattering problem.

(a)

(b)

Figure 8: Timings per DOF for AIRG with fixed sparsity, Falgout-CLJP and with
varying GMRES polynomial order in a 2D pure scattering problem. The $\times$
is AIRG with $m=1$, o is $m=2$, $\otimes$ is $m=3$, $\square$ is $m=4$,
$\diamond$ is $m=5$.

###  7.2 Additively preconditioned iterative method

In this section we show the performance of the iterative method from Section 3
in the scattering limit. The results in this section are not designed to show
the performance of a standard DSA method with different scattering ratios,
optical cell lengths, etc; we appeal to the wealth of literature on the topic.
Instead we wish to show that the additive combination of our preconditioners
is effective and that multigrid methods can be used to invert these operators
scalably. As discussed, this means forming the streaming/removal operator
$\textrm{{M}}_{\Omega}$, a CG diffusion operator
$\textrm{{D}}_{\textrm{diff}}$ and the streaming/removal components
$\textrm{{B}}_{\Omega}$ and/or $\textrm{{C}}_{\Omega}$, all of which can be
done scalably. We use 1 V-cycle of AIRG with the same drop/strong tolerances
as in Section 7.1.1 to apply $\textrm{{M}}_{\Omega}^{-1}$ and 1 V-cycle of
boomerAMG (with default options) to apply $\textrm{{D}}_{\textrm{diff}}^{-1}$
per outer GMRES iteration.

(a)

(b)

Figure 9: The 10 biggest and smallest eigenvalues (by real, imaginary parts
and magnitude) (dots) and field of values (solid lines) of different
operators, with the red the equivalent for $\textrm{{A}}_{\textrm{ff}}$ with
CF splitting by Falgout-CLJP. Computed on the third refined spatial grid with
level one angular refinement.

For our iterative method to be effective, 1 V-cycle of both methods must
reduce the error by a fixed amount with space/angle refinement (which is
equivalent to a solve with a fixed tolerance taking a fixed amount of work).
We can assume that multigrid methods such as boomerAMG can invert the
diffusion operator with fixed work (we also scale the diffusion operator by
its inverse diagonal prior to use), but in Section 7.1.1 we only showed that
AIRG can invert the streaming operator with fixed work in the solve, rather
than the streaming/removal operator. Thankfully the removal term results in a
better conditioned matrix given extra term on the (block) diagonals; the
streaming limit is the most difficult to solve. Fig. 9 shows (part of) the
spectrum of the streaming operator vs the streaming/removal operator for the
2D source problem with the third refined grid, level one angular refinement
and total cross-section of 10.0. We can see that the smallest eigenvalues of
the streaming/removal operator are (slightly) further from the origin.

The convergence of AIRG relies on the convergence of our GMRES polynomial
approximations to $\textrm{{A}}_{\textrm{ff}}^{-1}$. We can also see in Fig. 9
that the eigenvalues of $\textrm{{A}}_{\textrm{ff}}$ are more compact than
that of the full operators, confirming that the CF splitting is helping
produce a better conditioned $\textrm{{A}}_{\textrm{ff}}$ in both cases. We
know that our operators are non-normal, so the spectrum does not completely
determine the convergence of our GMRES polynomials [63]. Given this, Fig. 9
also plots the field of values (a.k.a., the numerical range) of our operators,
given by

$\mathcal{F}(\textrm{{A}})=\\{\textrm{{x}}^{*}\textrm{{A}}\textrm{{x}}\,|\,\textrm{{x}}\textrm{{x}}^{*}=1,\textrm{{x}}\in\mathbb{C}^{n}\\},$
(43)

which is a convex set that contains the eigenvalues. To give some insight into
the convergence of the GMRES polynomials, we define $\mu$ to be the distance
from the origin, or

$\mu=\min_{z\in\mathcal{F}(\textrm{{A}})}|z|.$ (44)

[80] show that (see also [63]) if the field of values doesn’t contain the
origin, $\beta\in(0,\pi/2)$ such that $\cos(\beta)=\mu/||\textrm{{A}}||$ and
the Hermitian part of A, namely $(\textrm{{A}}+\textrm{{A}}^{*})/2$ is
positive definite then the residual at step $m$ is bounded by

$||\textrm{{r}}^{m}||\leq||\textrm{{r}}^{0}||\left(2+\frac{2}{\sqrt{3}}\right)(2+\gamma_{\beta})\gamma_{\beta}^{m},$
(45)

where

$\gamma_{\beta}=2\sin\left(\frac{\beta}{4-2\beta/\pi}\right).$ (46)

We confirmed numerically that none of our operators or their fine-fine sub-
matrices touch the origin (and hence the field of values are all in the right-
half of the complex plane) and that their Hermitian parts are positive
definite.

For the $\textrm{{A}}_{\textrm{ff}}$ component of the streaming operator on
the top grid, pictured in Fig. 9 we found that with $m=4$, (45) gives
$||\textrm{{r}}^{m}||\leq 9.41||\textrm{{r}}^{0}||$, while for the
$\textrm{{A}}_{\textrm{ff}}$ component of the streaming/removal operator we
have $||\textrm{{r}}^{m}||\leq 9.40||\textrm{{r}}^{0}||$ (the disk bound in
[81] gives a similar conclusion). These bounds are not particularly tight, but
they do indicate that a 3rd order GMRES polynomial should result in a smaller
residual for the streaming/removal operator and hence we would expect AIRG to
perform better.

We should note that as $\mu\rightarrow 0$, $\beta$ gets closer to $\pi/2$ and
the asymptotic convergence factor, $\gamma_{\beta}^{m}\rightarrow 1$. In
general this means the further the minimum field of values is from the origin,
the better the convergence; this is also demonstrated by considering the disk
bound in [81], given by
$|\delta/c|=(1-\cos(\beta))/(1+\cos(\beta))<\gamma_{\beta}$, where $\delta$
and c are the radius and centre of a disk, respectively, that covers
$\mathcal{F}(\textrm{{A}})$. This helps explain why using GMRES polynomials to
approximate $\textrm{{A}}_{\textrm{ff}}$ can be effective even when GMRES
polynomial preconditioning of the full operators may not be; Fig. 9a shows
that the field of values for the streaming operator almost touches the origin,
with $\mu\approx 3.2\times\text{10}^{\text{-5}}$ and hence single-level GMRES
polynomial preconditioning would likely not be effective in this problem (this
is backed by numerical experiments; we find considerable growth in the
iteration count with refinement). This hints at the importance of combining
GMRES polynomials with a reduction multigrid.

CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{mf}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 23 | 3.0 | 1.6 | 33 | 61 | 16.5
591 | 1.6$\times$ 10${}^{\text{4}}$ | 24 | 4.0 | 1.0 | 36 | 67 | 17.2
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 25 | 4.1 | 1.7 | 38 | 69 | 17.2
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 26 | 4.2 | 2.4 | 39 | 72 | 17.2
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 26 | 4.4 | 2.8 | 39 | 73 | 17.3
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 26 | 4.6 | 3.2 | 39 | 73 | 17.5
Table 11: Results from using additive preconditioning on a pure scattering problem with total and scattering cross-section of 10.0 in 2D. The cycle and operator complexity listed are for AIRG on $\textrm{{M}}_{\Omega}$ with CF splitting by Falgout-CLJP. CG nodes | Angle lvl. | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{mf}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---|---
2313 | 1 | 6.3$\times$ 10${}^{\text{4}}$ | 25 | 4.1 | 1.7 | 38 | 69 | 17.2
2313 | 2 | 2.5$\times$ 10${}^{\text{5}}$ | 27 | 4.0 | 1.4 | 40 | 71 | 16.5
2313 | 3 | 1$\times$ 10${}^{\text{6}}$ | 28 | 4.1 | 1.4 | 41 | 73 | 16.2
Table 12: Results from using additive preconditioning on a pure scattering problem with total and scattering cross-section of 10.0 in 2D with angle refinement. The cycle and operator complexity listed are for AIRG on $\textrm{{M}}_{\Omega}$ with CF splitting by Falgout-CLJP. CG nodes | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{mf}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---
97 | 2.4$\times$ 10${}^{\text{3}}$ | 18 | 3.6 | 1.9 | 27 | 50 | 17.1
591 | 1.6$\times$ 10${}^{\text{4}}$ | 18 | 4.0 | 2.4 | 27 | 50 | 17.2
2313 | 6.3$\times$ 10${}^{\text{4}}$ | 18 | 4.3 | 2.8 | 28 | 51 | 17.4
9166 | 2.5$\times$ 10${}^{\text{5}}$ | 19 | 4.5 | 3.1 | 29 | 54 | 17.5
35784 | 9.9$\times$ 10${}^{\text{5}}$ | 19 | 4.7 | 3.3 | 30 | 55 | 17.6
150063 | 4.2$\times$ 10${}^{\text{6}}$ | 19 | 4.8 | 3.5 | 30 | 55 | 17.7
Table 13: Results from using additive preconditioning on a pure scattering problem with total and scattering cross-section of 1.0 in 2D. The cycle and operator complexity listed are for AIRG on $\textrm{{M}}_{\Omega}$ with CF splitting by Falgout-CLJP. CG nodes | Angle lvl. | NDOFs | $n_{\textrm{its}}$ | CC | Op. Complx | WUs${}^{\textrm{mf}}$ | WUs${}^{\textrm{DG}}$ | Memory
---|---|---|---|---|---|---|---|---
2313 | 1 | 6.3$\times$ 10${}^{\text{4}}$ | 18 | 4.3 | 2.8 | 28 | 51 | 17.4
2313 | 2 | 2.5$\times$ 10${}^{\text{5}}$ | 18 | 4.3 | 2.8 | 27 | 49 | 16.7
2313 | 3 | 1$\times$ 10${}^{\text{6}}$ | 19 | 4.3 | 2.8 | 29 | 51 | 16.5
Table 14: Results from using additive preconditioning on a pure scattering
problem with total and scattering cross-section of 1.0 in 2D with angle
refinement. The cycle and operator complexity listed are for AIRG on
$\textrm{{M}}_{\Omega}$ with CF splitting by Falgout-CLJP.

We see in Table 11 our additively preconditioned iterative method is effective
with a total and scatter cross-section of 10, with the iteration count growing
from 23 to a plateau of 26 with spatial refinement. The work is very close to
constant, with fixed iteration count and slight growth in the cycle
complexity, even though we used AIRG with fixed sparsity. As such we didn’t
investigate using AIRG without fixed sparsity, as the fixed sparsity was
sufficient to give plateauing work. This helps confirm the observations above,
namely that AIRG is more effective on the streaming/removal operator.

Comparing to the results in Section 7.1.2, it uses roughly 73 DG WUs, compared
to around 135 when using either AIRG or lAIR as a preconditioner on the full
matrix. It also uses less memory at approximately 18 copies of the angular
flux, even though we have to store the diffusion operator and several extra
temporary vectors. Compared to the pure streaming problem, this method
requires approximately 4.1$\times$ more work; we can see in Table 11 that this
work largely comes from computing the matrix-free matvec with scattering, with
26 iterations at the highest spatial refinement requiring 39
WUs${}^{\textrm{mf}}$. The split of work is 26 WUs${}^{\textrm{mf}}$ in the
matvec required by the outer GMRES, 11 WUs${}^{\textrm{mf}}$ to apply the
additive preconditioners and around 2 WUs${}^{\textrm{mf}}$ to compute the
source and $\bm{\Theta}$. Similarly, Table 13 shows a (lower) constant
iteration count with a lower total and scatter cross-section of 1.0.

Given the streaming/removal operator is easier to solve, lAIR performs better
when used additively to invert $\textrm{{M}}_{\Omega}$, compared with just the
streaming operator. With a total and scatter cross-section of 10.0, we find
both distance 1 and 2 lAIR give 30, 33, 31, 34, 36 and 39 iterations with
spatial refinement, with similar grid and operator complexities to AIRG.
Increasing the number of FCF smooths from 1 to 3 results in an iteration count
with less growth, namely 30, 25, 25, 26, 26 and 27 iterations, but the cycle
complexity at the finest level of refinement is large at 10.7, compared with
AIRG at 4.6. We also know from Section 7.1.1 that the iteration count of lAIR
grows in the streaming limit, so we do not test lAIR any further as part of
our additive method.

Tables 12 and 14 show that the additive method with AIRG and angular
refinement perform well, as the iteration count and the work required is very
close to constant. Importantly we can also see the memory use is fixed. These
results show that as might be expected, using an (inconsistent) CG DSA can
form an effective preconditioner in scattering problems when used with an
outer GMRES iteration. Importantly the combination of a single V-cycle of AIRG
used to apply the streaming/removal operator and a single V-cycle of a
traditional multigrid on the diffusion operator can be used additively and
results in almost constant work with spatial and angular refinement on
unstructured grids.

##  8 Conclusions

This paper presented a new reduction multigrid based on approximate ideal
restrictors (AIR) combined with GMRES polynomials (AIRG) with excellent
performance in advection-type problems. Matrix polynomial methods have been
used for many years in multilevel methods but we believe we are the first to
use GMRES polynomials in this fashion. Reduction multigrids and LDU methods in
particular benefit from using GMRES polynomials, as the improved conditioning
of $\textrm{{A}}_{\textrm{ff}}$, when compared to A, can allow the formation
of good approximate inverses with low polynomial orders. This allowed us to
easily build both approximate ideal restrictors, approximate ideal
prolongators (without the need to compute near-nullspace vectors) and perform
F-point smoothing (without the need to compute additional dampening
parameters).

GMRES polynomials share many advantages with other polynomial methods; in
particular their coefficients can be computed very simply; low-order
polynomials don’t require additional work to ensure stability (like in [61,
62]); explicitly forming approximate matrix inverses is simple and only
involves matrix-matrix products or if desired; the polynomials can be applied
matrix-free; their application is highly parallel with their setup able to use
communication-avoiding techniques; and they also work well across a range of
symmetric and asymmetric problems.

When applied to the time independent Boltzmann Transport Equation (BTE) we
could solve pure streaming problems (i.e., in the pure advection limit) on
unstructured spatial grids with space/angle refinement with fixed memory use.
The time-independent streaming limit is the most challenging to solve and we
found we could either get fixed work in the solve and growth in the setup, or
by introducing fixed sparsity into the matrix-powers of our GMRES polynomials,
we found fixed work in the setup with growth in the solve. We found good
performance from using between first to fourth order GMRES polynomials on each
level of our multigrid. Fixing the sparsity of our third-order ($m=4$) GMRES
polynomials resulted in a fixed FLOP count in the setup, and building an
implementation of a matmatmult $\textrm{{A}}=\textrm{{B}}\textrm{{C}}$ where
the three matrices share the same sparsity would reduce the implementation
costs of our setup for second order polynomials and higher. With fixed
sparsity we found at most 20% growth in the work to solve with either 6 levels
of spatial refinement or three levels of uniform angular refinement.

A balance must be struck between the scalability of the solve vs expense of
the setup, but we believe this is the first method to show scalable solves
with a stable spatial discretisation that doesn’t feature lower-triangular
structure in the streaming limit of the BTE. We did not spend much effort
tweaking parameters and we have found that that we can get better performance
in these problems with standard AMG tweaks such as level specific drop
tolerances.

We also compared AIRG to two different reduction multigrids and found
performance advantages; one where sparse approximation inverses (SAIs) are
used to approximate $\textrm{{A}}_{\textrm{ff}}^{-1}$, and the lAIR
implementation in hypre. We used ParaSails in hypre to form SAIs with the
fixed sparsity of $\textrm{{A}}_{\textrm{ff}}$ and found almost identical
convergence behavior to fixed sparsity AIRG (and hence the same solve time) in
the streaming limit. We found however that the setup of the SAIs took twice as
long as our GMRES polynomials with $m=4$. For lAIR, we could not find a set of
parameters that resulted in fixed work in the solve. Our further
investigations suggest the combination of distance 3 or 4 lAIR plus (only)
F-point smooths are required to get scalable results with lAIR, but this is
not practical given the setup/communication costs.

In comparison to distance 1 or 2 lAIR, AIRG took roughly two to three times
less work to solve. Timing the setup showed that computing Z with our third-
order GMRES polynomial approximations cost between that of computing Z with
distance 1 and 2 lAIR. The total time of our setup at the highest level of
spatial refinement matched that of distance 2 lAIR. The total time (setup plus
solve) for AIRG was roughly $3\times$ less than lAIR, though implementation
differences make this comparison difficult. We then investigated using AIRG
and lAIR on the full matrix formed with scattering. Forming this matrix cannot
be done scalably with angular refinement, but we showed that AIRG is
applicable in the diffuse limit, performing about as well as lAIR.

We then built an iterative method that used the additive combination of two
preconditioners applied to the angular flux; 1 V-cycle of AIRG was used to
invert the streaming/removal operator and 1 V-cycle of boomerAMG was used to
invert a CG diffusion operator. The streaming/removal operator is easier to
solve than the streaming operator and hence we found the work in the solve
plateaued with fixed sparsity AIRG. Using distance 1 or 2 lAIR to invert the
streaming/removal operator resulted in cycle complexities over twice that of
fixed sparsity AIRG. Given the performance shown here it would be worth
investigating the use of AIRG as part of a standard DG FEM source iteration;
preliminary work reveals AIRG performs similarly when used with DG streaming
or streaming/removal operators.

The only remaining consideration is how we can apply this method with our
previously developed angular adaptivity and the parallel performance, which we
will investigate in future work. AIRG should be performant in parallel, as the
entire multigrid hierarchy can be applied with only matrix-vector products
(i.e., no reductions). The CF splitting algorithm used has a parallel
implementation available in hypre. The GMRES polynomial coefficients on each
level must be computed once during the setup and can be trivially stored for
multiple solves. Furthermore given our use of low-order GMRES polynomials in
AIRG, we found a single step method based on a QR factorisation of the Krylov
basis could be used to generate these coefficients stably. In parallel we
could therefore use a tall-skinny QR and generate the coefficients of a
polynomial of order $m-1$ with $m$ matvecs and a single all-reduce on each
level. For zero and first order polynomials, there is no other communication
required. For second order and higher, the remainder of the GMRES polynomial
setup uses $m-2$ matmatmults and matmatadds to compute matrix-powers. If we
impose the aforementioned fixed sparsity we only need to communicate the
required off-processor rows of $\textrm{{A}}_{\textrm{ff}}$ in the matmatmults
once in order to compute those matrix powers, regardless of the order of the
polynomial.

Given the results in this paper, we believe the combination of our low-memory
sub-grid scale discretisation, AIRG with low-order GMRES polynomials, and an
iterative method that additively preconditions with the streaming/removal
operator and an inconsistent CG DSA forms an excellent method for solving
transport problems on unstructured grids in both the streaming and scattering
limit.

## Acknowledgments

The authors would like to acknowledge the support of the EPSRC through the
funding of the EPSRC grants EP/R029423/1 and EP/T000414/1.

## References

## References

  * Warsa et al. [2004] J. S. Warsa, T. A. Wareing, J. E. Morel,  Krylov Iterative Methods and the Degraded Effectiveness of Diffusion Synthetic Acceleration for Multidimensional SN Calculations in Problems with Material Discontinuities,  Nuclear Science and Engineering 147 (2004) 218–248.
  * Ramone et al. [1997] G. L. Ramone, M. L. Adams, P. F. Nowak,  A Transport Synthetic Acceleration Method for Transport Iterations,  Nuclear Science and Engineering 125 (1997) 257–283. Publisher: Taylor & Francis _eprint: https://doi.org/10.13182/NSE97-A24274.
  * Adams and Larsen [2002] M. L. Adams, E. W. Larsen,  Fast iterative methods for discrete-ordinates particle transport calculations,  Progress in Nuclear Energy 40 (2002) 3–159.
  * Dargaville et al. [2020a] S. Dargaville, A. G. Buchan, R. P. Smedley-Stevenson, P. N. Smith, C. C. Pain,  Scalable angular adaptivity for Boltzmann transport,  Journal of Computational Physics 397 (2020a).
  * Dargaville et al. [2020b] S. Dargaville, R. P. Smedley-Stevenson, P. N. Smith, C. C. Pain,  Goal-based angular adaptivity for Boltzmann transport in the presence of ray-effects,  Journal of Computational Physics 421 (2020b).
  * Manteuffel et al. [1994] T. Manteuffel, S. McCormick, J. Morel, S. Oliveira, G. Yang,  A parallel version of a multigrid algorithm for isotropic transport equations,  SIAM Journal on Scientific Computing 15 (1994) 474–493.
  * Lansrud [2005] B. D. Lansrud, A spatial multigrid iterative method for two-dimensional discrete-ordinates transport problems, Book, Texas A&M University, 2005. Accepted: 2005-08-29T14:39:54Z Artwork Medium: electronic Interview Medium: electronic.
  * Kanschat and Ragusa [2014] G. Kanschat, J. Ragusa,  A Robust Multigrid Preconditioner for $S_n$DG Approximation of Monochromatic, Isotropic Radiation Transport Problems,  SIAM Journal on Scientific Computing 36 (2014) A2326–A2345.
  * Densmore et al. [2016] J. D. Densmore, D. F. Gill, J. M. Pounders,  Cellwise Block Iteration as a Multigrid Smoother for Discrete-Ordinates Radiation-Transport Calculations,  Journal of Computational and Theoretical Transport 0 (2016) 1–26.
  * Nowak [1988] P. F. Nowak, A Coupled Synthetic and Multigrid Acceleration Method for Two-Dimensional Transport Calculations., Ph.D. thesis, University of Michigan, 1988.
  * Morel and Manteuffel [1991] J. E. Morel, T. A. Manteuffel,  An Angular Multigrid Acceleration Technique for Sn Equations with Highly Forward-Peaked Scattering,  Nuclear Science and Engineering 107 (1991) 330–342.
  * Manteuffel et al. [1995] T. Manteuffel, S. Mccormick, J. Morel, S. Oliveira, G. Yang,  A fast multigrid algorithm for isotropic transport problems i: Pure scattering,  SIAM J. Sci. Comp 16 (1995) 601–635.
  * Manteuffel et al. [1996] T. Manteuffel, S. McCormick, J. Morel, G. Yang,  A fast multigrid algorithm for isotropic transport problems. II: With absorption,  SIAM Journal on Scientific Computing 17 (1996) 1449–1474.
  * Pautz et al. [1999] S. D. Pautz, J. E. Morel, M. L. Adams,  An angular multigrid acceleration method for SN equations with highly forward-peaked scattering,  in: Proc. of Int. Conf. Mathematics and Computation, Reactor Physics and Environmental Analysis in Nuclear Application, Madrid, Spain, volume 1, pp. 647–656.
  * Chang et al. [2007] B. Chang, T. Manteuffel, S. McCormick, J. Ruge, B. Sheehan,  Spatial multigrid for isotropic neutron transport,  SIAM Journal on Scientific Computing 29 (2007) 1900–1917.
  * Lee [2010a] B. Lee,  A novel multigrid method for sn discretizations of the mono-energetic boltzmann transport equation in the optically thick and thin regimes with anisotropic scattering, part i,  SIAM Journal on Scientific Computing 31 (2010a) 4744–4773.