DTIC ADA603963: Grid Effect on Spherical Shallow Water Jets Using Continuous and Discontinuous Galerkin Methods

20.5 MB·English

by  Defense Technical Information Center

Quarterly Journal of the Royal Meteorological Society Q. J. R. Meteorol. Soc. 00: 1–16 (2013)
Grid effect on spherical shallow water jets using
continuous and discontinuous Galerkin methods
∗
S. Marras , M. A. Kopera, F. X. Giraldo
Department of Applied Mathematics, Naval Postgraduate School
∗Correspondence to: 833 Dyer Road, Bldg. 232, Spanagel 249A, Monterey, CA 93943 U.S.A.
[email protected]; [email protected]
Element-based Galerkin methods are not constrained by the logical
structure of the computational grid. In this paper, the behavior of the
continuous and discontinuous Galerkin solutions of the shallow water
equationson thesphereis analyzedfor threedifferentgridsof ubiquitoususe
in atmospheric modeling: (A) hexahedron (or cubed sphere), (B) reduced
latitude-longitude, and (C) icosahedron. Conforming and nonconforming
mesh configurations are adopted. The nonconforming grids are based on
a quad-tree structure. The analysis is performed on the mid-latitude zonal
flow problem suggested by Galewsky et al. [An initial-value problem for testing
numerical models of the global shallow-water equation, Tellus 2004; 56A:429-440].
This test is sufficiently simple in its implementation, yet sufficiently complex
to capture some important modes that arise on the global scales of real
atmospheric dynamics. Because the inviscid solution on certain grids shows
a high sensitivity to the resolution, the viscous counterpart of the governing
equations are also solved and the results are compared. The study shows
that not only the resolution, but the grid alignment with respect to the flow
is important in order to obtain an accurate solution. This is especially true if
the governing equations are not regularized by the addition of a sufficiently
large, fully artificial, diffusion mechanism; we give a brief (not exhaustive)
analysis of the effect of viscosity on the solution quality. The flexibility and
accuracyofelement-basedGalerkinmethodsongeneralsphericalgeometries
is illustrated; we particularly emphasize the excellent properties of the
reduced Lat-Lon configuration in comparison with the cubed-sphere on the
one hand, and with the more regular and uniform icosahedral grid on the
other.
KeyWords: ReducedLat-Longrids;IcosahedralGrid;Cubed-Sphere;ShallowWaterEquations;Galerkin
Methods; Spectral Element; Discontinuous Galerkin; Barotropic Jet; Galewsky Flow; Grid Generation on
theSphere;TransfiniteInterpolation.
Received ...
1. Introduction from now on) and discontinuous Galerkin (DG) solutions
of the shallow water equations (SWE) on a set of common
gridsusedinglobalcirculationmodels(GCMs).Specifically,
As computational power increases, meteorologists demand
we study 1) the cubed-sphere (Hex) (Sadourny 1972),
greater resolution from global atmospheric models utilizing
2) the icosahedron (Ico) (Sadourny et al. 1968), and 3)
spherical domains. Driven by the need to select the
a reduced latitude-longitude (Lat-Lon) geometry (Phillips
most proper computational grid on the sphere for the
1957). Because SWE capture many of the essential features
Nonhydrostatic Unified Model of the Atmosphere (NUMA)
of GCMs while eliminating the vertical structure of the
(Kelly and Giraldo 2012; Giraldo et al. 2013), in this paper
atmosphere, the results from this study will be directly
we analyze the continuous Galerkin (or spectral elements
applicable to the fully three-dimensional model NUMA.
The logically structured Hex and Lat-Lon grids can
†Pleaseensurethatyouusethemostuptodateclassfile,availablefrom
be advantageous for a finite difference based solver, and
theQJRMSHomePageat
http://onlinelibrary.wiley.com/journal/10.1002/(ISSN)1477-870X are generally usable by grid point methods such as finite
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls[Version: 2013/10/14v1.1]
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
2013 2. REPORT TYPE 00-00-2013 to 00-00-2013
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
Grid effect on spherical shallow water jets using continuous and
5b. GRANT NUMBER
discontinuous Galerkin methods
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 Postgraduate School,Department of Applied REPORT NUMBER
Mathematics,Monterey,CA,93943
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

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 16
unclassified unclassified unclassified Report (SAR)
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
2 Marras et al.
and spectral element (FE, SE), discontinuous Galerkin with constant coefficient ν =1 105m2s−1 as suggested
×
(DG), or finite volume (FV) methods. Furthermore, Lat- by Galewsky et al. (2004). The positive effect of diffusion
Lon is the standard discretization for GCMs based on is clear for the low resolution results on the cubed-sphere
spectraltransformmethodsduetothefastFouriertransform and icosahedron; on the contrary, we will show how this
operation along the longitude direction (Hogan et al. 1991). is not completely true for the solution on the reduced
Icosahedral grids or, more generally, unstructured tri- and Lat-Lon grid. It must be pointed out that the diffusion
quad-based tessellations∗ of the sphere are geometrically correction that was adopted is fully artificial and has no
flexible, but not all numerical methods are able to handle physical meaning. The sole reason for its use in this study
them. (A review of different grid generation techniques is to be able to analyze the solution when the grid is
in atmospheric modeling can be found in the book by not of sufficient resolution to properly resolve the flow.
Behrens (2006).) Recently, different grids to be used Partially based on this study, but mostly on the analysis
with finite volume discretization techniques have been made by practitioners in computational fluid dynamics for
analyzedin,e.g.,Qaddouri et al.(2012);Weller et al.(2012) more than three decades, we strongly believe that simple
or Peixoto and Barros (2013), with a review appearing Laplace operators (of any order) for stabilization purposes
in Staniforth and Thuburn (2012). In the context of should be avoided for the sake of accuracy. If viscosity is
Galerkin methods, however, not much analysis has been indeed necessary, techniques similar to those described in
conducted with respect to the solution dependence on these Malm et al.(2013)orMarras et al.(2013b,a),amongothers,
computational grids. Rather, given a specific grid, we are may be valuable options to consider for the solution of the
likely to find a model that is developed around it and is shallowwaterandEulerequations.Researchinthisdirection
optimized to minimize the error in a specific configuration. is still on-going; especially so in atmospheric applications.
In the case of the solution of shallow water problems, Finally, to assess the ability of the algorithm to handle
examples of this approach are the works by Taylor et al. conforming and nonconforming grids with refinement, we
(1997) (SE on cubed-sphere), Lanser et al. (2000) (FV on a complete the study by building six statically adapted grids.
reducedLat-Longeometry),Giraldo(2001)(SEonthequad- We constructed a fixed high-resolution grid core in the
basedicosahedron),Giraldo et al.(2002)(DGonatriangle- region where the zonal jet is most likely to develop, and
based icosahedron), or Nair et al. (2007) (DG on the cubed coarsened the rest of the domain with a 1- and 2-level de-
sphere). One step forward in grid comparison was made refinement approach to estimate the error when the total
by Giraldo and Warburton (2005) who compared different numberofgridpointsisdrasticallyreduced.TheDGsolution
triangle-based unstructured grids amongst each other and, algorithm on the refined grids is based on the procedure by
subsequently,againstthecubedsphere.Theyconcludedthat Kopera and Giraldo (2013) and references therein.
the accuracy of the solution is clearly a function of the The paper is organized as follows. In Sec. 2, we report on
combinationofthegridandtheproblemthatisbeingsolved, theconstructionofhigh-ordergridsonthesphere.Themodel
and it does not depend on either the numerical method or equations are described in Sec. 3 whereas the numerical
method of solution is described in Sec. 4. We describe the
the grid alone. A comparison between dynamically adaptive
testsandanalyzetheresultsinSec.5.Wefinallysummarize
and uniform meshes was performed by Müller et al. (2013)
our findings in Sec. 6.
for triangle based discontinuous Galerkin methods.
In2008,St-Cyret al.compared twoshallowwatersolvers
2. High order grid generation on the sphere
based on spectral elements on the one hand, and finite
volumes on the other (see St-Cyr et al. (2008)). The two
Grid generation on the sphere is a relatively easy task.
models weredesigned toexecute ondifferent grids:acubed-
However, finding the grid that can suit different numerical
sphere (SE solver) and on a Lat-Lon grid (FV model). In
methods is not. Because element-based Galerkin methods
their study, they assert that the spectral element method
havetheadvantageofbeingflexiblewithrespecttothegrid,
is unable to capture the solution at certain resolutions.
weanalyzehowdifferentmeshesonthespherecanaffectthe
In this paper we show that the inaccuracies that they
spectral element and discontinuous Galerkin solution of the
encountered are not linked to the numerical method per se,
shallow water equations.
butrather,totheunfortunatecombinationoftheresolution
Astandardapproachtohighordergridgenerationisbased
andalignmentofthemeshwithrespecttothedynamicsthat
on the construction of a linear grid first (i.e. a grid with
was being resolved in that specific study. We will confirm
straight edges whose extrema are the only points that lie on
their conclusion that a coarse resolution is responsible for
the geometry), and then to populate it with the high order
a wavenumber-4 signal that destroys the solution on the
internal points. The population step is performed element-
hexahedral grid; nevertheless, the simple change of grid will
wise in the logical space obtained by a proper projection
disproof the responsibility of SE and, rather, show that
from the physical space. Once the higher order grid points
the instabilities are merely an undesired effect of the grid
havebeenbuiltontheplane,abackwardprojectionontothe
configuration.
physical geometry moves the new high order elements onto
This study is performed by solving the two nonlinear thesphere.Thisprocessissuchthatthehighorderelements
zonal flows proposed by Galewsky et al. (2004). At low approximate the surface faithfully. The projection from the
resolution, the flows on the cubed-sphere and icosahedron sphere to the logical space may differ from grid to grid.
completely diverge from the true solutions in the case of (i) The gnomonic projection by Ronchi et al. (1996) was used
the equilibrium, and (ii) the barotropic instability tests. As to build the high order grids. In what follows, we describe
the resolution is increased, the effect of the irregular grid how the RLL, Hex, and Ico grids were constructed.
geometryispartiallycamouflagedandeventuallydisappears.
Insearchofthecoarsestallowableresolutiontobeusedwhile 2.1. Reduced Lat-Lon grid (RLL)
still preserving accuracy, we added an artificial diffusion
Since the effort of Phillips (1957) to reduce the singularity
problem at the poles of a global latitude-longitude grid,
∗Thekeywordsgrid,mesh,andtessellationwillbeusedinterchangeably
throughoutthepaper. different types of reduced Lat-Lon grid (RLL from now
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
Grid effect on spherical shallow water jets 3
on) have been used. Partially based on the composite TobuildeachinterfacesurfacebyTFIweneedinformation
methods of Starius (1977) and Browning et al. (1989), from its four boundary edges. With reference to the
Lanser et al. (2000) solved the shallow water equations schematic of Fig. 2, the region delimited by edges 1,2,3,4
on a Lat-Lon grid away from the poles combined with a is built in two steps: (a) given the (λ,ϕ) coordinates of the
stereographic grid at the poles. In this paper, we build cornerpoints1-2-3-4,buildedges3and4andsubdividethem
and test a modified, high order version of Lanser’s grids. into N 1D linear elements, where N is the user-
elat,int elat,int
The interfaces between the reduced Lat-Lon area and the definednumberofelementsalongthelatitudedirectioninthe
polar caps are obtained by a transfinite interpolation (TFI) interfaceregion;(b)computethegridpointsintheinteriorof
(Gordon and Hall 1973; Eriksson 1984). The way this is thepatchusingtheplanarandlineartransfiniteinterpolation
done will be described shortly. To build the RLL grid we defined by the Boolean sum
used a multiblock approach (see, e.g., Thompson (1987))
that consists of building a set of independent, simply X(ξˆ,ηˆ)=U+V U V, (2)
− ⊗
connectedsurfacegridsthatareeventuallypatchedtogether
where, given the arrays ξˆ=0,...,1 and ηˆ=0,...,1 in
to form the global mesh. The main RLL region and the
computational space along the local (uˆ,vˆ) directions, the
polarcapsarethefirstblockstobebuilt.ThemainLat-Lon
univariate interpolations and tensor products ( ) are
region is composed of four faces obtained from the master ⊗
computed as
face γ =[ π/4 λ π/4] [ϕ ϕ ϕ ] and from
1 min max
− ≤ ≤ × ≤ ≤
its rotation with respect to the z-axis of the sphere. The
variables λ and ϕ indicate longitude and latitude. The RLL
U(ξˆ,ηˆ)=(1 ξˆ)X(0,ηˆ)+ξˆX(1,ηˆ)
band in the northern and southern hemispheres is delimited −
byϕ andϕ .Therotationmatrixfromγ tothethree V(ξˆ,ηˆ)=(1 ηˆ)X(ξˆ,0)+ηˆX(ξˆ,1)
min max 1 −
remaining faces γ ,γ ,γ is obtained by the combined effect U V(ξˆ,ηˆ)=(1 ξˆ)(1 ηˆ)X(0,0)+(1 ξˆ)ηˆX(0,1)
2 3 4
⊗ − − −
of a translation on the xy-plane and a rotation around the (1 ηˆ)ξˆX(1,0) ξˆηˆX(1,1). (3)
z-axis; this transformation is given by − − −
In Equations (2), X is the array of the (λ,ϕ) coordinates
1 0 0 0 cos(λ ) sin(λ ) 0 0
rot rot of the grid points in the interior of the surface,
[T] = 0 1 0 0−sin(λrot) cos(λrot) 0 0. given the known values of the surface boundary edges
0 0 1 0 0 0 1 0 stored in X(0,ηˆ),X(1,ηˆ),X(ξˆ,0),X(ξˆ,1), and corners
   X(0,0),X(1,1),X(1,0),X(0,1).
l m n 1 0 0 0 1
  
  
Thefirstmatrixproducesatranslation[lmn]=[ 1 1
− − −
1]along(x,y,z)andisfollowedbytherotationsλ =
rot,i=2,3,4
(π/2,π,3π/2)givenbytheactionofthesecondmatrix.The
coordinatesoneachrotatedplanarfacearesimplyafunction
of (x,y,z) and are computed as
γ1
[xyz 1] =[xyz 1] [T], i=2,3,4. (1)
γi γ1
Once rotated, each gridded face is projected onto the
sphere. The polar caps are built in the same way, except
thatthestartingpointisthesquareplanedefinedwithinthe
region γ =[ π/4 λ π/4] [ π/4 λ π/4]. At this
9
− ≤ ≤ × − ≤ ≤
pointwehavealineargridmadeof6disjointregularpatches.
The top view of this is shown in Fig. 1. Eight new surfaces Figure 2.Detailofoneinterfacepatchinthenorthernhemisphere.
must now be built to fill the ungridded gaps.
Thehigh-orderLegendre-Gauss-Lobatto(LGL)pointsare
added to the linear grid by projecting the linear elements
ontotheauxiliarygnomonicspace(ontheplane),andback-
projecting the new high order quads onto the sphere. The
total number of elements N and grid points N of the high
e p
order conforming grid are given by:
N =
e
4N N Lateral Lat-Lon
elon elat (4)
+2N N Polar caps
elon elon
+8N N Transition zone,
elon elat,int
and
N =
p
4N p(N p+1) Lateral
elon elat (5)
+2(N p+1)(N p+1) Polar
elon elon
+2(4N p(N p+1) 8N p) Transition.
elon elat,int − elon
In (4) and (5), N is the number of elements of order N
Figure 1.Linear grid. Top x−y view of the disjoint lateral Lat-Lon elon
and top patches in the multi-block RLL grid. TFI is used to connect along the longitude direction of one face in the main Lat-
themtogetherthroughaninterfacesurfacegrid. Lon band that, by construction, coincides with the number
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
4 Marras et al.
of longitude elements in the interface and polar regions. 2.3. Quad-based icosahedral grid
N and N are, respectively, the number of latitude
elat elat,int
elements of the Lat-Lon region and of the interface patches. The quad-based icosahedral grid of order N is derived from
The number of latitude elements in the polar caps is also aninitialicosahedronof20triangles.Forbetterpropertiesof
given by N . thehighordergrid,everytriangleismappedontoagnomonic
elon
Fig. 3 shows an example of a high-order conforming RLL space by the transformation formulas (Giraldo (2001))
tessellation.
x=rtanλ′, (6a)
y=rtanϕ′secλ′, (6b)
where r is the radius of the sphere and, given the centroids
(λ ,ϕ ) of each triangle, we define
c c
cosϕsin(λ λ )
λ′ =atan − c , (7a)
cosϕ sinϕ+cosϕ cosϕcos(λ λ )
(cid:20) c c − c (cid:21)
ϕ′ =asin[cosϕ sinϕ sinϕ cosϕcos(λ λ )]. (7b)
c c c
− −
After mapping, the triangles are subdivided into smaller
ones by a Lagrange polynomial of order n .
I
The number of quadrilateral elements and grid points of
the final grid are then given by
N =6(NT 2)p2+2
p p −
Figure 3.ConformingRLLgrid. Ne =6(NpT −2)
where the number of points NT of the original triangular
p
2.2. Cubed sphere (or Hexahedral) grid can be found in Giraldo (2001).
Fig. 5 shows the 8th order icosahedral grid. As for the
Thecubedsphere(see,e.g.Ronchi et al.(1996);Taylor et al. previous cases, the elements are curved and lie on the
(1997, 2013)) is derived from the projection of a cube onto spherical surface.
the sphere. As such, it consists of 6 faces that are then
subdivided into as many elements as necessary. Like RLL,
we built this grid in a multi-block fashion by going through
the same steps described above. The fundamental difference
is that the hexahedral grid (Hex, from now on) has only 6
faces and they are all equal. An example is shown in Fig.
4. The internal LGL points are omitted for simplicity. By
construction, the elements that are furthest from the center
of each face are smaller than those at the face center; in
addition, their distorted shape prove to be a concern when
the resolution is not sufficiently high. We will discuss this
issue shortly.
Figure 5.ConformingIcogrid.
3. Shallow water equations on the sphere
3.1. Equations of motion
Under the assumption of a shallow depth, the motion of
an incompressible fluid may be described by the non-linear
shallowwaterequations.Theshallowwaterassumptionholds
aslongasthecharacteristichorizontalextensionofthefluid
is much larger than its depth. The vector representation
Figure 4.ConformingHexgrid. of the shallow-water equations in the Cartesian coordinates
x=(x,y,z) is given, in flux form, by
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
Grid effect on spherical shallow water jets 5
where F˘(q) indicates the numerical flux of the advection
∂φ
+ (φu)=0 (8a) and diffusion terms across the boundary Γ. In the case of
∂t ∇· spectral elements, where the basis functions are continuous
∂φu across neighboring elements, and given that the problem
+ (φu u)= φ φ f(x φu)+µx+δν 2(φu)
∂t ∇· ⊗ − ∇ − × ∇ is solved on a closed manifold, the boundary integral
Γ
(8b) vanishes. On the other hand, in the case of the DG method,
where φ=gh is the geopotential height (g and h are the theboundaryintegralremainsbecauseofthenon-zerointeRr-
modulus of the acceleration of gravity and the vertical element fluxes. Furthermore, the surface integrals for DG
height of the fluid), ν 2 is the artificial viscosity term of aredefinedelement-wiseoverΩ ratherthangloballyoverΩ.
viscouscoefficientν =1∇ 105m2s−1 thatisde-activatedby ThismakesEquations(11)ande(12)asystemofN equations
the switch δ=0, u=(u×,v,w)T is the velocity vector, = that are coupled via the flux of Equation (12). Boeth options
(∂x,∂y,∂z)T,andf =2ωz/r2 istheCoriolisparameterw∇ith are available within the same code used for this study and
Earth’s angular velocity ω=7.292 10−5s−1 and mean both methods share most of the numerical machinery.
×
radius r=6.37122 106m. For the fluid particles to remain
Theintegralsabovearesolvedelement-wiseonthediscrete
×
on the spherical surface defined in 3D cartesian space, the polyhedral approximation Ωh. The discrete domain Ωh is
normal component of the velocity field with respect to
defined by the union of N non-overlapping quadrilateral
e
the sphere is removed by the term µx in the momentum elements Ωh. For every element we define a bijective
e
equation, where µ is the Lagrange multiplier that serves
transformation : (x,y,z) (ξ,η,ζ) that maps the
this purpose. The construction of the Lagrange multiplier FΩhe →
physical coordinate system, (x,y,z) to the reference system
is described in Sec. 3.2. (ξ,η,ζ) and is such that the reference element Ωˆh(ξ,η)=
e
[ 1,1] [ 1,1] lies on the spherical surface. ζ is thus the
3.2. Forcing the fluid onto the sphere by a Lagrange − × −
spherical radius vector whose discrete values identify a
multiplier
spherical shell (e.g. ζ =1 is the sphere of unit radius.)
The Jacobian matrix of the transformation has com-
WeusedtheapproachofCoté(1988).Specifically,attheend
ponents Ji =∂ i and determinant J =∂ x G, where
ofeverytimestepn+1,themomentumoftheflowmustbe ξF | | ζ ·
G=∂ x ∂ x is the surface conforming component of the
correctedfortheparticlestoremainonthesphere.Ifcandu ξ × η
transformation (Song and Wolf 1999; Giraldo 2001).
denote the constrained and unconstrained values of φu, the
Giventhedefinitionofthereferenceelement,theelement-
new momentum is corrected according to
wise solution can be approximated by the expansion
(φu)n+1 =(φu)n+1+µx. (9)
c u
For the fluid to remain on the sphere, the condition (N+1)2
u·x=0 must hold (i.e., the velocity field and the position qN(x,t)|Ωe = ψk(FΩ−e1(x))qN(xk,t), e=1,...,Ne
ivnecetqouramtiounst(b9)e.oFrotrhtohgiosntaol)b,ehternucee,itmhepvlyailnugetohfattheφmucu·ltxip=lie0r Xk=1 (13)
is found to be µ= φun+1 x/r2. where(N +1)2isthenumberofcollocationpointswithinthe
TheGalerkindis−cretizuatio·n of(8)isdescribed inthenext elementoforderN,andψk aretheinterpolationpolynomials
section. evaluatedatpointk.Thebasisfunctionsψk areconstructed
as the tensor product of the one-dimensional functions
4. Numerical method hi(ξ(x)) and hj(η(x)) as:
4.1. Spectral elements and discontinuous Galerkin ψ =h (ξ(x)) h (η(x)),
k i j
⊗
approximation
i,j =0,...,N.h (ξ(x))andh (η(x))arethebasisfunctions
i j
∀
To solve the shallow water equations by element-based associated with the N LGL points ξi and ηj, respectively,
Galerkin methods on a domain Ω, we proceed by defining given by the roots of
the weak form of (8) that we first write in compact form as
(1 ξ2)P′ (ξ)=0,
∂q − N
+ F(q)=S(q), (10)
∂t ∇· where P′ (ξ) is the derivative of the Nth-order Legendre
N
whereq=[φ,φu]T isthearrayofthesolutionvariables,and polynomial. Given these definitions, the one-dimensional
Lagrange polynomials h (ξ) are
F and S are the flux and source terms that can be easily i
identified in the equations (8) above. Given a basis function
1 (1 ξ2)P′ (ξ)
ψthatbelongstotheSobolevspaceH1inthecaseofspectral h (ξ)= − N .
i
elements,andtoL2 inthecaseofDG,theweakformof(10) −N(N +1)(ξ ξi)PN(ξi)
−
is given by the integral
The functions h (η) are computed in the same way.
j
The same expansion above is used to construct the
∂q
ψ + F(q) S(q) dΩ=0, (11) derivatives. We write:
∂t ∇· −
ZΩ (cid:20) (cid:21)
that, after integration by parts of the flux, becomes
(N+1)2
∂q (x,t) ∂q (x ,t)
N = ψ ( −1(x)) N k ,
∂q ∂t |Ωe k FΩe ∂t
ψ dΩ ψ F(q)dΩ ψS(q)dΩ+ k=1
X
ZΩ ∂t −ZΩ∇ · −ψnZΩF˘(q)dΓ=0, (12) ∂qN∂(xx,t)|Ωe = (N+1)2 ∂ψk(F∂Ωx−e1(x))qN(xk,t). (14)
·
ZΓ Xk=1
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
6 Marras et al.
The integrals are computed by numerical quadrature on where R(q) indicates the right-hand sides of equations
the reference element as follows: (16,17), K is the number of stages of the RK method, and
q0 =qn,qK =qn+1. All the computations were executed
with a maximum Courant-Friedrichs-Lewy number (CFL,
q(x,t)dx= q(ξ,t)J(ξ)dξ Courant et al.(1967))approximately0.5.Thisvalueensures
ZΩhe ZΩˆhe | | ≈ the stability of the scheme given that the time step ∆t is
Nlgl such that
q(ξ ,η ,t)J(ξ ,η )ωξωη, (15)
i,j=1 i j | i j | i j ∆t CFL min[u χ +φ√χ χ]−1.
X ≤ x∈Ωh | · | ·
wheretheGaussianquadratureweightsω arecomputedas
i,j
The local grid distortion, χ, is defined as the vector
2
2 1
ωi,j = N(N +1)(cid:18)PN(ξi,ηj)(cid:19) . χ= |∆ξxξ| + |∆ηxη|,|∆ξyξ| + |∆ηyη|,|∆ξzξ| + |∆ηzη|
(cid:18) (cid:19)
In the case of spectral elements, the substitution of the
expansions (13,14) into the weak form (12) yields the semi- where ∆ξ and ∆η are the local average grid size across the
discrete (in space) matrix problem surface.
Intheabsenceofdissipation(e.g.,artificialdiffusion),the
∂q =DTF(q)+S(q) (16) high order numerical solution of advective problems may be
∂t characterizedbyspurioushighfrequencymodes.Topreserve
where, for the global mass and differentiation matrices, M full stability, along with respecting the CFL condition, we
b
and D, we have that D=M−1D. The global matrices applied the standard spatial Boyd-Vandeven filter (Boyd
are obtained from their local counterparts, Me and De, 1996) at every time step. The filter constant was chosen to
by means of the direct sbtiffness summation operation that reduce by 5% the highest modes only.
maps the local degrees of freedom of an element Ωh to the
e
corresponding global degrees of freedom in Ωh, and adds 5. Tests
theelementvaluesintheglobalsystem.Byconstruction,M
is diagonal (assuming inexact integration), with an obvious To evaluate the behaviour of the solution on the three
advantage if explicit time integration is used. grid families, the highly nonlinear zonal dynamics test by
ConcerningthediscontinuousGalerkinapproximation,the Galewsky et al.(2004)wasrunatdifferentresolutionsusing
problem at hand is solved only locally, and unlike the case elementsoforder8.BasedontheanalysisbyGaberšek et al.
of spectral elements, the flux integral in Equation (12) must (2012),8istheorderthatwearelikelytousewhenrunning
be discretized as well. The element-wise counterpart of the operational (weather prediction) simulations with NUMA.
matrix problem (16) is then written as: The test consists of (i) an unperturbed mid-latitude zonal
jet on top of a geostrophically balanced geopotential height,
andof(ii)itsperturbedcounterpartwheretheperturbation
∂qe
= (Ms,e)TF˘(qe)+(De)TF(qe)+S(qe), (17) istriggeredbyabumpinthegeopotential.Ineithercase,the
∂t −
initialzonalvelocityandheightareafunctionofthelatitude
where we obtcain Ms,e =(Meb)−1Ms,e from the element ϕ (measured in radians unless otherwise specified) as
boundary matrix, Ms,e, and the element mass matrix,
Me. Out of varioucs possible choices for the definition of
0 if ϕ ϕ or ϕ ϕ
the numerical flux F˘(q) in Equation (17), we selected the u(ϕ)= ≤ 0 ≥ 1
Rusanov flux that, for the inviscid part of the flux, is ( umenaxexp (ϕ−ϕ0)1(ϕ−ϕ1) if ϕ0 <ϕ<ϕ1,
constructed as h i (19)
where u =80ms−1 is the maximum zonal velocity
1 max
F˘(q)= 2 FN(qNL)+FN(qNR)−|λ|(qNR−qNL)n , (18) and (ϕ0,ϕ1)=(π/7,π/2−ϕ0) are the lowest and highest
latitudes that delimit the jet. At the jet mid-point, where
(cid:2) (cid:3)
where λ = n u +√φ is the maximum wave speed of ϕ=π/4,thenon-dimensionalparameteren =exp[ 4/(ϕ1
the sha|llo|w |wa·ter| equations. By identifying an intrinsic ϕ0)2]normalizesthejetmagnitudetoumax.From−theiniti−al
tangentialdirectionineachinter-elementboundaryedge,qL zonal flow given by (19), the height gφ is found by solving
N
and qR denote the solution in the left (L) and right (R)
N
elements, respectively. For a detailed description of the DG
ϕ tan(ϕ′)
algorithmthatwasusedforthisstudy,thereaderisreferred gφ(ϕ)=gφ ru(ϕ′) f +u(ϕ′) dϕ′, (20)
to Kopera and Giraldo (2013). 0− r
Z (cid:20) (cid:21)
4.2. Time integration
where the integral on the right-hand side is computed by
numerical quadrature. As in the reference, ϕ is chosen
The solution of the systems of ordinary differential 0
to give a global mean layer depth of 10 km. The height
equations (16,17) is advanced in time by a third-order
perturbation of test (ii) is the bump whose shape is given
strong stability preserving Runge-Kutta method (SSP-
by
RK3) (Cockburn and Shu 2001; Spiteri and Ruuth 2002)
that yields the solution at step n+1 as
φ′(λ,ϕ)=φˆcos(ϕ)exp[ λ2/α2]exp[ (ϕ ϕ )2/β2]
2 1
− − −
qk =αkqn+αkqk−1+βk∆tR(qk−1),fork=1,...,K, for π<λ<π, (21)
0 1 −
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
Grid effect on spherical shallow water jets 7
Table 1.Equilibrium. Normalized L2 error of (φ,us,vs) at day 6
where λ is the longitude, φˆ=120m, ϕ =π/4, α=1/3 using CG. The error is computed with respect to the analytical
2 initialfieldsgivenbyequations(19)and(20),andisreportedagainst
and β =1/15. The definition of this bump is such that the
thetotalnumberofgridpoints.
perturbation is forced to zero at the poles.
The flows in (i) and (ii) are simulated along a time CG, L
2
span of 144 hours. In the absence of any perturbation, the N points 100000 50000 25000 12000
balanced field and zonal flow are expected to maintain their
RLL
initial state indefinitely. The analytical initial field is used
φ 4.878e-8 1.802e-7 2.011e-5 1.466e-4
tocomputetheL errornormsofthenumericalsolutionsas
2 u 5.910e-5 1.586e-4 5.162e-3 2.701e-2
s
time evolves.
v 5.910e-5 1.586e-4 5.167e-3 2.706e-2
s
Although its solution seems trivial, this test is especially
Hex
demanding for a shallow water model because of the large
φ 6.099e-5 4.896e-4 5.440e-3 5.460e-3
sensitivityofthesolutiontothegridresolutionandgeometry.
u 6.679e-3 5.346e-2 5.458e-1 5.059e-1
s
Unlikeotherteststhatmaybemoreforgivinginthisrespect
v 6.679e-3 5.346e-2 5.458e-1 5.059e-1
s
(e.g. the steady-state nonlinear zonal geostrophic flow by
Ico
Williamson et al. (1992)), when the grid is not sufficiently
φ 2.230e-5 2.593e-4 1.057e-2 1.263e-2
fine, grid related oscillations may spoil the equilibrium with
u 2.040e-3 2.313e-2 8.035e-1 9.163e-1
s
devastating consequences. This is clearly visible in Fig. 6,
v 2.040e-3 2.313e-2 8.035e-1 9.163e-1
s
where the initially unperturbed and balanced zonal flow
completely loses its initial state when the test is executed
Table 2.Equilibrium. As table 1, but the error is reported against
on the low resolution hexahedral and icosahedral grids. All
theequivalentresolution∆λ.
the runs were executed at the equivalent resolutions ∆λ
∆ϕ= 0.625◦,1.25◦,2.5◦,5.0◦ and for a varying number o≈f
CG, L
{ } 2
grid points Np = 12000,25000,50000,100000 . The reason RLL 0.625◦ 1.25◦ 2.5◦ 5.0◦
{ }
for the two sets of runs will be clarified in Remark 1.
N points 180000 49000 16000 3000
Below, the keywords Equilibrium and Perturbed jet are
φ 1.377e-8 5.017e-7 1.106e-4 2.056e-3
used to identify, respectively, problems (i) and (ii).
u 2.185e-5 3.167e-4 2.553e-2 4.694e-1
s
v 2.185e-5 3.167e-4 2.555e-2 4.694e-1
s
Remark 1: equivalent resolutions A few words are in Hex 0.625◦ 1.25◦ 2.5◦ 5.0◦
order regarding the issue of equivalent resolution when grids N points 124000 31000 9600 1500
of a very different nature are compared against each other. φ 5.109e-5 4.721e-3 6.336e-3 1.038e-2
InthecaseofspectralelementsthatrelyonLGLnodes,the u 5.601e-3 4.709e-1 4.832e-1 6.918e-1
s
distancebetweentwoconsecutivegridpointschangeswithin v 5.601e-3 4.709e-1 4.832e-1 6.918e-1
s
the same element. Because of this, the equivalent angular
resolution along a latitude circle subdivided into N
e,λ
elements of order N is 2π/(Ne,λN). Given this definition, that are not aligned with the characteristic flow (which is
we add that two resolutions are approximately equivalent likely to be the case in realistic simulations), we should not
when, within the region of major interest (e.g., where the trust the solution with high fidelity.
dynamics is most likely to develop from a certain initial In Tables 1 and 2 we report the normalized L error
2
state), the size of the angular resolution of two different norms of the geopotential height and zonal velocity against
grids are comparable. In other words, when we say that the number of grid points and, respectively, the equivalent
the resolution is, e.g., ∆λ=0.625◦, it means that the mean resolution ∆λ. These errors were computed using spectral
angular distance between two meridians in the proximity elements on conforming grids; from here on, we shall
of the jet measures 0.625◦. This definition makes the refer to spectral elements as Continuous Galerkin (CG) to
measurement on the icosahedral grid not straightforward distinguishitfromthediscontinuousGalerkin(DG)method.
because of the irregular orientation of the elements. To Wementionedabovethatthestructureoftheicosahedral
obviate this problem, the error norms in the analysis below gridissomewhatirregularanditishencenotstraightforward
are also compared with respect to the total number of grid to find a direct relationship between the number of grid
points. points and the equivalent resolution. For this reason, we
report the errors computed on the Ico grid in Table 3
againsttheequivalentresolutionsandagainstthesubdivision
5.1. Equilibrium parameters [N6N5N4N3N2N1] that are defined by the
user to construct the grid. We plot the errors of Tables 1-3
Fig. 6 shows the vorticity field Λ after 6 days for the in Figs. 7 and 8.
unperturbed jet problem on the three different conforming Let us focus on the curves with continuous lines as
grid types (rows) and N = 25000,50000,100000 grid they trace the error on the uniform conforming grids. The
p
{ }
points (from left to right). For two out of the three grid conclusions drawn from the vorticity contours of Fig. 6 are
geometries, the effect of the misalignment is still important confirmed by these error estimates. Furthermore, from Fig.
for as many as 25000 grid points, which approximately 8itisinterestingtoseethattheL erroroftheHexsolution
2
corresponds to an equivalent resolution ∆ϕ=1◦ on the reaches its peak value at ∆λ 1.25◦. This leaves very little
≈
cubed-sphere. The solution on the RLL grid can maintain room for resolution options when the cubed-sphere is used
the jet with approximately 80% fewer grid points than withoutviscosity.Theresultsreportedsofarwerecomputed
the hexahedral grid and 50% fewer than the icosahedral withCGonconforminggrids.However,wealsoranthesame
grid. This is not reason enough to draw conclusions on the tests using DG on every grid. As expected, the errors are
superiority of the RLL grid; rather, it should only suggest practically identical to those obtained by CG; the solution
that, unless a sufficiently high resolution is used with grids accuracy for CG and DG should be identical because the
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
8 Marras et al.
Hex 25k Hex 50k Hex 100k
Rll 25k Rll 50k Rll 100k
Ico 25k Ico 50k Ico 100k
Figure 6.Equilibrium: Vorticity field (×10−5 s−1) after 6 days. The solutions on the hexahedral (top row), on the RLL (middle), and on the
icosahedralgrids(bottomrow)areplottedfor25000,50000,and100000gridpoints(fromlefttoright).Thecolorscalerangesbetween−12×10−5
s−1 (darkblue)and16×10−5 s−1.
Table 3.Equilibrium. As table 2, but on the icosahedron. Because of the geometrical structure of the icosahedron, the correspondence
between the equivalent resolution and the number of subdivision indicated by N6,...,N1 is not as straightforward as it is for Hex and
RLL.
Geometry N6 N5 N4 N3 N2 N1
Approximate ∆λ 0.625◦ 1.25◦ 2.5◦ 5.0◦
φ 5.064e-6 2.299e-5 2.593e-4 6.157e-3 1.261e-2 1.033e-2
u 4.381e-4 2.039e-3 2.313e-2 5.169e-1 9.166e-1 7.774e-1
s
v 4.381e-4 2.039e-3 2.313e-2 5.169e-1 9.166e-1 7.774e-1
s
solutionisrelativelysmooth(differencesshouldonlyemerge Table4.Equilibrium.Asintable2,butusingDG.
betweenthetwomethodsinthepresenceofdiscontinuities).
For comparison, we report the values of the DG normalized DG, L2
error norms in Table 4. We will show DG results on non- RLL 0.625◦ 1.25◦ 2.5◦ 5.0◦
conforming grids with static refinement in Sec. 5.3. N points 180000 49000 16000 3000
φ 1.384e-8 6.556e-7 1.203e-4 2.151e-3
InFig.9weplotthetimeevolutionof φ during6days
k kL2 u 2.189e-5 3.261e-4 2.636e-2 4.794e-1
for 100000 point grids. After the gravity wave adjustment s
v 2.189e-5 3.260e-4 2.640e-2 4.794e-1
during the first 6 hours (see Galewsky et al. (2004)), the s
Hex 0.625◦ 1.25◦ 2.5◦ 5.0◦
regularityofthegridgeometrystillplaysanimportantrolein
N points 124000 31000 9600 1500
theevolutionoftheerrorevenataresolutionthat,inglobal
φ 6.649e-5 4.948e-3 6.397e-3 1.540e-2
mode, is typically considered fairly high in an operational
u 7.291e-3 4.933e-1 4.972e-1 9.369e-1
setting. s
v 7.291e-3 4.933e-1 4.972e-1 9.369e-1
Sofar,theviscoustermintheequationshasbeenignored s
(i.e. δ=0). In the quest of diminishing the negative effects
of the low resolution on certain grids, we turned viscosity
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls
Grid effect on spherical shallow water jets 9
10-1 101
RLL RLL
RLL-1L RLL-1L
10-2 RLL-2L 100 RLL-2L
Hex Hex
10-3
ror Hex-1L ror 10-1 Hex-1L
r r
e e
L 2 10-4 Hex-2L L 2 Hex-2L
d Ico d Ico
ze ze 10-2
mali 10-5 Ico-1L mali Ico-1L
r Ico-2L r Ico-2L
o o
Φ n 10-6 u ns 10-3
10-7 10-4
10-8103 104 105 106 10-5103 104 105 106
Total number of points Total number of points
Figure 7.Equilibrium at day 6: Normalized L2 error norm of φ (left) and us (right) against the number of grid points. The error is computed
withrespecttotheinitialconditionofthebalancedproblemoneverygrid.
10-1 101
10-2
100
10-3
ror ror 10-1
er 10-4 RLL er RLL
2 2
L L
RLL-1L RLL-1L
d d
ze 10-5 RLL-2L ze 10-2 RLL-2L
ali ali
m Hex m Hex
or 10-6 or
Φ n Hex-1L u ns 10-3 Hex-1L
10-7 Hex-2L Hex-2L
Ico Ico
10-4
10-8 Ico-1L Ico-1L
Ico-2L Ico-2L
10-9 10-5
0.625 1.25 2.5 5.0 0.625 1.25 2.5 5.0
Resolution (Degrees) Resolution (Degrees)
Figure 8.Equilibriumatday6:NormalizedL2 errornormofφ(left)andus (right)againsttheresolutionmeasuredindegreesfortheRLLand
Hexgrids.Theerroriscomputedwithrespecttotheinitialconditionofthebalancedproblemoneverygrid.
on and recomputed the normalized L error norms that that was well behaved in the inviscid case has deteriorated.
2
we report in Table 5. We compare these errors against the This behavior is expected for two main reasons: (1) the
inviscid estimates using 12000 and 25000 grid points and viscous and inviscid equations are, mathematically, not the
point out the instances when artificial viscosity was either samesetofequations,andphysicallyrepresenttwodifferent
beneficialorharmfultotheerrormeasure.Tointerpretthese dynamics;(2)viscosityissolargethatitisdampingallofthe
values, we need to take a parallel look at the contour plots modestriggeredbythegridontheonehandbut,inthesame
of Figures 10 and 11. In the case of RLL, where the grid is way,dampsthesolutionalsowhereitshouldnot.Athorough
alignedwiththecharacteristicflow,itwasshownabovethat analysis of the diffusion mechanisms that are currently used
the resolution does not significantly affect the solution. The inatmosphericproblemsgoesbeyondthescopeofthispaper.
opposite, however, occurs in the case of the other two grids. Nevertheless,itshouldbeemphasizedthatusinganarbitrary
By adding viscosity, what was a bad inviscid solution seems diffusion operator may not actually improve the solution
to improve (see the flow pattern after 6 days in Figs. 10, 11, quality even though the final solution may appear smooth,
and the error measure in Table 5); in contrast, the solution as we have shown for the RLL grid. This simple result
(cid:13)c 2013RoyalMeteorologicalSociety Preparedusingqjrms4.cls