Implicit Filtered PN for High-Energy Density Thermal Radiation Transport using Discontinuous Galerkin Finite Elements

2.3 MB·

by  Vincent M. Laboure

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.:[phone]
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