DTIC ADA572213: A Spectral Element Semi-Lagrangian (SESL) Method for the Spherical Shallow Water Equations

30 Pages·0.66 MB·English

by  Defense Technical Information Center

JournalofComputationalPhysics190(2003)623–650
www.elsevier.com/locate/jcp
A spectral element semi-Lagrangian (SESL) method
for the spherical shallow water equations
F.X. Giraldo a,*, J.B. Perot b, P.F. Fischer c
aNavalResearchLaboratory,Monterey,CA93943,USA
bUniversityofMassachusetts,Amherst,MA01003,USA
cArgonneNationalLaboratory,Argonne,IL60439,USA
Received29September2001;receivedinrevisedform23April2003;accepted27May2003
Abstract
Aspectralelementsemi-Lagrangian(SESL)methodfortheshallowwaterequationsonthesphereispresented.The
sphere is discretized using a hexahedral grid although any grid imaginable can be used as long as it is comprised of
quadrilaterals. The equationsare written in Cartesian coordinates to eliminatethe pole singularitywhich plaguesthe
equationsinsphericalcoordinates.Inapreviouspaper[Int.J.Numer.MethodsFluids35(2001)869]weshowedhow
to construct an explicit Eulerian spectral element (SE) model on the sphere; we now extend this work to a semi-La-
grangianformulation.ThenoveltyoftheLagrangianformulationpresentedisthatthehighorderSEbasisfunctions
are used as the interpolation functions for evaluating the values at the Lagrangian departure points. This makes the
method not only high order accurate but quite general and thus applicable to unstructured grids and portable to
distributed memory computers. The equations are discretized fully implicitly in time in order to avoid having to in-
terpolatederivativesatdeparturepoints.ByincorporatingtheCoriolistermsintotheLagrangianderivative,theblock
LUdecompositionoftheequationsresultsinasymmetricpositive-definitepseudo-Helmholtzoperatorwhichwesolve
using the generalized minimum residual method (GMRES) with a fast projection method [Comput. Methods Appl.
Mech. Eng. 163 (1998) 193]. Results for eight test cases are presented to confirm the accuracy and stability of the
method. These results show that SESL yields the same accuracy as an Eulerian spectral element semi-implicit (SESI)
whileallowing for time-steps 10timesas large andbeing upto 70%moreefficient.
(cid:1) 2003Elsevier Science B.V. Allrightsreserved.
Keywords: Cubic gnomonic; Flux-corrected transport; Hexahedral grid; Shallow water equations; Semi-Lagrangian; Semi-implicit;
Spectralelementmethod;Sphericalgeometry
1. Introduction
In [9] a method obtained by fusing the the spectral element method with the semi-Lagrangian method
was introduced and applied to the advection–diffusion equation. A stability analysis revealed that this
*Correspondingauthor.
E-mailaddress:[email protected](F.X.Giraldo).
0021-9991/$-seefrontmatter (cid:1) 2003ElsevierScienceB.V.Allrightsreserved.
doi:10.1016/S0021-9991(03)00300-0
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
APR 2003 2. REPORT TYPE 00-00-2003 to 00-00-2003
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
A spectral element semi-Lagrangian (SESL) method for the spherical
5b. GRANT NUMBER
shallow water equations
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Naval Research Laboratory,Monterey,CA,93943 REPORT NUMBER
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
A spectral element semi-Lagrangian (SESL) method for the shallow water equations on the sphere is
presented. The sphere is discretized using a hexahedral grid although any grid imaginable can be used as
long as it is comprised of quadrilaterals. The equations are written in Cartesian coordinates to eliminate
the pole singularity which plagues the equations in spherical coordinates. In a previous paper [Int. J.
Numer. Methods Fluids 35 (2001) 869] we showed how to construct an explicit Eulerian spectral element
(SE) model on the sphere; we now extend this work to a semi-Lagrangian formulation. The novelty of the
Lagrangian formulation presented is that the high order SE basis functions are used as the interpolation
functions for evaluating the values at the Lagrangian departure points. This makes the method not only
high order accurate but quite general and thus applicable to unstructured grids and portable to distributed
memory computers. The equations are discretized fully implicitly in time in order to avoid having to
interpolate derivatives at departure points. By incorporating the Coriolis terms into the Lagrangian
derivative, the block LU decomposition of the equations results in a symmetric positive-definite
pseudo-Helmholtz operator which we solve using the generalized minimum residual method (GMRES)
with a fast projection method [Comput. Methods Appl. Mech. Eng. 163 (1998) 193]. Results for eight test
cases are presented to confirm the accuracy and stability of the method. These results show that SESL
yields the same accuracy as an Eulerian spectral element semi-implicit (SESI) while allowing for time-steps
10 times as large and being up to 70% more efficient.
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE Same as 28
unclassified unclassified unclassified Report (SAR)
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
624 F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650
hybrid method yields exponential convergence for increasing basis function order while experiencing nei-
therdispersionnordissipationerrors;thatis,providedthatthesolutionissmooth.Inthatwork,however,
we only showed results for the 1D and 2D advection–diffusion equations. In a more recent paper [11], we
showed how to apply the spectral element method to the shallow water equations on the sphere using
Cartesianratherthansphericalcoordinates.Althoughusing Cartesian coordinatesseemscounter-intuitive
this coordinate system avoids the singularities at the poles encountered by spherical coordinates. This
formulation of the equations is not dependent on the grid as is the case with all currently existing models
[15,17,20,32]. In that paper [11] we showed that this approach yields exponential convergence for the
shallow water equations; however, in that work the time integration scheme used was an explicit Eulerian
method (third order Adams–Bashforth).
Although explicit time-integration methods for atmospheric simulations are quite efficient, their main
disadvantage is that small time-steps must be observed in order to maintain stability. The reason for this
prohibitivelysmalltime-stepisduetothefastmovinggravitywaves.Thesewavesrequireasmalltime-step
whileonlycarryingaverysmallpercentofthetotalenergyinthesystem.Inordertoamelioratethisrather
stringenttime-steprestriction researchershavetriedvarious approachessuchasusingalargerdifferencing
stencilforthegravitywavetermstherebyeffectivelyreducingtheCourantnumber[24,35],anddiscretizing
thegravitywavetermsimplicitlyintime[19].Thefirststrategy,thatisusingalargerdifferencingstencilfor
the gravity wave terms, is typically done to avoid the inf–sup (Babuska–Brezzi) condition [18] associated
with the Stokes problem. However, discretizing the gravity wave terms implicitly in time is a much more
effective way of increasing the time-step; this is the approach we follow in this paper.
After the gravity wave terms have been successfully discretized the next set of terms responsible for
controlling the maximum time-step are the advection terms. In order to use increasingly larger time-steps
researchers have turned to Lagrangian methods for treating these recalcitrant terms [26,33]. By rewriting
theequationsintermsoftheLagrangianderivativethetroublesomeadvectiontermsareabsorbedintothe
Lagrangian derivative. Thus the equations in this form are now discretized in time along characteristics
which results in a much more stable numerical method due to the virtual disappearance of the Courant–
Friedrichs–Lewy (CFL) condition. We say virtual because any time-step is possible provided that the
Lipschitz condition is not violated and that the departure points of the semi-Lagrangian method are
computed accurately.
In the current paper we combine a semi-Lagrangian method with high order basis functions and
extend this hybrid method to the solution of the spherical shallow water equations. Thus we are es-
sentially extracting the highlights of our previous papers [9,11] and fusing them into a powerful technique
for solving the shallow water equations on the sphere. The allure of this method is that it achieves the
same order of accuracy obtained with exponentially high order Eulerian methods [11] while using much
larger time-steps. It should be mentioned that this is not the only possible characteristic method. There
are other variants that are certainly worth considering especially the characteristic method of Maday et
al. [22] which we hope to compare with SESL in future work. In this paper we only consider SESL
because it represents the only true characteristic high order method which fuses both space and time into
a single formulation.
The remainder of this paper is organized as follows. Section 2 describes the governing equations of
motion used to test our numerical method. In Section 3 we describe the discretization of the governing
equations.Thisincludesthespatialdiscretizationbythespectralelementmethod,thetimediscretizationby
thesemi-Lagrangianmethod,andthesolutionoftheresultingimplicitequations.InSection4wedescribe
thetessellationofthesphereintoquadrilateralelementswhichareusedtoconstructthediscreteoperators.
Finally,inSection5wepresentconvergenceratesofthespectralelementsemi-Lagrangian(SESL)method
and compare it with an Eulerian spectral element semi-implicit (SESI) method. This then leads to some
conclusions about the feasibility of this approach for constructing future atmospheric models and a dis-
cussion on the direction of future work.
F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650 625
2. Shallow water equations
Theshallowwaterequationsareasystemoffirstordernonlinearhyperbolicequationswhichgovernthe
motionofaninviscidincompressiblefluidinashallowdepth.Thepredominantfeatureofthistypeoffluid
isthatthecharacteristiclengthofthefluidisfargreaterthanitsdepthwhichisanalogoustothemotionof
air in the atmosphere. For this reason these equations are typically used as a first step toward the con-
structionofeithernumericalweatherprediction(NWP)orclimatemodels;theformertypeofmodelisthe
final goal of our research.
The spherical shallow water equations in Cartesian advection form are
ou
þu(cid:1)$uu¼(cid:2)u$(cid:1)u; ð1Þ
ot
ou 2xz
þu(cid:1)$u¼(cid:2) ðx(cid:3)uÞ(cid:2)$ðuþusÞ(cid:2)lx; ð2Þ
ot a2
whereuisthegeopotentialheight(u¼gh,wheregisthegravitationalconstantandhistheverticalheight
ofthefluid),us isthesurfacetopography(e.g.,mountains),uistheCartesianwindvelocityvectorðu;v;wÞ,
and ðx;aÞ represent the rotation of the earth and its radius, respectively.
Thetermlx,wherex¼ðx;y;zÞisthepositionvectorofthegridpoints,isafictitiousforceintroducedto
constrain the fluid particles to remain on the surface of the sphere. By switching from spherical (2D) to
Cartesian (3D) coordinates we have allowed the fluid particles an additional degree of freedom which will
manifest itself in the fluid particles flying off the sphere. In order to prevent this from happening we in-
troduce the Lagrange multiplier l (see [4]); we address the calculation of this coefficient in Section 3.4.
The shallow water equations in Cartesian form have received significant attention recently (see
[10–12,16,31]). It should be mentioned that using these equations in Cartesian form on the sphere intro-
ducesnoapproximationswhatsoever;theequationsarecompletelyequivalenttotheequationsinspherical
coordinates as shown by Swarztrauber et al. [31]. Although the formulation we use in our paper may
appear slightly different it can be shown to be identical to that of [16,31].
3. Discretization
Inthis section wedescribethe discretization oftheshallowwaterequations.InSection3.1, wedescribe
the spatial discretization by the spectral element method including: the choice of basis functions and the
mappingfromphysical3DCartesian spacetothelocal2Dcomputationalspace.InSection3.2wediscuss
theimplicittimediscretizationbythesemi-Lagrangianmethodincludingthediscretizationofthetrajectory
equation.Inaddition, we describetheinterpolationoperatorsrequiredby thesemi-Lagrangianmethod as
wellastheconstructionofmonotoneinterpolationschemes.Finally,inSection3.4wediscussthesolution
of the coupled implicit equations.
3.1. Spatial discretization by the spectral element method
3.1.1. Basis functions and integration
To define the local operators which shall be used to construct the global approximation of the solution
webeginbydecomposingthesphericaldomainXintoN non-overlappingquadrilateralelementssuchthat
e
[Ne
X¼ X :
e
e¼1
626 F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650
Toperformdifferentiationandintegrationoperations,weintroducethenon-singularmappingx¼WðnÞ
whichdefinesatransformationfromthephysicalCartesiancoordinatesystemx¼ðx;y;zÞdefinedinX to
e
the reference coordinate system n¼ðn;g;fÞ defined in each element such that ðn;gÞ lies on the spherical
surface; where ðn;gÞ2½(cid:2)1;þ1(cid:4)2 in each element.
Associated with the local mapping, W, is the transformation Jacobian, J ¼ox=on, and the determinant
ox ox ox
jJj¼ (cid:1)G; G ¼ (cid:3) ;
of on og
where G represents the surface conforming component of the mapping (see [11] for further details).
We can now use this mapping to define the local representation of the solution, q¼ðu;uÞ, and the
approximationofoperationssuchasdifferentiationandintegration.Forsimplicity,weassumeftobeunity
in what remains and denote n¼ðn;gÞ.
Thesimplestructureofthereferenceelement,I,spannedbyn2½(cid:2)1;1(cid:4)2,makesitnaturaltorepresentthe
local element-wise solution q by an Nth order polynomial in n as
ðNþ1Þ2
X
q ðxÞ¼ w ðxÞq ðx Þ;
N k N k
k¼1
where x represents ðN þ1Þ2 grid points and w ðxÞ reflects the associated multivariate Lagrange polyno-
k k
mial.ThelogicalsquarestructureofIsimplifiesmattersinthatwecanexpresstheLagrangepolynomialby
a tensor-product as
w ðxÞ¼hðnðxÞÞhðgðxÞÞ;
k i j
where i;j¼0;...;N.
For the grid points ðn;gÞ we choose the Legendre–Gauss–Lobatto (LGL) points, given as the tensor-
i j
product of the roots of
ð1(cid:2)n2ÞP0ðnÞ¼0;
N
whereP ðnÞistheNthorderLegendrepolynomial.Thischoicesimplifiestheconstructionofthealgorithm
N
because the LGL points are also used as the sampling points in the Gaussian quadrature rule required by
the numerical integration which we shall describe shortly. The one-dimensional Lagrange polynomials,
hðnÞ, are
i
1 ð1(cid:2)n2ÞP0ðnÞ
hðnÞ¼(cid:2) N ;
i NðN þ1Þ ðn(cid:2)nÞP ðnÞ
i N i
and likewise for hðgÞ.
j
ThechoiceoftheLGLpointsenablesthestraightforwardapproximationoflocalelementintegrals,i.e.,
Z Z N
X
qðxÞdx¼ qðnÞJðnÞdn’ xðnÞxðgÞqðn;gÞJðn;gÞ;
i j i j i j
Xe I i;j¼0
whereJ representsthelocalJacobianforthetransformationbetweenX andI,andxðnÞandxðgÞarethe
e i j
Gaussian quadrature weights,
2 (cid:1) 1 (cid:2)2
xðnÞ¼ ;
i NðN þ1Þ P ðnÞ
N i
associated with the one-dimensional LGL quadrature.
F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650 627
To simplify the description of the numerical algorithm, let us define the following local element oper-
ators: let
Z
Me ¼ wðxÞwðxÞdx
ij i j
Xe
represent the mass matrix and
Z
De ¼ wðxÞ$wðxÞdx
ij i j
Xe
the differentiation matrix where i;j¼1;...;ðN þ1Þ2. Note that D¼ðDx;Dy;DzÞ is a vector of matrices
corresponding to the three spatial directions for D.
3.1.2. Satisfying the equations globally
To satisfy the equations globally requires assembling the global solution by virtue of an element-wise
construction. This element-wise construction is based on the summation of the local element matrices
to form their global representation. This summation procedure is known as the global assembly
or direct stiffness summation. Let us represent this global assembly procedure by the summation
operator
^Ne
e¼1
with the mapping ði;eÞ!ðIÞ where i¼1;...;ðN þ1Þ2 are the local element grid points, e¼1;...;N are
e
thespectralelementscoveringthesphericalshell,andI ¼1;...;N aretheglobalgridpoints.Applyingthe
p
global assembly operator to the local element matrices results in the following global matrices:
^Ne
M ¼ Me
e¼1
for the mass matrix and
^Ne
D¼ De
e¼1
for the differentiation matrix.
With these operators defined and by denoting the global grid vector for the grid points as x , the
G
geopotentialasu ,andthewindvelocity asu wecannowwritethesemi-discreteapproximation toEqs.
G G
(1) and (2) as follows:
ou
M GþuTDu ¼(cid:2)u DTu ; ð3Þ
ot G G G G
(cid:1) (cid:2)
MouGþuTDu ¼(cid:2)M 2xzGðx (cid:3)u Þ (cid:2)D(cid:3)u þus (cid:4)(cid:2)Mðlx Þ: ð4Þ
ot G G a2 G G G G G
Itshouldbenotedthatthemassmatrix,M,isdiagonalandtherebytrivialtoinvert.Thediagonalproperty
ofthismatrixisduetothedualroleoftheLGLpointswhichareusedbothasgridpointsandquadrature
points.
628 F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650
3.2. Time discretization by the semi-Lagrangian method
Thetime-steprestrictionoftheshallowwaterequationsisdominatedbytheadvectionandgravitywave
p
terms.Infact,thecharacteristicwavespeedofthissystemisC ¼U þ ffiuffiffiffi.However,thegravitywaveterm
although traveling at far greater speeds than the wind velocity only carries a fraction of the total energy.
Thus, explicit time integration schemes must use very small time-steps in order to maintain stability while
muchlargervaluescouldbeusedwithoutdiminishingaccuracy.Typically,inordertoenlargethetime-step,
semi-implicit schemes [19] have been used whereby the advection and Coriolis terms are discretized ex-
plicitly in time while the remaining terms are discretized implicitly. The generalized leapfrog method
unþ1(cid:2)un(cid:2)1 ¼(cid:2)ðu(cid:1)$uÞn(cid:2)U(cid:3)h$(cid:1)unþ1þð1(cid:2)hÞ$(cid:1)un(cid:2)1(cid:4); ð5Þ
2Dt
unþ1(cid:2)un(cid:2)1 2xz (cid:6) (cid:7)
¼(cid:2)ðu(cid:1)$uÞn(cid:2) ðx(cid:3)uÞn(cid:2) h$ðuþusÞnþ1þð1(cid:2)hÞ$ðuþusÞn(cid:2)1 ð6Þ
2Dt a2
has been the most popular method for constructing semi-implicit algorithms for the shallow water and
atmospheric equations (see for example the operational models in [14,17,29,30] and the shallow water
modelin[34]).NotethatuhasbeenlinearizedaboutsomemeanstateUinordertoavoidhavingtosolvea
nonlinearproblemandwehaveomittedtheLagrangemultiplierforthetimebeing.InEqs.(5)and(6)the
parameterhdeterminesthetypeofschemeusedinconjunctionwiththeleapfrogmethodfortheadvection
terms. For example, h¼0;1;1 yields the Euler method (explicit), the Crank–Nicholson method (implicit),
2
and the backward Euler method (implicit), respectively. For h¼1 the scheme is second order accurate in
2
time and for all other values it is only first order accurate.
Tofurtherenlargethetime-stepwecanincorporatetheadvectiontermsintotheLagrangianderivative;
this is the semi-Lagrangian method [26]. Integrating Eqs. (5) and (6) along characteristics results in
unþ1(cid:2)uuen(cid:2)1 ¼(cid:2)U(cid:3)h$(cid:1)unþ1þð1(cid:2)hÞ$(cid:1)uun(cid:2)1(cid:4);
2Dt e
unþ1(cid:2)uun(cid:2)1 2xzz (cid:6) (cid:7)
e ¼(cid:2) eðxx(cid:3)uuÞn(cid:2) h$ðuþusÞnþ1þð1(cid:2)hÞ$ðuuþuusÞn(cid:2)1 ;
2Dt a2 e e e e
where it should be understood that the values at times n and n(cid:2)1 are not the values at the Eulerian grid
points but rather those at the Lagrangian departure points found by integrating the trajectory equation
dx
¼u; ð7Þ
dt
where the Lagrangian derivative is
d o
¼ þu(cid:1)$:
dt ot
We have marked the variables that need to be interpolated at the departure points with a tilde (~). The
statement of the semi-Lagrangian method is: find qqðn;n(cid:2)1Þ using Eq. (7) such that at time nþ1 these
e
characteristics arrive at the Eulerian grid points (arrival points).
Although the Crank–Nicholson method has been the most popular method used with the semi-La-
grangianmethoditrequiresthecomputationofnotonlytheprimitivevariablesatthedeparturepointsbut
their derivatives as well. This step in the semi-Lagrangian procedure requires a very careful treatment;
solvingthesetermsindiscriminatelyresultsinseverelossofconservation.Forthisreason,manyoperational
F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650 629
NWP models use off-centering, h>1=2, to combat the inaccuracies of computing derivatives at the de-
parturepoints.Theproblemwithoff-centeringisthattheschemeisonlyfirstorderaccurate.Modelswhich
perform off-centering in either an Eulerian or semi-Lagrangian context include the operational NWP and
climate models in [14,17,29,30].
Galerkin methods pose some difficulties for semi-Lagrangian methods. The problem stems from the C0
basisfunctionstypicallyusedinGalerkinmethodswhichtranslatesintothederivativesbeingdiscontinuous
acrosselements.Forthisreasonmanysemi-Lagrangianimplementationsconstructtheright-handside(i.e.,
thetermsthatmustbeevaluatedatthedeparturepoints)atthearrivalpointsandtheninterpolatetheentire
right-hand side vector to the departure points. This, however, will conserve neither mass nor energy. This
leads us to the following conjecture.
Conjecture 1. Semi-Lagrangian methods with C0 basis functions should only be used for homogeneous
Lagrangian equations of the type
du
¼0:
dt
However, if there are source terms with derivatives on the right-hand side then they should be either
evaluated at arrival points only or Ck basis functions should be used, where k represents the maximum
orderofthederivativesinthesystem (k ¼1fortheshallow waterequations).Wereservetheproofofthis
conjecture for a future study.
The problem with using Ck functions is that the number of unknowns increases from 4 to 16 per grid
point due to the first order derivatives that need to be carried. Xiu and Karniadakis [37] also noted some
problemswithevaluatingderivativesatthedeparturepointsandthusdecidedtouseabackwarddifference
formula (BDF) which gives rise to a fully implicit scheme.
Beforediscretizingtheequationsalongthecharacteristicsletuswritetheshallowwaterequationsinthe
Lagrangian form which we shall use to solve the equations. Using the idea of Rochas [27] we include the
Coriolis term inside the Lagrangian derivative and write the equations as follows:
du
¼(cid:2)U$(cid:1)u; ð8Þ
dt
d
ðuþ2xðk(cid:3)xÞÞ¼(cid:2)$ðuþusÞ(cid:2)lx; ð9Þ
dt
dx
¼u;
dt
wherek¼ð0;0;1Þistheunitvectorinthez-direction.Wecannowdiscretizetheseequationsusingasecond
order BDF to yield
unþ1(cid:2)4uunþ1uun(cid:2)1 2
3e 3e ¼(cid:2) U$(cid:1)unþ1; ð10Þ
Dt 3
uunþ1(cid:2)4uunþ1uun(cid:2)1 2
b 3b 3b ¼(cid:2) $ðuþusÞnþ1; ð11Þ
Dt 3
where
uu ¼uþ2xðk(cid:3)xÞ
b
630 F.X.Giraldoetal./JournalofComputationalPhysics190(2003)623–650
represents the absolute velocity in an inertial reference frame. Note that Eqs. (10) and (11) do not require
the interpolation of derivatives at the departure points; for this reason we have elected to use a BDF.
Moving the implicit terms to the left-hand side and the explicit terms to the right-hand side results
in
2
unþ1þ DtU$(cid:1)unþ1 ¼bu; ð12Þ
3
2
unþ1þ Dt$unþ1 ¼bu; ð13Þ
3
where the right-hand side vectors are defined as follows:
4 1
bu ¼ uun(cid:2) uun(cid:2)1; ð14Þ
3e 3e
4 1 2
bu ¼ ðuuþ2xðk(cid:3)xxÞÞn(cid:2) ðuuþ2xðk(cid:3)xxÞÞn(cid:2)1(cid:2) Dtð$usÞnþ1(cid:2)2xðk(cid:3)xÞnþ1: ð15Þ
3 e e 3 e e 3
We now describe the computation of the departure point values.
3.2.1. Trajectory equation
Inordertoevaluatethevariables qqðn;n(cid:2)1Þ requiretheaccuratecomputationofthedeparturepointswhich
e
are found from Eq. (7). The simplest second order scheme for discretizing this equation is the two-step
Runge–Kutta scheme:
Dt
xxnþ1=2 ¼xnþ1(cid:2) uðxnþ1;tnþ1Þ;
e 2
xxn ¼xnþ1(cid:2)Dtu(cid:3)xxnþ1=2;tnþ1=2(cid:4);
e e
Dt
xxn(cid:2)1=2 ¼xn(cid:2) uðxxn;tnÞ;
e 2 e
xxn(cid:2)1 ¼xn(cid:2)Dtu(cid:3)xxn(cid:2)1=2;tn(cid:2)1=2(cid:4);
e e
where the velocity field is extrapolated as follows:
uðx;tnþ1Þ¼2uðx;tnÞ(cid:2)uðx;tn(cid:2)1Þ;
u(cid:3)x;tnþ1=2(cid:4)¼3uðx;tnÞ(cid:2)1uðx;tn(cid:2)1Þ;
2 2
u(cid:3)x;tn(cid:2)1=2(cid:4)¼1uðx;tnÞþ1uðx;tn(cid:2)1Þ:
2 2
However, we do not have to stop at second order but instead can construct higher order trajectory ap-
proximations. For example we can construct the fourth order Runge–Kutta approximation
xxn ¼xnþ1(cid:2)Dt(cid:8)u(cid:3)xxð1Þ;tnþ1(cid:4)þ2u(cid:3)xxð2Þ;tnþ1=2(cid:4)þ2u(cid:3)xxð3Þ;tnþ1=2(cid:4)þu(cid:3)xxð4Þ;tn(cid:4)(cid:9); ð16Þ
e 6 e e e e