Solving nonlinear discretizations of SN transport calculations

Annals of Nuclear Energy xxx (xxxx) xxx

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Solving nonlinear discretizations of SN transport calculations Michael W. Hackemack Naval Nuclear Laboratory, P.O. Box 1072, Schenectady, NY 12301-1072, USA

a r t i c l e

i n f o

Article history: Received 20 July 2019 Received in revised form 14 September 2019 Accepted 21 September 2019 Available online xxxx Keywords: Neutron transport Nonlinear spatial discretization Newton-like iterations Diffusion acceleration

a b s t r a c t In this paper, we investigate iterative strategies to converge transport solutions using nonlinear discretization schemes, specifically the adaptive-weighted diamond-difference method and a set-to-zero fixup approach for the linear discontinuous method compatible with arbitrary polytopal cells. The nonlinear nature of these schemes precludes the use of linear methods, such as highly effective Krylov iterative solvers like GMRES. Instead, nonlinear Newton-like solvers such as Broyden’s method, nonlinear Krylov acceleration, and the Jacobian-free Newton-Krylov method are analyzed in conjunction with physics-based diffusion acceleration schemes. The Richardson method (i.e., source iteration) is also analyzed for comparison. Numerical experiments (both monoenergetic and multigroup) demonstrate effective reduction in transport work by the use of appropriately chosen iterative schemes, including choosing sufficiently consistent diffusion acceleration. These results are analogous to linear transport operators, and they can be used as a point of reference in determining solution strategies for other nonlinear transport discretization schemes. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction The linearized Boltzmann transport equation, which we will simply call the transport equation, describes the transport of neutral particles (i.e., neutrons and photons) through background media. It has a large phase-space consisting of three spatial dimensions, 2 angular dimensions, 1 energy dimension, and 1 temporal dimension (for unsteady problems) (Lewis and Miller, 1984). The discrete ordinates (SN ) method is perhaps the most widely utilized angular discretization for the transport equation. However, it is well known that unrefined spatial discretizations or inadequate Legendre polynomial expansion of the anisotropic scattering source can produce negative flux solutions (Lathrop, 1969). While these negative flux solutions are mathematically valid for a given phase-space discretization, they can be extremely undesirable. Negative fluxes within regions of interest can give poor or even non-physical quantities of interest (e.g., reaction rates within a detector region). Additionally, common acceleration schemes may be adversely affected by negative fluxes (Smith, 1983; Adams and Larsen, 2002). Over the years, several methods have been developed to either mitigate or eliminate flux negativities. These methods can be divided into two broad categories: (1) non-negative nonlinear moment methods and (2) ad hoc ‘‘fixup” methods.

E-mail address: [email protected]

Discretizations that can be defined as non-negative nonlinear moment methods have inherently smooth nonlinearities, with examples including the linear exponential characteristic methods (EC) (Minor and Mathews, 1995) and the linear exponential discontinuous finite element method (ED) (Wareing, 1997). Due to their smoothness, these methods can be efficiently solved with Newton’s methods. However, both have significant drawbacks. Characteristic methods are expensive, especially on non-orthogonal meshes, and cannot be applied to curvilinear coordinate systems. Conversely, exponential discontinuous finite element methods can be applied to curvilinear geometries, but their accuracy can suffer in fine mesh limits relative to analogous finite element schemes (e.g., the linear discontinuous finite element method (LD) (Hill, 1975)). On the other hand, ad hoc fixup methods, which we now simply refer to as fixup methods, differ from the non-negative nonlinear moment methods by having non-smooth closure relations. Due to their ad hoc nature, many of these types of methods have been developed. The well-known diamond difference method can be fixed with a set-to-zero scheme (Lathrop, 1969) or an adaptiveweighted closure (Alcouffe, 1993). Likewise, the LD method has had many fixup methods developed, which can be separated into two broad categories: consistent (Maginot et al., 2012) and inconsistent (Hamilton et al., 2009; Hackemack and Becker, 2017) moment methods. Consistency refers to the method’s ability to continue satisfying the moment equations upon fixup (the balance equation naturally must still be preserved). The methods in

https://doi.org/10.1016/j.anucene.2019.107080 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

2

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

Hackemack and Becker (2017) ensure positivity of the cell averages on polytopal mesh cells. The consistent set-to-zero (CSZ) method (Maginot et al., 2012) was derived for 1D and 2D quadrilateral meshes and ensures positivity everywhere in a cell. In a similar manner, the bilinear consistent set-to-zero (BCSZ) method extends the CSZ method to bilinear discontinuous finite elements on 2D quadrilaterals (Maginot et al., 2017). Finally, it is noted that the methodology in Hamilton et al. (2009) can be used to fixup other polytopal-compatible discontinuous finite element method (DFEM) discretizations, such as those in Bailey (2008) and Hackemack and Ragusa (2018). The linear SN transport equations have historically been solved via source iteration (which is simply the fixed-point form of Richardson iteration with a unity relaxation parameter) (Lewis and Miller, 1984; Adams and Larsen, 2002). Each source iteration constitutes performing a ubiquitously-known transport sweep, which requires the matrix-free inversion of an independent block lower-triangular system for each direction (Baker and Koch, 1998; Hawkins et al., 2012). In this work, it is assumed that the mesh topology and boundary conditions do not produce any cycles in the sweep graphs. Since the discretizations employed in this work are dependent on the flux solution, their transport sweeps are inherently nonlinear. Therefore, the transport operators and the system of equations to solve are nonlinear. Unfortunately, this precludes the direct application of well-known Krylov iterative methods (e.g., GMRES or BiCGStab), which have been shown to be highly effective linear solvers (Guthrie et al., 1999; Patton and Holloway, 2002). It is well known that Richardson iterations without any preconditioning (i.e., acceleration) for transport problems with opticallythick and diffusive regions can converge prohibitively slowly (Adams and Larsen, 2002). To accelerate these difficult transport problems, there are two broad approaches that can be applied: (1) use nonlinear fixed-point and Newton-like iterative methods and (2) apply some form of physics-based preconditioning (e.g., diffusion acceleration). Combinations of these two approaches can also be utilized. Besides Richardson iteration, other nonlinear iterative methods include Broyden’s method (Broyden, 1965), Anderson acceleration (Anderson, 1965), nonlinear Krylov acceleration (NKA; Carlson and Miller, 1998), and Jacobian-Free NewtonKrylov (JFNK; Knoll and Keyes, 2004). These methods have been applied to nonlinear fixup methods in the works of Hackemack (2018) and Fichtl et al. (2010). They have also been extensively studied for nonlinear eigenvalue transport problems (Knoll et al., 2011; Gill et al., 2011; Calef et al., 2013). It has been shown that physics-based transport preconditioners provide effective acceleration to slowly converging transport problems (Adams and Larsen, 2002). In general, these methods accelerate convergence by projecting information from a lower-order dimensional space onto the iterating transport solution. Of these methods, diffusion acceleration methods have become popular because of their low cost and effectiveness. There are two general categories of diffusion acceleration schemes: linear diffusion synthetic acceleration (DSA) and nonlinear diffusion acceleration (NDA). Linear DSA acts by additively correcting the error between two successive transport iterations (Reed, 1971; Alcouffe, 1977; Wang and Ragusa, 2010); variants for the fission source have also been developed (Suslov, 2003; Barbu and Adams, 2019). If a linear transport operator were used, then source iteration with DSA is exactly a preconditioned Richardson scheme. On the other hand, the NDA class of methods, with the coarse mesh finite difference (CMFD) method being a popular variant (Smith, 1983), corrects the high-order transport solution with multiplicative prolongations of equivalently-derived low-order solutions (Smith and Rhodes, 2002). Unfortunately, the NDA method becomes unstable and diverges for optically thick mesh cells just like cell-based

DSA methods (Reed, 1971). To stabilize NDA, a partial-current nonlinear diffusion acceleration (pNDA) formulation can be used which effectively under-relaxes the scheme (Cho et al., 2003). A lot of the corresponding analysis of the NDA methods has occurred on Cartesian mesh cells utilizing the same spatial basis between the transport and diffusion problems (e.g., a cell-based finitedifference or finite-volume method). Other work has extended these methods to unstructured grids (Kim and DeHart, 2011) as well as using different spatial bases between the high-order and low-order problems (Lee and Joo, 2014). Using combinations of advanced iterative solvers with diffusion preconditioning has also been investigated. Linear transport problems solved with GMRES using consistent and partially-consistent DSA schemes were analyzed by Warsa et al. (2004). They showed that the use of a Krylov method stabilized the preconditioned Richardson iterations that suffered degraded effectiveness (or even diverged) because of either strong material discontinuities or the partial-consistency of the DSA scheme. While this is not exactly applicable to our work due to the nonlinear transport operators, the use of the Broyden, Anderson, NKA, and JFNK methods act as ‘‘nonlinear Krylov” schemes which can stabilize divergent diffusion acceleration schemes (both DSA and NDA schemes). This behavior is demonstrated with numerical results later in this work. The use of a NDA scheme, even with linear discretization schemes, makes the fixed-point transport operator nonlinear. JFNK (Knoll et al. (2011)) and Anderson acceleration (Willert et al., 2014) methods have been applied to accelerate the NDA iterative scheme. Finally, Bruss et al. (2014) utilized a linear S2 synthetic acceleration (S2 SA) scheme to precondition the 1D CSZ equations. We utilize an analogous methodology in this work whereby linear DSA operators precondition our employed nonlinear transport operators. The primary goals of this paper are to extend the work of Hackemack (2018) to form a more complete analysis of iterative techniques for the solution of nonlinear SN discretizations. This extension will be done in several ways. First, we analyze the efficacy of physics-based diffusion preconditioning with the chosen nonlinear iterative methods: Richardson, Broyden, NKA, and JFNK. Second, we extend our analysis to discretization strategies compatible for unstructured polytopal meshes. Third, we analyze the iterative performance to a generalized multigroup framework along with appropriate multigroup diffusion preconditioning. Finally, we correct some results from Hackemack (2018) stemming from an incorrect implementation of the Anderson method. The NKA and Anderson methods are mathematically equivalent (using unity weighting factors in the Anderson method) except for their condition monitoring (Calef et al., 2013), and as such should provide identical iterative results. Therefore, NKA (but not Anderson) is analyzed in this work. The remainder of this work is organized as follows. Section 2 presents the SN transport equation along with brief descriptions of the utilized nonlinear discretization methods. Next, Section 3 details the transport iterative operations, including physics-based diffusion preconditioning, to be used for the nonlinear iterative solvers described in Section 4. Finally, numerical results are given in Section 5 as well as some concluding remarks in Section 6. 2. The transport equation 2.1. The discretized multigroup SN equations Given an angular quadrature set with M directions, n oM ~m wm ; X , the angular flux solution of group g traveling in m¼1   ~m , of the multigroup SN transport direction m; wm;g ð~ r Þ  wg ~ r; X equation with isotropic scattering and with G number of groups

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

3

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

within an open, convex spatial domain ~ r 2 D, with boundary @D, is given by





~m  r ~ þ rt;g ð~ X rÞ wm;g ð~ rÞ ¼

0 G X rgs ;g ð~rÞ

g 0 ¼1

4p

/g0 ð~ rÞ þ Q m;g ð~ r Þ;

ð1Þ

0

/g ð~ rÞ ¼

M   X ~  dXwg ~ wm wm;g ð~ r; X rÞ:

ð2Þ

m¼1

4p

The boundary conditions are defined for any point on the inflow n o ~ ~ ~ boundary: ~ r b 2 @D m ¼ @D such that Xm  nb < 0 , where nb is the outward unit normal on the boundary at ~ r b . Two types of boundary conditions are employed: incident boundaries (Cinc ) and reflecting boundaries (Cref ), such that @D ¼ Cinc [ Cref . Up to this point, Eq. (1) has been discretized in terms of angle and energy only, with no loss of generality in the spatial domain. Next, we lay down a mesh, Th 2 Rd , over the spatial domain, where d is the spatial dimension (d ¼ 2 in this work). This mesh consists of non-overlapping spatial elements to form a complete S union over the entire spatial domain: D ¼ K2Th K. Likewise, the faces of the mesh are formed as a union of the interior faces, Cint , and the boundary faces: C ¼ Cint [ Cinc [ Cref . Without providing a specific discretization method at this time (details are provided in Section 2.2), we discretize Eq. (1) using a full upwind scheme (Lewis and Miller, 1984). Eqs. (1) and (2) can be written in operator form as

LðUÞW ¼ MSU þ Q ;

U ¼ DW;

ð3Þ

where

L¼

the discrete collision and streaming operator;

S¼

the discrete scattering operator;

M ¼ the angular moment-to-discrete operator; D¼ Q ¼

the angular discrete-to-moment operator; the discrete angular fixed source;

W ¼ the discrete angular flux; U ¼ the moment-based angular flux: Eq. (3) can be rearranged to only be in terms of the angular moments:

  I  DLðUÞ1 MS U ¼ DLðUÞ1 Q :

 XXm;j  þ 
   D  wm;g;K;j  wm;g;K;j þ rt;g;K wm;g;K ¼ Q m;g;K ; j

where rt;g is the total cross section of group g; rgs g is the scattering cross section from group g 0 into group g; Q m;g is the distributed source of group g in direction m which can include contributions from fixed or fission sources (i.e., eigenvalue problems), and /g is the scalar flux of group g given by

Z

sian meshes. For a mesh cell K with span of DK;j in dimension j, the WDD method yields a balance equation of

ð4Þ

ð5Þ

K;j


m;g;K and Q m;g;K are the cell average flux and source, respecwhere w þ tively, and w m;g;K;j and wm;g;K;j are the average inflow and outflow fluxes in dimension j, respectively. To eliminate the outgoing fluxes in Eq. (5), standard closure is given by


m;g;K ¼ w

1 aj wþ þ w ; 1 þ aj m;g;K;j 1 þ aj m;g;K;j

ð6Þ

where aj 2 ½0; 1 are the dimensional weighting factors. Inserting Eq. (6) into Eq. (5) gives the average cell flux:

Q m;g;K þ
m;g;K ¼ w

X j

rt;g;K þ

X   1 þ aj  Dm;j wm;g;K;j K;j

X j

X  1 þ aj  Dm;j  K;j

:

ð7Þ

Once the average cell flux is computed with Eq. (7), then the average outflow fluxes are computed with Eq. (6). If the values of aj ; 8j remain constant, then WDD is a linear method. Indeed aj ¼ 0 and aj ¼ 1 correspond to the standard step and diamond closures, respectively. Instead, AWDD is a two-pass predictor-corrector scheme. A predictor estimate of the cell aver
p , is computed with the diamond closure (i.e., age flux, w m;g;K aj ¼ 1) in Eq. (7). Next, corrected dimensional weighting factors, acj , are computed with

1   
p   w m;g;K A; ¼ min @1;   p 
 2 wm;g;K;j  w m;g;K 0

acj

ð8Þ


p which is determined from w m;g;K as a sufficient condition to ensure non-negative outgoing fluxes (Alcouffe, 1993). Finally, these corrected dimensional weighting factors are used in Eqs. (7) and (6) to compute corrected cell and face average fluxes, respectively. 2.2.2. Linear discontinuous with flux fixup on polytopes The second nonlinear spatial discretization we consider in this work is a set-to-zero fixup method applied to a formulation of the linear discontinuous finite element discretization (which we denote as LD-STZ) developed for polytopal mesh elements (Hackemack and Becker, 2017). As is standard with discontinuous finite element schemes, the angular flux of Eq. (1) is expanded over mesh cell K as

wm;g ð~ rÞ 

dþ1 X ^ m;g;K;j bK;j ð~ w r Þ;

ð9Þ

j¼1

The operation DL1 ðUÞ is known as the transport sweep because it sweeps a given source across the spatial mesh for all discrete ordinates to obtain a new estimate of the flux moments. Due to the nonlinearity of the transport discretization, the flux solution arising from the sweep is dependent on the solution itself making the sweep operator inherently nonlinear. Next, we provide details of the nonlinear discretizations chosen for this work.

where, for a 2D problem, the expansion coefficients and basis functions with compact support on cell K are given by n oT   ~ ^ m;g;K;a ; w ^ m;g;K;x ; w ^ m;g;K;y and bK  bK;a ; bK;x ; bK;y T , respecwm;g;K  w

2.2. Nonlinear discretizations

1
xK ¼ jKj

2.2.1. Adaptive-weighted diamond difference The first nonlinear spatial discretization we consider in this work is the adaptive-weighted diamond-difference method (AWDD) (Alcouffe, 1993). It is a nonlinear predictor-corrector variant of the weighted diamond-difference method (WDD) for Carte-

tively. The average and x-moment functions are defined as bK;a ¼ 1 and bK;x ¼ ðx 
xK Þ=DxK , respectively, where the centroid and halfthickness quantities are given by

Z K

ffi vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z u u 1 2 drx and DxK ¼ t dr ðx  x
K Þ ; jKj

ð10Þ

K

respectively. The y-moment functions and quantities are computed analogously. To ease the implementation on polytopal mesh cells, a facewise coordinate system is also utilized. Analogous expressions

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

4

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

for face-wise expansion coefficients, ~ wm;g;f , and basis functions, bf , can easily be computed. The face-wise average and x-moment   basis functions are defined as bf ;a ¼ 1 and bf ;x ¼ x 
xf =Dxf , respectively, and use


xf ¼

1 jf j

Z

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 1 Z  2 u
dsx and Dxf ¼ t ds x  xf : jf j

f

ð11Þ

f

Mapping between cell K and one of its faces, f, can be performed through geometric means with a transform matrix:

2

1

6 ð
xf 
xK Þ 6 Tf ;K ¼ 6 DxK 4

ðyf yK Þ DyK

0

0

3

7 0 7 7: 5

Dxf D xK

ð12Þ

Dyf DyK

0

Therefore, cell-to-face mapping for the fluxes and basis functions is given by ~ wm;g;f ¼ TTf;K ~ wm;g;K and bK ¼ Tf ;K bf , respectively. Using all of the cell-wise and face-wise quantities so far defined, the LD weak form for cell K is given by P P ~ m;g;K ; ~ ð Fm;f ;K TTf;K  M1 Fm;f ;K ~ w"m;g;f þ Q K Gm;K þ rt IÞ wm;g;K ¼ ~ m ~ X nf >0

~m ~ X nf <0

ð13Þ where I is the identity matrix, ~ w"m;g;f are the inflow fluxes taken from the upstream cell across face f, and we use the following elementary matrix terms:

Fm;f ;K MK Gm;K

jf j ~  jKj j Xm  ~ nf j M1 K Tf ;K Mf ; R T 1  jKj dr bK bK ; K  R ~ 1 ~ bK bT ;  jKj dr Xm  r K

 jf1j

R

T

ds bf bf :

f

~m  ~ X nf > 0), an alternate transformation matrix is applied when

negatives arise. Using the set-to-zero concept, we can define a fixup transformation matrix as Tf ;K;Fix  0, where 0 is an all-zero matrix. Splitting the outflow faces between those that stay positive (f 2 pos) and those that become negative (f 2 neg), Eq. (13) can be rewritten as

Fm;f ;K TTf;K f 2pos

X

þ Fm;f ;K TTf;K;Fix f 2neg M1 K Gm;K ¼

X ~m ~ X nf <0

þ rt;g I



So far, we have provided the space-angle-energy discretizations for the multigroup SN transport equation along with the introduction of the transport sweep in Eq. (4). This section provides the details for the iterative operations required for the different solvers presented in Section 4. Consider the system of equations: FðxÞ ¼ 0, where F is nonlinear in this work (F can also be referred to as a nonlinear residual). There are many problems, the transport equation included, that can then be recast into a fixed-point form: x ¼ GðxÞ, where FðxÞ ¼ GðxÞ  x. A fixed-point iteration for the ð‘ þ 1Þ iterate can thus be cast into the following form:

  xð‘þ1Þ ¼ G xð‘Þ ;

ð16Þ

where convergence is guaranteed if G is a contraction mapping in a neighborhood around the solution and xð0Þ is some initial approximation to the solution within that neighborhood (Kelley, 1995).     The corresponding residual can be written as F xð‘Þ ¼ G xð‘Þ  xð‘Þ . These iterative operations will be expanded in detail for the transport and diffusion preconditioning operators in Sections 3.1 and 3.2, respectively. 3.1. Transport sweep operations

Eq. (13) represents the standard LD method with a system of d þ 1 equations and unknowns that is solved during a transport sweep. It is also in a form that is easily extensible for a fixup method. By mon^ m;g;f ;a , for outflow faces (i.e., itoring the average face value, w

X

3. Iterative operations

ð14Þ

K

Mf

Eq. (15) is repeatedly used until no more outgoing faces become negative. It is important to note that once an outgoing flux is flagged negative, it remains that way in that transport sweep to prevent infinite alternating flagging of faces between positive and negative. While not necessary for this work, this sweep iterative process can still be used with sweep graphs containing cycles if the appropriate upstream fluxes are lagged (Wareing et al., 2001).

~ wFix m;g;K

ð15Þ

~ m;g;K : w"m;g;f þ Q Fm;f ;K ~

In the limit of all outgoing faces being flagged negative, along with a positive Q m;g;K;a , it is mathematically guaranteed for Eq. (15) to yield a non-negative average solution: wFix m;g;K;a P 0. Determining the negative outflow faces (f 2 neg) for Eq. (15) constitutes a within-sweep iterative process. During a transport sweep for a given cell/angle, Eq. (13) is used to determine ~ wm;g;K along with the positivity/negativity of the average values of all outgoing faces. Faces with negative values are flagged and Eq. (15) calculates the fixup flux ~ wFix . The average outgoing face values are m;g;K

then determined and flagged if negative, and Eq. (15) is again used.

In Eq. (4), the operation L1 is ubiquitously known as the transport sweep. The sweep acts by tracking particles along the angular abscissa through the mesh, beginning at the incident boundaries. This way, all necessary particle information streaming into a cell from its upwind neighbors is known. Therefore, the cell’s local calculation can then be performed before sending information to its downstream neighbors. This constitutes an exact inversion of L in a matrix-free manner. Consider a subset ð g 2 GÞ of the G energy groups which we assume to be contiguous in this work (e.g., the set of thermal/ upscatter groups), though we note that this restriction is not necessary. Eq. (1) can be written for the ð‘ þ 1Þ source iteration of a group in the G subset as





‘þ1Þ ~m  r ~ þ rt;g ð~ X r Þ wðm;g ð~ rÞ ¼

X rgs 0 ;g ð~ r Þ ð‘Þ /g0 ð~ r Þ þ Q m;g ð~ r Þ; 4p g 0 2G

ð17Þ

where it is noted that Q m;g now contains inscattering contributions from groups outside of G. Eq. (17) constitutes a Jacobi iteration of the g 2 G flux moments from the previous iterate. Using the moment-only form of the discretized transport equation in Eq. (4), the ð‘ þ 1Þ fixed-point iteration of Eq. (16) can be written as







UðG‘þ1Þ ¼ Gtrans UðG‘Þ ¼ DL1 UðG‘Þ G

h

i ð‘Þ Q G þ MSG UG :

ð18Þ

The operators and fluxes have the G subscript to simply denote that we are restricted to the g 2 G subset of groups. Likewise, the corresponding ð‘Þ iterate residual is given by

    ð‘Þ ð‘Þ ð‘Þ Ftrans UG ¼ Gtrans UG  UG ;  h i ¼ DL1 UðG‘Þ Q G þ MSG UðG‘Þ  UðG‘Þ ; G ð‘þ1Þ G

¼U

ð19Þ

ð‘Þ G ;

U

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

5

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

which is utilized in the nonlinear iterative solvers presented in Section 4 for the numerical cases in Section 5 that do not use diffusion preconditioning. 3.2. Physics-based preconditioning It is well known that the fixed-point iteration scheme of Eq. (18) can lead to undesirably slow convergence for the SN transport equation for problems that are optically thick and diffusive (i.e., large material regions that are scattering dominated). Diffusion preconditioning is a standard methodology to accelerate the convergence of these slowly converging transport problems. In this work, we will analyze two broad categories of diffusion preconditioning: linear diffusion synthetic acceleration (DSA) and nonlinear diffusion acceleration (NDA), with general details provided in Sections 3.2.1 and 3.2.2, respectively. Additionally, more detailed descriptions of the diffusion discretizations, their boundary conditions, and mapping operators are given in Appendix A. Diffusion preconditioning, like most transport acceleration methods, operates by reducing the transport problem to an approximate problem of a lower-order dimensional space. The solutions from these approximate low-order problems then correct or inform the convergence of the high-order transport iterations. For the g 2 G subset of groups used in Section 3.1, a transportequivalent diffusion operator is of the form

~ ~ r J g þ rt;g ug 

X g 0 2G

rgs ;g ug0 ¼ Qg ; g 2 G 0

ð20Þ

where u is the diffusion solution, J is the current, and Q is the diffusion source term (Roberts and Forget, 2014). We analyze a full multigroup diffusion representation on the same spatial mesh as the transport since it should be the most optimal diffusion-only preconditioning strategy based on minimizing the transport work. Other diffusion preconditioning strategies operate on coarsened space/energy grids and may reduce the overall computational work by drastically reducing the diffusion cost (Roberts and Forget, 2014; Hanuš et al., 2017). The diffusion forms analyzed in this work differ by the treatment of the current and source terms as well as the employed spatial discretization. We begin with a change of iterate notation by denoting the flux after the sweep in Eq. (18) as Uð‘þ1=2Þ and the flux after the diffusion preconditioning step as Uð‘þ1Þ . The diffusion problem used to accelerate the ð‘ þ 1=2Þ transport iterate is given by

  X 0 ð‘þ1=2Þ Dðg‘þ1=2Þ þ rt;g uðg‘þ1=2Þ  rgs ;g ug0 ¼ Qðg‘þ1=2Þ ;

ð21Þ

g 0 2G

where the current term is replaced with Dð‘þ1=2Þ that operates on the diffusion scalar flux. The source term has a functional form dependent on the type of diffusion preconditioning performed (i.e., DSA or NDA). It can be written in a general manner as

  Qðg‘þ1=2Þ ¼ R Q ; Uð‘þ1=2Þ ; Uð‘Þ ;

ð22Þ

where R is a restriction operator that maps information from the high-order transport space onto the low-order diffusion space. Once the diffusion solution of Eq. (21) is computed, the transport solution which we designated as Uð‘þ1Þ can be determined by

  /ðg‘þ1Þ ¼ P /ðg‘þ1=2Þ ; uðg‘þ1=2Þ ;

g 2 G;

ð23Þ

where P is a prolongation operator that maps information from the low-order diffusion space onto the high-order transport space. Eq. (23) constitutes a new fixed-point iteration of Eq. (16) and can be written as





UðG‘þ1Þ ¼ Gdiff Gtrans UðG‘Þ



ð24Þ

to show that both DSA and NDA can be seen to operate in a fixedpoint manner on the transport iterate resulting from Eq. (18). Finally, since Eq. (24) provides a valid fixed-point expression, a corresponding expression for the nonlinear residual can be given by

     ð‘Þ ð‘Þ ð‘Þ ð‘þ1Þ ð‘Þ Fdiff UG ¼ Gdiff Gtrans UG  UG :  UG ¼ UG

ð25Þ

Eq. (25) is therefore utilized in the nonlinear iterative solvers presented in Section 4 for the numerical cases in Section 5 that use diffusion preconditioning. 3.2.1. Diffusion synthetic acceleration We first provide details for the utilized DSA methodology which is similar to Bruss et al. (2014), except they used a linear S2 SA scheme to accelerate the nonlinear CSZ equations. These DSA methods solve for an approximation of the error between the (yet unknown) transport solution and the current iterate. Using a Fick’s law relation between the current and scalar flux, the diffusion term of Eq. (21) is given as

~; ~  Dg r Dg;DSA ¼ r

ð26Þ

where Dg is the diffusion coefficient of group g and is computed from the material properties of the problem: Dg ¼ 1=3rt;g . Because DSA methods solve for the transport error, the source term of Eq. (22) becomes a residual of the scattering source and is given by

Qðg‘þ1=2Þ ¼

X





Þ ð‘Þ rgs ;g RDSA /ðg‘þ1=2  /g0 ; 0 0

ð27Þ

g 0 2G

where the DSA restriction operator, RDSA , will be described shortly. Once a solution of Eq. (21) is obtained, the low-order DSA error is prolongated back onto the transport solution by an additive correction:

  /ðg‘þ1Þ ¼ /ðg‘þ1=2Þ þ P DSA uðg‘þ1=2Þ :

ð28Þ

In this work, we analyze two broad classes of DSA discretizations: cell-based DSA (C-DSA) schemes and vertex-based DSA (VDSA) schemes. Depending on the employed discretization strategy, the V-DSA could be considered a discrete (support at each vertex) or nodal (support at each degree of freedom) scheme. For both AWDD and LD-STZ, the employed C-DSA scheme is from Kim and DeHart (2011), where the nonlinear terms associated with NDA are simply removed. The discretization is a finite volume method in space with a constant flux representation, thereby corresponding to the average flux within a cell. Therefore, the restriction and prolongation operators are 1-to-1 for the AWDD method, but ^ g;K;a , for the only comprise the average expansion coefficient, / LD-STZ method. The work of Lee and Joo (2014) provides details on the restriction/prolongation mapping between DFEM transport discretizations and the cell-based diffusion acceleration methods we use in this work. Due to the known instability of C-DSA schemes (Reed, 1971), we also analyze V-DSA methods. Different schemes are used for the AWDD and LD-STZ methods. For the AWDD method, we employ the 2D equivalent of the reduced-dimensional scheme of Bando et al. (1985). It is a straightforward discretization scheme with a 2d þ 1 stencil, and the restriction and prolongation operators map cell-wise quantities to adjoining vertices (and vice versa). For the LD-STZ method, we use the modified interior penalty (MIP) scheme (Wang and Ragusa, 2010; Turcksin and Ragusa, 2014), developed for DFEM transport discretizations (which includes LD). Since LD-STZ and MIP use the same spatial bases, the restriction and prolongation operators are 1-to-1. Finally, we again state

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

6

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

that further details for the employed C-DSA and V-DSA methods are provided in Appendix A.

where k  k2 is the L2 -norm of the vector. Using Eq. (34) and the Sherman-Morrison formula, the Broyden update is given by

3.2.2. Nonlinear diffusion acceleration Next, we provide details for the nonlinear NDA methods which solve for an approximation of the solution at the current transport iterate. To ensure equivalence between the transport and diffusion solutions, the NDA diffusion operator is given by

Dxð‘Þ ¼

  ð‘þ1=2Þ ~D ~  Dg r b ð‘þ1=2Þ ; Dg;NDA ¼ r g

ð29Þ

bg where D is the vector-wise coefficient of a drift term. Unlike Dg , which is computed solely from the material properties of the b g is derived using information from the higher-order problem, D transport iterate. It has a functional form of ð‘þ1=2Þ

~ð‘þ1=2Þ þ Dg r ~ uðg‘þ1=2Þ b ð‘þ1=2Þ ¼ J D ; g ð‘þ1=2Þ

ug

ð30Þ

where the diffusion current and flux have been restricted from the ð‘ þ 1=2Þ transport iterate. Because NDA methods solve for the transport solution, the source term of Eq. (22) is the restriction of the transport source onto the diffusion phase-space:

  Qðg‘þ1=2Þ ¼ RNDA Q m;g :

ð31Þ

In a similar manner, following the solution of Eq. (21), the loworder NDA solution is prolongated directly back onto the transport solution by

/ðg‘þ1Þ ¼ P NDA



 uðg‘þ1=2Þ :

4. Nonlinear iterative solvers With the residual and fixed-point operations defined in Section 3, we can now provide some summary details for the iterative solvers used in this work: the Richardson, Broyden, NKA, and JFNK methods. Each of these methods can be written in the general form:

xð‘þ1Þ ¼ xð‘Þ þ Dxð‘Þ ;

ð33Þ

where the differences obviously lie in the determination of the update vector, Dxð‘Þ . The Richardson method is the simplest iterative   scheme with Dxð‘Þ ¼ F xð‘Þ using the typical relaxation parameter of unity. The remaining iterative methods are more complex and utilize a dimensional subspace to compute the update vector.

The next nonlinear iterative scheme is Broyden’s method (Broyden, 1965). It is a quasi-Newton, multisecant method that operates by forming an approximate Jacobian and was used as comparative analysis in Calef et al. (2013) for the k-eigenvalue transport problem. We choose to use the Broyden implementation given by Kelley (1995), which utilizes a subspace of maximum size n. If a space of successive updates is given by fsi g, where si ¼ DxðiÞ is used for notational ease and brevity, then the approximation of the inverse Jacobian for xð‘Þ is

¼

‘1 Y j¼‘nð‘Þþ1

Iþ

sjþ1 sTj ksj k22

:

!

;

ð34Þ

ð35Þ

4.2. Nonlinear Krylov acceleration Another class of multisecant methods is Anderson mixing (Anderson, 1965) with the NKA method being one particular variant (Carlson and Miller, 1998; Calef et al., 2013). The method differs from Broyden’s method by forming an approximate inverse Jacobian directly using two vector spaces: 1) the difference between solution states, v i ¼ xði1Þ  xðiÞ , and 2) the difference     between two nonlinear residuals, wi ¼ F xðiÞ  F xði1Þ . The NKA update is given by

"

Dxð‘Þ ¼

‘ X

#

ð‘Þ

zi

vi

" # ‘ X   ð‘Þ þ F xð ‘ Þ  zi wi ;

i¼‘nþ1

ð36Þ

i¼‘nþ1

ð‘Þ

where the zi terms are computed to minimize

  ‘ X ð ‘ Þ ð ‘ Þ ~ z ¼ minn F x  yi wi y2R i¼‘nþ1

ð37Þ

2

~T w ~ . Unlike using a Cholesky decomposition of the reduced matrix: w Broyden, NKA utilizes a reduced memory approach where the oldest ~ are removed when the vector spaces reach the vectors in ~ v and w maximum allowed subspace size. Because NKA uses two spaces to form its update, both the total number of global vectors and operations are 2n þ Oð1Þ. 4.3. Jacobian-free Newton-Krylov The final nonlinear iterative method employed in this work is the JFNK method (Knoll and Keyes, 2004; Knoll et al., 2011). It is an inexact Newton’s method where the inversion of the Jacobian   at iteration ‘; J xð‘Þ , is performed to arbitrary precision using a Krylov method (typically GMRES). This means that JFNK consists of ‘outer’ Newton steps requiring nested ‘inner’ Krylov iterations. In the JFNK method, the standard Newton iterate, given by

    J xð‘Þ Dxð‘Þ ¼ F xð‘Þ ;

4.1. Broyden’s method

J 1 ‘

2 ð‘Þ 1  sT‘1 J 1 ‘ Fðx Þ=ks‘1 k2

The employed Broyden implementation has similarities to the restarted GMRES method. Using a maximum subspace of size n, both methods utilize n þ Oð1Þ global vectors which get cleared and restarted every n iterations. This leads to the variable subspace size, nð‘Þ, in Eq. (34). Also, both methods utilize n þ Oð1Þ global operations (e.g., dot product) to generate the update vector. While small compared to a transport sweep, the global operations can become more noticeably time consuming for large parallel problems.

ð32Þ

In this work, we analyze both NDA and pNDA acceleration schemes where the methodologies are taken from Kim and DeHart (2011). They are the nonlinear analogues to the C-DSA scheme already described (with the same restriction and prolongation operators). Further details are provided in Appendix A.

 ð‘Þ  J 1 ‘ F x

ð38Þ

is solved without forming the Jacobian. Instead, only the action of the Jacobian on a vector is required. For the ‘-th Newton iteration with a solution of xð‘Þ , the action of the Jacobian on a vector v is approximated by a first-order Taylor series expansion,

     ð‘Þ  F xð‘Þ þ ‘ v  F xð‘Þ J x v ; ‘



ð39Þ

where  is the perturbation parameter. From Knoll and Keyes (2004), we use the simple choice of ‘ :

‘ ¼

pffiffiffiffiffiffiffiffiffiffiffi X  mach N  ð‘Þ 1þ j xi j ; N jjv jj2 i¼1

ð40Þ

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

7

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

where N is the system dimension and

mach

indicates the machine

precision limit (generally on the order of 1016 ). While not arising in the numerical experiments of this work, if the evaluation of   F xð‘Þ is precision limited, then another effective formula (Pernice pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and Walker, 1998) is ‘ ¼ mach ð1þ jj xð‘Þ jjÞ= jj v jj. Solving Eq. (38) tightly leads to quadratic convergence of the JFNK method in a neighborhood containing the solution. However, this can leading to ‘‘over-solving” in the early Newton iterations and be detrimental to performance (Knoll and Keyes, 2004). Therefore, we set the convergence criteria of each GMRES inner solve proportional to the current Newton residual,

      jj J xð‘Þ Dxð‘Þ þ F xð‘Þ jj2 6 g jj F xð‘Þ jj2 ;

ð41Þ

void, and scattering regions. The geometry and material properties are given in Fig. 1 and Table 2, respectively. Each boundary is a vacuum condition, and a S16 level-symmetric quadrature is used. Negative fluxes arise in both of the absorbing regions (i.e., Region (III)), though they eventually disappear on impracticably fine meshes. Example solutions are presented on a log scale in Figs. 2 and 3 for the AWDD and LD-STZ methods, respectively, to demonstrate the nonlinear discretization methods’ ability to preserve physical, strictly-positive flux solutions. The AWDD method is plotted against the WDD method with the standard diamond closure, and the LD-STZ method is plotted against the LD method with no flux fixup. White space indicates a mesh cell with negative average solutions. The known highly-oscillatory behavior of the

where the forcing term, g, is some value less than unity. For this work, g was kept constant, but techniques exist to adaptively select g to improve efficiency and the overall robustness of the JFNK iterative process (Pawlowski et al., 2006). Finally, to also improve performance, we initialize the JFNK solver by first performing a single Richardson iteration. This is accounted for when presenting the iteration results in Section 5. 4.4. Summary of Methods In this section, we have presented four nonlinear iterative schemes that utilize a Newton-like update methodology as given in Eq. (33), where each differs in their means of generating the update vectors. Table 1 provides a summary of each method’s mode in handling its subspace vectors along with the total number of possible vectors and global operations for a given space of size n. Broyden and JFNK (its linear GMRES solver) are similar in using a restarted methodology and require n þ Oð1Þ vectors and operations. NKA differs by using a reduced memory approach as well as requiring 2n þ Oð1Þ vectors and operations. The simple Richardson method is included in Table 1 for comparison.

Fig. 1. Configuration of the 2D shielding problem where (X) denotes the region number.

Table 2 Material properties for the 2D shielding problem. Region

rt

rs

Q

(I) (II) (III) (IV)

1.0 0.001 100.0 10.0

1.0 0.0 0.0 9.999

1.0 0.0 0.0 0.0

5. Numerical results In this section, we present a pair of computational examples to analyze the performance of all iterative methods defined in this work. All results were generated using a prototyping code written in MATLAB. A method listed as ‘None’ corresponds to no physicsbased preconditioning being applied. The AWDD method was used on strictly Cartesian grids and the LD-STZ method was used on polygonal grids. All of the polygonal grids were generated in a two-step procedure. First, a polygonal grid using Centroidal Voronoi Tessellation (CVT) is generated with the PolyMesher software (Talischi et al., 2012). Then, the material boundaries of the problem were preserved by adding cut lines through the domain and splitting intersected mesh cells. Simulations are terminated when   kFðxÞk2 6 108 kF xð0Þ k2 , where the initial solution, xð0Þ , is simply an all-zero vector. 5.1. Monoenergetic shielding-type problem The first numerical problem we investigate is a 1-group, heterogeneous 2D shielding-type problem with source, absorbing, nearTable 1 Storage mode and number of global vectors and maximum number of global all reduce operations (per iteration) for a given space of size n. Method

Mode

# Vectors

Operations

Richardson Broyden NKA JFNK

N/A Restart Reduce Restart

1 nþ2 2n þ 3 nþ3

0 nþ1 2n þ 1 nþ1

Fig. 2. Example solutions plotted on a log scale of the WDD method on a Cartesian grid using the (a) diamond closure and the (b) AWDD closure. White space indicates a cell with a negative solution.

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

8

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

to converge the Richardson method using all of the different preconditioning strategies on some coarse meshes. The AWDD method was analyzed on both a 10  5 and 20  10 Cartesian mesh, and the LD-STZ method was analyzed on both a 64 and 256 cell (pre-split) polygonal mesh. We can see how the C-DSA and NDA acceleration methods are unstable and diverge on meshes with cells having large optical thicknesses. Concurrently, the VDSA and pNDA schemes are stable. Next, the C-DSA and NDA acceleration methods are analyzed on the coarsest meshes (i.e., 10  5 Cartesian grid for AWDD and 64 cell polygonal mesh for LD-STZ) using the Broyden, NKA, and JFNK iterative solvers. A subspace of size 10 was consistently used. The number of required transport sweeps are given in Table 4, and it is easy to see that the solution no longer diverges for all cases. Also, the number of sweeps are comparable to those in Table 3 using the V-DSA and pNDA schemes. Next, we analyze the sensitivity of the performance of the JFNK method with varying values of the forcing term g (see Eq. (41)) and Krylov subspace size. It was mentioned that using a looser forcing term prevents ‘‘true” quadratic convergence of the Newton iterations. However, it can lead to an overall computational savings by not over-converging the inner GMRES iterations, especially in the early Newton iterations. Values of 101 ; 102 ; 103 ; 104 , and

Fig. 3. Example average solutions plotted on a log scale of the LD method on a polygonal grid using (a) no fixup method and the (b) STZ fixup method. White space indicates a cell with a negative average solution.

diamond closure (Petrovic´ and Haghighat, 1996) is evident in Fig. 2a. On a final note, it can clearly be seen from Fig. 3 how the polygonal mesh cells were split along cut lines to preserve the material boundaries. There are several iterative analyses to perform and present in this example, listed as the following. 1. Analyze the known instability of the Richardson method with cell-based DSA and NDA schemes; concurrently, show how the nonlinear iterative methods can stabilize the C-DSA and NDA preconditioning schemes analogously to linear Krylov methods preconditioned with DSA (Warsa et al., 2004). 2. Analyze the sensitivity of the performance of the JFNK method with varying values of the forcing term and Krylov subspace size. 3. Analyze the performance of the possible iterative strategies including all combinations of solvers and preconditioning methodologies. We first analyze the known instability of the Richardson method with cell-based DSA and NDA schemes on unresolved meshes. Table 3 provides the number of transport sweeps needed Table 3 Number of sweeps needed for the Richardson method for both discretizations using all preconditioning strategies on coarse meshes. AWDD Method None C-DSA V-DSA NDA pNDA ⁄

Solution diverged.

LD-STZ

10  5

20  10

64

256

138 Div⁄ 26 Div 17

151 Div 22 26 20

165 Div 25 Div 30

165 Div 23 19 22

105 are used for g along with subspace sizes of 5, 10, 20, and 30 for the inner GMRES iterations. The numbers of required Newton iterations and sweeps for the JFNK method without any diffusion preconditioning are provided in Tables 5 and 6 for the AWDD and LD-STZ methods, respectively. We can see for this problem that decreasing values of g leads to increasing sweep counts, no matter the Krylov subspace size due to the plateauing of the required number of Newton iterations. Therefore, a value of 101 for g is used for the remainder of the analysis in this problem. We conclude this numerical example by analyzing the solver performance of all possible iterative combinations including solver method, preconditioning scheme, and varying subspace size on a model that does not suffer any divergent behavior. This corresponds to analyzing the AWDD method on a 40  20 mesh and the LD-STZ method on a 1,024 cell (pre-split) polygonal mesh. For the Broyden, NKA, and JFNK methods, subspace sizes of 5, 10, 20, and 30 were utilized. The number of required sweeps are given in Tables 7 and 8 for the AWDD and LD-STZ methods, respectively. Selected convergence histories (as denoted in the tables) are plotted in Fig. 4. There Table 4 Number of sweeps needed for the Broyden(10), NKA(10), and JFNK(10) methods with the C-DSA and NDA preconditioning schemes on the coarsest analyzed meshes. AWDD Method C-DSA NDA

LD-STZ

Broy.

NKA

JFNK

Broy.

NKA

JFNK

21 17

16 15

19 22

27 20

22 19

23 33

Table 5 Number of Newton iterations and sweeps needed using JFNK without any diffusion preconditioning with varying subspace size and forcing term for the AWDD discretization on a 10  5 Cartesian mesh. X(Y) denotes X Newton iterations and Y transport sweeps.

g Subspace 5 10 20 30

10 7 7 7 7

1

(28) (28) (28) (28)

10 5 5 5 5

2

(33) (31) (31) (31)

103 4 4 4 4

(40) (35) (35) (35)

104 4 4 4 4

(51) (47) (47) (47)

105 4 4 4 4

(72) (69) (55) (55)

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

9

Table 6 Number of Newton iterations and sweeps needed using JFNK without any diffusion preconditioning with varying subspace size and forcing term for the LD-STZ discretization on a 64 cell (pre-split) polygonal mesh. X(Y) denotes X Newton iterations and Y transport sweeps.

g Subspace 5 10 20 30

101 7 7 7 7

(33) (33) (33) (33)

102 5 5 5 5

(40) (38) (38) (38)

103 4 4 4 4

105

104

(44) (39) (39) (39)

4 4 4 4

(58) (52) (50) (50)

4 4 4 4

(98) (88) (76) (76)

Table 7 Number of sweeps needed for the AWDD discretization on a 40  20 Cartesian mesh for the various solution methods. Method

a

Subspace

None

C-DSA

V-DSA

NDA

Richardson

–

156

28

20

19a

pNDA 22

Broyden

5 10 20 30

34 31 31 32

13 12 12 12

11a 11 11 11

11 11 11 11

12 11 11 11

NKA

5 10 20 30

29 27 24 24

12 12 12 12

10a 10 10 10

11 11 11 11

11 11 11 11

JFNK

5 10 20 30

29 29 29 29

13 13 13 13

11 11 11 11

10a 10 10 10

14 14 14 14

Error plotted in Fig. 4a.

Table 8 Number of sweeps needed for the LD-STZ discretization on a 1,024 cell (pre-split) polygonal mesh for the various solution methods.

a

Method

Subspace

None

C-DSA

V-DSA

NDA

Richardson

–

165

17a

22

18

pNDA 21

Broyden

5 10 20 30

35 36 32 30

13 12 12 12

12a 13 13 13

12 12 12 12

13 12 12 12

NKA

5 10 20 30

38 27 25 25

12 12 12 12

12 12 12 12

12 11a 11 11

12 12 12 12

JFNK

5 10 20 30

27 27 27 27

12a 12 12 12

13 13 13 13

13 12 12 12

13 13 13 13

Error plotted in Fig. 4b.

are several iterative behaviors to note. First, the Richardson method with diffusion acceleration and the Broyden, NKA, and JFNK methods without diffusion acceleration provide comparable reduction in the number of sweeps. Additional reduction in the number of required transport sweeps is realized by using the Broyden, NKA, and JFNK methods with any diffusion acceleration scheme, and it is not affected by increasing subspace size. This shows that the necessary iterative work is dominated by the diffusion operators for this problem. Finally, the Broyden and NKA methods without diffusion acceleration have reductions in the required transport sweeps with increasing subspace size (at the expense of memory). This is not true for JFNK because the number of sweeps per Newton iteration is already small by using g ¼ 101 .

Fig. 4. Normalized L2 -norm of the residual error of select solver permutations as identified in Tables 7 and 8.

5.2. Multigroup detector problem In this second numerical problem, we analyze the solution algorithms on a heterogeneous, multigroup detector problem. The materials for this problem include 252Cf, graphite, 303 stainless steel, borated water, and BF3. This problem primarily differs from the previous one by being a multigroup problem where the diffusive material (graphite) has significant thermal neutron scattering. The geometry and component materials are given in Fig. 5 and Table 9, respectively. The source resides in the bottom-left corner and the detector in the top-right with layers of graphite, steel, and borated water between them. The bottom and left boundaries are reflective, and the top and right boundaries are vacuum. This problem uses 12 energy groups (7 fast and 5 thermal) with the cross sections generated at room temperature. The graphite material leads to strong thermal neutron upscattering in the 5

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

10

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx Table 10 Total interaction rate (rxn=cm3  s) of the detector region of the different discretization methods with increasing mesh refinement.

Table 9 Component materials for the multigroup detector problem along with their colors used in Fig. 5.

Source Moderator Structural Absorber Detector

Material 252

Cf graphite 303-SS Borated H2O BF3

Diamond

AWDD

LD

LD-STZ

1 2 3 4 5 6

1.476e9 1.023e7 2.624e5 6.544e2 7.976e2 9.822e2

6.997e3 1.622e3 1.109e3 1.051e3 1.040e3 1.039e3

1.726e3 6.860e2 1.002e3 1.026e3 1.032e3 1.039e3

1.935e3 6.844e2 9.945e2 1.026e3 1.032e3 1.039e3

for the coarsest spatial grids (a 12  12 Cartesian mesh and a 64 cell (pre-split) polygonal mesh for the AWDD and LD-STZ methods, respectively). Only iterative results from the thermal group set are presented since the focus of this problem is the solution strategies in a multigroup setting. Tables 11 and 12 give the required number of sweeps for the thermal group set for the AWDD and LD-STZ methods, respectively. Selected convergence histories (as denoted in the tables) are then plotted in Fig. 6. For both discretizations, the C-DSA and NDA acceleration schemes are unstable with the Richardson iterative method but stabilized using the Broyden, NKA, and JFNK methods. Also, we again see that the Richardson method with stable diffusion acceleration (V-DSA and pNDA) as well as unaccelerated Broyden/NKA/JFNK provide a significant reduction in the required transport sweeps. Like the monoenergetic problem, AWDD solved with Broyden/NKA/JFNK and any

Fig. 5. Configuration of the 2D detector problem.

Component

Mesh Level

Table 11 Number of sweeps needed for the thermal groups subset for the AWDD discretization on a 12  12 Cartesian mesh for the various solution methods.

Color Method

Pink Yellow Blue White Green

thermal groups. No upscattering occurs into any of the 7 fast groups. The 252Cf source emits neutrons with the Watt spectrum (Fröhner, 1990), and again, a S16 level-symmetric quadrature is used. A specific iterative strategy is employed to minimize the overall computational work. The energy groups are split into 8 subsets. Each of the seven fast groups are in their own subset, and the 5 thermal groups are collected into a subset. An outer iteration consisting of a Gauss-Seidel sweep through the group subsets is performed. Within each subset, Jacobi iterations are performed using the prescribed iterative strategy. Using this strategy, only 1 outer iteration through the group subsets is needed since all of the upscattering is accounted for in the thermal groups’ Jacobi iterations. First, we provide the total interaction rate of the detector region in Table 10 for the different discretization schemes under increasing mesh refinement. The coarsest meshes are a 12  12 Cartesian mesh and a 64 cell (pre-split) polygonal mesh for the WDD and LD methods, respectively. The AWDD, LD, and LD-STZ schemes all converge to the same value after a few iterations and always remain positive. We note that for the coarser meshes the LD method does yield mesh cells with negative reaction rates, just not in the detector region. On the other hand, the diamond-difference method on the coarser meshes yields incredibly inaccurate reaction rates (even becoming negative), again demonstrating its known oscillatory behavior (Petrovic´ and Haghighat, 1996). Next, we provide some iterative convergence analysis for this problem. Just like the monoenergetic problem, some low-fidelity analysis was performed on the forcing factor for the JFNK method, with g ¼ 101 again selected. For this problem, we present results

a b

Subspace

None

C-DSA

V-DSA

NDA

pNDA

Richardson

–

465

Diva

46b

Div

56

Broyden

5 10 20 30

73 75 60 56

28 27 25 25

25 24 24 24

21 20 19 19

18 19 16b 16

NKA

5 10 20 30

65 40 36 35

25 21 20 20

22 21 20 20

18 16 16 16

17 16b 18 18

JFNK

5 10 20 30

63 57 57 57

32 30 29 29

29 30 29 29

22 22 22 22

20b 20 20 20

Solution diverged. Error plotted in Fig. 6a.

Table 12 Number of sweeps needed for the thermal groups subset for the LD-STZ discretization on a 64 cell (pre-split) polygonal mesh for the various solution methods. Method

a b

Subspace

None

C-DSA

V-DSA

NDA

pNDA

Richardson

–

539

Diva

17b

Div

54

Broyden

5 10 20 30

90 73 59 51

46 41 38 38

12b 12 12 12

58 55 54 47

40 40 37 36

NKA

5 10 20 30

70 43 36 34

38 35 34 34

12b 12 12 12

53 37 34 31

41 42 42 39

JFNK

5 10 20 30

69 52 52 52

49 39 39 39

13b 13 13 13

53 49 49 49

30 32 32 32

Solution diverged. Error plotted in Fig. 6b.

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

11

selected nonlinear discretizations used and were analyzed on Cartesian and polygonal meshes, respectively. Due to the transport operator being nonlinear from these nonlinear discretization methods, powerful linear solvers such as GMRES and BiCGStab are unusable. Instead, different fixed-point and Newton-like iterative schemes were utilized, including the Richardson, Broyden, NKA, and JFNK methods. These methods were further preconditioned with physics-based diffusion acceleration: linear diffusion synthetic acceleration and nonlinear diffusion acceleration. Two variants of DSA were considered: cell-based DSA and vertex-based DSA. Also, an under-relaxing variant of NDA, the partial-current NDA method, was analyzed. Numerical results demonstrated the necessary considerations for desirable solution methodologies. The C-DSA and NDA schemes accelerating the Richardson method were demonstrated to be unstable on meshes with optically thick cells, but stabilized with the Broyden, NKA, and JFNK methods. Either using Richardson with diffusion acceleration or Broyden, NKA, and JFNK without diffusion acceleration provided appreciable reductions in the required number of transport sweeps, but the best reduction in transport sweeps was realized using Broyden/NKA/JFNK with diffusion acceleration. For the AWDD method, all diffusion acceleration schemes provided similar iterative performance benefits. However, for the LD-STZ, the cellbased DSA and NDA schemes, while still stable using Broyden, NKA, and JFNK, provided little additional iterative benefits for the more optically-diffusive multigroup problem, due to the greater inconsistency between the transport and diffusion discretizations. Instead, the more consistent V-DSA scheme continued providing excellent iterative performance. All of these results are analogous to linear transport operators. Acknowledgments The submitted manuscript has been authored by a contractor of the U.S. Government under contract No. DOE89233018CNR000004. Accordingly, the U.S. Government retains a non-exclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes. Appendix A. Diffusion acceleration operators Here, we provide further details of the employed diffusion operators. Restricting ourselves to the monoenergetic case for simplicity, a low-order diffusion approximation of the transport equation can be written in analytical form as Fig. 6. Normalized L2 -norm of the residual error of select solver permutations as identified in Tables 11 and 12.

diffusion acceleration methodology yields reduced sweep numbers of similar amounts. Conversely, due to the increased diffusivity and optical thicknesses of this problem along with the greater inconsistencies between SN and diffusion, the LD-STZ method using Broyden/NKA/JFNK with cell-based diffusion acceleration provides little additional reduction in the number of required sweeps. However, the greater consistency of the V-DSA scheme leads to a significant reduction in the number of sweeps for all solver methods, thereby demonstrating the need to consider the employed diffusion acceleration strategy analogously to linear transport operators.

~ ~ r J þ rr u ¼ Q;

ð42Þ

where rr is the removal cross section. The fundamental difference between linear and nonlinear diffusion acceleration is how to handle the current, ~ J, and the source term, Q. The analytical form of the current for DSA and NDA are

~ ~ u and ~ ~u þ D b u; J DSA ¼ Dr J NDA ¼ Dr

ð43Þ

b u is a drift term with D b a vector quantity respectively, where D computed to maintain equivalence with the transport operator. The source term, Q, is either the transport scattering residual or source, respectively (Eqs. (27) and (31), respectively), restricted to the P0 angular space.

6. Conclusions

A.1. Cell-based diffusion operators

In this work, we investigated the effectiveness of iterative solution strategies for solving SN transport problems utilizing nonlinear discretizations. The AWDD and LD-STZ fixup methods were the

We first provide discretization details for the cell-based diffusion operators: C-DSA, NDA, and pNDA. Using a constant finite volume approach, the discretized Eq. (42) can be written for cell K as

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

12

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

Xjf j J þ rr;K uK ¼ QK ; j K j f ;K f 2@K

ð44Þ

where j f j and j K j are the measures of face f and cell K, respectively (length and area, respectively, in 2D). Depending on the diffusion method, the current becomes

  ef u 0  u J DSA;f ;K ¼  D f ;K f ;K ;     ef u 0  u b J NDA;f ;K ¼  D f ;K f ;K þ D f uf ;K 0 þ uf ;K ;     ef u 0  u bþ b J pNDA;f ;K ¼  D f ;K f ;K  D f uf ;K 0  D f uf ;K ;

ð45aÞ ð45bÞ ð45cÞ

where uf ;K 0 denotes the cell flux connected to cell K through face f,

ef ¼ D

8
Df ;K Df ;K 0

f ;K Lf ;K 0 þDf ;K 0 Lf ;K

bDf ;K bLf ;K þDf ;K

; f 2 Cint

ð46Þ

f 2 Cinc [ Cref

;

and where Lf ;K is the distance between the centroids of face f and cell K and b is an albedo term (b ¼ 0 for f 2 Cref and b ¼ 0:5 for f 2 Cinc ). For the NDA method, the equivalence diffusion coefficient, b f , in Eq. (45b) can be directly computed with J D being the

A.2. Vertex-based finite-difference operator The next diffusion discretization we present is the V-DSA scheme used for the AWDD method in this work and is taken from Bando et al. (1985). On a 2D Cartesian grid where cells are indexed at full integer values (i; j) and vertices are half-range indices (i þ 1=2; j þ 1=2), the equation for the (i þ 1=2; j þ 1=2) vertex is given by



 Diþ1;jþ1=2  uiþ3=2;jþ1=2  uiþ1=2;jþ1=2 hiþ1  Di;jþ1=2  þ uiþ1=2;jþ1=2  ui1=2;jþ1=2 hi  Diþ1=2;jþ1  uiþ1=2;jþ3=2  uiþ1=2;jþ1=2  hjþ1  Diþ1=2;j  uiþ1=2;jþ1=2  uiþ1=2;j1=2 þ hj ^ r;iþ1=2;jþ1=2 uiþ1=2;jþ1=2 þr

NDA;f ;K

restriction of the transport net current onto face f. In a similar manb  , can be comner, the pNDA equivalence diffusion coefficients, D f puted with the transport restricted partial currents:





1e b þu ; D f uf ;K 0  uf ;K þ D f ;K f 2  1e  b u 0 : ¼ D f uf ;K 0  uf ;K þ D f ;K f 2

J þf ;K ¼ 

iþ1 X 1X V i0 ;j0 Q i0 ;j0 4 0 0 jþ1

¼

where

Diþ1=2;jþ1 ¼

ð‘þ1=2Þ

ð48Þ

m¼1

and

( ð‘þ1Þ /K

¼

ð‘þ1=2Þ

/K

þ uK

uðK‘þ1=2Þ ;

ð‘þ1=2Þ

; C  DSA ðpÞNDA

ð49Þ

respectively. LD-STZ differs by mapping the average transport flux ^ K;a ) to the diffusion flux. Therefore, the cell-based difin cell K (i.e., / fusion source and prolongation for the ‘-th LD-STZ transport sweep are

ð‘þ1=2Þ

QK

8   ^ ð‘þ1=2Þ  / ^ ð‘Þ ; C  DSA > > K;a < rs;K /K;a ¼ M P > > wm Q m;K ; ðpÞNDA :

1X h 0D 0; 2 0 j i;j

ð53bÞ

iþ1 1X h 0D 0 ; 2 0 i i ;jþ1

ð53cÞ

iþ1 1X h 0D 0 ; 2 0 i i ;j

ð53dÞ

iþ1 X 1X rr;i0 ;j0 V i0 ;j0 ; 4 0 0

ð53eÞ

jþ1

It is easy to see that Eq. (45c) was derived by using  JpNDA;f ;K ¼ J þ f ;K  J f ;K . The boundary conditions are naturally tied to the face currents. If face f is on a reflective boundary, then Jf ;K ¼ 0. Conversely, incident boundary conditions are already specified by Eqs. (45) and (46), where uf ;K 0 ¼ 0. The restriction and prolongation operations differ between AWDD and LD-STZ. For AWDD, the mappings are 1-to-1. Therefore, the cell-based diffusion source and prolongation for the ‘-th AWDD transport sweep are

QK

ð53aÞ

j ¼j

ð47aÞ ð47bÞ

  8 ð‘þ1=2Þ ð‘Þ > r  /K ; C  DSA s;K /K > < M ¼ X > > wm Q m;K ; ðpÞNDA :

1X h 0D 0; 2 0 j iþ1;j jþ1

Diþ1;jþ1=2 ¼ Di;jþ1=2 ¼

J f ;K

ð52Þ

i ¼i j ¼j

j ¼j

i ¼i

Diþ1=2;j ¼

i ¼i

r^ r;iþ1=2;jþ1=2 ¼

jþ1

i ¼i j ¼j

and hi ¼ hiþ1=2  hi1=2 and V i;j ¼ hi hj . The first four terms of Eq. (52) represent the current contributions to vertex (i þ 1=2; j þ 1=2) from adjoining faces. Therefore, if an adjoining boundary face to a vertex is reflective, then its current contribution is 0. This is analogous to setting the corresponding D to 0. For an incident boundary condition, the corresponding diffusion term in Eq. (52) becomes: D=h ! L=2, where L is the length of the face in 2D. The spatial bases are obviously different between the cell-based AWDD transport solution and the vertex-based diffusion solution. Therefore, mapping occurs between a cell and it’s adjacent vertices. The source term, Q i;j , in Eq. (52) is already a cell-based quantity and is simply ð‘þ1=2Þ

Q i;j

  ð‘þ1=2Þ ð‘Þ ¼ rs;i;j /i;j  /i;j ;

ð54Þ

for the ‘-th AWDD transport sweep. Once uð‘þ1=2Þ is computed, the DSA update for cell (i; j) is

ð50Þ ð‘þ1Þ

/i;j

m¼1

ð‘þ1=2Þ

¼ /i;j

iþ1 X 1X Þ uði0‘þ1=2 : 1=2;j0 1=2 4 0 0 jþ1

þ

ð55Þ

i ¼i j ¼j

and

( ^ ð‘þ1Þ ¼ / K;a

^ ð‘þ1=2Þ þ uð‘þ1=2Þ ; C  DSA / K;a K

u

respectively.

ð‘þ1=2Þ ; K

ðpÞNDA

ð51Þ

A.3. Modified interior penalty DFEM form The final diffusion discretization we present is the MIP method, which is the V-DSA scheme used for LD-STZ method in this work

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

(Turcksin and Ragusa, 2014). Defining a test function v, the MIP weak form is given by

XR

dr

K2Th K

þ

rr uv þ

XR

f 2Cint f

þ

XR

K2Th K

~u  r ~v dr Dr

  ds jf sutsv t þ fD@ n ugsv t þ sutfD@ n v g ds



f 2Cinc f

¼

XR

XR

K2Th K

jf uv  12 v D@ n u  12 uD@ n v



ð56Þ

dr Qv ;

~ and the jump and average terms of the discontinwhere @ n ¼ ~ nr uous quantities across an interior face are given by

sv t ¼ v þ  v 

and fv g ¼

1 þ ðv þ v  Þ; 2

ð57Þ

respectively. The left/right values along a face f are given by the trace: v  ¼ lims!0 v ð~ r þ s~ nÞ in which the face normal vector is arbitrarily chosen in a direction across f (it is required to be the outgoing normal for boundary faces). The penalty term is given by   jf ¼ max 14 ; jIPf , where

jIPf ¼

8

 Dþ D > > < 2c hþf þ hf f 2 Cint ; f

> D > : c hf f

f

ð58Þ

f 2 Cinc ;

c is a positive constant taken to be 8 in this work (Turcksin and 

Ragusa, 2014) and h is the length of the cell in the  direction orthogonal to face f. It is easy to see how the MIP boundary conditions are handled from Eq. (56). Because we utilize the same spatial mesh for the transport and diffusion operators and LD and LD-STZ have the same spatial basis, the restriction and prolongation operations using MIP are 1-to-1. This means that the diffusion source and DSA update of cell K for the ‘-th transport sweep are succinctly given by

  ð‘þ1=2Þ ð‘Þ ~ ð‘þ1=2Þ ¼ rs;K ~ /K ~ /K ; Q K

ð59Þ

and ð‘þ1Þ ð‘þ1=2Þ ~ ~ ðK‘þ1=2Þ ; /K ¼~ /K þu

ð60Þ

respectively. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, athttps://doi.org/10.1016/j.anucene.2019. 107080. References Adams, M.L., Larsen, E.W., 2002. Fast iterative methods for discrete-ordinates particle transport calculations. Prog. Nucl. Energy 40 (1), 3–159. Alcouffe, R.E., 1977. Diffusion synthetic acceleration methods for the diamonddifferenced discrete-ordinates equations. Nucl. Sci. Eng. 64 (2), 344–355. Alcouffe, R.E., 1993. An adaptive weighted diamond differencing method for threedimensional XYZ geometry. Trans. Am. Nucl. Soc. 68, 206–212. Anderson, D.G., 1965. Iterative procedures for nonlinear integral equations. J. ACM 12 (4), 547–560. Bailey, Teresa S., 2008. The Piecewise Linear Discontinuous Finite Element Method Applied to the RZ and XYZ Transport Equations PhD thesis. A&M University, Texas. Baker, R.S., Koch, K.R., 1998. An Sn algorithm for the massively parallel CM-200 computer. Nucl. Sci. Eng. 128 (3), 312–320. Bando, M., Yamamoto, T., Saito, Y., Takeda, T., 1985. Three-dimensional transport calculation method for eigenvalue problems using diffusion synthetic acceleration. Nucl. Sci. Technol. 22 (10), 841–850.

13

Barbu, A.P., Adams, M.L., 2019. A linear diffusion-acceleration method for keigenvalue transport. In: M&C 2019 – International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering. Broyden, C.G., 1965. A class of methods for solving nonlinear simultaneous equations. Math. Comput. 19 (92), 577–593. Bruss, D.E., Morel, J.E., Ragusa, J.C., 2014. S2 SA preconditiong for the Sn equations with strictly nonnegative spatial discretization. J. Comput. Phys. 273 (15), 706– 719. Calef, M.T., Fichtl, E.D., Warsa, J.S., Berndt, M., Carlson, N.N., 2013. Nonlinear Krylov acceleration applied to a discrete ordinates formulation of the k-eigenvalue problem. J. Comput. Phys. 238 (1), 188–209. Carlson, N.N., Miller, K., 1998. Design and application of a gradient-weighted moving finite element code I: in one dimension. SIAM J. Scientific Comput. 19 (3), 728–765. Cho, N.Z., Lee, G.S., Park, C.J., 2003. On a new acceleration method for 3D whole core transport calculations. In: Proc. Annual Meeting of Atomic Energy Society of Japan. Fichtl, E.D., Warsa, J.S., Densmore, J.D., 2010. The Newton-Krylov method applied to negative-flux fixup in SN transport calculations. Nucl. Sci. Eng. 165 (3), 331–341. Fröhner, F.H., 1990. Evaluation of 252 Cf prompt fission neutron data from 0 to 20 MeV by Watt Spectrum fit. Nucl. Sci. Eng. 106 (3), 345–352. Gill, D.F., Azmy, Y.Y., Warsa, J.S., Densmore, J.D., 2011. Newton‘s method for the computation of k-eigenvalues in SN transport applications. Nucl. Sci. Eng. 168 (1), 37–58. Guthrie, B., Holloway, J.P., Patton, B.W., 1999. GMRES as a multi-step transport sweep accelerator. Transport Theory Stat. Phys. 28 (1), 83–102. Hackemack, M.W., 2018. Fixed-point acceleration of nonlinear spatial discretizations. In: PHYSOR 2018: Reactor Physics Paving the Way towards More Efficient Systems. Hackemack, M.W., Becker, T.L., 2017. Alternate LD closures to ensure cell-averaged positivity on unstructured meshes. In: Transactions of the American Nuclear Society, pp. 706–709. Hackemack, M.W., Ragusa, J.C., 2018. Quadratic serendipity discontinuous finite element discretization for SN transport on arbitrary polygonal grids. J. Comput. Phys. 374 (1), 188–212. Hamilton, S., Benzi, M., Warsa, J., 2009. Negative flux fixups in discontinuous finite element SN transport. In: 2009 International Conference on Mathematics, Computational Methods & Reactor Physics (M&C 2009). Hanuš, M., Ragusa, J.C., Hackemack, M.W., 2017. A study of various thermal upscattering acceleration schemes for massively parallel transport sweeps. In: M&C 2017 – International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering. Hawkins, W.D., Smith, T., Adams, M.P., Rauchwerger, L., Amato, N.M., Adams, M.L., 2012. Efficient massively parallel transport sweeps. Trans. Am. Nucl. Soc. 107, 477–481. Hill, T.R., 1975. ONETRAN: a discrete ordinates finite element code for the solution of the one-dimensional multigroup transport equation (Technical Report LA5990-MS). Los Alamos Scientific Laboratory. Kelley, C.T., 1995. Iterative methods for linear and nonlinear equations. Soc. Ind. Appl. Math. Kim, K.-S., DeHart, M.D., 2011. Unstructured partial- and net-current based coarse mesh finite difference acceleration applied to the extended step characteristics method in NEWT. Ann. Nucl. Energy 38, 527–534. Knoll, D.A., Keyes, D.E., 2004. Jacobian-free Newton-Krylov methods: a survey of approaches and applications. J. Comput. Phys. 193 (2), 357–397. Knoll, D.A., Park, H., Newman, C., 2011. Acceleration of k-eigenvalue/criticality calculations using the Jacobian-free Newton-Krylov method. Nucl. Sci. Eng. 167 (2), 133–140. Knoll, D.A., Park, H., Smith, K.S., 2011. Application of the Jacobian-Free NewtonKrylov method to nonlinear acceleration of transport source iteration in slab geometry. Nucl. Sci. Eng. 167 (2), 122–132. Lathrop, K.D., 1969. Spatial differencing of the transport equation: positivity vs. accuracy. J. Comput. Phys. 4, 475–498. Lee, D.W., Joo, H.G., 2014. Coarse mesh finite difference acceleration of discrete ordinate neutron transport calculation employing discontinuous finite element method. Nucl. Eng. Technol. 46 (6), 783–796. Lewis, E.E., Miller, W.F., 1984. Computational Methods of Neutron Transport. John Wiley and Sons Inc, New York, NY. Maginot, P.G., Morel, J.E., Ragusa, J.C., 2012. A nonnegative momentpreserving spatial discretization scheme for the linearized Boltzmann transport equation in 1-D and 2-D Cartesian geometries. J. Comput. Phys. 231 (20), 6801–6826. Maginot, P.G., Ragusa, J.C., Morel, J.E., 2017. Nonnegative methods for bilinear discontinuous differencing of the SN equations on quadrilaterals. Nucl. Sci. Eng. 185 (1), 53–69. Minor, B., Mathews, K.A., 1995. Exponential characteristic spatial quadrature for discrete ordinates radiation transport with rectangular cells. Nucl. Sci. Eng. 120 (3), 165–186. Patton, B.W., Holloway, J.P., 2002. Application of preconditioned GMRES to the numerical solution of the neutron transport equation. Ann. Nucl. Energy 29 (2), 109–136. Pawlowski, R.P., Shadid, J.N., Simonis, J.P., Walker, H.F., 2006. Globalization techniques for Newton-Krylov methods and applications to the fully coupled solution of the Navier-Stokes equations. SIAM Rev. 48 (4), 700–721. Pernice, M., Walker, H.F., 1998. NITSOL: a Newton solver for nonlinear systems. SIAM J. Scientific Comput. 19 (1), 302–318.

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080

14

M.W. Hackemack / Annals of Nuclear Energy xxx (xxxx) xxx

Petrovic´, B.G., Haghighat, A., 1996. Analysis of inherent oscillations in multidimensional SN solutions of the neutron transport equation. Nucl. Sci. Eng. 124 (1), 31–62. Reed, W.H., 1971. The effectiveness of acceleration techniques for iterative methods in transport theory. Nucl. Sci. Eng. 45 (3), 245–254. Roberts, J.A., Forget, B., 2014. Multigroup diffusion preconditioners for multiplying fixed-source transport problems. J. Comput. Phys. 274 (1), 455–472. Smith, K.S., 1983. Nodal method storage reduction by non-linear iteration. Transactions of the American Nuclear Society, vol. 44, p. 265. Smith, K.S., Rhodes III, J., 2002. Full-core 2-D LWR core calculations with CASMO-4E. In: Proc. Int. Conf. New Frontiers of Nuclear Technology: Reactor Physics, Safety, and High-Performance Computing, pp. 7–10. Suslov, I.R., 2003. An algebraic collapsing acceleration method for acceleration of the inner (scattering) iterations in long characteristics transport theory. In: Proceedings of the International Conference on Supercomputing in Nuclear Applications SNA 2003. Talischi, C., Paulino, G.H., Pereira, A., Menezes, I., 2012. PolyMesher: a generalpurpose mesh generator for polygonal elements written in Matlab. Struct. Multidiscip. Optim. 45 (3), 309–328.

Turcksin, B., Ragusa, J.C., 2014. Discontinuous diffusion synthetic acceleration for S_n) transport on 2D arbitrary polygonal meshes. J. Comput. Phys. 274 (1), 356– 369. Wang, Y., Ragusa, J.C., 2010. Diffusion synthetic acceleration for high-order discontinuous finite element S_n) transport schemes and application to locally refined unstructured meshes. Nucl. Sci. Eng. 166 (2), 145–166. Wareing, T.A., 1997. An exponential discontinuous scheme for discrete-ordinate calculations in Cartesian geometries. In: Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Applications. Wareing, T.A., McGhee, J.M., Morel, J.E., Pautz, S.D., 2001. Discontinuous finite element SN methods on three-dimensional unstructured grids. Nucl. Sci. Eng. 138 (3), 256–268. Warsa, J.S., Wareing, T.A., Morel, J.E., 2004. Krylov iterative methods and the degraded effectiveness of diffusion synthetic acceleration for multidimensional SN calculations in problems with material discontinuities. Nucl. Sci. Eng. 147 (3), 218–248. Willert, J., Taitano, W.T., Knoll, D.A., 2014. Leveraging Anderson Acceleration for improved convergence of iterative solutions to transport systems. J. Comput. Phys. 273 (15), 278–286.

Please cite this article as: M. W. Hackemack, Solving nonlinear discretizations of SN transport calculations, Annals of Nuclear Energy, https://doi.org/ 10.1016/j.anucene.2019.107080