GitHub

SpeedyTransforms

SpeedyTransforms is a submodule that has been developed for SpeedyWeather.jl which is technically independent (SpeedyWeather.jl however imports it) and can also be used without running simulations. It is just not put into its own respective repository for now.

The SpeedyTransforms are based on RingGrids and LowerTriangularMatrices to hold data in either grid-point space or in spectral space. So you want to read these sections first for clarifications how to work with these. We will also not discuss mathematical details of the Spherical Harmonic Transform here, but will focus on the usage of the SpeedyTransforms module.

The SpeedyTransforms module also implements the gradient operators $\nabla, \nabla \cdot, \nabla \times, \nabla^2, \nabla^{-2}$ in spectral space. Combined with the spectral transform, you could for example start with a velocity field in grid-point space, transform to spectral, compute its divergence and transform back to obtain the divergence in grid-point space. Examples are outlined in Gradient operators.

SpeedyTransforms assumes a unit sphere

The operators in SpeedyTransforms generally assume a sphere of radius $R=1$. For the transforms themselves that does not make a difference, but the gradient operators div, curl, ∇, ∇², ∇⁻² omit the radius scaling. You will have to do this manually. Also note that the meridional derivate expects a $\cos^{-1}(\theta)$ scaling. More in Gradient operators.

Example transform

Lets start with a simple transform. We could be using SpeedyWeather but to be more verbose these are the modules required to load

using SpeedyWeather.RingGrids
using SpeedyWeather.LowerTriangularMatrices
using SpeedyWeather.SpeedyTransforms

As an example, we want to transform the $l=m=1$ spherical harmonic from spectral space in alms to grid-point space.

alms = zeros(LowerTriangularMatrix{ComplexF64}, 6, 6)     # spectral coefficients
alms[2, 2] = 1                                           # only l=1, m=1 harmonic
alms
6×6 LowerTriangularMatrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

Now gridded is the function that takes spectral coefficients alms and converts them a grid-point space map

map = gridded(alms)
128-element, 8-ring FullGaussianGrid{Float64}:
 -0.19278869685896918
 -0.17811353112752462
 -0.13632219488509478
 -0.07377704023518313
  0.0
  0.07377704023518313
  0.13632219488509478
  0.17811353112752462
  0.19278869685896918
  0.17811353112752462
  ⋮
  0.17811353112752462
  0.19278869685896918
  0.17811353112752462
  0.13632219488509478
  0.07377704023518313
  0.0
 -0.07377704023518313
 -0.13632219488509478
 -0.17811353112752462

By default, the gridded transforms onto a FullGaussianGrid unravelled here into a vector west to east, starting at the prime meridian, then north to south, see RingGrids. We can visualize map quickly with a UnicodePlot via plot (see Visualising RingGrid data)

import SpeedyWeather.RingGrids: plot    # not necessary when `using SpeedyWeather`
plot(map)
       8-ring FullGaussianGrid{Float64}     
       ┌────────────────┐  0.7
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └────────────────┘ -0.7
        0     ˚E     360

Yay! This is the what the $l=m=1$ spherical harmonic is supposed to look like! Now let's go back to spectral space with spectral

alms2 = spectral(map)
6×6 LowerTriangularMatrix{ComplexF64}:
 0.0+0.0im           0.0+0.0im  0.0+0.0im  …  0.0+0.0im           0.0+0.0im
 0.0+0.0im           1.0+0.0im  0.0+0.0im     0.0+0.0im           0.0+0.0im
 0.0+0.0im           0.0+0.0im  0.0+0.0im     0.0+0.0im           0.0+0.0im
 0.0+0.0im   3.75734e-16+0.0im  0.0+0.0im     0.0+0.0im           0.0+0.0im
 0.0+0.0im           0.0+0.0im  0.0+0.0im     0.0+0.0im           0.0+0.0im
 0.0+0.0im  -1.63058e-17+0.0im  0.0+0.0im  …  0.0+0.0im  -2.47957e-17+0.0im

Comparing with alms from above you can see that the transform is exact up to a typical rounding error from Float64. 

alms ≈ alms2
true

YAY! The transform is typically idempotent, meaning that either space may hold information that is not exactly representable in the other but the first two-way transform will remove that so that subsequent transforms do not change this any further. However, also note here that the default FullGaussianGrid is an exact grid, inexact grids usually have a transform error that is larger than the rounding error from floating-point arithmetic.

Transform onto another grid

While the default grid for SpeedyTransforms is the FullGaussianGrid we can transform onto other grids by specifying Grid too

map = gridded(alms, Grid=HEALPixGrid)
plot(map)
       7-ring HEALPixGrid{Float64}     
       ┌──────────────┐  0.7
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └──────────────┘ -0.7
        0    ˚E    360

which, if transformed back, however, can yield a larger transform error as discussed above

spectral(map)
6×6 LowerTriangularMatrix{ComplexF64}:
 0.0+0.0im           0.0+0.0im          0.0+0.0im  …  0.0+0.0im  0.0+0.0im
 0.0+0.0im       1.01215-1.4295e-17im   0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im   3.42299e-17+0.0im          0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im     0.0185292-1.63433e-17im  0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im  -6.28844e-17+0.0im          0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im    -0.0164851+1.42335e-17im  0.0+0.0im  …  0.0+0.0im  0.0+0.0im

On such a coarse grid the transform error (absolute and relative) is about $10^{-2}$, this decreases for higher resolution. The gridded and spectral functions will choose a corresponding grid-spectral resolution (see Matching spectral and grid resolution) following quadratic truncation, but you can always truncate/interpolate in spectral space with spectral_truncation, spectral_interpolation which takes trunc = $l_{max} = m_{max}$ as second argument

spectral_truncation(alms, 2)
3×3 LowerTriangularMatrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im

Yay, we just chopped off $l > 2$ from alms which contained the harmonics up to degree and order 5 before. If the second argument in spectral_truncation is larger than alms then it will automatically call spectral_interpolation and vice versa. Also see Interpolation on RingGrids to interpolate directly between grids. If you want to control directly the resolution of the grid gridded is supposed to transform onto you have to provide a SpectralTransform instance. More on that now.

The SpectralTransform struct

Both spectral and gridded create an instance of SpectralTransform under the hood. This object contains all precomputed information that is required for the transform, either way: The Legendre polynomials, pre-planned Fourier transforms, precomputed gradient, divergence and curl operators, the spherical harmonic eigenvalues among others. Maybe the most intuitive way to create a SpectralTransform is to start with a SpectralGrid, which already defines which spectral resolution is supposed to be combined with a given grid.

using SpeedyWeather
spectral_grid = SpectralGrid(Float32, trunc=5, Grid=OctahedralGaussianGrid, dealiasing=3)
SpectralGrid:
├ Spectral:   T5 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid:       12-ring OctahedralGaussianGrid{Float32}, 360 grid points
├ Resolution: 1190km (average)
└ Vertical:   8-level SigmaCoordinates

(We using SpeedyWeather here as SpectralGrid is exported therein). We also specify the number format Float32 here to be used for the transform although this is the default anyway. From spectral_grid we now construct a SpectralTransform as follows

S = SpectralTransform(spectral_grid)
SpectralTransform{Float32}:
├ Spectral:   T5, 7x6 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       12-ring OctahedralGaussianGrid{Float32}
├ Truncation: dealiasing = 3 (cubic)
└ Legendre:   recompute polynomials false (1.22 KB)

Note that because we chose dealiasing=3 (cubic truncation) we now match a T5 spectral field with a 12-ring octahedral Gaussian grid, instead of the 8 rings as above. So going from dealiasing=2 (default) to dealiasing=3 increased our resolution on the grid while the spectral resolution remains the same. The SpectralTransform also has options for the recomputation or pre-computation of the Legendre polynomials, fore more information see (P)recompute Legendre polynomials.

Passing on S the SpectralTransform now allows us to transform directly on the grid defined therein.

map = gridded(alms, S)
plot(map)
       12-ring OctahedralGaussianGrid{Float32}     
       ┌────────────────────────┐  0.7
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └────────────────────────┘ -0.7
        0         ˚E         360

Yay, this is again the $l=m=1$ harmonic, but this time on a slightly higher resolution OctahedralGaussianGrid as specified in the SpectralTransform S. Note that also the number format was converted on the fly to Float32 because that is the number format we specified in S! And from grid to spectral

alms2 = spectral(map, S)
7×6 LowerTriangularMatrix{ComplexF32}:
   2.05564f-9+0.0im          0.0+0.0im          …         0.0+0.0im
          0.0+0.0im          1.0-2.14344f-8im             0.0+0.0im
   -1.3654f-9+0.0im          0.0+0.0im                    0.0+0.0im
          0.0+0.0im  -2.66275f-8-6.93167f-11im            0.0+0.0im
 -3.22443f-10+0.0im          0.0+0.0im                    0.0+0.0im
          0.0+0.0im  -4.58858f-8-8.19443f-9im   …  6.58519f-9+4.79414f-9im
   1.85859f-9+0.0im          0.0+0.0im                    0.0+0.0im

As you can see the rounding error is now more like $10^{-8}$ as we are using Float32 (the OctahedralGaussianGrid is another exact grid). Note, however, that the returned LowerTriangularMatrix is of size 7x6, not 6x6 what we started from. The underlying reason is that internally SpeedyWeather uses LowerTriangularMatrixs of size $l_{max} + 2 \times m_{max} + 1$. One $+1$ on both degree and order for 0-based harmonics versus 1-based matrix sizes, but an additional $+1$ for the degrees which is required by the meridional derivative. For consistency, all LowerTriangularMatrixs in SpeedyWeather.jl carry this additional degree but only the vector quantities explicitly make use of it. See Meridional derivative for details.

For this interface to SpeedyTransforms this means that on a grid-to-spectral transform you will get one more degree than orders of the spherical harmonics by default. You can, however, always truncate this additional degree, say to T5 (hence matrix size is 6x6)

spectral_truncation(alms2, 5, 5)
6×6 LowerTriangularMatrix{ComplexF32}:
   2.05564f-9+0.0im          0.0+0.0im          …         0.0+0.0im
          0.0+0.0im          1.0-2.14344f-8im             0.0+0.0im
   -1.3654f-9+0.0im          0.0+0.0im                    0.0+0.0im
          0.0+0.0im  -2.66275f-8-6.93167f-11im            0.0+0.0im
 -3.22443f-10+0.0im          0.0+0.0im                    0.0+0.0im
          0.0+0.0im  -4.58858f-8-8.19443f-9im   …  6.58519f-9+4.79414f-9im

spectral_truncation(alms2, 5) would have returned the same, a single argument is then assumed equal for both degrees and orders. Alternatively, you can also pass on the one_more_degree=false argument to the SpectralTransform constructor

S = SpectralTransform(spectral_grid, one_more_degree=false)
SpectralTransform{Float32}:
├ Spectral:   T5, 6x6 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       12-ring OctahedralGaussianGrid{Float32}
├ Truncation: dealiasing = 3 (cubic)
└ Legendre:   recompute polynomials false (1.07 KB)

As you can see the 7x6 LowerTriangularMatrix in the description above dropped down to 6x6 LowerTriangularMatrix, this is the size of the input that is expected (otherwise you will get a BoundsError).

Power spectrum

How to take some data and compute a power spectrum with SpeedyTransforms you may ask. Say you have some global data in a matrix m that looks, for example, like

m
96×47 Matrix{Float32}:
 7.35118  6.925    5.76322  4.14488  …  -0.684485   -1.21956   -1.16158
 7.33802  6.83566  5.58725  3.92778     -0.818366   -1.26571   -1.15598
 7.32059  6.73733  5.3982   3.69483     -0.951692   -1.30922   -1.14778
 7.29945  6.63187  5.19874  3.44815     -1.08086    -1.34837   -1.13663
 7.27521  6.52119  4.99173  3.19035     -1.20213    -1.38143   -1.12217
 7.24846  6.40724  4.7802   2.92455  …  -1.31168    -1.40673   -1.10407
 7.21979  6.29201  4.56733  2.65426     -1.40578    -1.42268   -1.08202
 7.18979  6.17745  4.35639  2.38331     -1.48088    -1.42779   -1.05573
 7.15904  6.06551  4.1507   2.11576     -1.53379    -1.42078   -1.02498
 7.12808  5.95806  3.9536   1.85578     -1.5618     -1.40053   -0.989571
 ⋮                                   ⋱               ⋮         
 7.20324  7.07189  6.39528  5.15423      0.0830802  -0.927487  -1.14657
 7.24704  7.1283   6.43223  5.14501      0.061883   -0.935579  -1.15128
 7.28361  7.1645   6.43852  5.10502      0.0210026  -0.951535  -1.15556
 7.31299  7.18108  6.41562  5.03697  …  -0.0383737  -0.974899  -1.15936
 7.33533  7.17884  6.36516  4.94323     -0.114805   -1.005     -1.16253
 7.35084  7.15876  6.28885  4.82571     -0.206605   -1.04095   -1.16492
 7.35979  7.12197  6.18852  4.68596     -0.31183    -1.08169   -1.1663
 7.36255  7.06977  6.0661   4.52517     -0.428269   -1.12598   -1.1664
 7.35953  7.0036   5.92361  4.34445  …  -0.55342    -1.17244   -1.16494

You hopefully know which grid this data comes on, let us assume it is a regular latitude-longitude grid, which we call the FullClenshawGrid. Note that for the spectral transform this should not include the poles, so the 96x47 matrix size here corresponds to

We now wrap this matrix therefore to associate it with the necessary grid information

map = FullClenshawGrid(m)
plot(map)
                      47-ring FullClenshawGrid{Float32}                   
       ┌────────────────────────────────────────────────────────────┐  10 
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └────────────────────────────────────────────────────────────┘ -8  
        0                           ˚E                           360

Now we transform into spectral space and call power_spectrum(::LowerTriangularMatrix)

alms = spectral(map)
power = SpeedyTransforms.power_spectrum(alms)

Which returns a vector of power at every wavenumber. By default this is normalized as average power per degree, you can change that with the keyword argument normalize=false. Plotting this yields

using UnicodePlots
k = 0:length(power)-1
lineplot(k, power, yscale=:log10, ylim=(1e-15, 10), xlim=extrema(k),
    ylabel="power", xlabel="wavenumber", height=10, width=60)
               ┌────────────────────────────────────────────────────────────┐ 
           10¹ │⠤⠤⢄⣀⡀⠀⣀⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠈⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢱⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
   power       │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣇⣀⣀⣀⣀⣀⣀⣀⣀⣠⣀⣀⣀⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⣀│ 
         10⁻¹⁵ │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│ 
               └────────────────────────────────────────────────────────────┘ 
               ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀31⠀ 
               ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀wavenumber⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀

The power spectrum of our data is about 1 up to wavenumber 10 and then close to zero for higher wavenumbers (which is in fact how we constructed this fake data). Let us turn this around and use SpeedyTransforms to create random noise in spectral space to be used in grid-point space!

Example: Creating k^n-distributed noise

How would we construct random noise in spectral space that follows a certain power law and transform it back into grid-point space? Define the wavenumber $k$ for T31, the spectral resolution we are interested in. (We start from 1 instead of 0 to avoid zero to the power of something negative). Now create some normally distributed spectral coefficients but scale them down for higher wavenumbers with $k^{-2}$

k = 1:32
alms = randn(LowerTriangularMatrix{Complex{Float32}}, 32, 32)
alms .*= k.^-2
32×32 LowerTriangularMatrix{ComplexF32}:
    0.0557274-1.0116im       …          0.0+0.0im
     0.248452-0.0583573im               0.0+0.0im
    -0.059435+0.119851im                0.0+0.0im
   0.00129244-0.0822943im               0.0+0.0im
   -0.0198537+0.0167453im               0.0+0.0im
   0.00467321+0.00395513im   …          0.0+0.0im
  0.000128168-0.00397645im              0.0+0.0im
  -0.00367459+0.00159226im              0.0+0.0im
   0.00369835-0.00992157im              0.0+0.0im
     0.013639+0.011443im                0.0+0.0im
             ⋮               ⋱  
 -0.000105731-0.00132082im              0.0+0.0im
 -0.000388599-0.00218143im              0.0+0.0im
   3.70741f-5-0.00105248im   …          0.0+0.0im
  -0.00209573+4.62188f-5im              0.0+0.0im
  0.000906028-0.00146972im              0.0+0.0im
 -0.000153108+0.000693383im             0.0+0.0im
  -0.00011329-0.000303503im             0.0+0.0im
 -0.000842798-0.000456834im  …          0.0+0.0im
 -0.000913645-0.000345001im     0.000435727+0.000786138im

Awesome. For higher degrees and orders the amplitude clearly decreases! Now to grid-point space and let us visualize the result

map = gridded(alms)
plot(map)
                      48-ring FullGaussianGrid{Float32}                   
       ┌────────────────────────────────────────────────────────────┐  0.4
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
       └────────────────────────────────────────────────────────────┘ -0.5
        0                           ˚E                           360

You can always access the underlying data in map via map.data in case you need to get rid of the wrapping into a grid again!

(P)recompute Legendre polynomials

The spectral transform uses a Legendre transform in meridional direction. For this the Legendre polynomials are required, at each latitude ring this is a $l_{max} \times m_{max}$ lower triangular matrix. Storing precomputed Legendre polynomials therefore quickly increase in size with resolution. One can recompute them to save memory, but that uses more arithmetic operations. There is therefore a memory-compute tradeoff.

For a single transform, there is no need to precompute the polynomials as the SpectralTransform object will be garbage collected again anyway. For low resolution simulations with many repeated small transforms it makes sense to precompute the polynomials and SpeedyWeather.jl does that automatically anyway. At very high resolution the polynomials may, however, become prohibitively large. An example at T127, about 100km resolution

spectral_grid = SpectralGrid(trunc=127)
SpectralTransform(spectral_grid, recompute_legendre=false)
SpectralTransform{Float32}:
├ Spectral:   T127, 129x128 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       192-ring OctahedralGaussianGrid{Float32}
├ Truncation: dealiasing = 2 (quadratic)
└ Legendre:   recompute polynomials false (3.23 MB)

the polynomials are about 3MB in size. Easy that is not much. But at T1023 on the O1536 octahedral Gaussian grid, this is already 1.5GB, cubically increasing with the spectral truncation T. recompute_legendre=true (default false when constructing a SpectralTransform object which may be reused) would lower this to kilobytes

SpectralTransform(spectral_grid, recompute_legendre=true)
SpectralTransform{Float32}:
├ Spectral:   T127, 129x128 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       192-ring OctahedralGaussianGrid{Float32}
├ Truncation: dealiasing = 2 (quadratic)
└ Legendre:   recompute polynomials true (33.60 KB)

Gradient operators

SpeedyTransforms also includes many gradient operators to take derivatives in spherical harmonics. These are in particular $\nabla, \nabla \cdot, \nabla \times, \nabla^2, \nabla^{-2}$. However, the actually implemented operators are, in contrast to the mathematical Derivatives in spherical coordinates due to reasons of scaling as follows. Let the implemented operators be $\hat{\nabla}$ etc.

\[\hat{\nabla} A = \left(\frac{\partial A}{\partial \lambda}, \cos(\theta)\frac{\partial A}{\partial \theta} \right) = R\cos(\theta)\nabla A\]

So the zonal derivative omits the radius and the $\cos^{-1}(\theta)$ scaling. The meridional derivative adds a $\cos(\theta)$ due to a recursion relation being defined that way, which, however, is actually convenient because the whole operator is therefore scaled by $R\cos(\theta)$. The curl and divergence operators expect the input velocity fields to be scaled by $\cos^{-1}(\theta)$, i.e.

\[\begin{aligned} \hat{\nabla} \cdot (\cos^{-1}(\theta)\mathbf{u}) &= \frac{\partial u}{\partial \lambda} + \cos\theta\frac{\partial v}{\partial \theta} = R\nabla \cdot \mathbf{u}, \\ \hat{\nabla} \times (\cos^{-1}(\theta)\mathbf{u}) &= \frac{\partial v}{\partial \lambda} - \cos\theta\frac{\partial u}{\partial \theta} = R\nabla \times \mathbf{u}. \end{aligned}\]

And the Laplace operators omit a $R^2$ (radius $R$) scaling, i.e.

\[\hat{\nabla}^{-2}A = \frac{1}{R^2}\nabla^{-2}A , \quad \hat{\nabla}^{2}A = R^2\nabla^{2}A\]

Example: Geostrophy

Now, we want to use the following example to illustrate their use: We have $u, v$ and want to calculate $\eta$ in the shallow water system from it following geostrophy. Analytically we have

\[-fv = -g\partial_\lambda \eta, \quad fu = -g\partial_\theta \eta\]

which becomes, if you take the divergence of these two equations

\[\zeta = \frac{g}{f}\nabla^2 \eta\]

Meaning that if we start with $u, v$ we can obtain the relative vorticity $\zeta$ and, using Coriolis parameter $f$ and gravity $g$, invert the Laplace operator to obtain displacement $\eta$. How to do this with SpeedyTransforms?

Let us start by generating some data

using SpeedyWeather

spectral_grid = SpectralGrid(trunc=31, nlev=1)
forcing = SpeedyWeather.JetStreamForcing(spectral_grid)
drag = QuadraticDrag(spectral_grid)
model = ShallowWaterModel(; spectral_grid, forcing, drag)
simulation = initialize!(model);
run!(simulation, period=Day(30))

Now pretend you only have u, v to get vorticity (which is actually the prognostic variable in the model, so calculated anyway...).

u = simulation.diagnostic_variables.layers[1].grid_variables.u_grid
v = simulation.diagnostic_variables.layers[1].grid_variables.v_grid
vor = SpeedyTransforms.curl(u, v) / spectral_grid.radius

Here, u, v are the grid-point velocity fields, and the function curl takes in either LowerTriangularMatrixs (no transform needed as all gradient operators act in spectral space), or, as shown here, arrays of the same grid and size. In this case, the function actually runs through the following steps

RingGrids.scale_coslat⁻¹!(u)
RingGrids.scale_coslat⁻¹!(v)

S = SpectralTransform(u, one_more_degree=true)
us = spectral(u, S)
vs = spectral(v, S)

vor = curl(us, vs)
33×32 LowerTriangularMatrix{ComplexF32}:
        0.0+0.0im         0.0+0.0im        …        0.0+0.0im
    19.8125+0.0im    0.304942+0.11123im             0.0+0.0im
    51.9837+0.0im     6.09093+1.10217im             0.0+0.0im
     210.51+0.0im     23.5655-11.0583im             0.0+0.0im
    107.815+0.0im     5.01137+14.7836im             0.0+0.0im
   -69.3147+0.0im     12.5271+7.19151im    …        0.0+0.0im
   -88.8129+0.0im  -0.0505834-46.8792im             0.0+0.0im
    -299.85+0.0im    -29.9557+25.7812im             0.0+0.0im
   -6.80406+0.0im       35.25-18.9829im             0.0+0.0im
    95.4603+0.0im    -1.95749+7.63725im             0.0+0.0im
           ⋮                               ⋱  
   0.483588+0.0im     3.20616-2.6084im              0.0+0.0im
 -0.0494366+0.0im    -2.33938+0.926119im   …        0.0+0.0im
  -0.144139+0.0im     1.72886-0.344334im            0.0+0.0im
  -0.117126+0.0im    0.331543+2.25201im             0.0+0.0im
   0.450046+0.0im    -0.83997+0.465839im            0.0+0.0im
   -1.36186+0.0im    0.519825-0.564119im            0.0+0.0im
   0.647759+0.0im   -0.387562-0.0816269im  …        0.0+0.0im
  -0.305471+0.0im   -0.022625-0.301176im      0.0869728-0.0291413im
        0.0+0.0im         0.0+0.0im                 0.0+0.0im

(Copies of) the velocity fields are unscaled by the cosine of latitude (see above), then transformed into spectral space, and the returned vor requires a manual division by the radius. We always unscale vector fields by the cosine of latitude before any curl, or div operation, as these omit those.

One more degree for spectral fields

The SpectralTransform in general takes a one_more_degree keyword argument, if otherwise the returned LowerTriangularMatrix would be of size 32x32, setting this to true would return 33x32. The reason is that while most people would expect square lower triangular matrices for a triangular spectral truncation, all vector quantities always need one more degree (= one more row) because of a recursion relation in the meridional gradient. So as we want to take the curl of us, vs here, they need this additional degree, but in the returned lower triangular matrix this row is set to zero.

One more degree for vector quantities

All gradient operators expect the input lower triangular matrices of shape $(N+1) \times N$. This one more degree of the spherical harmonics is required for the meridional derivative. Scalar quantities contain this degree too for size compatibility but they should not make use of it. Use spectral_truncation to add or remove this degree manually.

Example: Geostrophy (continued)

Now we transfer vor into grid-point space, but specify that we want it on the grid that we also used in spectral_grid. The Coriolis parameter for a grid like vor_grid is obtained, and we do the following for $f\zeta/g$.

vor_grid = gridded(vor, Grid=spectral_grid.Grid)
f = SpeedyWeather.coriolis(vor_grid)
fζ_g = spectral_grid.Grid(vor_grid .* f ./ model.planet.gravity)

The last line is a bit awkward for now, as the element-wise multiplication between two grids escapes the grid and returns a vector that we wrap again into a grid. We will fix that in future releases. Now we need to apply the inverse Laplace operator to $f\zeta/g$ which we do as follows

fζ_g_spectral = spectral(fζ_g, one_more_degree=true)
η = SpeedyTransforms.∇⁻²(fζ_g_spectral) * spectral_grid.radius^2
η_grid = gridded(η, Grid=spectral_grid.Grid)

Note the manual scaling with the radius $R^2$ here. We now compare the results

plot(η_grid)
                   48-ring OctahedralGaussianGrid{Float32}                 
       ┌────────────────────────────────────────────────────────────┐  5e⁹ 
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘ 
       └────────────────────────────────────────────────────────────┘ -1e¹⁰
        0                           ˚E                           360

Which is the interface displacement assuming geostrophy. The actual interface displacement contains also ageostrophy

plot(simulation.diagnostic_variables.surface.pres_grid)
                   48-ring OctahedralGaussianGrid{Float32}                 
       ┌────────────────────────────────────────────────────────────┐ 2 000
    90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
    ˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
       │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│ 
   -90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘ 
       └────────────────────────────────────────────────────────────┘ -600 
        0                           ˚E                           360

Strikingly similar! The remaining differences are the ageostrophic motions but also note that the mean is off. This is because geostrophy only use/defines the gradient of $\eta$ not the absolute values itself. Our geostrophic $\eta_g$ has by construction a mean of zero (that is how we define the inverse Laplace operator) but the actual $\eta$ is some 1400m higher.

Functions and type index

SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
S = SpectralTransform(  alms::AbstractMatrix{Complex{NF}};
                        recompute_legendre::Bool=true,
                        Grid::Type{<:AbstractGrid}=DEFAULT_GRID)

Generator function for a SpectralTransform struct based on the size of the spectral coefficients alms and the grid Grid. Recomputes the Legendre polynomials by default.

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
SpectralTransform(
    ::Type{NF},
    Grid::Type{<:SpeedyWeather.RingGrids.AbstractGrid},
    lmax::Int64,
    mmax::Int64;
    recompute_legendre,
    legendre_shortcut,
    dealiasing
) -> SpectralTransform

Generator function for a SpectralTransform struct. With NF the number format, Grid the grid type <:AbstractGrid and spectral truncation lmax, mmax this function sets up necessary constants for the spetral transform. Also plans the Fourier transforms, retrieves the colatitudes, and preallocates the Legendre polynomials (if recompute_legendre == false) and quadrature weights.

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
S = SpectralTransform(  map::AbstractGrid;
                        recompute_legendre::Bool=true)

Generator function for a SpectralTransform struct based on the size and grid type of gridded field map. Recomputes the Legendre polynomials by default.

source
SpeedyWeather.RingGrids.get_nlat_half — Function
get_nlat_half(trunc::Integer) -> Any
get_nlat_half(trunc::Integer, dealiasing::Real) -> Any

For the spectral truncation trunc (e.g. 31 for T31) return the grid resolution parameter nlat_half (number of latitude rings on one hemisphere including the Equator) following a dealiasing parameter (default 2) to match spectral and grid resolution.

source
SpeedyWeather.SpeedyTransforms.UV_from_vor! — Method
UV_from_vor!(
    U::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    V::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    vor::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    S::SpectralTransform{NF<:AbstractFloat}
)

Get U, V (=(u, v)*coslat) from vorticity ζ spectral space (divergence D=0) Two operations are combined into a single linear operation. First, invert the spherical Laplace ∇² operator to get stream function from vorticity. Then compute zonal and meridional gradients to get U, V.

source
SpeedyWeather.SpeedyTransforms.UV_from_vordiv! — Method
UV_from_vordiv!(
    U::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    V::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    vor::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    div::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    S::SpectralTransform{NF<:AbstractFloat}
) -> Complex

Get U, V (=(u, v)*coslat) from vorticity ζ and divergence D in spectral space. Two operations are combined into a single linear operation. First, invert the spherical Laplace ∇² operator to get stream function from vorticity and velocity potential from divergence. Then compute zonal and meridional gradients to get U, V.

source
SpeedyWeather.SpeedyTransforms._divergence! — Method
_divergence!(
    kernel,
    div::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    u::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    v::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    S::SpectralTransform{NF<:AbstractFloat}
)

Generic divergence function of vector u, v that writes into the output into div. Generic as it uses the kernel kernel such that curl, div, add or flipsign options are provided through kernel, but otherwise a single function is used.

source
SpeedyWeather.SpeedyTransforms.curl! — Method
curl!(
    curl::LowerTriangularMatrix,
    u::LowerTriangularMatrix,
    v::LowerTriangularMatrix,
    S::SpectralTransform;
    flipsign,
    add
)

Curl of a vector u, v written into curl, curl = ∇×(u, v). u, v are expected to have a 1/coslat-scaling included, then curl is not scaled. flipsign option calculates -∇×(u, v) instead. add option calculates curl += ∇×(u, v) instead. flipsign and add can be combined. This functions only creates the kernel and calls the generic divergence function _divergence! subsequently with flipped u, v -> v, u for the curl.

source
SpeedyWeather.SpeedyTransforms.curl — Method
curl(
    u::LowerTriangularMatrix,
    v::LowerTriangularMatrix
) -> Any

Curl (∇×) of two vector components u, v of size (n+1)xn, the last row will be set to zero in the returned LowerTriangularMatrix. This function requires both u, v to be transforms of fields that are scaled with 1/cos(lat). An example usage is therefore

RingGrids.scale_coslat⁻¹!(u_grid)
RingGrids.scale_coslat⁻¹!(v_grid)
u = spectral(u_grid)
v = spectral(v_grid)
vor = curl(u, v)
vor_grid = gridded(div)
source
SpeedyWeather.SpeedyTransforms.curl — Method
curl(
    u::SpeedyWeather.RingGrids.AbstractGrid,
    v::SpeedyWeather.RingGrids.AbstractGrid
) -> Any

Curl (∇×) of two vector components u, v on a grid. Applies 1/coslat scaling, transforms to spectral space and returns the spectral curl, which is scaled with the radius of the sphere. Divide by radius for unscaling.

source
SpeedyWeather.SpeedyTransforms.divergence! — Method
divergence!(
    div::LowerTriangularMatrix,
    u::LowerTriangularMatrix,
    v::LowerTriangularMatrix,
    S::SpectralTransform;
    flipsign,
    add
)

Divergence of a vector u, v written into div, div = ∇⋅(u, v). u, v are expected to have a 1/coslat-scaling included, then div is not scaled. flipsign option calculates -∇⋅(u, v) instead. add option calculates div += ∇⋅(u, v) instead. flipsign and add can be combined. This functions only creates the kernel and calls the generic divergence function _divergence! subsequently.

source
SpeedyWeather.SpeedyTransforms.divergence — Method
divergence(
    u::LowerTriangularMatrix,
    v::LowerTriangularMatrix
) -> Any

Divergence (∇⋅) of two vector components u, v which need to have size (n+1)xn, the last row will be set to zero in the returned LowerTriangularMatrix. This function requires both u, v to be transforms of fields that are scaled with 1/cos(lat). An example usage is therefore

RingGrids.scale_coslat⁻¹!(u_grid)
RingGrids.scale_coslat⁻¹!(v_grid)
u = spectral(u_grid, one_more_degree=true)
v = spectral(v_grid, one_more_degree=true)
div = divergence(u, v)
div_grid = gridded(div)

Both div and div_grid are scaled with the radius.

source
SpeedyWeather.SpeedyTransforms.divergence — Method
divergence(
    u::SpeedyWeather.RingGrids.AbstractGrid,
    v::SpeedyWeather.RingGrids.AbstractGrid
) -> Any

Divergence (∇⋅) of two vector components u, v on a grid. Applies 1/coslat scaling, transforms to spectral space and returns the spectral divergence, which is scaled with the radius of the sphere. Divide by radius for unscaling.

source
SpeedyWeather.SpeedyTransforms.get_recursion_factors — Method
get_recursion_factors(  ::Type{NF}, # number format NF
                        lmax::Int,  # max degree l of spherical harmonics (0-based here)
                        mmax::Int   # max order m of spherical harmonics
                        ) where {NF<:AbstractFloat}

Returns a matrix of recursion factors ϵ up to degree lmax and order mmax of number format NF.

source
SpeedyWeather.SpeedyTransforms.get_truncation — Function
get_truncation(nlat_half::Integer) -> Any
get_truncation(nlat_half::Integer, dealiasing::Real) -> Any

For the grid resolution parameter nlat_half (e.g. 24 for a 48-ring FullGaussianGrid) return the spectral truncation trunc (max degree of spherical harmonics) following a dealiasing parameter (default 2) to match spectral and grid resolution.

source
SpeedyWeather.SpeedyTransforms.gridded! — Method
gridded!(   map::AbstractGrid,
            alms::LowerTriangularMatrix,
            S::SpectralTransform)

Spectral transform (spectral to grid) of the spherical harmonic coefficients alms to a gridded field map. The spectral transform is number format-flexible as long as the parametric types of map, alms, S are identical. The spectral transform is grid-flexible as long as the typeof(map)<:AbstractGrid. Uses the precalculated arrays, FFT plans and other constants in the SpectralTransform struct S.

source
SpeedyWeather.SpeedyTransforms.gridded — Method
gridded(
    alms::AbstractArray{T<:Complex{NF}, 2};
    recompute_legendre,
    Grid,
    kwargs...
) -> Any

Spectral transform (spectral to grid space) from spherical coefficients alms to a newly allocated gridded field map. Based on the size of alms the grid type grid, the spatial resolution is retrieved based on the truncation defined for grid. SpectralTransform struct S is allocated to execute gridded(alms, S).

source
SpeedyWeather.SpeedyTransforms.gridded — Method
gridded(
    alms::AbstractMatrix,
    S::SpectralTransform{NF};
    kwargs...
) -> Any

Spectral transform (spectral to grid space) from spherical coefficients alms to a newly allocated gridded field map with precalculated properties based on the SpectralTransform struct S. alms is converted to a LowerTriangularMatrix to execute the in-place gridded!.

source
SpeedyWeather.SpeedyTransforms.roundup_fft — Method
m = roundup_fft(n::Int;
                small_primes::Vector{Int}=[2, 3, 5])

Returns an integer m >= n with only small prime factors 2, 3 (default, others can be specified with the keyword argument small_primes) to obtain an efficiently fourier-transformable number of longitudes, m = 2^i * 3^j * 5^k >= n, with i, j, k >=0.

source
SpeedyWeather.SpeedyTransforms.spectral! — Method
spectral!(  alms::LowerTriangularMatrix,
            map::AbstractGrid,
            S::SpectralTransform)

Spectral transform (grid to spectral space) from the gridded field map on a grid<:AbstractGrid to a LowerTriangularMatrix of spherical harmonic coefficients alms. Uses FFT in the zonal direction, and a Legendre Transform in the meridional direction exploiting symmetries. The spectral transform is number format-flexible as long as the parametric types of map, alms, S are identical. The spectral transform is grid-flexible as long as the typeof(map)<:AbstractGrid. Uses the precalculated arrays, FFT plans and other constants in the SpectralTransform struct S.

source
SpeedyWeather.SpeedyTransforms.spectral — Method
spectral(
    map::AbstractMatrix;
    Grid,
    kwargs...
) -> LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat

Converts map to grid(map) to execute spectral(map::AbstractGrid; kwargs...).

source
SpeedyWeather.SpeedyTransforms.spectral — Method
spectral(
    map::SpeedyWeather.RingGrids.AbstractGrid,
    S::SpectralTransform{NF}
) -> LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat

Spectral transform (grid to spectral) map to grid(map) to execute spectral(map::AbstractGrid; kwargs...).

source
SpeedyWeather.SpeedyTransforms.spectral — Method
spectral(
    map::SpeedyWeather.RingGrids.AbstractGrid{NF};
    recompute_legendre,
    one_more_degree
) -> LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat

Converts map to Grid(map) to execute spectral(map::AbstractGrid; kwargs...).

source
SpeedyWeather.SpeedyTransforms.spectral_interpolation — Method
alms_interp = spectral_interpolation(   ::Type{NF},
                                        alms::LowerTriangularMatrix,
                                        ltrunc::Integer,
                                        mtrunc::Integer
                                        ) where NF

Returns a spectral coefficient matrix alms_interp that is alms padded with zeros to interpolate in spectral space. If trunc is smaller or equal to the implicit truncation in alms obtained from its size than spectral_truncation is automatically called instead, returning alms_trunc, a coefficient matrix that is smaller than alms, implicitly setting higher degrees and orders to zero.

source
SpeedyWeather.SpeedyTransforms.spectral_smoothing! — Method
spectral_smoothing!(A::LowerTriangularMatrix, c; power=1)

Smooth the spectral field A following A = (1-(1-c)∇²ⁿ) with power n of a normalised Laplacian so that the highest degree lmax is dampened by multiplication with c. Anti-diffusion for c>1.

source
SpeedyWeather.SpeedyTransforms.spectral_smoothing — Method
A_smooth = spectral_smoothing(A::LowerTriangularMatrix, c; power=1)

Smooth the spectral field A following A_smooth = (1-c*∇²ⁿ)A with power n of a normalised Laplacian so that the highest degree lmax is dampened by multiplication with c. Anti-diffusion for c<0.

source
SpeedyWeather.SpeedyTransforms.spectral_truncation! — Method
spectral_truncation!(alms::AbstractMatrix, ltrunc::Integer, mtrunc::Integer)

Truncate spectral coefficients alms in-place by setting (a) the upper right triangle to zero and (b) all coefficients for which the degree l is larger than the truncation ltrunc or order m larger than the truncaction mtrunc.

source
SpeedyWeather.SpeedyTransforms.spectral_truncation! — Method
spectral_truncation!(alms::LowerTriangularMatrix, ltrunc::Integer, mtrunc::Integer)

Truncate spectral coefficients alms in-place by setting all coefficients for which the degree l is larger than the truncation ltrunc or order m larger than the truncaction mtrunc. Similar to spectral_truncation!(::AbstractMatrix, ...) but skips the upper triangle which is zero by design for LowerTriangularMatrix.

source
SpeedyWeather.SpeedyTransforms.spectral_truncation — Method
alms_trunc = spectral_truncation(alms, trunc)

Returns a spectral coefficient matrix alms_trunc that is truncated from alms to the size (trunc+1)². alms_trunc only contains those coefficient of alms for which m, l ≤ trunc, and l ≥ m are zero anyway. If trunc is larger than the implicit truncation in alms obtained from its size than spectral_interpolation is automatically called instead, returning alms_interp, a coefficient matrix that is larger than alms with padded zero coefficients.

source
SpeedyWeather.SpeedyTransforms.ϵlm — Method
ϵ = ϵ(l, m)

Recursion factors ϵ as a function of degree l and order m (0-based) of the spherical harmonics. ϵ(l, m) = sqrt((l^2-m^2)/(4*l^2-1)) with default number format Float64.

source
SpeedyWeather.SpeedyTransforms.ϵlm — Method
ϵ = ϵ(NF, l, m)

Recursion factors ϵ as a function of degree l and order m (0-based) of the spherical harmonics. ϵ(l, m) = sqrt((l^2-m^2)/(4*l^2-1)) and then converted to number format NF.

source
SpeedyWeather.SpeedyTransforms.∇! — Method
∇!(
    dpdx::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    dpdy::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    p::LowerTriangularMatrix{Complex{NF<:AbstractFloat}},
    S::SpectralTransform{NF<:AbstractFloat}
) -> Tuple{LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat, LowerTriangularMatrix{Complex{NF}} where NF<:AbstractFloat}

Applies the gradient operator ∇ applied to input p and stores the result in dpdx (zonal derivative) and dpdy (meridional derivative). The gradient operator acts on the unit sphere and therefore omits the radius scaling

source
SpeedyWeather.SpeedyTransforms.∇²! — Method
∇²!(    ∇²alms::LowerTriangularMatrix,
        alms::LowerTriangularMatrix,
        S::SpectralTransform;
        add::Bool=false,
        flipsign::Bool=false,
        inverse::Bool=false)

Laplace operator ∇² applied to the spectral coefficients alms in spherical coordinates. The radius R is omitted in the eigenvalues which are precomputed in S. ∇²! is the in-place version which directly stores the output in the first argument ∇²alms.

Keyword arguments

  • add=true adds the ∇²(alms) to the output
  • flipsign=true computes -∇²(alms) instead
  • inverse=true computes ∇⁻²(alms) instead

Default is add=false, flipsign=false, inverse=false. These options can be combined.

source
SpeedyWeather.SpeedyTransforms.∇² — Method
∇²(alms::LowerTriangularMatrix, S::SpectralTransform) -> Any

Laplace operator ∇² applied to input alms, using precomputed eigenvalues from S. The Laplace operator acts on the unit sphere and therefore omits the 1/radius^2 scaling

source
SpeedyWeather.SpeedyTransforms.∇² — Method
∇²(alms::LowerTriangularMatrix) -> Any

Returns the Laplace operator ∇² applied to input alms. The Laplace operator acts on the unit sphere and therefore omits the 1/radius^2 scaling

source
SpeedyWeather.SpeedyTransforms.∇⁻²! — Method
∇⁻²!(   ∇⁻²alms::LowerTriangularMatrix,
        alms::LowerTriangularMatrix,
        S::SpectralTransform;
        add::Bool=false,
        flipsign::Bool=false)

Calls ∇²!(∇⁻²alms, alms, S; add, flipsign, inverse=true).

source
SpeedyWeather.SpeedyTransforms.∇⁻² — Method
∇⁻²(
    ∇²alms::LowerTriangularMatrix,
    S::SpectralTransform
) -> Any

InverseLaplace operator ∇⁻² applied to input alms, using precomputed eigenvalues from S. The Laplace operator acts on the unit sphere and therefore omits the radius^2 scaling

source
SpeedyWeather.SpeedyTransforms.∇⁻² — Method
∇⁻²(∇²alms::LowerTriangularMatrix) -> Any

Returns the inverse Laplace operator ∇⁻² applied to input alms. The Laplace operator acts on the unit sphere and therefore omits the radius^2 scaling

source