Application of a Galerkin Finite Element Scheme to ...

AIAA 2010-690

48th AIAA Aerospace Sciences Meeting
Including the
New Horizons Forum and Aerospace 4 - 7 January 2010, Orlando, Florida

Application of a Galerkin Finite Element Scheme to Atmospheric Buoyant and Gravity Driven Flows Simone Marras∗ , Mariano Vazquez† , Oriol Jorba‡ Romain Aubry† , Guillaume Houzeaux† and Jos´e M. Baldasano§ Barcelona Supercomputing Center BSC-CNS, c/ J. Girona 29, 08034 - Spain

The application of a new finite element (FE) technique for the solution of stratified, non-hydrostatic, low-Mach number flows is introduced in the context of mesoscale atmospheric modeling. In this framework, a Compressible Variational Multiscale (VMS-C) finite element algorithm to solve the conservative form of the Euler equations coupled to the conservation of potential temperature was developed. This methodology is new in the fields of Computational Fluid Dynamics for compressible flows and in Numerical Weather Prediction (NWP), and we mean to show its ability to maintain stability in the solution of thermal, gravity-driven flows in a stratified environment. This effort is justified by the advantages offered by a Galerkin finite element algorithm when massive parallel efficiency is a constraint, which is indeed becoming the paradigm for both CFD and NWP practitioners. The algorithm is validated against the standard test cases specifically designed to test the dynamical core of new atmospheric models. In the context of buoyant and gravity flows three tests are selected among those presented in the literature: the warm rising smooth anomaly, and two versions of the density current evolution from a cold disturbance defined by different initial conditions. The reference quantitative and qualitative values are taken from the literature and from the output obtained with the Weather Research and Forecasting model (WRF-ARW), a state-of-the-art research NWP model. Keywords: Numerical Weather Prediction, Low Mach, Compressible Flows, High Performance Computing

I.

Introduction

n the context of Numerical Weather Prediction (NWP), the new high-performance computing (HPC) Igroups paradigm based on the use of thousands of cheap processors is driving the majority of NWP research to review and modify the numerical methods their models are based on . 1–5

The ability of taking full advantage of such powerful computer architectures is linearly dependent on the numerics behind the discretized mathematical model. The limitations of the commonly used finite difference schemes (FD) or those of global methods that employ spectral transforms (ST)6 (this method has dominated global NWP and climate models for the past two decades7 ), appear clearly when efficiency on massively parallel machines is the final goal (e.g., large communication footprint of the parallel finite differences is a major limitating factor). These limitations can be circumvented by means of more global methods such as the finite volume technique (FV), or by the local domain decomposition property of finite element methods (FE). In the past two decades this direction has been taken extensively by many after the pioneering work of Ref. 8, 9, and others10–13 in the 60s, 70s, and 80s, who studied variational formulations for the treatment of stratified flows. The recent developments of Ref. 14 and 15 within the operational global spectral model of the European Centre for Medium-Range Weather Forecasts (ECMWF) and the finite difference Met Office Unified Model (UM) respectively, are examples of the successful application of a FE approach to the ∗ Ph.D.

Candidate, BSC-CNS and Universitat Politecnica de Catalunya. [email protected] - AIAA Student Member Research Scientist, Dept. of Computer Applications, BSC-CNS. Contact: [email protected] ‡ Senior Research Scientist, Dept. of Earth Sciences, BSC-CNS. Contact: [email protected] § Head of the Earth Sciences Department, BSC-CNS

† Senior

1 of 10

Institute Aeronautics andPublished Astronautics Copyright © 2010 by Simone Marras, Mariano Vazquez, OriolAmerican Jorba, Romain Aubry, of Guillaume Houzeaux. by the American Institute of Aeronautics and Astronautics, Inc., with permi

vertical discretization of two operational and research NWP models; moreover, the Canadian Environmental Multiscale model (GEM)16–18 , based on the efforts of Ref. 19 and 20, was indeed developed on a FE basis with the parallel efficiency in mind16 . From the field of computational fluid dynamics (CFD) other schemes began to appear as applied to NWP problems. Extensive work on the discontinuous Galerkin (DG) or Spectral Element techniques for the solution of the Euler and Navier-Stokes equations comes, for example, with 21–25, who applied CFD stabilization methodologies for the treatment of stratified atmospheres in the very low Mach regime. Similarly, and with a mind set to unstructured, dynamically adaptive meshes, the studies on FV models such as the NWP software OMEGA,26 and the NICAM,27, 28 are worth mentioning. With more interaction between NWP scientists and CFD developers, the two fields may take full advantage of each other achievements when stability and accuracy in the solution of low Mach, stratified flows are the main concern. In this article we introduce the use of a Compressible Variational Multiscale (VMS-C) finite element method to yield stable solutions in the context of low-Mach atmospheric flows, focusing on the performance of this technique when run on linear continuous quadrilateral elements (Q1). The algorithm is validated against the standard test cases specifically designed to test the dynamical core of new atmospheric models.29 In the context of buoyant and gravity flows we selected three tests among those presented by other authors: the warm rising smooth anomaly,30 the density current test,31 and a second version of density current (density current II) that comes from the analysis shown in Ref. 32 and that is useful for a more quantitative analysis. After this introduction, the mathematical problem and its numerical formulation are described in section II; the validating results are given in section III, and the algorithm parallel performance in section IV. Conclusions follow in section V.

II.

Governing Equations

The solution of a particular set of governing equations used for a dry, stratified, non-rotating atmosphere is described in this section; this set is a specific manipulation of the compressible Euler equations where energy is formulated in terms of potential temperature. The reasons that justify this particular choice are here given, and the formulation that is most commonly adopted in atmospheric models is also shown. Firstly, given that all physical processes that drive the dynamics of dry atmospheres can be considered to be adiabatic, viscous effects are depreciable, and hence the main reason for using an inviscid scheme. Secondly, in spite of the low peak Mach numbers characteristic of atmospheric motion, the stratification is such that all thermodynamic quantities vary as a function of pressure along the vertical coordinate. When the energy equation is formulated in terms of potential temperature per unit volume (Θ = θρ), the governing set is expressed as follows: ∂Φ ∂Fi (Φ) ∂Li (Φ) + + = Si , ∂t ∂xi ∂xi

(1)

where Φ = (ρ, Uj , ρ θ) is the vector of the conserved variables (density, linear momentum, and potential temperature per unit volume). The array of the fluxes Fi , the pressure term Li , and the source term Si read:       Ui 0 0       Fi =  ui Uj  , Li =  δij P  , Si =  −ρg  0 ui ρ θ 0 with the equation of state expressed by:

p = p0

'

Rρθ p0

( ccp v

.

(2)

Above, g = 9.81 m s−2 is the acceleration of gravity, cp = 1004 Jkg −1 K −1 and cv = 717 Jkg −1 K −1 are the coefficients of specific heat at constant air pressure and volume respectively, R = 287 J K −1 kg −1 the gas constant, and po = 105 P a a reference pressure. Because the equations are written on a 2D, non-rotating system of reference, Coriolis effects are not accounted for. Table 1 displays the characteristics of some meteorological models based on this set of equations. 2 of 10 American Institute of Aeronautics and Astronautics

Table 1. Existing atmospheric operational and research NWP models based on set (1).

Origin

Type†

UK, Met. Office

NH, mesoscale regional and global

Fully compr., FD semi-impl., semi-Lagrang. terrain following

Germany, COSMO Proj.

NH, mesoscale regional

Fully compr., FD RK split-expl, semi-impl. terrain-following Schar coords. is an option

USA, NCAR and collaborations

NH, mesoscale regional

Fully compr., FD semi-Impl., split Expl. Terrain-foll.

USA

NH, mesoscale

Fully compr., FV

Model Unified Model (UM)33

Lokal Model (LM)34

WRF-ARW35

Ahmad solver36*

Equations and Numerics•

* This

is not a NWP model. However, given that the Ahmad model solves the same set of equations within the same environment, this solver is worth mentioning. † NH stands for non-hydrostatic. • RK for Runge-Kutta; FD for Finite Differences; FV for Finite Volumes.

II.A.

Discrete Formulation

The numerical technique that we developed discretizes space in a Finite Element manner, and uses a Finite Difference scheme for time discretization. The projection on a FE space W of problem (1) yields the following compact discrete equations: ) ) ) ∂Φ ∂Li ∂Fi Ψ dΩ = Ψ dΩ + Ψ dΩ ∀Ψ ∈ W (3) ∂t ∂xi ∂xi which after integration by parts of some of the fluxes (e.g. F ), yields ) ) ) ) ∂ ∂Ψ ∂Li dΩ − Fi dΩ + ΨFi ni d∂Ω ΨΦdΩ = Ψ ∂t ∂xi ∂xi

∀Ψ ∈ W.

(4)

The last of the right hand side terms can be used to impose Neumann-like boundary conditions on the fluxes, with ∂Ω being the domain boundary and ni its exterior normal unit vector. The interpolation space where the variational form solution is to be found is (practically) the same as W , explicitly time independent. For that reason, the time derivative is taken out of the space integral, with the additional hypothesis that Ω does not change with time. The discretization is done by chosing compact support functions as the interpolators, defined from tesselations of the spacial domain. In this paper, we introduce the Variational Multiscale Method for Compressible flows (VMS-C)37 to mesoscale meteorological problems. The VMS stabilization technique was pioneered by Hughes38 and is based on the explicit resolution of the coarser scales while modeling the finer. The results here presented are still preliminary, although very promising, showing the accuracy of the VMS-C scheme for very low-Mach number flows. The stabilization terms of the VMS-C approach are of the form: S α = L∗ (Ψ) τ L(Φ), and are added to every equation of (4), where L(·) and L∗ (·) are respectively the space operator coming from (4) and its adjoint. τ is the stabilization matrix, which, in this case, is a diagonal matrix whose coefficients are the critical time steps for each equation.

III.

Validating Benchmarks

The VMS-C FE algorithm is tested in the context of atmospheric dynamics by means of standard benchmarks accepted by NWP practitioners. The ability to reproduce thermal transport and buoyant effects is 3 of 10 American Institute of Aeronautics and Astronautics

tested by imposing a thermal perturbation of proper definition within a base atmosphere in initial thermal and hydrostatic equilibrium. In all tests presented here, an initial elliptical pool of warm or cold air is introduced in a uniform potential temperature field, and allowed to evolve by means of buoyant and gravity effects. Results are now presented for (1) the warm rising bubble test of Ahmad et al.36 after Wiker and Skamarock,30 (2) the density current test of Straka et al.,31 and (3) the version of the density current described by Haertel, Johnson and Tulich.32 III.A.

Initial Field

The initial atmosphere is an ideal adiabatic environment in hydrostatic equilibrium and defined by a set of distributions for ρ, P, π, and θ, whose evolution is a combination of the mean value and its perturbation ∆ ¯ + ∆θ(x, t). as: ρ(x, t) = ρ¯(z) + ∆ρ(x, t), P (x, t) = P¯ (z) + ∆P (x, t), π(x, t) = π ¯ (z) + ∆π(x, t), θ(x, t) = θ(z) Every overlined quantity is in hydrostatic equilibrium throughout the run and does not change with time; the perturbation is the cause of the deviation from the initial equilibrium of every thermodynamic field.

III.B.

Smooth Warm Bubble

The test originally defined by Robert39 is here reproduced within the same domain [0.0; 20.0]×[0.0; 10.0] km2 as in Wiker and Skamarock.30 The environment is set to a uniform potential temperature θ¯ = 300 K with no mean flow and in initial hydrostatic equilibrium, and is perturbed by a smooth, linear distribution as in Ref.:36 * 0 if R > ro ¯ θ=θ+ (5) ∆Θ [1.0 − R/ro ] if R ≤ ro Above, ∆Θ = 2 K is the + intensity of the circular anomaly centered in (xc , zc ) = (10.0, 2.0) km; ro = 2.0 km is its radius, and R = (x − xc )2 + (z − zc )2 its analytical outline. The warm pool of air triggers an upward motion whose evolution is represented in figure 1. Initial and boundary conditions: the initial velocity is set to zero everywhere, and solid walls are imposed at the top, bottom and lateral boundaries. Results are shown in figures 1 and 2. Runs were made on structured grids of quadrangular elements and uniform resolution of 50 m. No spurious oscillations were observed despite the fully inviscid solution. The qualitative and quantitative comparison of the results displayed in figure 1 against those presented by Ahmad et al.36 is encouraging. III.C.

Density Current I

The initial base atmosphere is characterized by a uniform θ¯ = 300 K within a domain 50.0 km wide and 10.0 km deep. The field is perturbed by a cosine distribution of θ¯ centered in (xc , zc ) = (25.0 km, 3.0 km), with radius ro = (4.0 km, 2.0 km) and amplitude ∆Θ = −15 K: * 0 , - if r > 1 θ = θ¯ + (6) 1+cos(πc R) ∆Θ if r ≤ 1 2 ./ 02 / 02 x−xc c R= + z−z is the analytical boundary of the cold pool. xr zr Initial and boundary conditions: the initial velocity is set to zero everywhere; the domain lateral boundaries are open to the flow (pressure is imposed), while the top and bottom walls are solid and do not allow normal flow. The heavy pool of air gains a downward motion that causes the development of the propagating front. Once the front propagation is triggered, inertia is such that the top layers of the front develop into the Kelvin-Helmholtz vortical structures shown in figure 3. As it appears in table 2, the algorithm presents a smaller diffusivity at the finer computational mesh and all solutions agree in the value of the front speed. A good qualitative analogy appears between the profiles that we present and those of reference (e.g. Ref. 31). 4 of 10 American Institute of Aeronautics and Astronautics

Figure 1. Evolution of the thermal perturbations. From bottom left counterclockwise: initial distribution, after 300 s, after 600 s, and at the final step after 1020 s. The plot is a detail extracted from the full domain that extends within [0.0 : 20.0] × [0.0 : 10.0] km2 . A mesh resolution of 50 m in x and z was used.

Table 2. Density Current I: comparative results of potential temperature and front location obtained with different models (and schemes) at 900s.

Model VMS 50 m VMS 100 m VMS 150 m WRF-ARW 50 m (FD) WRF-ARW 100 m (FD) WRF-ARW 150 m (FD) DG22 f-wave (FV)36

III.D.

∆θmin -15.03 -8.204 -7.575 -7.58 -7.46 -7.35 -8.90 -9.82

∆θmax 0.00051 0.00058 0.00022 0.0 0.0 0.0 0.00012 0.00892

Front Location [m] 15352 15285 15508 15500 15800 15650 14767 14975

Density Current II

Similarly to that described in the previous paragraph, here we consider the characterization of a gravity current as reported in Ref. 32. Haertel, Johnson and Tulich simulate an idealized thunderstorm outflow initialized as an instantaneous cooling of the lower layers of an isentropic atmosphere (a description of the physics behind this phenomenon is given in their paper). The difference that this case has with respect to the one presented in Ref. 40 and Ref. 41 is in the geometrical definition of the cold region: instead of the linear stratification or cosine-law distribution of the potential temperature as in the Droegemeier and Straka cases 5 of 10 American Institute of Aeronautics and Astronautics

Figure 2. x-velocity (top) and y-velocity (bottom) after 1020 seconds for the evolution of the thermal bubble of figure 1.

respectively, the cold pool is here given by a uniform potential temperature θ = 297 K within a semi-elliptic region centered in the middle of the domain and extending over 10 km in the x, and 2 km in the z direction. The outer atmosphere is isentropic with reference value of θ = 300 K. The simple geometrical shape that develops upstream of the moving leading edge justifies the introduction of this second type of density current. In fact, the resulting shape makes it is easy to verify the kinematic relation that Benjamin42 gave between the velocity of the leading edge c and the depth d of the current a few meters upstream. If Density Current I provides important information on the quality of the simulation when compared to other models, Benjamin’s law can be the bridge towards a quantitative evaluation of the test. Benjamin’s equation is the following: + c = k g$ d , (7)

where g’ is the product of the acceleration of gravity and the ratio of the perturbation density ∆ρ$ to the environment ρ, and k is a function of the ratio of the depth of the + density current to the depth of the domain in which it is advancing; being this ratio defined by α, k = (α2 − 3α + 2)/(1 + α). In physical thunderstorm outflows values between 0.71 and √ 1.25 are usually observed for k (e.g. Ref. 43,44). The largest k-value theoretically obtainable is limited to 2 when α → 0, decreasing for increasing values of α. At equal geometry, in Ref. 45 it is observed that theoretical studies and physical outflows can be characterized by different values of k partly because of surface friction and mixing effects that may differ in the two cases. Figure 5 shows the evolution of the simulated outflow. If we take the largest absolute value of the perturbation density ∆ρ$ = ρ − ρ$ inside the perturbed flow, and coherently, the largest value of density within the atmosphere, the entries displayed in table 3 are thus obtained for Eq.(7): giving a theoretical leading edge speed of c = 14.05 ms−1 . The difference between the 6 of 10 American Institute of Aeronautics and Astronautics

(a) t = 300 s

(b) t = 600 s

(c) t = 900 s Figure 3. Density current I at different time steps and resolutions: structured grid of quadrangular elements with dx = dz = 50 m (left columns); with dx = dz = 100 m (center); and with dx = dz = 150 m (right column). All runs are inviscid. Contours from 283.39 K to 299.5 K.

Figure 4. x-velocity (top) and y-velocity (bottom) after 900 s for the nonlinear density current test. dx = dz = 100 m.

Table 3. Entries for Benjamin’s law

Time

d

600 s

800 m

α=

d Hdomain

0.125

ρ

"∆ρ" "

k = f (α)

g " = f (g, ∆ρ" , ρ)

1.16

0.02

1.208

0.16

numerical speed measured in our simulation and that obtained by Benjamin’s law is O(1.5ms−1 ), while as much as O(5ms−1 ) is the difference with that computed in Ref. 32. However, the two domains are slightly different in height and so are the depths of the currents, and hence a difference was expected.

7 of 10 American Institute of Aeronautics and Astronautics

Figure 5. Density current II: Initial state (top left), and profiles with vector field of the density current of Ref.32 Evolution after 300 s (top right), 600 s (bottom left), and 900 s (bottom right). Vector fields are shown for direct comparison with the reference case. We consider the front to be identified by the last grid cell on the bottom boundary where θ < Tref ; for this simulation the front reaches 9350 m at 600 s, with a speed of 15.66 ms−1 ; this value is not captured in this image because, being the ∆t computed at each time step, the time of the image is not the nominal 600 s.

IV.

Parallel Performance

The parallel performance of the algorithm was tested on a massively parallel architecturea up to 5000 processors. Automatic domain partitioning and MPI were used together with OpenMP for opening threads in elementary loops. The strong scalability is plotted in figure 6 for a three-dimensional, compressible boundary layer benchmark. Average # elements per CPU 13.2k 6.6k 3.3k

1600

1.7k

5000

Explicit NS, MPI Explicit NS, MPI+OpenMP Ideal

Speed up

Speed up

5k

Explicit NS, MAR Ideal

4000

1200

800

400

0

Average # elements per CPU 24k 12k

3000 2000 1000

0

400

800

1200

1600

0

0

1000

2000

3000

4000

5000

Figure 6. Strong scalability for a Euler compressible explicit 3D benchmark run on Marenostrum. Left, comparison between hybrid MPI+OpenMP and MPI schemes, up to 1600 cores. Right, MPI scheme up to 500 cores. On the x-axes is the number of cores involved (top scale) and the average FEM elements per core (bottom scale). The y-axis is the speed-up. In both cases, the ideal linear scalability is also plotted. For an extensive explanation on scalability strategies used in our code see.46

V.

Conclusions and Recommendations

We presented a new algorithm for the solution of the compressible Euler equations in the low Mach, stratified environment typical of atmospheric dynamics. The algorithm is an extension of the VMS technique developed by Hughes to compressible flows (VMS-C), and has been adapted to the treatment of thermally perturbed, stratified atmospheres in the context of nonhydrostatic atmospheric dynamics. The approach at hand, originally developed for CFD applications, showed to give encouraging results in terms of stability of the solution, accuracy, and parallel performance, despite being in its preliminary form. Further work is currently being carried for the solution on 9-node quadratic elements. Future work is planned for an extension to three dimensions and for the implementation of a fully implicit compressible solver; computational speed and large time steps are an issue when the application is that of numerical weather prediction, and a fully compressible implicit solver may enhance the system in this respect. a MareNostrum

at the Barcelona Supercomputing Center BSC-CNS.


Acknowledgments The authors are indebted to Dr. Francis X. Giraldo from the Naval Postgraduate School for fruitful and clarifying discussions, and Dr. John Zevenbergen from the Delft University of Technology, for his comments and corrections that helped improve this manuscript.

References 1 Giraldo, F. and Restelli, M., “A conservative semi-implicit discontinuous galerkin method for the Navier-Stokes equations in nonhydrostatic mesoscale modeling,” SIAM J. Sci. Comp., Vol. 31, 2009, pp. 2231–2257. 2 Drake, J., Foster, I., Michalakes, J., Toonen, B., and Worley, P., “Design and performance of a scalable parallel community climate model,” Parallel Comput., Vol. 21, 1995, pp. 1571–1592. 3 Hammond, S., Loft, R., Dennis, J., and Sato, R., “Implementation and performance issues of a massively parallel atmospheric model,” Parallel Comput., Vol. 21, 1995, pp. 1593–1610. 4 Barros, S., Dent, D., Isaksen, L., and Robinson, G., “The IFS model overview and parallel stategies,” Parallel supercomputing in atmospheric science: sixth ECMWF workshop on the use of Parallel processors in meteorology., No. 1, 1995, pp. 303–318. 5 Henderson, T., Baillie, C., Benjamin, S., Govett, M., Hart, L., Maroquin, A., Rodriguez, B., Black, T., Bleck, R., Carr, G., and Middlecoff, J., “Progress toward demonstrating operational capability of massively parallel processors at the Forecast Systems Laboratory,” Parallel supercomputing in atmospheric science: sixth ECMWF workshop on the use of Parallel processors in meteorology., No. 1, 1995, pp. 162–176. 6 Thomas, S., Loft, R., and Dennis, J., “Parallel implementation issues: global versus local methods,” Computing in Science and Engineering, Vol. 4, 2002, pp. 26–31. 7 Washington, W. and Parkinson, C., An introduction to three-dimensional climate modeling, University Science Books, 2nd ed., 2005. 8 Holmstrom, I., “On a method for parametric representation of the state of the atmosphere,” Tellus, Vol. 15, 1963, pp. 127–149. 9 Simons, T., “A three-dimensional spectral prediction equation,” J. Atmos. Sci., Vol. 127, 1968, pp. 1–27. 10 Cullen, M., “A finite element method for a non-linear initial value problem,” IMA J. Appl. Math., Vol. 13, 1974, pp. 233– 247. 11 Eliasen, E., Machenhauer, B., and Rasmussen, E., “On a numerical method for integration of the hydro-dynamical equations with a spectral representation of the horizontal fields,” Tech. Rep. 2, Institute of Meteorology, Univ. of Copenhagen, 1970. 12 Francis, P., “The possible use of Laguerre polynomials for representing the vertical structure of numerical models of the atmosphere,” Quarterly Journal of the Royal Meteorological Society, Vol. 98, 1972, pp. 662–667. 13 Burridge, D., Steppeler, J., and Struffing, R., “Finite element schemes for the vertical discretization of the ECMWF forecast model using linear elements,” Tech. Rep. 54, ECMWF, Sheffild Park, Reading - UK, 1986. 14 Untch, A. and Hortal, M., “A finite-element scheme for the vertical discretization of the semi-Lagrangian version of the ECMWF forecast model,” Q.J.R. Meteorol. Soc., Vol. 130, 2004, pp. 1505–1530. 15 Davies, T., Cullen, M., Malcolm, A., Mawson, M., Staniforth, A., White, A., and Wood, N., “A new dynamical core for the Met Office’s global and regional modelling of the atmosphere,” Quart. J. Roy. Meteor. Soc., Vol. 131, 2006, pp. 1759–1782. 16 Cote, J., Gravel, S., Methot, A., Patoine, A., Roch, M., and Staniforth, A., “The operational CMC-MRB global environmental multiscale (GEM) model. Part I: design considerations and formulation,” Mon. Wea. Rev., Vol. 126, 1998, pp. 1373–1395. 17 Cote, J., Desmarais, J., Gravel, S., Methot, A., Patoine, A., Roch, M., and Staniforth, A., “The Operational CMCMRB Global Environmental Multiscale (GEM) Model. Part II: results,” Mon. Wea. Rev., Vol. 126, 1998, pp. 1397–1418. 18 Yeh, K., Cote, J., Gravel, S., Methot, A., Patoine, A., Roch, M., and Staniforth, A., “The CMC-MRB Global Environmental Multiscale (GEM) Model. Part III: nonhydrostatic formulation,” Mon. Wea. Rev., Vol. 130, 2002, pp. 339–356. 19 Staniforth, A., “The application of the finite-element method to meteorological simulations - a review,” Int. J. Num. Meth. Fluids, Vol. 4, 1984, pp. 1–12. 20 Beland, M., Cote, J., and Staniforth, A., “The accuracy of a finite-element vertical discretization scheme for primitive equation models: Comparison with a finite-difference scheme,” Mon. Wea. Rev., Vol. 111, 1983, pp. 2298–2318. 21 Lauter, M., Giraldo, F., Handorf, D., and Dethloff, K., “A discontinuous Galerkin method for the shallow water equations using spherical triangular coordinates,” J. Comp. Phys., Vol. 227, 2008, pp. 10226–10242. 22 Giraldo, F. and Restelli, M., “A study of spectral element and discontinuous Galerkin methods for the Navier-Stokes equations in nonhydrostatic mesoscale atmospheric modeling: Equation sets and test cases,” J Comp Phys, Vol. 227, 2008, pp. 3849–3877. 23 Giraldo, F. and Rasmond, T., “A scalable Spectral Element Eulerian Atmospheric Model (SEE-AM) for numerical weather prediction: dynamical core tests,” Mon. Wea. Rev., Vol. 132, 2004, pp. 133–153. 24 Bonaventura, L., “A semi-implicit, semi-Lagrangian scheme using the height coordinate for a nonhydrostatic and fully elastic model of atmospheric flows,” J. Comp. Phys., Vol. 158, 2000, pp. 186–213. 25 Dawson, C., Westerink, J., Feyen, J., and Pothina, D., “Continuous, discontinuous and coupled discontinuous-continuous Galerkin Finite Element Methods for the shallow water equations,” Int. J. Numer. Meth. Fluids, Vol. 52, 2006, pp. 63–88. 26 Bacon, D., Ahmad, N., Boybeyi, Z., Dunn, T., Hall, M., Lee, C., Sarma, R., and Turner, M., “A dynamically adaptive weather and dispersion model: the operational multiscale environment model with grid adaptivity (OMEGA),” Monthly Weather Review , Vol. 128, 2000, pp. 2044–2075.


27 Tomita, H. and Satoh, M., “A new dynamical framework of nonhydrostatic global model using the icosahedral grid,” Fluid Dynamics Research, Vol. 34, 2004, pp. 357–400. 28 Satoh, M., Matsuno, T., Tomita, H., Miura, H., Nasuno, T., and Iga, S., “Nonhydrostatic icosahedral atmospheric model (NICAM) for global cloud resolving simulations,” J. Comp. Phys., 2008, pp. 3486–3514. 29 Skamarock, W., Doyle, J., Clark, P., and Wood, N., “A standard test set for nonhydrostatic dynamical cores of NWP models,” 2004, AMS NWP-WAF Conference Poster P2.17. 30 Wicker, L. and Skamarock, W., “A time-splitting scheme for the elastic equations incorporating second-order RungeKutta time differencing,” Mon. Weather Rev., Vol. 126, 1998, pp. 1992–1999. 31 Straka, J., Wilhelmson, R., Wicker, L., Anderson, J., and Droegemeier, K., “Numerical solution of a nonlinear density current: a benchmark solution and comparisons,” Int. J. Num. Meth. in Fluids, Vol. 17, 1993, pp. 1–22. 32 Haertfel, P., Johnson, R., and Tulich, S., “Some simple simulations of thunderstorm outflows,” J. Atmosph. Sci., Vol. 58, 2001, pp. 504–516. 33 Malcolm, A., “Evaluation of the proposed New Unified Model Scheme vs the Current Unified Model Scheme on the shallow water equations,” Tech. Rep. 180, MetOffice, Reading, 1996. 34 Schattler, U. and Doms, G., “The Nonhydrostatic Limited-Area Model LM (Lokal-Modell) of DWD - Part I: Implementation Documentation,” Tech. rep., Deutscher Wetterdienst, Meteorologische Analyse und Modellierung, Germany, 1999. 35 Skamarock, W., Klemp, J., Dudhia, J., Gill, D., Barker, D., Wang, W., and Powers, J., “A description of the Advanced Research WRF Version 2,” Tech. Rep. NCAR/TN 468+STR, NCAR, 2007. 36 Ahmad, N. and Lindeman, J., “Euler solutions using flux-based wave decomposition,” Int. J. Numer. Meth. Fluids, Vol. 54, 2007, pp. 47–72. 37 Vazquez, M., Houzeaux, G., and Moraques, M., “Variational multiscale stabilization of computational compressible flows,” Unpublished. 38 Hughes, T., “Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods,” Comp. Methods Appl. Mech. and Eng., Vol. 127, 1995, pp. 387–401. 39 Robert, A., “Bubble convection experiments with a semi-implicit formulation of the Euler equations,” J. Atmosph. Sci., Vol. 50, 1993, pp. 1865–1873. 40 Carpenter, R., Droegemeier, K., Woodward, P., and Hance, C., “Application of the piecewise parabolic method (PPM) to meteorological modeling,” Mon. Wea. Rev., Vol. 118, 1990, pp. 586–612. 41 Droegemeier, K. and Wilhelmson, R., “Numerical simulation of thunderstorm outflow dynamics. Part 1: outflow sensitivity experiments and turbulence dynamics,” J. Atmosph. Sci., Vol. 44, 1987, pp. 1180–1210. 42 Benjamin, T. B., “Gravity currents and related phenomena,” J. Fluid Mech., Vol. 31, 1968, pp. 209–248. 43 Mueller, C. K. and Carbone, R. E., “Dynamics of a thunderstorm outflow,” J. Atmos. Sci., Vol. 44, 1987, pp. 1879–1898. 44 CharbaJ., “Application of a gravity current model to analysis of squall line gust front,” Mon. Wea. Rev., Vol. 102, 1974, pp. 140–156. 45 Klemp, J. B., Rotunno, R., and Skamarock, W. C., “On the dynamics of gravity currents in a channel,” J. Fluid Mech., Vol. 269, 1994, pp. 169–198. 46 Houzeaux, G., V´ azquez, M., and Aubry, R., “A parallel incompressible Navier-Stokes solver for large scale supercomputers,” J. Comput. Phys., 2009, In Press.
