Energy and Society
Energy and Society PHYSICS 190
Popular in Course
Popular in Physics 2
This 28 page Class Notes was uploaded by Yvonne McKenzie on Friday October 30, 2015. The Class Notes belongs to PHYSICS 190 at University of Massachusetts taught by Staff in Fall. Since its upload, it has received 11 views. For similar materials see /class/232315/physics-190-university-of-massachusetts in Physics 2 at University of Massachusetts.
Reviews for Energy and Society
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/30/15
Available online at www5ciencedirectcom SCIENCEDIRECT COMPUTATIONAL PHYSICS EL EVIER Journal of Computational Physics 190 2003 623 650 wwwelseviercomlocatejcp A spectral element semiLagrangian SESL method for the spherical shallow water equations FX Giraldo 3 JB Perot b PF Fischer C a Naval Research Laboratory Monterey CA 93943 USA b University of Massachusetts Amherst MA 01003 USA C Argonne National Laboratory Argonne IL 60439 USA Received 29 September 2001 received in revised form 23 April 2003 accepted 27 May 2003 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 La grangian 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 in terpolate 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 de nite 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 con rm 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 el cient 2003 Elsevier Science BV All rights reserved Keywords Cubic gnomonic Flux corrected transport Hexahedral grid Shallow water equations Semi Lagrangian Semi implicit Spectral element method Spherical geometry 1 Introduction In 9 a method obtained by fusing the the spectral element method with the semiLagrangian method was introduced and applied to the advection diffusion equation A stability analysis revealed that this Corresponding author E mail address giraldonrlmrynavymil FX Giraldo 0021 9991 see front matter 2003 Elsevier Science BV All rights reserved doi101016S0021 99910300300 0 624 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 hybrid method yields exponential convergence for increasing basis function order while experiencing nei ther dispersion nor dissipation errors that is provided that the solution is smooth In that work however we only showed results for the ID 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 Cartesian rather than spherical coordinates Although using Cartesian coordinates seems counterintuitive 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 15172032 In that paper II 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 timeintegration methods for atmospheric simulations are quite ef cient their main disadvantage is that small timesteps must be observed in order to maintain stability The reason for this prohibitively small timestep is due to the fast moving gravity waves These waves require a small timestep while only carrying a very small percent of the total energy in the system In order to ameliorate this rather stringent timestep restriction researchers have tried various approaches such as using a larger differencing stencil for the gravity wave terms thereby effectively reducing the Courant number 2435 and discretizing the gravity wave terms implicitly in time 19 The first strategy that is using a larger differencing stencil for 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 timestep 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 timestep are the advection terms In order to use increasingly larger timesteps researchers have turned to Lagrangian methods for treating these recalcitrant terms 2633 By rewriting the equations in terms of the Lagrangian derivative the troublesome advection terms are absorbed into the 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 timestep is possible provided that the Lipschitz condition is not violated and that the departure points of the semiLagrangian method are computed accurately In the current paper we combine a semiLagrangian 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 911 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 I I while using much larger timesteps 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 This includes the spatial discretization by the spectral element method the time discretization by the semiLagrangian method and the solution of the resulting implicit equations In Section 4 we describe the tessellation of the sphere into quadrilateral elements which are used to construct the discrete operators Finally in Section 5 we present convergence rates of the spectral element semiLagrangian SESL method and compare it with an Eulerian spectral element semiimplicit SESD 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 EX Giralda et al Jaumal of Computational Physics 190 2003 623 650 625 2 Shallow Water equations The shallow water equations are a system of first order nonlinear hyperbolic equations which govern the motion of an inviscid incompressible uid in a shallow depth The predominant feature of this type of uid is that the characteristic length of the uid is far greater than its depth which is analogous to the motion of air in the atmosphere For this reason these equations are typically used as a first step toward the con struction of either numerical weather prediction NWP or climate models the former type of model is the nal goal of our research The spherical shallow water equations in Cartesian advection form are a g0uVgou chu I at an 2032 s auVa 7xgtltu Vcpcp ix 2 where q is the geopotential height 42 gh where g is the gravitational constant and h is the vertical height of the uid of is the surface topography eigi mountains u is the Cartesian wind velocity vector u 1 w and co a represent the rotation of the earth and its radius respectively The term ax where x x yz is the position vector of the grid points is a fictitious force introduced to constrain the uid particles to remain on the surface of the sphere By switching from spherical 2D to Cartesian 3D coordinates we have allowed the uid particles an additional degree of freedom which will manifest itself in the uid particles ying off the sphere In order to prevent this from happening we in troduce the Lagrange multiplier a see 4 we address the calculation of this coe icient in Section 34 The shallow water equations in Cartesian form have received signi cant attention recently see 10 12163l It should be mentioned that using these equations in Cartesian form on the sphere intro duces no approximations whatsoever the equations are completely equivalent to the equations in spherical coordinates as shown by Swarztrauber et a1 31 Although the formulation we use in our paper may appear slightly different it can be shown to be identical to that of 1631 3 Discretiza on In this section we describe the discretization of the shallow water equations In Section 31 we describe the spatial discretization by the spectral element method including the choice of basis functions and the mapping from physical 3D Cartesian space to the local 2D computational space In Section 32 we discuss the implicit time discretization by the semiLagrangian method including the discretization of the trajectory equation In addition we describe the interpolation operators required by the semiLagrangian method as well as the construction of monotone interpolation schemes Finally in Section 34 we discuss the solution of the coupled implicit equations 3 Spatial discretization by the spectral element method 31 Basis functions and integration To define the local operators which shall be used to construct the global approximation of the solution we begin by decomposing the spherical domain 2 into Ne nonoverlapping quadrilateral elements such that Na Qert 21 626 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 To perform differentiation and integration operations we introduce the nonsingular mapping x TE which de nes a transformation from the physical Cartesian coordinate system x x yz de ned in 22 to the reference coordinate system E Em E de ned in each element such that Em lies on the spherical surface where E n E 1 12 in each element Associated with the local mapping T is the transformation Jacobian J fix6E and the determinant 6x 6x 6x J G G X l l 5C 55 57 where G represents the surface conforming component of the mapping see 11 for further details We can now use this mapping to de ne the local representation of the solution 11 4214 and the approximation of operations such as differentiation and integration For simplicity we assume E to be unity in what remains and denote E E n i The simple structure of the reference element spanned by E E 1 llz makes it natural to represent the local elementwise solution 11 by an Nth order polynomial in E as N1 2 lt tinX Z WkxqNxk where xk represents N 12 grid points and 11x re ects the associated multivariate Lagrange polyno mial The logical square structure of simpli es matters in that we can express the Lagrange polynomial by a tensorproduct as 1PM hx xhnx where ij 0H Ni For the grid points E m we choose the Legendre Gauss Lobatto LGL points given as the tensor product of the roots 0 1 ZP1 v 0 where PNE is the Nth order Legendre polynomial This choice simpli es the construction of the algorithm 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 onedimensional Lagrange polynomials hE are W 1 1 5 PM NN1m and likewise for h The choice of the LGL points enables the straightforward approximation of local element integrals ie 9 gm dx mm d Zwltcgtwltmgtqlt woman whereJ represents the local Jacobian for the transformation between 22 and and cuE and 0300 are the Gaussian quadrature weights ma associated with the onedimensional LGL quadrature EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 627 To simplify the description of the numerical algorithm let us de ne the following local element oper ators et M ixxwx dx represent the mass matrix and D wxwwxx dx the differentiation matrix where ij 1V N 1 Note that D UDVU is a vector of matrices corresponding to the three spatial directions for D 312 Satisfying the equations globally To satisfy the equations globally requires assembling the global solution by virtue of an elementwise construction This elementwise 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 sti ness summation Let us represent this global assembly procedure by the summation operator 3 with the mapping 12 gt I where i l t t t N l2 are the local element grid points z l t t t Ne are the spectral elements covering the spherical shell and I l t t t Np are the global grid points Applying the global assembly operator to the local element matrices results in the following global matrices Na M1Me for the mass matrix and for the differentiation matrix With these operators defined and by denoting the global grid vector for the grid points as xg the geopotential as ch and the wind velocity as m we can now write the semidiscrete approximation to Eqs 1 and 2 as follows a M quot D l G CPGDTuG7 3 OMG I 2cozG s M uGDuG M a2 xG gtlt uG Dch 42G Mprt 4 It should be noted that the mass matrix M is diagonal and thereby trivial to invert The diagonal property of this matrix is due to the dual role of the LGL points which are used both as grid points and quadrature points 628 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 32 Time discretization by the semiLagrangian method The timestep restriction of the shallow water equations is dominated by the advection and gravity wave terms In fact the characteristic wave speed of this system is C U p However the gravity wave term 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 timesteps in order to maintain stability while much larger values could be used without diminishing accuracy Typically in order to enlarge the timestep semiimplicit 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 n1 pan 2A uVlt0quotltP9V39quot 19Vquot 1gt 5 uquot12AtuLl ltu Wquot ltx X nquot WW Wm lt1 6N NH 6 has been the most popular method for constructing semiimplicit algorithms for the shallow water and atmospheric equations see for example the operational models in l4l72930 and the shallow water model in 34 Note that g has been linearized about some mean state D in order to avoid having to solve a nonlinear problem and we have omitted the Lagrange multiplier for the time being In Eqs 5 and 6 the parameter 6 determines the type of scheme used in conjunction with the leapfrog method for the advection terms For example 6 0 1 yields the Euler method explicit the Crank Nicholson method implicit and the backward Euler method implicit respectively For 6 the scheme is second order accurate in time and for all other values it is only first order accurate To further enlarge the timestep we can incorporate the advection terms into the Lagrangian derivative this is the semiLagrangian method 26 Integrating Eqs 5 and 6 along characteristics results in nlanrl T lt1gt6Vu 11 6V quot 1 quotHI W71 ZWEN 71 5nl N Ngnil T7xxu9Vltltpltp lt1 6gtvltltpltpgt where it should be understood that the values at times n and n l are not the values at the Eulerian grid points but rather those at the Lagrangian departure points found by integrating the trajectory equation dx 7 d u where the Lagrangian derivative is d 6 a a ll Vi We have marked the variables that need to be interpolated at the departure points with a tilde i The statement of the semiLagrangian method is find 39quot using Eq 7 such that at time n1 these characteristics arrive at the Eulerian grid points arrival points Although the Crank Nicholson method has been the most popular method used with the semiLa grangian method it requires the computation of not only the primitive variables at the departure points but their derivatives as well This step in the semiLagrangian procedure requires a very careful treatment solving these terms indiscriminately results in severe loss of conservation For this reason many operational EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 629 NWP models use offcentering 6 gt 12 to combat the inaccuracies of computing derivatives at the de parture points The problem with offcentering is that the scheme is only first order accurate Models which perform offcentering in either an Eulerian or semiLagrangian context include the operational NWP and climate models in l4l72930 Galerkin methods pose some dif culties for semiLagrangian methods The problem stems from the CO basis functions typically used in Galerkin methods which translates into the derivatives being discontinuous across elements For this reason many semiLagrangian implementations construct the righthand side ie the terms that must be evaluated at the departure points at the arrival points and then interpolate the entire righthand side vector to the departure points This however will conserve neither mass nor energy This leads us to the following conjecture Conjecture 1 SemiLagrangian methods with CO basis functions should only be used for homogeneous Lagrangian equations of the type 1 4 0 dt However if there are source terms with derivatives on the righthand side then they should be either evaluated at arrival points only or CE basis functions should be used where k represents the maximum order of the derivatives in the system k l for the shallow water equations We reserve the proof of this conjecture for a future study The problem with using CIt 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 problems with evaluating derivatives at the departure points and thus decided to use a backward difference formula BDF which gives rise to a fully implicit scheme Before discretizing the equations along the characteristics let us write the shallow water equations in the 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 dw E lt1 Vu 8 d S 50 2wltk x x vltqz a gt x lt9 E d 7quot where k 0 0 l is the unit vector in the z direction We can now discretize these equations using a second order BDF to yield nlananrl 2 lt1iVu 1 10 At 3 lt nl nl nrl 2 V snl 11 A 3 gm 7 where iiu2cukgtltx 630 EX Giralda et a Jaumal 0f Camputatianal Physics 190 2003 623 650 represents the absolute velocity in an inertial reference frame Note that Eqs 10 and l 1 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 lefthand side and the explicit terms to the righthand side results in 2 lt0quot1 Atlt1 Vuquot1 m 12 1 2 n1 u uquot Ath2 b7 13 where the righthand side vectors are defined as follows 4 l b Nquot Nquot 1 14 3 p 3 p 1174 N 71 1 n71 2 571l n1 b 75u2mkgtltx u2cukgtlt x Atch 2mkgtltx t 15 We now describe the computation of the departure point values 32 Trajectory equation In order to evaluate the variables 39quot require the accurate computation of the departure points which are found from Eq 7 The simplest second order scheme for discretizing this equation is the twostep Runge Kutta scheme N At xnlZ xnl 7quot xHerH7 in xnl Atu nl27trxlZ7 N At N xanZ xn 714x717 33an xn Atu x4nrl27trxrlZ7 where the velocity field is extrapolated as follows uxtquot1 2uxtquot uxtquot 17 ux KM 2ux t lux t H 7 2 7 2 7 7 1 1 ux7 tquot 12 514x t 514x tquot 1 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 At 6 xv n1 Waxy 2uEZtquot1Z 2ult3gttn1Z quotimam 16 EX Giralda et al Jaumal of Computational Physics 190 2003 623 650 631 where the trial points are EU 17 At Ea 117ult x17tn17 N At N x3 11 7quotformation27 334 11 Atu nl27tnlZ and the velocity eld needs to be extrapolated to suf ciently high order We have included a tilde in the trial points EU for k l t t t 4 to remind the reader that the value of the velocity field 14 must be inter polated at these trial points since they will generally not lie on a grid point 322 Interpolation and uxcorrected transport Once we determine the departure point and the element in which it lies we must then construct an in terpolation operator using its local information the departure point search is described in Section 42 Let us assume that the element 2d claims the departure point xd Then the interpolated value E2 is simply given as ad Z th dh17ld w x010 where 42 are the grid point values of the element 5 and the departure point value xd has been mapped to computational space id nd Thus for high order N we cannot expect to have a monotonic interpolation Therefore just as the spectral element method is not naturally monotonic than neither is the high order semiLagrangian method The lack of monotonicity does not have serious repercussions except when very large gradients such as shocks are encountered One of the test cases that we use in this paper contains a shock Case 1C and for this case we shall use a monotonic version of SESL The simplest way to make a high order interpolation monotonic is to use a linear interpolation operator In other words instead of using all of the N 1 grid points in an element we can just use its vertices However using these values indiscriminately throughout the computation will result in a very diffused solution Therefore we apply the monotonicity constraint based on whether the full polynomial interpo lation generates values that are larger or smaller than the values of the element 2d This FCTtype method uxcorrected transport see 1 can be written as follows ad minlcpmymaxwm 10 where N is the interpolation using the full polynomial N N am Z th d xndmp x o J o and pmX and 42mm are the maximum and minimum values in the element 5 respectively In other words this method now keeps the interpolated value from exceeding the extrema of the element 2d and thus resulting in a monotonic interpolation operator 632 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 33 Space time discretization In all previous uid dynamics applications whereby Lagrangian schemes are combined with Galerkin methods either spectral or finite elements the two methods are applied virtually independentlyi As an example when spectral methods are used only adhoc cubic Lagrange polynomials are constructed to perform the interpolation at the departure points 25 In the finite element method only low order finite element methods linear are used in conjunction with cubic interpolation 5 In these approaches the interpolation operator is constructed independently of the spatial discretization methods note that the two previously mentioned works both use cubic interpolation even though they use different order accuracy spatial discretization methods In our approach we use the high order basis functions of the spectral el ement method as the interpolation functions and so the interpolation operators are naturally and inherently contained within the spectral element method Since the appearance of 9 in which we introduced SESL for the advection diffusion equation two papers have appeared which have contributed to the further development of SESLtype methods In the paper by Falcone and Ferretti 7 they show a rigorous proof for the approximation error of the advection equation with Lipschitz continuous solutions In that paper the approximation error for Lagrangian methods using high order interpolation was shown to AxN1 K 0ltAt At 17 where K is the order of accuracy used to solve the trajectory equation dx m quot andN is the order of the interpolation functionsi Therefore if either high order K is used with low order N or viceversa then the overall accuracy of the scheme is compromised In all previous semiLagrangian im plementations researchers used both low order K and N i They then noticed that when they increased one or the other they did not obtain increased accuracy This led to the incorrect conclusion that it is not bene cial to increase these parameters beyond a certain value however no one tried to increase both values In the paper by Xiu and Karniadakis 37 they show numerically that the error of the semiLagrangian method combined with the spectral element method is indeed controlled by Eq 17 they show this for the 2D advection diffusion equation and argue that the same behavior holds for the Navier Stokes equations at high Reynolds number In this paper we show that the error in SESL is bounded by Eq 17 Besides giving a bound on the error for high order Lagrangian methods the result of Falcone and Ferretti can also be used to obtain the optimal timestep to use with a given order of accuracy spatial discretization The optimal timestep is obtained when the first and second terms in Eq ID are balanced which occurs for the following timestep 1 AtK A i 18 Following Malevsky and Thomas 21 we can write the timestep as a function of the grid spacing as follows At Ax 19 Substituting Eq 19 into Eq 18 results in the optimal timestep a 7 N l TKr EX Giralda et a Jaumal 0f Camputatianal Physics 190 2003 623 650 633 34 Implicit formulation of the discrete system In this section we describe the implicit formulation for the discrete equations Discretizing Eqs 12 and 13 by the spectral element method results in the following discrete system M422 MDTu39gI B 20 Mute we B More lt21 where xl 2 3At B M273 and B MFG We next need to compute the Lagrange multiplier which will constrain the uid particles to the sphere Using Eq 21 we obtain the grid point values of the momentum equation as follows quot251 M7131 Wichprgl ygl 22 Since we need to ensure that the velocity field remains tangential to the sphere we require xG uG 0 Taking the scalar product of Eq 22 with x39gl and rearranging yields 1 n I1 7 n n e 3M 1ltxG1Bugt M We m t lt23 Writing 1 Pq q 3 81 1 24 where the projection matrix P is given by 1 a2 x2 xy xz P 2 xy 12 yZ yz 25 d xz Z Z yz a z we can now write the constrained equations as M422 MDTu39gI B 26 Mug 1 1130ng PB 27 The projection matrix P constrains any vector quantity 1 to be tangential to the sphere Eqs 26 and 27 can now be written in the matrix form M APD m quot1 7 PB 28 MiDT M 42G B I Applying a block LU decomposition yields the following matrix system M APD m quot1 7 PB 29 0 M AzdiDTM IPD 42G B I MiDTPBu The twostep solution process is solve the pseudoHelmholtz equation 634 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 M xtzdiDTM IPD ltng B MiDTPBu 30 for p31 and then use this solution to solve u gl M IPBE AM IPDCPgR 31 This formulation only requires the solution of an Np gtlt Np matrix system corresponding to Eq 30 Fur thermore the associated matrix is symmetric positivedefinite and thus solvable by conjugate gradient CG methods We use GMRES with a fast projection method see 6 to solve this matrix only because this is what we happen to have developed at the moment The solution of Eq 31 is completely trivial due to the diagonal nature of M and thus only requires the computation of a 3Np vector 4 Grid generation on the sphere The elementbased interpolation operator used by SESL permits the use of any type of grid Defining the order of the SE basis functions automatically defines the order of the interpolation stencil as well This independence from the grid allows for the use of grids other than the typical latitude longitude grids used by spectral models We only use hexahedral grids in this paper but to show the method s independence on the grid we refer the reader to the papers 1112 in which we used icosahedral grids for the explicit Eulerian method The choice of using hexahedral grids in this paper is predicated on that fast searching algorithms can be constructed quite easily on this class of grids we are currently investigating similarly fast search algorithms for the icosahedral gridsi 41 Hexahedral grids Hexahedral aka cubic gnomonic grids are constructed by mapping the six faces of a hexahedron onto the surface of a sphere Each of the six faces of the hexahedron are then subdivided into the desired number of quadrilateral elements We begin by constructing a spectral element grid on the gnomonic space Gt G is defined by the square region X TL47 7L4Z in a 2D Cartesian space see 28 This region is divided into the elements and inside each element we construct the LGL grid points Upon constructing this grid we then map the gnomonic quot X to the Lui A quot spherical 1 Ag via AG X 32 7 arcsin tan Y CPG 1taan tanZ Y i It should be noted that 11G only gives the spherical coordinates of one of the six faces of the hexahedroni Therefore we have to rotate this face to the six faces of the hexahedron by the rotation mapping R 33 l x AC arctan cos 4 sm G cos ch cos AG cos pC sin ch sin pC q arcsinsin ch cos pC cos ch cos AG sin goo7 where the centroids Amcpc of the six faces are located at Amgoc c ln20 for c lui4 and 157425 Om2 1674 s Own2 This approach results in the construction of the hexahedral grid with the following properties Np 6nHNZ 2 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 635 Ne 6011027 where Np and Ne are the number of points and elements comprising the grid The parameter nH refers to the number of quadrilateral elements in each direction contained in each of the six faces of the hexahedron and N is the order of the Legendre cardinal functions of each element Table 1 gives some grid con gurations for the hexahedral grid Fig 1 shows sample hexahedral grids for N 4 8 and 16 with nH 4 42 Search algorithms The semiLagrangian method requires finding the location of the departure points EWH Unfortunately for quadrilateral elements there is no general way of finding the element which claims a specific departure point in triangles this is possible through the use of the natural coordinates such as in Giraldo 810 and Xiu and Karniadakis 37 Therefore the searching algorithms for quadrilateral elements must be speci cally tailored to the grid In the case of hexahedral grids we know that the grid is constructed by mapping each of the six faces of the hexahedron onto a local 2D planar Cartesian space Therefore the first step in determining the location of a departure point is to determine which of the six faces of the hexahedron claims the departure point We determine this face by checking which face centroid is closest to the departure point Then we map this face to the rotatedgnomonic space via Eq 33 and the inverse of Eq 32 This gives the departure point and the elements contained in this face in terms of the local 2D Cartesian coordinates X Because on this gnomonic projection each face is a perfect square then the element which claims the departure point can be obtained by using the following integer operations Xd E 4 1d l 1ntT 34 Table 1 The number of grid points and elements for the hexahedral grid as a function of nH and N m N Np N 1 4 98 6 1 8 386 6 1 16 1538 6 1 32 6146 6 4 4 1538 96 4 8 6146 96 4 16 24578 96 4 32 98306 96 Fig 1 A hexahedral grid for m 4 and a N 4 b N 8 and c N 16 636 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 jd lam 11 where idjd are the indices of the quadrilateral element which claims the departure point Xd Yd is the departure point in gnomonic space and AX AY nZnH The next step is to compute the departure point in the local computational space We require the local computational coordinate Ed because the Legendre cardinal functions are only defined in terms of these coordinates Using the indices of the element which claims the departure point given by Eq 34 we can obtain the departure points in computational space as follows 5dw17 35 Ara Y1 1 lt nd T 7 where X1 id 1AX 36 Y1 Ud 1AY This search algorithm requires a maximum of eight oating point operations and two integer operations regardless of the number of elements in the grid or the order of the polynomial basis N i 5 Results For the numerical experiments we use the normalized L1 L2 and LDo error norms M 9 Wm Pl d9 L1 f9 l l exaal 19 M f9 em Zd L2 f9 ngactdg 7 max q q PHLW l l max lq mal to judge the accuracy of the numerical scheme The following two conservation metrics are also used M chde E chpuZvZwZcpZd2 f9 4721 d9 7 f9 exactltugxact Dgxact ngact gogxact d9 7 where M measures the conservation property of the mass and E measures the conservation of the total energy In addition we show plots as a function of Courant number for the timestep studies The Courant number we use is the maximum Courant number of all the LGL cells To compute the Courant number the elements are decomposed into their LGL grid points and these grid points form a fine grid which we refer to as the LGL cells The velocities and grid spacings are then defined at the centers of these cells Using these definitions the Courant number is then defined as FX Giraldo er al Journal of Computational Physics 190 2003 623 650 637 CM 6 Courant number maX Ve E 1 Ne AS LGL where C U for Case 1 U g0 for Cases 2 34 5 and 6 where C is the characteristic speed U u u and AS Ax2 A322 A22 is the grid spacing Eight test cases are used to judge the performance of SESL Cases 1A 2 3 5 and 6 correspond to the Williamson et al 36 standard test case suite Case 4 has been used as a test case for the shallow water equations in 10 122324 Cases 1A 1B and 1C involve the geopotential equation only whereas the re mainder of the test cases concern the full shallow water equations In addition Cases 1 3 have analytic solutions and are used to determine the accuracy of SESL quantitatively Cases 4 6 on the other hand do not have analytic solutions and are thus only used to determine the accuracy of the scheme qualitatively For many of the test cases in this section we compare SESL with the spectral element semiimplicit method de ned in Eqs 5 and 6 which we refer to as SESI the value 6 12 is used unless otherwise speci ed We also use a fully eXplicit spectral element leapfrog model for comparisons which we refer to as SEMLF 51 Case IA Advection of a cosine wave with Steadyslate velocity Case 1A concerns the solid body rotation of a cosine wave The velocity eld remains unchanged throughout the computation For this case the velocity eld is de ned as u us sinl vs sin6cos 2 v uscosl vs sinBsin 2 w 128 cos 6 where us and US are the zonal and meridional velocity components in spherical coordinates given in 36 Results are obtained after one full revolution which corresponds to an integration of 12 days Fig 2 shows that SESL converges algebraically for this test case We cannot eXpect eXponential con vergence because the derivative at the base of the cosine hill is nonsmooth Normalized Error 4 8 16 32 N Fig 2 Case 1A The normalized g0 error norms for SESL as a function of N after 12 days using rm 1 638 EX Giraldo er al Journal of Computational Physics 190 2003 623 650 Recall that Falcone and Ferretti 7 proved that the error of Lagrangian high order methods is governed by Eq 17 In this section we show that the error in SESL is governed by Eq 17 we hope to show this for the atmospheric equations in future work Fig 3 shows the normalized g0 L error norm for the K 2 and K 4 trajectory approximations using the nH 1 and N 16 grid for timesteps yielding Courant numbers from 01 to about 50 The K 2 scheme is the second order Runge Kutta method given by Eq 15 while the K 4 scheme is the fourth order Runge Kutta method given by Eq 16 The error of the K 2 scheme begins to increase at a Courant number of 25 This is due to the error being dominated by MK On the other hand the error for the K 4 scheme remains level throughout for Courant numbers below 25 Thus by simply increasing the order of K we have been able to sustain the same level of high accuracy for an order of magnitude greater time step 52 Case 1B Advection of a cosine wave with limedependent velocity This case is similar to Case 1A except that the velocity eld varies with time in the following manner u uS sin1 vs sin6cos x1 1 COS 12 gays v 11S cos1 vs sin6sin11 COS 12 gays 2 t w vscos61cos12 gays This case was introduced in 3 and the exact solution can be found in 10 Fig 4 shows that SESL con verges rapidly even for ow with nonconstant velocity To show that the Falcone and Ferretti result also holds for ows with timedependent velocity elds we plot in Fig 5 the normalized g0 L2 error norm as a function of timestep for nH 1 and N 16 Once again by increasing the order of K we are able to maintain the same level of high accuracy for a much larger time step from a Courant number of 4 for K 2 to a Courant number of 10 for K 4 10 NormalizedL Error 2 O 01 I I I39M 100 Courant Number Fig 3 Case 1A The normalized g0 L2 error norm for SESL using RK2 solid and RK4 dashed for the trajectory computation as a function of Courant number after 12 days using rm 1 and N 16 FX Giraldo ct al Journal Of Computational Physics 190 2003 623 650 639 Normalized Error N Fig 4 Case 1B The normalized g0 error norms for SESL as a function of N after 12 days using rm 1 10 NormalizedL Error 2 O 01 1 10 100 Courant Number Fig 5 Case 1B The normalized g0 L2 error norm for SESL using RK2 solid and RK4 dashed for the trajectory computation as a function of Courant number after 12 days using rm 1 and N 16 53 Case 1C Advection of a circular column In this case we use the steadystate velocity eld from Case 1A but the geopotential is now a circular column Unlike a cosine wave the circular column represents a discontinuous function and is the 2D analog of a 1D square wave For ows containing large gradients andor discontinuities it is well known that high order methods produce oscillations near the discontinuities which contaminate the entire ow eld ie Gibbs phenomena We illustrate this test case to show the ability of the semiLagrangian method to handle nonsmooth ows In the atmosphere severe fronts can form and any method that can handle this test case should be adequate for handling the strongest atmospheric fronts Fig 6 shows the algebraic convergence rate of SESL for ows with sharp discontinuities In Fig 7 we plot the normalized g0 L2 error norm as a function of timestep for nH l and N 16 The K 2 scheme begins to lose accuracy at a Courant number of 25 Meanwhile the K 4 scheme retains the same accuracy up to Courant numbers of 50 The results for Case I con rm that using large K permits increasing the timestep without increasing the error However this comes at a considerable cost especially for the full nonlinear equations The number of departure point searches per grid point scales like 9K and the number of iterations required to achieve convergence in the elliptic solver increases signi cantly with increasing timestep due to the growth of the condition number We found the most ef cient Courant numbers for the shallow water equations to be in 640 EX Giraldo er al Journal of Computational Physics 190 2003 623 650 101 Normalized Error 0 Fig 6 Case 1C The normalized g0 error norms for SESL as a function of N after 12 days using rm 1 10 u ICI 39D G E a E O Z 2 RK2 RK4 1 1 39 39 001 1 10 100 Courant Number Fig 7 Case 1C The normalized g0 L2 error norm for SESL using RK2 solid and RK4 dashed for the trajectory computation as a function of Courant number after 12 days using rm 1 and N 16 the range of 10 20 which can be achieved most ef ciently and with suf cient accuracy using K 2 Thus for the remainder of the paper we use SESL with K 2 To see how the spectral element and the spectral element semiLagrangian method handle sharp dis continuities let us compare g0 along the Equator for SESI and SESL using the grid nH 4 and N 16 Figs 8 and 9 show these slices after one complete revolution of the circular column Figs 8a and b show the I I l I I I O 90 1 80 270 360 O 90 180 270 36C 3 7h b 7r Fig 8 Case 1C The mass go along the Equator after 12 days using rm 4 and N 16 for a SESI and b SESI with ltering FX Girala o et al Journal of Computational Physics 190 2003 623 650 641 100 100 50 O 90 10 270 360 O 90 10 270 36C Fig 9 Case 1C The mass go along the Equator after 12 days using nH 4 and N 16 for a SESL and b SESL with FCT results for SESI with and without ltering using a Courant number of 05 This is the maximum timestep allowed by SESI for this test case In Fig 8a it is quite evident that the large gradients at the front and lee of the discontinuous wave are causing severe oscillations In fact the entire domain is contaminated by Gibbs phenomena In order to eliminate these aliasing errors high order methods typically employ some form of ltering In Fig 8b we plot the result for SESI using the Boyd Vandeven lter see 212 The lter has eliminated a substantial amount of the numerical noise however there are still some undershoots in the front and lee of the wave and some disturbing oscillations at the top of the wave The most disturbing aspect of this solution is that the correct extrema go 6 0 100 are violated Thus if this equation repre sented moisture in the atmosphere we would overpredict the amount of rain expected in a certain region of the world and obtain negative rain rates in others Figs 9a and b show the results for SESL and SESLFCT using a Courant number of 5 with no ltering SESL has eliminated almost all of the oscillations although small undershoots and overshoots still occur Fig 9b shows that SESL with the FCTtype monotone interpolation does not introduce any undershoots or overshoots 54 Case 2 Steadyslate nonlinear zonal geostrophic ow This case is a steadystate solution to the nonlinear shallow water equations The equations are geo strophically balanced and remain so for the duration of the integration The velocity eld thus remains constant throughout the computation The geopotential height g0 undergoes a solid body rotation but since the initial height eld is given as a constant band in the zonal direction and the ow eld is purely zonal then the solution remains the same throughout the time integration The velocity eld is the same as that used in Case 1 Williamson et a1 36 recommend that the error be computed after ve days of integration Fig 10 illustrates the normalized g0 error norms as a function of polynomial order N for SESL The errors are decreasing quite rapidly Fig 11 compares the normalized g0 L2 error norm of SESI and SESL as a function of the timestep for the grid nH 4 and N 16 This gure shows that SESL clearly retains its accuracy for signi cantly large Courant numbers while the accuracy of SESI diminishes quickly and eventually blows up for Courant numbers larger than 25 Note that SESL can use a Courant number 10 times larger than SESI Courant number of 25 and yield the same accuracy as SESI for a Courant number of 25 55 Case 3 Steadyslate nonlinear zonal geostrophic ow with compact support This case is another steadystate solution to the nonlinear shallow water equations The equations remain geostrophically balanced for the duration of the integration The initial velocity eld is zero everywhere except 642 EX Giraldo er al Journal of Computational Physics 190 2003 623 650 Normalized Error Fig 10 Case 2 The normalized g0 error norms for SESL as a function of N after ve days using rm 1 Error O 2 Normalized L SESL 01 I I H m1 00 Courant Number Fig 11 Case 2 The normalized g0 L2 error norm for SESI dashed and SESL solid as a function of Courant number after ve days using rm 4 and N 16 in a very small isolated region This isolated region or jet encapsulates the ow and limits the geopotential height eld to remain Within a very con ned circular region The results are reported for a ve day integration as suggested in 36 Fig 12 illustrates the normalized g0 error norms as a function of polynomial order N for SESL SESL converges exponentially for this test case Fig 13 compares the normalized g0 L2 error norm of SESI and Normalized Error N Fig 12 Case 3 The normalized g0 error norms for SESL as a function of N after ve days using rm 1 FX Girala o et al Journal of Computational Physics 190 2003 623 650 643 10 3E 4 Error l O 39139 I I39ll 2 Normalized L SESL 939 I I I 01 1 10 100 Courant Number Fig 13 Case 3 The normalized g0 L2 error norm for SESI dashed and SESL solid as a function of Courant number after ve days using nH 4 and N 16 SESL as a function of the timestep for the grid nH 4 and N 16 This gure shows that SESL clearly retains its accuracy for signi cantly large Courant numbers while the accuracy of SESI diminishes quickly and eventually blows up for Courant numbers larger than 2 SESL can achieve the same order of accuracy for Courant numbers of about 10 56 Case 4 Dancing high low waves This case was introduced in McDonald and Bates 23 and is not an analytic solution to the shallow water equations The initial geopotential height is comprised of two large waves with the left wave being the low wave and the right wave being the high wave when viewed from the north pole The waves rotate clockwise in a swirling dancelike fashion Fig 14 shows the contours of go and the zonal and meridional velocities 113123 for SEMLF SESI and SESL after 10 days The Eulerian models SEMLF and SESI are run using their maximum allowable timesteps The following Courant numbers are used for the three models 05 for SEMLF Courant number 14 for SESI 140 for SESL For SESI the advection terms which are integrated explicitly limits the timestep to only three times that of the purely explicit leapfrog scheme SEMLF SESL uses a timestep 10 times larger than SESI and almost 30 times larger than SEMLF The results show almost no differences between the three models In ad dition for 30 day integrations not shown SEMLF conserves mass exactly and energy within 008 SESI conserves mass within 0005 and energy within 010 and SESL conserves mass within 002 and energy within 075 After 30 days all three models yield almost identical results with go only varying by less than 085 between the models 5 7 Case 5 Zonal ow over an isolated mountain This case uses the same initial conditions as Case 2 with the addition of a conical mountain at xi 6 180 30 where xi is the zonal direction and 6 the meridional direction Due to the zonal ow impinging on the mountain various wave structures form and propagate throughout the sphere Fig 15 shows the contours of go as and Us for SEMLF SESI and SESL after 15 days which is the maximum timelength suggested in 36 The following Courant numbers are used for these results 644 EX Giraldo er al Journal of Computational Physics 190 2003 623 650 Fig 14 Case 4 Contours of go left us center and vs right for a SEM LF b SESI and c SESL after 10 days using nH 4 and N 16 05 for SEMLF Courant number 15 for SESI 120 for SESL SESI uses a timestep three times larger than SEMLF while SESL uses a timestep eight times larger than SESI and 24 times larger than SEMLF Once again the results show almost no differences between the models especially between the two implicit models SESI and SESL SEMLF shows some differences but the underlying wave structures of the ow are clearly evident in all three models For 30 day integrations the conservation metrics are as follows SEMLF conserves mass exactly and energy within 070 SESI conserves mass within 003 and energy within 057 and SESL conserves mass within 003 and energy within 055 After 30 days the models appear slightly different although g0 only varies by less than 03 between the models SESI and SESL yield virtually identical solutions 58 Case 6 Rossby Haurwitz wave Although Rossby Haurwitz waves are not analytic solutions to the shallow water equations they have become a de facto test case In a nondivergent barotropic model the initial geopotential height eld undergoes a solid body rotation in a counterclockwise direction when viewed from the north pole Although this case does not have an analytic solution it is wellknown that the initial wave structure of the Rossby Haurwitz wave should remain intact for the duration of the timeintegration Thus we know that for at least 30 days the waves should persist without any loss in shape or mass For this reason this case represents the best test for judging the performance of the model for long timeintegrations FX Giraldo er al Journal of Computational Physics 190 2003 j 623 650 645 Fig 15 Case 5 Contours of go left us center and vs right for a SEM LF b SESI and c SESL after 15 days using nH 4 and N 16 Fig 16 shows the contours of g0 us and Us for SEMLF SESI and SESL after 14 days The following Courant numbers are used for these results 032 for SEMLF Courant number 103 for SESI 103 for SESL SESI uses a timestep three times larger than SEMLF while SESL uses a timestep 10 times larger than SESI and 30 times larger than SEMLF Once again the results show almost no differences between the models For 30 day integrations the conservation metrics are as follows SEMLF conserves mass exactly and energy within 037 SESI conserves mass within 003 and energy within 057 and SESL conserves mass within 003 and energy within 040 After 30 days the results of the three models are virtually indistinguishable with no variation in g0 However for SESI to run for 30 days required offcentering to a value of 6 08 The most important result is that all three models predict the same rate of rotation This is an extremely important result because it con rms that the implicit models are not damping the Rossby waves Rossby waves are the most important waves for largescale meteorological processes and are often referred to as planetary waves Rossby Haurwitz waves are the Rossby waves on the sphere Thus any numerical method being considered for constructing atmospheric models must handle these waves properly 646 EX Giraldo er al Journal of Computational Physics 190 2003 623 650 Fig 16 Case 6 Contours of go left us center and vs right for a SEM LF b SESI and c SESL after 14 days using nH 4 and N 16 5 9 Computational cost In this section we compare the computational costs of the three models SEMLF SESI and SESL We shall use the explicit model SEMLF as the reference with which to compare the Eulerian implicit model SESI and the Lagrangian implicit model SESL Fig 17 illustrates the ratio of computer time required by SESI relative to SEMLF dashed line and SESL relative to SEMLF solid line as a function of N for nH 1 These time ratios are based on a 24 h integration for Case 5 with the models running the following Courant numbers 040 for SEMLF Courant number 150 for SESI 1200 for SESL We chose Case 5 because it is the only test case with topography and thereby is the most realistic test case of the entire suite SEMLF is anywhere between 125 and 4 times more ef cient than SESI In contrast SEM LF is only more ef cient than SESL for N gt 8 In addition SESL is more ef cient than SESI for N g 12 This study seems to indicate that SESL should be used with N g 12 which is in agreement with our previous experience of cost versus accuracy see 1112 where we found the range 8 g N g 16 to be optimal for the construction of atmospheric models see 13 In Fig 18 we plot the time ratio comparing SESI and SESL to SEMLF as a function of ml with N 8 and using the same Courant numbers given above SESI uses a timestep almost four times larger than SEMLF while SESL uses a timestep eight times larger than SESI and 32 times larger than SEMLF FX Giraldo er al Journal Of Computational Physics 190 2003 623 650 647 SESISEM LF 35 SESLSEM LF Time Ratio 4 1392 1396 2394 32 N Fig 17 The ratio of computational cost for SESLSEM LF dashed and SESLSESI solid for various N for Case 5 using nH 1 The ratios re ect the amount of computer time required to complete a 24 h time integration SESISEM LF SESLSEMLF Time Ratio L l 2 4 8 16 Fig 18 The ratio of computational cost for SESISEM LF dashed and SESLSEM LF solid for various nH for Case 5 using N 8 The ratios re ect the amount of computer time required to complete a 24 h time integration Fig 18 shows that SESL is as high as 70 more ef cient than SESI for nH l6 SESL is 30 less ef cient than SEMLF for nH lt4 However for nH gt 4 SESL becomes increasingly more ef cient than SEMLF with increasing nH For nH l6 SESL is 7 more ef cient than SEMLF This is very good news because ml is the parameter which controls the grid resolution for constant N Thus based on the results outlined in Fig 18 we can expect SESL to continue to increase its performance with increasing grid res olution We certainly expect to be able to increase the performance of SESL beyond the current level To further increase the performance of SESL requires increasing the performance of SESI Increasing the perfor mance of SESI will increase the performance of SESL because both models use the same solver and preconditioner While we did not see a performance increase with SESI over SEMLF Thomas et a1 34 have shown that for the shallow water equations Eulerian semiimplicit methods can be constructed to be more ef cient than explicit methods on parallel computers We hope to learn from their experiences and use similar techniques in order to achieve the same kind of performance In this paper we present the results for SESL to show proofofconcept only and in future work we shall report on developments in the solvers preconditioners and search algorithms in the context of a 3D atmospheric model on parallel computers 648 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 6 Conclusions In this paper we present a high order Lagrangian method for solving advectiondominated uid ows which may be wellsuited for constructing future atmospheric models The highlight of this method is that it uses the high order basis functions of the spectral element method as the interpolation functions required by the Lagrangian method to interpolate values at the Lagrangian departure points This elementbased compact interpolation renders the spectral element semiLagrangian SESL method quite general In addition to high order accuracy and generality locality is another property inherent to SESL High order accuracy means that SESL can achieve exponential convergence for smooth ows Generality means that SESL can be used on adaptive unstructured grids Locality means that SESL can be ported to distributed memory computers quite easily In this paper we have shown the high order accuracy property We have only shown results on hexahedral grids but it is immediately obvious that any type of quadrilateral grid can be used generality property The parallel implementation locality property of SESL is the topic of future work The results obtained for all eight cases show that SESL is not only high order accurate but possibly very ef cient due to the large timesteps that it can use The results show that SESL yields good convergence for increasing order of the basis functions In addition the results show that SESL with an FCTtype inter polation scheme is capable of resolving nonsmooth ows such as those encountered in atmospheric fronts The results obtained by SESL with FCT were shown to be better than those obtained by the spectral el ement method with ltering We compare SESL with an Eulerian semiimplicit method and show that these two methods yield results that are almost identical even when SESL uses timesteps 8 10 times larger These timesteps are 20 30 times larger than that of the explicit Eulerian model SEMLF The performance comparisons revealed that SESL is 70 more ef cient than an Eulerian semiimplicit model but is only 7 more ef cient than an Eulerian explicit model We hope to develop faster solvers preconditioners and search algorithms in order to further increase the performance of SESL and we will report our progress in future work Acknowledgements The first author FXG was supported by the Of ce of Naval Research through program el ement PE0602435N The second author JBP was supported in part by the Air Force Of ce of Scienti c Research F49620001003 and the Of ce of Naval Research N0001401l0483 The effort of the third author PFF was funded by the US Department of Energy Of ce of Ad vanced Scienti c Computing under the SciDAC program Terascale Simulation Tools and Tech niques TSTT Center References 1 R Bermejo A Staniforth The conversion of semi Lagrangian advection schemes to quasi monotone schemes Monthly Weather Review 120 1992 2622 2632 2 J Boyd The erfc log lter and the asymptotics of the Euler and Vandeven sequence accelerations Houston Journal of Mathematics 1996 3 G Chukapalli Weather and climate numerical algorithms a uni ed approach to an efficient parallel implementation PhD Dissertation University of Toronto 1997 4 J C3te A Lagrange multiplier approach for the metric terms of semi Lagrangian models on the sphere Quarterly Journal of the Royal Meteorological Society 114 1988 1347 1352 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 649 5 J Cote M Roch A Staniforth L Fillion A variable resolution semi Lagrangian nite element global model of the shallow water equations Monthly Weather Review 121 1993 231 243 6 PF Fischer Projection techniques for iterative solution of Ax b with successive right hand sides Computer Methods in Applied Mechanics and Engineering 163 1998 193 204 7 M Falcone R Ferretti Convergence analysis for a class of high order semi Lagrangian advection schemes SIAM Journal of Numerical Analysis 35 1998 909 940 8 FX Giraldo Lagrange Galerkin methods on spherical geodesic grids Journal of Computational Physics 136 1997 197 213 9 FX Giraldo The Lagrange Galerkin spectral element method on unstructured quadrilateral grids Journal of Computational Physics 147 1998 114 146 10 FX Giraldo Lagrange Galerkin methods on spherical geodesic grids the shallow water equations Journal of Computational Physics 160 2000 336 368 1 FX Giraldo A spectral element shallow water on spherical geodesic grids International Journal for Numerical Methods in Fluids 35 2001 869 901 12 FX Giraldo JS Hesthaven T Warburton Nodal high order discontinuous Galerkin methods for the spherical shallow water equations Journal of Computational Physics 181 2002 499 525 13 FX Giraldo TE Rosmond A scalable spectral element Eulerian atmospheric model SEE AM for fNWP dynamical core test Monthly Weather Review 2003 in press 14 JJ Hack BA Boville BP Briegleb JT Kiehl PJ Rasch DL Williamson Description of the NCAR community climate model CCM2 NCAR Technical Note NCARTN 382STR National Center for Atmospheric Research Climate Modeling ection PO Box 3000 Boulder CO 80307 1992 15 R Heikes DA Randall Numerical integration of the shallow water equations on a twisted icosahedral grid Part I Basic design and results of tests Monthly Weather Review 123 1995 1862 1880 16 T Heinze Hense The shallow water equations on the sphere and their Lagrange Galerkin solution Meteorology and Atmospheric Physics 81 2002 129 137 17 TF Hogan TE Rosmond The description of the navy global operational prediction system s spectral forecast model Monthly Weather Review 119 1991 1786 1815 18 M Iskandarani DB Haidvogel JP Boyd Staggered spectral element model with application to the oceanic shallow water equations International Journal for Numerical Methods in Fluids 20 1995 393 414 19 M Kwizak AJ Robert A semi implicit scheme for grid point atmospheric models of the primitive equations Monthly Weather Review 99 1971 32 36 20 SJ Lin RB Rood An explicit ux form semi Lagrangian shallow water model on the sphere Quarterly Journal of the Royal Meteorological Society 123 1997 2477 2498 21 AV Malesky SJ Thomas Parallel algorithms for semi Lagrangian advection International Journal for Numerical Methods in Fluids 25 1997 455 473 22 Y Maday AT Patera EM Ronquist An operator integration factor splitting method for time dependent problems application to incompressible uid ow Journal of Scienti c Computing 5 1990 263 292 23 A McDonald JR Bates Semi Lagrangian integration of a gridpoint shallow water model on the sphere Monthly Weather Review 117 1989 130 137 24 B Neta FX Giraldo IM Navon Analysis of the Turkel Zwas scheme for the 2D shallow water equations in spherical coordinates Journal of Computational Physics 133 1997 102 122 25 H Ritchie Semi Lagrangian advection on a Gaussian grid Monthly Weather Review 115 1987 608 619 26 A Robert A semi Lagrangian and semi implicit numerical integration scheme for the primitive meteorological equations Journal of the Meteorological Society ofJapan 60 1982 319 325 27 M Rochas ARPEGE Documentation Part 2 Chapter 6 available from Meteo France 1990 28 C Ronchi R Iacono PS Paolucci The cubed s here a new method for the solution of partial differential equations in spherical geometry Journal of Computational Physics 124 1996 93 114 29 JG Sela Spectral modeling at the National Meteorological Center Monthly Weather Review 108 1980 1279 1292 30 AJ Simmons DM Burridge M Jarraud C Girard W Wergen The ECMWF medium range prediction models development of the numerical formulations and the impact of increased resolution Meteorology and Atmos heric Physics 40 1989 28 60 31 RN Swarztrauber DL Williamson JB Drake The Cartesian method for solving partial differential equations in spherical geometry Dynamics of Atmospheres and Oceans 27 1997 679 706 32 M Taylor J Tribbia M Iskandarani The spectral element method for the shallow water equations on the sphere Journal of Computational Physics 130 1997 92 108 33 C Temperton A Staniforth An efficient two time level semi Lagrangian semi implicit integration scheme Quarterly Journal of the Royal Meteorological Society 113 1987 1025 1039 34 SJ Thomas RD Loft JM Dennis Parallel implementation issues global versus local methods Computing in Science and Engineering 4 2002 26 31 650 EX Giralda et al Jaumal 0f Camputatianal Physics 190 2003 623 650 35 E Turkel G Zwas Explicit large time step schemes for the shallow water equations in R VichneVetsky RS Stepleman Eds Advances in Computer Methods for Partial Differential Equations IMACS Lehigh University 1979 p 65 36 DL Williamson JB Drake JJ Hack R Jakob PN Swarztrauber A standard test set for numerical approximations to the shallow water equations in spherical geometry Journal of Computational Physics 102 1992 211 224 37 D Xiu GE Karniadakis A semi Lagrangian high order method for the Navier Stokes equations Journal of Computational Physics 172 2001 658 684
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'