Positivity-preserving discontinuous Galerkin methods with Lax-Wendroff time discretizations

3.7 MB·

by  Scott A. Moe

Positivity-preserving discontinuous Galerkin methods with
Lax-Wendroff time discretizations
ScottA.Moe1,JamesA.Rossmanith2,andDavidC.Seal3
6
1
1UniversityofWashington,DepartmentofAppliedMathematics,Seattle,WA98195,USA([email protected])
0
2 2IowaStateUniversity,DepartmentofMathematics,396CarverHall,Ames,IA50011,USA
([email protected])
n
a 3U.S.NavalAcademy,DepartmentofMathematics,121BlakeRoad,Annapolis,MD21402,USA
J
([email protected])
9
2

]
A Abstract
N Thisworkintroducesasingle-stage,single-stepmethodforthecompressibleEulerequationsthatisprovablypositivity-
. preserving and can be appliedon both Cartesian and unstructured meshes. Thismethod is the first case of a single-
h
stage,single-stepmethodthatissimultaneouslyhigh-order,positivity-preserving,andoperatesonunstructuredmeshes.
t
a Time-steppingisaccomplishedviatheLax-Wendroffapproach,whichisalsosometimescalledtheCauchy-Kovalevskaya
m
procedure, where temporal derivatives in a Taylor series in time are exchanged for spatial derivatives. The Lax-
[ WendroffdiscontinuousGalerkin(LxW-DG)methoddevelopedinthisworkisformulatedsothatitlookslikeafor-

wardEulerupdatebutwithahigh-ordertime-extrapolatedflux. Inparticular,thenumericalfluxusedinthisworkisa
1
v linearcombinationofalow-orderpositivity-preservingcontributionandahigh-ordercomponentthatcanbedamped
5 toenforcepositivityofthecellaveragesforthedensityandpressureforeachtimestep. Inadditiontothisfluxlim-
4 iter, a moment limiter is applied that forces positivity of the solution at finitely many quadrature points within each
1
cell. Thecombinationofthefluxlimiterandthemomentlimiterguaranteespositivityofthecellaveragesfromone
8
time-step to the next. Finally, a simple shock capturing limiter that uses the same basic technology as the moment
0
. limiterisintroducedinordertoobtainnon-oscillatoryresults. Theresultingschemecanbeextendedtoarbitraryorder
1
withoutincreasingthesizeoftheeffectivestencil. Wepresentnumericalresultsinoneandtwospacedimensionsthat
0
6 demonstratetherobustnessoftheproposedscheme.
1
:
v
i
X 1 Introduction
r
a 1.1 Governingequations
The purpose of this work is to develop a positivity-preserving version of the Lax-Wendroff discontinuous
Galerkin method for the compressible Euler equations on unstructured meshes. The compressible Euler
equationsformasystemofhyperbolicconservationlawthatcanbewrittenasfollows:
   
ρ ρ(cid:126)u
 ρ(cid:126)u  +∇x· ρ(cid:107)(cid:126)u(cid:107)2+p =0. (1)
E (E+p)(cid:126)u
,t
The conserved variables are the mass density, ρ, the momentum density, M(cid:126) =ρ(cid:126)u, and the energy density,
E;theprimitivevariablesarethemassdensity,ρ,thefluidvelocity,(cid:126)u,andthepressure, p. TheenergyE is
1
relatedtotheprimitivevariablesthroughtheequationofstate,
p 1
E = + ρ(cid:107)(cid:126)u(cid:107)2, (2)
γ−1 2
wheretheconstantγistheratioofspecificheats(aka,thegasconstant).
ThecompressibleEulerequationsareanimportantmathematicalmodelinthestudyofgasesandplasma.
Attemptsatnumericallysolvingthetheseequationshasledtoaplethoraofimportanthistoricaladvancesin
thedevelopmentofnumericalanalysisandscientificcomputing(seee.g.,[15,22,28,29,32]).
1.2 DiscontinuousGalerkinspatialdiscretization
The focus of this work is on high-order discontinuous Galerkin (DG) methods, which were originally de-
velopedforgeneralhyperbolicconservationlawsbyCockburn,Shu,etal. inseriesofpapers[10–14]. The
purposeofthissectionistosetthenotationusedthroughoutthepaperandtobrieflydescribetheDGspatial
discretization.
Let Ω⊂Rd be a polygonal domain with boundary ∂Ω. The domain Ω is discretized via a finite set of
non-overlapping elements, Ti, such that Ω=∪Ni=1Ti. Let PMD(cid:0)Rd(cid:1) denote the set of polynomials from Rd
toRwithmaximalpolynomialdegreeM . LetWh denotethebrokenfiniteelementspaceonthemesh:
D
Wh:=(cid:110)wh∈[L∞(Ω)]ME : wh(cid:12)(cid:12)Ti∈(cid:2)PMD(cid:3)ME,∀Ti∈Th(cid:111), (3)
wherehisthemeshspacing. Theaboveexpressionmeansthatwh∈WhhasM components,eachofwhich
E
when restricted to some element T is a polynomial of degree at most M and no continuity is assumed
i D
acrosselementedges(orfacesin3D).
TheapproximatesolutiononeachelementT attimet =tn isoftheform
i
(cid:12) ML(MD)
qh(tn,x(ξξξ))(cid:12) = ∑ Q((cid:96))nϕ((cid:96))(ξξξ), (4)
(cid:12) i
Ti (cid:96)=1
where M is the number of Legendre polynomials and ϕ((cid:96))(ξξξ) : Rd (cid:55)→ R are the Legendre polynomials
L
definedonthereferenceelementT intermsofthereferencecoordinatesξξξ∈T . TheLegendrepolynomials
0 0
areorthonormalwithrespecttothefollowinginnerproduct:
(cid:40)
1 (cid:90) 1 if k=(cid:96),
ϕ(k)(ξξξ)ϕ((cid:96))(ξξξ)dξξξ= (5)
|T0| T0 0 if k(cid:54)=(cid:96).
We note that independent of h, d, M , and the type of element, the lowest order Legendre polynomial is
D
alwaysϕ(1)≡1. ThismakesthefirstLegendrecoefficientthecellaverage:
1 (cid:90) (cid:12)
Q(1)n= qh(tn,x(ξξξ))(cid:12) ϕ(1)(ξξξ)dξξξ=:qn. (6)
i |T0| T0 (cid:12)Ti i
1.3 Timestepping
Themostcommonapproachfortime-advancingDGspatialdiscretizationsisviaexplicitRunge-Kuttatime-
stepping; the resulting combination of time and space discretization is often referred to as the “RK-DG”
method [10]. The primary advantage for this choice of time stepping is that explicit RK methods are easy
toimplement,theycanbeconstructedtobelow-storage,andasubclassofthesemethodshavetheso-called
2
strong stability preserving (SSP) property [23], which is important for defining a scheme that is provably
positivity-preserving. However, there are no explicit Runge-Kutta methods that are SSP for orders greater
thanfour[26,36].
ThemaindifficultywithRunge-Kuttamethodsisthattheytypicallyrequiremanystages;andtherefore,
manycommunicationsareneededpertimestep. Onedirectconsequenceofthecommunicationrequiredat
each RK stage is that is difficult to combine RK-DG with locally adaptive mesh refinement strategies that
simultaneouslyrefineinbothspaceandtime.
The key piece of technology required in locally adaptive DG schemes is local time-stepping (see e.g.,
Dumbser et al. [17]). Local time-stepping is easier to accomplish with a single-stage, single-step (Lax-
Wendroff) method than with a multi-stage Runge-Kutta scheme. For these reasons there is interest from
discontinuous Galerkin theorists and practitioners in developing single-step time-stepping techniques for
DG(seee.g.,[16,18–21,34,42,43]),aswellashybridmultistagemultiderivativealternatives[38].
In this work, we construct a numerical scheme that uses a Lax-Wendroff time discretization that is
coupled with the discontinuous Galerkin spatial discretization. In subsequent discussions in this paper
we demonstrate the advantages of switching to single-stage and single-stage time-stepping in regards to
enforcingpositivityonarbitrarymeshes.
1.4 Positivitypreservation
Insimulationsinvolvingstrongshocks,high-orderschemes(i.e. morethanfirst-order)forthecompressible
Eulerequationsgenerallycreatenonphysicalundershoots(belowzero)inthedensityand/orpressure. These
undershootstypicallycausecatastrophicnumericalinstabilitiesduetoalossofhyperbolicity. Moreover,for
manyapplicationsthesepositivityviolationsexistevenwhentheequationsarecoupledwithwell-understood
totalvariationdiminishing(TVD)ortotalvariationbounded(TVB)limiters. Thechiefgoalofthelimiting
scheme developed in this work is to address positivity violations in density and pressure, and in particular,
toaccomplishthistaskwithahigh-orderscheme.
In the DG literature, the most widely used strategy to maintain positivity was developed by Zhang and
Shuinaseriesofinfluentialpapers[47–49]. ThebasicstrategyofZhangandShuforapositivity-preserving
RK-DGmethodcanbesummarizedasfollows:
Step0. Writethecurrentsolutionintheform
qh(cid:12)(cid:12) =q +θ(cid:16)qh(cid:12)(cid:12) −q(cid:17), (7)
Ti i Ti i
whereθisyet-to-be-determined. θ=1representstheunlimitedsolution.
(cid:12)
Step1. Find the largest value of θ, where 0 ≤ θ ≤ 1, such that qh(cid:12) satisfies the appropriate positivity
Ti
conditionsatsomeappropriatelychosenquadraturepointsandlimitthesolution.
Step2. Findthelargeststabletime-stepthatguaranteesthatwithaforwardEulertimethatthecellaverage
ofthenewsolutionremainspositive.
Step3. Rely on the fact that strong stability-preserving Runge-Kutta methods are convex combinations of
forwardEulertimesteps;andtherefore,thefullmethodpreservesthepositivityofcellaverages(under
someslightlymodifiedmaximumallowabletime-step).
For a Lax-Wendroff time discretization, numerical results indicate that the limiting found in Step 1 is
insufficient to retain positivity of the solution, even for simple 1D advection. Therefore, the strategy we
pursueinthisworkwillstillcontainanequivalentStep1;however,inplaceofStep2andStep3above,we
3
will make use of a parameterized flux, sometimes also called a flux corrected transport (FCT) scheme, to
maintainpositivecellaveragesaftertakingasingletimestep. Indoingso,weavoidintroducingadditional
time step restrictions that often appear (e.g., in Step 2. above) when constructing a positivity-preserving
schemebasedonRunge-Kuttatimestepping.
This idea of computing modified fluxes by combining a stable low-order flux with a less robust high-
order flux is relatively old, and perhaps originates with Harten and Zwas and their self adjusting hybrid
scheme [25]. The basic idea is the foundation of the related flux corrected transport (FCT) schemes of
Boris, Bookandcollaborators[2–5], wherefluxesareadjustedinordertoguaranteethataveragevaluesof
theunknownareconstrainedtoliewithinlocallydefinedupperandlowerbounds. Thisfamilyofmethods
isusedinanextensivevarietyofapplications,rangingfromseismologytometeorology[27,44,46,50]. A
thoroughanalysisofsomeoftheearlymethodsisconductedin[41]. Identicaltomodernmaximumprinciple
preserving(MPP)schemes,FCTcanbeformulatedasaglobaloptimizationproblemwherea“worstcase”
scenario assumed in order to decouple the previously coupled degrees of freedom [1]. Here we do not
attempttouseFCTtoenforceanysortoflocalbounds(inthesenseofdevelopingashock-capturinglimiter),
instead we leverage these techniques in order to retain positivity of the density and pressure associated to
qh(tn,(cid:126)x); such approaches have recently received renewed interest in the context of weighted essentially
non-oscillatory(WENO)methods[6,8,9,30,39,45].
Tosummarize,ourlimitingschemedrawsonideasfromthetwoaforementionedfamiliesoftechniques
that are well established in the literature. First, we start with the now well known (high-order) pointwise
limiting developed for discontinuous Galerkin methods [47–49], and second, we couple this with the very
large family of flux limiters [6, 8, 39]. (developed primarily for finite-difference (FD) and finite-volume
(FV)schemes).
1.5 Anoutlineoftheproposedpositivity-preservingmethod
ThecompressibleEulerequations(1)canbewrittencompactlyas
q +∇·F(q)=0, in Ω⊂Rd, (8)
,t
wheretheconservedvariablesareq=(ρ,M,E)andthefluxfunctionis
 M(cid:126) ·(cid:126)n 
(cid:16) (cid:17)
(cid:126)F·(cid:126)n= M(cid:126) ·(cid:126)n (cid:126)u+p(cid:126)n, (9)
 
(cid:126)u·(cid:126)n(E+p)
whereM(cid:126) =ρ(cid:126)uandE = p +1ρ(cid:107)(cid:126)u(cid:107)2.
γ−1 2
The basic positivity limiting strategy proposed in this work is summarized below. Some important
detailsareomittedhere,butweelaborateonthesedetailsinsubsequentsections.
Step0. Oneachelementwewritethesolutionas
(cid:12) ML(MD)
qh(tn,(cid:126)x((cid:126)ξ))(cid:12) :=qn+θ ∑ Q(k)(t)ϕ(k)((cid:126)x), (10)
(cid:12) i i
Ti k=2
whereθisyet-to-be-determined. θ=1representstheunlimitedsolution.
Step1. Assumethatthissolutionispositiveinthemean. Thatis,weassumeforallithatρ >0and
i
1(cid:107)M(cid:107)2
i
p :=(γ−1)E − >0. (11)
i i 2 ρ
i
4
AconsequenceoftheseassumptionsisthatE >0.
i
Step2. Findthelargestvalueofθ,where0≤θ≤1,suchthatthedensityandpressurearepositiveatsome
suitablydefinedquadraturepoints. Thisstepiselaborateduponin§3.2.
Step3. Construct time-averaged fluxes through the Lax-Wendroff procedure. That is, we start with the
exactdefinitionofthetime-averageflux:
1 (cid:90) tn+1 1 (cid:90) ∆t
Fn((cid:126)x):= F(q(t,(cid:126)x))dt = F(q(tn+s,(cid:126)x))ds (12)
∆t tn ∆t 0
andapproximatethisviaaTaylorseriesexpansionarounds=0:
Fn((cid:126)x):=F(q(tn,(cid:126)x))+∆t dF(q(tn,(cid:126)x))+∆t2d2F(q(tn,(cid:126)x))=Fn((cid:126)x)+O(∆t3). (13)
T 2! dt 3! dt2
Alltimederivativesinthisexpressionarereplacedbyspatialderivativesusingthechainruleandthe
governingPDE(8). Theapproximatetime-averagedflux(13)isfirstevaluatedatsomeappropriately
chosen set of quadrature points – in fact, the same quadrature points as used in Step 2 – and then,
using appropriate quadrature weights, summed together to define a high-order flux, at both interior
and boundary quadrature points. Thanks to Step 2, all quantities of interest used to construct this
expansionarepositiveateachquadraturepoint.
Step4. Timestepthesolutionsothatcellaveragesareguaranteedtobepositive. Thatis,weupdatethecell
averagesviaaformulaoftheform
∆t
Q(1)n+1=Q(1)n− ∑ (cid:126)Fh∗·(cid:126)n , (14)
i i |T| e e
i e∈Ti
where(cid:126)n isanoutward-pointing(relativetoT)normalvectortoedgeewiththepropertythat(cid:107)(cid:126)n (cid:107)is
e i e
thelengthofedgee∈T,andthenumericalfluxonedge,(cid:126)Fh∗,isaconvexcombinationofahigh-order
i e
flux,FH,andalow-orderfluxFL:
e e
(cid:126)Fh∗:=θFH+(1−θ)FL. (15)
e e e
Thelow-orderflux,FL,isbasedonthe(approximate)solutiontotheRiemannproblemdefinedbycell
e
averagesonly,andthe“high-order”flux,FH,isconstructedafterintegratingviaGaussianquadrature
e
the(approximate)Riemannsolutionsatquadraturepointsalongtheedgee∈T:
i
M
1 Q
FH= ∑ω FH, (16)
e 2 k ek
k=1
whereω aretheGaussianquadratureweightsforquadraturewithM pointsandFH arethenumer-
k Q ek
icalfluxesateachoftheM quadraturepoints. Notethatthissumhasonlyasinglesummandinthe
Q
one-dimensionalcase. Theselectionofθisdescribedinmoredetailin§3.2. Thisstepguaranteesthat
thesolutionretainspositivity(inthemean)forasingletimestep.
Step5. Applyashock-capturinglimiter. Thepositivity-preservinglimiterisdesignedtopreservepositivity
of the solution, but it fails at reducing spurious oscillations, and therefore a shock-capturing limiter
needs to be added. There are many choices of limiters available; we use the limiter recently devel-
oped in [31] because of its ability to retain genuine high-order accuracy, and its ability to push the
polynomialordertoarbitrarydegreewithoutmodifyingtheoverallscheme.
5
Step6. Repeatallofthesestepstoupdatethesolutionforthenexttimestep.
Each step of this process is elaborated upon throughout the remainder of this paper. The end result is
thatourmethodisthefirstschemetosimultaneouslyobtainallofthefollowingproperties:
• High-orderaccuracy. Theproposedmethodisthird-orderinspaceandtime,andcanbeextendedto
arbitraryorder.
• Positivity-preserving. The proposed limiter is provably positivity-preserving for the density and
pressure,atafinitesetofpointvalues,fortheentiresimulation.
• Single-stage, single-step. We use a Lax-Wendroff discretization for time stepping the PDE, and
thereforeweonlyneedonecommunicationpertimestep.
• Unstructuredmeshes. BecauseweusethediscontinuousGalerkinmethodforourspatialdiscretiza-
tionandallofourlimitersaresufficientlylocal,weareabletorunsimulationswithDG-FEMonboth
Cartesianandunstructuredmeshes.
• No additional time-step restrictions. Because we do not rely on a SSP Runge-Kutta scheme, we
do not have to introduce additional time-step restrictions to retain positivity of the solution. This
differentiatesusfrompopularpositivity-preservinglimitersbasedonRKtimediscretizations[48].
1.6 Structureofthepaper
The remainder of this paper has the following structure. The Lax-Wendroff DG (LxW-DG) method is
described in §2, where we view the scheme as a method of modified fluxes. The positivity-preserving
limiter is described in §3, where the discussion of the limiter is broken up into two parts: (1) the moment
limiter(§3.1)and(2)theparameterizedfluxlimiter(§3.2). In§4wepresentnumericalresultsonseveraltest
casesin1D,2DCartesian,and2Dunstructuredmeshes. Finallyweclosewithconclusionsandadiscussion
offutureworkin§5.
2 TheLax-WendroffdiscontinuousGalerkinscheme
2.1 Thebasescheme: Amethodofmodifiedfluxes
TheLax-WendroffdiscontinuousGalerkin(LxW-DG)method[34]servesasthebaseschemeforthemethod
developedinthiswork. ItistheresultofanapplicationoftheCauchy-Kovalevskayaproceduretohyperbolic
PDE: we start with a Taylor series in time, then we replace all time derivatives with spatial derivatives via
the PDE. Finally, a Galerkin projection discretizes the overall scheme, where a single spatial derivative is
reservedforthefluxesinordertoperformintegration-by-parts.
We review the Lax-Wendroff DG scheme for the case of a general nonlinear conservation law that is
autonomousinspaceandtimeinmultipledimensions[34]. Thecurrentpresentationillustratesthefactthat
Lax-Wendroff schemes can be viewed as a method of modified fluxes, wherein higher-order information
about the PDE is directly incorporated by simply redefining the fluxes that would typically be used in an
“Eulerstep.”
Weconsideragenericconservationlawoftheform
q +∇·F(q)=0, (17)
,t
6
where the matrix ∂F ·nˆ is diagonalizable for every unit length vector nˆ and q in the domain of interest.
∂q
Formalintegrationof(17)overaninterval[tn,tn+1]resultsinanexactupdatethrough
q(t+∆t,(cid:126)x)=q(t,(cid:126)x)−∆t∇·F(q(t,(cid:126)x)), (18)
wherethetime-averaged flux[7]isdefinedas
1 (cid:90) tn+∆t
F(q(t,(cid:126)x)):= F(q(t,(cid:126)x))dt. (19)
∆t tn
Moreover,aTaylorexpansionofFandachangeofvariablesyields
1 (cid:90) ∆t(cid:18) 1 (cid:19)
F(q)= F(qn)+τF(qn) + τ2F(qn) +··· dτ
,t ,t,t
∆t 2
0 (20)
1 1
=F(qn)+ ∆tF(qn) + ∆t2F(qn) +···,
,t ,t,t
2! 3!
which can be inserted into (18). In a numerical discretization of (18), the Taylor series in (20) is truncated
afterafinitenumberofterms.
Remark1. IfF≈F(q(tn,(cid:126)x)),then(18)reducestoaforwardEulertimediscretizationforhyperboliccon-
servationlaw(17). Thisfactwillallowustoincorporatepositivity-preservinglimitersintotheLax-Wendroff
fluxconstruction.
Thisobservationallowsustoincorporatethepositivity-preservinglimitersthatarepresentedin§3.1and
§3.2,becauseweviewtheLxW-DGmethodasamethodofmodifiedfluxes.
2.2 Constructionofthetime-averagedflux
Wenowdescribehowtocomputethetemporalderivativeterms:
F(qn) , F(qn) , F(qn) , ... (21)
,t ,t,t ,t,t,t
that are required to define the time-averaged flux in (20). This discussion is applicable to high-order finite
difference methods, finite volume methods (e.g., ADER), as well as discontinuous Galerkin finite element
methods.
Asingleapplicationofthechainruletocomputethetimederivativeofthefluxfunctionyields
∂F
=F(cid:48)(q)·q =−F(cid:48)(q)·(∇·F), (22)
,t
∂t
wherethefluxJacobianis
∂F
F(cid:48)(q) := i, 1≤i,j≤M. (23)
ij
∂q
j
The matrix-vector products in (22) can be compactly written using the Einstein summation convention
(whererepeatedindicesareassumedtobesummedover),whichproducesavectorwhoseith-componentis
∂F ∂F ∂q ∂F
i i j i
= =− (∇·F) . (24)
∂t ∂q ∂t ∂q j
j j
7
Asecondderivativeof(22)yields
∂2F = ∂ (cid:0)−F(cid:48)(q)∇·F(cid:1)=F(cid:48)(cid:48)(q)·(∇·F(q),∇·F(q))+F(cid:48)(q)·∇(cid:0)F(cid:48)(q)(∇·F)(cid:1), (25)
∂t2 ∂t
whereF(cid:48)(cid:48)(q)istheHessianwithelementsgivenby
∂2F (cid:18) ∂2f ∂2g (cid:19)
i i i
F := = , . (26)
ijk
∂q ∂q ∂q ∂q ∂q ∂q
j k j k j k
Equations (22) and (25) are generic formulae; the equalities are appropriate for any two-dimensional hy-
perbolic system, and similar identities exist for three dimensions. The first product in the right hand side
of (25) is understood as a Hessian-vector product. Scripts that compute these derivatives, as well as the
matrix,andHessianvectorproductsthatarenecessarytoimplementathird-orderLax-Wendroffschemefor
multidimensionalEulerequationscanbefoundintheopensourcesoftwareFINESS[37].
Finally, these two time derivatives are sufficient to construct a third-order accurate method by defining
thetime-averagedfluxthrough
1 1
Fn(q):=F(qn)+ ∆tF(qn) + ∆t2F(qn) , (27)
T 2! ,t 3! ,t,t
andthenupdatingthesolutionthrough
qn+1((cid:126)x)=qn((cid:126)x)−∆t∇·Fn (28)
T
inplaceof(18).
2.3 Fully-discreteweakformulation
The final step is to construct a fully discrete version of (28). The LxW-DG scheme follows the following
process[24,38]:
Step1. Ateachquadraturepointevaluatethenumericalflux, F(qn), andthenintegratethisnumericalflux
againstbasisfunctionstoobtainaGalerkinexpansionofFh insideeachelement.
Step2. UsingtheGalerkinexpansionsofqhandFh,evaluateallrequiredspatialderivativestoconstructthe
timeexpansionFn in(27)ateachquadraturepoint.
T
Step3. Multiply (28) by a test function ϕ((cid:96)), integrate over a control element T, and apply the divergence
i
theoremtoyield
(cid:90) (cid:90) (cid:90) (cid:73)
qn+1ϕ((cid:96))dx= qnϕ((cid:96))dx−∆t ∇ϕ((cid:96))·F (qn)dx+∆t ϕ((cid:96))F (qn)·nˆds, (29)
T T
Ti Ti Ti ∂Ti
wherenˆ istheoutwardpointingunitnormaltoelementT,whichreducesto
i
∆t (cid:90) ∆t (cid:73)
Q((cid:96))n+1=Q((cid:96))n− ∇ϕ((cid:96))·Fh dx+ ϕ((cid:96))Fh∗·nˆds (30)
i i |T| T |T| T
i Ti i ∂Ti
(cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125)
Interior Edges
byorthogonalityofthebasisϕ. Inpractice, boththeinteriorandedgeintegralsareapproximatedby
appropriate numerical quadrature rules. The flux values, Fh∗, in the edge integrals still need to be
T
defined.
8
ρBC ρb
M(cid:126)BC M(cid:126)b
EBC Eb
−→ ←−
∂Ω
Fig.1: Interior,qb,andexterior,qBC,solutionvaluesoneithersideoftheboundary∂Ω.
Step4. Along each edge solve Riemann problems at each quadrature point by using the left and right
interfacevalues. Inthisworkweusethewell-knownLax-Friedrichsflux:
(cid:16) (cid:17) 1(cid:104) (cid:16) (cid:16) (cid:17) (cid:16) (cid:17)(cid:17) (cid:16) (cid:17)(cid:105)
Fh∗ qh,qh ·nˆ = nˆ· F qh +F qh −s qh −qh , (31)
T − + 2 T + T − + −
where s is an estimate of the maximum global wave speed, qh is the approximate solution evaluated
−
on the element boundary on the interior side of T, and qh is the approximate solution evaluated on
i +
theelementboundaryontheexteriorsideofT.
i
2.4 Boundaryconditions
Inordertoachievehigh-orderaccuracyattheboundariesofthecomputationaldomain, acarefultreatment
ofthesolutionineachboundaryelementisrequired. Inparticular,allsimulationsinthisworkrequireeither
reflective(hardsurface)ortransparent(outflow)boundaryconditions.
Supposethesolutiontakesonthevalue
qb=(cid:0)ρb,M(cid:126)b,Eb(cid:1) (32)
at a quadrature point xb on the boundary ∂Ω, and we wish to define a boundary value, qBC, on the exterior
side of the boundary ∂Ω that yields a flux with one of two desired boundary conditions. This is depicted
in Figure 1. Let ˆt and nˆ be the unit tangent and unit normal vectors to the boundary at the boundary point
xb, respectively. In both the reflective and transparent boundary conditions, we enforce continuity of the
tangentialcomponents:
M(cid:126)BC·ˆt=M(cid:126)b·ˆt. (33)
The only difference between the two types of boundary conditions we consider lie in the normal direction.
Weset
M(cid:126)BC·nˆ =∓M(cid:126)b·nˆ, (34)
wheretheminussigncorrespondstothereflectiveboundaryconditionandtheplussigncorrespondstothe
transparentboundarycondition.
FromthiswecaneasilywriteoutthefullboundaryconditionsatthepointxBC:
ρBC  ρb 
(cid:16) (cid:17) (cid:16) (cid:17)
M(cid:126)BC= M(cid:126)b·ˆt ˆt∓ M(cid:126)b·nˆ nˆ, (35)
 
EBC Eb
whereagaintheminussigncorrespondstothereflectiveboundaryconditionandtheplussigncorresponds
tothetransparentboundarycondition.
9
In order to achieve high-order time accuracy we need apply the above boundary conditions to the time
derivativesontheboundary:
ρB,tC (cid:16) (cid:17) ρb,t(cid:16) (cid:17)  ρB,t,Ct  (cid:16) (cid:17) ρb,t,(cid:16)t (cid:17) 
M(cid:126),BtC= M(cid:126),bt·ˆt ˆt∓ M(cid:126),bt·nˆ nˆ, M(cid:126),Bt,Ct= M(cid:126),bt,t·ˆt ˆt∓ M(cid:126),bt,t·nˆ nˆ, (36)
E,BtC E,bt E,Bt,Ct E,bt,t
wherealltimederivativesmustbereplacedbyspatialderivativesusingthePDE.
3 Positivitypreservation
The temporal evolution described in the previous section fails to retain positivity of the solution, even in
the simple case of linear advection with smooth solutions that are near zero. In fact, in extreme cases the
projection of the initial conditions can fail to retain positivity of the solution due to the Gibbs phenomena.
Thepositivity-preservinglimiterwepresentfollowsatwostepprocedure:
Step1. Limitthemomentsintheexpansionsothatthesolutionispositiveateachquadraturepoint.
Step2. Limitthefluxessothatthecellaveragesretainpositivityafterasingletimestep.
Wenowdescribethefirstofthesetwosteps.
3.1 Positivityatinteriorquadraturepointsviamomentlimiters
Thissectiondescribesaprocedurethatimplementsthefollowing: ifthecellaveragesarepositive,thenthe
solution is forced to be positive at a preselected and finite collection of quadrature points. We select only
thequadraturepointsthatareactuallyusedinthenumericalupdate;thisincludesinternalGaussquadrature
pointsaswellasface/edgeGaussquadraturepoints. Unlikeotherpositivitylimitingschemes[47,49]inthe
SSP Runge-Kutta framework, this step is not strictly necessary to guarantee positivity of the cell average
atthenexttime-step;however,themainreasonforapplyingthelimiteratquadraturepointsistoguarantee
thateachtermintheupdateisphysical,whichwillreducethetotalamountofadditionallimitingofthecell
average updated needed in Section 3.2. The process to maintain positivity at quadrature points is carried
out in a series of three simple steps. Because this part of the limiter is entirely local, we focus on a single
elementT ,andthereforedropthesubscriptforeaseofnotation.
k
3.1.1 Step0: Assumepositivityofthecellaverages.
We assume that the cell averages for the density satisfies ρn ≥ε , where ε >0 is a cutoff parameter that
0 0
guarantees hyperbolicity of the system. In this work, we set ε =10−12 in all simulations. Furthermore,
0
weassumethatthecellaverageforthepressuresatisfies pn>0,wherethe(average)pressure pn isdefined
throughtheaveragesoftheotherconservedquantitiesin(11). Notethatthesetwoconditionsaresufficient
n
toimplythattheaverageenergyE ispositive.
3.1.2 Step1: Enforcepositivityofthedensity.
In this step, we enforce positivity of the density at each quadrature point x ∈T . Because ρn ≥ε , there
m k 0
existsa(maximal)valueθρ∈[0,1]suchthat
ML(MD)
ρθ :=ρn+θ ∑ ρ((cid:96))nϕ((cid:96))(x )≥ε (37)
m m 0
(cid:96)=2
10