Dimension Splitting and a Long Time-Step Multi-Dimensional Scheme for Atmospheric Transport

6.5 MB·

by  Yumeng Chen

DimensionsplittingforAdvection 1
Dimension Splitting and a Long Time-Step Multi-Dimensional
Scheme for Atmospheric Transport
Yumeng Chen,a Hilary Weller,b∗ Stephen Pringc and James Shawb
aUniversitätHamburg
7
1 bMeteorology,UniversityofReading
0
2 cMetOffice,UK

n
a ∗Correspondenceto:E-mail:<[email protected]>
J

4
2

] Dimensionally split advection schemes are attractive for atmospheric modelling due to their
A
efficiencyandaccuracyineachspatialdimension.Accuratelongtimestepscanbeachieved
N
.
h without significant cost using the flux-form semi-Lagrangian technique. The dimensionally
t
a
splitschemeusedinthispaperisconstructedfromtheone-dimensionalPiecewiseParabolic
m
[ Method and extended to two dimensions using COSMIC splitting. The dimensionally split

1
scheme is compared with a genuinely multi-dimensional, method of lines scheme with
v
7
implicittime-steppingwhichisstableforverylargeCourantnumbers.
0
9
Two-dimensional advection test cases on Cartesian planes are proposed that avoid the
6
0
complexities of a spherical domain or multi-panel meshes. These are solid body rotation,
.
1
0 horizontal advection over orography and deformational flow. The test cases use distorted
7
1 mesheseithertorepresentslopingterrainortomimicthedistortionsofacubedsphere.
:
v
Such mesh distortions are expected to accentuate the errors associated with dimension
i
X
r splitting,however,thedimensionallysplitschemeisveryaccurateonorthogonalmeshesand
a
accuracydecreasesonlyalittleinthepresenceoflargemeshdistortions.Thedimensionally
splitschemealsolosessomeaccuracywhenlongtime-stepsareused.Themulti-dimensional
scheme is almost entirely insensitive to mesh distortions and asymptotes to second-order
accuracy at high resolution. As is expected for implicit time-stepping, phase errors occur
whenusinglongtime-stepsbutthespatiallywellresolvedfeaturesareadvectedatthecorrect
speedandthemulti-dimensionalschemeisalwaysstable.
A naive estimate of computational cost (number of multiplies) reveals that the implicit
scheme is the most expensive, particularly for large Courant numbers. If the multi-
dimensional scheme is used instead with explicit time-stepping, the Courant number is
restrictedtolessthanonebutthecostbecomessimilartothedimensionallysplitscheme.
(cid:13)c 0000RoyalMeteorologicalSociety Preparedusingqjrms4.cls
KeyWords: Multi-dimensional;advection;stability;accuracy;longtime-steps
Received...
QuarterlyJournaloftheRoyalMeteorologicalSociety Q.J.R.Meteorol.Soc.00:2–24(0000)
1. Introduction
Weatherandclimatemodelsarebeingdevelopedonquasi-uniformmeshesinordertobetterexploitmoderncomputers(eg ?Lauritzen
et al. 2014; ?; Skamarock and Gassmann 2011; K.K. et al. 2015) and so accurate and efficient transport (or advection) schemes on
non-orthogonalmeshesarerequired.Thereisanabundanceofdesirablepropertiesofadvectionschemes,including:
1. Inherentlocalconservationoftheadvectedquantity.
2. StabilityinthepresenceoflargeCourantnumbers.
3. AccuracyinthepresenceoflargeCourantnumbers.
4. Highorder-accuracy.
5. Lowcomputationalcost,goodparallelscalingandmulti-tracerefficiency.
6. Lowphaseanddispersionerrors(advectionofallwavenumbersoftheadvectedquantityatclosetothecorrectspeed).
7. Lowdiffusionerrors(maintainingamplitudeofallwavenumbersoftheadvectedquantity).
8. Boundedness,monotonicity,positivityandmaintainingcorrelationsbetweenmultipleadvectedtracers.
Four(important)propertiesarelistedtogetherinthefinalitembecausetheywillnotbeaddressedhere.Inthispaper,weaddressthe
issueofconservative,accurate,efficientadvectionschemesforlogicallyrectangular,non-orthogonalmesheswhicharestableinthe
presenceoflargeCourantnumbers.Theseschemeswouldbeparticularlyrelevantforcubed-spheremeshesandforterrainfollowing
meshes.Anothernovelaspectofthispaperisthat,forsimplicity,wetestadvectionschemesentirelyonplanarmeshesratherthanon
thesphere,proposingtestcasestochallengeadvectionschemesonnon-orthogonalmesheswithouttheneedtoimplementmeshesin
sphericalgeometry.
Dimensionallysplitschemes(operatingseparatelyineachspatialdimension)areattractiveforatmosphericmodellingduetotheir
efficiencyandhighaccuracyineachspatialdimension(eg. LinandRood1996a;Leonardetal.1996a;BrassingtonandSanderson
1999; Putman and Lin 2007; K.K. et al. 2015). Inherent conservation is guaranteed by using the flux-form semi-Lagrangian (FFSL
or forward in time) technique (eg. Colella and Woodward 1984a) which integrates the dependent variable over a swept distance
upstreamofeveryfaceinordertocalculatefluxesinandoutofcells.Accuratelongtimestepscanbeachievedwithoutsignificantcost
bycalculatingcumulativemassfluxesalongthedomain(?).Thisisdonewithoutremappingwhichcanbeanexpensiveprocedure,
finding a conservative map between fields on overlapping meshes. One-dimensional schemes can be used with operator splitting
to create dimensionally split, second-order accurate schemes (eg Leonard et al. 1996a) on logically rectangular, multidimensional
meshes. Dimensionally split schemes have been found to give good accuracy on non-orthogonal meshes such as the cubed-sphere
(Putman and Lin 2007; K.K. et al. 2015) with special treatment over cube edges. Putman and Lin (2007) use the average of two
one-sidedschemesatcubeedgeswhereasK.K.etal.(2015)createghostcellsoutsideeachcubepanelboundary.Withoutthisspecial
treatment,dimensionallysplitschemesdonotaccountproperlyformeshdistortionanderrorsmayoccur.Accuracybetweensecond
andfourthorderwasfoundinpracticeonarangeoftestcasesonthecubed-spherebyK.K.etal.(2015).
The same problems occur when using dimensionally split schemes over orography since terrain following layers become non-
orthogonaloverorographyandspecialtreatmentcannotbeusedeverywherewherethereisorography.Acommonsolutionhasbeen
to make the terrain following layers as smooth as possible, reducing non-orthogonality (eg ?) or to use floating Lagrangian vertical
co-ordinates(?).Dimensionsplittingmayaccountforsomeoftheerrorsoverorographyreportedby?althoughtheydonotcitethisas
areasonforerrors.
Dimensionsplittingerrorsondistortedmeshescanbeeliminatedbyusinggenuinelymulti-dimensionaladvectionschemes.These
can be either FFSL (ie swept area, eg. Lashley 2002; Lipscomb and Ringler 2005; Miura 2007; Thuburn et al. 2014), method of
lines (MOL, discretising space and time separetely, eg. Weller et al. 2009; Skamarock and Gassmann 2011; K.K. et al. 2015; ?) or
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
DimensionsplittingforAdvection 3
conservativesemi-Lagrangian(withconservativere-mapping,egIskeandKaser2004;Lauritzenetal.2010;Zerroukatetal.2004).
TheFFSLandMOLmulti-dimensionalschemeshavenotpreviouslybeenextendedtoworkwithCourantnumberssignificantlylarger
thanone.FFSLmulti-dimensionalschemescouldbeextendedtohandlelargetime-stepsbyintegratingtheupstreamsweptvolumeover
alargeupstreamvolume,interactingwithalargenumberofupstreamcells.Howeverthecostwouldbeproportionaltothetime-step
sincethelargertheCourantnumber,themorecellstheupstreamsweptvolumewouldneedtooverlapwith.Thistechniquetherefore
offersnoadvantageoverusingasmallertime-step.MOLmulti-dimensionalschemescanbeextendedtoworkwithCourantnumbers
larger than one by using implicit time-stepping. This will increase the computational cost per tracer advected since the solution of
amatrixequationwouldbeneededforeveryadvectedtracer.Twootherdisadvantagesofimplicittime-steppingarethelargephase
errorswhenlongtime-stepsareused(eg??)andthedifficultyofachievingmonotonicity.Thisisincontrasttosemi-Lagrangianor
FFSL schemes which maintain accuracy with long time-steps (???) although monotonocity with long time-steps is still challenging
(?).Conservativesemi-Lagrangiannaturallyextendstolongtime-stepsbuttheconservativeremappingiscomplicatedandexpensive,
particularlyonnon-rectangularmeshesandwillnotbeinvestigatedhere.Lauritzenetal.(2014)describedhowtheFFSLtechnique
withalongtime-stepcanbemadeequivalenttotheconservativesemi-Lagrangian.
Itisthereforenotclearwhatapproachshouldbetakenforachievinglongtime-stepswhenadvectingmultipletracersondistorted
meshes. In this paper, we show the effect of dimension splitting errors using a FFSL dimensionally split scheme on a number of
test-caseswhichusedistortedmeshesandcomparewithagenuinelymulti-dimensionalimplicitMOLschemeusinglargeandsmall
Courantnumbers.
Thetheoreticalpropertiesofdimensionallysplitadvectionschemesareoftentestedonuniform,orthogonalmeshes(egLeonardetal.
1996a).Onthecubed-sphere,specialtreatmentisneededforthecubeedges(egLinandRood1996a;K.K.etal.2015).Developinga
transportschemetotheextentthatitcanbeusedonamulti-panelcubed-spherewithspecialtreatmentofcubeedgesisaconsiderable
undertaking. Hence there is a need for more challenging advection test cases which are simpler to implement, without the need for
spherical meshes. We therefore propose some modifications of existing test cases to use distorted meshes, or distorted co-ordinate
systems,onalogicallyrectangular,two-dimensionalplane.
Thelongtime-steppermitting,dimensionallysplitschemeandthelongtime-steppermittingmulti-dimensionalschemearedefined
insection2.Insection3wepresentresultsofthreeadvectiontestcasesondistortedmeshesintwo-dimensionalplanesusingCourant
numbers above and below one. These are the solid body rotation test case of Leonard et al. (1996a) modified to use a mesh (or
co-ordinate system) with distortions similar to a cubed-sphere (section 3.1), the horizontal advection test case over orography (?),
examiningsensitivitytotime-step,resolutionandmountainheight,allonthemaximallydistortedbasicterrainfollowingmesh(section
3.2)andamodificationofthedeformationalflowtestcaseofLauritzenetal.(2012b)foraperiodicrectangularplane(section3.3).
Someestimatesaremadeofcomputationalcostinsection3.4andfinalconclusionsaredrawninsection4.
2. TransportSchemes
Wepresenttwoconservativeadvectionschemessuitableforlongtime-steps(stableforCourantnumberssignificantlylargerthanone)
forsolvingthelinearadvectionequation:
∂φ
+∇·uφ =0 (1)
∂t
wherethedependentvariableφ isadvectedbyvelocityfieldu(x,t).Thedimensionallysplitschemeisthepiecewiseparabolicmethod
(PPMColellaandWoodward1984a)whichusestheflux-formsemi-Lagrangianapproachextendedtolongtime-stepsfollowing?with
COSMICsplitting(Leonardetal.1996a)toextendPPMtotwodimensions.Themulti-dimensionalschemeusesthemethodoflines
approach (treating space and time independently). The second-order accurate spatial discretistaion of Weller and Shahrokhi (2014)
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
4 Y.Chen,H.Weller,S.PringandJ.Shaw
is combined with implicit, Crank-Nicholson time-stepping to allow long time-steps. Neither scheme has monotonicity or positivity
preservation.
Thecodefordimensionallysplittingschemeisavailableathttps://github.com/yumengch/COSMIC-splitting.The
multi-dimensionalschemeisimplementedusing?andisavailableathttps://github.com/AtmosFOAM/.
2.1. One-dimensionalPPMwithLongTime-steps
Wedescribealongtime-stepversionofPPMforsolvingtheone-dimensionaladvectionequation:
∂φ ∂uφ
+ =0. (2)
∂t ∂x
Colella and Woodward (1984a) defined PPM with monotonicity constraints and for variable resolution but for simplicity (and for
comparisonwiththemulti-dimensionalscheme)wewilldefinePPMwithoutmonotonicityconstraintsandforafixedresolution,∆x,
andtime-step,∆t.Withtheserestrictions,PPMshouldbefourth-orderaccurateinonedimensionforvanishingtime-step.Wedefine
(n)
thedependentvariable,φi ,tobethemeanvalueofφ incelliattime-levelnwherexi=i∆xandt=n∆t.SincePPMisaflux-form
(n)
finitevolumemethod,φ isupdatedusing:
i
u φ −u φ
φ(n+1)=φ(n)+X (φ)=φ(n)− i+1/2 i+1/2 i−1/2 i−1/2 (3)
i i C i ∆x
where X (φ) is the conservative advection operator for φ in the x direction. The fluxes, u φ , are found by integrating a
C i±1/2 i±1/2
piecewisepolynomial, p,alongthedistancetravelledineachtime-stepupwindofcellboundaryx .Thepolynomialisdefinedin
i±1/2
eachcell,i,suchthat:
ˆ
1 xi+1/2
φi= pi(x)dx (4)
∆x
xi−1/2
by
(cid:18) (cid:18) (cid:19)(cid:19)
1(cid:0) (cid:1)
pi(x)=pi−1/2+ξ pi+1/2−pi−1/2+(1−ξ)6 φi−2 pi−1/2+pi+1/2 (5)
whereξ =(x−x )/∆xand
i−1/2
7 1
pi+1/2= 12(φi+φi+1)−12(φi+2+φi−1) (6)
Inordertocopewithlongtime-steps,wefollow?anddividetheCourantnumberintoasignedintegerpart,c ,andaremainder,c .
N r
Thedeparturepointoflocationx isthuscomputedasx =x −u ∆tandthedeparturecell,forcelledgei−1,isi =i−c −1
i−1 d i−1 i−1 2 d N
2 2 2
foru >0andi =i−c foru <0.Thenforu >0thefluxthroughx betweentimesn∆t and(n+1)∆t is:
i−1 d N i−1 i−1 i−1/2
2 2 2
ˆ (cid:32) ˆ (cid:33)
1 xi−1/2 1 xi−1/2−cN
u φ = p(x)dx= M −M + p(x)dx (7)
i−21 i−12 ∆t xd ∆t i−1/2 i−cN−1/2 xi−1/2−cN−cr∆x
whereM isthecumulativemassfromthestartpointtopositionx :
i−1/2 i−1/2
M =∑∆xφ . (8)
i−1/2 k
k<i
Thisdeparturepointcalculationassumesthatthevelocityisuniformonthecomputationalmeshwhichhasafirst-ordererrorwhich
couldbeparticularlydamagingforlongCourantnumbers,whenthewrongdeparturecellcouldbefound.
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
DimensionsplittingforAdvection 5
ThevelocityisderivedfromastreamfunctionandtheJacobianoftheco-ordinatetransform:
   
u Ψ
y
 =J  (9)
   
v −Ψ
x
Forstability,thetime-stepisrestrictedbythedeformationalCourantnumber:
(cid:18)(cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12) (cid:12)(cid:19)
(cid:12)∂u(cid:12) (cid:12)∂u(cid:12) (cid:12)∂u(cid:12) (cid:12)∂v(cid:12) (cid:12)∂v(cid:12) (cid:12)∂v(cid:12) (cid:12)∂w(cid:12) (cid:12)∂w(cid:12) (cid:12)∂w(cid:12)
cd=∆tmax (cid:12)(cid:12)∂x(cid:12)(cid:12),(cid:12)(cid:12)∂y(cid:12)(cid:12),(cid:12)(cid:12)∂z(cid:12)(cid:12),(cid:12)(cid:12)∂x(cid:12)(cid:12),(cid:12)(cid:12)∂y(cid:12)(cid:12),(cid:12)(cid:12)∂z(cid:12)(cid:12),(cid:12)(cid:12)∂x(cid:12)(cid:12),(cid:12)(cid:12)∂y(cid:12)(cid:12),(cid:12)(cid:12)∂z(cid:12)(cid:12) (10)
(?)suchthatc ≤1.
d
2.2. COSMICSplitting
COSMICoperatorsplitting(Leonardetal.1996a)allowssinglestage,one-dimensionalschemessuchasPPMtobeusedstablyintwo
ormoredimensionswhilstretainingconservation,constancypreservationandsecond-orderaccuracy(onorthogonalmeshes).Aswe
arenowconsideringtwospatialdimensions,wedefineφij,uijandvij,thevaluesofφ andthevelocitycomponents,uandvincell(i,j)
wherex=i∆xandy= j∆y.COSMICsplittingusesbothadvectiveandconservativeadvectionoperatorsinthexandydirections:
1 1
XC(φ)=− (ueφe−uwφw) YC(φ)=− (vnφn−vsφs) (11)
∆x ∆y
φij φij
XA(φ)=Xc(φ)+ (ue−uw) YA(φ)=Yc(φ)+ (vn−vs) (12)
∆x ∆y
whereφn=φi,j+1/2,φs=φi,j−1/2,φe=φi+1/2,j,φw=φi−1/2,j,vn=vi,j+1/2,vs=vi,j−1/2,ue=ui+1/2,janduw=ui−1/2,jarethevalues
of φ, u and v at the cell boundaries. If COSMIC is being used to extend PPM to two spatial dimensions then φn,s,e,w are calculated
fromequation7.AssumingC-gridstaggering,v ,v ,u andu aredependentvariables.Insteadofusingcellcenteredvelocity(Lin
n s e w
andRood1996a)orupwindvelocity(Leonardetal.1996a),theadvectiveoperatorsarecalculatedinasimilarmannerto?.
Meshdistortionscanbeincludedintheadvectionequationwithaco-ordinatetransformwithJacobianJ:
∂|J|−1φ ∂|J|−1uφ ∂|J|−1vφ
+ + =0. (13)
∂t ∂x ∂y
TheoperatorsX ,Y X andY arethen:
c c A A
XC(φ)=−∆1x(cid:0)|J|−e1ueφe−|J|−w1uwφw(cid:1) YC(φ)=−∆1y(cid:0)|J|−n1vnφn−|J|−s1vsφs(cid:1) (14)
XA(φ)=Xc(φ)+φ∆ixj(cid:0)|J|−e1ue−|J|−w1uw(cid:1) YA(φ)=Yc(φ)+φ∆iyj(cid:0)|J|−n1vn−|J|−s1vs(cid:1) (15)
whicharecombinedtoupdateφij ineachcellby:
(cid:18) |J| (cid:16) (cid:17)(cid:19) (cid:18) |J| (cid:16) (cid:17)(cid:19)
φi(jn+1)=φinj+|J|ijXC φi(jn)+ 2ijYA φi(jn) +|J|ijYC φi(jn)+ 2ijXA φi(jn) . (16)
where|J| isthedeterminantofJacobianatthecenterofcell(i,j),and∆xand∆ygivethecellsizeinuniformcomputationaldomain.
ij
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
6 Y.Chen,H.Weller,S.PringandJ.Shaw
φ
φ φ φ
y down φ
φ φ
φ φ f x
y φdown
φ φ
up up
φ φ φf x
φ φ
φ
φ
φ
φ
φ
φ
φ
Figure1.Stencilsofupwind-biasedcellsforinterpolatingontofacefusinga2dcubicpolynomialfortwodifferentmeshstructures.
2.3. Multi-dimensionalMethodofLines(MOL)Scheme
The MOL scheme uses the finite volume method for arbitrary meshes and is implemented in (?). This uses a cubic upwind spatial
discretisation (Weller and Shahrokhi 2014; ?) combined with Crank-Nicholson in time. Although the interpolation uses a cubic
polynomial,cellcentrevaluesareapproximatedascellaveragevaluesandfacecentrevaluesareapproximatedasfaceaveragevalues
sothemethodislimitedto2ndorderaccuracyinspace.TheadvectionschemeusesGauss’sdivergencetheoremtoapproximatethe
divergencetermoftheadvectionequation:
1
∇·uφ ≈ ∑φ u ·S (17)
f f f
V
f∈c
whereV isthecellvolume,thesummationisoverallfaces, f,ofcellc,φ isthevalueofthedependentvariable,φ interpolatedonto
f
face f,u isthevelocityatface f andS isthevectornormaltoface f withmagnitudeequaltotheareaofface f (iethefacearea
f f
vector).Theadvectionschemeusesafittoa2d(or3d)polynomialusinganupwind-biasedstencilofcells(fig1)inordertointerpolate
fromknowncellvaluesontoaface.In2d,thecubicpolynomialis:
φ =a+bx+cy+dx2+exy+fy2+gx3+hx2y+ixy2 (18)
omittingthey3term,wherexisthedirectionnormaltoacellfaceandyisperpendiculartox.(They3termisomittedbecauseitcannot
besetwithastencilthatisnarrowinthedirectionoftheflow.)Coefficientsatoiaresetfromaleastsquaresfittothecelldatainthe
stencil.Theleast-squaresprobleminvolvesa9×mmatrixsingularvaluedecomposition(wheremisthesizeofthestencil)forevery
faceandforbothorientationsofeachface.Howeverthisispurelyageometriccalculationandisthereforeapre-processingactivity
sincethemeshisfixed.Thisgeneratesasetofweightsforcalculatingφ fromthecellvaluesinthestencil,leavingmmultipliesfor
f
eachfaceforeachcalloftheadvectionoperator.Thestencilsarefoundforthree-dimensional,arbitrarilystructuredmeshesbyfinding
theface(s)closesttoupwindofthefaceweareinterpolatingonto,takingthetwocellseithersideoftheupwindface(s)andthentaking
thevertexneighboursofthosecentralcells(fig1).Foreachfacetherearetwopossiblestencilsdependingontheupwinddirection.
Bothofthestencilsarestoredandtheinterpolationweightsforbothstencilsarecalculated.
Inordertoensurethatthefitisaccurateinthecellseithersideofface f andtoensurethatthevaluesintheseadjacentcellshave
thestrongestcontroloverφ ,rowsassociatedwiththesevaluesintheleastsquaresfitmatrixareweightedafactorof1000relativeto
f
theotherrows(followingLashley2002).Thisdoesnotaffecttheorderofaccuracy.Mathematically,anarbitrarilylargevalueofthe
weightcanbeusedtoensurethatthefitgoesexactlythroughtheupwindanddownwindcell.Howeverifavaluetoolargeisused,the
singularvalueproblembecomesillconditioned.Wedonotuseanyoftheotherstabilisationproceduresasdescribedby?.Thevalue
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
DimensionsplittingforAdvection 7
φ isthencalculatedasahigher-ordercorrectiontofirst-orderupwind:
f
φf =φup+∑wcφc (19)
c
where w are the weights for each cell of the stencil calculated from the least squares fit, with w reduced by 1 to make the fit a
c up
correctiononupwind.
2.3.1. ImplicitSolution,MatrixSolversandTolerances
ThetrapezoidalimplicitorCrank-Nicholsontime-steppingleadstoamatrixequationwhichneedstobesolvedtofindalltheφsatthe
nexttime-step.Inordertoensurethatthematrixisdiagonallydominantforarbitrarytime-steps,thecubicinterpolationappliedisa
deferredcorrectiononfirst-orderupwindsothatonlythecoefficientcorrespondingtotheupwindcellsareincludedinthematrix.This
meansthatmorethanoneimplicitsolvesareneededpertime-stepsothatthehigher-ordertermsaresolvedtobesecond-orderaccurate
intime.IftheCourantnumberislessthanorclosetoone,weusetwoimplicitsolvespertime-step.Consequently,assumingthatthe
velocityfieldandmeshareconstantintime,thetime-steppingschemeisdefinedas:
φ(cid:48)−φn 1 (cid:18) (cid:19)
∆t = −2V ∑ φunp+φu(cid:48)p+2∑wcφcn uf·Sf (20)
f∈c c
φn+1−φn 1 (cid:18) (cid:19)
∆t = −2V ∑ φunp+φunp+1+∑wcφcn+∑wcφc(cid:48) uf·Sf. (21)
f∈c c c
ForlargerCourantnumbersweusefourimplicitsolvespertime-stepalthoughsensitivitytothischoicehasnotbeeninvestigated.
Ifthisimplicitschemeisappliedonalogicallyrectangular,two-dimensionalmeshwithhorizontalandverticalCourantnumbers
c andc ,thenthediagonalcoefficientsofthematrixwouldbe1+c /2+c /2and,assumingtwoupwinddirections,therewouldbe
x z x z
exactlytwooff-diagonalelements,−c /2and−c /2.Consequentlythematrixisverysparse,asymmetricanddiagonallydominantfor
x z
alltime-steps.ItissolvedusingtheOpenFOAMbi-conjugategradientsolverusingDILUpre-conditioningtoatoleranceof10−8every
iteration.Sensitivitytothesolverorsolvertolerancehavenotbeeninvestigated.Informationaboutthenumberofsolveriterationsfor
differenttestcasesingiveninsection3.4.
3. ResultsofTestCasesinPlanarGeometry
Inordertomaketestcasesassimpleaspossiblewithouttheneedforincorporatingsphericalgeometry,multi-panelmeshesornon-
rectangularcells,ourcomputationaldomainconsistsofaperiodictwo-dimensionalplanewithdeformationsintheco-ordinatesystem
(ormesh)mimickingthekindofdistortionswhichareproducedbyacubed-spheremesh.Wealsouseatwo-dimensionalverticalslice
testcaseoverorographyusingterrainfollowingco-ordinates(orterrainfollowingmeshes).
A range of test cases are undertaken on uniform and distorted meshes using the dimensionally split scheme and the multi-
dimensional scheme in order to evaluate the influences of mesh distortions, the validity of using a dimensionally split scheme on a
distortedmeshandtheschemes’accuracyandstabilityforlongtime-steps.
3.1. SolidBodyRotation
The solid body rotation test case of Leonard et al. (1996a) is used to compare the accuracy of the dimensionally split and multi-
dimensional schemes on orthogonal and non-orthogonal meshes. We define this test case on a domain that is 104×104m2. The
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
8 Y.Chen,H.Weller,S.PringandJ.Shaw
velocityisdefinedbynumericallydifferentiatingthestreamfunctionwhichcanbedefinedatmeshvertices,x,by:
ψ(x,t)=A|x−xc|2 (22)
where xc is the centre of the domain, and A=5π/3000s−1 so that the angular velocity is 2A. The initial tracer takes a Gaussian
distributioninordertoensurethatalladvectionschemesachievetheirtheoreticalorderofaccuracy:
(cid:32) (cid:33)
1|x−x |2
φ(x)=exp − φ (23)
2 r2
φ
wherex =x +r jistheinitialcentreofthetracerdistribution,r =2500m,r =500mandjistheunitvectorintheydirection.
φ c cφ cφ φ
Theanalyticsolutionhasthesametracerdistributionbutwiththecentreofthetracerat:
 
cos(π/2+2At)
x =x +r   (24)
φ c cφ 
sin(π/2+2At)
andtheGaussianrotatesanti-clockwiseexactlyonerevolutionin600s.
The solid body rotation test case is performed on a uniform, orthogonal mesh and on a non-orthogonal mesh on a plane with
non-orthogonalitysimilartothatofacubed-spheremesh.Forthedimensionallysplitscheme,non-orthogonalityisachievedusingthe
co-ordinatetransform:

(cid:16) (cid:17)
X=x Y =ym 1+2yym−−ff fory≥ f (25)
(cid:16) (cid:17)
ym 1+y−ff fory< f
wherey =5000mand f istheequationforaypositionofuniformY halfwayupthedomain.Inordertocreateanglesof120ointhe
m
mesh,similartoacubed-sphere, f isgivenby:

(cid:16) (cid:17)
ym 1+ √1 −√x forx≤xm
f = 2 3 3 . (26)
(cid:16) (cid:17)
ym 1− √1 +x−√xm forx>xm
2 3 3
where x =5000 m. For a 50×50 mesh this gives the x and y co-ordinate locations as shown in figure 2. The multi-dimensional
m
schememodelusesCartesianco-ordinatesandadistortedmeshratherthananon-orthogonalco-ordinatesystemonaCartesianmesh.
Howeverthisdoesnotaffectthenumericalresultsassumingthattheco-ordinatetransformsareimplementedinaconsistentwaytothe
distortedmeshinCartesianco-ordinates.
Forthedimensionallysplitscheme,bi-periodicboundaryconditionsareapplied. Forthemulti-dimensionalscheme,itwasmore
straightforwardtoimposefixedvalueboundaryconditionsofφ =0wherethevelocityisintothedomainandzeronormalgradient
wherethevelocityisoutofthedomain.Howeverφ remainsalmostzeroneartheboundariessotheseboundaryconditionsdonotaffect
theaccuracy.
Resultsofthistestcaseontheorthogonalandnon-orthogonalmeshesof100×100cellswith∆t=1sareshowninfigure3forboth
advectionschemes(whichgivesamaximumCourantnumberclosetoone).Thecontoursshowthetracervalueevery100secondsand
thecoloursshowtheerrorssummedevery100seconds.Thedimensionallysplitschemeoutperformsthemulti-dimensionalscheme
onbothmeshesduetothehigher-orderaccuracyofthesplitscheme.Thedimensionallysplitschemeintroducesasmallerrorat300
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
DimensionsplittingforAdvection 9
10000
9000
8000
7000
6000
5000
4000
3000
2000
1000
0
0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000
Figure2.Twonon-orthogonalmeshonaplanewith50×50meshpointsandanglesupto120odesignedtohavesimilaranglestoacubed-spheremesh.
secondswherethetracergoesthroughthechangeindirectionofthemeshwhichwouldbeamelioratedifwewereusingmonotonicity
constraints.Thesecond-order,multi-dimensionalschemeshowsphaselagbuterrorsarealmostentirelyinsensitivethemeshdistortions.
The multi-dimensional and dimensionally split schemes take very different approaches to handling large Courant numbers. The
multi-dimensional scheme uses implicit time-stepping whereas the dimensionally split scheme uses a flux-form semi-Lagrangian
approach,integratingoveralineofcellsinordertocalculatethefluxacrossaface.Implicitschemesareknowntosufferfromphase
errors (eg ??) for long time-steps whereas the accuracy of semi-Lagrangian is less sensitive to time-step (?). Therefore we present
results of both schemes on orthogonal and non-orthogonal meshes for time-steps 10 times those used in figure 3 (∆t =10s) giving
maximumCourantnumbersofaround10using100×100cellsinfigure4.Theerrorofthemulti-dimensionalschemeisagainmuch
largerthanthedimensionallysplitschemeonbothmeshes.ThedimensionallysplitschemeisaccurateatlargeCourantnumbersdespite
thefirst-ordercalculationofdeparturepoints.However,thedimensionallysplitschemeintroducesoscillationsonthenon-orthogonal
mesh,particularlywherethemeshchangesdirectionwhereasforthemulti-dimensionalscheme,errorsarenotstronglyaffectedbythe
non-orthogonality.Figure4clearlyshowsphaseerrorsoftheimplicittime-steppingbut,despitethelargeCourantnumber,thewell
resolvedpartoftheprofileispropagatingatclosetothecorrectspeed.Dispersionanalysis(?)showsthathighfrequencyoscillations
(whicharepoorlyresolvedintimeandspace)willbesloweddramaticallybutfastmovingfeatureswhicharewellresolvedinspace
willpropagateatamuchmorerealisticspeed,supportingtheresultsshowninfigure4.
Inordertocompareconvergencewithresolutionofthedifferentnumericalmethods,weusethe(cid:96) and(cid:96) errornormsdefinedin
2 ∞
theusualway:
(cid:115)ˆ ˆ
(cid:30)
(cid:96)2 = (φ−φT)2dV φT2 dV (27)
V V
(cid:14)
(cid:96)∞ = max|φ−φT| max|φT| (28)
where φT is the analytic solution and the integrations and maxima are over the whole domain, with volume V. Figure 5 shows
convergencewithresolutionofthe(cid:96) and(cid:96) errormeasuresformeshesof50×50,100×100,200×200and400×400cellswithtime-
2 ∞
stepsscaledinordertomaintainamaximumCourantnumberof1(∆t=2, 1, 0.5, 0.25s)orscaledtoachieveamaximumCourant
numberof10(∆t=20, 10, 5, 2.5s).Theerrornormsarecalculatedatt=500seconds,whenthetracerhasmade5/6ofonerevolution
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls
10 Y.Chen,H.Weller,S.PringandJ.Shaw
Figure3.Resultsofthesolidbodyrotationtestcaseonorthogonalandnon-orthogonalmeshesevery100secondsformeshesof100×100cellsand∆t=1sgivinga
maximumCourantnumberofaround1.Thecontoursshowthetracervaluefrom0.1every0.1to0.9andthecoloursshowtheerrors.Thedashedcontourshowstheinitial
conditions.Theminandmaxtextineachplotgivestheminimumandmaximumerrors.
inordertoavoiderrorcancellation.Onbothorthogonalandnon-orthogonalmeshes,themulti-dimensionalschemehassecondorder
convergence once errors are low enough to avoid error saturation (stable errors are bounded at around one). With a large Courant
number, both schemes are less accurate. For the multi-dimensional scheme, this is due to phase errors of the implicit time-stepping
whereasforthemulti-dimensionalschemethefirst-ordererrorsincalculatingthedeparturepointandtrajectorywillbeemergingand
therecouldalsobesignificanterrorsfromthesecond-orderCOSMICsplitting.Ontheorthogonalmesh,thedimensionallysplitscheme
hasthirdorderconvergencefortheCourantnumberclosetooneandsecondorderforthelargerCourantnumber(consistentwiththe
results of Colella and Woodward 1984a; Leonard et al. 1996a). On the non-orthogonal mesh, the dimensionally split scheme has
secondorderconvergesforbothCourantnumbers.
Thedifferenttime-steppingschemesofthetwomodelsalsoaffectaccuracy.The(cid:96) and(cid:96) errormeasuresasafunctionoftime-step
2 ∞
for meshes of 100×100 cells are shown in figure 6. The dimensionally split scheme, which uses flux-form semi-Lagrangian time-
stepping,haserrorsreducingastime-stepincreases,uptoaCourantnumberof2(∆t=2s)whereasthemulti-dimensionalscheme,
(cid:13)c 0000RoyalMeteorologicalSociety
Preparedusingqjrms4.cls