An Accurate Globally Conservative Subdomain Discontinuous Least-squares Scheme for Solving Neutron Transport Problems
AAn Accurate Globally Conservative SubdomainDiscontinuous Least-squares Scheme for SolvingNeutron Transport Problems
Weixiong Zheng a, ∗ , Ryan G. McClarren b , Jim E. Morel b a Nuclear Engineering, University of California, Berkeley, Berkeley, CA 94709 b Nuclear Engineering, Texas A&M University, College Station, TX 77843-3133
Abstract
In this work, we present a subdomain discontinuous least-squares (SDLS) schemefor neutronics problems. Least-squares (LS) methods are known to be inaccu-rate for problems with sharp total-cross section interfaces. In addition, theleast-squares scheme is known not to be globally conservative in heterogeneousproblems. In problems where global conservation is important, e.g. k-eigenvalueproblems, a conservative treatment must be applied. We, in this study, proposean SDLS method that retains global conservation, and, as a result, gives highaccuracy on eigenvalue problems. Such a method resembles the LS formulationin each subdomain without a material interface and differs from LS in that anadditional least-squares interface term appears for each interface. The scalarflux is continuous in each subdomain with continuous finite element method(CFEM) while discontinuous on interfaces for every pair of contiguous subdo-mains. SDLS numerical results are compared with those obtained from othernumerical methods with test problems having material interfaces. High accu-racy of scalar flux in fixed-source problems and k eff in eigenvalue problems aredemonstrated. Keywords:
Global conservation; least-squares; discrete ordinates ∗ Corresponding author
Email addresses: [email protected] (Weixiong Zheng), [email protected] (Ryan G. McClarren), [email protected] (Jim E. Morel)
Preprint submitted to Nuclear Science and Engineering July 8, 2018 a r X i v : . [ phy s i c s . c o m p - ph ] N ov . Introduction Neutral particle transport problems are governed by a first-order, hyperbolicequation. As a result particle transport problems, especially those that arestreaming dominated, can have sharp gradients along the characteristic lines.Such problems require spatial discretizations that allow for discontinuities inspace. Moreover, in regions where there is a strong scattering, discontinuous fi-nite element method (DFEM) have been shown to properly preserve the asymp-totic diffusion limit of the transport equation [1]. These two facts have led todiscontinuous Galerkin finite elements being widely accepted in transport cal-culations. Nevertheless, the discontinuous finite element methods have moredegrees of freedom (DoFs) than their continuous counterparts, especially in 3-D. Another approach to solving transport problems involves forming a second-order transport operator, and solving the resulting equations using continuousfinite elements. Second-order transport problems based on the parity of theequations and the self-adjoint angular flux (SAAF) equation are well-known[2,3, 4, 5, 6]. The resulting equations are symmetric, but are ill-posed in voidregions and ill-conditioned in near voids.More recently, Hansen, et al. [7] derived a second-order form that is equiva-lent to minimizing the squared residual of the transport operator; it is thereforecalled the least-squares transport equation (LSTE). This method can also beformed by multiplying the transport equation by the adjoint transport operatoror by applying least-squares finite elements in space to the transport equation.The resulting equations are well-posed in void and are symmetric positive def-inite (SPD). The least-squares method has been continuously investigated inthe applied mathematics and computational fluid dynamics communities fordecades [8, 9, 10, 11, 12, 13]. However, the method generally has low accuracyin problems with interfaces between optically thin and thick materials withoutlocal refinement near the interface[14, 16, 15, 17]. Additionally, the standardleast-squares method does not have particle conservation (except in the limit2s the number of zones goes to infinity), which would induce undesirably large k -eigenvalue errors in neutronics simulations [18], unless conservative accelera-tion schemes, such as nonlinear diffusion acceleration described in [19, 20, 21],are applied to regain the conservation.In the process of deriving the least-squares equations, the symmetrizationof the transport equation converts a hyperbolic equation to an elliptic equationand causes loss of causality. For example, the presence of a strong absorberdownstream of a source can influence the solution upstream of the source. Inremedying this deficiency, it is desirable to add back some asymmetry to theoperator, but do so in such a way that still preserves the positive properties ofthe least-squares method.To this end we develop a least-squares method that minimizes the square ofthe transport residual over certain regions of the entire problem domain. In eachof these subregions the interaction cross-sections are slowly varying functions ofspace. We allow the solution between these regions to be discontinuous. The re-sulting scheme essentially solves least-squares independently in each subregion,and the regions are coupled through a sweep-like procedure . Furthermore, withthe correct weighting in the least-squares procedure, we can construct a methodthat is globally conservative.The reminder of the paper is organized as follows: we start off reviewing thegoverning equation and the ordinary least-squares finite element weak formula-tion derived from the minimization point of view in Sec. 2. In Sec. 3, we propose anew method that is based on a least-squares formulation over contiguous blocksin the problem base on a novel subdomain-discontinuous functional. We alsoderive the corresponding weak formulations in this section. Next, we furthertheoretically demonstrate the method, unlike standard least-squares methods, In this work, we will utilize the discrete-ordinates method for angular discretization.When solving the discrete-ordinates equations with first-order discretization techniques, suchas discontinuous Galerkin method, the procedure is such that the transport equation is solvedfrom upstream cells to downstream cells, a procedure called a “sweep”.
2. Revisit least-squares discretization of transport equation
We will consider steady, energy-independent transport problems with isotropicscattering in this work. The complications of energy dependence and anisotropicscattering can be readily incorporated into our method. The steady transportequation used to describe neutral particles of a single speed is expressed inoperator form as: Lψ = q s , (1a)where the transport operator L is defined as the sum of the streaming operator Ω · ∇ ( · ) and total collision operator σ t : L = Ω · ∇ ( · ) + σ t , (1b)and q s represents the total volumetric source defined as the sum of the scatteringsource Sψ and a fixed volumetric source q ( (cid:126)r, Ω): q s = Sψ + q. (1c)In these equations ψ ( (cid:126)r, Ω) is the angular flux of neutral particles with units ofparticles per unit area per unit time, where (cid:126)r ∈ D ⊂ R , D stands for problemdomain and Ω ∈ S is a point on the unit sphere representing a direction oftravel for the particles. S is the scattering operator. For isotropic scattering, itis defined as Sψ = (cid:90) π d Ω (cid:48) σ s (cid:16) Ω · Ω (cid:48) (cid:17) ψ ( (cid:126)r, Ω (cid:48) ) = σ s φ π , (2)4here σ s (cid:16) Ω · Ω (cid:48) (cid:17) is differential scattering cross section defined as σ s (cid:16) Ω · Ω (cid:48) (cid:17) = σ s π (3)and φ is the scalar flux defined as φ = (cid:90) π d Ω ψ. (4)Additionally, σ s is the scattering cross section.The boundary conditions for Eq. (1) specify the angular flux ψ inc on theboundary for incoming directions: ψ ( (cid:126)r, Ω) = ψ inc ( (cid:126)r, Ω) for r ∈ ∂ D , (cid:126)n · Ω < . (5)We discretize the angular component of the transport equation using discrete-ordinates (S N ) method [22]. Therein, we use a quadrature set { w m , Ω m } , con-taining weights w m and quadrature collocation points Ω m , for the angular spaceto obtain the set of equations: Ω m · ∇ ψ m + σ t ψ m = q s , ψ m = ψ (cid:16) Ω m (cid:17) , m = 1 , · · · , M, (6)with M being the total number of angles in the quadrature. Additionally, theangular integration for generic function f is defined as: (cid:90) π d Ω f ( Ω) = M (cid:88) m =1 w m f ( Ω m ) . (7)For instance, the scalar flux is expressed as φ ( (cid:126)r ) = (cid:90) π d Ω ψ ( (cid:126)r, Ω) ≈ M (cid:88) m =1 w m ψ m ( (cid:126)r ) . (8) In order to employ CFEM to solve the transport equation, Hansel, et al. [7]derived a least-squares form of transport equation by multiplying the transportequation by adjoint streaming and removal operator, i.e. L † ( Lψ − q s ) = 0 , (9)5here L † = − Ω · ∇ ( · ) + σ t . (10)The corresponding weak form of the problem in Eq. (9) with CFEM is: given afunction space V , find ψ ∈ V such that ∀ v ∈ V : (cid:73) π d Ω (cid:90) D dV vL † ( Lψ − q s ) = 0 . (11)The property of adjoint operators allows us to express Eq. (11) as (cid:73) π d Ω (cid:90) D dV Lv ( Lψ − q s ) − (cid:73) π d Ω (cid:90) ∂ D ds (cid:126)n · Ω v ( Lψ − q s ) = 0 . (12)Additionally, we require the transport equation to be satisfied on the boundaryas well: Lψ − q s = 0 . (13)Therefore, the weak form can be expressed as: (cid:73) π d Ω (cid:90) D dV Lv ( Lψ − q s ) = 0 . (14)On the other hand, one could also derive the weak form from defining aleast-squares functional of transport residual:Γ LS = 12 (cid:73) π d Ω (cid:90) D dV ( Lψ − q s ) . (15)Minimizing Eq. (15) in a discrete function space V leads to Eq. (14) as well . In the derivation above, the incident boundary conditions have been ignored.Though it is possible to impose a strong boundary condition [5], a weak bound-ary condition is chosen in this work instead. To weakly impose the boundary The procedure of the minimization is to find a stationary “point” in a specific functionspace, see Ref. [23] for details. LS = 12 (cid:73) π d Ω (cid:90) D dV ( Lψ − q s ) + 12 (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds λ (cid:0) ψ − ψ inc (cid:1) , (16)where λ is a freely chosen Lagrange multiplier. The resulting LS weak formula-tion with boundary condition is then expressed as: (cid:73) π d Ω (cid:90) D dV Lv ( L − S ) ψ + (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds λvψ = (cid:73) π d Ω (cid:90) D dV Lvq + (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds λvψ inc , (17)where the weight function v is any element of the trial space.If we choose λ = σ t (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) , then in non-void situations the LS scheme isglobally conservative in homogeneous media (i.e., where σ t is spatially indepen-dent in the whole domain). This can be seen by taking the weight function v to be 1 in Eq. (17) and performing the integration over angle to get σ t B = 0 , (18a) B := (cid:90) ∂ D ds j out − (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) ψ inc + (cid:90) D dV ( σ a φ − Q ) , (18b) j out := (cid:90) (cid:126)n · Ω > d Ω (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) ψ and Q := (cid:90) π d Ω q. (18c)Here B is the global balance: it states that the outgoing current, j out , plus theabsorption σ a φ is equal to the total source plus the incoming current. When σ t = 0, B being any value does not disturb σ t B = 0. In general, what we observeis B (cid:54) = 0. Therefore, conservation is lost.7 . A least-squares discretization allowing discontinuity on subdomaininterface Given that LS is conservative when the cross-sections are constant and non-void, we propose to reformulate the problem to solve using least-squares insubdomains where cross-sections are constant. On the boundary of these sub-domains we connect the angular fluxes using an interface condition that allowsthe fluxes to be discontinuous. The result is a discretization that can be solvedusing a procedure analogous to transport sweeps where the sweeps are oversubdomains, instead of mesh zones.In the context of a minimization problem, we can define a functional in thefollowing form:Γ
SDLS = 12 (cid:88) D i (cid:73) π d Ω (cid:90) D i dV ( L i ψ i − q si ) + 12 (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ( ψ i − ψ inc ) + 12 (cid:88) D i (cid:88) F i , j (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ( ψ i − ψ j ) , (19)where F i , j is the interface between D i and any contiguous subdomain D j . Ac-cordingly, the variational problem turns to: find ψ i in a polynomial space V such that ∀ v i ∈ V , (cid:88) D i (cid:73) π d Ω (cid:90) D i dV L i v i ( L i − S i ) ψ i + (cid:88) F i , j σ ti (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ( ψ i − ψ j ) + (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ i = (cid:88) D i (cid:73) π d Ω (cid:90) D i dV L i v i q i (20)+ (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ inc . Compared with ordinary LS method as illustrated in Eq. (17), SDLS doesnot enforce the continuity on subdomain interface. That presents possibility ofcombining the least-squares method and transport sweeps. For a given direction,
Ω, Eq. (20) can be written as a block-lower triangular system if no re-entering8nterface manifests. Therein, each block is LS applied to a subdomain.The inver-sion of this system requires solving a LS system for each subdomains, connectedvia boundary angular fluxes.
A favorable SDLS property is both subdomain conservation and global con-servation are preserved. The demonstration is similar to the derivation in Sec.2.3. Taking v i = 1, L i v i simplifies to L i v i = σ ti . (21)Accordingly, the SDLS weak form in Eq. (20) is transformed to (cid:88) D i (cid:73) π d Ω (cid:90) D i dV σ ti ( L i − S i ) ψ i + (cid:88) F i , j (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ( ψ i − ψ j ) + (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ i = (cid:88) D i (cid:73) π d Ω (cid:90) D i dV σ ti q i (22)+ (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ inc . Denote the boundary of subdomain D i by ∂ D i3 , the first integral in Eq. (22)can be expressed as: (cid:88) D i (cid:73) π d Ω (cid:90) D i dV σ ti ( L i − S i ) ψ i = − (cid:88) D i (cid:90) ∂ D i ds (cid:90) (cid:126)n i · Ω < d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ψ i + (cid:88) D i (cid:90) ∂ D i ds (cid:90) (cid:126)n i · Ω > d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ψ i + (cid:88) D i (cid:90) D i dV σ ti σ ai φ i . (23)Note that (cid:88) D i (cid:88) F i , j (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ψ i + (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ i − (cid:88) D i (cid:90) ∂ D i ds (cid:90) (cid:126)n i · Ω < d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ψ i = 0 . (24) Note that ∂ D i could either be on the problem boundary or interior interfaces. (cid:88) D i (cid:90) D i dV σ ti σ ai φ i + (cid:88) ∂ D i (cid:90) (cid:126)n i · Ω > d Ω (cid:90) F i , j ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ i = (cid:88) D i (cid:88) F i , j (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) σ ti ψ j + (cid:88) D i (cid:90) D i dV σ ti Q i (25)+ (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) (cid:126)n i · Ω < d Ω (cid:73) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ inc . Define the incoming currents from contiguous subdomain D j and problem bound-ary for subdomain D i and outgoing current as j ini , j ( (cid:126)r ) := (cid:90) (cid:126)n i · Ω < d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ i , (26a) j ini , b ( (cid:126)r ) := (cid:90) (cid:126)n i · Ω < d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ inc (26b)and j outi ( (cid:126)r ) := (cid:90) (cid:126)n i · Ω > d Ω (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) ψ i , (26c)respectively. It is then straightforward to get (cid:88) D i (cid:90) D i dV σ ti σ ai φ i + (cid:90) ∂ D i ds σ ti j outi = (cid:88) D i (cid:88) F i , j (cid:90) F i , j ds σ ti j ini , j + (cid:88) D i ∩ ∂ D(cid:54) = ∅ (cid:90) ∂ D∩D i ds σ ti j ini , b + (cid:88) D i (cid:90) D i dV σ ti Q i . (27)Assuming σ ti is spatially uniform within D i , we can further transform Eq. (27)to (cid:88) D i σ ti (cid:90) D i dV σ ai φ i + (cid:90) ∂ D i ds j outi − (cid:88) D i σ ti (cid:88) F i , j (cid:90) F i , j ds j ini , j − (cid:88) D i ∩ ∂ D(cid:54) = ∅ σ ti (cid:90) ∂ D∩D i ds j ini , b − (cid:88) D i σ ti (cid:90) D i dV Q i = 0 . (28)10r equivalently (cid:88) D i σ ti B i = 0 (29)For all nonzero σ ti , in order to make Eq. (29) true, one must have subdomain-wise conservation: B i ≡ , ∀D i ⊂ D (30)Additionally, this implies that there is the global conservation: (cid:88) D i B i = 0 (31) CFEM-SAAF is globally conservative, yet, not compatible with void andpotentially ill-conditioned in near-void situations. Therefore, efforts have beenput in alleviating CFEM-SAAF in void [17, 18, 21, 24]. We will specially give abrief review of the treatment developed in [18, 24].Therein, Laboure, et al. developed a hybrid method compatible with voidand near-void based on SAAF and LS. Denote non-void and void/uniform near-void subdomains by D n and D v , respectively . Then the hybrid formulation ispresented as: (cid:73) π d Ω (cid:90) D dV (cid:16) τ Ω · ∇ vψ + σ t vψ − (1 − σ t τ ) Ω · ∇ vψ (cid:17) + (cid:90) (cid:126)n · Ω > d Ω (cid:90) ∂ D ds (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) vψ = (cid:73) π d Ω (cid:90) D dV (cid:16) τ Ω · ∇ v + v (cid:17) q s + (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) vψ inc , (32)where τ = /σ t , (cid:126)r ∈ D n /c (cid:126)r ∈ D v . (33) In our work, D v is defined with σ t < .
01 cm − . is a freely chosen constant and usually set to be 1. One notices that in non-voidsubdomain, the formulation resembles CFEM-SAAF. At the same time, globalconservation is preserved. Therefore, the method is named “SAAF-ConservativeLS (SAAF-CLS)”.Consider a homogeneous near-void problem. The weak form in Eq. (32) canbe transformed to (cid:73) π d Ω (cid:90) D v dV Lv ( Lψ − q s ) + (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D v ds σ t (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) v (cid:0) ψ − ψ inc (cid:1) + (cid:73) π d Ω (cid:90) D v dV cv ( Lψ − q s ) + (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D v ds c (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) v (cid:0) ψ − ψ inc (cid:1) = 0 . (34)In void/near-void, CLS is conservative:( σ t + c ) B = 0 . (35)We therefore modify SDLS such that void treatment based on Eq. (34) isincorporated: (cid:88) D i ⊂D n (cid:73) π d Ω (cid:90) D i dV L i v i ( L i − S i ) ψ i + (cid:88) F i , j σ ti (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ( ψ i − ψ j ) + (cid:88) D i ∩ ∂ D(cid:54) = ∅D i ⊂D n (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ i + (cid:88) D i ⊂D v (cid:73) π d Ω (cid:90) D i dV ( cv i + L i v i ) ( L i − S i ) ψ i + (cid:88) F i , j (cid:90) (cid:126)n i · Ω < d Ω (cid:90) F i , j ds ( c + σ ti ) (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ( ψ i − ψ j ) + (cid:88) D i ∩ ∂ D(cid:54) = ∅D i ⊂D v (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds ( c + σ ti ) (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ i (36)= (cid:88) D i ⊂D n (cid:73) π d Ω (cid:90) D i dV L i v i q i + (cid:88) D i ∩ ∂ D(cid:54) = ∅D i ⊂D n (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds σ ti (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ inc . + (cid:88) D i ⊂D v (cid:73) π d Ω (cid:90) D i dV ( cv i + L i v i ) q i + (cid:88) D i ∩ ∂ D(cid:54) = ∅D i ⊂D v (cid:90) (cid:126)n i · Ω < d Ω (cid:90) ∂ D∩D i ds ( c + σ ti ) (cid:12)(cid:12)(cid:12) (cid:126)n i · Ω (cid:12)(cid:12)(cid:12) v i ψ inc . . Numerical results The implementation is carried out by the
C++ finite element library deal.II [25]. The Bi-conjugate gradient stabilized method [26] is used as linear solverfor SDLS and CFEM-SAAF-CLS while LS and CFEM-SAAF are solved usingconjugate gradient method. Symmetric successive overrelaxation [27] is usedas preconditioner for all calculations with relaxation factor fixed at 1.4. Inall tests, we also include results from solving the globally conservative self-adjoint angular flux (SAAF) equation with CFEM as a comparison[2, 21]. Inall problems, piecewise linear polynomial basis functions are used. In the 1Dexamples, the angular quadrature is the Gauss quadrature. In the 2D test, thequadrature is a Gauss-Chebyshev quadrature with azimuthal point number oneach polar level specified in a way similar to level-symmetric quadrature. SeeRef. [28] for further details.
The first problem is Reed’s problem in 1D slab geometry [29] designed to testspatial differencing accuracy and stability. The material properties are listed inTable 1. Specifically, a void region is set at x ∈ (3 , in angleare presented in Figure 1a with zoomed results for x ∈ (4 . , .
25) cm in Figure1b. In this figure 32 cells are used for the CFEM-LS (green line), CFEM-SAAF-CLS (blue line) and SDLS (red lines) calculations. For SDLS, the interfacesare set at x = 3 , , x ∈ (3 ,
5) is flat and therelative error is lower than 3 × − . By defining the relative balance as: B rel = | B | (cid:90) (cid:126)n · Ω < d Ω (cid:90) ∂ D ds (cid:12)(cid:12)(cid:12) (cid:126)n · Ω (cid:12)(cid:12)(cid:12) ψ inc + (cid:90) D dV Q , we found that the SDLS is globally conservative as B rel = 5 . × − . On theother hand, LS solution is distorted in void and the balance is poor ( B rel = 3 . × − ). Meanwhile, we notice that CFEM-SAAF-CLS gives a better solution inthe graph norm than LS in most regions of the problem. Yet, solution in voidregion is heavily affected by the thick absorber ( x ∈ (5 ,
6) cm) as continuity isenforced at x = 5 cm. Table 1: Material configuration for Reed’s problem. x [cm] (0,2) (2,3) (3,5) (5,6) (6,8) σ t [cm − ] 1 1 0 5 50 σ s [cm − ] 0 . . Q φ Ref.LSCFEM-SAAF-CLSSDLS (a) Reed’s problem result comparison. x [cm] (zoomed) φ Ref.LSCFEM-SAAF-CLSSDLS (b) Zoomed results x ∈ (4 . , .
25) cm.
Figure 1: Reed’s problem results from different methods. 32 cells are used.
Another test problem is a 1D slab pure absorber problem firstly proposed byZheng et al. [38]. There is a unit isotropic incident angular flux on left boundaryof the slab. No source appears in the domain. In this problem σ t = . − x < − otherwise . [cm], 16 cells φ S DLSLS (a) 16 cells. x [cm], 40 cells φ S DLSLS (b) 40 cells.
Figure 2: Two-region absorption results comparison with LS, SAAF and SDLS.
The results in Figure 2 compare coarse solutions using different methods.The reference solution was computed with first order S using diamond differ-ence using 2 × spatial cells. The LS and SAAF solutions in the thin region( x < x > DoF counts -4 -3 -2 -1 R e l a t i v e E rr o r s f o r Lea k age LSCFEM-SAAFSDLS
Figure 3: Leakage errors vs. DoF counts per direction at the right boundary.Table 2: Relative global balances with different methods (absolute values).
Number of Cells 20 40 80 160 320LS 8 . × − . × − . × − . × − . × − CFEM-SAAF 5 . × − . × − . × − . × − . × − SDLS 2 . × − . × − . × − . × − . × − k -eigenvalue problem A one-group k -eigenvalue problem is also tested in 1D slab geometry. Anabsorber region in x ∈ (0 , .
3) cm is set adjacent to a multiplying region in x ∈ (0 . , .
5) cm. Material properties are presented in Table 3. Reflective BCsare imposed on both sides of the slab. The configuration is to mimic the impactfrom setting a control rod near fuel.A reference is provided by CFEM-SAAF using 20480 spatial cells. With thepresence of the strong absorber, the scalar flux has a large gradient near x = 0 . x = 0 . Table 3: Material configuration for k -eigenvalue problem. x [cm] (0,0.5) (0.5,1.5) σ t [cm − ] 5 1 σ s [cm − ] 0 0 . νσ f . x [cm] φ SDLSCFEM-SAAFLS
Figure 4: Scalar flux from different methods.
We also examine the k eff in Figure 5. We observe that SDLS converges tothe reference k eff in the graph norm using only 10 spatial cells in Figure 5a. Yet,CFEM-SAAF slowly converges to the reference and LS still has a large errorwith 160 cells. Figure 5b shows k eff error change with respect to DoF counts and17 DoF counts per direction0.70.80.911.1 k e ff Ref.SDLSCFEM-SAAFLS (a) k eff values from different methods. DoF counts per direction10 -2 -1 k e ff ab s . e rr o r s [ p c m ] SDLSCFEM-SAAFLSh Ref. (b) Absolute errors in pcm.
Figure 5: k-eigenvalue results. illustrates the benefit of adding the subdomain interface discontinuity. With 5cells (the points with smallest DoF counts), SDLS has over an order lower of k eff error than CFEM-SAAF. When refining the spatial mesh, SDLS shows roughlysecond-order convergence and has an error three orders of magnitude lower thanCFEM-SAAF with 160 cells (the points with largest DoF counts) although bothhave the global conservations. The last test is a modified one-group iron-water shielding problem [30] usedto test accuracy of numerical schemes in relatively thick materials (see the con-figuration in Figure 6a). Material properties listed in Table 4 are from thethermal group data. S is used in the angular discretization. The reference isSDLS using 1200 × ×
120 cells), we see SDLS graphically agrees with CFEM-SAAFand LS in the line-out illustrated in Figure 7a. We further examine the absorp-tion rate errors in iron. As shown in Figure 7b, LS and CFEM-SAAF presentssimilar spatial convergence rates. However, LS presents lower accuracy thanCFEM-SAAF with the presence of material interface between iron and water.18n the other hand, by setting two interfaces between iron and water and intro-ducing extra DoFs on the subdomain interfaces, SDLS converges to the referencesolution much earlier than the other methods.
Table 4: Material cross sections in one-group iron-water test.
Materials σ t [cm − ] σ s [cm − ]Water 3 . . . . (a) Problem configuration. (b) Scalar flux distribution from reference. Figure 6: Results from the iron water problem.
In addition, we present a timing comparison of different methods in Fig-ure 8 using the same error data employed in Figure 7b. Both CFEM-SAAFand SDLS are more accurate given a specific computing time. For a specifichigh error level ( > × − ) obtained using coarse meshes, the SDLS is moretime-consuming than CFEM-SAAF. However, for low error levels ( < × − )achieved by refining the mesh, SDLS cost much less total computing time thanCFEM-SAAF. On the other hand, when the total computing time increases,SDLS error drops faster than both CFEM-SAAF and LS. Overall, the linearsolver presents reasonable efficiency for SDLS despite the fact that the system19 x [cm], y=10.1cm -6 -5 -4 -3 -2 -1 φ SAAF: 120x120SDLS: 120x120LS: 120x120CFEM-SAAF Ref. (a) Line-out plot for x = 10 . DoF counts per direction [s] -6 -5 -4 -3 -2 -1 R e l a t i v e e rr o r s CFEM-SAAFSDLSLSh Ref.h Ref. (b) Relative iron absorption errors.
Figure 7: Results from the iron water problem. is non-SPD. -2 -1 Total computing time [s] -6 -5 -4 -3 -2 -1 R e l a t i v e e rr o r s CFEM-SAAFSDLSLS
Figure 8: Relative errors of absorption rate in iron vs total computing time. . Concluding remarks and further discussions In this work, we proposed a subdomain discontinuous least-squares dis-cretization method for solving neutral particle transport. It solves a least-squares problem in each region of the problem assuming spatially uniform crosssections within each region and couples the contiguous regions with discon-tinuous interface conditions. We demonstrated that our formulation preservesconservation in each subdomain and gives smaller numerical error than eitherSAAF or standard least squares methods.Though SDLS is angularly discretized with S N method in this work, P N would be applicable in the angular discretization as well. Therein, LSP N [31,32, 33] is applied in each subdomain with Mark-type boundary condition usedas interface condition [34]. Further, since SDLS allows the discontinuity on thesubdomain interface, different angular schemes, e.g. S N and P N can be used indifferent subdomains. In fact, the angular coupling is desirable for SAAF in thecode suite Rattlesnake [35] at Idaho National Laboratory. An initial imple-mentation of the coupling has been developed in [36] through enforcing strongcontinuity of angular flux on coupling subdomain interfaces. And later, theinterface discontinuity is applied to SAAF to allow an effectively improved cou-pling scheme reported in [37, 38] and implemented in
Rattlesnake to improvethe neutronics calculations.
We restrict the method to the configuration that total cross section is uni-form within subdomains. In cases total cross sections varies smoothly withinsubdomains, one could simply weight the subdomain functional with reciprocaltotal cross section while weighting the subdomain interface functional with 1.The rationale is such that weak form resembles the conservative CFEM-SAAFin the corresponding subdomain with spatially varying cross sections (see Ref.[6] for demonstration). 21he other interesting direction is to explore the effects of re-entrant curvedsubdomain interfaces, which is not considered in current study. The efficacy ofthe method and according solving techniques in these cases are not known andworthy of investigation.
Acknowledgements
W. Zheng is thankful to Dr. Bruno Turksin from Oak Ridge National Labo-ratory for the suggestions on multi-D implementation. Also, Zheng appreciatesDr. Yaqi Wang from Idaho National Laboratory for his suggestions to make thenotation simple and clear. He wants to give appreciation to Hans Hammer fromTexas A&M Nuclear Engineering for running reference tests as well. Finally, wewish to thank the anonymous reviewers with their invaluable advice to improvethe paper to current shape.This project was funded by Department of Energy NEUP research grantfrom Battelle Energy Alliance, LLC- Idaho National Laboratory, Contract No:C12-00281.