Spherical Geodesic Grid Research Articles

A finite element method for angular discretization of the radiation transport equation on spherical geodesic grids

Discrete ordinate (SN) and filtered spherical harmonics (FPN) based schemes have been proven to be robust and accurate in solving the Boltzmann transport equation but they have their own strengths and weaknesses in different physical scenarios. We present a new method based on a finite element approach in angle that combines the strengths of both methods and mitigates their disadvantages. The angular variables are specified on a spherical geodesic grid with functions on the sphere being represented using a finite element basis. A positivity-preserving limiting strategy is employed to prevent non-physical values from appearing in the solutions. The resulting method is then compared with both SN and FPN schemes using four test problems and is found to perform well when one of the other methods fail.

SELEN<sup>4</sup> (SELEN version 4.0): a Fortran program for solving the gravitationally and topographically self-consistent sea-level equation in glacial isostatic adjustment modeling

Abstract. We present SELEN4 (SealEveL EquatioN solver), an open-source program written in Fortran 90 that simulates the glacial isostatic adjustment (GIA) process in response to the melting of the Late Pleistocene ice sheets. Using a pseudo-spectral approach complemented by a spatial discretization on an icosahedron-based spherical geodesic grid, SELEN4 solves a generalized sea-level equation (SLE) for a spherically symmetric Earth with linear viscoelastic rheology, taking the migration of the shorelines and the rotational feedback on sea level into account. The approach is gravitationally and topographically self-consistent, since it considers the gravitational interactions between the solid Earth, the cryosphere, and the oceans, and it accounts for the evolution of the Earth's topography in response to changes in sea level. The SELEN4 program can be employed to study a broad range of geophysical effects of GIA, including past relative sea-level variations induced by the melting of the Late Pleistocene ice sheets, the time evolution of paleogeography and of the ocean function since the Last Glacial Maximum, the history of the Earth's rotational variations, present-day geodetic signals observed by Global Navigation Satellite Systems, and gravity field variations detected by satellite gravity missions like GRACE (the Gravity Recovery and Climate Experiment). The “GIA fingerprints” constitute a standard output of SELEN4. Along with the source code, we provide a supplementary document with a full account of the theory, some numerical results obtained from a standard run, and a user guide. Originally, the SELEN program was conceived by Giorgio Spada (GS) in 2005 as a tool for students eager to learn about GIA, and it has been the first SLE solver made available to the community.

An Adaptive Wavelet Collocation Method for Solution of the Convection-Dominated Problem on a Sphere

In this paper, we present a numerical simulation of a convection-dominated problem using adaptive wavelet collocation method. This is based on second-generation spherical wavelets on a recursively refined spherical geodesic grid. The power lies in the fact that they only require small number of coefficients to represent this problem more accurately, which allows compression and efficient computation. In contrast with other approximate schemes, this approach reduced the number of collocation points used by factors of many orders of magnitude when compared to uniform grids of equivalent resolution. Numerical results demonstrate the efficiency of our approach.

Approximation of the differential operators on an adaptive spherical geodesic grid using spherical wavelets

In this work, a new adaptive multi-level approximation of surface divergence and scalar-valued surface curl operator on a recursively refined spherical geodesic grid is presented. A hierarchical finite volume scheme based on the wavelet multi-level decomposition is used to approximate the surface divergence and scalar-valued surface curl operator. The multi-level structure provides a simple way to adapt the computation to the local structure of the surface divergence and scalar-value surface curl operator so that the high resolution computations are performed only in regions where singularities or sharp transitions occur. This multi-level approximation of the surface divergence operator is then used in an adaptive wavelet collocation method (AWCM) to solve two standard advection tests, solid-body rotation and divergent flow on the sphere. In contrast with other approximate schemes, this approach can be extended easily to other curved manifolds by considering appropriate coarse approximation to the desired manifold (here we used the icosahedral approximation to the sphere at the coarsest level) and using recursive surface subdivision.

A Dynamic Adaptive Wavelet Method for Solution of the Schrodinger Equation

In this paper, we present a dynamically adaptive wavelet method for solving Schrodinger equation on one-dimensional, two-dimensional and on the sphere. Solving one-dimensional and two-dimensional Schrodinger equations are based on Daubechies wavelet with finite difference method on an arbitrary grid, and for spherical Schrodinger equation is based on spherical wavelet over an optimal spherical geodesic grid. The method is applied to the solution of Schrodinger equation for computational efficiency and achieve accuracy with controlling spatial grid adaptation — high resolution computations are performed only in regions where a solution varies greatly (i.e., near steep gradients, or near-singularities) and a much coarser grid where the solution varies slowly. Thereupon the dynamic adaptive wavelet method is useful to analyze local structure of solution with very less number of computational cost than any other methods. The prowess and computational efficiency of the adaptive wavelet method is demonstrated for the solution of Schrodinger equation on one-dimensional, two-dimensional and on the sphere.

Multilevel approximation of the gradient operator on an adaptive spherical geodesic grid

This work presents a new adaptive multilevel approximation of the gradient operator on a recursively refined spherical geodesic grid. The multilevel structure provides a simple way to adapt the computation to the local structure of the gradient operator so that high resolution computations are performed only in regions where singularities or sharp transitions occur. This multilevel approximation of the gradient operator is used to solve the linear spherical advection equation for both time-independent and time-dependent wind field geophysical test cases. In contrast with other approximation schemes, this approach can be extended easily to other curved manifolds by choosing an appropriate coarse approximation and using recursive surface subdivision. The results indicate that the adaptive gradient calculation and the solution of spherical advection equation accurate, efficient and free of numerical dispersion.

MAGNETOHYDRODYNAMIC MODELING OF SOLAR SYSTEM PROCESSES ON GEODESIC GRIDS

This report describes a new magnetohydrodynamic numerical model based on a hexagonal spherical geodesic grid. The model is designed to simulate astrophysical flows of partially ionized plasmas around a central compact object, such as a star or a planet with a magnetic field. The geodesic grid, produced by a recursive subdivision of a base platonic solid (an icosahedron), is free from control volume singularities inherent in spherical polar grids. Multiple populations of plasma and neutral particles, coupled via charge-exchange interactions, can be simulated simultaneously with this model. Our numerical scheme uses piecewise linear reconstruction on a surface of a sphere in a local two-dimensional "Cartesian" frame. The code employs HLL-type approximate Riemann solvers and includes facilities to control the divergence of magnetic field and maintain pressure positivity. Several test solutions are discussed, including a problem of an interaction between the solar wind and the local interstellar medium, and a simulation of Earth's magnetosphere.

Integration of barotropic vorticity equation over spherical geodesic grid using multilevel adaptive wavelet collocation method

In this paper, we present the multilevel adaptive wavelet collocation method for solving non-divergent barotropic vorticity equation over spherical geodesic grid. This method is based on multi-dimensional second generation wavelet over a spherical geodesic grid. The method is more useful in capturing, identifying, and analyzing local structure [1] than any other traditional methods (i.e. finite difference, spectral method), because those methods are either full or partial miss important phenomena such as trends, breakdown points, discontinuities in higher derivatives of the solution. Wavelet decomposition is used for interpolation and adaptive grid refinement on different levels.

Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Reconstructions for Forward-in-Time Schemes

Abstract The finite-volume transport scheme of Miura, for icosahedral–hexagonal meshes on the sphere, is extended by using higher-order reconstructions of the transported scalar within the formulation. The use of second- and fourth-order reconstructions, in contrast to the first-order reconstruction used in the original scheme, results in significantly more accurate solutions at a given mesh density, and better phase and amplitude error characteristics in standard transport tests. The schemes using the higher-order reconstructions also exhibit much less dependence of the solution error on the time step compared to the original formulation. The original scheme of Miura was only tested using a nondeformational time-independent flow. The deformational time-dependent flow test used to examine 2D planar transport in Blossey and Durran is adapted to the sphere, and the schemes are subjected to this test. The results largely confirm those generated using the simpler tests. The results also indicate that the scheme using the second-order reconstruction is most efficient and its use is recommended over the scheme using the first-order reconstruction. The second-order reconstruction uses the same computational stencil as the first-order reconstruction and thus does not create any additional parallelization issues.

Fourier Analysis of Multigrid Methods on Hexagonal Grids

This paper applies local Fourier analysis to multigrid methods on hexagonal grids. Using oblique coordinates to express the grids and a dual basis for the Fourier modes, the analysis proceeds essentially the same as for rectangular grids. The framework for one- and two-grid analyses is given and then applied to analyze the performance of multigrid methods for the Poisson problem on a hexagonal grid. Numerical results confirm the analysis. Uniform hexagonal grids provide an approximation to spherical geodesic grids; numerical results for the latter show similar performance. While the analysis is similar to that for rectangular grids, the results differ somewhat: full weighting is superior to injection for restriction, Jacobi relaxation performs about as well as Gauss–Seidel relaxation, and underrelaxation is not required for good performance. Also, coarse-fine or four-color ordering (both analogues of red-black ordering on the rectangular grid) improves the performance of Jacobi relaxation, with the latter achieving a smoothing factor of approximately 0.25. An especially simple compact fourth-order discretization works well, and the full multigrid algorithm produces the solution to the level of truncation error in work proportional to the number of unknowns.

An Adaptive Multilevel Wavelet Solver for Elliptic Equations on an Optimal Spherical Geodesic Grid

An adaptive multilevel wavelet solver for elliptic equations on an optimal spherical geodesic grid is developed. The method is based on second-generation spherical wavelets on almost uniform optimal spherical geodesic grids. It is an extension of the adaptive multilevel wavelet solver [O. V. Vasilyev and N. K.-R. Kevlahan, J. Comput. Phys., 206 (2005), pp. 412–431] to curved manifolds. Wavelet decomposition is used for grid adaption and interpolation. A hierarchical finite difference scheme based on the wavelet multilevel decomposition is used to approximate the Laplace–Beltrami operator. The optimal spherical geodesic grid [Internat. J. Comput. Geom. Appl., 16 (2006), pp. 75–93] is convergent in terms of local mean curvature and has lower truncation error than conventional spherical geodesic grids. The overall computational complexity of the solver is $O(\mathcal{N})$, where $\mathcal{N}$ is the number of grid points after adaptivity. The accuracy and efficiency of the method is demonstrated for the spherical Poisson equation. Although the present paper considers the sphere, the strength of this new method is that it can be extended easily to other curved manifolds by choosing an appropriate coarse approximation and using recursive surface subdivision.

Semi‐Lagrangian advection on a spherical geodesic grid

AbstractA simple and efficient numerical method for solving the advection equation on the spherical surface is presented. To overcome the well‐known ‘pole problem’ related to the polar singularity of spherical coordinates, the space discretization is performed on a geodesic grid derived by a uniform triangulation of the sphere; the time discretization uses a semi‐Lagrangian approach. These two choices, efficiently combined in a substepping procedure, allow us to easily determine the departure points of the characteristic lines, avoiding any computationally expensive tree‐search. Moreover, suitable interpolation procedures on such geodesic grid are presented and compared. The performance of the method in terms of accuracy and efficiency is assessed on two standard test cases: solid‐body rotation and a deformation flow. Copyright © 2007 John Wiley & Sons, Ltd.

Interpolation on spherical geodesic grids: A comparative study

Most operational models in atmospheric physics, meteorology and climatology nowadays adopt spherical geodesic grids and require “ad hoc” developed interpolation procedures. The author does a comparison between chosen representatives of linear, distance-based and cubic interpolation schemes outlining their advantages and drawbacks in this specific application field. Numerical experiments on a standard test problem, while confirming a good performance of linear and distance-based schemes in a single interpolation step, also show their minor accuracy with respect to the cubic scheme in the more realistic simulation of advection of a meteorological field.

An Incremental Remapping Transport Scheme on a Spherical Geodesic Grid

Abstract Weather and climate models contain equations for transporting conserved quantities such as the mass of air, water, ice, and associated tracers. Ideally, the numerical schemes used to solve these equations should be conservative, spatially accurate, and monotonicity-preserving. One such scheme is incremental remapping, previously developed for transport on quadrilateral grids. Here the incremental remapping scheme is reformulated for a spherical geodesic grid whose cells are hexagons and pentagons. The scheme is tested in a shallow-water model with both uniform and varying velocity fields. Solutions for standard shallow-water test cases 1, 2, and 5 are obtained with a centered scheme, a flux-corrected transport (FCT) scheme, and the remapping scheme. The three schemes are about equally accurate for transport of the height field. For tracer transport, remapping is far superior to the centered scheme, which produces large overshoots, and is generally smoother and more accurate than FCT. Remapping has a high startup cost associated with geometry calculations but is nearly twice as fast as FCT for each added tracer. As a result, remapping is cheaper than FCT for transport of more than about seven tracers.

Climate modeling with spherical geodesic grids

A new approach to climate simulation uses geodesic grids generated from an icosahedron and could become on attractive alternative to current models. We implemented an atmospheric general circulation model using a geodesic discretization of the sphere. Our model uses the message-passing interface and runs efficiently on massively parallel machines.

A Potential Enstrophy and Energy Conserving Numerical Scheme for Solution of the Shallow-Water Equations on a Geodesic Grid

Using the shallow water equations, a numerical framework on a spherical geodesic grid that conserves domain-integrated mass, potential vorticity, potential enstrophy, and total energy is developed. The numerical scheme is equally applicable to hexagonal grids on a plane and to spherical geodesic grids. This new numerical scheme is compared to its predecessor and it is shown that the new scheme does considerably better in conserving potential enstrophy and energy. Furthermore, in a simulation of geostrophic turbulence, the new numerical scheme produces energy and enstrophy spectra with slopes of approximately K−3 and K−1, respectively, where K is the total wavenumber. These slopes are in agreement with theoretical predictions. This work also exhibits a discrete momentum equation that is compatible with the Z-grid vorticity-divergence equation.

A spectral element shallow water model on spherical geodesic grids

SUMMARY The spectral element method for the two-dimensional shallow water equations on the sphere is presented. The equations are written in conservation form and the domains are discretized using quadrilateral elements obtained from the generalized icosahedral grid introduced previously (Giraldo FX. Lagrange‐ Galerkin methods on spherical geodesic grids: the shallow water equations. Journal of Computational Physics 2000; 160: 336‐368). The equations are written in Cartesian co-ordinates that introduce an additional momentum equation, but the pole singularities disappear. This paper represents a departure from previously published work on solving the shallow water equations on the sphere in that the equations are all written, discretized, and solved in three-dimensional Cartesian space. Because the equations are written in a three-dimensional Cartesian co-ordinate system, the algorithm simplifies into the integration of surface elements on the sphere from the fully three-dimensional equations. A mapping (Song Ch, Wolf JP. The scaled boundary finite element method—alias consistent infinitesimal finite element cell method—for diffusion. International Journal for Numerical Methods in Engineering 1999; 45: 1403‐1431) which simplifies these computations is described and is shown to contain the Eulerian version of the method introduced previously by Giraldo (Journal of Computational Physics 2000; 160: 336‐368) for the special case of triangular elements. The significance of this mapping is that although the equations are written in Cartesian co-ordinates, the mapping takes into account the curvature of the high-order spectral elements, thereby allowing the elements to lie entirely on the surface of the sphere. In addition, using this mapping simplifies all of the three-dimensional spectral-type finite element surface integrals because any of the typical two-dimensional planar finite element or spectral element basis functions found in any textbook (for example, Huebner et al. The Finite Element Method for Engineers. Wiley, New York, 1995; Karniadakis GE, Sherwin SJ. Spectral:hp Element Methods for CFD. Oxford University ( ˘

Modeling the Atmospheric General Circulation Using a Spherical Geodesic Grid: A New Class of Dynamical Cores

This paper documents the development and testing of a new type of atmospheric dynamical core. The model solves the vorticity and divergence equations in place of the momentum equation. The model is discretized in the horizontal using a geodesic grid that is nearly uniform over the entire globe. The geodesic grid is formed by recursively bisecting the triangular faces of a regular icosahedron and projecting those new vertices onto the surface of the sphere. All of the analytic horizontal operators are reduced to line integrals, which are numerically evaluated with second-order accuracy. In the vertical direction the model can use a variety of coordinate systems, including a generalized sigma coordinate that is attached to the top of the boundary layer. Terms related to gravity wave propagation are isolated and an efficient semi-implicit time-stepping scheme is implemented. Since this model combines many of the positive attributes of both spectral models and conventional finite-difference models into a single dynamical core, it represents a distinctively new approach to modeling the atmosphere’s general circulation. The model is tested using the idealized forcing proposed by Held and Suarez. Results are presented for simulations using 2562 polygons (approximately 4.5° × 4.5°) and using 10 242 polygons (approximately 2.25° × 2.25°). The results are compared to those obtained with spectral model simulations truncated at T30 and T63. In terms of first and second moments of state variables such as the zonal wind, meridional wind, and temperature, the geodesic grid model results using 2562 polygons are comparable to those of a spectral model truncated at slightly less than T30, while a simulation with 10 242 polygons is comparable to a spectral model simulation truncated at slightly less than T63. In order to further demonstrate the viability of this modeling approach, preliminary results obtained from a full-physics general circulation model that uses this dynamical core are presented. The dominant features of the DJF climate are captured in the full-physics simulation. In terms of computational efficiency, the geodesic grid model is somewhat slower than the spectral model used for comparison. Model timings completed on an SGI Origin 2000 indicate that the geodesic grid model with 10 242 polygons is 20% slower than the spectral model truncated at T63. The geodesic grid model is more competitive at higher resolution than at lower resolution, so further optimization and future trends toward higher resolution should benefit the geodesic grid model.

Lagrange–Galerkin Methods on Spherical Geodesic Grids: The Shallow Water Equations

The weak Lagrange–Galerkin finite element method for the 2D shallow water equations on the sphere is presented. This method offers stable and accurate solutions because the equations are integrated along the characteristics. The equations are written in 3D Cartesian conservation form and the domains are discretized using linear triangular elements. The use of linear triangular elements permits the construction of accurate (by virtue of the second-order spatial and temporal accuracies of the scheme) and efficient (by virtue of the less stringent CFL condition of Lagrangian methods) schemes on unstructured domains. Using linear triangles in 3D Cartesian space allows for the explicit construction of area coordinate basis functions thereby simplifying the calculation of the finite element integrals. The triangular grids are constructed by a generalization of the icosahedral grids that have been typically used in recent papers. An efficient searching strategy for the departure points is also presented for these generalized icosahedral grids which involves very few floating point operations. In addition a high-order scheme for computing the characteristic curves in 3D Cartesian space is presented: a general family of Runge–Kutta schemes. Results for six test cases are reported in order to confirm the accuracy of the scheme.

Trajectory Calculations for Spherical Geodesic Grids in Cartesian Space

Abstract This paper shows how to obtain accurate and efficient trajectory calculations for spherical geodesic grids in Cartesian space. Determination of the departure points is essential to characteristic-based methods that trace the value of a function to the foot of the characteristics and then either integrate or interpolate at this location. In this paper, the departure points are all computed in relation to the spherical geodesic grids that are composed of a disjoint set of unstructured equilateral triangles. Interpolating and noninterpolating trajectory calculation approaches are both illustrated and the accuracy of both methods are compared. The noninterpolating method of McGregor results in the most accurate trajectories. The challenge in using McGregor’s method on unstructured triangular grids lies in the computation of the derivatives required in the high-order terms of the Taylor series expansion. This paper extends McGregor’s method to unstructured triangular grids by describing an accurate an...