Meshless discretization of the discrete-ordinates transport equation with integration based on Voronoi cells
MMeshless discretization of the discrete-ordinates transport equationwith integration based on Voronoi cells
Brody R. Bassett a,(cid:63) , J. Michael Owen a [email protected], [email protected] a Lawrence Livermore National Laboratory7000 East Avenue, Livermore, CA, 94550 (cid:63)
Corresponding author
Abstract
The time-dependent radiation transport equation is discretized using the meshless-local Petrov-Galerkinmethod with reproducing kernels. The integration is performed using a Voronoi tessellation, which creates apartition of unity that only depends on the position and extent of the kernels. The resolution of the integrationautomatically follows the particles and requires no manual adjustment. The discretization includes streamline-upwind Petrov-Galerkin stabilization to prevent oscillations and improve numerical conditioning. The angularquadrature is selectively refineable to increase angular resolution in chosen directions. The time discretizationis done using backward Euler. The transport solve for each direction and the solve for the scattering source areboth done using Krylov iterative methods. Results indicate first-order convergence in time and second-orderconvergence in space for linear reproducing kernels.
Keywords : radiation transport, meshless local Petrov-Galerkin (MLPG), streamline-upwind Petrov-Galerkin(SUPG), Voronoi tessellation 1 a r X i v : . [ phy s i c s . c o m p - ph ] A ug Introduction
Meshless methods have a rich variety of applications in hydrodynamic modeling, particularly in situations thatare challenging to mesh initially or where complex flows make it difficult to maintain a good meshed description.Examples include problems with large deformations, fractures, unstable flows, and a variety of astrophysical prob-lems (where such methods originated) [1]. However, many astrophysical problems in the high energy density physicsregime also require a thermal radiation transfer treatment [2], and to date there has been much less work on meshfreetreatments for radiation transport. One fundamental choice one must make is what sort of angular representationof the radiation is appropriate for the physics problem at hand. The few prior meshfree thermal radiative transfertreatments have generally focused on radiation diffusion and similar approximations [3, 4, 5, 6, 7], wherein theangular distribution of the radiation is neglected. This is appropriate for situations dominated by scattering andabsorption (such as deep in stellar interiors), but many interesting problems include both transparent and opaqueregions. To model these problems accurately, a more complicated angular discretization of the radiation transportequation is needed. The goal of this work is to develop a discrete ordinates radiation transport implementationthat is compatible with a method such as smoothed particle hydrodynamics [8] and variations of the reproducingkernel particle method [9, 10]. For this paper, the transport does not include radiation hydrodynamics effects andthe nonlinear emission terms in thermal radiative transfer.The discrete ordinates radiation transport equation has been solved previously using meshless methods, includingfor collocation methods [11, 12, 13, 14, 15, 16] and the meshless local Petrov-Galerkin (MLPG) method [17, 18]. Thecollocation methods often use the second-order form of the radiation transport equation, such as the self-adjointangular flux equation [19], while the MLPG methods use a background mesh for the integration.This paper extends the past MLPG discretizations of the radiation transport equation in several ways. First, re-producing kernels [9] are used instead of moving least squares, which reduces the number of linear solves needed toevaluate the kernels. Using these RK functions, higher-than-second-order convergence is demonstrated for certainchoices of RK correction order. Second, the implementation is time-dependent, which adds additional complexityand means that the integration, which for for some problems may be performed at each time step, cannot requiremanual adjustment and needs to be efficient. Third, the integration is done using a Voronoi decomposition, whichmeets both of these criteria. Fourth, the angular quadrature can be refined, which permits consideration of prob-lems that require high angular resolution in specific directions. Finally, the implementation of the code has beendone inside a code that already includes radiation hydrodynamics with diffusion [20, 21] and reproducing kernelhydrodynamics [10], which should permit future consideration of meshless radiation hydrodynamics.In many cases, including the meshless Galerkin approach [22], the integration is done using a background mesh [23].In the original MLPG paper, the authors recommend using circular (2D) or spherical (3D) domains of integrationto retain a truly meshless method [24]. Integration can also be performed by introducing a quadrature into thelens-shaped intersection of two kernels [25] or by reducing the dimensionality of the integrals [26]. One method ofavoiding any evaluations outside of the kernel centers is by using nodal integration, which may require stabilizationfor the derivatives [27] and cannot be further refined to increase accuracy [28], although using a Voronoi tessellationto create the volumes can increase accuracy [29]. By using a quadrature within a Voronoi tessellation to integratethe kernels, the resolution of the integrals follows the resolution of the particles and that the integration meshdepends only on the location of the particles and the boundary surfaces of the problem. This, in effect, makes theintegration invariant to the rotation or translation of the points within the domain.The paper is structured as follows. In Sec. 2, the smoothed particle hydrodynamics (SPH) and reproducing kernels(RK) are introduced, along with methods for calculating derivatives. Next, the integration method using a Voronoitessellation is discussed in Sec. 3 and transformations from local to global coordinates are derived in App. A. Thetime-dependent transport equation is introduced and discretized using these kernels and integration methods inSec. 4. The discretization is then tested in Sec. 5 for a purely absorbing problem, two manufactured solutions, apurely-scattering problem with disparate cross sections, and a problem with a source in a void far from a strongabsorber whose solution is derived in App. B. Finally, conclusions and future work are discussed in Sec. 6.
This section introduces the smoothed particle hydrodynamics (SPH) and reproducing kernel (RK) functions. Hereand throughout this paper, Greek superscripts represent dimensional indices, with repeated indices representingsummation, e.g. x α is the component α of the vector x and X αβ is the component α, β of the tensor X . The2otation ∂ αx denotes the partial derivative with respect to x α . Subscripts denote evaluation of a function at adiscrete point, e.g. f j = f ( x j ) , unless otherwise noted. In SPH, spatial fields are interpolated from the discrete values at individual points using functions referred to askernels, which are symmetric functions (typically of vector distance) with compact support, i.e., they fall to zero atsome range. The support of the kernel is determined by a smoothing parameter, which translates from physical spaceto the reference space for the kernel. This smoothing parameter can either be a scalar h , as in standard SPH, or asymmetric tensor H αβ , as in adaptive smoothed particle hydrodynamics (ASPH) [30]. This smoothing parameteris generally allowed to vary point to point, so H αβi in general. Note that when using SPH, the interpolation kernel isradially symmetric, while under ASPH this is not necessarily true. Because the ASPH H αβ tensor is symmetric, thecorresponding smoothing scale isocontours around a point are elliptical (2D) or ellipsoidal (3D). Using the tensorsmoothing parameter (which has units of inverse length), the transformed distance vector in ASPH reference space η ( x ) and the scaled distance χ ( η ) are η α = H αβ x β , (1) χ = √ η α η α . (2)Note that SPH is a simply a special case of ASPH, wherein H αβ = h − δ αβ , so in the SPH case Eq. (1) reduces to η α = x α /h. (3)With these conventions in mind the kernel equations can be written in terms of H αβ and apply equally to SPH orASPH. The kernel W ( x ) and its derivatives can be defined in terms of the base kernel in reference space W b ( χ ) as W = W b , (4) ∂ γx W = ∂ αη χ∂ γx η α ∂ χ W b . (5)The derivative equation can be simplified by inserting the derivatives of χ and η , ∂ αη χ = η α χ , (6) ∂ γx η α = H αγ , (7)which results in ∂ γx W = η α χ H αγ ∂ χ W b . (8)As the kernels approximate delta functions ( W ( x ) → δ ( x ) as H αβ → ∞ ), they can be used in interpolation, f ( x ) = (cid:90) V ∂ ( x − x (cid:48) ) f ( x (cid:48) ) dV (cid:48) ≈ (cid:90) V W ( x − x (cid:48) ) f ( x (cid:48) ) dV (cid:48) . (9)The kernels are normalized such that formally the volume integral is unity, (cid:90) V W ( x − x (cid:48) ) dV (cid:48) = 1 , (10)though this property is only approximately true in the discrete case for SPH. The interpolant can be discretized fora set of these kernels with discrete positions x j with associated volumes V j , f ( x ) = (cid:88) j V j W j ( x ) f j , (11)where W i ( x ) = W ( x − x i ) (12)and f denotes the discrete interpolant of f . 3 .2 Introduction to reproducing kernels In general, the standard SPH kernels cannot reproduce even a constant solution exactly, (cid:88) j V j W j ( x ) (cid:54) = const . (13)The SPH kernels can be augmented with RK functions [9], which permit exact interpolation of functions up to acertain polynomial order. Interpolation with RK functions U i ( x ) works the same as in SPH, f ( x ) = (cid:88) j V j U j f j , (14)with the caveat that the RK functions have the property that (cid:88) j V j U j ( x ) f j = f ( x ) , f ( x ) ∈ P n (15)where P n is the space of polynomials with degree less than or equal to n . These functions and their derivatives aredefined in terms of the SPH functions as U i = P (cid:62) i CW i , (16) ∂ γ U i = (cid:0) ∂ γx P (cid:62) i C + P i ∂ γx C (cid:1) W i + P (cid:62) i C∂ γx W i , (17)where P ( x ) is the polynomial basis vector, P ( x ) = (cid:2) , x α , x α x β , · · · (cid:3) (cid:62) , (18) P i = P ( x − x i ) , (19) P ij = P ( x j − x i ) , (20)and C ( x ) = (cid:2) C ( x ) , C ( x ) , C ( x ) , · · · (cid:3) (cid:62) is a corrections vector of the same size as the polynomial vector withcoefficients C k (i.e. the component k of C ) to be determined. Suppose that F is a vector of arbitrary coefficients.The RK method calculates C such that the reproducing kernels can exactly represent F (cid:62) P . The term F (cid:62) P i isinterpolated as F (cid:62) P i = (cid:88) j V j F (cid:62) P ji U j = F (cid:62) (cid:88) j V j P ji P (cid:62) j CW j . (21)This equation must be true for each component of F , P i = (cid:88) j V j P ji P (cid:62) j CW j , (22)and can be evaluated at the point x i to produce simple conditions for the coefficients, (cid:88) j V j P ji P (cid:62) ji W ji C i = G, (23)where W ji = W j ( x i ) and G = [1 , , , · · · ] (cid:62) . The matrix for this linear system and its derivatives can be writtenexplicitly as M i = (cid:88) j V j P ji P (cid:62) ji W ji , (24) ∂ γ M i = (cid:88) j V j (cid:2)(cid:0) ∂ γ P ji P (cid:62) ji + P ji ∂ γ P (cid:62) ji (cid:1) W ji + P ji P (cid:62) ji ∂ γ W ji (cid:3) . (25)4n terms of these matrices, the linear systems to solve for C i and its derivatives can be written as M i C i = G, (26) ∂ γ M i C i + M i ∂ γ C i = 0 . (27)By first solving for C i and then ∂ γ C i , the only matrix that needs to be inverted is M i , C i = M − i G, (28) ∂ γ C i = − M − i ∂ γ M i C i , (29)which lets us reuse its factorization. In this section, the methodology for creating a Voronoi tessellation is introduced, the meshless integration process isdescribed, and the connectivity for the weak-form kernels is derived in terms of a similar strong-form connectivity.
As discussed in the introduction (Sec. 1), there have been several methods developed to integrate radial basisfunctions. For these results, the problem is decomposed using what is essentially a Voronoi tessellation constructedusing the PolyClipper library [31], with one line segment (1D), polygon (2D), or polyhedron (3D) per meshfreepoint. Each cell is then further decomposed into triangles (2D) or tetrahedra (3D), with surfaces defined by points(1D), line segments (2D) or triangles (3D). It is worth pointing out that the decomposition is not truly the Voronoi.Rather the decomposition begins with an initial polytope for each point that encompasses the finite kernel extent ofthat point (i.e., the space over which its kernel value is non-zero), which is progressively clipped by planes halfwaybetween the point in question and each neighbor point it interacts with. In the end, this results in a tiling of spacewith these polytopes per point that exactly constructs a partition of unity in space for all points that overlap.Note that the topological connection between the polytopes for each point is not computed, but only a uniquepolygon or polyhedron for each point independently. Figure 1a shows a cartoon of this process. The goal is toconstruct the Voronoi-like polygon for the central red point, which has a set of neighbor points it overlaps (in blue),and a non-zero kernel extent represented as the gray region. The starting polygon for this point is the boundingsurface of this gray region, which is progressively clipped by planes half-way between the central red point and eachof its neighbors. In the end all that is left is the central light red polygon, which is the unique volume closer tothe red point than any of its neighbors. To facilitate simple integration quadratures, these polytopes for each pointare further broken down into triangles (in 2D) and tetrahedra (in 3D). Fig. 2 shows this procedure for the polygongenerated in Fig. 1, where the cross markers denote the centroid of each of the sub-triangles in the polygon. Notethe centroid of the polygon does not necessarily coincide with the original point used to construct it.
The set of integration quadratures over the base shapes represents a contiguous, non-overlapping quadrature thatcovers the domain. A Gauss-Legendre quadrature is used for integration of the line segments, while symmetricquadrature rules as described in Ref. [32] are used to integrate the triangles and tetrahedra. Appendix A presentsinformation on how the integrals are transformed from physical to reference space.At each integration point, all the functions whose support includes the integration point must be evaluated. TheRK functions (Sec. 2.2) are expensive to evaluate relative to a standard SPH kernel and the evaluation of one suchfunction depends on the values of all other functions at that point. As such, a large amount of computation can besaved by making the quadrature the outermost loop in the code, as shown in Alg. 1. For each quadrature point,all functions whose support includes the integration quadrature point are evaluated. Then these values are used toperform each integral.It is worth comparing this approach with prior background integration methodologies, wherein a traditional back-ground mesh (often some sort of orthogonal Cartesian grid aligned with the lab frame) is placed independently ofthe meshless points, even if the point locations inform characteristics of the background mesh. The approach in this5ection, using a Voronoi tessellation for the integration, produces an integration mesh whose properties are solelydetermined by the volume and relative positioning of the meshless points, which makes it invariant to rotation ortranslation. This integration does not meet the strictest of meshless criteria, in which no mesh is allowed [33], butdoes meet a looser criterion in that all information can be derived directly from the meshless points. The geometryof each point’s unique volume is constructed based solely on the positions of surrounding points without storing,evolving, or specifying anything except the the point positions and kernel extents.
For the following sections, to simplify the notation, volume and surface integrals will be written as (cid:104) f, g (cid:105) = (cid:90) V f gdV, (30) ( f, g ) = (cid:90) S f gdS. (31)For many applications, such as hydrodynamics and diffusion, SPH and RK can be used to directly discretize theequations via collocation [8]. The equation, in this example a α ∂ αx f + bf = 0 , (32)is first integrated by parts, ( U i , n α a α f ) − (cid:104) ∂ αx U i , a α f (cid:105) + (cid:104) U i , bf (cid:105) = 0 , (33)(with n denoting the surface normal), the surface term is discarded, and interpolants [Eq. (14)] and the deltafunction property [Eq. (9)] are used to simplify the equation to − (cid:88) j V j a αj f j ∂ αx i U ji + b i f i = 0 . (34)While this form of the equation has the advantage of simplicity, it depends on a one-point quadrature rule for theintegration, (cid:90) V U i f ( x ) ≈ f i , (35)and generally throws away surface terms. Another option, and the one used in this paper, is to insert a basisfunction expansion, f ( x ) = (cid:88) j V j U j g j (36)(with coefficients g j ), (cid:88) j V j [( U i , n α a α U j ) − (cid:104) ∂ αx U i , a α U j (cid:105) + (cid:104) U i , bU j (cid:105) ] g j = 0 , (37)and perform the integrals directly using a quadrature like the one described in Sec. 3.2, (cid:104) ∂ αx U i , a α U j (cid:105) ≈ (cid:88) k w k ∂ αx k U ik a αk U jk , (38)where w k are the weights of a quadrature spanning the integration volume. For functions with compact support,this is the Meshless-Local Petrov-Galerkin (MLPG) method. Two types of connectivity are used in evaluating and storing the MLPG integrals. In a standard SPH code, theconnectivity is the sets of points whose evaluation is nonzero at the center of the other, so points i and j areneighbors if the support of W i includes the point x j or vice versa, or χ ( η ( x i − x j )) ≤ r, (39)6here r is the dimensionless support radius of the kernel. For MLPG, the connectivity (or sets of points for whichthe bilinear integrals are nonzero) is determined by whether points have overlapping support, so points i and j areneighbors if the support regions of W i and W j intersect. The overlap connectivity is a subset of the set of pointsfor which χ ( η ( x i − x j )) ≤ r. (40)To show this, suppose that there is a point x k that is in the support radius of W i and W j . It follows from Eq. (39)and the triangle inequality in Euclidean space that χ ( η ( x i − x j )) ≤ χ ( η ( x i − x k )) + χ ( η ( x k − x j )) ≤ r .Standard SPH connectivity information can be used to both create the overlap connectivity and calculate whichfunctions are nonzero at each integration point, which is similar to standard SPH connectivity. If the radius ofthe MLPG kernels is doubled (or the smoothing length altered to produce a similar effect), then by Eqs. (39) and(40), the standard SPH connectivity of the doubled-radius will include all overlap neighbors. As each cell producedby the Voronoi tessellation is completely contained within the support of its associated MLPG point, this samedouble-radius SPH connectivity for i and j will include all integration points in the cell i for which the kernel W j is nonzero. This is why the same connectivity can be used for both the overlap of two kernels and the overlap of akernel with a Voronoi cell associated with a kernel in Alg. 1. In this section, the radiation transport equation is introduced and discretized using MLPG with RK functions.Then, the iterative methods that are used for the solution of the discretized equation are described. Finally, theangular quadrature with selective refinement is introduced.
The gray radiation transport equation, which is the transport equation integrated over all energies, is ∂ t ψ + Ω α ∂ αx ψ + σ t ψ = 14 π σ s φ + q, (41)with the boundary condition ψ = ψ b , x ∈ ∂V, Ω α n α > , (42)and initial condition ψ = ψ init , t = 0 , where Ω is the radiation propagation direction, ψ is the angular flux, φ = (cid:82) π ψd Ω is the scalar flux, σ s is thescattering cross section, σ a is the absorption cross section, σ t = σ a + σ s is the total cross section, q is a source thatmay include physics such as thermal emission, ψ b is the incoming flux at the boundary, ψ init is the initial angularflux, and ∂V denotes the boundary of the domain. The discrete-ordinates approximation evaluates this equationat discrete angles that are ordinates of a quadrature over a unit sphere. Denoting these ordinates as Ω m for theangular index m , integrals over the unit sphere become (cid:90) π f d Ω ≈ (cid:88) m w m f m , (43)where w m are the weights of the quadrature. The transport equation evaluated at the discrete ordinate m becomes ∂ t ψ m + Ω αm ∂ αx ψ m + σ t ψ m = 14 π σ s φ + q m , (44)with the scalar flux φ = (cid:80) m w m ψ m . The backward Euler (or fully-implicit) method is used to discretize in time, c ∆ t ψ m + Ω αm ∂ αx ψ m + σ t ψ m = 1 c ∆ t ψ n − m + 14 π σ s φ + q m , (45)where all variables are evaluated at time index n except where noted otherwise as a superscript.A standard Galerkin approach to transport would be to multiply the transport equation by U i and then integrateover the support of U i . For streamline-upwind Petrov-Galerkin (SUPG) stabilization, the transport equation is7nstead multiplied by U i + τ i Ω α ∂ αx U i , where τ is a proportionality constant with unit length that is chosen to beconstant for each trial function. Performing this operation, the transport equation becomes c ∆ t (cid:104) U i , ψ m (cid:105) + 1 c ∆ t τ i Ω αm (cid:104) ∂ αx U i , ψ m (cid:105) + ( U i , Ω αm n α ψ m ) Ω α n α > − Ω αm (cid:104) ∂ αx U i , ψ m (cid:105) + τ i Ω αm Ω βm (cid:10) ∂ αx U i , ∂ βx ψ m (cid:11) + (cid:104) U i , σ t ψ m (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , σ t ψ m (cid:105) = 1 c ∆ t (cid:10) U i , ψ n − m (cid:11) + 1 c ∆ t τ i Ω αm (cid:10) ∂ αx U i , ψ n − m (cid:11) + (cid:0) U i , | Ω αm n α | ψ bm (cid:1) Ω α n α < + 14 π (cid:104) U i , σ s φ (cid:105) + 14 π τ i Ω αm (cid:104) ∂ αx U i , σ s φ (cid:105) + (cid:104) U i , q m (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , q m (cid:105) . (46)The surface integral term produced by integration by parts has been split into known (incoming) and unknown(outgoing) parts and the boundary condition [Eq. (42)] has been applied. Inserting a basis function expansion forthe scalar and angular flux [as in Eq. (36)], φ = (cid:88) j V j U j Φ j , (47a) ψ = (cid:88) j V j U j Ψ j , (47b)the equation becomes (cid:88) j V j (cid:20) c ∆ t (cid:104) U i , U j (cid:105) + 1 c ∆ t τ i Ω αm (cid:104) ∂ αx U i , U j (cid:105) + ( U i , Ω αm n α U j ) Ω α n α > − Ω αm (cid:104) ∂ αx U i , U j (cid:105) + τ i Ω αm Ω βm (cid:10) ∂ αx U i , ∂ βx U j (cid:11) + (cid:104) U i , σ t U j (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , σ t U j (cid:105) (cid:3) Ψ m,j = (cid:88) j V j (cid:20) c ∆ t (cid:104) U i , U j (cid:105) + 1 c ∆ t τ i Ω αm (cid:104) ∂ αx U i , U j (cid:105) (cid:21) Ψ n − m,j + (cid:0) U i , | Ω αm n α | ψ bm (cid:1) Ω α n α < + 14 π (cid:88) j V j [ (cid:104) U i , σ s U j (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , σ s U j (cid:105) ] Φ j + (cid:104) U i , q m (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , q m (cid:105) . (48)Note that, in general, Φ i (cid:54) = φ i and Ψ i (cid:54) = ψ i . Once Eq. (48) is solved for Φ and Ψ , the solution at each point mustbe recovered using Eqs. (47).The effect of the SUPG stabilization is to reduce oscillations and make the system of equations easier to solveiteratively while not affecting global particle balance [18]. For the results in Sec. 5, the stabilization parameter isset to be τ i = h i k , where k is approximately the number of points across the kernel radius. This results in a τ i that is approximatelyequal to the spacing between the points.Because RK permits interpolation, the initial condition can be set using initial values of the angular flux instead ofneeding to interpolate coefficients, or Ψ i (cid:12)(cid:12)(cid:12)(cid:12) t =0 = ψ init (cid:12)(cid:12)(cid:12)(cid:12) x = x i . This is what is done for the results in Sec. 5. The transport equation can be written in operator form as L Ψ = T Ψ n − + MS Φ + r, (49)or in terms of Φ as (cid:0) I − DL − MS (cid:1) Φ = DL − (cid:0) T Ψ n − + r (cid:1) , (50)8ith the operators defined as ( L Ψ) m,i = (cid:88) j V j (cid:20) c ∆ t (cid:104) U i , U j (cid:105) + 1 c ∆ t τ i Ω αm (cid:104) ∂ αx U i , U j (cid:105) + ( U i , Ω αm n α U j ) Ω α n α > − Ω αm (cid:104) ∂ αx U i , U j (cid:105) + τ i Ω αm Ω βm (cid:10) ∂ αx U i , ∂ βx U j (cid:11) + (cid:104) U i , σ t U j (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , σ t U j (cid:105) (cid:3) Ψ m,j , (51) (cid:0) T Ψ n − (cid:1) m,i = (cid:88) j V j (cid:20) c ∆ t (cid:104) U i , U j (cid:105) + 1 c ∆ t τ i Ω αm (cid:104) ∂ αx U i , U j (cid:105) (cid:21) Ψ n − m,j , (52) ( MS Φ) i,m = 14 π (cid:88) j V j [ (cid:104) U i , σ s U j (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , σ s U j (cid:105) ] Φ j , (53) ( D Ψ) i = (cid:88) m w m Ψ m,i , (54) ( r ) i,m = (cid:0) U i , | Ω αm n α | ψ bm (cid:1) Ω α n α < + (cid:104) U i , q m (cid:105) + τ i Ω αm (cid:104) ∂ αx U i , q m (cid:105) . (55)Note that using this notation, Φ = D Ψ . The notation L − denotes the linear inverse of the L operator, which isblock diagonal in angle. With the first-flight source defined as b ψ = L − (cid:0) T Ψ n − + r (cid:1) , (56) b φ = D b ψ , (57)the equation can be simplified to (cid:0) I − DL − MS (cid:1) Φ = b φ . (58)Equation (58) can be solved directly using a matrix-free linear solver such as GMRES or iteratively using a methodsuch as fixed point iteration. Once the scattering source is converged, the angular flux is recovered for use in thenext time step by performing an additional solve using the converged scalar flux, Ψ = L − MS Φ + b ψ . (59)For the results in Sec. 5, the L − operation is performed using two packages from Trilinos [34], the Belos packagefor GMRES and the Ifpack2 package for the ILUT preconditioner. The ILUT factorizations for each angle areprecomputed at the start of the time step and reused to minimize computation. The iterations to converge thescattering source [the solution of Eq. (58)] are also performed using GMRES from Belos without a preconditioner.For a discussion on preconditioners for the scattering iterations, see Sec. 6. For the problems in Sec. 5, the Gauss-Legendre quadrature is used in 1D, while the LDFE (linear discontinuousfinite element) quadrature [35] is used in 2D and 3D. The LDFE quadrature is hierarchal, meaning that each octantof the unit sphere can be further subdivided into four ordinates, which can themselves be subdivided and so on.This can be used to produce a high density of angular ordinates in chosen directions object to prevent ray effects,which as shown in Sec. 5.4. Given a goal quadrature, the refined angular discretization keeps all ordinates from thegoal quadrature that hit the object and combines the ordinates that do not hit the object inasmuch as is possible(Alg. 2). This significantly reduces the number of angles needed for a given number of rays from a small source tohit a distant object.
In this section, five problems are considered to test the discretization described in Sec. 4 with the integration inSec. 3. The first two problems use the method of manufactured solutions. The third problem considers a purelyabsorbing medium, while the fourth problem considers a purely scattering medium. The final problem shows apossible application of the code to simulate an asteroid absorbing a large quantity of radiation from a distantsource. All the problems use a kernel sampling radius of 4 neighbors (i.e., the equivalent smoothing scale is 4times the local particle spacing) and an RK order of one, as described in Sec. 3, except for the purely absorbingproblem, which also explores other combinations of RK order and neighbors. Note that for uniformly spaced points,9 sampling radius of 4 neighbors implies a total number of overlap neighbors for each point of 8 (1D), 50 (2D), and268 (3D).The first two sections use the relative error as a measure for convergence. This is defined as (cid:15) rel = (cid:80) i (cid:12)(cid:12)(cid:12) φ numeric i − φ analytic i (cid:12)(cid:12)(cid:12)(cid:80) i φ analytic i (60)for the numeric and analytic scalar fluxes at the MLPG centers i . The method of manufactured solutions works by selecting a solution for ψ , solving for a source q by insertingthis solution into the continuous transport equation [Eq. (41)], assigning the boundary source ψ b to be equalto the solution, and then calculating a numerical solution using the discretized transport equation [Eq. (48)] forcomparison to the original solution.The spatial and time convergence of the manufactured problems is considered in 1D, 2D, and 3D. In 1D, dueto the low cost of integration and because the integration cells are not subdivided, the integration is performedusing a 64-point Gauss-Legendre quadrature. In 2D, the integration of the subcell triangles is performed usinga tenth-order symmetric quadrature with 25 points, while in 3D, the integration of the surface triangles and thesubcell tetrahedra is performed using a third-order symmetric quadrature with 8 points. For more information onthe integration quadratures, see Sec. 3.The spatially-dependent results are run for several cases between d and d points (for the dimension d ) withouttime dependence. For the time-dependent case, the simulation is run until t = 1 with time steps between 0.001 and1.0 and d points. The minimum time step is increased to 0.01 in 2D and 0.1 in 3D. In each case the points arelaid down in a spatially uniform lattice configuration. For the non-uniform cases, the point positions are randomlyperturbed by up to 0.2 times the point distance in each dimension, or x α = x α ± ∆ x α γ α , where − . ≤ γ α ≤ . is randomly generated for each point and dimension independently and ∆ x α is the point spacing for the givendimension. The first manufactured solution, ψ sinusoidal = 1 + 12 π (cid:89) α cos ( π ( x α + t )) , (61)is designed to test convergence of the discretized transport equation [Eq. (48)]. The solution is chosen such thatthe manufactured source never becomes negative, which would be unphysical. To ensure that the integration of thecross sections [Sec. 3] works correctly, the scattering and absorption opacities are also chosen to have sinusoidalvalues that are out of phase with the solution and one another, σ a = 1 + 23 (cid:89) α cos (3 x α ) , (62) σ s = 1 + 34 (cid:89) α sin (2 x α ) . (63)The domain is − ≤ x α ≤ . For the steady-state case, the manufactured solution is fixed at t = 0 .The spatial convergence results in Fig. 3a indicate second-order convergence in 1D, 2D, and 3D, as expected forlinear RK corrections. The 3D results eventually plateau around 96 points. It is likely that the difference betweenthe numeric and analytic solution is reaching the accuracy limit of the third-order quadrature in 3D, which appearsto reduce the convergence order to first-order. The inclusion of spatially-dependent cross sections does not appearto hinder convergence.When the point positions are randomly perturbed (Fig. 3c), the convergence rate stays the same, with a caveat.The algorithm for calculating kernel extents is designed for hydrodynamics and requires that the average level ofsupport is above a certain threshold, not the support for each point. It is possible that in 2D and 3D, the solution10eaches the accuracy of the poorly-supported RK calculation and stops converging for certain points. Before thisoccurs, the convergence is second-order. When the calculation is run with a higher kernel extent (not pictured), thequantitative behavior is similar but the error levels off at a lower value. The other difference between dimensions isthe integration quadrature, which is coarser and less accurate between dimensions. This could be contributing tothe leveling off of the error.The temporal convergence rate is first-order (Fig. 4a), as expected from the backward Euler time discretization.The 2D and 3D results have similar or lower error than the 1D results for a similar number of points, which may bedue to the higher level of connectivity in 2D and 3D. For this problem, the time discretization error even with thesmallest time step considered (0.001) is similar to the spatial discretization error with d points. It is expected thatto increase the accuracy of a time-dependent simulation at the point where the temporal and spatial discretizationerrors are similar, the time step would need to be decreased as the distance between points squared. The second manufactured problem represents a wave traveling from the origin outward, ψ wave = 1 + 1 t + 6 exp (cid:16) −
10 [ | x | − t ] (cid:17) , with constant opacities of σ a = 0 . and σ s = 2 . and a domain of − ≤ x α ≤ . For the steady-state case, themanufactured solution is fixed at t = 0 . .As in the sinusoidal case, the spatial convergence is second-order (Fig. 3b) and the temporal convergence is first-order (Fig. 4b). The magnitude of the error is similar in the wave and the sinusoidal case, and just as in thesinusoidal case, the error in 3D plateaus on the spatial convergence plot, probably due to the integration error. Theperturbed version of the steady-state problem (Fig. 3d) has similar behavior to the sinusoidal case described above,with second-order convergence until reaching issues with either RK kernel support or integration. One challenge in transport is handling highly absorptive regions without incurring negative fluxes. In this problem,a single ray with
Ω = { , , } is incident on a purely absorbing slab with a domain ≤ x ≤ , which is modeledin 1D. First, a constant cross section of σ a = 5 is considered with a variable number of points between 8 and 64.Then, the number of points is held constant at and the cross section is varied between and . The points areagain placed uniformly.Convergence results are shown in Fig. 5 for a few cases of the RK order and the number of neighbors, which is thenumber of other points across a kernel radius for the base connectivity used to create the overlap connectivity. Thenumber of neighbors for the reduced-radius kernels should be at least one higher than the RK order to prevent thesystem in Eq. (24) from being singular, which means that the number of neighbors for the original kernels (whichis the number reported here) should be two times the RK order plus one. For zeroth-order RK corrections, thesolution converges with approximately second-order accuracy. For first-order corrections, the solution convergeswith approximately second-order accuracy for 4 neighbors and between second and third order for 6 neighbors.With second-order corrections and 6 neighbors, the convergence order is between third and fourth. In general, fora smooth solution, the expected convergence order is one greater than the RK order.Results for the case with a constant number of points and a changing cross section are shown in Fig. 6. For thisproblem, in which the primary gradient is at the edge of the problem with the lowest point density, the MLPGapproach requires around one point per mean free path of the material to avoid negativities, which is reflected in theresults. The solution begins to exhibit negativities at σ a = 32 for the case with 6 neighbors, while for 4 neighbors,the negativities show up for σ a = 64 and above. Note that the error is an absolute error, since the normalizationwould otherwise skew the results. As in the convergence study, an increasing number of neighbors and RK orderdecrease the solution error.Based on these results, it may be tempting to use kernels with large radii and high RK order for other problems toincrease solution accuracy. One issue with this is computational cost, which increases significantly in 2D and 3Das the function radii increase. Since the RK order is limited for kernels with small radii, this also limits the RKorder. For instance, in 3D, the number of neighbors increases as r , where r is the kernel radius, so moving from 4neighbors across the kernel radius to 6 will more than triple the cost. This also increases the difficulty of solvingthe transport system (the L − operation). Another issue is negativities, which are present for all the 32-point11imulations with 6 neighbors but not for any 32-point simulations with 4 neighbors. These negativities can becomeamplified in time-dependent problems, where a negative absorption becomes an unphysical source of particles. Formore discussion on negative fluxes, see Sec. 6. This problem is described in Ref. [36] as a test of diffusion synthetic acceleration (DSA). Results for the MLPG codewith the Krylov iteration as described in Sec. 4.2 are compared to those from a code based on the discontinuousfinite element method (DFEM) with acceleration based on a variable Eddington factor (VEF), as described in Ref.[37].The geometry for the problem is shown in Fig. 7, with σ s = 200 in the wall and σ s = 0 . in the pipe. Theabsorption cross section is zero in both regions. The relatively coarse angular discretization with 16 ordinatesis identical between the MLPG and DFEM codes. The MLPG results have a spatial discretization with 716,800equally-spaced points (or 160 by 160 points per unit area), while the DFEM results are calculated on a mesh with1,335,296 elements, with a higher density of elements placed near the pipe-wall boundary. The MLPG points andDFEM mesh at quarter resolution (or 16 times fewer points) is shown in Fig. 8. The units for the problem are setsuch that the speed of light is c = 1 . The problem is run until t = 20 with a fixed time step of ∆ t = 0 . .A comparison of the crooked pipe results at t = 10 and t = 20 is shown in Fig. 9. The propagation speed of theradiation appears is nearly identical between the two codes. For the t = 10 plot, the radiation would have traveled10 unit distance at most in the 10 unit time (since c = 1 ). The minimum path the radiation could take to reach theplane at x = 2 from the source at x = − . is . unit distance. Depending on the direction, the actual distance theradiation would need to travel to reach the plane x = 2 would be at between 7 and 12 unit distance. The strongestvisible ray is at Ω α = 1 / √ , which would have reached x = 2 at t = 9 . if scattering and corners were neglected.As the radiation appears to have just reached x = 2 at t = 10 , the calculated time of arrival is close to the distancethe radiation would have traveled in that time.The results show significant ray effects, but because the two codes use the same angular quadrature, the effectsappear to be the same. Before reaching the crooked part of the problem, there are no significant differences visiblebetween the two solutions. After the radiation has gone around the obstacle, the MLPG solution is higher inmagnitude, which is visible in the t = 10 plot near the right edge of the obstacle or in the t = 20 plot at the exitingsurface of the pipe. Part of this could be due to the higher resolution along the pipe-wall interface in the DFEMsimulation, which could affect the scattering rate at the interface. The contours of the solution near the interfacesalso line up very closely, except at the wall edge at x = 0 . , where the MLPG solution reaches a further throughthe wall, which could again be due to the lower resolution near the interface.As mentioned before, this problem has been used as a test of DSA. The Krylov solution procedure [for the solutionof the MLPG system in Eq. (58)] works for this case without preconditioning, with an average of 64 iterations toconverge. As the ILUT factorization of the matrices representing the L − operation is performed once and stored,this doesn’t increase the total simulation time by nearly 64 times more than a single iteration, as the factorizationis a far larger cost than a single solve. Within each scattering iteration, the transport GMRES solver converges to atolerance of − in 55 iterations on average. The DFEM solution, however, required only one transport solve andtwo VEF solves per time step, which if applied to the MLPG solution could open up more cost-effective methodsfor solving for the scattering source.For general reference, the RK transport code takes 15,943 seconds, or 79 seconds per time step, to run the crookedpipe problem with 16 ordinates and 716,800 points on 288 processors, which equates to 2,488 points or 39,808unknowns per processor on average. This includes the time for integration, computation of the ILUT preconditionersfor all 16 directions, and convergence of the solution and scattering source. With an effective preconditioner andpossibly avoiding ILUT decompositions, the solve time would decrease significantly. The need for appropriatepreconditioning is discussed further in Sec. 6. One motivation for combining radiation transport with a smoothed particle hydrodynamics code is for a planetarydefense application: the deflection of an asteroid due to radiation from a standoff nuclear burst. In this scenario,the absorbed radiation energy ablates the surface of the asteroid, causing material to blow off and alter the orbitof the asteroid via momentum conservation. This problem is a simplified version of that scenario, a spherical rock“asteroid” that absorbs radiation from a distant point source in 2D. This problem is similar to the purely-absorbing12ersion of Kobayashi benchmarks [38], which also have features that are difficult to angularly resolve and a smallsource emitting particles into a void. The asteroid has an absorption cross section of σ a = 10 . , while the mediumsurrounding the asteroid has an absorption cross section of σ a = 0 . . Neither material includes scattering. Theasteroid has a radius of 35 and is centered at the origin. The point source is located a distance of 70 away from thesurface of the asteroid. The asteroid can be modeled by a shell, since almost all of the radiation is absorbed at thesurface of the asteroid.The shell of the asteroid is set to be 20 mean free paths thick. The distance between points is set to be 0.2 at theinside of the shell of the asteroid, 0.1 at the outside of the shell, 1.0 halfway between the source and the asteroid, and0.1 near the source. The initial angular quadrature of 4,096 ordinates is refined as described in Sec. 4.3 down to 550ordinates, of which 480 hit the asteroid. The problem is run with a single time step large enough for the radiationto propagate throughout the domain. Afterward, the numeric solution is compared to the analytic solution, whichis derived in App. B. The solution points and integration mesh for this problem are shown in Fig. 10. Note thatthe integration mesh is further broken down into triangles for use with standard quadratures (Sec. 3).The analytic and numeric solutions to the asteroid problem are shown in Fig. 11. The most obvious difference atfirst glance is the large areas at the top of the numeric solution where the solution is close to zero. These are areaswhere, by design, the angular refinement has not put a sufficient number of angles to resolve the solution. Theseshould not affect the solution at the asteroid, since there, the solution should be sufficiently resolved in the angulardomain. While there are 480 rays that hit the asteroid from the point source, the radiation is still not angularlyuniform in the region that is resolved, as can be seen by the more intense ray that hits the asteroid around thepoint { , } .Near the y = 0 plane, where there is approximately one point per mean free path in the direction the solution ischanging, the contours of the analytic and numeric solutions line up very well, with the exception of the aforemen-tioned oscillations. This agrees with the purely absorbing results (Sec. 5.2), in which the solutions with around onepoint per mean free path showed higher accuracy and few oscillations compared to those with more than one pointper mean free path.There are two connected difficulties in this problem, which are ray effects and negativities. Oscillations can beseen toward the inner surface of the asteroid at all positions, but the oscillations are by far the worst where theradiation is traveling nearly parallel to the surface of the asteroid. At these points, the solution will change from thevacuum solution just outside of the surface to nearly zero inside of the surface. This causes oscillations that lead tonegativities. The SUPG stabilization does a good job of handling the oscillations that may develop in the directionof radiation propagation (near the y = 0 plane), but more consideration is needed to prevent the oscillations thatdevelop perpendicular to the radiation propagation direction or when the solution changes discontinuously.This problem is designed to stress the code and show opportunities for future work. If the results needed tobe accurate, the source could be analytically calculated just before it hits the asteroid and inserted as a boundarysource there, which would reduce much of the need for a refined and specialized quadrature. The negatives, however,would persist, which is something that would need to be addressed before this calculation would work well in atime-dependent or thermal radiative transfer scenario, as discussed in Sec. 6. The MLPG discretization in this paper simplifies the process running a problem with meshless transport. Thefully-implicit time differencing is stable for large time steps. The SUPG stabilization works to prevent oscillationsand increase the efficiency of inverting the transport matrix. The integration with a Voronoi diagram is robust andfollows the resolution of the meshless points without user input. The SUPG stabilization and Voronoi integrationadd complexity to the code, which could be a barrier to entry, but reduces the need for specialized solvers orrepeated adjustments of a background mesh.The RK functions used in the discretization permit higher-than-second-order convergence, but practically, the radiiof the kernels should often be minimized to reduce negativities and computation cost, which constrains the RKorder. With the fully-implicit time discretization and first-order RK corrections, the results are consistent withsecond-order convergence in space and first-order convergence in time. For a purely absorbing problem with anincoming source, higher-order convergence is achieved in space with second-order RK corrections and larger kernelradii. The two manufactured solutions show the capability to represent spatially-dependent cross sections andconverge in 1D, 2D, and 3D.There are at least three issues that remain to be resolved. The first is negativities. In the purely-absorbing13roblem, around one point is needed per mean free path to avoid negativities. In the asteroid problem, negativitiesare difficult to avoid due to ray effects, as the SUPG stabilization applies numerical diffusion only in the directionof radiation propagation. A negative flux fixup method such as the zero-and-rescale approach [39] may work formeshless transport, but care would need to be taken to rescale a quantity that should always be positive and not theexpansion coefficients. While this could inhibit the effects of oscillations on time-dependent problems and perhapskeep them from growing, it would be much more difficult to remove oscillations entirely.The second issue is preconditioning. In the crooked pipe problem, the solution converged when using only GMRESto converge the scattering source and agreed well with a DFEM solution, but this required many iterations toachieve. Combining the Krylov solve with a method such as DSA [40] could significantly reduce the number ofiterations needed to converge the scattering source. The MLPG transport equation without SUPG should havethe diffusion limit, similar to a high-order DFEM discretization [41], but it is not apparent whether the same istrue with SUPG. The process of deriving DSA for the SUPG system should give information on whether it has thediffusion limit and if not, what changes may be made to the discretized system to ensure it has the diffusion limit.The third issue is ensuring proper support for the RK kernels. As shown in the manufactured problems, randomlyperturbing the point positions can lead to a limit on convergence. Ensuring that every point has the needed supportindividually through a more robust calculation may resolve these issues.The current meshless discretization works well for problems in which the solution does not go negative. Once anegative flux fixup treatment is applied, the addition of additional physics to the transport discretization such asthermal radiative transfer and radiation hydrodynamics would be more achievable. In a radiation hydrodynamicssimulation, where the meshless topology is constantly changing, the consistently discretized transport with a Voronoiintegration approach eliminates mapping to and from a mesh for radiation transport and is far cheaper than placinga non partition-of-unity quadrature for each kernel.
Acknowledgements
This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore NationalLaboratory under Contract DE-AC52-07NA27344. This document was prepared as an account of work sponsoredby an agency of the United States government. Neither the United States government nor Lawrence LivermoreNational Security, LLC, nor any of their employees makes any warranty, expressed or implied, or assumes any legalliability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, orprocess disclosed, or represents that its use would not infringe privately owned rights. Reference herein to anyspecific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does notnecessarily constitute or imply its endorsement, recommendation, or favoring by the United States government orLawrence Livermore National Security, LLC. The views and opinions of authors expressed herein do not necessarilystate or reflect those of the United States government or Lawrence Livermore National Security, LLC, and shall notbe used for advertising or product endorsement purposes.
References [1] Ted Belytschko, Yury Krongauz, Daniel Organ, Mark Fleming, and Petr Krysl. Meshless methods: an overviewand recent developments.
Computer methods in applied mechanics and engineering , 139(1):3–47, 1996.[2] John I Castor.
Radiation hydrodynamics . Cambridge University Press, 2004.[3] Stuart C Whitehouse and Matthew R Bate. Smoothed particle hydrodynamics with radiative transfer in theflux-limited diffusion approximation.
Monthly Notices of the Royal Astronomical Society , 353(4):1078–1094,2004.[4] Stuart C Whitehouse, Matthew R Bate, and Joe J Monaghan. A faster algorithm for smoothed particlehydrodynamics with radiative transfer in the flux-limited diffusion approximation.
Monthly Notices of theRoyal Astronomical Society , 364(4):1367–1377, 2005.[5] Serge Viau, Pierre Bastien, and Seung-Hoon Cha. An implicit method for radiative transfer with the diffusionapproximation in smooth particle hydrodynamics.
The Astrophysical Journal , 639(1):559, 2006.146] Lucio Mayer, Graeme Lufkin, Thomas Quinn, and James Wadsley. Fragmentation of gravitationally unstablegaseous protoplanetary disks with radiative transfer.
The Astrophysical Journal Letters , 661(1):L77, 2007.[7] Margarita Petkova and Volker Springel. An implementation of radiative transfer in the cosmological simulationcode gadget.
Monthly Notices of the Royal Astronomical Society , 396(3):1383–1403, 2009.[8] Joe J Monaghan. Smoothed particle hydrodynamics.
Reports on progress in physics , 68(8):1703, 2005.[9] Wing Kam Liu, Sukky Jun, and Yi Fei Zhang. Reproducing kernel particle methods.
International journal fornumerical methods in fluids , 20(8-9):1081–1106, 1995.[10] Nicholas Frontiere, Cody D Raskin, and J Michael Owen. CRKSPH–a conservative reproducing kernelsmoothed particle hydrodynamics scheme.
Journal of Computational Physics , 332:160–209, 2017.[11] Hamou Sadat. On the use of a meshless method for solving radiative transfer with the discrete ordinatesformulations.
Journal of Quantitative Spectroscopy and Radiative Transfer , 101(2):263–268, 2006.[12] Hamou Sadat, Cheng-An Wang, and Vital Le Dez. Meshless method for solving coupled radiative andconductive heat transfer in complex multi-dimensional geometries.
Applied Mathematics and Computation ,218(20):10211–10225, 2012.[13] Manuel Kindelan, Francisco Bernal, Pedro González-Rodríguez, and Miguel Moscoso. Application of therbf meshless method to the solution of the radiative transport equation.
Journal of Computational Physics ,229(5):1897–1908, 2010.[14] LH Liu and JY Tan. Least-squares collocation meshless approach for radiative heat transfer in absorbing andscattering media.
Journal of Quantitative Spectroscopy and Radiative Transfer , 103(3):545–557, 2007.[15] JM Zhao, JY Tan, and LH Liu. A second order radiative transfer equation and its solution by meshless methodwith application to strongly inhomogeneous media.
Journal of Computational Physics , 232(1):431–455, 2013.[16] S Kashi, A Minuchehr, A Zolfaghari, and B Rokrok. Mesh-free method for numerical solution of the multi-groupdiscrete ordinate neutron transport equation.
Annals of Nuclear Energy , 106:51–63, 2017.[17] LH Liu and JY Tan. Meshless local Petrov-Galerkin approach for coupled radiative and conductive heattransfer.
International journal of thermal sciences , 46(7):672–681, 2007.[18] Brody Bassett and Brian Kiedrowski. Meshless local Petrov–Galerkin solution of the neutron transport equationwith streamline-upwind petrov–galerkin stabilization.
Journal of Computational Physics , 377:1–59, 2019.[19] JE Morel and JM McGhee. A self-adjoint angular flux equation.
Nuclear Science and Engineering , 132(3):312–325, 1999.[20] Brody R Bassett, J Michael Owen, and Thomas A Brunner. Efficient smoothed particle radiation hydrody-namics i: Thermal radiative transfer. arXiv preprint arXiv:2001.11606 , 2020.[21] Brody R Bassett, J Michael Owen, and Thomas A Brunner. Efficient smoothed particle radiation hydrody-namics ii: Radiation hydrodynamics. arXiv preprint arXiv:2001.11608 , 2020.[22] Shaofan Li and Wing Kam Liu. Meshfree and particle methods and their applications.
Appl. Mech. Rev. ,55(1):1–34, 2002.[23] Vinh Phu Nguyen, Timon Rabczuk, Stéphane Bordas, and Marc Duflot. Meshless methods: a review andcomputer implementation aspects.
Mathematics and computers in simulation , 79(3):763–813, 2008.[24] Satya N Atluri and Tulong Zhu. A new meshless local petrov-galerkin (mlpg) approach in computationalmechanics.
Computational mechanics , 22(2):117–127, 1998.[25] Donat Racz and Tinh Quoc Bui. Novel adaptive meshfree integration techniques in meshless methods.
Inter-national journal for numerical methods in engineering , 90(11):1414–1434, 2012.[26] Amir Khosravifard and Mohammad Rahim Hematiyan. A new method for meshless integration in 2d and 3dgalerkin meshfree methods.
Engineering Analysis with Boundary Elements , 34(1):30–40, 2010.1527] Jiun-Shyan Chen, Cheng-Tang Wu, Sangpil Yoon, and Yang You. A stabilized conforming nodal integrationfor galerkin mesh-free methods.
International journal for numerical methods in engineering , 50(2):435–466,2001.[28] JX Zhou, JB Wen, HY Zhang, and L Zhang. A nodal integration and post-processing technique based onvoronoi diagram for galerkin meshless methods.
Computer methods in applied mechanics and engineering ,192(35-36):3831–3843, 2003.[29] MA Puso, JS Chen, E Zywicz, and W Elmer. Meshfree and finite element nodal integration methods.
Inter-national Journal for Numerical Methods in Engineering , 74(3):416–446, 2008.[30] J Michael Owen, Jens V Villumsen, Paul R Shapiro, and Hugo Martel. Adaptive smoothed particle hydrody-namics: Methodology. ii.
The Astrophysical Journal Supplement Series , 116(2):155, 1998.[31] J Michael Owen. Polyclipper. https://github.com/LLNL/PolyClipper , 2020.[32] Freddie D Witherden and Peter E Vincent. On the identification of symmetric quadrature rules for finiteelement methods.
Computers & Mathematics with Applications , 69(10):1232–1241, 2015.[33] Satya N Atluri, H-G Kim, and J Ya Cho. A critical assessment of the truly meshless local petrov-galerkin(mlpg), and local boundary integral equation (lbie) methods.
Computational mechanics , 24(5):348–372, 1999.[34] Michael A Heroux, Roscoe A Bartlett, Vicki E Howle, Robert J Hoekstra, Jonathan J Hu, Tamara G Kolda,Richard B Lehoucq, Kevin R Long, Roger P Pawlowski, Eric T Phipps, et al. An overview of the trilinosproject.
ACM Transactions on Mathematical Software (TOMS) , 31(3):397–423, 2005.[35] Joshua J Jarrell and Marvin L Adams. Discrete-ordinates quadrature sets based on linear discontinuous finiteelements. In
International Conference on Mathematics and Computational Methods Applied to Nuclear Scienceand Engineering (M&C 2011), Rio de Janeiro, RJ, Brazil , 2011.[36] RP Smedley-Stevenson, AW Hagues, and J Kphzi. A benchmark for assessing the effectiveness of diffusionsynthetic acceleration schemes. In
Proceedings of the Joint International Conference on Mathematics andComputation, Supercomputing in Nuclear Applications and the Monte Carlo Method (M&C 2015), Nashville,TN , 2015.[37] Ben C Yee, Samuel S Olivier, Terry S Haut, Milan Holec, Vladimir Z Tomov, and Peter G Maginot. Aquadratic programming flux correction method for high-order dg discretizations of sn transport.
Journal ofComputational Physics , page 109696, 2020.[38] Keisuke Kobayashi, Naoki Sugimura, and Yasunobu Nagaya. 3D radiation transport benchmark problems andresults for simple geometries with void region.
Progress in Nuclear Energy , 39(2):119–144, 2001.[39] Steven Hamilton, Michele Benzi, and Jim Warsa. Negative flux fixups in discontinuous finite element sntransport. In
International Conference on Mathematics, Computational Methods and Reactor Physics (M&C2009), American Nuclear Society, LaGrange Park, Illinois, USA . Citeseer, 2009.[40] James S Warsa, Todd A Wareing, and Jim E Morel. Krylov iterative methods and the degraded effectiveness ofdiffusion synthetic acceleration for multidimensional sn calculations in problems with material discontinuities.
Nuclear science and engineering , 147(3):218–248, 2004.[41] TERRY S Haut, BS Southworth, PETER G Maginot, and VLADIMIR Z Tomov. Dsa preconditioning for dgdiscretizations of s _{ N } transport and high-order curved meshes. arXiv preprint arXiv:1810.11082 , 2018. A Integral transformations
This section describes transformations of volume and surface integrals from a reference element (a line segment in1D, a triangle in 2D, and a tetrahedron in 3D) to an element in physical space. For information on the meshlessintegration methods that produce these elements, see Sec. 3.16 .1 Volume integrals
The integrals described in Sec. 3 need to be mapped from reference space R to physical space V , (cid:90) V f ( x ) dV = (cid:90) R f ( x ( ξ )) | J | dR, (64)where x and ξ are the coordinates in physical and reference space, respectively, and J is the Jacobian determinantof the transformation, J αβ = ∂x β ∂ξ α . (65)In discrete form, the quadrature is mapped using (cid:90) V f ( x ) dV ≈ (cid:88) m w m | J m | f ( x ( ξ m )) , (66)where w m are the weights of the reference quadrature, which effectively converts the ordinates to x ( ξ m ) and theweights to w m | J m | for the integration.The line quadrature is assumed to have the bounds − ≤ ξ ≤ , while the triangular and tetrahedral quadratureshave the bounds ≤ α ξ α ≤ , where α represents a vector of ones. For a line with points p and p in physicalspace, the mapping to reference space is x ( ξ ) = p + ( p − p ) ( ξ + 1)2 . (67)For a triangle with points p , p , and p or a tetrahedron with an additional point p , the mapping is x ( ξ ) = p + ( p α − p ) ξ α . (68) A.2 Surface integrals
For surface integrals, the mapping is between reference space T with coordinates ξ and physical space S withcoordinates x , (cid:90) S f ( x ) dS = (cid:90) T f ( x ( ξ )) KdT, (69)where K is a differential surface element. For a triangle mapped from a 2D reference element to a 3D surfaceelement, this term is K = | ∂ ξ x × ∂ ξ x | , (70)while for a line mapped from a 1D to a 2D line element, the differential element is K = | ∂ ξ x | , (71)with the line and surface integrals mapped as in Eqs. (67) and (68), with the exception that the p vectors haveone more element than the ξ ones (e.g. in 3D, the p vectors have three elements, while the ξ , which is a surfaceparameterization, has only two). In discrete form, the surface integrals are (cid:90) S f ( x ) dS ≈ (cid:88) m w m K m f m . (72) B Analytic solution to asteroid problem
In this section, the analytic solution for the asteroid problem described in Sec. 5.4 is derived. In 2D, the “asteroid”is actually an infinite cylinder and the “point source” is a line source. The transport equation for this problem canbe written in cylindrical geometry as µ r ∂ r ( rψ ) + σ a ψ = δ ( r ) , (73)17here r is the distance from the point source and µ is the cosine of the angle between the x-y plane and the radiationpropagation direction. The cross section is σ a = (cid:40) σ a, background , r < d asteroid ,σ a, asteroid , otherwise , (74)where d asteroid is the distance from the source to the asteroid for a given evaluation point. The solution to thisequation is ψ = 1 r exp (cid:18) − µ ω (cid:19) , (75)where ω = (cid:90) r σ a dr (cid:48) , (76)which can be calculated using only the distance travelled in each of the asteroid and the background along withtheir respective cross sections. Integrating this equation over all forward angles, ≤ µ ≤ (which is also the onlyset of angles for which the solution is nonzero), the scalar flux solution is φ = 1 r (cid:20) exp ( − ω ) − ω (cid:90) ∞ ω ω (cid:48) exp ( − ω (cid:48) ) dω (cid:48) (cid:21) . (77)The remaining integral is the exponential integral, which can be evaluated directly in many mathematical softwarepackages.To find the distance from the source to the asteroid, let the position of the source in the xy plane be x src and theevaluation point be x eval . The parametric equation for the line connecting these is x = x src + ( x eval − x src ) s, (78)where s = 0 is located at the source and s = 1 is at the evaluation point. Inserting this equation into the equationfor the asteroid centered at the origin, | x | = r asteroid , (79)and solving for s results in the two intercept locations, s int = − (cid:96) ± (cid:112) (cid:96) − (cid:96) (cid:96) (cid:96) , (80)where (cid:96) = x α src x α src − r , (81) (cid:96) = x αsrc ( x α eval − x α src ) , (82) (cid:96) = ( x α eval − x α src ) ( x α eval − x α src ) . (83)If (cid:96) − (cid:96) (cid:96) < , then the ray from the source to the evaluation point does not travel through the asteroid. Theevaluation point may be before the ray intersects with the asteroid, inside the asteroid, or after the ray has exitedthe asteroid. Given the distances traveled in each of the asteroid and the background, the integral in Eq. (76) canbe evaluated, which permits evaluation of either the angular or scalar flux [Eq. (75) and (77), respectively].18 lgorithm 1 Meshless integration algorithm shown for the example functions (cid:104) ∂ αx U i , f U j (cid:105) and ( U i , n α g α U j ) . c l a s s B i l i n e a r K e r n e l D K e r n e l : p u b l i c B i l i n e a r I n t e g r a l f u n c a d d T o I n t e g r a l ( b a s i s , d b a s i s , o r d i n a t e , w e i g h t ) : s e t numBasis t o number o f b a s i s f u n c t i o n s ( n e i g h b o r s ) f o r ( i = 0 ; i < numBasis ; ++i ) f o r ( j = 0 ; j < numBasis ; ++j ) i n t e g r a l ( i , j ) += w e i g h t ∗ c o e f f i c i e n t ( o r d i n a t e ) ∗ b a s i s [ i ] ∗ b a s i s [ j ] S p a r s e M a t r i x i n t e g r a l c l a s s B i l i n e a r S u r f a c e K e r n e l K e r n e l : p u b l i c B i l i n e a r I n t e g r a l f u n c a d d T o S u r f a c e I n t e g r a l ( b a s i s , o r d i n a t e , weight , normal ) : s e t numBasis t o number o f b a s i s f u n c t i o n s ( n e i g h b o r s ) f o r ( i = 0 ; i < numBasis ; ++i ) f o r ( j = 0 ; j < numBasis ; ++j ) i n t e g r a l ( i , j ) += w e i g h t ∗ normal . dot ( c o e f f i c i e n t ( o r d i n a t e ) ) ∗ b a s i s [ i ] ∗ b a s i s [ j ] S p a r s e M a t r i x i n t e g r a l f u n c p e r f o r m I n t e g r a t i o n ( v o l u m e I n t e g r a l s , s u r f a c e I n t e g r a l s ) : f o r ( i = 0 ; i < numPoints ; ++i ) s e t n e i g h b o r s t o n e i g h b o r s o f p o i n t i s e t c e l l t o t h e v o r o n o i t e s s e l a t i o n f o r t h e p o i n t i decompose c e l l i n t o s u b c e l l s f o r ( c = 0 ; c < n u m S u b c e l l s ; ++c ) s e t o r d i n a t e s and w e i g h t s t o volume q u a d r a t u r e f o r s u b c e l l s [ c ] f o r ( q = 0 ; q < numOrdinates ; ++q ) s e t b a s i s / d b a s i s t o RK e v a l u a t i o n s / d e r i v a t i v e s a t o r d i n a t e s [ q ] f o r eachn e i g h b o r f o r each v o l u m e I n t e g r a l i n t e g r a l . a d d T o I n t e g r a l ( b a s i s , d b a s i s , o r d i n a t e s [ q ] , w e i g h t s [ q ] ) decompose c e l l s u r f a c e i n t o s u b s u r f a c e s f o r ( s = 0 ; s < n u m S u b s u r f a c e s ; ++s ) s e t o r d i n a t e s and w e i g h t s t o s u r f a c e q u a d r a t u r e f o r s u b s u r f a c e s [ s ] s e t normal t o normal f o r s u b s u r f a c e s [ s ] f o r ( q = 0 ; q < numOrdinates ; ++q ) s e t b a s i s t o RK e v a l u a t i o n s f o r each n e i g h b o r p o i n t a t o r d i n a t e s [ q ] f o r each s u r f a c e I n t e g r a l : i n t e g r a l . a d d T o S u r f a c e I n t e g r a l ( b a s i s , o r d i n a t e s [ q ] , w e i g h t s [ q ] , normal [ s ] ) Algorithm 2
Method for refinement of the angular quadrature to resolve a region of interest. f u n c g e t R e f i n e d Q u a d r a t u r e ( minRule , g o a l R u l e , a n g l e H i t s O b j e c t ) : i n i t i a l i z e o r d i n a t e s and w e i g h t s t o be empty g e t q u a d r a t u r e f o r t h e g o a l r u l e f o r ( r u l e = minRule ; r u l e < g o a l R u l e ; ++r u l e ) i n i t i a l i z e o r d i n a t e s and w e i g h t s f o r t h i s q u a d r a t u r e f o r ( i = 0 ; i < numPoints f o r t h i s r u l e ; ++i ) g e t t h e i n d i c e s from t h e g o a l q u a d r a t u r e t h a t t h i s p o i n t r e p r e s e n t s a sg o a l P o i n t s i f not a n g l e H i t s O b j e c t ( g o a l P o i n t s [ j ] ) f o r any such g o a l q u a d r a t u r e i n d e x j : add t h i s p o i n t t o t h e o r d i n a t e s and w e i g h t s remove t h i s p o i n t from t h e g o a l q u a d r a t u r e and w e i g h t s add t h e r e m a i n i n g g o a l q u a d r a t u r e o r d i n a t e s and w e i g h t s t o t h e o r d i n a t e s andw e i g h t s r e t u r n t h e o r d i n a t e s and w e i g h t s
12 16 24 32 48 64 96 128number of points10 r e l a t i v e e rr o r dim = 1dim = 2dim = 3 2 nd order (a) Manufactured sinusoidal solution, uniform. r e l a t i v e e rr o r dim = 1dim = 2dim = 3 2 nd order (b) Manufactured wave solution, uniform. r e l a t i v e e rr o r dim = 1dim = 2dim = 3 2 nd order (c) Manufactured sinusoidal solution, perturbed. r e l a t i v e e rr o r dim = 1dim = 2dim = 3 2 nd order (d) Manufactured wave solution, perturbed. Figure 3: Spatial convergence of the manufactured solutions at steady-state in 1D, 2D, and 3D, with the dottedline indicating second-order convergence. For two of the cases, the point positions are randomly perturbed. Theperturbation may cause insufficient support for the RK functions or lower integration accuracy, leading to a lack ofconvergence. r e l a t i v e e rr o r dim = 1dim = 2dim = 3 1 st order (a) Manufactured sinusoidal solution. r e l a t i v e e rr o r dim = 1dim = 2dim = 3 1 st order (b) Manufactured wave solution. Figure 4: Temporal convergence of the manufactured solutions in 1D, 2D, and 3D, with the dotted line indicatingfirst-order convergence. 21
12 16 24 32 48 64number of points10 r e l a t i v e e rr o r neighb, RK ord4, 04, 16, 06, 16, 2 conv.order2 nd rd th Figure 5: Relative error of the numeric solution to the purely absorbing problem as the number of points is increasedfor various combinations of the number of neighbors and RK correction order, with the dotted lines indicating second,third, and fourth-order convergence. a b s o l u t e e rr o r neighb, RK ord4, 04, 16, 06, 16, 2 Figure 6: Absolute error of the numeric solution to the purely absorbing problem for various combinations of thenumber of neighbors, RK correction order, and absorption cross section.22igure 7: Geometry of the cooked pipe problem. The scattering opacities are σ s = 200 in the black region and σ s = 0 . in the gray region.Figure 8: The MLPG points (top) and DFEM mesh (bottom) for the crooked pipe problem, at 1/4 resolution forvisibility. Note that the DFEM mesh has additional resolution near the thin-thick boundary.23 a) t = 10 (b) t = 20= 20