Table Of ContentA conservative semi-Lagrangian HWENO method for the Vlasov
equation∗
Xiaofeng Cai† Jianxian Qiu‡ Jingmei Qiu§
February 2, 2016
6
1
0
2 Abstract: In this paper, we present a high order conservative semi-Lagrangian (SL) Hermite weighted
n
a essentially non-oscillatory (HWENO) method for the Vlasov equation based on dimensional splitting [Cheng
J
0 and Knorr, Journal of Computational Physics, 22(1976)]. The major advantage of HWENO reconstruction,
3
compared with the original WENO reconstruction, is compact. For the split one-dimensional equation, to
]
A
ensure local mass conservation, we propose a high order SL HWENO scheme in a conservative flux-difference
N
.
h form, following the work in [J.-M. Qiu and A. Christlieb, Journal of Computational Physics, v229(2010)].
t
a
m Besidesperformingdimensionalsplittingfortheoriginal2Dproblem, wedesignapropersplittingforequations
[
of derivatives to ensure local mass conservation of the proposed HWENO scheme. The proposed fifth order
1
v
SL HWENO scheme with the Eulerian CFL condition has been tested to work well in capturing filamentation
4
7
0 structures without introducing oscillations. We introduce WENO limiters to control oscillations when the time
0
0
stepping size is larger than the Eulerian CFL restriction. We perform classical numerical tests on rigid body
.
2
0 rotation problem, and demonstrate the performance of our scheme via the Landau damping and two-stream
6
1
instabilities when solving the Vlasov-Poisson system.
:
v
i
X Keywords: Conservative semi-Lagrangian scheme; HWENO reconstruction; Vlasov-Poisson system, Lan-
r
a dau damping, Two-stream instability.
1 Introduction
This paper focuses on a high order conservative semi-Lagrangian scheme with high order HWENO recon-
struction for the Vlasov-Poisson (VP) simulations based on dimensional splitting. The VP system, arise from
∗ResearchwaspartiallysupportedbyNSFCgrants91230110,11328104,11571290,91530107andNSFDMS-1217008andDMS-
1522777.
†SchoolofMathematicalSciences,XiamenUniversity,Xiamen,Fujian,361005,P.R.China. E-mail: [email protected].
‡SchoolofMathematicalSciencesandFujianProvincialKeyLaboratoryofMathematicalModeling&High-PerformanceScien-
tificComputing,XiamenUniversity,Xiamen,Fujian,361005,P.R.China. E-mail: [email protected].
§DepartmentofMathematics,UniversityofHouston,Houston,77204. E-mail: [email protected].
1
collisionless plasma applications, reads as following,
∂f
+v·∇ f +E(t,x)·∇ f =0, (1.1)
∂t x v
and
E(t,x)=−∇ φ(t,x), −∆ φ(t,x)=ρ(t,x), (1.2)
x x
where x and v are coordinates in phase space (x,v) ∈ R3×R3, E is the electric field, φ is the self-consistent
electrostaticpotentialandf(t,x,v)isprobabilitydistributionfunctionwhichdescribestheprobabilityoffinding
a particle with velocity v at position x at time t. The probability distribution function couples to the long
(cid:82)
rangefieldsviathechargedensity,ρ(t,x)= f(t,x,v)dv−1,wherewetakethelimitofuniformlydistributed
R3
infinitely massive ions in the background. Equations (1.1) and (1.2) have been nondimensionalized so that all
physical constants are one.
PopularmethodsinfusionsimulationsincludeLagrangian,semi-LagrangianandEulerianmethods. Popular
lagrangian methods include the particle-in-cell (PIC) [2, 15, 18], Lagrangian particle methods [3, 12]; Eulerian
methodsincludeweightedessentiallynon-oscillatory(WENO)coupledwithFouriercollocation[38], continuous
finite element methods [37, 36], Runge-Kutta discontinuous Galerkin methods [1, 11, 17, 8]. Each method
has its own advantages and limits. For example, Lagrangian methods are well known for its reasonable low
computational cost for high dimensional problems. However, it suffers from statistical noise due to the initial
sampling of macro-particles. Eulerian methods offer a good alternative to overcome this lack of precision, but
they suffer from ’the curse of dimensionality’ and the CFL time step restriction. Compared with the Eulerian
approach, theSLmethodsisrelieffromtheCFLtimesteprestriction, becauseinformationisbeingpropagated
along characteristics.
AmongSLschemes,theschemebasedondimensionalsplitting,introducedbyChengandKnorroriginally[7],
isverypopularSLschemewithhighordercubicsplineinterpolationwasproposedin[31]. Apositivitypreserving
and flux conservative finite volume SL scheme with ENO reconstruction for the VP system is proposed in [14]
and the scheme for the guiding center Vlasov model is proposed in [10]. A conservative finite different semi-
LagrangianschemewithWENOreconstructionisproposedin[25]; laterthealgorithmisgeneralizedtovariable
coefficient case [26, 27] and maximum principle preserving limiter [33]. In the finite element discontinuous
Galerkinframework,thereareSLdiscontinuousGalerkinschemeswithpositivtypreservinglimiters[28,29]and
hybrid SL finite element-finite difference methods in [16]. High order propagation methods based on Hermite
2
interpolationareproposedin[22,13,5,4,35]. HWENOschemewasintroducedin[23]andfurtherdevelopedin
[39,21]forhyperbolicconservationlaws. Besidestheoriginalequation,onealsoevolvesequationsofderivatives
in the WENO fashion. Hence their reconstruction stencils are more compact than the original WENO scheme
[19], given the same order of approximation. A similar technique, called the CIP (Constrained Interpolation
Profile/Cubic Interpolated Propagation) scheme [34], has also been proposed for the VP system [22]. In [5],
besides the function values themselves, the gradients of the function are evolved in a semi-Lagrangian fashion.
In [35], a semi-Lagrangian HWENO is proposed without evolving the gradients. The proposed method in
this paper successfully couples the semi-Lagrangian framework with HWENO method. Compared with those
earlier work, our method achieves local mass conservation, has relatively compact stencil and is able to capture
under-resolve solution structures without numerical artifacts.
In this paper, we design a conservative SL HWENO scheme based on dimensional splitting. Its stencil is
more compact than the regular WENO scheme. We follow the idea in [25] to express the SL update in a flux
difference conservative form. There are several new ingredients we developed in this paper in order to ensure
the effectiveness and robustness of the proposed scheme: firstly, we design a SL HWENO scheme in a flux
differenceformfor1Dproblem; secondly, weapplyaspecialsplittingsimilartothesplittingintheCIPmethod
[22] for mass conservation for high dimensional problems; thirdly, we introduce WENO limiters for controlling
oscillations.
The paper is organized as follows. In Section 2, we introduce a conservative form of high order SL HWENO
method for 1D transport problems. In Section 3, we introduce a conservative SL HWENO for VP system by
a special form of splitting, and introduce WENO limiters. In Section 4, we present our numerical results for
basic test problems, such as linear advection and rigid body rotation, and for the VP simulations. Concluding
remarks are given in Section 5.
2 Conservative SL HWENO method for 1D transport problem
Inthissection,weintroducetheSLHemiteinterpolationinaflux-differenceformfor1Dtransportproblems.
Then we incorporate the HWENO mechanism into the flux function reconstruction procedure to realize a non-
oscillatory capturing of shocks.
3
2.1 The SL Hemite interpolation in a flux-difference form for 1D transport prob-
lem
In this section, we consider a 1D transport problem
f +vf =0, f(x,t=0)=f (x), on [a,b], (2.1)
t x 0
where v is a constant. For the Hermite method, we also consider the evolution equation for the solution’s
.
derivative g =f . For the linear transport problem (2.1), g satisfies the same linear transport equation
x
g +vg =0.
t x
We discretize the domain [a,b] as
a=x <x <···<x =b,
1 3 N+1
2 2 2
withtheuniformgridpointsx =a+(i−1)∆xandthecellsize∆x=x −x . WeletI =[x ,x ]and
i 2 i+1 i−1 i i−1 i+1
2 2 2 2
I = [x ,x ]. In the HWENO approach, the numerical solutions associated with each grid point are point
i+1 i i+1
2
values fn and derivatives gn. Here the subscript i means the solution at the grid point x and the superscript
i i i
n means the solution at time level tn. To design a SL HWENO scheme, we update {fn+1,gn+1}N from the
i i i=1
corresponding solutions at time tn.
For the linear problem (2.1) with constant characteristic speed, the solutions fn+1 and gn+1 can be obtain
i i
by shifting the information at tn in the SL framework. We define the amount of shift v∆t as xshift. There are
∆x
threecasesofxshift: shifttotherightbysomeamountlessthanhalfacell(xshift∈[0,1]),shifttotheleftby
2
some amount less than half a cell (xshift∈[−1,0]) and shift a distance greater than half a cell (|xshift|> 1).
2 2
Toillustratetheidea,weonlyconsiderHermiteinterpolationsforxshift∈[0,1],whiletheoneforxshift∈
2
[−1,0] is mirror symmetric with respect to x of the previous interpolations. In the case when |xshift| > 1,
2 i 2
whole grid shifting is carried out and followed by a final update based on the procedure for xshift ∈ [−1,1].
2 2
In the following, we present the Hermite interpolation with cubic polynomials. Higher order schemes will be
discussed later.
1. Theunderlyingfunctionattn canbeapproximatedbyaHermite-typereconstruction,basedonthestencil
{fn ,fn,gn ,gn},
i−1 i i−1 i
f(cid:101)n (ξ) = fn−gn∆xξ+(cid:0)2gn∆x−3fn+3fn +gn ∆x(cid:1)ξ2
i−1 i i i i i−1 i−1
2
+(cid:0)2fn−gn∆x−2fn −gn ∆x(cid:1)ξ3,
i i i−1 i−1
4
where ξ(x)= xix−−1−xixi ∈[0,1],x∈Ii−21.
2. fn+1 and gn+1 can be obtained by tracing the characteristic back to time t = tn and evaluating the
i i
interpolant f(cid:101)in−1(ξ) at the foot of characteristics x(cid:63) =xi−v∆t,
2
fn+1 = fn−ξ ((3fnξ −2fnξ2)−(3fn ξ −2fn ξ2))
i i 0 i 0 i 0 i−1 0 i−1 0
−gn∆xξ +(2gn∆x+gn ∆x)ξ2+(−gn∆x−gn ∆x)ξ3 (2.2)
i 0 i i−1 0 i i−1 0
(cid:18) 6fn−6fn (cid:19)
gn+1 = gn+ −4gn+ i i−1 −2gn ξ
i i i ∆x i−1 0
(cid:18) 6fn−fn (cid:19)
+ − i i−1 +3gn+3gn ξ2, (2.3)
∆x i i−1 0
where ξ0 = xxi−(cid:63)−1−xxii.
For the above linear scheme (2.2) and (2.3), we have the following mass conservation result.
PROPOSITION 1. If (cid:80)N gn ≡ 0 and with periodic boundary condition, then the scheme (2.2) and (2.3)
i=1 i
conserve the total mass, i.e., (cid:80)N fn+1 ≡(cid:80)N fn and (cid:80)N gn+1 ≡0.
i=1 i i=1 i i=1 i
Proof:
N N
(cid:88) (cid:88)
fn+1 = [fn−ξ ((3fnξ −2fnξ2)−(3fn ξ −2fn ξ2))
i i 0 i 0 i 0 i−1 0 i−1 0
i=1 i=1
−gn∆xξ +(2gn∆x+gn ∆x)ξ2+(−gn∆x−gn ∆x)ξ3]
i 0 i i−1 0 i i−1 0
N
(cid:88)
= fn−ξ ((3fnξ −2fnξ2)−(3fnξ −2fnξ2))
i 0 N 0 N 0 0 0 0 0
i=1
N N
(cid:88) (cid:88)
+ gn ∆xξ2− gn ∆xξ3
i−1 0 i−1 0
i=1 i=1
N
(cid:88)
= fn (periodic boundary condition),
i
i=1
(cid:88)N (cid:88)N (cid:20) (cid:18) 6fn−6fn (cid:19)
gn+1 = gn+ −4gn+ i i−1 −2gn ξ
i i i ∆x i−1 0
i=1 i=1
(cid:18) 6fn−fn (cid:19) (cid:21)
+ − i i−1 +3gn+3gn ξ2
∆x i i−1 0
(cid:88)N 6fn −6fn 6fn −fn (cid:88)N (cid:88)N
= gn+ N 0 ξ − N 0 ξ2+2 gn ∆xξ2−3 gn ∆xξ3
i ∆x 0 ∆x 0 i−1 0 i−1 0
i=1 i=1 i=1
N
(cid:88)
= gn (periodic boundary condition).
i
i=1
Hence the proposition is proved.
5
In order to guarantee (cid:80)N g0 ≡ 0 in the assumption of the proposition, we introduce a sliding average
i=1 i
function h(x) in [30] which satisfies
1 (cid:90) x+∆2x
f(x)= h(ξ)dξ,
∆x
x−∆x
2
then
(cid:18) (cid:18) (cid:19) (cid:18) (cid:19)(cid:19)
1 ∆x ∆x
g(x)=f(x) = h x+ −h x− .
x ∆x 2 2
Thus (cid:80)N g0 = (cid:80)N (h0 −h0 ) ≡ 0 where h0 ≈ h(x± ∆x) can be obtained by reconstruction from
i=1 i i=1 i+1 i−1 i±1 2
2 2 2
{f0} .
j j
In the following, we will adopt a matrix notation for presentation of the Hermite interpolation. The matrix
A will denote the interpolation matrix. We use A(i,j) to denote the element at the ith row and jth column,
A(i,:) to denote the ith row of A, and A(:,j) to denote the jth column of A.
We rewrite the scheme (2.2) and (2.3) into a flux difference form, in order to ensure local mass conserva-
tion, especially when the nonlinear HWENO mechanism is applied. In order to do so, we propose to update
(cid:16) (cid:17)
{fn,hn } insteadof{fn,gn} ,observingthatgn canberecoveredfrom{hn } bygn = hn −hn /∆x.
i i+1 i i i i i i+1 i i i+1 i−1
2 2 2 2
Specifically, (2.2) can be rewritten in the following flux difference form using the new {hn } ,
i+1 i
2
fn+1 = fn−ξ ((3fnξ −2fnξ2)−(3fn ξ −2fn ξ2))
i i 0 i 0 i 0 i−1 0 i−1 0
−gn∆xξ +(2gn∆x+gn ∆x)ξ2+(−gn∆x−gn ∆x)ξ3
i 0 i i−1 0 i i−1 0
= fn−ξ (fn(3ξ −2ξ2)−fn (3ξ −2ξ2))+gn∆x(−ξ +2ξ2−ξ3)+gn ∆x(ξ2−ξ3)
i 0 i 0 0 i−1 0 0 i 0 0 0 i−1 0 0
= fn−ξ (fn(3ξ −2ξ2)−fn (3ξ −2ξ2))
i 0 i 0 0 i−1 0 0
(cid:16) (cid:17) (cid:16) (cid:17)
− hn −hn ξ (1−2ξ +ξ2)− hn −hn ξ (−ξ +ξ2)
i+1 i−1 0 0 0 i−1 i−3 0 0 0
2 2 2 2
(cid:110)(cid:104) (cid:105)
= fn−ξ fn(3ξ −ξ2)+hn (1−2ξ +ξ2)+hn (−ξ −ξ2)
i 0 i 0 0 i+1 0 0 i−1 0 0
2 2
(cid:104) (cid:105)(cid:111)
− fn (3ξ −ξ2)+hn (1−2ξ +ξ2)+hn (−ξ −ξ2)
i−1 0 0 i−1 0 0 i−3 0 0
2 2
= fi−ξ0(f(cid:98)in+1(ξ0)−f(cid:98)in−1(ξ0)).
2 2
where
f(cid:98)in−1(ξ0)=(fin−1,hni−1,hni−3)·C3L·(1,ξ0,ξ02)(cid:48)
2 2 2
with
0 3 −2
C3L = 1 −2 1 .
0 −1 1
6
We update gn+1 by
i
hn+1−hn+1
gn+1 = i+12 i−12 (2.4)
i ∆x
where
hn+1 =(fn ,hn ,hn )·DL·(1,ξ ,ξ2)(cid:48) (2.5)
i−1 i−1 i−1 i−3 3 0 0
2 2 2
with
0 6 −6
D3L = 1 −4 3 .
0 −2 3
REMARK 1. The flux-difference form for the SL finite difference scheme was originally proposed in [25].
There are two main advantages to work with the flux difference form:
1. The flux difference form can ensure local mass conservation.
2. We can design a nonlinear HWENO mechanism for the flux reconstructions, see discussions in the next
subsection.
In order to work with the flux difference form, it is crucial to work with the {hn } instead of the original
i+1 i
2
derivative function g =f .
x
REMARK 2. We observe that DL(:,k)=kCL(:,k), k =1,2,3.
3 3
REMARK 3. The case presented here is for the third order scheme. Similar procedure can be used to obtained
higher order scheme, e.g. the fifth order case with HWENO is presented in the next subsection.
2.2 HWENO reconstruction for flux functions
In general, high order fixed stencil reconstruction of numerical fluxes performs well when the solution is
smooth. However, around discontinuities, oscillations will be introduced. In this subsection, a nonlinear SL
HWENOprocedureisintroducedforreconstructingthefluxf(cid:98)n (ξ). Byadaptivelyassigningnonlinearweights
i−1
2
to neighboring candidate stencils, the nonlinear HWENO reconstruction preserves high order accuracy of the
linear scheme around smooth regions of the solution, while producing a sharp and essentially non-oscillatory
capture of discontinuities.
WeadopttheideaoftheHWENOreconstruction[23,21]intotheproposedconservativeSLframework. We
present a fifth order HWENO reconstruction as an example. Similar procedure can be generalized to higher
order case.
7
Our discussion will be focused on the case of xshift ∈ [−1,1]. As before, the case of |xshift| > 1 will be
2 2 2
handledwithawholegridshiftfollowedbythecaseofxshift∈[−1,1]toaccountforthefractionalremainder.
2 2
When xshift ∈ [0,1], the fifth order conservative SL method the {fn ,fn ,fn,fn ,gn ,gn } is the
2 i−2 i−1 i i+1 i−2 i+1
following,
(cid:18) (cid:19)
8 19 2 1
fn+1 = fn+ − fn +fn − fn + gn ∆x− gn ∆x ξ
i i 27 i−2 i−1 27 i+1 9 i+1 9 i−2 0
(cid:18) (cid:19)
1 2 7 19 25
+ − gn ∆x− gn ∆x− fn− fn +fn + fn ξ2
18 i−2 9 i+1 4 i 108 i−2 i−1 27 i+1 0
(cid:18) (cid:19)
1 1 1 5 3 1
+ gn ∆x− gn ∆x+ fn+ fn − fn + fn ξ3
6 i−2 6 i+1 4 i 12 i−2 4 i−1 12 i+1 0
(cid:18) (cid:19)
1 2 3 19 1 23
+ gn ∆x+ gn ∆x+ fn+ fn − fn − fn ξ4
18 i−2 9 i+1 4 i 108 i−2 2 i−1 54 i+1 0
(cid:18) (cid:19)
1 1 1 13 1 13
+ − gn ∆x− gn ∆x− fn− fn + fn + fn ξ5,
18 i−2 18 i+1 4 i 108 i−2 4 i−1 108 i+1 0
hn −hn
Using the flux difference form for g function, gn = i+12 i−12 in tn, then
i ∆x
fn+1 =fn−ξ ((fn ,fn ,f ,fn ,hn ,hn ,hn ,hn )·BL·(1,ξ ,ξ2,ξ3,ξ4)(cid:48))
i i 0 i−2 i−1 i i+1 i−5 i−3 i+1 i+3 5 0 0 0 0
2 2 2 2
where
8 19 − 5 − 19 13
27 108 12 108 108
−1 −1 3 1 −1
0 7 −41 −23 14
B5L = −12971 −−422157 −11412 251344 −−1411038 .
19 118 −61 −181 118
29 −182 −61 218 −181
9 9 6 9 18
−2 2 1 −2 1
9 9 6 9 18
Then fn+1 can be written in the flux difference form,
i
fn+1 =fn−ξ ((fn ,fn,fn ,hn ,hn )·CL−(fn ,fn ,fn,hn ,hn )·CL)·(1,ξ ,ξ2,ξ3,ξ4)(cid:48) (2.6)
i i 0 i−1 i i+1 i−3 i+3 5 i−2 i−1 i i−5 i+1 5 0 0 0 0
2 2 2 2
where
− 8 − 19 5 19 − 13
27 108 12 108 108
19 89 −1 − 35 7
CL = 2179 −10285 − 31 21308 −5143
5 217 127 −112 −541 1108
9 18 6 18 18
−2 2 1 −2 1
9 9 6 9 18
And we have the flux difference form the derivative gn+1 in tn+1,
i
hn+1−hn+1
gn+1 = i+12 i−12
i ∆x
where
hn+1 =(fn ,fn ,fn,hn ,hn )·DL·(1,ξ ,ξ2,ξ3,ξ4)(cid:48) (2.7)
i−1 i−2 i−1 i i−5 i+1 5 0 0 0 0
2 2 2
8
where DL satisfies DL(:,k)=kCL(:,k), k =1,··· ,5.
5 5 5
When xshift∈[−1,0], we update {fn+1,hn+1} by the following formulas,
2 i i+1 i
2
fin+1 =fin+ξ0(f(cid:98)in+1(ξ0)−f(cid:98)in−1(ξ0)), (2.8)
2 2
where the flux function
f(cid:98)in−1(ξ)=(fin−1,fin,fin+1,hi−3,hi+3)·C5R·(1,ξ0,ξ02,ξ03,ξ04)(cid:48), (2.9)
2 2 2
where
19 −25 − 1 23 − 13
27 27 12 54 108
19 89 −1 − 35 7
CR = −278 −10189 53 11908 −5143 . (2.10)
5 −227 1208 112 −1028 1108
9 9 6 9 18
1 1 −1 − 1 1
9 18 6 18 18
And
hn+1 =(fn ,fn,fn ,hn ,hn )·DR·(1,ξ ,ξ2,ξ3,ξ4)(cid:48) (2.11)
i−1 i−1 i i+1 i−3 i+3 5 0 0 0 0
2 2 2
where DR satisfies DR(:,k)=kCR(:,k), k =1,··· ,5.
5 5 5
In the following, we illustrate the corresponding HWENO reconstruction of flux functions. We only discuss
the HWENO reconstruction for the flux f(cid:98)n and hn+1 when xshift ∈ [0,1]. When xshift ∈ [−1,0], the
i−1 i−1 2 2
2 2
flux f(cid:98)in−1 and hni−+11 could be reconstructed symmetrically with respect to xi. From equations (2.6) and (2.7),
2 2
the stencil {fn ,fn ,fn,hn ,hn } is used to construct the flux f(cid:98)n (ξ) and hn+1. It is composed of the
i−2 i−1 i i−5 i+1 i−1 i−1
2 2 2 2
information from three potential stencils
S ={hn ,fn ,fn }, S ={fn ,fn ,fn}, S ={fn ,fn,hn }. (2.12)
1 i−5 i−2 i−1 2 i−2 i−1 i 3 i−1 i i+1
2 2
Intuitively, in regions where the function is smooth, we want to use information from S ,S and S in an
1 2 3
optimal way, to obtain a fifth order approximation. On the other hand, around a big jump, we only want to
use the information from the relatively smooth stencil. Following [25], we only use the HWENO mechanism in
adaptively reconstructing the coefficients in front of the constant 1 in the equation for f(cid:98)n and hn+1, while
i−1 i−1
2 2
leaving coefficients for ξ ,ξ2,ξ3,ξ4 unchanged. We can observe that the first column of matrix CL is the same
0 0 0 0 5
as that of DL. Thus we only consider the HWENO procedure for constructing f(cid:98)n ,
5 i−1
2
1. Compute the linear weights, γ ,γ and γ , such that
1 2 3
(fn ,fn ,fn,hn ,hn )·CL(:,1)
i−2 i−1 i i−5 i+1 5
2 2
1 5 1
=γ (fn ,fn ,hn )·(−2,2,1)(cid:48)+γ (fn ,fn ,fn)·(− , , )(cid:48)
1 i−2 i−1 i−52 2 i−2 i−1 i 6 6 3
1 5 1
+γ (fn ,fn,hn )·( , ,− )(cid:48),
3 i−1 i i+21 4 4 2
9
where(fn ,fn ,hn )·(−2,2,1)(cid:48), (fn ,fn ,fn)·(−1,5,1)(cid:48) and(fn ,fn,hn )·(1,5,−1)(cid:48) arethird
i−2 i−1 i−5 i−2 i−1 i 6 6 3 i−1 i i+1 4 4 2
2 2
orderreconstructionsoffluxesfromthreestencilsS ,S andS respectively. Fromequation(2.6),γ = 1,
1 2 3 1 9
γ = 4 and γ = 4.
2 9 3 9
2. We compute the smoothness indicator, denoted by β , for each stencil S , which measures how smooth
j j
the function p (x) is in the target cell I . The smaller this smoothness indicator β , the smoother the
j i j
function p (x) is in the target cell. We use the same recipe for the smoothness indicator as in [19],
j
(cid:88)2 (cid:90) (cid:18) ∂ (cid:19)2
β = ∆x2l−1 p (x) dx.
j ∂xl j
l=1 Ii
Intheactualnumericalimplementationthesmoothnessindicatorsβ arewrittenoutexplicitlyasquadratic
j
forms of the points {fn,hn } in the stencil,
i i+1 i
2
13(cid:18) 9 3 3 (cid:19)2 (cid:18)31 9 13 (cid:19)2
β = − fn + hn + fn + fn − hn − fn ,
1 3 4 i−2 2 i−25 4 i−1 4 i−2 2 i−25 4 i−1
β = 13(cid:0)−fn +2fn −fn(cid:1)2+(cid:18)−3fn+2fn − 1fn (cid:19)2,
2 12 i−2 i−1 i 2 i i−1 2 i−2
13(cid:18) 9 3 3 (cid:19)2 (cid:18)5 1 3 (cid:19)2
β = − fn+ fn + hn + fn+ fn − hn .
3 3 4 i 4 i−1 2 i+12 4 i 4 i−1 2 i+21
3. We compute the nonlinear weights based on the smoothness indicators.
ω γ
ω = j , j =1,2,3, ω = k
j (cid:80)3k=1ωk k (cid:15)+βk
where (cid:15) is a small number to prevent the denominator to becoming zero. In our numerical tests we take
(cid:15) to be 10−6.
4. Compute numerical fluxes constructed in HWENO fashion. Define the matrix C(cid:101)L and D(cid:101)L as,
5 5
1 5 1 1 5 1
D(cid:101)5L(:,1)=C(cid:101)5L(:,1)=ω1·(−2,2,0,1,0)+ω2·(−6,6,3,0,0)+ω3·(0,4,4,0,−2)
C(cid:101)L(:,k)=CL(:,k), D(cid:101)L(:,2)=kCL(:,k), k =2,··· ,5.
5 5 5 5
The updated numerical flux is computed using C(cid:101)L and D(cid:101)L, i.e.,
5 5
f(cid:98)in−1(ξ0)=(fin−2,fin−1,fin,hni−5,hni+1)·C(cid:101)5L·(1,ξ0,ξ02,ξ03,ξ04)(cid:48), (2.13)
2 2 2
hni−+11 =(fin−2,fin−1,fin,hni−5,hni+1)·D(cid:101)5L·(1,ξ0,ξ02,ξ03,ξ04)(cid:48). (2.14)
2 2 2
10
