Regular ArticleA Spectral Method for Polar Coordinates
Abstract
A new set of polynomial functions that can be used in spectral expansions of C∞ functions in polar coordinates (r, φ) is defined by a singular Sturm-Liouville equation. With the use of the basis functions, the spectral representations remain analytic at the pole despite the coordinate singularity because the pole condition is exactly satisfied at the origin for ag azimuthal modes, not just a few of the gravest modes (which is the usual case). Based on recurrence relations, fast and stable numerical operators for 1/r, r(d/dr), the Laplacian and Helmholtz operators and their inverses are developed. Although the spacings in the azimuthal direction of the collocation points near the origin are small (i.e., α 1/M2, where M is the number of radial modes), the explicit numerical method for Euler's equation is not stiff at the origin. Namely, the CFL number σ is O(1) where the grid size in is defined as π/M (i.e., the maximum allowable timestep is proportional to 1/M, not 1/M2).
References (0)
Cited by (93)
Gyroscopic polynomials
2023, Journal of Computational PhysicsGyroscopic alignment gives rise to highly spatially anisotropic columnar structures that in combination with complex domain boundaries pose challenges for efficient numerical discretizations and computations. We define gyroscopic polynomials to be three-dimensional polynomials expressed in a coordinate system that conforms to rotational alignment. We remap the original domain with radius-dependent boundaries onto a right cylindrical or annular domain to create the computational domain in this coordinate system. We find the volume element expressed in gyroscopic coordinates leads naturally to a hierarchy of orthonormal bases. We build the bases out of Jacobi polynomials in the vertical and generalized Jacobi polynomials in the radial. Because these coordinates explicitly conform to flow structures found in rapidly rotating systems the bases represent fields with a relatively small number of modes. We develop the operator structure for one-dimensional semi-classical orthogonal polynomials as a building block for differential operators in the full three-dimensional cylindrical and annular domains. The differentiation operators of generalized Jacobi polynomials generate a sparse linear system for discretization of differential operators acting on the gyroscopic bases. This enables efficient simulation of systems with strong gyroscopic alignment.
A gyroscopic polynomial basis in the sphere
2022, Journal of Computational PhysicsStandard spectral codes for full sphere dynamics utilize a combination of spherical harmonics and a suitable radial basis to represent fluid variables. These basis functions have a rotational invariance not present in geophysical flows. Gyroscopic alignment - alignment of dynamics along the axis of rotation - is a hallmark of geophysical fluids in the rapidly rotating regime. The Taylor-Proudman theorem, resulting from a dominant balance of the Coriolis force and the pressure gradient force, yields nearly invariant flows along this axial direction. In this paper we tailor a coordinate system to the cylindrical structures found in rotating spherical flows. This “spherindrical” coordinate system yields a natural hierarchy of basis functions, composed of Jacobi polynomials in the radial and vertical direction, regular throughout the ball. We expand fluid variables using this basis and utilize sparse Jacobi polynomial algebra to implement all operators relevant for partial differential equations in the spherical setting. We demonstrate the representation power of the basis in three eigenvalue problems for rotating fluids.
Tensor calculus in spherical coordinates using Jacobi polynomials. Part-I: Mathematical analysis and derivations
2019, Journal of Computational Physics: XCitation Excerpt :We also show how our basis sets map to other past works. Some authors use rescaled Jacobi polynomials in one form or another, but do not allow the parameters to change as derivatives act [3,29,32,51,64]. In §6 we provide conclusions for this Part-I of our two-part series.
This paper presents a method for accurate and efficient computations on scalar, vector and tensor fields in three-dimensional spherical polar coordinates. The method uses spin-weighted spherical harmonics in the angular directions and rescaled Jacobi polynomials in the radial direction. For the 2-sphere, spin-weighted harmonics allow for automating calculations in a fashion as similar to Fourier series as possible. Derivative operators act as wavenumber multiplication on a set of spectral coefficients. After transforming the angular directions, a set of orthogonal tensor rotations put the radially dependent spectral coefficients into individual spaces each obeying a particular regularity condition at the origin. These regularity spaces have remarkably simple properties under standard vector-calculus operations, such as gradient and divergence. We use a hierarchy of rescaled Jacobi polynomials for a basis on these regularity spaces. It is possible to select the Jacobi-polynomial parameters such that all relevant operators act in a minimally banded way. Altogether, the geometric structure allows for the accurate and efficient solution of general partial differential equations in the unit ball.
Tensor calculus in spherical coordinates using Jacobi polynomials. Part-II: Implementation and examples
2019, Journal of Computational Physics: XWe present a simulation code which can solve a broad range of partial differential equations in a full sphere. The code expands tensorial variables in a spectral series of spin-weighted spherical harmonics in the angular directions and a scaled Jacobi polynomial basis in the radial direction, as described in Vasil et al. (2018; hereafter, Part-I). Nonlinear terms are calculated by transforming from the coefficients of the spectral series to the value of each quantity on the physical grid, where it is easy to calculate products and perform other local operations. The expansion makes it straightforward to solve equations in tensor form (i.e., without decomposition into scalars). We propose and study several unit tests which demonstrate the code can accurately solve linear problems, implement boundary conditions, and transform between spectral and physical space. We then run a series of benchmark problems proposed in Marti et al. (2014), implementing the hydrodynamic and magnetohydrodynamic equations. We are able to calculate more accurate solutions than reported in Marti et al. (2014) by running at higher spatial resolution and using a higher-order timestepping scheme. We find the rotating convection and convective dynamo benchmark problems depend sensitively on details of timestepping and data analysis. We also demonstrate that in low resolution simulations of the dynamo problem, small changes in a numerical scheme can lead to large changes in the solution. To aid future comparison to these benchmarks, we include the source code used to generate the data, as well as the data and analysis scripts used to generate the figures.
Tensor calculus in polar coordinates using Jacobi polynomials
2016, Journal of Computational PhysicsSpectral methods are an efficient way to solve partial differential equations on domains possessing certain symmetries. The utility of a method depends strongly on the choice of spectral basis. In this paper we describe a set of bases built out of Jacobi polynomials, and associated operators for solving scalar, vector, and tensor partial differential equations in polar coordinates on a unit disk. By construction, the bases satisfy regularity conditions at for any tensorial field. The coordinate singularity in a disk is a prototypical case for many coordinate singularities. The work presented here extends to other geometries. The operators represent covariant derivatives, multiplication by azimuthally symmetric functions, and the tensorial relationship between fields. These arise naturally from relations between classical orthogonal polynomials, and form a Heisenberg algebra. Other past work uses more specific polynomial bases for solving equations in polar coordinates. The main innovation in this paper is to use a larger set of possible bases to achieve maximum bandedness of linear operations. We provide a series of applications of the methods, illustrating their ease-of-use and accuracy.
A fully spectral methodology for magnetohydrodynamic calculations in a whole sphere
2016, Journal of Computational PhysicsWe present a fully spectral methodology for magnetohydrodynamic (MHD) calculations in a whole sphere. The use of Jones–Worland polynomials for the radial expansion guarantees that the physical variables remain infinitely differentiable throughout the spherical volume. Furthermore, we present a mathematically motivated and systematic strategy to relax the very stringent time step constraint that is present close to the origin when a spherical harmonic expansion is used for the angular direction. The new constraint allows for significant savings even on relatively simple solutions as demonstrated on the so-called full sphere benchmark, a specific problem with a very accurately-known solution. The numerical implementation uses a 2D data decomposition which allows it to scale to thousands of cores on present-day high performance computing systems. In addition to validation results, we also present three new whole sphere dynamo solutions that present a relatively simple structure.