$\mathcal{C}^{(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}}$:

$\mathcal{C}^{(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}}\coloneqq\left(V(\mathcal{C}),\left\\{\varphi_{{\mathbf{X}}}(R)\colon{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell},R\in\mathcal{C}[{{\mathbf{X}}}]\right\\}\right).$

This is again a $k$-ary coherent configuration, as it suffices to check this
on the union of at most $k+1$ fibers, and there, $\mathcal{C}$ is
combinatorially isomorphic to
$\mathcal{C}^{(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}}$. Because
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\colon\mathcal{C}\to\mathcal{C}^{(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}}$
is an $\ell$-local isomorphism, and $\mathcal{C}$ is $\ell$-locally
determined, there exists some combinatorial isomorphism inducing
$(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}$. Thus,
$(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}\in\mathbb{C}_{\ell}(\mathcal{C})$.

Conversely, assume
$\mathbb{L}_{\ell}(\mathcal{C})=\mathbb{C}_{\ell}(\mathcal{C})$ and let
$\varphi=(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}\colon\mathcal{C}\to\mathcal{D}$
be an $\ell$-local isomorphism. We need to find a combinatorial isomorphism
$\hat{\varphi}$ that induces $(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}$. For
this, we pick for every fiber $X\in F(\mathcal{C})$ some extension
$X\in{{\mathbf{Y}}}_{X}\in\binom{F(\mathcal{C})}{\ell}$ and define
$\psi_{X}\coloneqq\varphi_{{{\mathbf{Y}}}_{X}}|_{X}\colon X\to X^{\varphi}$.
We combine these maps into a bijection $\psi\colon V(\mathcal{C})\to
V(\mathcal{D})$, which is a combinatorial isomorphism on every fiber, and
furthermore induces
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}$ on
every fiber.

Now, consider the configuration
$\mathcal{D}^{\psi^{-1}}\coloneqq(V(\mathcal{C}),\\{\psi^{-1}(R)\colon
R\in\mathcal{D}\\})$ that we get by pulling back the structure of
$\mathcal{D}$ through $\psi$. This turns $\psi$ into a combinatorial
isomorphism $\mathcal{D}^{\psi^{-1}}\to\mathcal{D}$. In particular,
$\mathcal{D}^{\psi^{-1}}$ is an a $k$-ary coherent configuration. Furthermore,
the family
$(\psi^{-1}\circ\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}$
is an $\ell$-local permutation of $\mathcal{C}$ that fixes every basis
relation inside each fiber. Thus,
$(\psi^{-1}\circ\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\in\mathbb{L}_{\ell}(\mathcal{C})=\mathbb{C}_{\ell}(\mathcal{C})$.
But because also
$(\psi|_{{\mathbf{X}}})_{{\mathbf{X}}}\in\mathbb{C}_{\ell}(\mathcal{C})$, and
$\mathbb{C}_{\ell}(\mathcal{C})$ is closed under composition by Lemma 5.4,
this implies
$(\varphi_{{\mathbf{X}}})_{{\mathbf{X}}}\in\mathbb{C}_{\ell}(\mathcal{C})$.

Before we continue, we need one Lemma from [15]:

######  Lemma 5.8 ([15]).

For every finite group $\Gamma$, there exists an edge-colored directed graph
$G_{\Gamma}$ with the following properties:

  1. 1.

$V(G_{\Gamma})=\Gamma$,

  2. 2.

$\operatorname{Aut}(G_{\Gamma})\cong\Gamma$, where $\Gamma$ acts on
$G_{\Gamma}$ by left-multiplication.

We now show that we can check the second condition of Lemma 5.6 in polynomial
time.

######  Lemma 5.9.

For every $c\in\mathbb{N}$, there exists an algorithm that for every $k$-ary
coherent configuration $\mathcal{C}$ of order $n$ with abelian fibers of size
at most $c$ runs in time $O_{k,c,\ell}(n^{O(\ell)})$, and either

  1. 1.

outputs an $\ell$-local isomorphism
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\colon\mathcal{C}\to\mathcal{D}$
which is not induced by a combinatorial isomorphism, or

  2. 2.

correctly returns that $\mathcal{C}$ is $\ell$-locally determined.

######  Proof 5.10.

By Lemma 5.6, it suffices to decide whether
$\mathbb{L}_{\ell}(\mathcal{C})=\mathbb{C}_{\ell}(\mathcal{C})$. As in the
proof of Lemma 3.17, we first compute in time $O_{k,c,\ell}(n^{O(\ell)})$ a
(vertex- and edge-colored directed) graph $G_{\mathcal{C}}^{\ell}$ of order
$O_{k,c,\ell}(n^{O(\ell)})$ and bounded color multiplicity such that
$\operatorname{Aut}(G_{\mathcal{C}}^{\ell})\cong\mathbb{L}_{\ell}(\mathcal{C})$.
Lemma 3.16 then allows us to compute a generating set $S$ of
$\mathbb{L}_{\ell}(\mathcal{C})$ in time $O_{k,c,\ell}(n^{O(\ell)})$. Then, it
only remains to check that $S\subseteq\mathbb{C}_{\ell}(\mathcal{C})$, which
we can again reduce to graph isomorphism for graphs of bounded color
multiplicity as in the proof of Lemma 3.19.

We define the vertices of the graph $G_{\mathcal{C}}^{\ell}$ by

$V(G_{\mathcal{C}}^{\ell})=\mathop{\dot{\bigcup}}\limits_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])$

and color the different parts in different colors. For every
${{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}$, the set
$\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])$ has size at most $(c!)^{\ell}$, thus
$G_{\mathcal{C}}^{\ell}$ has color multiplicity at most $(c!)^{\ell}$ and
order at most $n^{\ell}\cdot(c!)^{\ell}\in O_{c,\ell}(n^{\ell})$.

To define the edges of $G_{\mathcal{C}}^{\ell}$, we equip every color class
$\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])$ of $G_{\mathcal{C}}^{\ell}$ with the
structure of the graph $G_{\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])}$ from
Lemma 5.8 such that
$\operatorname{Aut}(G_{\mathcal{C}}^{\ell}[\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])])\cong\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])$.
Furthermore, we connect permutations
$\varphi_{{\mathbf{Y}}}\in\mathbb{P}(\mathcal{C}[{{\mathbf{Y}}}])$ and
$\varphi_{{\mathbf{Z}}}\in\mathbb{P}(\mathcal{C}[{{\mathbf{Z}}}])$ with an
edge if they induce the same algebraic isomorphism on
$\mathcal{C}[{{\mathbf{Y}}}\cap{{\mathbf{Z}}}]$.

We claim that

$\displaystyle\mathbb{L}_{\ell}(\mathcal{C})$
$\displaystyle\to\operatorname{Aut}(G_{\mathcal{C}}^{\ell}),$
$\displaystyle(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}}$
$\displaystyle\mapsto(\psi_{{\mathbf{Y}}}\mapsto\varphi_{{\mathbf{Y}}}\circ\psi_{{\mathbf{Y}}})$

is an isomorphism of groups.

We start by showing that the map is well-defined. As each of the maps
$\psi_{{\mathbf{Y}}}\mapsto\varphi_{{\mathbf{Y}}}\circ\psi_{{\mathbf{Y}}}$ is
an automorphism on every color class, it remains to argue that they preserve
edges between color classes. For this, consider an edge
$\psi_{{\mathbf{Y}}}\psi_{{\mathbf{Z}}}$. By the definition of the edge set,
this means that for all basis relations
$R\in\mathcal{C}[{{\mathbf{Y}}}\cap{{\mathbf{Z}}}]$ we have
$\psi_{{\mathbf{Y}}}(R)=\psi_{{\mathbf{Z}}}(R)$, and we need to show that also
$(\varphi_{{\mathbf{Y}}}\circ\psi_{{\mathbf{Y}}})(R)=(\varphi_{{\mathbf{Z}}}\circ\psi_{{\mathbf{Z}}})(R)$.
But this follows from Equation (2). Injectivity of the map is clear, hence it
remains to argue for surjectivity. Due to the structure of the graphs on the
color classes, every automorphism of
$\operatorname{Aut}(G_{\mathcal{C}}^{\ell})$ is of the form
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}$
with $\varphi_{{\mathbf{X}}}\in\mathbb{P}(\mathcal{C}[{{\mathbf{X}}}])$ for
all ${{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}$. Applying this
automorphism to the clique
$\\{\operatorname{id}_{{\mathbf{X}}}\colon{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}\\}\subseteq
V(G_{\mathcal{C}}^{\ell})$ yields that also
$\\{\varphi_{{\mathbf{X}}}\colon{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}\\}$
is a clique in $G_{\mathcal{C}}^{\ell}$. But this means precisely that
$\varphi_{{\mathbf{Y}}}$ and $\varphi_{{\mathbf{Z}}}$ induce the same
algebraic action for all choices of fibers
${{\mathbf{Y}}},{{\mathbf{Z}}}\in\binom{F(\mathcal{C})}{\ell}$ and thus that
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\in\mathbb{L}_{\ell}(\mathcal{C})$.

Because the graph $G_{\mathcal{C}}^{\ell}$ has bounded color multiplicity, we
can compute a generating set $S$ of
$\operatorname{Aut}(G_{\mathcal{C}}^{\ell})\cong\mathbb{L}_{\ell}(\mathcal{C})$
in the required runtime via Lemma 3.16. Now, fix an arbitrary colored variant
$\mathfrak{C}$ of $\mathcal{C}$. Let
$\varphi=(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\in
S$ be arbitrary. There is a unique colored variant $\mathfrak{C}^{\varphi}$ of
$\mathcal{C}$ such that all of the maps $\varphi_{{\mathbf{X}}}$ become color-
preserving when seen as maps
$\mathfrak{C}[{{\mathbf{X}}}]\to\mathfrak{C}^{\varphi}[{{\mathbf{X}}}]$. An
isomorphism $\hat{\varphi}\colon\mathfrak{C}\to\mathfrak{C}^{\varphi}$ is
precisely a combinatorial permutation of $\mathcal{C}$ that induces $\varphi$.
Thus, $\varphi\in\mathbb{C}_{\ell}(\mathcal{C})$ if and only if
$\mathfrak{C}\cong\mathfrak{C}^{\varphi}$. The latter can be decided in time
$O_{k,c,\ell}(n^{O(\ell)})$ by Lemma 3.16. The algorithm tests whether
$\mathfrak{C}\cong\mathfrak{C}^{\varphi}$ for every $\varphi\in S$. If this is
the case, then $S\subseteq\mathbb{C}_{\ell}(\mathcal{C})$. Otherwise, there is
a
$\varphi=(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\in
S$ such that $\mathfrak{C}\not\cong\mathfrak{C}^{\varphi}$. As argued before,
then
$(\varphi_{{\mathbf{X}}})_{{{\mathbf{X}}}\in\binom{F(\mathcal{C})}{\ell}}\colon\mathcal{C}\to\mathcal{C}^{\varphi}$
is an $\ell$-local isomorphism witnessing non-identification.

Finally, the main approximation theorem is implied by Corollary 5.3 and Lemma
5.9. *

##  6 Conclusion

We have shown that when the dimension $k$ is part of the input, the
$k$-$\mathrm{WL}$-equivalence problem and the $k$-$\mathrm{WL}$-identification
problem are co-NP-hard and NP-hard respectively. On the other hand, when the
dimension $k$ is fixed, the equivalence problem is trivially solvable in
polynomial time, and we have shown that the identification problem is solvable
in polynomial time over graphs of color multiplicity $5$. There is no obvious
reason why this should not generalize to arbitrary graphs. In particular, our
NP-hardness proof reduced from the problem of deciding whether the tree-width
of a given graph is at most $k$. For fixed $k$ this problem can be solved even
in linear time [7].

As a first step to generalizing this, we proved an approximation result for
graphs with bounded abelian color classes. It would be interesting to
determine whether one can drop the restriction to abelian color classes or
improve the approximation algorithm to an exact algorithm.

## References

  * [1] Markus Anders and Pascal Schweitzer.  Parallel computation of combinatorial symmetries.  In 29th Annual European Symposium on Algorithms, ESA 2021, September 6-8, 2021, Lisbon, Portugal (Virtual Conference), volume 204 of LIPIcs, pages 6:1–6:18. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2021.  doi:10.4230/LIPICS.ESA.2021.6.
  * [2] Vikraman Arvind, Johannes Köbler, Gaurav Rattan, and Oleg Verbitsky.  Graph isomorphism, color refinement, and compactness.  Comput. Complex., 26(3):627–685, 2017.  doi:10.1007/S00037-016-0147-6.
  * [3] Albert Atserias and Elitza N. Maneva.  Sherali-Adams relaxations and indistinguishability in counting logics.  SIAM J. Comput., 42(1):112–137, 2013.  doi:10.1137/120867834.
  * [4] László Babai.  Graph isomorphism in quasipolynomial time [extended abstract].  In Proceedings of the 48th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2016, Cambridge, MA, USA, June 18-21, 2016, pages 684–697. ACM, 2016.  doi:10.1145/2897518.2897542.
  * [5] László Babai and Ludek Kucera.  Canonical labelling of graphs in linear average time.  pages 39–46, 1979.  doi:10.1109/SFCS.1979.8.
  * [6] László Babai.  Monte-Carlo algorithms in graph isomorphism testing.  Université de Montréal Technical Report, DMS, (79-10), 1979.
  * [7] Hans L. Bodlaender.  A linear-time algorithm for finding tree-decompositions of small treewidth.  SIAM J. Comput., 25(6):1305–1317, 1996.  doi:10.1137/S0097539793251219.
  * [8] Hans L. Bodlaender.  A partial _k_ -arboretum of graphs with bounded treewidth.  Theor. Comput. Sci., 209(1-2):1–45, 1998.  doi:10.1016/S0304-3975(97)00228-4.
  * [9] Hans L. Bodlaender, Édouard Bonnet, Lars Jaffke, Dušan Knop, Paloma T. Lima, Martin Milanič, Sebastian Ordyniak, Sukanya Pandey, and Ondřej Suchý.  Treewidth is NP-complete on cubic graphs.  In 18th International Symposium on Parameterized and Exact Computation, IPEC 2023, September 6-8, 2023, Amsterdam, The Netherlands, volume 285 of LIPIcs, pages 7:1–7:13. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2023.  doi:10.4230/LIPICS.IPEC.2023.7.
  * [10] Béla Bollobás.  Distinguishing vertices of random graphs.  North-holland Mathematics Studies, 62:33–49, 1982.  doi:10.1016/S0304-0208(08)73545-X.
  * [11] Jin-yi Cai, Martin Fürer, and Neil Immerman.  An optimal lower bound on the number of variables for graph identification.  volume 12, pages 389–410, 1992.  doi:10.1007/BF01305232.
  * [12] G. Chen and I. Ponomarenko.  Lectures on Coherent Configurations.  Central China Normal University Press, 2019.  A draft is available at http://www.pdmi.ras.ru/~inp/ccNOTES%20(13.12.23).pdf.
  * [13] Anuj Dawar and David Richerby.  The power of counting logics on restricted classes of finite structures.  In Computer Science Logic, 21st International Workshop, CSL 2007, 16th Annual Conference of the EACSL, Lausanne, Switzerland, September 11-15, 2007, Proceedings, volume 4646 of Lecture Notes in Computer Science, pages 84–98. Springer, 2007.  doi:10.1007/978-3-540-74915-8_10.
  * [14] Zdenek Dvorák.  On recognizing graphs by numbers of homomorphisms.  J. Graph Theory, 64(4):330–342, 2010.  doi:10.1002/JGT.20461.
  * [15] Robert Frucht.  Herstellung von Graphen mit vorgegebener abstrakter Gruppe.  Compositio Mathematica, 6:239–250, 1939.  URL: http://www.numdam.org/item/CM_1939__6__239_0/.
  * [16] Frank Fuhlbrück, Johannes Köbler, and Oleg Verbitsky.  Identifiability of graphs with small color classes by the Weisfeiler-Leman algorithm.  SIAM J. Discret. Math., 35(3):1792–1853, 2021.  doi:10.1137/20M1327550.
  * [17] Merrick Furst, John Hopcroft, and Eugene M. Luks.  A subexponential algorithm for trivalent graph isomorphism.  Technical report, USA, 1980.
  * [18] Leslie M. Goldschlager.  The monotone and planar circuit value problems are log space complete for P.  SIGACT News, 9(2):25–29, 1977.  doi:10.1145/1008354.1008356.
  * [19] Martin Grohe.  Equivalence in finite-variable logics is complete for polynomial time.  volume 19, pages 507–532, 1999.  doi:10.1007/S004939970004.
  * [20] Martin Grohe.  Descriptive Complexity, Canonisation, and Definable Graph Structure Theory, volume 47 of Lecture Notes in Logic.  Cambridge University Press, 2017.  doi:10.1017/9781139028868.
  * [21] Martin Grohe, Moritz Lichter, Daniel Neuen, and Pascal Schweitzer.  Compressing CFI graphs and lower bounds for the Weisfeiler-Leman refinements.  In 64th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2023, Santa Cruz, CA, USA, November 6-9, 2023, pages 798–809. IEEE, 2023.  doi:10.1109/FOCS57990.2023.00052.
  * [22] Martin Grohe and Julian Mariño.  Definability and descriptive complexity on databases of bounded tree-width.  In Database Theory - ICDT ’99, 7th International Conference, Jerusalem, Israel, January 10-12, 1999, Proceedings, volume 1540 of Lecture Notes in Computer Science, pages 70–82. Springer, 1999.  doi:10.1007/3-540-49257-7_6.
  * [23] Martin Grohe and Daniel Neuen.  Canonisation and definability for graphs of bounded rank width.  ACM Trans. Comput. Log., 24(1):6:1–6:31, 2023.  doi:10.1145/3568025.
  * [24] Martin Grohe and Martin Otto.  Pebble games and linear equations.  J. Symb. Log., 80(3):797–844, 2015.  doi:10.1017/JSL.2015.28.
  * [25] Lauri Hella.  Logical hierarchies in PTIME.  Inf. Comput., 129(1):1–19, 1996.  doi:10.1006/INCO.1996.0070.
  * [26] Neil Immerman and Eric S. Lander.  Describing Graphs: A First-Order Approach to Graph Canonization, pages 59–81.  Springer New York, New York, NY, 1990.  doi:10.1007/978-1-4612-4478-3_5.
  * [27] Tommi A. Junttila and Petteri Kaski.  Engineering an efficient canonical labeling tool for large and sparse graphs.  In Proceedings of the Nine Workshop on Algorithm Engineering and Experiments, ALENEX 2007, New Orleans, Louisiana, USA, January 6, 2007. SIAM, 2007.  doi:10.1137/1.9781611972870.13.
  * [28] Tommi A. Junttila and Petteri Kaski.  Conflict propagation and component recursion for canonical labeling.  In Theory and Practice of Algorithms in (Computer) Systems - First International ICST Conference, TAPAS 2011, Rome, Italy, April 18-20, 2011. Proceedings, volume 6595 of Lecture Notes in Computer Science, pages 151–162. Springer, 2011.  doi:10.1007/978-3-642-19754-3_16.
  * [29] Sandra Kiefer, Ilia Ponomarenko, and Pascal Schweitzer.  The Weisfeiler-Leman dimension of planar graphs is at most 3.  J. ACM, 66(6):44:1–44:31, 2019.  doi:10.1145/3333003.
  * [30] Sandra Kiefer, Pascal Schweitzer, and Erkal Selman.  Graphs identified by logics with counting.  volume 23, pages 1:1–1:31, 2022.  doi:10.1145/3417515.
  * [31] Tuukka Korhonen.  A single-exponential time 2-approximation algorithm for treewidth.  In 62nd IEEE Annual Symposium on Foundations of Computer Science, FOCS 2021, Denver, CO, USA, February 7-10, 2022, pages 184–192. IEEE, 2021.  doi:10.1109/FOCS52979.2021.00026.
  * [32] Ludek Kucera.  Canonical labeling of regular graphs in linear average time.  In 28th Annual Symposium on Foundations of Computer Science, Los Angeles, California, USA, 27-29 October 1987, pages 271–279. IEEE Computer Society, 1987.  doi:10.1109/SFCS.1987.11.
  * [33] Brendan D. McKay and Adolfo Piperno.  Practical graph isomorphism, II.  J. Symb. Comput., 60:94–112, 2014.  doi:10.1016/J.JSC.2013.09.003.
  * [34] Christopher Morris, Martin Ritzert, Matthias Fey, William L. Hamilton, Jan Eric Lenssen, Gaurav Rattan, and Martin Grohe.  Weisfeiler and Leman go neural: Higher-order graph neural networks.  In The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019, pages 4602–4609. AAAI Press, 2019\.  doi:10.1609/AAAI.V33I01.33014602.
  * [35] Daniel Neuen and Pascal Schweitzer.  Benchmark graphs for practical graph isomorphism.  In 25th Annual European Symposium on Algorithms, ESA 2017, September 4-6, 2017, Vienna, Austria, volume 87 of LIPIcs, pages 60:1–60:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2017.  doi:10.4230/LIPICS.ESA.2017.60.
  * [36] Daniel Neuen and Pascal Schweitzer.  An exponential lower bound for individualization-refinement algorithms for graph isomorphism.  In Proceedings of the 50th Annual ACM SIGACT Symposium on Theory of Computing, STOC 2018, Los Angeles, CA, USA, June 25-29, 2018, pages 138–150. ACM, 2018.  doi:10.1145/3188745.3188900.
  * [37] Sang-il Oum and Paul D. Seymour.  Approximating clique-width and branch-width.  J. Comb. Theory, Ser. B, 96(4):514–528, 2006.  doi:10.1016/J.JCTB.2005.10.006.
  * [38] Neil Robertson and Paul D. Seymour.  Graph minors. XIII. The disjoint paths problem.  J. Comb. Theory, Ser. B, 63(1):65–110, 1995.  doi:10.1006/JCTB.1995.1006.
  * [39] Kyoungah See and Sung Y. Song.  Association schemes of small order.  Journal of Statistical Planning and Inference, 73(1):225–271, 1998\.  doi:10.1016/S0378-3758(98)00064-0.
  * [40] Paul D. Seymour and Robin Thomas.  Graph searching and a min-max theorem for tree-width.  J. Comb. Theory, Ser. B, 58(1):22–33, 1993.  doi:10.1006/JCTB.1993.1027.
  * [41] B. Weisfeiler and A. Leman.  The reduction of a graph to canonical form and the algebra which appears therein.  Nauchno-Technicheskaya Informatsia, Seriya 2, 9:12–16, 1968.  An english translation due to Grigory Ryabov is available at https://www.iti.zcu.cz/wl2018/pdf/wl_paper_translation.pdf.
  * [42] Faried Abu Zaid, Erich Grädel, Martin Grohe, and Wied Pakusa.  Choiceless polynomial time on structures with small abelian colour classes.  In Mathematical Foundations of Computer Science 2014 - 39th International Symposium, MFCS 2014, Budapest, Hungary, August 25-29, 2014. Proceedings, Part I, volume 8634 of Lecture Notes in Computer Science, pages 50–62. Springer, 2014.  doi:10.1007/978-3-662-44522-8_5.


# Exploring One-Cell Inversion Method
for Transient Transport on GPU

J. P. Morgan1,2, Ilham Variansyah1,3, Todd S. Palmer1,3, and Kyle E.
Niemeyer1,2
1Center for Exascale Monte Carlo Neutron Transport
2School of Mechanical https://cement-psaap.github.io/ Industrial  and
Manufacturing Engineering
Oregon State University  204 Rogers Hall  Corvallis  OR 97331
3School of Nuclear Science and Engineering
Oregon State University  Merryfield Hall  Corvallis  OR 97331
<EMAIL_ADDRESS>
###### Abstract

To find deterministic solutions to the transient $S_{N}$ neutron transport
equation, iterative schemes are typically used to treat the scattering (and
fission) source terms. We explore the one-cell inversion iteration scheme to
do this on the GPU and make comparisons to a source iteration scheme. We
examine convergence behavior, through the analysis of spectral radii, of both
one-cell inversion and source iterations. To further boost the GPU parallel
efficiency, we derive a higher-order discretization method, simple corner
balance (in space) and multiple balance (in time), to add more work to the
threads and gain accuracy. Fourier analysis on this higher-order numerical
method shows that it is unconditionally stable, but it can produce negative
flux alterations that are critically damped through time. We explore a whole-
problem (in all angle and all cell) sparse linear algebra framework, for both
iterative schemes, to quickly produce performant code for GPUs. Despite one-
cell inversion requiring additional iterations to convergence, those
iterations can be done faster to provide a significant speedup over source
iteration in quadrature sets at or below $S_{128}$. Going forward we will
produce a two-dimensional implementation of this code to experiment with
memory and performance impacts of a whole-problem framework including methods
of synthetic acceleration and pre-conditioners for this scheme, then we will
begin making direct comparisons to traditionally implemented source iteration
in production code.

###### keywords:

deterministic transport, one-cell inversion, source iteration, gpu, transient,
dynamic

##  1 INTRODUCTION

To find deterministic solutions to the transient $S_{N}$ (where N is the
number of angles in a quadrature set) neutron transport equation, iterative
schemes are typically used to treat the scattering (and fission) source terms.
The source iteration (SI) method is commonly used to do this, often
accompanied by preconditioners or synthetic accelerators, where the
contribution to the solution from the scattering source is allowed to lag,
while the angular flux is solved in every ordinate via transport sweeps
through the spatial domain [1].

SI sweeps in Cartesian geometries are readily parallelized over the number of
angles, as the source term is known from the previous iteration, allowing the
angular flux in each ordinate to be computed independently. While any
parallelization is a boon to performance, a scheme that is embarrassingly
parallel over the dimension with the greatest number of degrees of freedom
—space— may be advantageous. Furthermore, many $S_{N}$ production codes that
implement SI use some kind of parallel algorithm that works in higher spatial
dimensions (e.g., PARTISN, which implements the Koch–Baker–Alcouffe or KBA
algorithm [2]). In this paper, we explore one-cell inversion (OCI), which is
inherently parallel over space on both CPUs and graphics processor units
(GPUs), and its application for transient slab geometry transport problems.

##  2 ONE-CELL INVERSION (OCI)

In OCI, all angular fluxes within a cell are computed in a single linear
algebra step, assuming that the angular fluxes incident on the surfaces of the
cell are known from a previous iteration. OCI allows for parallelizing over
the number of cells as each cell is solved independently of the others in a
parallel block Jacobi scheme (where each block corresponds to the problem in a
single cell in all angles). Previous explorations of OCI have primarily
focused on steady-state problems [3, 4].

Because there is no communication of information between cells within an
iteration, OCI can require more iterations to converge to a solution for some
of problems. Specifically, as cellular optical thickness goes down, OCI’s
relative performance degrades. Figure 1 illustrates this behavior, showing the
spectral radii of the two iteration schemes as a function of cell thickness
(in mean free path) and the scattering ratio. These values were computed
numerically from an infinite medium problem (via reflecting boundaries) using
steady-state calculations in $S_{4}$. Gauss Legendre quadrature was used in
all presented explorations. The smaller the spectral radius, the faster a
method is converging. The spectral radius for SI depends strongly on the
scattering ratio, and for problems that are many mean free paths in size, it
is nearly independent of cell optical thickness. The spectral radius of OCI
decreases substantially as the optical thickness of the cells increases.

Figure 1: Spectral radii ($\boldsymbol{\rho}$) of OCI (left) and SI (middle)
and the ratio between the two (right), where $\boldsymbol{\Sigma}$ is the
total cross section, $\boldsymbol{\Delta x}$ is the cell width, and
$\boldsymbol{\Sigma_{s}}$ is the scattering cross section

Since both dimensions in Fig. 1 are governed by relationships with the total
cross section ($\Sigma$), altering that value will impact convergence
behavior. As the scattering ratio decreases, both iterative schemes converge
more quickly. However, the spectral radius of OCI also decreases with
increasing optical thickness, which is an added benefit. When solving
optically thick and highly scattering problems, small increases in $\Sigma$
can drastically improve the relative performance of OCI in comparison with SI.

sectiontocBACKWARD EULER TRANSIENT IMPLEMENTATION If we use a backward Euler
scheme to treat the time derivative in the NTE, the total cross section is
altered so that

$\hat{\Sigma}=\Sigma+\frac{1}{v\Delta t},$ (1)

where $v$ is the velocity of the particles and $\Delta t$ is the time-step.
When a backward Euler discretization is implemented for the transient part of
the NTE, the convergence behavior is improved through increasing the number of
mean free paths and decreasing the scattering ratio, moving where a given
problem lies on Fig. 1. This algorithmic boon to convergence behavior is then
compounded with the parallelization benefits of one-cell inversion and could
make a backward Euler code with one-cell inversion preform better overall (in
decreasing wall time) on many-core machines than a that same code with a
source iteration implementation.

Consider a $10\text{\,}\mathrm{cm}$ thick slab wall where $\Sigma=$
$0.5\text{\,}{\mathrm{cm}}^{-1}$, $\Delta x=$ $0.5\text{\,}\mathrm{cm}$, and
$v=$ $2\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$ with initially no flux
within. Then, at $t=$ $0\text{\,}\mathrm{s}$ the wall is subjected to a
constant isotropic incident flux of magnitude 1 at $x=$
$0\text{\,}\mathrm{cm}$ with a vacuum condition at $x=L$. We initialize our
iterative methods assuming 0 angular or scalar flux everywhere, then iterate
until an error of $1\text{\times}{10}^{-9}$ between the current solution and
the last iteration’s. Fig. 2 shows the ratio of the number of iterations one-
cell inversion requires over source iterations, for this transient problem
across choices of scattering ratio and time-step.

Figure 2: Ratios of iterations required by one-cell inversion as compared to
source iterations at a constant $\boldsymbol{\Delta x=0.5}$ and
$\boldsymbol{\Sigma=0.5}$

While one-cell inversion require more iterations than source iterations for
many choices of $\Delta t$, as the time-step shrinks, one-cell inversion
improve convergence behavior faster than source iterations do, resulting in a
higher value in Fig. 2. Eventually it requires even fewer iterations than
source iterations for small enough time-steps (everywhere $N_{SI}/N_{OCI}>1$).
Note that this behavior will change with the statically chosen macroscopic
cross-section but improvements to this ratio will still happen as the time-
step shrinks.

##  3 HIGHER-ORDER METHODS

To further improve the GPU parallel performance, we investigate higher-order
discretization methods, particularly the robust, second-order accurate spatial
discretization method simple corner balance (SCB) [5] and the (also) robust,
second-order accurate time discretization method multiple balance (MB) [6]. In
coupling this higher temporal accuracy scheme with an iterative method that
can be efficiently solved, we hope to optimize the ratio of compute work to
communication work to better suit the numerical method for GPUs.

We start by deriving the discretized slab geometry transport equations—SCB in
space, MB in time, $S_{N}$ in angle—show that the resulting method MB-SCB is
unconditionally stable via Fourier analysis, and then derive the respective
iterative systems.

###  3.1 Derivation of MB-SCB

We begin with the time-dependent, mono-energetic, isotropic scattering slab
geometry SN transport equations with an isotropic source:

$\frac{1}{v}\frac{\partial\psi_{m}(x,t)}{\partial
t}+\mu_{m}\frac{\partial\psi_{m}(x,t)}{\partial
x}+\Sigma(x)\psi_{m}(x,t)=\frac{1}{2}\left(\Sigma_{s}(x)\sum\limits_{n=1}^{N}w_{n}\psi_{n}(x,t)+Q(x,t)\right)\;,\\\
\qquad m=1\ldots N\;,\qquad t>0\;,\qquad x\in[0,X]\;,$ (2)

where $\psi$ is the angular flux, $t$ is time, $x$ is location, $v$ is
velocity, $w_{m}$ is angular quadrature weight, $\mu_{m}$ is the angular
quadrature ordinate, $m$ is the quadrature index, and $Q$ is the isotropic
material source. The initial and boundary conditions are prescribed angular
flux distributions:

$\psi_{m}(x,0)=\psi_{m,0}(x),\qquad m=1\ldots N\;,$
$\psi_{m}(0,t)=\psi_{m,L}(t),\qquad\mu_{m}>0\;,$
$\psi_{m}(X,t)=\psi_{m,R}(t),\qquad\mu_{m}<0\;.$

We discretize these equations in time using the MB approach [6]:

$\frac{1}{v}\left(\frac{\psi_{m,k+1/2}(x)-\psi_{m,k-1/2}(x)}{\Delta
t}\right)+\mu_{m}\frac{\partial\psi_{m,k}(x)}{\partial
x}+\Sigma(x)\psi_{m,k}(x)\\\
=\frac{1}{2}\left(\Sigma_{s}(x)\sum\limits_{n=1}^{N}w_{n}\psi_{n,k}(x)+Q_{k}(x)\right)\;,$
(3) $\frac{1}{v}\frac{\psi_{m,k+1/2}(x)-\psi_{m,k}(x)}{\Delta
t/2}+\mu_{m}\frac{\partial\psi_{m,k+1/2}(x)}{\partial
x}+\Sigma(x)\psi_{m,k+1/2}(x)\\\
=\frac{1}{2}\left(\Sigma_{s}(x)\sum\limits_{n=1}^{N}w_{n}\psi_{n,k+1/2}(x)+Q_{k+1/2}(x)\right)\;,$
(4)

where $\Delta t$ is the time step size, $k$ indexes time-average quantities,
and $k\pm 1/2$ indexes time-edge quantities. Then, we discretize in space
using SCB, which involves a spatial integration over the right and left halves
of a spatial cell:

$\frac{\Delta
x_{j}}{2}\frac{1}{v}\left(\frac{\psi_{m,k+1/2,j,L}-\psi_{m,k-1/2,j,L}}{\Delta
t}\right)+\mu_{m}\left[\frac{\left(\psi_{m,k,j,L}+\psi_{m,k,j,R}\right)}{2}-\psi_{m,k,j-1/2}\right]\\\
+\frac{\Delta x_{j}}{2}\Sigma_{j}\psi_{m,k,j,L}=\frac{\Delta
x_{j}}{2}\frac{1}{2}\left(\Sigma_{s,j}\sum\limits_{n=1}^{N}w_{n}\psi_{n,k,j,L}+Q_{k,j,L}\right)\;,$
(5) $\frac{\Delta
x_{j}}{2}\frac{1}{v}\left(\frac{\psi_{m,k+1/2,j,R}-\psi_{m,k-1/2,j,R}}{\Delta
t}\right)+\mu_{m}\left[\psi_{m,k,j+1/2}-\frac{\left(\psi_{m,k,j,L}+\psi_{m,k,j,R}\right)}{2}\right]\\\
+\frac{\Delta x_{j}}{2}\Sigma_{j}\psi_{m,k,j,R}=\frac{\Delta
x_{j}}{2}\frac{1}{2}\left(\Sigma_{s,j}\sum\limits_{n=1}^{N}w_{n}\psi_{n,k,j,R}+Q_{k,j,R}\right)\;,$
(6) $\frac{\Delta
x_{j}}{2}\frac{1}{v}\left(\frac{\psi_{m,k+1/2,j,L}-\psi_{m,k,j,L}}{\Delta
t/2}\right)+\mu_{m}\left[\frac{\left(\psi_{m,k+1/2,j,L}+\psi_{m,k+1/2,j,R}\right)}{2}-\psi_{m,k+1/2,j-1/2}\right]\\\
+\frac{\Delta x_{j}}{2}\Sigma_{j}\psi_{m,k+1/2,j,L}=\frac{\Delta
x_{j}}{2}\frac{1}{2}\left(\Sigma_{s,j}\sum\limits_{n=1}^{N}w_{n}\psi_{n,k+1/2,j,L}+Q_{k+1/2,j,L}\right)\;,$
(7) $\frac{\Delta
x_{j}}{2}\frac{1}{v}\left(\frac{\psi_{m,k+1/2,j,R}-\psi_{m,k,j,R}}{\Delta
t/2}\right)+\mu_{m}\left[\psi_{m,k+1/2,j+1/2}-\frac{\left(\psi_{m,k+1/2,j,L}+\psi_{m,k+1/2,j,R}\right)}{2}\right]\\\
+\frac{\Delta x_{j}}{2}\Sigma_{j}\psi_{m,k+1/2,j,R}=\frac{\Delta
x_{j}}{2}\frac{1}{2}\left(\Sigma_{s,j}\sum\limits_{n=1}^{N}w_{n}\psi_{n,k+1/2,j,R}+Q_{k+1/2,j,R}\right)\;,$
(8)

where $\Delta x$ is the cell width, $j$ is the spatial index, $R$ is the right
hand side sub cell division, and $L$ is the left. These equations contain the
first of the two simple spatial closures—the angular flux at the cell midpoint
is a simple average of the two half-cell average quantities:

$\psi_{m,k}(x_{j})=\frac{\left(\psi_{m,k,j,L}+\psi_{m,k,j,R}\right)}{2}\;,$
(9)
$\psi_{m,k+1/2}(x_{j})=\frac{\left(\psi_{m,k+1/2,j,L}+\psi_{m,k+1/2,j,R}\right)}{2}\;.$
(10)

The second is an upstream prescription for the cell-edge angular flux:

$\psi_{m,k,j+1/2}=\begin{cases}\psi_{m,k,j,R},&\mu_{m}>0,\\\
\psi_{m,k,j+1,L},&\mu_{m}<0\;,\end{cases}$ (11)
$\psi_{m,k+1/2,j+1/2}=\begin{cases}\psi_{m,k+1/2,j,R},&\mu_{m}>0,\\\
\psi_{m,k+1/2,j+1,L},&\mu_{m}<0\;.\end{cases}$ (12)

###  3.2 Fourier Analysis of Numerical Method

To ensure that our combination of higher-order discretization schemes is still
unconditionally stable, we performed a Fourier analysis to find the slowest-
converging mode of the method. Assuming no scattering, no sources, and making
physical assumptions (only constant positive values of $\Delta x$, $\Delta t$,
$\Sigma$, and $v$) we derived SCB-MB’s eigenfunction and numerically solved it
[7].

We found that the MB-SCB scheme is unconditionally stable over $\mu\in[-1,1]$.
While experimenting with this method we did find that, under some conditions,
it can produce negative fluxes; however, the negative flux oscillations were
critically damped and dissipated with time.

###  3.3 Source Iteration (SI)

In the traditional SI, the scattering source is presumed known from a previous
iteration, which leads to the following set of equations to be solved in
transport “sweeps”. This means that new estimates of both the end of time-step
value of angular flux and time-averaged angular flux are computed together in
each cell.

For SCB in slab geometry, this means there is a local $4\times 4$ matrix to be
solved in each cell. For $\mu_{m}>0$, this equation has the form:

$\mathbf{A}^{+}\begin{bmatrix}\psi_{m,k,j,L}\\\ \psi_{m,k,j,R}\\\
\psi_{m,k+1/2,j,L}\\\ \psi_{m,k+1/2,j,R}\end{bmatrix}=\mathbf{b}^{+}\;.$ (13)

Here, the $\mathbf{A}^{+}$ matrix has the following structure and element
definitions:

$\mathbf{A}^{+}=\begin{bmatrix}\frac{\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}&\frac{\mu_{m}}{2}&\frac{\Delta x_{j}}{2v\Delta t}&0\\\
-\frac{\mu_{m}}{2}&\frac{\mu_{m}+\Delta x_{j}\Sigma_{j}}{2}&0&\frac{\Delta
x_{j}}{2v\Delta t}\\\ -\frac{\Delta x_{j}}{v\Delta t}&0&\frac{\Delta
x_{j}}{v\Delta t}+\frac{\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}&\frac{\mu_{m}}{2}\\\ 0&-\frac{\Delta x_{j}}{v\Delta
t}&-\frac{\mu_{m}}{2}&\frac{\Delta x_{j}}{v\Delta t}+\frac{\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}\\\ \end{bmatrix}$ (14)

and $\mathbf{b}^{+}$ is given by

$\mathbf{b}^{+}=\begin{bmatrix}\frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k,j,L}^{(l)}+Q_{k,j,L}\right)+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,L}+\mu_{m}\psi_{m,k,j-1,R}\\\ \frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k,j,R}^{(l)}+Q_{k,j,R}\right)+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,R}\\\ \frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k+1/2,j,L}^{(l)}+Q_{k+1/2,j,L}\right)+\mu_{m}\psi_{m,k+1/2,j-1,R}\\\
\frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k+1/2,j,R}^{(l)}+Q_{k+1/2,j,R}\right)\end{bmatrix}\;.$
(15)

For $\mu_{m}<0$, this equation has the form:

$\mathbf{A}^{-}\begin{bmatrix}\psi_{m,k,j,L}\\\ \psi_{m,k,j,R}\\\
\psi_{m,k+1/2,j,L}\\\ \psi_{m,k+1/2,j,R}\end{bmatrix}=\mathbf{b}^{-}\;.$ (16)

Here, the $\mathbf{A}^{-}$ matrix has the following structure and element
definitions:

$\mathbf{A}^{-}=\begin{bmatrix}\frac{-\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}&\frac{\mu_{m}}{2}&\frac{\Delta x_{j}}{2v\Delta t}&0\\\
-\frac{\mu_{m}}{2}&\frac{-\mu_{m}+\Delta x_{j}\Sigma_{j}}{2}&0&\frac{\Delta
x_{j}}{2v\Delta t}\\\ -\frac{\Delta x_{j}}{v\Delta t}&0&\frac{\Delta
x_{j}}{v\Delta t}+\frac{-\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}&\frac{\mu_{m}}{2}\\\ 0&-\frac{\Delta x_{j}}{v\Delta
t}&-\frac{\mu_{m}}{2}&\frac{\Delta x_{j}}{v\Delta t}+\frac{-\mu_{m}+\Delta
x_{j}\Sigma_{j}}{2}\\\ \end{bmatrix}$ (17)

and $\mathbf{b}^{-}$ is given by

$\mathbf{b}^{-}=\begin{bmatrix}\frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k,j,L}^{(l)}+Q_{k,j,L}\right)+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,L}\\\ \frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k,j,R}^{(l)}+Q_{k,j,R}\right)+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,R}-\mu_{m}\psi_{m,k,j+1,L}\\\ \frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k+1/2,j,L}^{(l)}+Q_{k+1/2,j,L}\right)\\\
\frac{\Delta
x_{j}}{4}\left(\Sigma_{s,j}\phi_{k+1/2,j,R}^{(l)}+Q_{k+1/2,j,R}\right)-\mu_{m}\psi_{m,k+1/2,j+1,L}\end{bmatrix}\;.$
(18)

After sweeping the mesh cells in the appropriate directions for each angle in
the quadrature set, the scalar flux vector can be updated via

$\begin{bmatrix}\phi_{k,j,L}\\\ \phi_{k,j,R}\\\ \phi_{k+1/2,j,L}\\\
\phi_{k+1/2,j,R}\end{bmatrix}=\sum\limits_{n=1}^{N}w_{n}\begin{bmatrix}\psi_{n,k,j,L}\\\
\psi_{n,k,j,R}\\\ \psi_{n,k+1/2,j,L}\\\ \psi_{n,k+1/2,j,R}\end{bmatrix}\;,$
(19)

and the source iteration can continue until this scalar flux vector ceases
changing between iterations. After converging, the simulation can move to the
next time step.

###  3.4 One-Cell Inversion (OCI)

In OCI, the scattering source is subtracted to the left-hand side, and the
information that comes from cells other than cell $j$ is assumed to be known
from a previous iteration. This means that all $4N$ angular fluxes ($N$
angles, $L$ and $R$, $k$ and $k+1/2$) are computed simultaneously in cell $j$.
For SCB in slab geometry, this means there is a local $4N\times 4N$ matrix to
be solved in each cell.

$\left(\mathbf{A}-\mathbf{S}\right)\mathbf{\Psi}=\mathbf{c}\;,$ (20)

where,

$\mathbf{A}=\left[{\begin{array}[]{cccc}\mathbf{A}_{1}&0&\cdots&0\\\
0&\mathbf{A}_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\mathbf{A}_{N}\\\ \end{array}}\right]\;.$ (21)

Here, the $\mathbf{A}_{m}$ matrix has the following structure and element
definitions:

$\mathbf{A}_{m}=\begin{cases}\mathbf{A}^{+}&\mu_{m}>0\\\
\mathbf{A}^{-}&\mu_{m}<0\\\ \end{cases}\;.$ (22)

The scattering source $\mathbf{S}$ is defined by

$S_{l,n}=\frac{\Delta x_{j}\Sigma_{s,j}}{4}w_{n}\;,$ (23)

$\mathbf{\Psi}$ is given by:

$\mathbf{\Psi}=\begin{bmatrix}\mathbf{\psi}_{1}\;\mathbf{\psi}_{2}\;\cdots\;\mathbf{\psi}_{N}\end{bmatrix}^{T}\;.$
(24)

where

$\mathbf{\psi}_{n}=\begin{bmatrix}\psi_{n,k,j,L}\;\psi_{n,k,j,R}\;\psi_{n,k+1/2,j,L}\;\psi_{n,k+1/2,j,R}\end{bmatrix}^{T}\;,$
(25)

and $\mathbf{c}$ is given by

$\mathbf{c}=\left[\mathbf{c}_{1}\;\mathbf{c}_{2}\;\cdots\;\mathbf{c}_{N}\right]^{T}\;,$
(26)

where

$\mathbf{c}_{m}=\begin{cases}\mathbf{c}^{+}&\mu_{m}>0\\\
\mathbf{c}^{-}&\mu_{m}<0\\\ \end{cases}\;,$ (27)
$\mathbf{c}^{+}=\begin{bmatrix}\frac{\Delta x_{j}}{4}Q_{k,j,L}+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,L}+\mu_{m}\psi^{(l)}_{m,k,j-1,R}\\\
\frac{\Delta x_{j}}{4}Q_{k,j,R}+\frac{\Delta x_{j}}{2v\Delta
t}\psi_{m,k-1/2,j,R}\\\ \frac{\Delta
x_{j}}{4}Q_{k+1/2,j,L}+\mu_{m}\psi^{(l)}_{m,k+1/2,j-1,R}\\\ \frac{\Delta
x_{j}}{4}Q_{k+1/2,j,R}\end{bmatrix}\;,$ (28)
$\mathbf{c}^{-}=\begin{bmatrix}\frac{\Delta x_{j}}{4}Q_{k,j,L}+\frac{\Delta
x_{j}}{2v\Delta t}\psi_{m,k-1/2,j,L}\\\ \frac{\Delta
x_{j}}{4}Q_{k,j,R}+\frac{\Delta x_{j}}{2v\Delta
t}\psi_{m,k-1/2,j,R}-\mu_{m}\psi^{(l)}_{m,k,j+1,L}\\\ \frac{\Delta
x_{j}}{4}Q_{k+1/2,j,L}\\\ \frac{\Delta
x_{j}}{4}Q_{k+1/2,j,R}-\mu_{m}\psi^{(l)}_{m,k+1/2,j+1,L}\end{bmatrix}\;.$ (29)

Where $l$ indexes solutions from the previous iteration. One cell inversion
iterations continue until this angular flux vector ceases changing between
iterations. After convergence, the time-step counter is incremented and the
within time-step process can be repeated.

##  4 GPU IMPLEMENTATIONS

Modern exascale-class supercomputers (e.g., El Capitan, Frontier, Aurora)
currently being deployed use heterogeneous architectures—with both CPUs and
accelerators like GPUs—to reach exascale performance, motivating the
exploration of the performance of these numerical and iterative schemes on the
GPU.

When we first implemented our MB-SCB-OCI scheme on a GPU using discrete dense
matrix solves in every cell we found that wall-clock performance was
drastically worse then even pure Python implementations on a single CPU core.
We profiled on our naive GPU kernels and found that due to a) non-optimum data
transfers between the GPU and host, b) naive scheduling of our compute solves
in every cell, and c) overhead in calling an individual dense solver for every
cell, performance was lackluster.

Therefore we implement our neutron transport deterministic methods in a sparse
linear algebra framework—representing a problem in all cells in all angles as
a single sparse linear algebra problem—when we solve on GPUs. This whole-
problem sparse linear algebra method will produce the same solutions, within
the iteration convergence criteria, as implementations that manually conduct a
sweep or a cell inversion for source iteration or one-cell inversion,
respectively. The underlying discretizations are the same.

A sparse linear algebra framework has the benefits of using off the shelf
libraries which are often written and optimized by the manufacturers of the
hardware themselves [8]. Further benefits of utilizing pre-written optimized
APIs via whole-problem sparse linear algebra are:

  1. 1.

Generalizes solve schemes to operations that can run using almost any linear
algebra API in many languages (e.g., PETSc [9], CuPy [10], Trilinos [11],
Scipy [12]) on almost any most back-end (e.g., x86/ARM CPUs, AMD/Nvidia/Intel
GPUs);

  2. 2.

Enables rapid experimentation or even run-time alteration of matrix solve
schemes (e.g., Gauss-Sidel, Krylov solvers);

  3. 3.

Can take advantage of off-the-shelf preconditioner libraries; and

  4. 4.

Enables easy implementation of performance improvements to the mathematical
methods as hardware and libraries are further optimized or developed.

Historically this whole-problem sparse matrix approach has not been commonly
implemented, especially for source iteration where sweeps can be optimized to
decrease memory footprint within a given architecture. In forming the matrices
themselves, we are introducing a bottleneck on the CPU, since all our
implementations form both the whole-problem versions of the coefficient and
right-hand-side matrices on the host then send them to the GPU. As we increase
the problem size or increase dimensions, this bottleneck will grow and
potentially become prohibitively large for the limited amount of memory on a
GPU.

However, one of the benefits of using one-cell inversion is that the problem
can be dynamically decomposed to split over multiple accelerators within an
iteration. While we believe this whole-problem approach is appropriate for
these initial explorations in comparing and contrasting iterative schemes on a
GPU, specifically written GPU kernels in device code that do not form the
whole-problem matrices and that are written for a given hardware accelerator
will likely outperform our implementation after sufficient optimization.

##  5 PERFORMANCE ANALYSIS & RESULTS

We developed a Python-based implementation of the MB-SCB numerical method
using both SI and OCI iterative schemes [7]. We examine performance on multi-
core CPUs (via Numba [13] and SciPy [12]) using the manual parallelization
scheme that an iterative method might allow (sweeps within an angle in all
cells for SI and parallel dense matrix solutions within a single cell in all
angles in OCI), and on GPUs using a whole-problem sparse matrix
implementation. All dense matrices were solved with LAPACK’s _gesv (LU
decomposition with partial pivoting and row interchanges) via SciPy, and all
sparse linear algebra systems were solved with cuSparse’s implementation of
GMRES via CuPy.

To initially verify our discretizations and iterative methods on all hardware
back-ends, we used a transient, source-free, pure absorber problem with an
analytical solution, the AZURV1 problem [14], and transient Monte Carlo
solutions from MC/DC [15]. Once confident in the correctness of our solutions,
we conducted run-time testing on a transient version of Reed’s Problem and a
transient homogeneous slab with varying mean free path, scattering ratio, time
step, and quadrature order.

This performance analysis was conducted on an Intel Core i7-10875H CPU with 16
logical cores at 3.4 GHz and an Nvidia RTX 2060 (mobile). The point of these
simulations is to compare the wall clock run-time behavior of these iterative
methods, not to provide a realistic reactor benchmark.

###  5.1 Transient Reed’s Problem

Figure 3 shows a time-dependent form of Reed’s Problem [16] that serves as an
initial performance benchmark and to further verify our implementations. An
initial condition is set at $t=$ $0\text{\,}\mathrm{s}$ with 0 angular flux
everywhere. Then, at $t>$$0\text{\,}\mathrm{s}$, the conditions shown in
Figure 3 are applied with $v=$
$4\text{\,}\mathrm{cm}\text{\,}{\mathrm{s}}^{-1}$, $\Delta t=$
$0.1\text{\,}\mathrm{s}$, $N_{\text{time}}=$ $5\text{\,}\mathrm{,}$ $Q$ is in
$\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{cm}}^{-3}$, $\Sigma$ and
$\Sigma_{s}$ are in $\text{\,}{\mathrm{cm}}^{-1}$, and $\Delta x$ is in
$\text{\,}\mathrm{cm}$ in $S_{64}$. Figure 4 at right shows the transient
solution through time as it approaches the known steady-state solution.

Figure 3: Multi-region transient Reed’s problem Figure 4: Solution at various
times as it approaches steady state ($\boldsymbol{t}$ is in (s))

Table 1 shows run-times of the two iterative schemes on both hardware backends
from the fifth and final time-step. These run-times include every operation in
one time step: building and moving the coefficient matrix once per time-step,
building and moving the right-hand-side vector once per iteration within a
time-step, calling the matrix solver, and resorting the raw matrix solution
vector back into our standard form if required.

On a CPU, OCI takes one order of magnitude longer in wall clock time than SI
(both threaded on 16 CPU cores). This inverts when going to the GPU with OCI
pulling an almost factor of two speedup over SI. However, SI threaded on 16
CPU cores requires the least wall clock time to complete this time-step.

Table 1 also shows how readily one-cell inversion can be implemented in a
sparse linear algebra framework. It seems that allowing the GPU purpose-built
libraries to decide when and how to solve the presented parallel block Jacobi
system is superior to even a 16 thread manual implementation on CPU, given the
2.4$\times$ speedup. As building and moving data to and from the GPU is
further optimized, we would expect the relative performance of one-cell
inversion to continue to be enhanced.

Table 1: Wall clock run-times of one time-step of the transient version of Reed’s problem. Iteration Scheme | Back-end | Wall Time (s)
---|---|---
OCI | GPU | 17.01
SI | GPU | 33.76
OCI | CPU | 41.5
SI | CPU | 4.85

###  5.2 GPU Explorations Across Quadrature Sets

As convergence—and thus wall clock time—depends on spatial cell optical
thickness, time step size, scattering ratio, and the number of angles in
quadrature, we performed a series of run-time tests to compare the GPU methods
directly to one another. We used a slab thickness of $L=$
$10\text{\,}\mathrm{cm}$, $\Sigma=$ $2\text{\,}{\mathrm{cm}}^{-1}$, with five
time steps ($N_{time}=$ $5\text{\,}$), an isotropic material source of $Q=$
$0.1\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{cm}}^{-3}$, an initial
condition of zero angular flux everywhere, and vacuum boundary conditions on
both the left and right sides. Figure 5 shows how performance varies with
scattering cross section ($c=\Sigma_{s}/\Sigma$), cell optical thickness
($\Sigma\;\Delta x$), and time-step size ($\Delta t$), measured via speedup
($t_{\text{wall clock}}^{\text{SI}}/t_{\text{wall clock}}^{\text{OCI}}$)
across a wide range of angular quadrature orders.

Figure 5: Wall clock speedup of OCI over SI across different problems and
quadrature orders.

In lower order quadratures, OCI implementation on GPU provides a significant
speedup with a maximum of 3.4$\times$ over SI. $S_{64}$ shows the best
performance, with runs regularly reaching over a factor of three speedup. As
we further increase the quadrature order, the speedup of OCI over SI
decreases, with no observed speedup at $S_{256}$. This points to our original
idea that parallelization over the most dominate domain will incur the
greatest performance boost. For this scheme when N (quadrature order) is
smaller than J (number of cells), OCI leads in performance. Then as N
approaches J, the benefit of OCI decreases.

Closer inspection using profilers shows that when the quadrature order is
increased by a factor of four, the SI run-time of the GPU-based GMRES matrix
solver increases linearly (i.e., four times the angles, four times the solve
time). However, in OCI implementation, runtime increases by a factor of seven.
This suggests that there is an optimal matrix size for OCI/GMRES that can be
tuned, as OCI allows for easy dynamic domain decomposition at any arbitrary
location in the system within an iteration.

##  6 CONCLUSIONS & FUTURE WORK

OCI solved as a whole matrix problem on GPUs provides significant speedup in
quadrature sets at or below $S_{128}$ over SI in the MB-SCB scheme. A whole-
matrix sparse linear algebra framework on GPUs has benefits that merit further
exploration when applied to transient deterministic neutron transport methods,
since it exploits access to the computational power locked in accelerators
found on modern exascale computing systems.

While both OCI and the whole-problem sparse linear algebra framework do have
significant downsides (primarily in convergence rate at higher ordinates and
memory footprint, respectively), we expect to continue to see performance
speedups as the implementation is further optimized. In the future, we plan to
implement this method in two dimensions using more advanced linear algebra
libraries (e.g., PETSc) in a multi-node enabled C++ framework. There we can
further compare and contrast the whole-problem implementation versus a
standard algorithm, as well as explore any memory footprint issues and begin
to make direct comparisons to production codes.

## ACKNOWLEDGEMENTS

This work was supported by the Center for Exascale Monte-Carlo Neutron
Transport (CEMeNT), a PSAAP-III project funded by the Department of Energy,
grant number DE-NA003967.

## References

  * [1] M. L. Adams and E. W. Larsen.  “Fast iterative methods for discrete-ordinates particle transport calculations.”  _Progress in nuclear energy_ , volume 40(1), pp. 3–159 (2002).
  * [2] R. S. Baker.  “An SN Algorithm for Modern Architectures.”  _Nuclear science and engineering_ , volume 185(1), pp. 107–116 (2017).
  * [3] D. Anistratov and Y. Y. Azmy.  “Iterative stability analysis of spatial domain decomposition based on block Jacobi algorithm for the diamond-difference scheme.”  _Journal of Computational Physics_ , volume 297(C), pp. 462–479 (2015).
  * [4] D. S. Hoagland and Y. Y. Azmy.  “Hybrid approaches for accelerated convergence of block-Jacobi iterative methods for solution of the neutron transport equation.”  _Journal of Computational Physics_ , volume 439 (2021).  URL https://www.osti.gov/biblio/1815296.
  * [5] M. L. Adams.  “Subcell balance methods for radiative transfer on arbitrary grids.”  _Transport Theory and Statistical Physics_ , volume 26(4-5), pp. 385–431 (1997).  URL https://doi.org/10.1080/00411459708017924.
  * [6] I. Variansyah, E. W. Larsen, and W. Martin.  “A Robust Second-Order Multiple Balance Method for Time-Dependent Neutron Transport Simulations.”  _The European Physical Journal Conferences_ , volume 247, p. 03024 (2021).
  * [7] J. P. Morgan and I. Variansyah.  “Therefore.” (2023).  URL https://doi.org/10.5281/zenodo.7806833.
  * [8] “Nvidia cuSPARSE.” (2022).  URL https://docs.nvidia.com/cuda/cusparse/.
  * [9] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, and et. al.  “PETSc Web page.”  https://petsc.org/ (2023).  URL https://petsc.org/.
  * [10] R. Okuta, Y. Unno, D. Nishino, S. Hido, and C. Loomis.  “CuPy: A NumPy-Compatible Library for NVIDIA GPU Calculations.”  In _Proceedings of Workshop on Machine Learning Systems (LearningSys) in The Thirty-first Annual Conference on Neural Information Processing Systems (NIPS)_ (2017).  URL http://learningsys.org/nips17/assets/papers/paper_16.pdf.
  * [11] T. Trilinos Project Team.  “The Trilinos Project Website.” (2020 (acccessed May 22, 2020)).  URL https://trilinos.github.io.
  * [12] P. Virtanen, R. Gommers, Oliphant, SciPy 1.0 Contributors, and et. al.  “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.”  _Nature Methods_ , volume 17, pp. 261–272 (2020).
  * [13] S. K. Lam, A. Pitrou, and S. Seibert.  “Numba: A LLVM-Based Python JIT Compiler.”  In _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_ , LLVM ’15. Association for Computing Machinery, New York, NY, USA (2015).  URL https://doi.org/10.1145/2833157.2833162.
  * [14] B. Ganapol, R. Baker, J. Dahl, and R. E. Alcouffe.  “Homogeneous Infinite Media Time-Dependent Analytical Benchmarks.”  In _International Meeting on Mathematical Methods for Nuclear Applications_ , volume 41(25). Salt Lake City, UT (2001).
  * [15] I. Variansyah, J. P. Morgan, K. E. Niemeyer, and R. G. McClarren.  “Development of MC/DC: a performant, scalable, and portable Python-based Monte Carlo neutron transport code.”  In _International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering_. Niagara Falls, Ontario, Canada (2023).
  * [16] J. Warsa.  “Analytical SN solutions in heterogeneous slabs using symbolic algebra computer programs.”  _Annals of Nuclear Energy_ , volume 29, pp. 851–874 (2002).