Implicit Filtered P for High-Energy Density Thermal N Radiation Transport using Discontinuous Galerkin (cid:73) Finite Elements 6 1 0 Vincent M. Labourec,∗, Ryan G.McClarrenc, Cory D.Haucka,b 2 y aComputational and Applied Mathematics Group, Oak Ridge National Laboratory, Oak a Ridge, TN 37831 USA M bDepartment of Mathematics, University of Tennessee Knoxville, TN 37996-1320 cNuclear Engineering Department, Texas A&M University College Station, TX 77843 2 1 ] h Abstract p - p In this work, we provide a fully-implicit implementation of the time-dependent, m filtered spherical harmonics (FPN) equations for non-linear, thermal radiative transfer. We investigate local filtering strategies and analyze the effect of the o c filter on the conditioning of the system, showing in particular that the filter . improves the convergence properties of the iterative solver. We also investi- s c gate numerically the rigorous error estimates derived in the linear setting, to si determine whether they hold also for the non-linear case. Finally, we simulate y a standard test problem on an unstructured mesh and make comparisons with h implicit Monte-Carlo (IMC) calculations. p [ Keywords: Radiation transport, Thermal radiative transfer, Spherical 2 harmonics, Spectral filtering, Fully implicit methods, Discontinuous Galerkin v 2 4 1. Introduction 2 8 The equations of thermal radiative transfer describe the movement of pho- 0 . tons through a material as well as the exchange of energy between the photon 1 radiation and the material. There are two equations: a radiation transport 0 equation that tracks the energy in the radiation field via an angular intensity 6 1 I and a temperature equation that tracks the internal energy of the material. : v i X (cid:73)Thismaterialisbased,inpart,uponworksupportedbytheNationalScienceFoundation r under Grant No. 1217170. The research of the third author is sponsored by the Office a of Advanced Scientific Computing Research; U.S. Department of Energy. The work was performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC underContractNo. De-AC05-00OR22725. ∗Correspondingauthor. Tel.:+19792248506 Email addresses: [email protected](VincentM.Laboure),[email protected] (RyanG.McClarren),[email protected](CoryD.Hauck) Preprint submitted to Elsevier May 16, 2016 The coupling of these two equations reflects the exchange of energy as photons are emitted and absorbed by the material. Various numerical methods for the radiation transport equation have been developedtosolvetheradiativetransferproblem. Thechallengehereisthatthe angular intensity is, in the most general setting, a function of six phase space variables (position, energy, and direction of propagation) plus time. The most common approaches are implicit Monte Carlo methods [1, 2], discrete ordinate methods [3], spectral approximations [4], finite element discretizations [5], and nonlinear moments methods [6, 7, 8]. In this paper, we focus on a variation of the spherical harmonics (or P ) N method. The P method is a spectral Galerkin method that approximates N the angular dependence of the radiation intensity using a finite expansion in spherical harmonics up to degree N. The result is a linear, hyperbolic system of time-dependent equations for the expansion coefficients, which can then be discretized with respect to space and time in a variety of ways. The P approach offers several benefits. Among these are spectral con- N vergence for smooth solutions and preservation of the rotational invariance of the transport equation.1 However, the method also poses challenges. Chief among these is that the P approximation of the angular intensity can be N highly-oscillatory and even negative when the underlying exact solution is not sufficiently smooth; this happens typically in regions where the material cross- section is small. In addition to being non-physical, the negative radiation en- ergy can cause the material temperature T to become negative, in which case the model for photon emission is not well-defined.2 In addition, the material cross-section may become negative and thereby introduce instabilities into the simulation. To address the problem of Gibbs phenomena in the P approximation, N McClarren and Hauck [9] applied filtering techniques to smooth out the an- gular dependency of the solution; they coined the name filtered P (FP ) for N N this method. While filtering does not ensure positivity of the approximation for I, it does suppress oscillations in the P approximation at a computa- N tional cost that is much lower than other closures that are robustly positive, such as positive P (PP ) closures [10, 11, 12] and entropy-based closures N N [7, 8]. In practice, the filtering approach has so far shown very promising re- sults[9,13,14,15,16]. However,themethodhasyettobeimplementedwithan implicit time-integration scheme, which is often preferred due to the extremely fast scales in the transport equation. Indeed, in an explicit scheme, the time step ∆t required by particle advection is bounded by ∆x/c, where ∆x charac- terizesthesizeofthespatialmeshandcisthespeedoflight.3 Suchacondition 1Roughlyspeaking,thesolutionoftheequationisunchangedwhenthespatialandangular variablesinphasespaceundergothesamerotation. 2ForT ≤0,theexpressionofthePlanckian–asshowninEq.4below–isnotintegrable withrespecttoenergy. 3Anexplicittreatmentoftheenergyexchangetermsmayrequireanevensmallertimestep. However,becausethesetermsarespatiallylocal,theyarerelativelyeasytotreatimplicitly. 2 is often too restrictive. Implicit methods, on the other hand, maintain stability with a much larger time step. However, each step requires the inversion of a large set of algebraic equations. In this paper, we detail an implicit implementation of the filtered P equa- N tionsusingDiscontinuousGalerkin(DG)FiniteElements. TheDGapproachis one of several possible spatial discretization methods. Other finite element ap- proaches for P include least-squares formulations [17, 18, 19], parity-based N formulations [20, 21], self-adjoint formulations [22], and streamlined-upwind Petrov-Galerkin methods [23, 24]. Discontinuous Galerkin methods were in- vented for transport problems in Ref.[25]. There it was observed that the discontinuous basis, while more expensive than a standard continuous approxi- mation, give better approximations to problems with non-smooth solutions. In additiontobeingrobustinstreamingregimes,wherenon-smoothsolutionstyp- icallyoccur,DGmethods(withasufficientlyrichbasisset)alsoperformwellin the diffusion limit [26, 27, 28].4 A semi-implicit discretization of the P equa- N tions with DG methods, which treats the flux terms explicitly, can be found in Ref.[31]. The remainder of this paper is organized as follows. The radiative transfer equations and the FP equations are presented in Section 2. The spatial dis- N cretizationoftheFP equationsispresentedinSection3. InSection4,weshow N the impact the filter has on the convergence properties of the iterative solver and then consider the error estimates derived in Ref.[32] for the linear setting. Finally, in Section 5, we test the method with different filtering strategies on the challenging benchmark problem known as the Crooked Pipe [33]. Because this problem is particularly hard to converge, we first show good agreement betweenourcodeandimplicitMonte-Carlocalculationsonasimplifiedversion. Wethenshowfortheharderproblemthatthefiltermitigatesdeficienciesinthe P solutions, especially for smaller values of N. N 2. Implicit Filtered P N We consider the grey (frequency integrated) form of the thermal radiative transfer equations, given by [34]: 1∂I σ +Ω(cid:126) ·∇(cid:126)I+σ (T)I =σ (T)B(T)+ sφ+Q, (1) c ∂t t a 4π ∂ (cid:0) (cid:1) E(T)=σ (T) φ−4πB(T) , (2) ∂t a alongwithappropriateinitialandboundaryconditions. Eq.1governstheangu- larintensityI((cid:126)r,Ω(cid:126),t)ofthephotonradiation,with(cid:126)r andΩ(cid:126) being, respectively, the spatial and angular coordinates and t being the time. Meanwhile, Eq.2 4Roughly speaking, this limit occurs when particle interactions with the surrounding medium isotropize the radiation field and the angular average of the photon distribution satisfiesamuchsimplerdiffusionequation[29,30]. 3 governstheevolutionofthematerialenergyE(T), whereT((cid:126)r,t)isthematerial temperature. The derivative C =E(cid:48)(T) is the material heat capacity; for cal- v culations, we assume it is independent of T, although the formulations do not require it. The constant c is the speed of light; σ , σ , and σ =σ +σ are the s a t s a scattering, absorption, and total macroscopic cross-sections, respectively, with unitsofinverselength;Qisthe(known)volumetricsource. Thescalarintensity φ((cid:126)r,t) is the integral of the specific intensity with respect to angle (cid:90) φ≡ IdΩ, (3) S2 where S2 is the unit sphere; the frequency-integrated Planckian blackbody source is given by (cid:90) ∞ 2hν3 1 acT4 B(T)≡ dν = , (4) c2 exp(hν)−1 4π 0 kT where a = (8π5k4)/(15h3c3) is the radiation constant, with h and k being the Planck and Boltzmann constants, respectively. This integral is only defined for T >0, which is one reason to maintain a positive material temperature. We assume that Eq.1 is defined over a bounded spatial domain D and we let S = D × S2. Boundary conditions for I must be specified for incoming data—that is, on the set ∂S− ={((cid:126)r,Ω(cid:126))∈∂D×S2 :(cid:126)n ((cid:126)r)·Ω(cid:126) <0}, (5) 0 where (cid:126)n ((cid:126)r) is the outward normal at a point(cid:126)r ∈∂D. 0 2.1. Fully-implicit radiation transfer ApplyingthebackwardEulermethodtodiscretizeEqs.1and2intimeleads to the following quasi-steady form of the radiative transfer system: σ Ω(cid:126) ·∇(cid:126)In+1+σ∗(Tn+1)In+1 =σ (Tn+1)B(Tn+1)+ sφn+1+Q∗, (6) t a 4π E(Tn+1)−E(Tn) (cid:16) (cid:17) =σ (Tn+1) φn+1−4πB(Tn+1) , (7) ∆t a where 1 1 (cid:90) tn+1 In σ∗ =σ + and Q∗ = Qdt+ . (8) t t c∆t ∆t c∆t tn Here and throughout, the superscript n indicates the discrete approximation of a time-dependent quantity at time tn. When a superscript is not specified, it is assumed that such approximations are evaluated at tn+1. The system 6-7 is nonlinear due to the Planckian term B and possibly the materialproperties. Inthiswork,wechooseafullynonlineartreatmentalthough ithasbeenshownthatexpandingB aboutTn andevaluatingthecross-sections at the previous time step also performs well [35, 31]. 4 2.2. Spherical Harmonics Expansion of the transport equation In the P equations, I is approximated by a finite spherical harmonic ex- N pansion:5 N (cid:96) In+1((cid:126)r,Ω(cid:126))≈Iˆ((cid:126)r,Ω(cid:126))=(cid:88) (cid:88) Im((cid:126)r)Rm(Ω(cid:126)), (9) (cid:96) (cid:96) (cid:96)=0m=−(cid:96) where,forvariablesµ∈[−1,1]andφ∈[0,2π)suchthatΩ(cid:126) =(cid:112)1−µ2cosϕ(cid:126)e + x (cid:112) 1−µ2sinϕ(cid:126)e +µ(cid:126)e , the real-form spherical harmonics are given by: y z √  2CmPm(µ)cos(mϕ), 0<m≤(cid:96)≤N  (cid:96) (cid:96) Rm(Ω(cid:126))= C0P0(µ), 0≤(cid:96)≤N . (10) (cid:96) √(cid:96) (cid:96)  2C|m|P|m|(µ)sin(|m|ϕ), 0<−m≤(cid:96)≤N (cid:96) (cid:96) (cid:113) HereCm = (2(cid:96)+1)((cid:96)−m)! isanormalizationconstantchosensuchthat(cid:82) RmRm(cid:48)dΩ= (cid:96) 4π ((cid:96)+m)! S2 (cid:96) (cid:96)(cid:48) δ δ , with δ being the Kronecker delta, and Pm denotes the associated (cid:96),(cid:96)(cid:48) m,m(cid:48) (cid:96),(cid:96)(cid:48) (cid:96) Legendre polynomial of degree (cid:96) and order m. Integrating Eq.6 in angle against Rm and applying the approximation in (cid:96) Eq.9 gives, for all ((cid:96),m)∈N ≡{((cid:96),m)∈N2 :0≤|m|≤(cid:96)≤N}, (cid:90) √ Ω(cid:126) ·∇(cid:126)Iˆ RmdΩ+σ∗Im−σ I0δ = 4πσ Bδ +Qm, (11) (cid:96) t (cid:96) s 0 (cid:96),0 a (cid:96),0 (cid:96) S2 whereQm =(cid:82) Q∗RmdΩ.InEq.11,theangularmomentsarecoupledtoeach (cid:96) S2 (cid:96) other only through the streaming operator Ω(cid:126) ·∇(cid:126). This coupling is expressed through the matrices (cid:90) A ≡ Ω(cid:126) ·(cid:126)e RRT dΩ , χ∈{x,y,z}, (12) χ χ S2 where R is the vector containing the spherical harmonics Rm. These matrices l can be evaluated using well-known recursion relations (see [32, 16] or [36] for the complex version) or exact quadrature rules. We collect the expansion coefficients Im into a vector I using a consistent (cid:96) ordering with a single index, and write Eq.11 as the following linear system: A ∂I +A ∂I +A ∂I +σ∗I−σ Φ=(cid:0)√4πσ B(cid:1)1+Q, (13) x∂x y∂y z∂z t s a where 1 = (1 0 ··· 0)T and Φ = (I0 0 ··· 0)T. In a slight abuse of 0 notation, we denote the jth component of I in the single index ordering by I = Im, where the map ((cid:96),m) ↔ j is a bijection between the two sets of j (cid:96) indices.6 The same convention will be used for Qm. (cid:96) 5EventhoughIˆdependsonbothnandN,weomitthesedependenciesinordertosimplify thenotation. 6Wehavethusassumedthat((cid:96),m)=(0,0)isassociatedtoj=1. 5 The solution vector I has (N +1)2 components, but in reduced geometries, thereareonlyP <(N+1)2thatarenotredundantortriviallyzero. IfI depends on only two spatial dimensions, then P = 1(N +1)(N +2); if I depends on 2 only one spatial dimension, then P =N +1 [36].7 2.3. Angular filtering Thepurposeofangularfilteringistoreduceunphysicaloscillationsthatcan arise from truncating the spherical harmonics expansion. It has been demon- strated theoretically in [37] and observed numerically in [11, 38, 39] that these oscillations may lead to negative solution values for the scalar intensity φ. In itsoriginalimplementation, thefiltersuppressesthembydampinghigher-order angular moments ((cid:96)>0) after each time step in the given temporal integration scheme and, in doing so, effectively mitigates negative scalar intensities in the P solution. Inaddition,thefilterwasconstructedinsuchawayastoconserve N energy, preserve rotational invariance, and maintain formal convergence of the solution as N goes to infinity. Radice et al.[16] later showed that with an appropriate modification of the filter strength, one can derive a modified set of equations. In effect, their for- mulation adds artificial scattering to the system, replacing Eq.13 by: A ∂I +A ∂I +A ∂I +σ∗I+σDI−σ Φ=(cid:0)√4πσ B(cid:1)1+Q, (14) x∂x y∂y z∂z t f s a where σ is a free parameter, (DI)m = f((cid:96),N)Im, and the filter function f is f (cid:96) (cid:96) given by: (cid:18) (cid:19) (cid:96) f((cid:96),N)≡−logρ . (15) filterType N +1 The Lanczos and spherical spline filters are considered [16]: sinζ 1 ρ (ζ)= ; ρ (ζ)= . (16) Lanczos ζ SSpline 1+ζ4 The variable σ in Eq.14 is a tuning parameter – henceforth called filter f strength – that may be spatially dependent. Strategies for determining a good localvalueofσ arediscussedinSection4.1. Inthiscontext,oneofthestrengths f of the reformulation in [16] is that — unlike the original implementation in [9] — the filter strength is independent of the size of the time step and the spatial mesh [16]. Thus the value of σ needs to be tuned only once, and this can be f done using relatively cheap simulations on coarse meshes. 3. Spatial Discretization We discretize Eq.14 in space, along with Eq.7 for the material tempera- ture,usingtheDiscontinuousGalerkinFiniteElementMethod(DGFEM).This 7Inpractice,wesolveonlyforthenontrivialmoments,butforsimplicitywemaintainthe notationsinEq.13evenwhenP <(N+1)2. 6 method is, by now, fairly standard. Thus the presentation here will be brief. Roughly speaking, the method relies on a piecewise polynomial approximation of the true solution. It requires the specification of a numerical flux at points of discontinuities, the effect of which is to add a stabilizing term to the usual variational form. For details, we refer the reader to the review in [40]. 3.1. Variational formulation Let T be a collection of open convex, polyhedral cells K ⊂ D such that ∪K = D, and let h > 0 be the size of the largest disk that can be inscribed inside any cell K. Let Γ be the set of interior facets: int Γ ={e:e=K ∩K for any K ,K ∈T, K (cid:54)=K }. (17) int 1 2 1 2 1 2 Let V be a finite-dimensional trial space of functions that are polynomial on each K ∈ T. For each j, we seek a function Ih ∈ V that approximates the j coefficient I =Im. Thus the approximation Ih of the vector-valued function I j l lives in the (Cartesian) product space VP. We also approximate T by Th ∈V. The formulation of the DGFEM is as follows: Find (Ih,Th) ∈ VP × V such that a((Ih,Th),v) = L(v) for all v ∈ VP+1, where for each ((u,θ),v) ∈ (VP ×V)×VP+1, P+1 P+1 (cid:88) (cid:88) a((u,θ),v)= a ((u,θ),v ) and L(v)= L (v ). (18) i i i i i=1 i=1 Here a and L denote (for each i, 1 ≤ i ≤ P) forms associated to the i-th i i equation of Eq.14, and to Eq.7 for i=P +1. They are derived by multiplying the corresponding equation by v and integrating over each cell K ∈ T. For i 1 ≤ i ≤ P, following an integration by parts, a can expressed as the sum of i four terms: a ((u,θ),v )=avol((u,θ),v )+aint(u,v )+aext(u,v )+aBC(u,v ). (19) i i i i i i i i i i Here avol is the volumetric contribution; aint is the contribution from interior i i facets; aext is a contribution from exterior facets (along the boundary of the i spatial domain ∂D); and aBC is a boundary contribution. i The term aBC accounts for any boundary conditions that express incoming i informationintermsofoutgoinginformation(suchasreflectiveboundaries)and will be discussed in Section 3.3. The remaining terms are8 avol((u,θ),v )=(cid:18)(cid:88)P (cid:90) (cid:0)(σ∗δ +σD )v −A ∂vi −A ∂vi −A ∂vi(cid:1)u dx(cid:19) i i t ij f ij i x,ij ∂x y,ij ∂y z,ij ∂z j j=1 D (cid:90) (cid:16) √ (cid:17) − σ u + 4πσ B v δ dx, s 1 a 1 i,1 D (20) 8Recall that we have assumed that v1 is the test function associated to the 0-th moment equationofEq.14andthatvP+1 isthetestfunctionassociatedtothetemperatureequation (Eq.7). 7 (cid:88)P (cid:88) (cid:18)(cid:90) (cid:16) 1 (cid:17) aint(u,v )= (cid:126)e ·(cid:126)nA v (cid:104)u (cid:105)+ |(cid:126)e ·(cid:126)n|M v u ds i i j=1e∈Γint e x x,ij(cid:74) i(cid:75) j 2 x x,ij(cid:74) i(cid:75)(cid:74) j(cid:75) (cid:90) (cid:16) 1 (cid:17) + (cid:126)e ·(cid:126)nA v (cid:104)u (cid:105)+ |(cid:126)e ·(cid:126)n|M v u ds e y y,ij(cid:74) i(cid:75) j 2 y y,ij(cid:74) i(cid:75)(cid:74) j(cid:75) (cid:90) (cid:16) 1 (cid:17) (cid:19) + (cid:126)e ·(cid:126)nA v (cid:104)u (cid:105)+ |(cid:126)e ·(cid:126)n|M v u ds , e z z,ij(cid:74) i(cid:75) j 2 z z,ij(cid:74) i(cid:75)(cid:74) j(cid:75) (21) (cid:88)P (cid:18)(cid:90) 1(cid:16) (cid:17) aext(u,v )= (cid:126)e ·(cid:126)n A +|(cid:126)e ·(cid:126)n |M v u ds i i 2 x 0 x,ij x 0 x,ij i j j=1 ∂D (cid:90) 1(cid:16) (cid:17) (22) + (cid:126)e ·(cid:126)n A +|(cid:126)e ·(cid:126)n |M v u ds 2 y 0 y,ij y 0 y,ij i j ∂D (cid:90) 1(cid:16) (cid:17) (cid:19) + (cid:126)e ·(cid:126)n A +|(cid:126)e ·(cid:126)n |M v u ds , 2 z 0 z,ij z 0 z,ij i j ∂D Here (cid:126)n is a unit vector normal to the interior facet; (cid:126)n is the outward9 unit 0 normal vector on the domain boundary; and the exact form of the dissipation matricesM , M andM dependsonthechoiceofnumericalflux. Inthispaper, x y z we use a global Lax-Friedrich flux: M =M =M =λI, (23) x y z with λ = 1. This form of numerical flux was chosen over the upwind flux, as was used in [36], because it generates significantly fewer non-zero terms in the variational formulation. The average operator (cid:104)·(cid:105) and jump operator · are defined at any facet for any variable ψ by: (cid:74)(cid:75) ψ++ψ− ψ ≡(ψ+−ψ−) , (cid:104)ψ(cid:105)≡ , (24) (cid:74) (cid:75) 2 with ψ+ and ψ− being defined with respect to the unit normal on the facet, cf.Fig.1. Figure1: Notationfordiscontinuousvariables,givenaunitnormalvector(cid:126)n. 9TheoutwarddirectionisdefinedwithrespecttoD. 8 For each i, 1≤i≤P, the linear form L is given by i (cid:90) L (v )= Q v dx+LBC(v ). (25) i i i i i i D Here LBC accounts for the boundary conditions and will discussed along with i aBC in Section 3.3. For i=P +1 (i.e.the terms associated to Eq.7), we have i (cid:90) (cid:18)E(θ) (cid:16)√ (cid:17)(cid:19) a ((u,θ),v )= −σ 4πu −4πB v dx, (26) P+1 P+1 ∆t a 1 P+1 D (cid:90) E(Tn) L (v )= v dx. (27) P+1 P+1 ∆t P+1 D 3.2. Mass matrix lumping For robustness in optically thick regions, it may be necessary to lump the matrices corresponding to the collision terms. This was demonstrated in [28] in the context of discontinuous Galerkin discretizations of discrete ordinate equa- tions. Inpractice,lumpingamatrixisdonebyreplacingitbyadiagonalmatrix whosei-thtermisthesumoftheelementsonthei-throwoftheoriginalmatrix. For the Crooked Pipe test problem (see Section 5) this lumping proved to be necessary to avoid non-physical instabilities in the solution. 3.3. Initial and boundary conditions Initial conditions for Ih and Th are derived by projecting the initial data for I and T onto VP and V, respectively. Boundary conditions are required for V , but not Th. Unfortunately, the conditions for Ih cannot be derived h directlyfromtheboundaryconditionsforI,sincetheformerrequirefullmoment information and the latter are specified only for incoming data. The boundary conditions that apply to our system are natural, i.e.they are imposed weakly in the variational form by adding appropriate terms to the forms a and L. We impose incoming Dirichlet and reflective boundary conditions on B and B , d r respectively, where B ∪B = ∂D. For χ ∈ {d,r}, we define B− = {((cid:126)r,Ω(cid:126)) ∈ d r χ B ×S2 :(cid:126)n ((cid:126)r)·Ω(cid:126) <0}. The boundary conditions can then expressed as: χ 0 (cid:40) ∀((cid:126)r,Ω(cid:126))∈B−, I((cid:126)r,Ω(cid:126))=g((cid:126)r,Ω(cid:126)), (28) d ∀((cid:126)r,Ω(cid:126))∈B−, I((cid:126)r,Ω(cid:126))=I((cid:126)r,Ω(cid:126) −2(Ω(cid:126) ·(cid:126)n )(cid:126)n ), (29) r 0 0 where g is given. Then aBC = aBC,d +aBC,r, where each term are described i i i below. IncomingDirichletboundary. Dirichletconditionsareimposedbysettingvalues to the incoming data of I on the outward side of the exterior facets. The j outgoing data is obtained by continuity, that is using the outgoing data of I j 9 on the inward side of the exterior facets. The numerical flux (still using a Lax-Friedrich flux) can then be defined on B as: d F(u,g)= (cid:126)ex·(cid:126)n0 (cid:16)A u+H⊕ u+g(cid:9) (cid:17)+ |(cid:126)ex·(cid:126)n0|(cid:0)u−H⊕u−g(cid:9)(cid:1) 2 x flux,x flux,x 2 + (cid:126)ey·(cid:126)n0 (cid:16)A u+H⊕ u+g(cid:9) (cid:17)+ |(cid:126)ey·(cid:126)n0|(cid:0)u−H⊕u−g(cid:9)(cid:1) 2 y flux,y flux,y 2 + (cid:126)ez·(cid:126)n0 (cid:16)A u+H⊕ u+g(cid:9) (cid:17)+ |(cid:126)ez·(cid:126)n0|(cid:0)u−H⊕u−g(cid:9)(cid:1) 2 z flux,z flux,z 2 (30) where we have defined the following half-range integrals for all χ∈{x,y,z}: (cid:90) (cid:90) g(cid:9) ≡ Ω(cid:126) ·(cid:126)e RgdΩ , H⊕ ≡ Ω(cid:126) ·(cid:126)e RRT dΩ, (31) flux,χ χ flux,χ χ S− S+ (cid:90) (cid:90) g(cid:9) ≡ RgdΩ , H⊕ ≡ RRT dΩ, (32) S− S+ where S±((cid:126)r) = {Ω(cid:126) ∈ S2 : ±(cid:126)n ((cid:126)r)·Ω(cid:126) > 0} for all (cid:126)r ∈ ∂D. If (cid:126)n ((cid:126)r) is colinear 0 0 to(cid:126)e ,(cid:126)e or(cid:126)e , the matrices H⊕ and H⊕ can be evaluated exactly using an x y z flux,χ (N +1)-point Gauss-Jacobi quadrature rule. If not, they can be derived using rotation matrices and then applying the quadrature. According to Eq.30, the boundary contribution to a (cf.Eq.19) is (cid:88)P (cid:18)(cid:90) 1(cid:16) (cid:17) aBC,d(u,v )= H⊕ (cid:126)e ·(cid:126)n −H⊕|(cid:126)e ·(cid:126)n | v u ds i i 2 flux,x,ij x 0 ij x 0 i j j=1 Bd (cid:90) 1(cid:16) (cid:17) + H⊕ (cid:126)e ·(cid:126)n −H⊕|(cid:126)e ·(cid:126)n | v u ds (33) 2 flux,y,ij y 0 ij y 0 i j Bd (cid:90) 1(cid:16) (cid:17) (cid:19) + H⊕ (cid:126)e ·(cid:126)n −H⊕|(cid:126)e ·(cid:126)n | v u ds , 2 flux,z,ij z 0 ij z 0 i j Bd while the boundary contribution to L is (cf.Eq.25) (cid:18)(cid:90) 1(cid:16) (cid:17) LBC(v )=− g(cid:9) (cid:126)e ·(cid:126)n −g(cid:9)|(cid:126)e ·(cid:126)n | v ds i i 2 flux,x,i x 0 i x 0 i Bd (cid:90) 1(cid:16) (cid:17) + g(cid:9) (cid:126)e ·(cid:126)n −g(cid:9)|(cid:126)e ·(cid:126)n | v ds (34) 2 flux,y,i y 0 i y 0 i Bd (cid:90) 1(cid:16) (cid:17) (cid:19) + g(cid:9) (cid:126)e ·(cid:126)n −g(cid:9)|(cid:126)e ·(cid:126)n | v ds . 2 flux,z,i z 0 i z 0 i Bd Eq.22 already accounts for the terms in Eq.30 that correspond to the inside of the exterior facet. 10