SHTns  3.4.6
Spatial data layouts and grids used by SHTns

The angular coordinates on a spherical shell are the co-latitude $ \theta $ and the longitude $ \phi $.

The spatial discretization is defined by a call to shtns_init; or to shtns_set_grid or shtns_set_grid_auto after shtns_create.

A field A (or one of its component in case of a vector field) is discretized on an ordered grid, which consists of

  • NPHI equally spaced nodes in longitude, spanning the range of angle between 0 (included) and $2\pi$/MRES (excluded),
  • NLAT gauss nodes or equally spaced nodes (depending on shtns_type) in latitude, spanning angles from 0 to $ \pi $.

There are constraints on the sizes NLAT and NPHI. The sampling theorem requires NPHI > 2*MMAX, with MMAX the highest Fourier mode to be resolved. Similarly NLAT > LMAX is required for a Gauss-Legendre quadrature, while on a regular grid NLAT > 2*LMAX is needed for the Féjer and Clenshaw-Curtis quadratures.

When dealing with non-linear equations, or more generally when analysing data that is not band-limited to LMAX and/or MMAX, anti-aliasing constraints exist. For a given non-linear order of the equation N (typically 2 for a quadratic non-linear term), the anti-aliasing conditions are :

  • NPHI > (N+1)*MMAX
  • NLAT > (N+1)*LMAX/2 for Gauss-Legendre quadrature (Gauss nodes, non-equispaced).
  • NLAT > (N+1)*LMAX for Féjer or Clenshaw-Curtis quadratures (equispaced nodes).

Furthermore, the FFT is more efficient for certain values of NPHI and/or NLAT. If you do not require specific NLAT or NPHI, you can let SHTns choose them for you, by using shtns_set_grid_auto with nphi and/or nlat set to zero.

Identifying a grid point with its indices it and ip in latitudinal and longitudinal direction respectively (indices start at 0), one has :

  • $ \theta $ = acos(ct[it]) where ct = shtns_info::ct is initialized by the call to shtns_init
  • $ \phi $ = ip * 2 $\pi$/(NPHI*MRES) and PHI_RAD and PHI_DEG can be used to get the phi angle corresponding to ip.
  • How to access $ A(\theta,\phi) $ depends on the data layout, but it is always a contiguous array of double precision floating point values (64 bit)

Currently, three data layouts are supported. Contiguous longitudes, Contiguous latitudes, and Native layout.

Contiguous longitudes

In this layout, increasing longitudes are stored next to each other for each latitude. That is $ A(\theta,\phi) $ = A[it*NPHI + ip] in C or A(ip,it) in Fortran. Use SHT_PHI_CONTIGUOUS to instruct shtns_init to use this spatial data layout :

shtns_cfg shtns_init(enum shtns_type flags, int lmax, int mmax, int mres, int nlat, int nphi)
Definition: sht_init.c:1514
use Contiguous longitudes
Definition: shtns.h:63
@ sht_gauss
use gaussian grid and quadrature. Allows on-the-fly or matrix-based algorithms.
Definition: shtns.h:53

will tell shtns to precompute everything for a gauss grid, and spatial data stored with longitude varying fastest.

Contiguous latitudes

In this layout, increasing latitude are stored next to each other for each longitude. That is $ A(\theta,\phi) $ = A[ip*NLAT + it] in C or A(it,ip) in Fortran. Use SHT_THETA_CONTIGUOUS to instruct shtns_init to use this spatial data layout.

Additionally, SHT_ALLOW_PADDING instructs shtns to optimize the layout to avoid cache bank conflicts. This can lead to significant performance boost (from 1% to 50% depending on the architecture). In that case, shtns_info::nlat_padded > shtns_info::nlat and shtns_info::nspat > shtns_info::nlat * shtns_info::nphi to reflect the data layout.

Native layout

The native way of storing spatial field data (which will help you achieve the best performance with SHTns) is the same as the Contiguous latitudes layout, except that it requires you to allocate slightly more memory for a field than the NLAT*NPHI double values. Namely NLAT*(NPHI/2+1)*2 doubles are required (instead of NLAT*NPHI) to be able tu use the in-place FFT of FFTW.

In Fortran this means you will allocate data as if the phi direction had (NPHI/2+1)*2 points instead of NPHI. The additional space, is located at the end of the useful data, and you don't need to worry about it.

To instruct shtns_init that your spatial data has been set up using this layout, use SHT_NATIVE_LAYOUT.

One must be careful when allocating spatial data, as the Fourier transform requires some additional space.
Precisely, the number of required double is NSPAT_ALLOC = NLAT*(NPHI/2+1)*2. The first NLAT*NPHI values are the spatial data.
In addition, spatial data must be allocated using the fftw_malloc function, for example using :
A = fftw_malloc( NSPAT_ALLOC * sizeof(double) );
#define NSPAT_ALLOC(shtns)
total number of 'doubles' required for a spatial field (includes FFTW reserved space).
Definition: shtns.h:114
It is not recommended to use the native layout with Fortran or Python.