Cubed-sphere Grid Research Articles

Parallel and fully implicit simulations of the thermal convection in the Earth’s outer core

A parallel and fully implicit method is developed for numerical simulations of the thermal convection of an incompressible fluid in a rotating spherical shell. The method is based on a pseudo-compressibility approach with dual time stepping in order to tackle with the numerical difficulty of the solenoidal condition in the incompressible Navier-Stokes equations. The numerical solution of the nonlinear governing equations is based on the method of lines, whereby the partial differential equations are transformed into ordinary differential equations by a suitable spatial discretization. A second-order cell-centered finite volume discretization based on a collocated cubed-sphere grid is used here. In order to relax the time step limit and accurately integrate the semi-discrete nonlinear ordinary differential equations in time, we employ a second-order explicit-first-step, single-diagonal-coefficient, diagonally implicit Runge–Kutta (ESDIRK) method with adaptive time stepping. The nonlinear algebraic system arising at each pseudo-time step is solved by a Newton–Krylov–Schwarz algorithm to achieve good load balance on modern supercomputers. The numerical results are in good agreement with the benchmark solutions and more accurate than those obtained with a fractional step approach in our previous research. Large-scale tests on the Tianhe-2 supercomputer indicate that the fully implicit solver can achieve good parallel scalabilities in both strong and weak senses.

Verification of the Time-Split Method for Higher-Order Diffusion in the Spectral Element Method Model on a Cubed-Sphere Grid

This study used the time-split method for higher-order diffusion in a global numerical weather prediction model of the spectral element method on a cubed-sphere grid, and verified its performance compared to those of the explicit and implicit timestep diffusion schemes. Because the spectral element method on a cubed-sphere grid involves intrinsic characteristics of eigenvalues of the Laplacian operator, allowable timestep sizes are more restricted for stable time integration in higher-order diffusion. For example, the 6th-order diffusion significantly limits the timestep size for adequate filtering of small-scale numerical noise during simulation. The time-split scheme can be an alternative to overcome the severe constraints on the timestep size and allow the model a timestep size as large as is potentially acceptable even in higher-order diffusion, while maintaining atmospheric kinetic energy compared to that of the explicit and implicit schemes. During analysis of the amplification factor and numerical simulation of the schemes, it was shown that the time-split scheme has a theoretical response comparable to that of the explicit and implicit schemes. The numerical simulations were conducted using the Korean Integrated Model (KIM) with a 6th-order diffusion operator of three different diffusion schemes at resolutions of 0.25- or 0.125-degree grid spacing. The time-split scheme yields a forecast performance similar to both the explicit and implicit timestep schemes while greatly reducing the computational time.

Symmetric equations on the surface of a sphere as used by model GISS:IB

Abstract. Standard vector calculus formulas of Cartesian three space are projected onto the surface of a sphere. This produces symmetric equations with three nonindependent horizontal velocity components. Each orthogonal axis has a velocity component that rotates around its axis (eastward velocity rotates around the north–south axis) and a specific angular momentum component that is the product of the velocity component multiplied by the cosine of axis' latitude. Angular momentum components align with the fixed axes and simplify several formulas, whereas the rotating velocity components are not orthogonal and vary with location. Three symmetric coordinates allow vector resolution and calculus operations continuously over the whole spherical surface, which is not possible with only two coordinates. The symmetric equations are applied to one-layer shallow water models on cubed-sphere and icosahedral grids, the latter being computationally simple and applicable to an ocean domain. Model results are presented for three different initial conditions and five different resolutions.

High-Order Finite-Volume Method with Block-Based AMR for Magnetohydrodynamics Flows

A high-order central essentially non-oscillatory (CENO) finite volume scheme combined with a block-based adaptive mesh refinement (AMR) algorithm is proposed for the solution of the ideal magnetohydrodynamics equations. The high-order CENO finite-volume scheme is implemented with fourth-order spatial accuracy within a flexible multi-block, body-fitted, hexahedral grid framework. An important feature of the high-order adaptive approach is that it allows for anisotropic refinement, which can lead to large computational savings when anisotropic flow features such as isolated propagating fronts and/or waves, shocks, shear surfaces, and current sheets are present in the flow. This approach is designed to handle complex multi-block grid configurations, including cubed-sphere grids, where some grid blocks may have degenerate edges or corners characterized by missing neighboring blocks. A procedure for building valid high-order reconstruction stencils, even at these degenerate block edges and corners, is proposed, taking into account anisotropic resolution changes in a systematic and general way. Furthermore, a non-uniform or heterogeneous block structure is used where the ghost cells of a block containing the solution content of neighboring blocks are stored directly at the resolution of the neighbors. A generalized Lagrange multiplier divergence correction technique is applied to achieve numerically divergence-free magnetic fields while preserving high-order accuracy on the anisotropic AMR grids. Parallel implementation and local grid adaptivity are achieved by using a hierarchical block-based domain partitioning strategy in which the connectivity and refinement history of grid blocks are tracked using a flexible binary tree data structure. Physics-based refinement criteria as well as the CENO smoothness indicator are both used for directing the mesh refinement. Numerical results, including solution-driven anisotropic refinement of cubed-sphere grids, are presented to demonstrate the accuracy and efficiency of the approach.

GEOS-Chem High Performance (GCHP v11-02c): a next-generation implementation of the GEOS-Chem chemical transport model for massively parallel applications

Abstract. Global modeling of atmospheric chemistry is a grand computational challenge because of the need to simulate large coupled systems of ∼100–1000 chemical species interacting with transport on all scales. Offline chemical transport models (CTMs), where the chemical continuity equations are solved using meteorological data as input, have usability advantages and are important vehicles for developing atmospheric chemistry knowledge that can then be transferred to Earth system models. However, they have generally not been designed to take advantage of massively parallel computing architectures. Here, we develop such a high-performance capability for GEOS-Chem (GCHP), a CTM driven by meteorological data from the NASA Goddard Earth Observation System (GEOS) and used by hundreds of research groups worldwide. GCHP is a grid-independent implementation of GEOS-Chem using the Earth System Modeling Framework (ESMF) that permits the same standard model to operate in a distributed-memory framework for massive parallelization. GCHP also allows GEOS-Chem to take advantage of the native GEOS cubed-sphere grid for greater accuracy and computational efficiency in simulating transport. GCHP enables GEOS-Chem simulations to be conducted with high computational scalability up to at least 500 cores, so that global simulations of stratosphere–troposphere oxidant–aerosol chemistry at C180 spatial resolution (∼0.5∘×0.625∘) or finer become routinely feasible.

A Scalable High-Performance I/O System for a Numerical Weather Forecast Model on the Cubed-Sphere Grid

The design and implementation of a high-performance Input/Output (I/O) library for the Korean Integrated Model (KIM, KIM-IO) is described in this paper. The KIM is a next-generation global operational model for the Korea Meteorological Administration (KMA). The horizontal discretization of KIM consists of the spectral-element method on the cubed-sphere grid. The KIM-IO is developed to be a consistent and efficient approach for input and output of essential data in this particular grid structure in a multiprocessing environment. The KIM-IO provides three main features, comprising the sequential I/O, parallel I/O, and I/O decomposition methods, and adopts user-friendly interfaces similar to the Network Common Data Form (NetCDF). The efficiency of the KIM-IO is verified using experiments to analyze the performance of its three features. The scalability is also verified by implementing the KIMIO in the KIM at a resolution of approximately 12 km using the 4th supercomputer of KMA. The experimental results show that both regular parallel I/O and sequential I/O undergo performance degradation with an increasing number of processes. However, the I/O decomposition method in the KIM-IO overcomes this degradation, leading to improvement in scalability. The results also indicate that with using the new I/O decomposition method, the KIM attains good parallel scalability up to Ο (100,000) cores.

The Korean Integrated Model (KIM) System for Global Weather Forecasting

The Korea Institute of Atmospheric Prediction Systems (KIAPS) began a national project to develop a new global atmospheric model system in 2011. The ultimate goal of this 9-year project is to replace the current operational model at the Korea Meteorological Administration (KMA), which was adopted from the United Kingdom’s Meteorological Office’s unified model (UM) in 2010. The 12-km Korean Integrated Model (KIM) system, consisting of a spectral-element non-hydrostatic dynamical core on a cubed sphere grid and a state-of-the-art physics parameterization package, has been launched in a real-time forecast framework, with initial conditions obtained via the advanced hybrid four-dimensional ensemble variational data assimilation (4DEnVar) over its native grid. A development strategy for KIM and the evolution of its performance in medium-range forecasts toward a world-class global forecast system are described. Outstanding issues in KIM 3.1 as of February 2018 are discussed, along with a future plan for operational deployment in 2020.

Structure of Eigenvalues in the Advection-Diffusion Equation by the Spectral Element Method on a Cubed-Sphere Grid

To fundamentally understand discretized equations and differential operators, this study investigates the eigenvalues of the spherical gradient matrix and spherical Laplacian matrix discretized by the spectral element method (SEM) on the cubed-sphere grid (CS) for solid-body rotation. The gradient matrix with prescribed wind for the solid-body rotation has small positive real numbers in the eigenvalues, which implies that the solutions of the advection equation for solid-body rotation can be unstable and distorted. However, the unstable model can be effectively controlled by the Laplacian matrix, and all eigenvalues of the Laplacian matrix are aligned along the negative part of the real axis. The maximum eigenvalues of the gradient matrix and Laplacian matrix exhibit the characteristic of a linear function of the resolution, which is important to determine the upper limit of the time-step size for stable time integration. This study suggests suitable e-folding time scale of the diffusion using the 6th order for applications of the advection-diffusion equation. Regarding the characteristics of the eigenvalues, the setting of the stable time step is discussed. A method to implement the diffusion operators in the third-order Runge-Kutta time integration scheme in the advection-diffusion model is also discussed.

Multi-resolution Hybrid Data Assimilation Core on a Cubed-sphere Grid (HybDA)

This study illustrates the characteristics of the data assimilation system at the Korea Institute of Atmospheric Prediction Systems (KIAPS), based on the cubed-sphere grid system. The most interesting feature is the use of spherical harmonic functions defined on cubed-sphere grid points, which makes it possible to control the allowable physical wavenumber for the analysis increments. The relevant computational costs and parallel scalability are represented. The multiple-resolution approach is a distinguishable aspect of this data assimilation system. The wavenumber, up to which the analysis is conducted, increases as the outer iteration progresses. This multiresolution strategy is based on an investigation into the change of spectral components of analysis increments. The multi-resolution outer-loop provides cost-effective analysis-improvement, by explicitly controlling the analysis increments entered into the observation operator. To utilize the high-resolution deterministic forecast as a background state, it is subtracted from the forecast ensemble, to produce ensemble forecast perturbation that is hybridized with static background error covariance. Based on the cycled analysis experiments, the higher-resolution deterministic forecast is shown to preserve the high-frequency feature of the analysis increment relative to the ensemble mean forecast.

Seasonal range test run with Global Eta Framework

Abstract. Global Eta Framework (GEF) is a global atmospheric model developed in general curvilinear coordinates and capable of running on arbitrary rectangular quasi-uniform spherical grids, using stepwise (Eta) representation of the terrain. In this study, the model is run on a cubed-sphere grid topology, in a version with uniform Jacobians (UJ), which provides equal-area grid cells, and a smooth transition of coordinate lines across the edges of the cubed-sphere. Within a project at the Brazilian Center for Weather Forecasts and Climate Studies (CPTEC), a nonhydrostatic version of this model is under development and will be applied for seasonal prediction studies. This note describes preliminary tests with the GEF on the UJ cubed-sphere in which model performance is evaluated in seasonal simulations at a horizontal resolution of approximately 25 km, running in the hydrostatic mode. Comparison of these simulations with the ERA-Interim reanalyses shows that the 850 hPa temperature is underestimated, while precipitation pattern is mostly underestimated in tropical continental regions and overestimated in tropical oceanic regions. Nevertheless, the model is still able to well capture the main seasonal climate characteristics. These results will be used as a control run in further tests with the nonhydrostatic version of the model.

Runge-Kutta Discontinuous Local Evolution Galerkin Methods for the Shallow Water Equations on the Cubed-Sphere Grid

AbstractThe paper develops high order accurate Runge-Kutta discontinuous local evolution Galerkin (RKDLEG) methods on the cubed-sphere grid for the shallow water equations (SWEs). Instead of using the dimensional splitting method or solving one-dimensional Riemann problem in the direction normal to the cell interface, the RKDLEG methods are built on genuinely multi-dimensional approximate local evolution operator of the locally linearized SWEs on a sphere by considering all bicharacteristic directions. Several numerical experiments are conducted to demonstrate the accuracy and performance of our RKDLEG methods, in comparison to the Runge-Kutta discontinuous Galerkin method with Godunov's flux etc.

Vertical discretization with finite elements for a global hydrostatic model on the cubed sphere

A formulation of Galerkin finite element with basis-spline functions on a hybrid sigma-pressure coordinate is presented to discretize the vertical terms of global Eulerian hydrostatic equations employed in a numerical weather prediction system, which is horizontally discretized with high-order spectral elements on a cubed sphere grid. This replaces the vertical discretization of conventional central finite difference that is first-order accurate in non-uniform grids and causes numerical instability in advection-dominant flows. Therefore, a model remains in the framework of Galerkin finite elements for both the horizontal and vertical spatial terms. The basis-spline functions, obtained from the de-Boor algorithm, are employed to derive both the vertical derivative and integral operators, since Eulerian advection terms are involved. These operators are used to discretize the vertical terms of the prognostic and diagnostic equations. To verify the vertical discretization schemes and compare their performance, various two- and three-dimensional idealized cases and a hindcast case with full physics are performed in terms of accuracy and stability. It was shown that the vertical finite element with the cubic basis-spline function is more accurate and stable than that of the vertical finite difference, as indicated by faster residual convergence, fewer statistical errors, and reduction in computational mode. This leads to the general conclusion that the overall performance of a global hydrostatic model might be significantly improved with the vertical finite element.

Development of a hybrid parallel MCV-based high-order global shallow-water model

Utilization of high-order spatial discretizations is an important trend in developing global atmospheric models. As a competitive choice, the multi-moment constrained volume (MCV) method can achieve high accuracy while maintaining similar parallel scalability to classical finite volume methods. In this work, we introduce the development of a hybrid parallel MCV-based global shallow-water model on the cubed-sphere grid. Based on a sequential code, we perform parallelization on both the process and the thread levels. To enable process-level parallelism, we first decompose the six patches of the cubed-sphere in a same 2-D partition and then employ a conflict-free pipe-flow communication scheme for overlapping the halo exchange with computations. To further exploit the heterogeneous computing capacity of an Intel Xeon Phi accelerated supercomputer, we propose a guided panel-based inner---outer partition to distribute workload among the CPUs and the coprocessors. In addition to the above, thread-level parallelism along with various optimizations is done on both the multi-core CPU and the many-core accelerator. Numerical experiments are carried out to validate the correctness of the optimized parallel code and examine its parallel performance. Test results show that both the CPU-only and the hybrid codes scale well to hundreds of processes in terms of both the strong and weak scaling. In particular, the hybrid code can achieve a speedup of $$2.56\times $$2.56× as compared to the CPU-only version. In the largest run on a $$9216\,\times \,9216\,\times \,6$$9216×9216×6 mesh (1.5 billion unknowns), the hybrid code sustains an aggregative performance of 26.5 Tflops with 486 processes (33,534 cores).

Implementation of a Semi-implicit Time Integration Scheme in Non-Hydrostatic Euler Equations

A semi-implicit time integration scheme is implemented in a non-hydrostatic Euler problem on the cubed-sphere grid. The semi-implicit time integration scheme is a 3-stage additive RungeKutta method, which is an Implicit-Explicit (IMEX) multi-stage single-step scheme. The system treats acoustic and gravity waves implicitly and advection explicitly. In the implicit part, we compute the linear system by defining a proper linear operator. Then, the numerical results are presented to compare the semi-implicit time integration scheme with the explicit RungeKutta time integration scheme in non-hydrostatic Euler equations. In terms of accuracy, we demonstrate that the proposed method performs better than other schemes.

Characteristics of a Spectral Inverse of the Laplacian Using Spherical Harmonic Functions on a Cubed-Sphere Grid for Background Error Covariance Modeling

Abstract In this study, a spectral inverse method using spherical harmonic functions (SHFs) represented on a cubed-sphere grid (SHF inverse) is proposed. The purpose of the spectral inverse method studied is to help with data assimilation. The grid studied is the one that results from a spectral finite element decomposition of the six faces of the cubed sphere on Gauss–Legendre–Lobatto (GLL) points with equiangular gnomonic projection. For a given discretization of the cube in this form, as the total wavenumber of the test functions increases, there comes a point at which the cube’s eigenstructure fails to be able to replicate the spherical harmonic functions. The authors call this point a limit wavenumber in using the SHF inverse. In common with the authors’ previous research, the allowable total wavenumber of the SHF inverse increases more effectively with an enhanced polynomial order. The use of the eigenvectors and eigenvalues of the Laplacian, discretized on the grid spacing used in this study, to the Poisson equation is compared with the benchmark set by using the spherical harmonics solution to the problem. In terms of accuracy, the SHF inverse is superior to a direct inverse of the Laplacian using eigendecomposition. The feasibility of SHF inverse in operational implementation is examined under a massive computational environment.

Analyzing the Adaptive Mesh Refinement (AMR) Characteristics of a High-Order 2D Cubed-Sphere Shallow-Water Model

Abstract Adaptive mesh refinement (AMR) is a technique that has been featured only sporadically in atmospheric science literature. This paper aims to demonstrate the utility of AMR for simulating atmospheric flows. Several test cases are implemented in a 2D shallow-water model on the sphere using the Chombo-AMR dynamical core. This high-order finite-volume model implements adaptive refinement in both space and time on a cubed-sphere grid using a mapped-multiblock mesh technique. The tests consist of the passive advection of a tracer around moving vortices, a steady-state geostrophic flow, an unsteady solid-body rotation, a gravity wave impinging on a mountain, and the interaction of binary vortices. Both static and dynamic refinements are analyzed to determine the strengths and weaknesses of AMR in both complex flows with small-scale features and large-scale smooth flows. The different test cases required different AMR criteria, such as vorticity or height-gradient based thresholds, in order to achieve the best accuracy for cost. The simulations show that the model can accurately resolve key local features without requiring global high-resolution grids. The adaptive grids are able to track features of interest reliably without inducing noise or visible distortions at the coarse–fine interfaces. Furthermore, the AMR grids keep any degradations of the large-scale smooth flows to a minimum.

A global non-hydrostatic dynamical core using the spectral element method on a cubed-sphere grid

A new global model with a non-hydrostatic (NH) dynamical core is developed. It employs the spectral element method (SEM) in the horizontal discretization and the finite difference method (FDM) in the vertical discretization. The solver includes a time-split third-order Runge-Kutta (RK3) time-integration technique. Pursuing the quasi-uniform and pole singularity-free spherical geometry, a cubed-sphere grid is employed. To assess the performance of the developed dynamical solver, the results from a number of idealized benchmark tests for hydrostatic and non-hydrostatic flows are presented and compared. The results indicate that the non-hydrostatic dynamical solver is able to produce solutions with good accuracy and consistency comparable to reference solutions. Further evaluation of the model with a full-physics package demonstrates its capability in reproducing heavy rainfall over the Korean Peninsula, which confirms that coupling of the dynamical solver and full-physics package is robust.

Spectral Transformation Using a Cubed-Sphere Grid for a Three-Dimensional Variational Data Assimilation System

Abstract Atmospheric numerical models using the spectral element method with cubed-sphere grids (CSGs) are highly scalable in terms of parallelization. However, there are no data assimilation systems for spectral element numerical models. The authors devised a spectral transformation method applicable to the model data on a CSG (STCS) for a three-dimensional variational data assimilation system (3DVAR). To evaluate the 3DVAR system based on the STCS, the authors conducted observing system simulation experiments (OSSEs) using Community Atmosphere Model with Spectral Element dynamical core (CAM-SE). They observed root-mean-squared error reductions: 24% and 34% for zonal and meridional winds (U and V), respectively; 20% for temperature (T); 4% for specific humidity (Q); and 57% for surface pressure (Ps) in analysis and 28% and 27% for U and V, respectively; 25% for T; 21% for Q; and 31% for Ps in 72-h forecast fields. In this paper, under the premise that the same number of grid points is set, the authors show that the use of a greater polynomial degree, np, produces better performance than use of a greater element count, ne, on equiangular coordinates in terms of the wave representation.

High-Order Finite-Volume Transport on the Cubed Sphere: Comparison between 1D and 2D Reconstruction Schemes

Abstract This paper presents two finite-volume (FV) schemes for solving linear transport problems on the cubed-sphere grid system. The schemes are based on the central-upwind finite-volume (CUFV) method, which is a class of Godunov-type method for solving hyperbolic conservation laws, and combines the attractive features of the classical upwind and central FV methods. One of the CUFV schemes is based on a dimension-by-dimension approach and employs a fifth-order one-dimensional (1D) Weighted Essentially Nonoscillatory (WENO5) reconstruction method. The other scheme employs a fully two-dimensional (2D) fourth-order accurate reconstruction method. The cubed-sphere grid system imposes several computational challenges due to its patched-domain topology and nonorthogonal curvilinear grid structure. A high-order 1D interpolation procedure combining cubic and quadratic interpolations is developed for the FV schemes to handle the discontinuous edges of the cubed-sphere grid. The WENO5 scheme is compared against the fourth-order Kurganov–Levy (KL) scheme formulated in the CUFV framework. The performance of the schemes is compared using several benchmark problems such as the solid-body rotation and deformational-flow tests, and empirical convergence rates are reported. In addition, a bound-preserving filter combined with an optional positivity-preserving filter is tested for nonsmooth problems. The filtering techniques considered are local, inexpensive, and effective. A fourth-order strong stability preserving explicit Runge–Kutta time-stepping scheme is used for integration. The results show that schemes are competitive to other published FV schemes in the same category.

High-order finite volume shallow water model on the cubed-sphere: 1D reconstruction scheme

A central-upwind finite-volume (CUFV) scheme for shallow-water model on a nonorthogonal equiangular cubed-sphere grid is developed, consequently extending the 1D reconstruction CUFV transport scheme developed by us. High-order spatial discretization based on weighted essentially non-oscillatory (WENO) is considered for this effort. The CUFV method combines the alluring features of classical upwind and central schemes. This approach is particularly useful for complex computational domain such as the cubed-sphere. The continuous flux-form spherical shallow water equations in nonorthogonal curvilinear coordinates are utilized. Fluxes along the element boundaries are approximated by a Kurganov–Noelle–Petrova scheme. A fourth-order strong stability preserving Runge–Kutta time stepping scheme for time integration is employed in the present work. The numerical scheme is evaluated with standard shallow water test suite, which accentuate accuracy and conservation properties. In addition, an efficient yet inexpensive bound preserving filter with an optional positivity filter is used to remove spurious oscillations and to achieve strictly positive definite numerical solution. To tackle the discontinuities that arise at the edges of the cubed-sphere grid, we utilize a high-order 1D interpolation procedure combining cubic and quadratic interpolations.The results with the high-order scheme is compared with the results for the same tests for various schemes available in literature. Since, the scheme presented here uses local-cell information, it is expected to be scalable to high number of processors count in a distributed node high-performance computer.