Acceleration of the Sn Equations with Highly Anisotropic Scattering using the Fokker-Planck Equation

TL;DR

This paper proposes using the Fokker-Planck equation as a preconditioner to accelerate the convergence of discrete ordinates methods in highly anisotropic scattering transport problems, where traditional acceleration techniques fail.

Contribution

It introduces a novel approach of employing the Fokker-Planck equation as a preconditioner to improve convergence in forward-peaked transport simulations.

Findings

Fokker-Planck preconditioning significantly speeds up convergence.

Standard acceleration methods are ineffective for highly anisotropic scattering.

The approach enhances computational efficiency in forward-peaked transport modeling.

Abstract

The discrete ordinates method can model forward-peaked transport problems accurately. However, convergence of discrete ordinates solution can become arbitrarily slow upon use of standard iterative procedures like source iteration and GMRES. Standard zeroth and first moment-based acceleration methods like nonlinear diffusion acceleration and diffusion synthetic acceleration are ineffective in accelerating such problems because these methods do not correct higher order Legendre-moments of angular flux. We explore the idea of using Fokker-Planck as a preconditioner to accelerate forward-peaked transport problems in this paper.

Tables

Table 1. Table 1: Comparison of Numerical and Theoretical Spectral Radii

Kernel/Parameter	${\rho }_{\mathrm{FPSA}}^{\mathrm{MPD}}$-FA	${\rho }_{\mathrm{FPSA}}^{\mathrm{MPD}}$-Measured	${\rho }_{\mathrm{FPSA}}^{\mathrm{WFD}}$-FA	${\rho }_{\mathrm{FPSA}}^{\mathrm{WFD}}$-Measured
SRK/$\eta =2.83\times {10}^{-5}$	0.4706	0.4706	0.2121	0.2120
EK/$\Delta ={10}^{-5}$	0.1932	0.1954	0.6246	0.6327
HGK/$g=0.9$	0.4304	0.4303	0.4177	0.4177

Table 2. Table 2: SRK - Number of Iterations

BC/Source	Restart	GMRES Iteration Count	SI Iteration Count
Vacuum/Distributed			$>$ 150000
	50	3305
	100	2445
	150	1875
	200	1540
Beam/Zero			$>$ 150000
	50	2602
	100	2200
	150	1895
	200	1735

Table 3. Table 3: SRK - Solver Runtime [s]

BC/Source	Restart	GMRES Runtime	SI Runtime
Vacuum/Distributed			$>$ 3000
	50	64.97
	100	50.41
	150	37.76
	200	32.71
Beam/Zero			$>$ 3000
	50	50.68
	100	43.99
	150	41.29
	200	36.11

Table 4. Table 4: SRK - Vacuum Boundaries/Unit Distributed Source - Number of Iterations

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	1487	9	6	14	10
	15/32	1499	9	7	14	12
Factorize	15/16		9	7	14	10
	15/32		9	8	14	12

Table 5. Table 5: SRK - Vacuum Boundaries/Unit Distributed Source - Runtime [s]

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	28.75	8.79	6.15	12.01	5.98
	15/32	28.12	36.76	18.36	54.13	19.77
Factorize	15/16		1.62	2.45	0.3373	0.2501
	15/32		2.75	4.73	0.4244	0.3392

Table 6. Table 6: SRK - Beam Source - Number of Iterations

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	1357	12	8	21	13
	15/32	1335	13	10	23	19
Factorize	15/16		12	9	21	13
	15/32		13	11	23	19

Table 7. Table 7: SRK - Beam Source - Runtime [s]

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	26.54	11.39	6.975	9.69	5.403
	15/32	27	42.59	20.50	58.3	20.26
Factorize	15/16		1.687	2.547	0.4475	0.3074
	15/32		2.927	5.061	0.611	0.4828

Table 8. Table 8: EK - Vacuum Boundaries/Unit Distributed Source - Number of Iterations

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	2217	7	8	9	15
	15/32	2256	10	9	17	19
Factorize	15/16		7	9	9	15
	15/32		10	10	17	19

Table 9. Table 9: EK - Vacuum Boundaries/Unit Distributed Source - Runtime [s]

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	51.73	10.32	13.34	11.91	18.49
	15/32	56.10	31.14	25.06	42.82	27.85
Factorize	15/16		2.3754	4.5586	0.3762	0.454
	15/32		4.155	8.463	0.6098	0.6424

Table 10. Table 10: EK - Beam Source - Number of Iterations

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	2086	9	10	12	35
	15/32	1932	14	12	24	28
Factorize	15/16		14	13	12	35
	15/32		14	13	24	28

Table 11. Table 11: EK - Beam Source - Runtime [s]

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	39.55	9.45	11.63	8.335	19.42
	15/32	36.14	24.89	18.57	34.74	24.87
Factorize	15/16		2.842	6.657	0.2929	0.6853
	15/32		2.898	6.799	0.6329	0.6585

Table 12. Table 12: HGK - Vacuum Boundaries/Unit Distributed Source - Number of Iterations

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	1150	10	7	28	18
	15/32	1461	14	12	32	31
Factorize	15/16		10	8	28	18
	15/32		14	13	32	31

Table 13. Table 13: HGK - Vacuum Boundaries/Unit Distributed Source - Runtime [s]

Invert FP	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	83.66	288.1	160.9	931.9	373.4
	15/32	102.5	1055	572.7	2618	1428
Factorize	15/16		6.390	12.27	2.193	1.385
	15/32		11.75	28.820	2.651	2.5778

Table 14. Table 14: HGK - Beam Source - Number of Iterations

Invert FP	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	597	16	12	29	17
	15/32	1634	21	19	32	31
Factorize	15/16		12	9	29	17
	15/32		17	16	32	31

Table 15. Table 15: HGK - Beam Source - Runtime [s]

FP-Solve	L/N	$\mathrm{GMRES}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{GMRES}}^{\mathrm{WFD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{MPD}}$	${\mathrm{FPSA}}_{\mathrm{SI}}^{\mathrm{WFD}}$
GMRES	15/16	42.29	479.9	331.2	1040	285.7
	15/32	115.1	1804	990.6	2408	1315
Factorize	15/16		6.795	12.41	2.131	1.316
	15/32		12.59	29.69	2.945	2.452

Equations

$$\mu \frac{\partial }{\partial z}\psi (z,\mu )+{\sigma }_t(z)\psi (z,\mu )={\int }_{-1}^{1}d{\mu }^{\prime }{\sigma }_s(z,{\mu }_0)\psi (z,{\mu }^{\prime })+q(z,\mu ),$$

$$\psi (z_0,\mu )={\psi }_L(\mu )\text{for μ>0}\text{and}\psi (z_1,\mu )={\psi }_R(\mu )\text{for μ<0}.$$

$${\sigma }_{s,l}(z)={\int }_{-1}^{1}d{\mu }_0{\sigma }_s(z,{\mu }_0)P_l({\mu }_0).$$

$$\mu \frac{\partial }{\partial z}\psi (z,\mu )+{\sigma }_t(z)\psi (z,\mu )=\sum _{l=0}^L\frac{2l+1}{2}{\sigma }_{s,l}(z)P_l(\mu ){\varphi }_l(z)+q(z,\mu ),$$

$${\varphi }_l(z)={\int }_{-1}^{1}d\mu P_l(\mu )\psi (z,\mu ).$$

$$\mu \frac{\partial \psi }{\partial z}+{\sigma }_a(z)\psi (z,\mu )-\frac{{\sigma }_{tr}(z)}{2}{\mathcal{L}}_{\text{FP}}\psi (z,\mu )=q(z,\mu ),$$

$${\mathcal{L}}_{\text{FP}}=\frac{\partial }{\partial \mu }(1-{\mu }^2)\frac{\partial }{\partial \mu }$$

$$L\psi (z,\mu )=S\psi (z,\mu )+q(z,\mu ),$$

$$L=\mu \frac{\partial }{\partial z}+{\sigma }_t(z)\text{and}S=\sum _{l=0}^L\frac{2l+1}{2}{\sigma }_{s,l}(z)P_l(\mu ){\int }_{-1}^{1}d\mu P_l(\mu ).$$

$$L{\psi }^{m+1}(z,\mu )=S{\psi }^m(z,\mu )+q(z,\mu ),$$

$$L{\psi }^{m+\frac{1}{2}}(z,\mu )=S{\psi }^m(z,\mu )+q(z,\mu ),$$

$${\psi }^{m+1}(z,\mu )={\psi }^{m+\frac{1}{2}}(z,\mu )+F^{-1}S({\psi }^{m+\frac{1}{2}}(z,\mu )-{\psi }^m(z,\mu ))$$

$$||{\varphi }_l^{m+1}(z)-{\varphi }_l^m(z)|{|}_{\infty }>tolerence$$

$$F=\mu \frac{\partial }{\partial z}+{\sigma }_a(z)-\frac{{\sigma }_{tr}(z)}{2}{\mathcal{L}}_{\text{FP}}$$

$$\mu \frac{\partial {\psi }_n}{\partial z}+{\sigma }_a(z){\psi }_n(z)-{\nabla }_n^2\psi (z)=q_n(z),$$

$${\nabla }_n^2\psi (z)={\gamma }_{n+1/2}{\dot{\psi }}_{n+1/2}(z)-{\gamma }_{n-1/2}{\dot{\psi }}_{n-1/2}(z),$$

$${\dot{\psi }}_{n+1/2}(z)=\frac{{\psi }_{n+1}(z)-{\psi }_n(z)}{{\mu }_{n+1}-{\mu }_n},$$

$${\gamma }_{n+1/2}={\gamma }_{n-1/2}+\nu {\mu }_n{w}_n,\text{with}{\gamma }_{n-1/2}=0.$$

$${\mathcal{L}}_{\text{FP}}{P}_l(\mu )=-l(l+1)P_l(\mu ).$$

$${\int }_{-1}^1\left[{\mathcal{L}}_{\text{FP}}{P}_l(\mu )\right]f(\mu )d\mu ={\int }_{-1}^1{P}_l(\mu )\left[{\mathcal{L}}_{\text{FP}}f(\mu )\right]d\mu .$$

$${\int }_{-1}^1{P}_l(\mu )\left[{\mathcal{L}}_{\text{FP}}f(\mu )\right]d\mu =-l(l+1){\int }_{-1}^1{P}_l(\mu )f(\mu )d\mu .$$

$$\begin{gathered} \sum _{n=1}^N{w}_n{P}_l({\mu }_n){\nabla }_n^2\psi (z) \\ =-l(l+1)\sum _{n=1}^N{w}_n{P}_l({\mu }_n){\psi }_n(z), \end{gathered}$$

$${\nabla }^2\Psi (z)=\mathbf{F}\Psi (z)$$

$$\mathbf{F}={\mathbf{V}}^{-1}\mathbf{LV},$$

$${\mathbf{V}}_{i,j}$$

$${\mathbf{L}}_{i,i}$$

$$\mathbf{H}\frac{\partial }{\partial z}\Psi (z)+{\sigma }_a(z)\Psi (z)-\frac{{\sigma }_{tr}(z)}{2}\mathbf{F}\Psi (z)=Q(z),$$

$$\mathbf{H}={\mathrm{diag}}_{n=1,N}({\mu }_n),$$

$${\nabla }^2\Psi (z)=\mathbf{W}\Psi (z)$$

$${\mathbf{W}}_{i,j}=\{\begin{array}{l}\begin{array}{rlrl}& \frac{1}{w_n}\left(\frac{{\gamma }_{n-1/2}}{{\mu }_n-{\mu }_{n-1}}\right),& j& =i-1,j>1,\\ & \frac{1}{w_n}\left(\frac{{\gamma }_{n+1/2}}{{\mu }_{n+1}-{\mu }_n}\right),& j& =i+1,j<N,\\ & -\frac{1}{w_n}\left(\frac{{\gamma }_{n-1/2}}{{\mu }_n-{\mu }_{n-1}}+\frac{{\gamma }_{n+1/2}}{{\mu }_{n+1}-{\mu }_n}\right),& j& =i,\end{array}\end{array}$$

$$\mathbf{H}\frac{\partial }{\partial z}\Psi (z)+{\sigma }_a(z)\Psi (z)-\frac{{\sigma }_{tr}(z)}{2}\mathbf{W}\Psi (z)=Q(z),$$

$${\psi }^{m+1}-\psi ={\psi }^{m+\frac{1}{2}}-\psi +F^{-1}S({\psi }^{m+\frac{1}{2}}-\psi +\psi -{\psi }^m).$$

$$\psi -{\psi }^{m+\frac{1}{2}}={ϵ}^{m+\frac{1}{2}},\psi -{\psi }^{m+\frac{1}{2}}={ϵ}^{m+\frac{1}{2}},\mathrm{and}\psi -{\psi }^m={ϵ}^m,$$

$${ϵ}^{m+1}={ϵ}^{m+\frac{1}{2}}-F^{-1}S({ϵ}^m-{ϵ}^{m+\frac{1}{2}}).$$

$${ϵ}_l^m=\underset{-1}{\overset{1}{\int }}d\mu P_l(\mu ){ϵ}^m=\sum _{n=1}^N{w}_n{P}_l({\mu }_n){ϵ}^m.$$

Full text

Acceleration of the $\text{S}_{N}$ Equations with Highly Anisotropic Scattering using the Fokker-Planck Equation

$\mathrm{Japan~{}K.~{}Patel,}^{1}~{}\mathrm{James~{}S.~{}Warsa,}^{2}~{}\mathrm{and~{}Anil~{}K.~{}Prinja}^{3}$

${}^{1}~{}\mathrm{Department~{}of~{}Mechanical~{}and~{}Aerospace~{}Engineering,~{}The~{}Ohio~{}State~{}University,~{}Columbus,~{}OH}$

${}^{2}~{}\mathrm{Transport~{}Methods~{}Group~{}(CCS-2),~{}Los~{}Alamos~{}National~{}Laboratory,~{}Los~{}Alalmos,~{}NM}$

${}^{3}~{}\mathrm{Department~{}of~{}Nuclear~{}Engineering,~{}The~{}University~{}of~{}New~{}Mexico,~{}Albuquerque,~{}NM}$

patel.3545@osu.edu, warsa@lanl.gov, prinja@unm.edu

Abstract: The discrete ordinates method can model forward-peaked transport problems accurately. However, convergence of discrete ordinates solution can become arbitrarily slow upon use of standard iterative procedures like source iteration and GMRES. Standard zeroth and first moment-based acceleration methods like nonlinear diffusion acceleration and diffusion synthetic acceleration are ineffective in accelerating such problems because these methods do not correct higher order Legendre-moments of angular flux. We explore the idea of using Fokker-Planck as a preconditioner to accelerate forward-peaked transport problems in this paper.

Keywords: radiation transport, highly anisotropic scattering, synthetic acceleration, Fokker-Planck

1 Introduction

Transport problems with forward-peaked scattering kernels are encountered in several applications related to plasma physics, radiation sheilding, medical physics, and astrophysics. Such problems have extremely small mean free paths and nearly singular differential scattering cross-sections in the forward direction. Use of discrete ordinates method with standard methods like source iteration and GMRES can become extremely inefficient due to these properties. Standard acceleration techniques like diffusion synthetic acceleration (DSA) (Alcouffe, 1977) and nonlinear diffusion acceleration (NDA) (Smith et al., 2011) are ineffective in accelerating such problems because they assume all moments higher than the zeroth moment are inconsequential to the convergence of the solution. We see, later in this paper, that such an assumption becomes invalid for forward-peaked problems.

Several innovations have been made to accelerate the convergence of such problems. Valougeorgis, Williams, and Larsen (Valougeorgis et al., 1988) presented their work on stability analysis of $\text{P}_{L}$ acceleration applied to anisotropic neutron transport problems. This paper presented an extremely valuable framework for theoretical development and testing of future acceleration methods for transport problems with anisotropic scattering. Khattab and Larsen presented their modified $\text{P}_{L}$ acceleration method that used a modified form of $\text{P}_{L}$ equations with modified scattering cross section moments in (Khattab and Larsen, 1991). Morel and Manteuffel presented their angular multigrid method for solution of problems with high anisotropy in (Morel and Manteuffel, 1991). The angular multigrid method proved to be effective in 1D with a maximum spectral radius of $0.6$. The method however had to be modified later to preserve stability for problems in higher spatial dimensions in (Pautz et al., 1998). Turcksin and Morel integrated diffusion synthetic acceleration and angular multigrid to develop their diffusion synthetic acceleration-angular multigrid method in (Trucksin, 2012).

Several approximations to the transport equation have also been derived to tackle forward-peaked problems. Most prominent of these is the Fokker-Planck approximation which is an asymptotic limit of the Boltzmann equation (Pomraning, 1992) as scattering angle and energy loss become diminishingly small (Fokker-Planck limit). Renormalization techniques can be applied for generating stable higher order approximations in this limit to obtain generalized-Fokker-Planck equations (Pomraning, 1996), (Prinja and Pomraning, 2001), and (Leakes and Larsen, 2001). The scattering kernel can be decomposed into smooth and singular parts (Caro and Ligou, 1983), (Landesman and Morel, 1989), (Aristova and Gol’din, 1998), (Dixon, 2015) to derive the Boltzmann-Fokker-Planck or Boltzmann-Fokker-Planck-like approximations.

In this paper, we primarily focus on the acceleration side of the solving forward-peaked problems. Problems with forward-peaked scattering kernels require acceleration of all slowly-converging Legendre-moments of angular flux with significant magnitudes. We develop and test a synthetic acceleration method - Fokker-Planck synthetic acceleration (FPSA) - where the lower-order approximation for the error-correction stage is obtained using asymptotic analysis (Bender and Orszag, 1978) in the Fokker-Planck limit. We call this approach to acceleration asymptotics-based acceleration.

We organize the remainder of this paper as follows. In the next section, we introduce the FPSA and describe how we discretize the angular Laplacian term of the Fokker-Planck equation - weighted finite difference (Morel, 1985) and moment preserving discretization (Warsa and Prinja, 2012). In section 3, we present angularly-continuous and angularly-discrete Fourier analyses for FPSA and contrast them. Thenafter, in section 4, we present an efficiency study for FPSA screened Rutherford kernel (SRK) (Pomraning, 1992), the exponential kernel (EK) (Prinja et al., 1992) and Henyey-Greenstein kernel (HGK) (Henyey and Greenstein, 1941). We conclude this paper with a summary in section 5.

2 Fokker-Planck Synthetic Acceleration

The ideas presented here can be extended to problems with time and energy dependence, in multi-dimensions, and in curvilinear coordinates, but our presentation takes place in the context of steady-state, monoenergetic, one-dimensional Cartesian coordinates. We use standard notation and assume cgs units (Lewis and Miller, 1984), such that for a domain $z\in[z_{0},z_{1}]$ the transport equation is

[TABLE]

The cross section $\sigma_{s}(\mu_{0})$ depends on the cosine of the laboratory frame scattering angle, $\mu_{0}={\hat{\Omega}}^{\prime}\!\cdot\!\hat{\Omega}$, for a particle traveling with incident direction ${\hat{\Omega}}^{\prime}$, and exiting after a scattering event in direction $\hat{\Omega}$. Typically, this dependence is represneted with an expansion in Legendre polynomials, whose expansion coefficients are

[TABLE]

Assuming the expansion is truncrated at order $L$, and using the addition theorem for the normalized spherical harmonics, Eq. (1a) becomes

[TABLE]

where the scalar flux moments are

[TABLE]

Under certain restrictions on the scattering cross section and its expansion, the Fokker-Planck equation is an asymptotic limit of the Boltzmann transport equation when scattering is highly forward-peaked (Pomraning, 1992). The slab geometry Fokker-Planck equation is

[TABLE]

and the momentum transfer (or transport) cross section is $\sigma_{tr}=\sigma_{s,0}-\sigma_{s,1}$.

2.1 Standard Solution and FPSA

We demonstrate the ideas behind source iteration and synthetic acceleration (Adams and Larsen, 2001) before describing FPSA. Source iteration is one of the simplest methods used to solve the Boltzmann equations. We begin by rewriting (3):

[TABLE]

Source iteration for the Boltzmann equation is then written as:

[TABLE]

where $m$ is the iteration index. Fourier analysis of source iteration (Valougeorgis et al., 1988) returns $\frac{\sigma_{s,l}}{\sigma_{t}}$ as the eigenvalues of the iteration matrix. Depending on the scattering kernel, the spectral radius can approach unity. This could make solution extremely expensive. In order to optimize the performance of such schemes, the solution must be accelerated. The eigenvalues also suggest that for forward-peaked kernels, where $\frac{\sigma_{s,l}}{\sigma_{t}}$ decay slowly, accelerating higher order moments with nonzero magnitudes becomes essential. In order to understand synthetic acceleration (Kopp, 1963) consider (7). We break the solution procedure as follows:

[TABLE]

Different choices of $F$ operator return different synthetic acceleration schemes. For example choosing $F$ as the diffusion operator returns DSA. For this paper, we choose $F$ as the Fokker-Planck operator:

[TABLE]

2.2 Discretization

In order to discretize the Boltzmann and FP equations, we use linear discontinuous finite element discretization (LD) in space (Warsa, 2014) and discrete ordinates ($\text{S}_{N}$ ) in angle (Lewis and Miller, 1984). Moreover in order to discretize $\mathcal{L}_{\text{FP}}$, we use weighted finite difference (WFD) (Morel, 1985) and moment preserving discretization (MPD) (Warsa and Prinja, 2012). For a more detailed presentation on spatial and angular discretization, we refer the readers to (Patel, 2016). For convenience, we briefly review angular discretization of FP equation based on (Morel, 1985) and (Warsa and Prinja, 2012). We use $\text{S}_{N}$ quadrature to discretize the FP equation (5) in angle by collocating the angular flux at the directions $\mu_{n}$

[TABLE]

for $n=1,\ldots,N$. A way to discretize the term $\nabla^{2}_{n}\psi(z)$, which denotes the discrete form of the angular Laplacian operator Eq. (5b) evaluated at angle $n$, has to be defined in terms of the $\text{S}_{N}$ quadrature points and weights.

The three-point, WFD scheme for the angular Laplacian has weights that identical to those used for the $\text{S}_{N}$ discretization of the one-dimensional spherical coordinates (Morel, 1985). The WFD scheme is given by the following expressions.

[TABLE]

where

[TABLE]

[TABLE]

and where the $\text{S}_{N}$ quadrature normalization is $\nu$.

To develop the MPD method, we first recognize that in one dimension the Legendre polynomials are eigenfunctions of the Fokker-Planck operator

[TABLE]

For a twice-differentiable function $f(\mu)$, integrating twice by parts shows that the angular Laplacian operator is self-adjoint with respect to the Legendre polynomials:

[TABLE]

Substituting (14) in (15), the following integral relationship is readily obtained

[TABLE]

We now evaluate this relationship with $\text{S}_{N}$ quadrature for the angular flux $\psi_{n}(z)$ to get

[TABLE]

for $\quad l=0,\ldots,N-1$. This defines an $(N\times N)$ operator for the vector of $N$ angular fluxes at the spatial location $z$, $\Psi(z)$, such that the result is the Fokker-Planck operator collocated at all $N$ quadrature points simultaneously. That is,

[TABLE]

The method is a similarity transformation that by definition preserves the moments of the flux up to order $N-1$. Using Eq. (18), the MPD method for the $\text{S}_{N}$ approximation to the Fokker-Planck equation is, in operator notation,

[TABLE]

where the operator $\mathbf{M}$ is

[TABLE]

and $Q(z)$ is a vector of source terms $q_{n}(z)$ for $n=1,\ldots,N$.

Notice that we may also write the WFD operator in the manner of Eq. (18a), that is,

[TABLE]

for $i,j=1,\ldots,N$. The WFD scheme for the $\text{S}_{N}$ approximation to the Fokker-Planck equation can then be written in operator notation as

[TABLE]

As observed in (Morel, 1985), we see that the WFD scheme results in a diagonally-dominant M-matrix such that the transport operator is inverse-positive (neglecting spatial discretization). Even though the MPD operator does not have a similar simple structure that allows us to show it so easily, we have observed numerically that it is in fact inverse-positive.

3 Fourier Analysis for FPSA

Now that we have some idea of what FPSA is and how we discretize equations, we analyze FPSA using angularly-continuous and angularly-discrete Fourier analysis in this section. The goal of fourier analysis is to theoretically deretmine how efficient our method is. We begin by determining the error eqution for FPSA. We assume constant material properties throughout this exercise. We also drop notation for $z$ and $\mu$ dependence henceforth for conveneince.

Consider Eq. 8b. Upon subtracting the exact solution $\psi$ from both sides and adding and subtracting $\psi$ from the scattering term, we have:

[TABLE]

Upon introduction of the following error definitions in Eq. (22):

[TABLE]

simplification, and rearrangement, we get the following error equation:

[TABLE]

Based on whether we analyze the equation via Legendre-moments of angular error or by discretizing angular error with $\text{S}_{N}$ quadrature, we get angularly-continuous or angularly-discrete analysis.

3.1 Angularly-Continuous Fourier Analysis

Angularly-continuous or $\text{P}_{L}$ -based Fourier analysis is inspired by (Valougeorgis et al., 1988). We begin by defining error moments:

[TABLE]

In terms of moments, the error equation is:

[TABLE]

We follow the following general steps:

Obtain an expression for $\epsilon_{l}^{m+\frac{1}{2}}$ by analying the predictor step. 2. 2.

Obtain an expression for $\int_{-1}^{1}d\mu P_{l}(\mu)F^{-1}S(\epsilon^{m}-\epsilon^{m+\frac{1}{2}})$ by analyzing the corrector step. 3. 3.

Combine results from previous steps to obtain the iteration matrix $I_{M}$ such that the error is written according to $[\epsilon_{l}^{m+1}]=I_{M}[\epsilon_{l}^{m}]$.

The spectral radius of $I_{M}$ determines the convergence rate of the iterative method (Hageman and Young, 1984). This is because of the following relation:

[TABLE]

Step 1: In order to proceed, we note that $\epsilon_{l}^{m+\frac{1}{2}}$ comes from the error moment equation of the predictor step Eq. (8a). We obtan that equation by subtracting the exact transport equation from Eq. (8a), and defining error according to Eq. (23):

[TABLE]

Now, we separate the error components into their angle and space dependent components by writing $\epsilon^{m+\frac{1}{2}}$ and $\epsilon^{m}$ as Fourier integral (Adams and Larsen, 2001):

[TABLE]

where, $\lambda$ is the wave number. Substituting this form of error into the error equation Eq. (28) returns:

[TABLE]

Simplifying the above equation and noting that Fourier modes, $e^{i\lambda\sigma_{t}z}$, are linearly independent for all $\lambda$, we obtain (Adams and Larsen, 2001):

[TABLE]

Dropping $\lambda$ and $\mu$ in notation of $\hat{\epsilon}$ in Eq. (31) for convenience, and using definitions of error-moments, we get:

[TABLE]

Then rearranging the equation and taking $n^{th}$ Legendre moment of Eq. (32), we obtain the following:

[TABLE]

Further rearrangement and use of definition of error moments returns:

[TABLE]

We note that Eq. (34) represents the following matrix equation:

[TABLE]

where, $[\hat{\epsilon}_{l}^{m}]$ is a vector of error-moments at iteration m, and

[TABLE]

is an iteration matrix. We multiply Eq. (35) by $e^{i\lambda\sigma_{t}z}$ and use Eq. (25) to get:

[TABLE]

Now that we have an equation for $\epsilon_{l}^{m+\frac{1}{2}}$, we move on to the next step.

Step 2: We begin from Eq. (24). We note that the correction $\upsilon^{m+1}=F^{-1}S(\epsilon^{m}-\epsilon^{m+\frac{1}{2}})$ comes from the solution of the following equation:

[TABLE]

Introducing the Fourier mode ansatz:

[TABLE]

Upon introduction of Eq. (25), and Eq. (39) in Eq. (38), we get:

[TABLE]

Simplifying Eq. (40), taking its Legendre moment, and using the orthogonality property of Legendre polynomials returns:

[TABLE]

Now, using the recurrence relation for Legendre polynomials on the first term of Eq. (41), expanding $\hat{\upsilon}_{\lambda}^{m+1}(\mu)$ in the third term of Eq. (41) using Legendre expansion, we get:

[TABLE]

Simple rearrangement of the third term in Eq. (42), followed by use of Legendre’s equation, and orthogonality property of Legendre polynomials returns:

[TABLE]

where,

[TABLE]

Eq. (43) can be written in matrix form as:

[TABLE]

where,

[TABLE]

Multiplying Eq. (45) with $e^{i\lambda z\sigma_{t}}$ and using Eq. (39) returns:

[TABLE]

Step 3: Combining Eqs. (26), (37), and (47), returns:

[TABLE]

Comparing Eq. (48) with Eq. (27) returns the iteration matrix $I_{M}=(A-B^{-1}XE)$. The spectral radius of $I_{M}$ is the spectral radius of FPSA.

3.2 FPSA as a Special Case of $\text{P}_{L}$ Acceleration

Upon carrying out similar analysis for $\text{P}_{L}$ acceleration (Valougeorgis et al., 1988), we find that the iteration matrix for $\text{P}_{L}$ acceleration has a similar form except, in this case, the definition of $B_{l,l}$ is slightly different:

[TABLE]

That for FPSA is rewritten as:

[TABLE]

When we equate the two equations, we see that FPSA is a special case of $\text{P}_{L}$ acceleration when:

[TABLE]

Another way of obtaining this equivalence relation is by noting that Legendre polynomials are eigenfunctions of both Boltzmann scattering operator and the Fokker-Planck operator as done by Morel in (Morel, 1981):

[TABLE]

and equating the eigenvalues of Fokker-Planck and the Boltzmann scattering operators:

[TABLE]

Simple rearrangement of Eq. (53) returns Eq. (51). We will call these cross-section moments $P_{L}$-equivalent cross-section moments in this paper.

According to the $S_{N}-P_{L}$ equivalence relation in slab geometry (Lewis and Miller, 1984), when $N=L+1$, $S_{N}$ and $\text{P}_{L}$ equations are equivalent. Taking this and Eq. (51) into account, we note that FPSA will converge in one iteration when it is analytically equivalent to $\text{P}_{L}$ acceleration. In other words, when the scattering cross-section moments are according to Eq. (51), FPSA will converge in one iteration. Moreover, it would be a valid to think that the convergence will be rapid in case the cross-section moments are close to those obtained from Eq. (51). However, in the case when we truncate scattering expansion arbitrarily and $N$ is no longer equal to $L+1$, the FPSA-$\text{P}_{L}$ acceleration equivalence will no longer hold. This is due to the inconsistent introduction of zero values for scattering crosssection moments with $N\geq l>L$ (Patel, 2016).

3.3 Angularly-Discrete Fourier Analysis

Angularly-discrete analysis is carried out by using $\text{S}_{N}$ quadrature to approximate angular error. We need angularly-discrete Fourier analysis to analyze FPSA because different discretizations of $\mathcal{L}_{\text{FP}}$ preserve different number of moments. While WFD only preserves zeroth and first Legendre moments of the angular flux (Morel, 1985), MPD preserves upto $N$ Legendre moments (Warsa and Prinja, 2012). The angluarly-continuous Fourier analysis ($\text{P}_{L}$ -based analysis) is moment-based and therefore requires the numerical implementation to preserve all relevant moments in order to get a consistent spectral radius measurement. Moreover, in case of continuous transport, the transport equation only limits to the Fokker-Planck equation when a ”sufficient” number of Legendre moments are used to represent the angular flux (Patel, 2016). This sufficient number of moments is scattering cross-section dependent. This, however, may not necessarily be true in the discrete $(S_{N})$ case. This creates a discrepancy between the angularly-continuous and angularly-discrete Fourier analyses when sufficient number of moments are not used to represent angular flux. Therefore, in order to verify convergence rates of the numerical implementation, irrespective of how many moments are used to represent angular flux, we introduce angularly-discrete analysis. First, we consider Fourier analysis for FPSA with WFD for $\mathcal{L}_{\text{FP}}$.

3.3.1 Analysis with WFD

We follow the following analogous steps to do angularly discrete Fourier analysis:

Obtain an expression for $\epsilon^{m+\frac{1}{2}}$. 2. 2.

Obtain an expression for $F^{-1}S(\epsilon^{m}-\epsilon^{m+\frac{1}{2}})$. 3. 3.

Obtain the overall matrix equation that is used to estimate the the spectral radius.

Step 1: Since we have already detailed angularly-continuous analysis, we skip furnishing the introduction of Fourier mode assumption and simplification steps here. We also ignore notation of $\mu$ and $z$ dependence of relevant quantities for convenience. Taking the $n^{th}$ Legendre moment of Eq. (32), and using orthogonality property of Legendre polynomial returns:

[TABLE]

Now we write each integral as a discrete weighted-sum using $\text{S}_{N}$ quadrature:

[TABLE]

Finally, we get the following matrix equation from Eq. (55):

[TABLE]

where,

[TABLE]

and,

[TABLE]

Here, $\hat{A}$ is the iteration matrix in angularly-discrete from. This returns $\frac{\sigma_{s,0}}{\sigma_{t}}$ as the spectral radius which is consistent with the angularly-continuous analysis (Patel, 2016). We multiply Eq. (59) by the relevant exponential from Fourier mode ansatz to get:

[TABLE]

Step 2: Upon introduction of Fourier mode assumption for $\upsilon$ in Eq. (38), taking Legendre moment of equation, using definition of error moments, carrying out the relevant spatial differentiation and simplifications, we get:

[TABLE]

Now we write each integral in the form of a weighted sum and the angular differential using the weighted difference formulation (Morel, 1989) to obtain:

[TABLE]

where, $a_{n},b_{n},$ and $c_{n}$ are according to Eq. (20b). From Eq. (61), we get the following matrix equation:

[TABLE]

where,

[TABLE]

with,

[TABLE]

[TABLE]

and,

[TABLE]

We obtain the following expression for $F^{-1}S(\epsilon^{m}-\epsilon^{m+\frac{1}{2}})$:,

[TABLE]

Step 3: Combining Eqs. (66), (59), and (24) returns:

[TABLE]

where, $I_{M}=\hat{A}-\hat{B}^{-1}\hat{C}\hat{D}$ is the iteration matrix and its spectral radius determines the convergence rate of FPSA with WFD. Now, we consider angularly-discrete analysis for FPSA with MPD.

3.3.2 Analysis with MPD

Angularly discrete Fourier analysis for FPSA with MPD is done in the same way as for FPSA with WFD. The only difference will be how the Fokker-Planck operator is represented in step 2. Introducing angularly discrete formulation for integrals and MPD formulation (Warsa and Prinja, 2012) for the angular Laplacian in Eq. (60) returns:

[TABLE]

From Eq. (61), we get the following matrix equation:

[TABLE]

where,

[TABLE]

[TABLE]

and,

[TABLE]

We have the following expression for $L^{-1}S(\epsilon^{m}-\epsilon^{m+\frac{1}{2}})$:

[TABLE]

Step 3: Just like with previous analyses for FPSA, we get:

[TABLE]

Thus the iteration matrix for FPSA with MPD is $A-B^{-1}CD$.

3.4 Comparison of Spectral Radii

In order to get a glimpse into how FPSA performs, we consider one problem with screened Rutherford kernel (SRK) (Dixon, 2015), exponential kernel (EK) (Prinja et al., 1992) and Henyey-Greenstein kernel (HGK) (Pomraning, 1992) each. We plot spectral radii of stand-alone $\text{S}_{N}$ and FPSA for each kernel in the figures that follow. We choose $L=15$, $N=16$. We note that the spectral radius reduces significantly upon introduction of FPSA. The spectral radius reduction is completely problem dependent. The spectral radius can potentially change with N, L, and how close the scattering cross-section moments of the problem are to the $P_{L}$-equivalent moments.

Next, we compare the numerically measured and theoretical (angularly-discrete Fourier analysis) spectral radii. We analyze convergence rates for three scattering kernels - SRK, EK, and HGK. We choose $L=15$, $N=16$. We use a slab of length, 100 cm, discretize it using 100 elements. We use vacuum boundaries for numerical measurements of spectral radius. The theoretical and numerically measured spectral radii have been presented in Table 1.

We obtain similar theoretical and measured spectral radii values for different scattering kernels with varying parameters. This indicates a relatively accurate analysis of the method.

4 Efficiency Study

In this section we will assess how the reduction in spectral radius results in reduction in runtime of source iteration (SI) and GMRES solves. We run all problems using MATLAB and track runtime using its tic-toc functionality. We place tic and toc before and after the solver function calls respectively. In other words, we do not include the stiffness matrix setup time in our calculation. We will only account for the solver runtime. Specifically, choose problems with $L=15$, and $N=16$, $32$. We use beam and vacuum boundaries. We have a unit distributed source for problems with vacuum boundaries and a unit beam source with the beam boundary. We do this for SRK with $\eta=2.83\times 10^{-5}$, and for EK with $\Delta=10^{-5}$. We solve the Fokker-Planck error equation (invert the preconditioner) using LU factorization via factorize object (Davis, 2009) in MATLAB, and GMRES.

First, we compare unpreconditioned SI and GMRES solves. In order to compare these solves, we choose $\eta=2.83\times 10^{-5}$, $L=15$, $N=16$, $H=1$cm, $K=100$, $tol=10^{-10}$. We do this to contrast source iteration and GMRES solves. Table 2 and 3 present this data. It is clear that GMRES is more suitable than source iteration for forward-peaked transport problems.

Next, we compare solution rutimes and iteration counts. We will compare these for unpreconditioned GMRES, FPSA-preconditioned SI, and FPSA-preconditioned GMRES solves. We do not include unpreconditioned source iteration in this study because its ineffectiveness for relevant problems has already been demonstrated in Table 2 and 3. We will arbitrarily choose our restart parameter for this study to be 150.

4.1 Screened Rutherford Kernel

We compare efficiency of FPSA for problems involving SRK in this section. We choose a slab of unit length discretized using hundred elements. We choose $\eta=2.83\times 10^{-5}$ and $\sigma_{a}=1$. Scattering cross-section moments come from SRK and their numerical values can be found in (Patel, 2016). Finally, $L=15$, and $N=16$ and $32$. Number of iterations and overall runtime data has been presented in Table 4, 5, 6, and 7.

We observe a significant decrease (almost three orders of magnitude compared to unpreconditioned GMRES and five orders of magnitude compared to SI) in the number of transport-sweeps required for convergence due to preconditioning. We also observe a decrease in overall solver runtimes due to preconditioning when FP-solve is done using LU factorization (by upto two orders of magnitude compared to unpreconditioned GMRES). The FP-solve, however, can be extremely expensive and render this preconditioner ineffective with respect to problem’s overall runtime if inefficient solvers are used. Here, the number of iterations required for one FP-solve using GMRES was of the same order as an unpreconditioned transport solve using GMRES. It is imperative that we find an effective preconditioner for FP-solves. We are looking into this. The potential, however, of using FP as preconditioner for transport solves is amply evident from the data presented in this section. Next, we look at efficiency data for problems with the exponential kernel.

4.2 Exponential Kernel

We calculate scattering cross-section moments using EK for $\Delta=10^{-5}$. The zeroth moment is calculated using SRK. The study is done using the same parameters as SRK except for scattering cross-section moments. Number of iterations and overall runtime data has been presented in Table 8, 9, 10, and 11.

We see similar behavior to what we saw in the case of SRK. The solver runtimes differ due to difference in rate at which FP-solve converges for this particular problem. Again, we note a significant decrease in number of iterations but a decrease in solver runtime strongly depends on the efficiency of the FP-solve.

4.3 Henyey-Greenstein Kernel

In this section, we let the asymmetry parameter, $g=0.9999$. The study is carried out in the same way as the previously for SRK and EK. For this section, we will choose $\sigma_{a}=0.00001$ $cm^{-1}$. The scattering cross-section moments are calculated using HGK. We will choose slab length of $50$ cm disretized using $200$ elements. Number of iterations and overall runtime data has been presented in Table 12, 13, 14, and 15.

We note that, just like for SRK and EK, preconditioned schemes have significantly less iteration counts. However depending on how the Fokker-Planck error equation is solved, the preconditioning may or may not be effective with respect to runtime reduction. Solving the FP equation with GMRES renders FPSA scheme unviable, however use of factorization reduces to overall runtime significantly.

5 Summary and Future Work

We ran several numerical experiments and assessed the speed-ups in iteration count and solver runtime. We saw that preconditioning transport solve using FP resulted in reduction in iteration count by upto three orders (when compared to unpreconditioned GMRES solves). The overall runtime, however, depended completely on how efficiently the FP preconditioner was solved. Direct factorization resulted in a runtime reduction by upto two orders of magnitude. We observed that FP can be a very effective preconditioner for transport solves with highly forward-peaked scattering. However, we must develop an effective solver for FP-solve itself in order to make this an attractive preconditioning method. In future, we would like to determine how do we optimize FP-solve. We would also like to test FPSA’s performance in energy dependent, multi-D settings. Moreover, we would also like to develop a nonlinear version of this method which would allow us to obtain a Fokker-Planck equation that is consistent with the relevant transport equation.

Acknowledgments

This information has been co-authored by an employee or employees of the Los Alamos National Security, LLC. (LANS), operator of the Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396 with the U.S. Department of Energy.