- Email: [email protected]
A two-mesh adaptive mesh refinement technique for SN neutral-particle transport using a higher-order DGFEM
Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A two-mesh adaptive mesh refinement technique for SN neutral-particle transport using a higher-order DGFEM Jean C. Ragusa ∗ , Yaqi Wang Department of Nuclear Engineering, Texas A&M University, College Station, TX, 77843-3133, USA
article
abstract
info
Article history: Received 3 April 2009 Received in revised form 7 December 2009
We present a two-mesh Adaptive Mesh Refinement technique for the SN neutral-particle transport equation in 2D. A high-order (up to order 4) Discontinuous Galerkin (DG) finite element discretization is employed with standard upwinding. A suitable Diffusion Synthetic Acceleration (DSA) scheme is derived to precondition the transport solves. The DSA scheme, also based on a DG discretization, is derived directly from the discretized highorder DG SN transport equation. Results are presented to demonstrate the superiority of AMR for DSA-accelerated SN transport solves. © 2009 Elsevier B.V. All rights reserved.
Keywords: Adaptive mesh refinement SN transport Discontinuous Galerkin finite element method Diffusion synthetic acceleration
1. Introduction The steady-state transport of mono-energetic neutral particle (e.g., neutrons, photons) through matter is described by the following integro-differential equation
E + σ (Er ) Ψ = E ·∇ Ω
Z 4π
E0 ·Ω E )Ψ (Er , Ω E 0 )dΩ 0 + Q (Er , Ω E ), σs (Er , Ω
in D × S 2 ,
(1)
with the general boundary condition
E ) = Ψ inc (Erb , Ω E)+ Ψ (Erb , Ω
Z E 0 ·Enb >0 Ω
E0 → Ω E )Ψ (Erb , Ω E 0 )dΩ 0 on Erb ∈ ∂ D − , β(Erb , Ω
(2)
E the direction of flight, Er the position, σ (Er ) the total cross section, where DR is the spatial domain, S 2 the unit sphere, Ω E0·Ω E ) the scattering cross section, σs (Er , Ω E0·Ω E ) the differential scattering cross section, Q (Er , Ω E ) the σs (Er ) = 4π dΩ σs (Er , Ω E ) the incoming flux, and β(Erb , Ω E0 → Ω E ) the boundary redistribution kernel (albedo). The incoming source term, Ψ inc (E rb , Ω domain boundary ∂ D − is defined by E · nE(Erb ) < 0 , ∂ D − = Erb ∈ ∂ D , Ω
(3)
E(Erb ) the outward unit normal vector at position Erb . For brevity, we only consider Dirichlet conditions hereafter (i.e., with n β = 0). Despite Eq. (1) being linear, its numerical solution is challenging, even on current computers, due to 1. the high dimensionality of the phase-space. In complete generality, the angular flux Ψ depends upon seven variables: three for the standard physical space, two for the angular variable, one variable for the particles’ energy, and one for time. After the application of an implicit time integrator and the multigroup approximation in energy, an equation of the form given in Eq. (1) is obtained;
∗
Corresponding author. Tel.: +1 979 862 2033. E-mail addresses: [email protected] (J.C. Ragusa), [email protected] (Y. Wang).
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.12.020
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3179
2. the presence of singularities in the transport solution [1] and the greatly disparate values of σ (E r ) and σs (E r ) in practical applications; 3. the change in the nature of the PDE, from a hyperbolic behavior in materials with small scattering ratio, i.e., c = σσs 1, to an elliptic behavior for highly diffusive materials, i.e., c ' 1, [2]. In radiation transport production codes, low-order spatial differencing techniques are often employed, either on structured [3] or unstructured [4] grids. In addition, while Adaptive Mesh Refinement (AMR) techniques are now widely used in many science and engineering disciplines [5–7], only a limited number of references are available regarding the applications of mesh adaptivity to radiation transport. The higher-order spatial adaptivity technique described here is applied to the SN angular approximation of Eq. (1). Patch-AMR techniques [8–10] for the SN transport equation have been devised and are based on a hierarchy of nested grids. Some of the drawbacks of patch-AMR include the fact that the physics are not followed as closely as possible (the extent of refined patch being often too large), leading to more unknowns than needed, and the need to converge inflow/outflow values in between nested grids (a feature that is not present in cellbased AMR). Oftentimes, a gradient-based error estimator is employed to drive the adaptive mesh refinement in Cartesian structured geometries but it is overly conservative for high-order spatial discretization schemes. In [11], a local refinement (cell-based AMR) technique is described for SN transport spatially discretized using low-order finite differences, where the value of the particle mean free path (MFP) in a given cell is used as a mesh refinement criterion. While this approach takes into account the size of potential internal layers at a given location in the domain, it does not account for the actual smoothness of the solution at these locations and is, therefore, far from optimal; for example, within optically thick areas, the solution may well be approximated by a smooth spatial representation on coarse meshes despite the smallness of the mean free path. In [12], a lower-order angular flux solution serves to compute a converged scattering source which is then used as a fixed source term for a higher-order transport solve. The differences between the two solutions is then used to prescribe local mesh refinements; their technique was demonstrated for one-group, 2D rectangular structured meshes. The error estimate employed in our present work is based on the difference between the numerical solutions computed on coarser and finer AMR meshes at any adaptivity cycle; the local error hereby obtained is used to prescribe the new coarser AMR mesh at the next adaptivity cycle, the finer AMR mesh being obtained by refinement of the coarser mesh. Additionally, preconditioning techniques are necessary to accelerate the transport solution in highly diffusive configurations; it is well established that the spatial discretization of the Diffusion Synthetic Acceleration (DSA) equations must be ‘‘consistent’’ in some sense with the one used in the SN transport solves in order to yield unconditionally stable and efficient DSA schemes [13,14,4,15,16]. With the exception of [17] for block structured AMR with linear finite elements, little work has been reported to devise DSA schemes (i) for higher-order spatial discretizations and, more importantly, (ii) for unstructured AMR meshes. Using a variational argument, we obtain here a DGFEM DSA procedure (since it is based on a DG method, it is suitable for AMR meshes) and we show its relation with the standard Interior Penalty (IP) DGFEM method for elliptic PDEs [18,19]. The work presented here describes (i) an AMR technique to accurately solve the SN transport equation using a high-order Discontinuous Galerkin Finite Element Method (DGFEM) in space and (ii) an efficient DSA preconditioner for AMR meshes. Section 2 briefly reviews the SN approximation and describes the DGFEM variational form of the SN transport equation. In Section 3, we present our modified IP DSA preconditioner. The AMR technique and error estimator employed are explained in Section 4. Numerical results are given in Section 5 and we provide conclusions in Section 6. 2. Numerical solution techniques 2.1. Space/angle discretization In this section, we briefly review the SN angular discretization for Eq. (1) and present a DGFEM high-order spatial E m , wm }, m = 1, . . . , M, consisting of M directions Ω E m and weights wm , discretization. Given an angular quadrature set {Ω the SN discrete ordinates [2] version of Eq. (1) is Na X n X E + σ Ψm (Er ) = Em · ∇ E m ), Ω σs,n (Er )Φn,k (Er ) + Qn,k (Er ) An,k (Ω
(4)
n=0 k=−n
for m = 1, . . . , M, supplemented with appropriate boundary conditions. The scattering source term in Eq. (4) has been manipulated by decomposing the differential scattering cross section on Legendre polynomials (up to anisotropy order Na ) and introducing the flux moments
Φn,k (Er ) =
Z 4π
E )Ψ (Er , Ω E) ≈ dΩ An,k (Ω
M X
E m )Ψm (Er ), wm An,k (Ω
(5)
m=1
where An,k represents the real spherical harmonic function of degree n and order k with the Schmidt semi-normalization.
3180
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
The spatial domain D is triangulated into an unstructured mesh Th = {K , ∪K ∈Th K = D } and we employ a DG discretization technique in space, with Pp (K ) denoting the space of polynomial of degree ≤ p on element K . Hierarchical E m , the DG basis functions (up to order 4) are chosen for the spatial discretization [20,21]. For a given streaming direction Ω scheme with standard upwinding reads: Find Ψm ∈ Vp (K ), such that ∀Ψm∗ ∈ Vp (K ),
E + σ )Ψm∗ Em · ∇ Ψm , (−Ω
K
Na X n X
E m ) σs,n Φn,k + Qn,k , Ψm∗ , + Ψm− , Ψm∗ ∂ K + − Ψm− , Ψm∗ ∂ K − = An,k (Ω K
(6)
n=0 k=−n
where ∂ K − is the inflow boundary, ∂ K + is the outflow boundary of element K . The traces Ψ + and Ψ − are defined with E m , i.e., on an inflow boundary, Ψ + (resp. Ψ − ) is the value of Ψ taken from within element respect to the particle direction Ω K (resp. from the upwind neighbor) and on an outflow boundary, Ψ + (resp. Ψ − ) is the value of Ψ taken from the downwind neighbor (resp. from within element K ). These definitions are expressed as follows:
E m) Ψm± ≡ lim Ψm (Er + sΩ
(7)
s→0±
(f , g )K ≡
Z
hf , g ie ≡
Z
fg
(8)
E m · nE(Er )|f g |Ω
(9)
K
e
Em R 0 ∂ K ± ≡ Er ∈ ∂ K , nE(Er ) · Ω
(10)
E is the outward unit normal vector of element K and e represents any edge (or face in 3D) of K . Since the finite where n E is not defined at the element interfaces. Here, the upwind element approximation used is discontinuous, the quantity ΩΨ E) defined as scheme is obtaining by explicitly employing a numerical flux H (f + , f − , n E · nEf − E) = Ω H (f + , f − , n
(11)
such that
Z
E) Ψm∗ = Ψm− , Ψm∗ ∂ K + − Ψm− , Ψm∗ ∂ K − . H (Ψm+ , Ψm− , n
∂K
(12)
Eq. (12) was used in deriving Eq. (6). Transport solves are based on Eq. (6). By integrating Eq. (6) by parts, an equivalent formulation is obtained:
E + σ )Ψm , Ψm∗ Em · ∇ (Ω
K
Na X n X
E m ) σs,n Φn,k + Qn,k , Ψm∗ , + [[Ψm ]], Ψm∗ ∂ K − = An,k (Ω K
(13)
n=0 k=−n
where the inter-element jump is given by [[Ψm ]] = Ψm+ − Ψm− . Eq. (13) will be useful in deriving the diffusion synthetic accelerator in Section 3. 2.2. Iterative procedure The system of equations formed by Eq. (6) needs not to be assembled for all elements but can be solved element by E m . The order in which the elements are element. This strategy is often referred to as a ‘‘transport sweep’’ in direction Ω solved constitutes the graph of the sweep (for brevity, we do not expand on situations where the graph can present some dependencies, in which case these dependencies are simply lagged within the iterative procedure; for additional details, refer to [4]). With appropriate numbering, the angular fluxes and flux moment unknowns are described by 9 and 8, respectively, and the transport sweeps over all directions can be represented in the following matrix form: L9 = M Σ8 + Q 8 = DΨ ,
(14)
where L is the block-diagonal matrix due to particle streaming and total interaction (one block = one direction), 6 is the scattering matrix which operates on the flux moments, M is the moment-to-direction matrix, D is the direction-to-moment matrix (hereafter, we assume MD = I for conciseness in the derivations of Section 3, although not all quadrature formulae necessarily satisfy this). We can recast Eq. (14) in terms of flux moments as follows: I − DL−1 M Σ 8 = DL−1 Q.
(15)
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3181
None of these global matrices are explicitly formed and L is inverted in a ‘‘matrix-free’’ fashion using transport sweeps. Eq. (15) is solved using GMRes. Alternatively, a Richardson iteration is employed, leading to the scattering Source Iteration (SI) technique
9(i+1) = L−1 M68(i) + Q 8(i+1) = D9(i+1) ,
(16)
where the previous flux moment iterate 8(i) is utilized to build the scattering source term. Both the SI scheme and the GMRes solution are known to be slow converging for highly diffusive problems and a Diffusion Synthetic Acceleration (DSA) [16] must be implemented to address this issue. This DSA technique can be utilized as a preconditioner for both the SI scheme or the GMRes solves and is discussed next. 3. Diffusion preconditioner for AMR meshes For highly diffusive configurations (i.e., materials with scattering ratios c close to 1 and optically thick domains), the standard iterative techniques to solve the transport equation can become quite ineffective due to their slow convergence properties and a DSA preconditioner is needed to accelerate the convergence of the process. The diffusion operator must be spatially discretized in a manner that is ‘‘consistent’’ with the discretization of the transport equation [13–16] to ensure unconditional stability (however, a diffusion-based synthetic accelerator is not necessarily effective for multidimensional problems containing highly diffusive regions with discontinuities in material properties, [22]). To obtain consistently discretized DSA equations, the zeroth and first angular moments of the discretized transport equation are taken, leading to a mixed DGFEM discretization of a diffusion equation. The mixed formulation results in a nonsymmetric PD saddle point problem [4] that is CPU expensive and partially consistent DSA schemes, motivated by the difficulties associated with the algebraic elimination of the vector unknowns to yield an elliptic diffusion equations, are more popular. With partial consistency, it is hoped that the reduction in complexity outweighs the degradation in effectiveness but few partially consistent schemes have been applied to AMR meshes [17]. Additionally, finite-element-based partially consistent schemes employ a Continuous Finite Element (CFE) discretization of the diffusion operator [23]. However, the algorithmic and implementation simplicity of employing an AMR-DG solution for the transport solves would be severely reduced if the DSA preconditioner solves were based on a CFE method requiring hanging nodes treatment. In our new approach, we have chosen to derive partially consistent DSA schemes employing a DGFEM discretization by directly deriving them from the DGFEM discretization of the SN transport equations. Our scheme belongs to the family of partially consistent DSA methods because in our derivation, only the zeroth angular moment of the discretized transport equation is used and we assume that the continuous first angular moment of the transport equation (Fick’s law) to be verified in order to eliminate the current unknowns EJ = [Φ1,−1 ; Φ1,0 ; Φ1,1 ]T . The resulting DGFEM discretization for the diffusion equation is remarkably similar to the Interior Penalty (IP) stabilization method for diffusion equations solved using a DGFEM. Due to the discontinuous nature of the DGFEM approximation, it is particularly well suited for meshes arising in AMR calculations, i.e., hanging nodes are seamlessly incorporated into the DG method. This property also leads to an easy implementation of high-order test functions. For brevity, we only outline the DSA derivation here. Step 1: Using Eq. (13), we have the following variational form the DGFEM SN transport formulation: h Find Ψm ∈ WD , ∀Ψm∗ ∈ WDh , m = 1, . . . , M such that: b(Ψ , Ψ ∗ ) = l(Ψ ∗ ) +
N X n X
σs,n Φn,k , Φn∗,k
D
,
(17)
n=0 k=−n
where the bilinear and linear forms are b(Ψ , Ψ ∗ ) =
M X
E + σ )Ψm , Ψm∗ Em · ∇ wm ( Ω
D
m=1
l(Ψ ∗ ) =
N X n X n=0 k=−n
+
M X
wm [[Ψm ]], Ψm∗+ E i
m=1
Qn,k , Φn∗,k
D
+
X
X
(18)
h
wm Ψminc , Ψm∗ e .
(19)
e∈∂ D Ω E m ·Enb <0
h WD = Ψ ∈ L2 (D ); Ψ |K ∈ V (K ), ∀K ∈ Th and Ehi = ∪K ∈Th ∂ K \ ∂ D is the set of all interior edges. Note that Eq. (17) is the zeroth angular moment of the discretized transport equation, hence the weighted summation over all directions in Eqs. (18) and (19).
Step 2: Assume a linear dependence in angle (i.e., P1 approximation) for the angular flux Ψm and the test function Ψm∗
Ψm =
1 4π
1 E m ) and Ψm∗ = E m) (φ + 3EJ · Ω (φ ∗ + 3EJ ∗ · Ω 4π
(20)
3182
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
E and with φ ≡ Φ0,0 the scalar flux and EJ the net current. Further assume that the continuous Fick’s law (EJ = −D∇φ E ∗ ) to be valid, yielding: EJ ∗ = D∇φ Ψm =
1 4π
E ·Ω E m) (φ − 3D∇φ
and
Ψm∗ =
1 4π
E ∗·Ω E m ). (φ ∗ + 3D∇φ
(21)
Step 3: Insert Eq. (21) in Eq. (17) and carry out the angular summation. After some considerable algebra, the application of angular quadrature properties, and further simplifications [24], the following Modified Interior Penalty (MIP) DGFEM formulation is obtained:
E ∇φ E ∗ + κeMIP [[φ]], [[φ ∗ ]] i + [[φ]], { D∂n φ ∗ } i + D∇φ, Eh Eh D 1 1 + { D∂n φ}}, [[φ ∗ ]] E i + κeMIP φ, φ ∗ ∂ D − φ, D∂n φ ∗ ∂ D − D∂n φ, φ ∗ ∂ D h 2 2 E ∗ + J inc , φ ∗ E1 , 3D∇φ , lMIP (φ ∗ ) = Q0 , φ ∗ D + Q ∂D bMIP (φ, φ ∗ ) =
σa φ, φ ∗
D
(22) (23)
D
with the incoming current given by: J inc =
X
E m · nE(Erb )|Ψminc . wm |Ω
(24)
E m ·En(Erb )<0 Ω
[[]] and { } represent the inter-element jumps and averages on interior edges, respectively. The absorption cross section 1 σa = σ − σs and the diffusion coefficient D = 3(σ −σ are not imposed or coerced but are naturally obtained from this s,1 ) derivation. The penalty coefficient κeMIP in Eq. (22) is given by:
1 κeMIP = max κeIP ,
(25)
4
where the value of 1/4 is the value obtained by strictly applying the procedure, κeIP =
c (p+ ) D+ + 2 h⊥
+
c (p− ) D− − 2 h⊥
is the standard IP
penalty value (see, e.g., [25]); the ± symbols in κ represent the values in the adjacent cells K and K sharing a given edge ± e, D is the diffusion coefficient, h± ⊥ is the length of cell K orthogonal to the edge, c (p) = p(p + 1), and p is the polynomial order. We note that Eq. (22) is SPD and can, therefore, be solved using a standard preconditioned Conjugate Gradient (CG) solver; we used SSOR as the preconditioner for CG. A Fourier analysis of the iteration operator (DSA+transport solves) shows that our modified scheme (MIP) yields unconditional stability when employed as a diffusion synthetic accelerator, in both the thin and thick diffusion limits, whereas the standard IP scheme (κ = κ IP ) or the original scheme (κ = 1/4) are only conditionally stable accelerators. This Fourier analysis is discussed after the derivation of the MIP-DG scheme. Steps 1–3 describe how a discretized diffusion operator, that is partially consistent with the discretized transport equation, is obtained. Next, we briefly present how this MIP-DGFEM diffusion is employed to accelerate transport sweeps. For a given as a preconditioner SI iteration, one transport solve yields 8(i+1/2) = DL−1 M ΣΦ(i) + Q (this is simply Eq. (16) in which the new iterate has −
IP e
+
now been labeled i + 1/2); next we define the error e(i+1/2) = 8 − 8(i+1/2) where 8 represents the converged solution. This error satisfies the following transport equation I − DL−1 M6 e(i+1/2) = DL−1 M6 8(i+1/2) − 8(i) ,
(26)
= 6 8(i+1/2) − 8(i) . In the error equation, Eq. (26), the external volumetric source is
or, equivalently, (DLM − 6) e(i+1/2) no longer present, the Dirichlet boundary conditions are now homogeneous (vacuum) and the new volumetric source is the scattering source due to the difference in successive iterates. Nonetheless, this error equation is just as difficult to solve as Eq. (16) and a simpler estimate for the error term e˜ (i+1/2) is sought by replacing the (DLM − 6) operator with the diffusion operator obtained through steps 1–3, resulting in the following equation: Ae˜ (i+1/2) = 6 8(i+1/2) − 8(i) ,
(27)
where A is the matrix notation for the MIP-DGFEM bilinear form of Eq. (22). Finally, the DSA-accelerated SI scheme (or preconditioned Richardson iteration) reads as follows:
8(i+1) = 8(i+1/2) + e˜ (i+1/2) = DL−1 M ΣΦ(i) + Q + A−1 6 8(i+1/2) − 8(i) ,
(28)
and the preconditioned GMRes solve is I + A−1 6
I − DL−1 M6 8 = I + A−1 6 DL−1 Q.
(29)
Fig. 1 presents the results of a Fourier Analysis for a 2D homogeneous medium with periodic boundary conditions; the mesh consists of two triangles obtained by division of a single square; the scattering ratio is chosen to be c = 0.9999. It
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3183
Fig. 1. Fourier analysis for the various DSA forms as a function of the mesh optical thickness, homogeneous infinite medium case (the right pane is a zoomed portion of the left pane).
is well known that spectral radius of SI process is also equal to c, hence SI can be arbitrarily slow to converge. In Fig. 1, we have graphed the spectral radius obtained with the DSA-accelerated SI schemes when using the following three penalty coefficients: κ = 1/4, κ = κ IP (the standard IP formulation), and κ = κ MIP (the Modified IP formulation). The results are presented as a function of the mesh optical thickness (expressed in MFP). We can note that the schemes employing either κ = 1/4 or κ = κ IP can have spectral radii greater than 1, and hence are only conditionally stable for some range of optical thickness, whereas the MIP scheme has maximum spectral radius of about 0.48 and thus always converges. An S4 angular quadrature was used for this Fourier analysis. 4. A two-mesh AMR approach for high-order SN transport In radiation transport and radiative transfer, lower-order spatial schemes are most commonly used. The AMR techniques employed with such schemes typically rely on gradient-based error estimates [8], with a few exceptions, e.g., [26], where an adjoint-based estimator is derived, but at the cost of computing an adjoint solution. Here, we wish to obtain error estimates for high-order spatial approximations and, therefore, cannot rely on simpler estimates such as the gradient of the numerical solution. We have chosen to follow the approach proposed in [7,27], where two adapted AMR meshes are utilized at each adaptivity cycle: a coarse AMR mesh Tkh and a fine AMR mesh Tkh/2 obtained by uniform refinement of Tkh ; in the preceding, k denotes the mesh adaptivity cycle index. 4.1. Error estimation procedure The error estimate at cycle k is simply the difference in the solutions obtained on Tkh/2 and Tkh . This is computed in the L2 norm, element by element. However, in steady-state transport algorithms, the angular fluxes Ψm are never kept (except on albedo boundaries and for elements causing dependencies in the graph of the sweep), hence our error estimate is based on angle-integrated quantities such as the scalar flux φ or the current EJ = [J x ; J y ; J z ]T , as shown in Eq. (30) for 2D:
i R h x y y x 2 2 ( J − J ) + ( J − J ) h h h/2 h/2 K R x µkK = ∀K ∈ Tkh . 2 + (J y )2 ( J ) h h D
(30)
We tested error estimates based on either the scalar flux or the current and did not find significant changes in the AMR results and grids. The motivation to employ an estimate based on the current comes from prior work on AMR for the diffusion equation [28] where the diffusion approximation for the current is a mathematically appropriate choice. For efficiency, we first solve for the transport solution on the coarser mesh and use the converged scattering source as an initial guess for the solution on the finer mesh, hence reducing the number of transport iterations needed on the finer mesh. The projection operation in Eq. (30) is defined as: h , such that, Find Πh Φ ∈ WD
Πh Φ , Φ ∗
D
= Φ, Φ∗
D
,
∀Φ ∗ ∈ WDh
(31)
This projection is easily carried out element by element, due to the discontinuous nature of the spatial discretization. The criterion for refinement is defined as follows: an element K of Tkh is selected for refinement if
µkK ≥ α max µkK 0 ,
(32)
K 0 ∈Tkh
where α is a user-defined fraction (we used α = 0.3). This criterion allows us to focus the computational effort on elements with the largest errors and tends to equi-distribute the spatial error. The global error, obtained by summing the error estimates over all elements, i.e.,
ηk =
sX K ∈Tkh
µkK ,
(33)
3184
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
(a) Same refinement level.
(b) Coarser upwind element.
(c) Finer upwind element.
Fig. 2. Edge-coupling cases for refinement level differences ≤ 1.
is used to terminate the adaptivity procedure (ηk < tolAMR ). The final result of the calculation is the last solution obtained on the finer AMR mesh. 4.2. Mesh coupling Once an element has been flagged for refinement, it becomes inactive and the child-elements are the new active cells. Appropriate refinement rules allow for an efficient visit of the tree structure. Any level of refinement difference is allowed between neighboring elements. The element-coupling algorithm is based on recursive calls to the function dealing with the 1-level refinement difference. The cases of 0- and 1-level refinement difference are shown in Fig. 2 and are dealt with as follows. Fig. 2(a) is the standard situation without refinement between neighbors, where the edge-coupling from element K 0 to K is simply achieved by LAB
Ci,j ΨK 0 , (34) 6 where Ci,j is a coupling matrix. In 2D problems, Ci,j is akin to a 1D mass matrix (see Eq. (6)) but contains several empty rows and columns because it applies to the whole solution vector of element K 0 (only the components of ΨK 0 belonging to edge K ∩ K 0 play a role in edge-coupling). It was our deliberate choice to achieve edge-coupling through whole-element operations, as this is necessary for the diffusion solves (where the edge integrals present gradient terms and hence require the entire element solution). A few templates for the Ci,j matrices need to be stored, depending on the local numbering of vertices i and j; thanks to the hierarchical nature of the basis functions used, the templates for any polynomial order are embedded in the single template for the highest order employed. Fig. 2(b) shows the situation when the contribution to element K results from the contribution of two smaller elements (1-level refinement difference). The contribution is simply computed for each element separately, yielding LAB
LBC
(35) PiT Ci,j ΨK 0 , 1 6 6 2 where Pk are cell prolongation matrices operating on the local solution vector and providing the solution vectors on the four child-elements. These matrices are independent of the element shape. Finally, the case shown in Fig. 2(c) requires a simple transposition of the case in Fig. 2(b), yielding the contribution from K 0 to K as follows: PiT1 Ci,j ΨK 0 + 2
LAB
Ci,j Pj1 ΨK 0 . (36) 6 The cases of multi-level irregularly, as shown in Fig. 3, are treated recursively using the 1-level difference templates. For example, for case (a) of Fig. 3, we have LAB
LBC T T LCD T T PiT Ci,j ΨK 0 + Pi Pi Ci,j ΨK 0 + Pi Pi Ci,j ΨK 0 . (37) 2 12 11 6 1 6 2 1 6 2 2 Other edge-coupling terms, appearing in the diffusion bilinear form Eq. (22), are treated in a similar fashion with a few additional coupling matrix templates. 5. Results 5.1. Example 1 In this example, we compare the AMR and uniform mesh refinement convergence rates, as a function of the number of degrees of freedom (DOFs) and as a function of the CPU time, for a homogeneous medium placed in vacuum. The domain
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3185
Fig. 3. Multi-irregularity, i.e., differences > 1 in refinement levels.
size is [10 × 10] cm2 ; the total cross section σt is varied from 1 to 100 cm−1 to model a wide range of domain optical thicknesses (from 10 MFP to 1000 MFP). The scattering ratio c = σs /σ is equal to 0.9. A uniform and isotropic volumetric source is present in the domain. All computations are performed using an S4 level-symmetric angular quadrature [29]. Fig. 4 displays the convergence rates obtained using the two-mesh AMR technique for polynomial orders 1 through 4 and compares it with a uniform mesh refinement approach. The comparison is performed in terms of both DOFs and CPU time. For the AMR runs, the DOFs are the ones from the finer AMR mesh, i.e., Th/2 . A reference solution, obtained using the AMR technique with tighter tolerance, more adaptivity cycles and polynomial order 4, serves an exact solution to compute the error displayed in Fig. 4. First, we see a clear benefit in employing at least quadratic polynomials, regardless of the refinement strategy. For optically thin domains, AMR provides savings in memory and CPU time of about one order of magnitude when compared to the uniform mesh approach. For thicker domains, the savings with AMR are even larger. Fig. 5 shows the adapted mesh after 14 refinement cycles. Clearly, the S4 singularity lines have been refined (with S4 , there are 3 singularities originating from each corner of the domain; the 45◦ singularity from the bottomleft corner and the 135◦ singularity from the top-right corner are exactly aligned with the mesh and thus not visible here). 5.2. Example 2 The geometry and material descriptions are shown in Fig. 6. The polynomial order used is 2 and the angular quadrature is S4 ; scattering is isotropic. Table 1 presents the numerical spectral radius (NSR) obtained with employing DSA as a kui+1 −ui k
preconditioner for SI. The NSR is defined as the ratio of successive differences, kui −ui−1 k2 . For this problem, we can note 2 that the NSR is about 0.34 and decreases as the AMR proceeds. Even though the mesh size is increased at each AMR cycle, the total number of CG iterations in DSA per AMR cycle does not vary significantly. Most of the acceleration work occurs at the inception of the AMR process, when the meshes are coarse. Consequently the CPU time spent in DSA is small, even for highly refined meshes. For comparison purposes, we also present results in which the solution from the previous adaptivity cycle is discarded, hence no bootstrapping for the next adaptivity cycle is performed. In this case, the NSR remains about constant and the number of DSA iterations grows significantly: the total number of CG iterations doubles each time the number of unknowns is quadrupled. Should a large initial mesh be provided, a significant amount of work would be spent in DSA iterations. We note that the AMR process (with bootstrapping as one would always prefer to utilize available information) can be quite effective at focusing most of the DSA work when the AMR meshes are still relatively coarse as was noted by the decrease in the NSR when the adaptivity index increased. The AMR meshes obtained at adaptivity cycles # 6, 12 and 18 are given in Fig. 7 where we note that the interfaces and their associated S4 singularities are predominantly refined. 6. Conclusions A two-mesh Adaptive Mesh Refinement technique for the SN neutral-particle transport equation in 2D has been devised using a high-order (up to order 4) Discontinuous Galerkin (DG) finite element discretization with standard upwinding. At a given mesh adaptivity iteration, the difference in the solutions obtained on two AMR meshes, a coarser one and a finer one, is used to drive the AMR process. A Diffusion Synthetic Acceleration (DSA) scheme has been derived to precondition the transport solves. The DSA scheme, also based on a DG discretization, is obtained directly from the discretized highorder DG SN transport equation. The diffusion form upon which this DSA scheme is based is closely related to the standard Interior Penalty method for elliptic equations. The penalty stabilization coefficient was adjusted to yield unconditionally stable acceleration.
3186
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
Fig. 4. Convergence rates for various domain thicknesses and various polynomial orders using the two-mesh AMR technique (blue lines) or uniform mesh refinement (red lines); top row: domain size = 10 MFP, bottom row: domain size = 1000 MFP, left column: convergence rate vs. DOFs, right column: convergence rate vs. CPU time.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. AMR mesh after 14 adaptivity cycles, polynomial order 3, error achieved = 3.20 × 10−6 , example #1.
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3187
Fig. 6. Geometry and material description for example # 2.
Table 1 DSA-SI AMR Results, Example # 2. With bootstrapping Cycle ID
Number of active cells
NSR
Number of DSA iterations
CPU time in DSA (sec)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
50 74 104 122 206 278 356 500 794 1112 1526 2174 2942 3716 5144 6848 8096 10328 12746
0.3157 0.2437 0.3408 0.3663 0.3205 0.3400 0.3105 0.3193 0.2917 0.2790 0.2268 0.2786 0.2022 0.2051 0.2063 0.1978 0.2022 0.2038 0.1401
75 75 65 60 90 91 112 108 93 117 71 142 65 67 67 38 28 38 19
0.01 0.01 0.02 0.02 0.05 0.07 0.11 0.15 0.20 0.35 0.30 0.84 0.53 0.71 1.07 0.85 0.78 1.35 0.89
50 74 104 122 206 278 356 500 794 1112 1526 2174 2942 3716 5144 6848 8096 10328 12746
0.3157 0.3181 0.3108 0.3372 0.3046 0.3204 0.3411 0.3422 0.3403 0.3430 0.3964 0.3471 0.3465 0.3713 0.3454 0.3583 0.3715 0.3232 0.3830
75 128 133 151 238 282 318 367 464 546 637 782 890 1005 1124 1286 1546 1617 2021
0.01 0.03 0.04 0.05 0.14 0.21 0.30 0.51 1.00 1.64 2.61 4.61 7.12 10.30 17.62 27.96 40.71 56.29 88.13
With reinitialization 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
3188
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
Fig. 7. Example #2: AMR meshes at adaptivity cycles # 6, 12, and 18.
References [1] R.B. Kellogg, Numerical analysis of the neutron transport equations, in: B. Hubbard (Ed.), Numerical Solution of Partial Differential Equations- III, SYNSPADE 1975, Academic Press, Inc, New York, 1976. [2] James J. Duderstadt, William R. Martin, Transport Theory, John Wiley & Sons Inc, 1979. [3] R.E. Alcouffe, R.S. Baker, J.A. Dahl, S.A. Turner, PARTISN Code Abstract, Proc. Int. Topl. Mtg. Advances in Reactor Physics and Mathematics and Computations into the Next Millenium, Pittsburgh, Pennsylvania, May 7–12, 2000, American Nuclear Society, 2000. [4] Jim E. Morel, James S. Warsa, An SN spatial discretization scheme for tetrahedral meshes, Nucl. Sci. Engr. 151 (2005) 157–166. [5] Ekkehard Ramm Rannacher, E. Rank, R. Schweizerhof, K. Stein, E. Wendland, W. Wittum, G. Wriggers, Peter Wunderlich, Walter Stein (Eds.), Errorcontrolled Adaptive Finite Elements in Solid Mechanics, Wiley, 2001. [6] Tomasz Plewa, Timur Linde, V. Gregory Weirs (Eds.), Adaptive Mesh Refinement – Theory and Applications, Springer, 2005. [7] L. Demkowicz, Computing with hp-Adaptive Finite Elements. Vol.1: One and Two Dimensional Elliptic and Maxwell Problems, Chapman & Hall/CRC, 2007. [8] J.P. Jessee, W.A. Fiveland, L.H. Howell, P. Colella, R.B. Pember, An Adaptive Mesh Refinement Algorithm for the Radiative Transport Equation, J. Comput. Phys. 139 (2) (1998) 380–398. [9] R.M. Shagaliev, A.V. Alekseev, I.M. Beliakov, A.V. Gichuk, A.A. Nuzhdin, V.Yu. Rezchikov, Different algorithms of 2d transport equation parallelization on random non-orthogonal grids, in: Frank Graziani (Ed.), Computational Methods in Transport, Granlibakken 2004, in: Lecture Notes in Computational Science and Engineering, vol. 48, Springer, 2006, pp. 235–254. [10] F. Ogando, P. Velarde, Development of a radiation transport fluid dynamic code under amr scheme, J. Quant. Spectrosc. Radiat. Transf. 71 (2001) 541–550. [11] C. Aussourd, A multidimensional AMR SN scheme, Nucl. Sci. Engr. 143 (2003) 281–290. [12] Jose I. Duo, Yousry Y. Azmy, Ludmil T. Zikatanov, A posteriori error estimator and AMR for discrete ordinates nodal transport methods, Ann. Nucl. Energy 36 (2009) 268–273. [13] R.E. Alcouffe, Diffusion synthetic acceleration methods for the diamond differenced discrete-ordinates equations, Nucl. Sci. Engr. 64 (344) (1977). [14] E.W. Larsen, Unconditionally stable diffusion-synthetic acceleration methods for slab geometry discrete ordinates equations. part I, Theory. Nucl. Sci. Engr. 82 (47) (1982). [15] M.L. Adams, W.R. Martin, Diffusion-synthetic acceleration of discontinuous finite-element transport iterations, Nucl. Sci. Engr. 111 (1992) 145–167. [16] Marvin L. Adams, Edward W. Larsen, Fast iterative methods for discrete ordinates particle transport calculations., Prog. Nucl. Energy 40 (1) (2002) 31–59. [17] Robert C. Ward, Randal S. Baker, Jim E. Morel, A diffusion synthetic acceleration method for block adaptive mesh refinement, Nucl. Sci. Engr. 152 (2) (2006) 164–179. [18] J. Douglas, T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin method, in: Computing Methods in Applied Sciences, in: Lecture Notes in Physics, vol. 58, Springer-Verlag, 1976, pp. 207–216. [19] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, L. Donatella Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [20] P. Solin, K. Segeth, I. Dolezel, Higher-Order Finite Element Methods, Chapman & Hall / CRC Press, 2003. [21] Yaqi Wang, Jean Ragusa, A high-order discontinuous Galerkin method for the SN transport equations on 2D unstructured triangular meshes, Ann. Nucl. Energy 36 (2009) 931–939. [22] J.S. Warsa, T.A. Wareing, J.E. Morel, Krylov iterative methods and degraded effectiveness of diffusion synthetic acceleration for multidimensional sn calculations in problems with material discontinuities, Nucl. Sci. Eng. 147 (2004) 218–248. [23] Todd Arlin Wareing, Asymptotic diffusion accelerated discontinuous finite element methods for transport problems, Los Alamos National Laboratory LA-1 2425-T (1992). [24] Yaqi Wang, Adaptive mesh refinement solution techniques for the multigroup SN transport equation using a higher-order discontinuous finite element method, Ph.D. dissertation, Texas A&M University, May 2009. [25] G. Kanschat, Discontinuous Galerkin Methods for Viscous Flow, Deutscher Universitätsverlag, Wiesbaden, 2007. [26] S. Richling, E. Meinkohn, N. Kryzhevoi, G. Kanschat, Radiative transfer with finite elements: I. basic method and tests, A&A 380 (2001) 776–788. [27] P. Solin, L. Demkowicz, A goal-oriented hp-adaptivity for elliptic problems, Comput. Methods Appl. Mech. Eng. 193 (2004) 449–468. [28] Y. Wang, W. Bangerth, J. Ragusa, Three-dimensional h-adaptivity for the multigroup neutron diffusion equations, Progr. Nucl. Energy 51 (2009) 543–555. [29] Elmer E. Lewis, Warren F. Miller Jr., Computational Methods of Neutron Transport, John Wiley and Sons, Inc, New York, NY, 1984.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A two-mesh adaptive mesh refinement technique for SN neutral-particle transport using a higher-order DGFEM Jean C. Ragusa ∗ , Yaqi Wang Department of Nuclear Engineering, Texas A&M University, College Station, TX, 77843-3133, USA
article
abstract
info
Article history: Received 3 April 2009 Received in revised form 7 December 2009
We present a two-mesh Adaptive Mesh Refinement technique for the SN neutral-particle transport equation in 2D. A high-order (up to order 4) Discontinuous Galerkin (DG) finite element discretization is employed with standard upwinding. A suitable Diffusion Synthetic Acceleration (DSA) scheme is derived to precondition the transport solves. The DSA scheme, also based on a DG discretization, is derived directly from the discretized highorder DG SN transport equation. Results are presented to demonstrate the superiority of AMR for DSA-accelerated SN transport solves. © 2009 Elsevier B.V. All rights reserved.
Keywords: Adaptive mesh refinement SN transport Discontinuous Galerkin finite element method Diffusion synthetic acceleration
1. Introduction The steady-state transport of mono-energetic neutral particle (e.g., neutrons, photons) through matter is described by the following integro-differential equation
E + σ (Er ) Ψ = E ·∇ Ω
Z 4π
E0 ·Ω E )Ψ (Er , Ω E 0 )dΩ 0 + Q (Er , Ω E ), σs (Er , Ω
in D × S 2 ,
(1)
with the general boundary condition
E ) = Ψ inc (Erb , Ω E)+ Ψ (Erb , Ω
Z E 0 ·Enb >0 Ω
E0 → Ω E )Ψ (Erb , Ω E 0 )dΩ 0 on Erb ∈ ∂ D − , β(Erb , Ω
(2)
E the direction of flight, Er the position, σ (Er ) the total cross section, where DR is the spatial domain, S 2 the unit sphere, Ω E0·Ω E ) the scattering cross section, σs (Er , Ω E0·Ω E ) the differential scattering cross section, Q (Er , Ω E ) the σs (Er ) = 4π dΩ σs (Er , Ω E ) the incoming flux, and β(Erb , Ω E0 → Ω E ) the boundary redistribution kernel (albedo). The incoming source term, Ψ inc (E rb , Ω domain boundary ∂ D − is defined by E · nE(Erb ) < 0 , ∂ D − = Erb ∈ ∂ D , Ω
(3)
E(Erb ) the outward unit normal vector at position Erb . For brevity, we only consider Dirichlet conditions hereafter (i.e., with n β = 0). Despite Eq. (1) being linear, its numerical solution is challenging, even on current computers, due to 1. the high dimensionality of the phase-space. In complete generality, the angular flux Ψ depends upon seven variables: three for the standard physical space, two for the angular variable, one variable for the particles’ energy, and one for time. After the application of an implicit time integrator and the multigroup approximation in energy, an equation of the form given in Eq. (1) is obtained;
∗
Corresponding author. Tel.: +1 979 862 2033. E-mail addresses: [email protected] (J.C. Ragusa), [email protected] (Y. Wang).
0377-0427/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2009.12.020
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3179
2. the presence of singularities in the transport solution [1] and the greatly disparate values of σ (E r ) and σs (E r ) in practical applications; 3. the change in the nature of the PDE, from a hyperbolic behavior in materials with small scattering ratio, i.e., c = σσs 1, to an elliptic behavior for highly diffusive materials, i.e., c ' 1, [2]. In radiation transport production codes, low-order spatial differencing techniques are often employed, either on structured [3] or unstructured [4] grids. In addition, while Adaptive Mesh Refinement (AMR) techniques are now widely used in many science and engineering disciplines [5–7], only a limited number of references are available regarding the applications of mesh adaptivity to radiation transport. The higher-order spatial adaptivity technique described here is applied to the SN angular approximation of Eq. (1). Patch-AMR techniques [8–10] for the SN transport equation have been devised and are based on a hierarchy of nested grids. Some of the drawbacks of patch-AMR include the fact that the physics are not followed as closely as possible (the extent of refined patch being often too large), leading to more unknowns than needed, and the need to converge inflow/outflow values in between nested grids (a feature that is not present in cellbased AMR). Oftentimes, a gradient-based error estimator is employed to drive the adaptive mesh refinement in Cartesian structured geometries but it is overly conservative for high-order spatial discretization schemes. In [11], a local refinement (cell-based AMR) technique is described for SN transport spatially discretized using low-order finite differences, where the value of the particle mean free path (MFP) in a given cell is used as a mesh refinement criterion. While this approach takes into account the size of potential internal layers at a given location in the domain, it does not account for the actual smoothness of the solution at these locations and is, therefore, far from optimal; for example, within optically thick areas, the solution may well be approximated by a smooth spatial representation on coarse meshes despite the smallness of the mean free path. In [12], a lower-order angular flux solution serves to compute a converged scattering source which is then used as a fixed source term for a higher-order transport solve. The differences between the two solutions is then used to prescribe local mesh refinements; their technique was demonstrated for one-group, 2D rectangular structured meshes. The error estimate employed in our present work is based on the difference between the numerical solutions computed on coarser and finer AMR meshes at any adaptivity cycle; the local error hereby obtained is used to prescribe the new coarser AMR mesh at the next adaptivity cycle, the finer AMR mesh being obtained by refinement of the coarser mesh. Additionally, preconditioning techniques are necessary to accelerate the transport solution in highly diffusive configurations; it is well established that the spatial discretization of the Diffusion Synthetic Acceleration (DSA) equations must be ‘‘consistent’’ in some sense with the one used in the SN transport solves in order to yield unconditionally stable and efficient DSA schemes [13,14,4,15,16]. With the exception of [17] for block structured AMR with linear finite elements, little work has been reported to devise DSA schemes (i) for higher-order spatial discretizations and, more importantly, (ii) for unstructured AMR meshes. Using a variational argument, we obtain here a DGFEM DSA procedure (since it is based on a DG method, it is suitable for AMR meshes) and we show its relation with the standard Interior Penalty (IP) DGFEM method for elliptic PDEs [18,19]. The work presented here describes (i) an AMR technique to accurately solve the SN transport equation using a high-order Discontinuous Galerkin Finite Element Method (DGFEM) in space and (ii) an efficient DSA preconditioner for AMR meshes. Section 2 briefly reviews the SN approximation and describes the DGFEM variational form of the SN transport equation. In Section 3, we present our modified IP DSA preconditioner. The AMR technique and error estimator employed are explained in Section 4. Numerical results are given in Section 5 and we provide conclusions in Section 6. 2. Numerical solution techniques 2.1. Space/angle discretization In this section, we briefly review the SN angular discretization for Eq. (1) and present a DGFEM high-order spatial E m , wm }, m = 1, . . . , M, consisting of M directions Ω E m and weights wm , discretization. Given an angular quadrature set {Ω the SN discrete ordinates [2] version of Eq. (1) is Na X n X E + σ Ψm (Er ) = Em · ∇ E m ), Ω σs,n (Er )Φn,k (Er ) + Qn,k (Er ) An,k (Ω
(4)
n=0 k=−n
for m = 1, . . . , M, supplemented with appropriate boundary conditions. The scattering source term in Eq. (4) has been manipulated by decomposing the differential scattering cross section on Legendre polynomials (up to anisotropy order Na ) and introducing the flux moments
Φn,k (Er ) =
Z 4π
E )Ψ (Er , Ω E) ≈ dΩ An,k (Ω
M X
E m )Ψm (Er ), wm An,k (Ω
(5)
m=1
where An,k represents the real spherical harmonic function of degree n and order k with the Schmidt semi-normalization.
3180
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
The spatial domain D is triangulated into an unstructured mesh Th = {K , ∪K ∈Th K = D } and we employ a DG discretization technique in space, with Pp (K ) denoting the space of polynomial of degree ≤ p on element K . Hierarchical E m , the DG basis functions (up to order 4) are chosen for the spatial discretization [20,21]. For a given streaming direction Ω scheme with standard upwinding reads: Find Ψm ∈ Vp (K ), such that ∀Ψm∗ ∈ Vp (K ),
E + σ )Ψm∗ Em · ∇ Ψm , (−Ω
K
Na X n X
E m ) σs,n Φn,k + Qn,k , Ψm∗ , + Ψm− , Ψm∗ ∂ K + − Ψm− , Ψm∗ ∂ K − = An,k (Ω K
(6)
n=0 k=−n
where ∂ K − is the inflow boundary, ∂ K + is the outflow boundary of element K . The traces Ψ + and Ψ − are defined with E m , i.e., on an inflow boundary, Ψ + (resp. Ψ − ) is the value of Ψ taken from within element respect to the particle direction Ω K (resp. from the upwind neighbor) and on an outflow boundary, Ψ + (resp. Ψ − ) is the value of Ψ taken from the downwind neighbor (resp. from within element K ). These definitions are expressed as follows:
E m) Ψm± ≡ lim Ψm (Er + sΩ
(7)
s→0±
(f , g )K ≡
Z
hf , g ie ≡
Z
fg
(8)
E m · nE(Er )|f g |Ω
(9)
K
e
Em R 0 ∂ K ± ≡ Er ∈ ∂ K , nE(Er ) · Ω
(10)
E is the outward unit normal vector of element K and e represents any edge (or face in 3D) of K . Since the finite where n E is not defined at the element interfaces. Here, the upwind element approximation used is discontinuous, the quantity ΩΨ E) defined as scheme is obtaining by explicitly employing a numerical flux H (f + , f − , n E · nEf − E) = Ω H (f + , f − , n
(11)
such that
Z
E) Ψm∗ = Ψm− , Ψm∗ ∂ K + − Ψm− , Ψm∗ ∂ K − . H (Ψm+ , Ψm− , n
∂K
(12)
Eq. (12) was used in deriving Eq. (6). Transport solves are based on Eq. (6). By integrating Eq. (6) by parts, an equivalent formulation is obtained:
E + σ )Ψm , Ψm∗ Em · ∇ (Ω
K
Na X n X
E m ) σs,n Φn,k + Qn,k , Ψm∗ , + [[Ψm ]], Ψm∗ ∂ K − = An,k (Ω K
(13)
n=0 k=−n
where the inter-element jump is given by [[Ψm ]] = Ψm+ − Ψm− . Eq. (13) will be useful in deriving the diffusion synthetic accelerator in Section 3. 2.2. Iterative procedure The system of equations formed by Eq. (6) needs not to be assembled for all elements but can be solved element by E m . The order in which the elements are element. This strategy is often referred to as a ‘‘transport sweep’’ in direction Ω solved constitutes the graph of the sweep (for brevity, we do not expand on situations where the graph can present some dependencies, in which case these dependencies are simply lagged within the iterative procedure; for additional details, refer to [4]). With appropriate numbering, the angular fluxes and flux moment unknowns are described by 9 and 8, respectively, and the transport sweeps over all directions can be represented in the following matrix form: L9 = M Σ8 + Q 8 = DΨ ,
(14)
where L is the block-diagonal matrix due to particle streaming and total interaction (one block = one direction), 6 is the scattering matrix which operates on the flux moments, M is the moment-to-direction matrix, D is the direction-to-moment matrix (hereafter, we assume MD = I for conciseness in the derivations of Section 3, although not all quadrature formulae necessarily satisfy this). We can recast Eq. (14) in terms of flux moments as follows: I − DL−1 M Σ 8 = DL−1 Q.
(15)
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3181
None of these global matrices are explicitly formed and L is inverted in a ‘‘matrix-free’’ fashion using transport sweeps. Eq. (15) is solved using GMRes. Alternatively, a Richardson iteration is employed, leading to the scattering Source Iteration (SI) technique
9(i+1) = L−1 M68(i) + Q 8(i+1) = D9(i+1) ,
(16)
where the previous flux moment iterate 8(i) is utilized to build the scattering source term. Both the SI scheme and the GMRes solution are known to be slow converging for highly diffusive problems and a Diffusion Synthetic Acceleration (DSA) [16] must be implemented to address this issue. This DSA technique can be utilized as a preconditioner for both the SI scheme or the GMRes solves and is discussed next. 3. Diffusion preconditioner for AMR meshes For highly diffusive configurations (i.e., materials with scattering ratios c close to 1 and optically thick domains), the standard iterative techniques to solve the transport equation can become quite ineffective due to their slow convergence properties and a DSA preconditioner is needed to accelerate the convergence of the process. The diffusion operator must be spatially discretized in a manner that is ‘‘consistent’’ with the discretization of the transport equation [13–16] to ensure unconditional stability (however, a diffusion-based synthetic accelerator is not necessarily effective for multidimensional problems containing highly diffusive regions with discontinuities in material properties, [22]). To obtain consistently discretized DSA equations, the zeroth and first angular moments of the discretized transport equation are taken, leading to a mixed DGFEM discretization of a diffusion equation. The mixed formulation results in a nonsymmetric PD saddle point problem [4] that is CPU expensive and partially consistent DSA schemes, motivated by the difficulties associated with the algebraic elimination of the vector unknowns to yield an elliptic diffusion equations, are more popular. With partial consistency, it is hoped that the reduction in complexity outweighs the degradation in effectiveness but few partially consistent schemes have been applied to AMR meshes [17]. Additionally, finite-element-based partially consistent schemes employ a Continuous Finite Element (CFE) discretization of the diffusion operator [23]. However, the algorithmic and implementation simplicity of employing an AMR-DG solution for the transport solves would be severely reduced if the DSA preconditioner solves were based on a CFE method requiring hanging nodes treatment. In our new approach, we have chosen to derive partially consistent DSA schemes employing a DGFEM discretization by directly deriving them from the DGFEM discretization of the SN transport equations. Our scheme belongs to the family of partially consistent DSA methods because in our derivation, only the zeroth angular moment of the discretized transport equation is used and we assume that the continuous first angular moment of the transport equation (Fick’s law) to be verified in order to eliminate the current unknowns EJ = [Φ1,−1 ; Φ1,0 ; Φ1,1 ]T . The resulting DGFEM discretization for the diffusion equation is remarkably similar to the Interior Penalty (IP) stabilization method for diffusion equations solved using a DGFEM. Due to the discontinuous nature of the DGFEM approximation, it is particularly well suited for meshes arising in AMR calculations, i.e., hanging nodes are seamlessly incorporated into the DG method. This property also leads to an easy implementation of high-order test functions. For brevity, we only outline the DSA derivation here. Step 1: Using Eq. (13), we have the following variational form the DGFEM SN transport formulation: h Find Ψm ∈ WD , ∀Ψm∗ ∈ WDh , m = 1, . . . , M such that: b(Ψ , Ψ ∗ ) = l(Ψ ∗ ) +
N X n X
σs,n Φn,k , Φn∗,k
D
,
(17)
n=0 k=−n
where the bilinear and linear forms are b(Ψ , Ψ ∗ ) =
M X
E + σ )Ψm , Ψm∗ Em · ∇ wm ( Ω
D
m=1
l(Ψ ∗ ) =
N X n X n=0 k=−n
+
M X
wm [[Ψm ]], Ψm∗+ E i
m=1
Qn,k , Φn∗,k
D
+
X
X
(18)
h
wm Ψminc , Ψm∗ e .
(19)
e∈∂ D Ω E m ·Enb <0
h WD = Ψ ∈ L2 (D ); Ψ |K ∈ V (K ), ∀K ∈ Th and Ehi = ∪K ∈Th ∂ K \ ∂ D is the set of all interior edges. Note that Eq. (17) is the zeroth angular moment of the discretized transport equation, hence the weighted summation over all directions in Eqs. (18) and (19).
Step 2: Assume a linear dependence in angle (i.e., P1 approximation) for the angular flux Ψm and the test function Ψm∗
Ψm =
1 4π
1 E m ) and Ψm∗ = E m) (φ + 3EJ · Ω (φ ∗ + 3EJ ∗ · Ω 4π
(20)
3182
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
E and with φ ≡ Φ0,0 the scalar flux and EJ the net current. Further assume that the continuous Fick’s law (EJ = −D∇φ E ∗ ) to be valid, yielding: EJ ∗ = D∇φ Ψm =
1 4π
E ·Ω E m) (φ − 3D∇φ
and
Ψm∗ =
1 4π
E ∗·Ω E m ). (φ ∗ + 3D∇φ
(21)
Step 3: Insert Eq. (21) in Eq. (17) and carry out the angular summation. After some considerable algebra, the application of angular quadrature properties, and further simplifications [24], the following Modified Interior Penalty (MIP) DGFEM formulation is obtained:
E ∇φ E ∗ + κeMIP [[φ]], [[φ ∗ ]] i + [[φ]], { D∂n φ ∗ } i + D∇φ, Eh Eh D 1 1 + { D∂n φ}}, [[φ ∗ ]] E i + κeMIP φ, φ ∗ ∂ D − φ, D∂n φ ∗ ∂ D − D∂n φ, φ ∗ ∂ D h 2 2 E ∗ + J inc , φ ∗ E1 , 3D∇φ , lMIP (φ ∗ ) = Q0 , φ ∗ D + Q ∂D bMIP (φ, φ ∗ ) =
σa φ, φ ∗
D
(22) (23)
D
with the incoming current given by: J inc =
X
E m · nE(Erb )|Ψminc . wm |Ω
(24)
E m ·En(Erb )<0 Ω
[[]] and { } represent the inter-element jumps and averages on interior edges, respectively. The absorption cross section 1 σa = σ − σs and the diffusion coefficient D = 3(σ −σ are not imposed or coerced but are naturally obtained from this s,1 ) derivation. The penalty coefficient κeMIP in Eq. (22) is given by:
1 κeMIP = max κeIP ,
(25)
4
where the value of 1/4 is the value obtained by strictly applying the procedure, κeIP =
c (p+ ) D+ + 2 h⊥
+
c (p− ) D− − 2 h⊥
is the standard IP
penalty value (see, e.g., [25]); the ± symbols in κ represent the values in the adjacent cells K and K sharing a given edge ± e, D is the diffusion coefficient, h± ⊥ is the length of cell K orthogonal to the edge, c (p) = p(p + 1), and p is the polynomial order. We note that Eq. (22) is SPD and can, therefore, be solved using a standard preconditioned Conjugate Gradient (CG) solver; we used SSOR as the preconditioner for CG. A Fourier analysis of the iteration operator (DSA+transport solves) shows that our modified scheme (MIP) yields unconditional stability when employed as a diffusion synthetic accelerator, in both the thin and thick diffusion limits, whereas the standard IP scheme (κ = κ IP ) or the original scheme (κ = 1/4) are only conditionally stable accelerators. This Fourier analysis is discussed after the derivation of the MIP-DG scheme. Steps 1–3 describe how a discretized diffusion operator, that is partially consistent with the discretized transport equation, is obtained. Next, we briefly present how this MIP-DGFEM diffusion is employed to accelerate transport sweeps. For a given as a preconditioner SI iteration, one transport solve yields 8(i+1/2) = DL−1 M ΣΦ(i) + Q (this is simply Eq. (16) in which the new iterate has −
IP e
+
now been labeled i + 1/2); next we define the error e(i+1/2) = 8 − 8(i+1/2) where 8 represents the converged solution. This error satisfies the following transport equation I − DL−1 M6 e(i+1/2) = DL−1 M6 8(i+1/2) − 8(i) ,
(26)
= 6 8(i+1/2) − 8(i) . In the error equation, Eq. (26), the external volumetric source is
or, equivalently, (DLM − 6) e(i+1/2) no longer present, the Dirichlet boundary conditions are now homogeneous (vacuum) and the new volumetric source is the scattering source due to the difference in successive iterates. Nonetheless, this error equation is just as difficult to solve as Eq. (16) and a simpler estimate for the error term e˜ (i+1/2) is sought by replacing the (DLM − 6) operator with the diffusion operator obtained through steps 1–3, resulting in the following equation: Ae˜ (i+1/2) = 6 8(i+1/2) − 8(i) ,
(27)
where A is the matrix notation for the MIP-DGFEM bilinear form of Eq. (22). Finally, the DSA-accelerated SI scheme (or preconditioned Richardson iteration) reads as follows:
8(i+1) = 8(i+1/2) + e˜ (i+1/2) = DL−1 M ΣΦ(i) + Q + A−1 6 8(i+1/2) − 8(i) ,
(28)
and the preconditioned GMRes solve is I + A−1 6
I − DL−1 M6 8 = I + A−1 6 DL−1 Q.
(29)
Fig. 1 presents the results of a Fourier Analysis for a 2D homogeneous medium with periodic boundary conditions; the mesh consists of two triangles obtained by division of a single square; the scattering ratio is chosen to be c = 0.9999. It
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3183
Fig. 1. Fourier analysis for the various DSA forms as a function of the mesh optical thickness, homogeneous infinite medium case (the right pane is a zoomed portion of the left pane).
is well known that spectral radius of SI process is also equal to c, hence SI can be arbitrarily slow to converge. In Fig. 1, we have graphed the spectral radius obtained with the DSA-accelerated SI schemes when using the following three penalty coefficients: κ = 1/4, κ = κ IP (the standard IP formulation), and κ = κ MIP (the Modified IP formulation). The results are presented as a function of the mesh optical thickness (expressed in MFP). We can note that the schemes employing either κ = 1/4 or κ = κ IP can have spectral radii greater than 1, and hence are only conditionally stable for some range of optical thickness, whereas the MIP scheme has maximum spectral radius of about 0.48 and thus always converges. An S4 angular quadrature was used for this Fourier analysis. 4. A two-mesh AMR approach for high-order SN transport In radiation transport and radiative transfer, lower-order spatial schemes are most commonly used. The AMR techniques employed with such schemes typically rely on gradient-based error estimates [8], with a few exceptions, e.g., [26], where an adjoint-based estimator is derived, but at the cost of computing an adjoint solution. Here, we wish to obtain error estimates for high-order spatial approximations and, therefore, cannot rely on simpler estimates such as the gradient of the numerical solution. We have chosen to follow the approach proposed in [7,27], where two adapted AMR meshes are utilized at each adaptivity cycle: a coarse AMR mesh Tkh and a fine AMR mesh Tkh/2 obtained by uniform refinement of Tkh ; in the preceding, k denotes the mesh adaptivity cycle index. 4.1. Error estimation procedure The error estimate at cycle k is simply the difference in the solutions obtained on Tkh/2 and Tkh . This is computed in the L2 norm, element by element. However, in steady-state transport algorithms, the angular fluxes Ψm are never kept (except on albedo boundaries and for elements causing dependencies in the graph of the sweep), hence our error estimate is based on angle-integrated quantities such as the scalar flux φ or the current EJ = [J x ; J y ; J z ]T , as shown in Eq. (30) for 2D:
i R h x y y x 2 2 ( J − J ) + ( J − J ) h h h/2 h/2 K R x µkK = ∀K ∈ Tkh . 2 + (J y )2 ( J ) h h D
(30)
We tested error estimates based on either the scalar flux or the current and did not find significant changes in the AMR results and grids. The motivation to employ an estimate based on the current comes from prior work on AMR for the diffusion equation [28] where the diffusion approximation for the current is a mathematically appropriate choice. For efficiency, we first solve for the transport solution on the coarser mesh and use the converged scattering source as an initial guess for the solution on the finer mesh, hence reducing the number of transport iterations needed on the finer mesh. The projection operation in Eq. (30) is defined as: h , such that, Find Πh Φ ∈ WD
Πh Φ , Φ ∗
D
= Φ, Φ∗
D
,
∀Φ ∗ ∈ WDh
(31)
This projection is easily carried out element by element, due to the discontinuous nature of the spatial discretization. The criterion for refinement is defined as follows: an element K of Tkh is selected for refinement if
µkK ≥ α max µkK 0 ,
(32)
K 0 ∈Tkh
where α is a user-defined fraction (we used α = 0.3). This criterion allows us to focus the computational effort on elements with the largest errors and tends to equi-distribute the spatial error. The global error, obtained by summing the error estimates over all elements, i.e.,
ηk =
sX K ∈Tkh
µkK ,
(33)
3184
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
(a) Same refinement level.
(b) Coarser upwind element.
(c) Finer upwind element.
Fig. 2. Edge-coupling cases for refinement level differences ≤ 1.
is used to terminate the adaptivity procedure (ηk < tolAMR ). The final result of the calculation is the last solution obtained on the finer AMR mesh. 4.2. Mesh coupling Once an element has been flagged for refinement, it becomes inactive and the child-elements are the new active cells. Appropriate refinement rules allow for an efficient visit of the tree structure. Any level of refinement difference is allowed between neighboring elements. The element-coupling algorithm is based on recursive calls to the function dealing with the 1-level refinement difference. The cases of 0- and 1-level refinement difference are shown in Fig. 2 and are dealt with as follows. Fig. 2(a) is the standard situation without refinement between neighbors, where the edge-coupling from element K 0 to K is simply achieved by LAB
Ci,j ΨK 0 , (34) 6 where Ci,j is a coupling matrix. In 2D problems, Ci,j is akin to a 1D mass matrix (see Eq. (6)) but contains several empty rows and columns because it applies to the whole solution vector of element K 0 (only the components of ΨK 0 belonging to edge K ∩ K 0 play a role in edge-coupling). It was our deliberate choice to achieve edge-coupling through whole-element operations, as this is necessary for the diffusion solves (where the edge integrals present gradient terms and hence require the entire element solution). A few templates for the Ci,j matrices need to be stored, depending on the local numbering of vertices i and j; thanks to the hierarchical nature of the basis functions used, the templates for any polynomial order are embedded in the single template for the highest order employed. Fig. 2(b) shows the situation when the contribution to element K results from the contribution of two smaller elements (1-level refinement difference). The contribution is simply computed for each element separately, yielding LAB
LBC
(35) PiT Ci,j ΨK 0 , 1 6 6 2 where Pk are cell prolongation matrices operating on the local solution vector and providing the solution vectors on the four child-elements. These matrices are independent of the element shape. Finally, the case shown in Fig. 2(c) requires a simple transposition of the case in Fig. 2(b), yielding the contribution from K 0 to K as follows: PiT1 Ci,j ΨK 0 + 2
LAB
Ci,j Pj1 ΨK 0 . (36) 6 The cases of multi-level irregularly, as shown in Fig. 3, are treated recursively using the 1-level difference templates. For example, for case (a) of Fig. 3, we have LAB
LBC T T LCD T T PiT Ci,j ΨK 0 + Pi Pi Ci,j ΨK 0 + Pi Pi Ci,j ΨK 0 . (37) 2 12 11 6 1 6 2 1 6 2 2 Other edge-coupling terms, appearing in the diffusion bilinear form Eq. (22), are treated in a similar fashion with a few additional coupling matrix templates. 5. Results 5.1. Example 1 In this example, we compare the AMR and uniform mesh refinement convergence rates, as a function of the number of degrees of freedom (DOFs) and as a function of the CPU time, for a homogeneous medium placed in vacuum. The domain
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3185
Fig. 3. Multi-irregularity, i.e., differences > 1 in refinement levels.
size is [10 × 10] cm2 ; the total cross section σt is varied from 1 to 100 cm−1 to model a wide range of domain optical thicknesses (from 10 MFP to 1000 MFP). The scattering ratio c = σs /σ is equal to 0.9. A uniform and isotropic volumetric source is present in the domain. All computations are performed using an S4 level-symmetric angular quadrature [29]. Fig. 4 displays the convergence rates obtained using the two-mesh AMR technique for polynomial orders 1 through 4 and compares it with a uniform mesh refinement approach. The comparison is performed in terms of both DOFs and CPU time. For the AMR runs, the DOFs are the ones from the finer AMR mesh, i.e., Th/2 . A reference solution, obtained using the AMR technique with tighter tolerance, more adaptivity cycles and polynomial order 4, serves an exact solution to compute the error displayed in Fig. 4. First, we see a clear benefit in employing at least quadratic polynomials, regardless of the refinement strategy. For optically thin domains, AMR provides savings in memory and CPU time of about one order of magnitude when compared to the uniform mesh approach. For thicker domains, the savings with AMR are even larger. Fig. 5 shows the adapted mesh after 14 refinement cycles. Clearly, the S4 singularity lines have been refined (with S4 , there are 3 singularities originating from each corner of the domain; the 45◦ singularity from the bottomleft corner and the 135◦ singularity from the top-right corner are exactly aligned with the mesh and thus not visible here). 5.2. Example 2 The geometry and material descriptions are shown in Fig. 6. The polynomial order used is 2 and the angular quadrature is S4 ; scattering is isotropic. Table 1 presents the numerical spectral radius (NSR) obtained with employing DSA as a kui+1 −ui k
preconditioner for SI. The NSR is defined as the ratio of successive differences, kui −ui−1 k2 . For this problem, we can note 2 that the NSR is about 0.34 and decreases as the AMR proceeds. Even though the mesh size is increased at each AMR cycle, the total number of CG iterations in DSA per AMR cycle does not vary significantly. Most of the acceleration work occurs at the inception of the AMR process, when the meshes are coarse. Consequently the CPU time spent in DSA is small, even for highly refined meshes. For comparison purposes, we also present results in which the solution from the previous adaptivity cycle is discarded, hence no bootstrapping for the next adaptivity cycle is performed. In this case, the NSR remains about constant and the number of DSA iterations grows significantly: the total number of CG iterations doubles each time the number of unknowns is quadrupled. Should a large initial mesh be provided, a significant amount of work would be spent in DSA iterations. We note that the AMR process (with bootstrapping as one would always prefer to utilize available information) can be quite effective at focusing most of the DSA work when the AMR meshes are still relatively coarse as was noted by the decrease in the NSR when the adaptivity index increased. The AMR meshes obtained at adaptivity cycles # 6, 12 and 18 are given in Fig. 7 where we note that the interfaces and their associated S4 singularities are predominantly refined. 6. Conclusions A two-mesh Adaptive Mesh Refinement technique for the SN neutral-particle transport equation in 2D has been devised using a high-order (up to order 4) Discontinuous Galerkin (DG) finite element discretization with standard upwinding. At a given mesh adaptivity iteration, the difference in the solutions obtained on two AMR meshes, a coarser one and a finer one, is used to drive the AMR process. A Diffusion Synthetic Acceleration (DSA) scheme has been derived to precondition the transport solves. The DSA scheme, also based on a DG discretization, is obtained directly from the discretized highorder DG SN transport equation. The diffusion form upon which this DSA scheme is based is closely related to the standard Interior Penalty method for elliptic equations. The penalty stabilization coefficient was adjusted to yield unconditionally stable acceleration.
3186
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
Fig. 4. Convergence rates for various domain thicknesses and various polynomial orders using the two-mesh AMR technique (blue lines) or uniform mesh refinement (red lines); top row: domain size = 10 MFP, bottom row: domain size = 1000 MFP, left column: convergence rate vs. DOFs, right column: convergence rate vs. CPU time.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 5. AMR mesh after 14 adaptivity cycles, polynomial order 3, error achieved = 3.20 × 10−6 , example #1.
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
3187
Fig. 6. Geometry and material description for example # 2.
Table 1 DSA-SI AMR Results, Example # 2. With bootstrapping Cycle ID
Number of active cells
NSR
Number of DSA iterations
CPU time in DSA (sec)
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
50 74 104 122 206 278 356 500 794 1112 1526 2174 2942 3716 5144 6848 8096 10328 12746
0.3157 0.2437 0.3408 0.3663 0.3205 0.3400 0.3105 0.3193 0.2917 0.2790 0.2268 0.2786 0.2022 0.2051 0.2063 0.1978 0.2022 0.2038 0.1401
75 75 65 60 90 91 112 108 93 117 71 142 65 67 67 38 28 38 19
0.01 0.01 0.02 0.02 0.05 0.07 0.11 0.15 0.20 0.35 0.30 0.84 0.53 0.71 1.07 0.85 0.78 1.35 0.89
50 74 104 122 206 278 356 500 794 1112 1526 2174 2942 3716 5144 6848 8096 10328 12746
0.3157 0.3181 0.3108 0.3372 0.3046 0.3204 0.3411 0.3422 0.3403 0.3430 0.3964 0.3471 0.3465 0.3713 0.3454 0.3583 0.3715 0.3232 0.3830
75 128 133 151 238 282 318 367 464 546 637 782 890 1005 1124 1286 1546 1617 2021
0.01 0.03 0.04 0.05 0.14 0.21 0.30 0.51 1.00 1.64 2.61 4.61 7.12 10.30 17.62 27.96 40.71 56.29 88.13
With reinitialization 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
3188
J.C. Ragusa, Y. Wang / Journal of Computational and Applied Mathematics 233 (2010) 3178–3188
Fig. 7. Example #2: AMR meshes at adaptivity cycles # 6, 12, and 18.
References [1] R.B. Kellogg, Numerical analysis of the neutron transport equations, in: B. Hubbard (Ed.), Numerical Solution of Partial Differential Equations- III, SYNSPADE 1975, Academic Press, Inc, New York, 1976. [2] James J. Duderstadt, William R. Martin, Transport Theory, John Wiley & Sons Inc, 1979. [3] R.E. Alcouffe, R.S. Baker, J.A. Dahl, S.A. Turner, PARTISN Code Abstract, Proc. Int. Topl. Mtg. Advances in Reactor Physics and Mathematics and Computations into the Next Millenium, Pittsburgh, Pennsylvania, May 7–12, 2000, American Nuclear Society, 2000. [4] Jim E. Morel, James S. Warsa, An SN spatial discretization scheme for tetrahedral meshes, Nucl. Sci. Engr. 151 (2005) 157–166. [5] Ekkehard Ramm Rannacher, E. Rank, R. Schweizerhof, K. Stein, E. Wendland, W. Wittum, G. Wriggers, Peter Wunderlich, Walter Stein (Eds.), Errorcontrolled Adaptive Finite Elements in Solid Mechanics, Wiley, 2001. [6] Tomasz Plewa, Timur Linde, V. Gregory Weirs (Eds.), Adaptive Mesh Refinement – Theory and Applications, Springer, 2005. [7] L. Demkowicz, Computing with hp-Adaptive Finite Elements. Vol.1: One and Two Dimensional Elliptic and Maxwell Problems, Chapman & Hall/CRC, 2007. [8] J.P. Jessee, W.A. Fiveland, L.H. Howell, P. Colella, R.B. Pember, An Adaptive Mesh Refinement Algorithm for the Radiative Transport Equation, J. Comput. Phys. 139 (2) (1998) 380–398. [9] R.M. Shagaliev, A.V. Alekseev, I.M. Beliakov, A.V. Gichuk, A.A. Nuzhdin, V.Yu. Rezchikov, Different algorithms of 2d transport equation parallelization on random non-orthogonal grids, in: Frank Graziani (Ed.), Computational Methods in Transport, Granlibakken 2004, in: Lecture Notes in Computational Science and Engineering, vol. 48, Springer, 2006, pp. 235–254. [10] F. Ogando, P. Velarde, Development of a radiation transport fluid dynamic code under amr scheme, J. Quant. Spectrosc. Radiat. Transf. 71 (2001) 541–550. [11] C. Aussourd, A multidimensional AMR SN scheme, Nucl. Sci. Engr. 143 (2003) 281–290. [12] Jose I. Duo, Yousry Y. Azmy, Ludmil T. Zikatanov, A posteriori error estimator and AMR for discrete ordinates nodal transport methods, Ann. Nucl. Energy 36 (2009) 268–273. [13] R.E. Alcouffe, Diffusion synthetic acceleration methods for the diamond differenced discrete-ordinates equations, Nucl. Sci. Engr. 64 (344) (1977). [14] E.W. Larsen, Unconditionally stable diffusion-synthetic acceleration methods for slab geometry discrete ordinates equations. part I, Theory. Nucl. Sci. Engr. 82 (47) (1982). [15] M.L. Adams, W.R. Martin, Diffusion-synthetic acceleration of discontinuous finite-element transport iterations, Nucl. Sci. Engr. 111 (1992) 145–167. [16] Marvin L. Adams, Edward W. Larsen, Fast iterative methods for discrete ordinates particle transport calculations., Prog. Nucl. Energy 40 (1) (2002) 31–59. [17] Robert C. Ward, Randal S. Baker, Jim E. Morel, A diffusion synthetic acceleration method for block adaptive mesh refinement, Nucl. Sci. Engr. 152 (2) (2006) 164–179. [18] J. Douglas, T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin method, in: Computing Methods in Applied Sciences, in: Lecture Notes in Physics, vol. 58, Springer-Verlag, 1976, pp. 207–216. [19] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, L. Donatella Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (5) (2002) 1749–1779. [20] P. Solin, K. Segeth, I. Dolezel, Higher-Order Finite Element Methods, Chapman & Hall / CRC Press, 2003. [21] Yaqi Wang, Jean Ragusa, A high-order discontinuous Galerkin method for the SN transport equations on 2D unstructured triangular meshes, Ann. Nucl. Energy 36 (2009) 931–939. [22] J.S. Warsa, T.A. Wareing, J.E. Morel, Krylov iterative methods and degraded effectiveness of diffusion synthetic acceleration for multidimensional sn calculations in problems with material discontinuities, Nucl. Sci. Eng. 147 (2004) 218–248. [23] Todd Arlin Wareing, Asymptotic diffusion accelerated discontinuous finite element methods for transport problems, Los Alamos National Laboratory LA-1 2425-T (1992). [24] Yaqi Wang, Adaptive mesh refinement solution techniques for the multigroup SN transport equation using a higher-order discontinuous finite element method, Ph.D. dissertation, Texas A&M University, May 2009. [25] G. Kanschat, Discontinuous Galerkin Methods for Viscous Flow, Deutscher Universitätsverlag, Wiesbaden, 2007. [26] S. Richling, E. Meinkohn, N. Kryzhevoi, G. Kanschat, Radiative transfer with finite elements: I. basic method and tests, A&A 380 (2001) 776–788. [27] P. Solin, L. Demkowicz, A goal-oriented hp-adaptivity for elliptic problems, Comput. Methods Appl. Mech. Eng. 193 (2004) 449–468. [28] Y. Wang, W. Bangerth, J. Ragusa, Three-dimensional h-adaptivity for the multigroup neutron diffusion equations, Progr. Nucl. Energy 51 (2009) 543–555. [29] Elmer E. Lewis, Warren F. Miller Jr., Computational Methods of Neutron Transport, John Wiley and Sons, Inc, New York, NY, 1984.