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.

Notation: Spectral resolution

There are different ways to describe the spectral resolution, the truncation wavenumber (e.g. T31), the maximum degree $l$ and order $m$ of the spherical harmonics (e.g. $l_{max}=31$, $m_{max} = 31$), or the size of the lower triangular matrix, e.g. 32x32. In this example, they are all equivalent. We often use the truncation, i.e. T31, for brevity but sometimes it is important to describe degree and order independently (see for example One more degree for spectral fields). Note also how truncation, degree and order are 0-based, but matrix sizes are 1-based.

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. Note, the $+1$ on both degree (first index) and order (second index) for 0-based harmonics versus 1-based matrix indexing, see Size of LowerTriangularArray. Create a LowerTriangularMatrix for T5 resolution, i.e. 6x6 matrix size

alms = zeros(LowerTriangularMatrix{ComplexF64}, 6, 6)     # spectral coefficients T5
alms[2, 2] = 1                                            # only l=1, m=1 harmonic
alms
21-element, 6x6 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 transform is the function that takes spectral coefficients alms and converts them a grid-point space map (or vice versa)

map = transform(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 transforms 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 transform

alms2 = transform(map)
21-element, 6x6 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 = transform(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

transform(map)
21-element, 6x6 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.48537e-17im  0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im   3.42299e-17-1.51822e-34im  0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im     0.0185292-1.7445e-17im   0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im  -6.28844e-17+2.78914e-34im  0.0+0.0im     0.0+0.0im  0.0+0.0im
 0.0+0.0im    -0.0164851+5.53707e-18im  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 transform function 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)
6-element, 3x3 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 you want to transform onto, use the keyword dealiasing (default: 2 for quadratic, see Matching spectral and grid resolution). But you can also provide a SpectralTransform instance to reuse a precomputed spectral transform. More on that now.

The SpectralTransform struct

The function transform only with arguments as shown above, will 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(NF=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-layer SigmaCoordinates
└ Device:     CPU using Array

(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, Array}:
├ Spectral:   T5, 7x6 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       12-ring OctahedralGaussianArray{Float32}
├ Truncation: dealiasing = 3 (cubic)
├ Legendre:   Polynomials 720 bytes, shortcut: linear
└ Memory:     for 8 layers (18.94 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.

Passing on S the SpectralTransform now allows us to transform directly on the grid defined therein. Note that we recreate alms to be of size 7x6 instead of 6x6 for T5 spectral resolution because SpeedyWeather uses internally One more degree for spectral fields meaning also that's the default when creating a SpectralTransform from a SpectralGrid. But results don't change if the last degree (row) contains only zeros.

alms = zeros(LowerTriangularMatrix{ComplexF64}, 7, 6)     # spectral coefficients
alms[2, 2] = 1                                            # only l=1, m=1 harmonic

map = transform(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 = transform(map, S)
27-element, 7x6 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). While 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)
21-element, 6x6 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, Array}:
├ Spectral:   T5, 6x6 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       12-ring OctahedralGaussianArray{Float32}
├ Truncation: dealiasing = 3 (cubic)
├ Legendre:   Polynomials 576 bytes, shortcut: linear
└ Memory:     for 8 layers (18.94 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).

SpectralTransform generators

While you can always create a SpectralTransform from a SpectralGrid (which defines both spectral and grid space) there are other constructors/generators available:

SpectralTransform(alms)
SpectralTransform{Float64, Array}:
├ Spectral:   T5, 7x6 LowerTriangularMatrix{Complex{Float64}}
├ Grid:       8-ring FullGaussianArray{Float64}
├ Truncation: dealiasing = 1.67 (linear)
├ Legendre:   Polynomials 936 bytes, shortcut: linear
└ Memory:     for 1 layers (1.62 KB)

Now we have defined the resolution of the spectral space through alms but create a SpectralTransform by making assumption about the grid space. E.g. Grid=FullGaussianGrid by default, dealiasing=2 and nlat_half correspondingly. But you can also pass them on as keyword arguments, for example

SpectralTransform(alms, Grid=OctahedralClenshawGrid, nlat_half=24)
SpectralTransform{Float64, Array}:
├ Spectral:   T5, 7x6 LowerTriangularMatrix{Complex{Float64}}
├ Grid:       47-ring OctahedralClenshawArray{Float64}
├ Truncation: dealiasing = 15 (>cubic)
├ Legendre:   Polynomials 5.26 KB, shortcut: linear
└ Memory:     for 1 layers (45.78 KB)

Only note that you don't want to specify both nlat_half and dealiasing as you would otherwise overspecify the grid resolution (dealiasing will be ignored in that case). This also works starting from the grid space

grid = rand(FullClenshawGrid, 12)
SpectralTransform(grid)
SpectralTransform{Float64, Array}:
├ Spectral:   T15, 16x16 LowerTriangularMatrix{Complex{Float64}}
├ Grid:       23-ring FullClenshawArray{Float64}
├ Truncation: dealiasing = 2 (quadratic)
├ Legendre:   Polynomials 13.13 KB, shortcut: linear
└ Memory:     for 1 layers (10.58 KB)

where you can also provide spectral resolution trunc or dealiasing. You can also provide both a grid and a lower triangular matrix to describe both spaces

SpectralTransform(grid, alms)
SpectralTransform{Float64, Array}:
├ Spectral:   T5, 7x6 LowerTriangularMatrix{Complex{Float64}}
├ Grid:       23-ring FullClenshawArray{Float64}
├ Truncation: dealiasing = 7 (>cubic)
├ Legendre:   Polynomials 2.66 KB, shortcut: linear
└ Memory:     for 1 layers (10.58 KB)

and you will precompute the transform between them. For more generators see the docstrings at ?SpectralTransform.

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}:
 -2.30471  -4.58838  -5.88513  -6.02969  …  2.13554  1.56549   0.462376
 -2.25735  -4.49369  -5.74758  -5.86137     2.14974  1.58759   0.471419
 -2.19786  -4.37567  -5.5771   -5.65183     2.16231  1.60132   0.474548
 -2.12658  -4.23539  -5.37615  -5.40592     2.17527  1.60686   0.471705
 -2.04389  -4.07409  -5.14761  -5.1293      2.18982  1.60421   0.462843
 -1.95024  -3.89313  -4.89455  -4.82807  …  2.20625  1.59315   0.447935
 -1.84612  -3.69395  -4.62017  -4.50832     2.22373  1.57325   0.42697
 -1.73207  -3.47806  -4.32768  -4.17586     2.24032  1.54387   0.39996
 -1.60865  -3.24698  -4.02011  -3.83582     2.25297  1.5042    0.366939
 -1.47647  -3.00223  -3.70027  -3.49243     2.25765  1.45328   0.327976
  ⋮                                      ⋱           ⋮         
 -2.16397  -4.31884  -5.49803  -5.50151     1.48697  0.924264  0.134974
 -2.22942  -4.44687  -5.67924  -5.71998     1.63243  1.03862   0.190968
 -2.28286  -4.55173  -5.82936  -5.90673     1.75631  1.14199   0.242646
 -2.32401  -4.63255  -5.94633  -6.05752  …  1.85914  1.23422   0.289701
 -2.35262  -4.68861  -6.0284   -6.16829     1.94208  1.31534   0.331866
 -2.36854  -4.71942  -6.07419  -6.23545     2.00688  1.38555   0.368911
 -2.37168  -4.72471  -6.08281  -6.2562      2.05581  1.44522   0.400638
 -2.36204  -4.70445  -6.0539   -6.22884     2.09157  1.49477   0.42688
 -2.33967  -4.65886  -5.98772  -6.153    …  2.11712  1.5347    0.447498

You hopefully know which grid this data comes on, let us assume it is a regular latitude-longitude grid, which we call the FullClenshawGrid (in analogy to the Gaussian grid based on the Gaussian quadrature). Note that for the spectral transform this should not include the poles, so the 96x47 matrix size here corresponds to 23 latitudes north and south of the Equator respectively plus the equator (=47).

We now wrap this matrix into a FullClenshawGrid (input_as=Matrix is required because all grids organise their data as vectors, see Creating data on a RingGrid) therefore to associate it with the necessary grid information like its coordinates

map = FullClenshawGrid(m, input_as=Matrix)

using CairoMakie
heatmap(map)

Random pattern

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

alms = transform(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
A = randn(Complex{Float32}, 32, 32)
A .*= k.^-2
alms = LowerTriangularArray(A)
528-element, 32x32 LowerTriangularMatrix{ComplexF32}
    -0.142038-0.993244im     …           0.0+0.0im
    -0.112602-0.103906im                 0.0+0.0im
     0.110879-0.0703427im                0.0+0.0im
   -0.0176798+0.0320909im                0.0+0.0im
    0.0245083+0.0514916im                0.0+0.0im
   0.00919612+0.0162823im    …           0.0+0.0im
   0.00710186+0.0137954im                0.0+0.0im
    0.0152184+0.000884688im              0.0+0.0im
  -0.00217564-0.009574im                 0.0+0.0im
  -0.00485983+0.0030024im                0.0+0.0im
             ⋮               ⋱  
  0.000124134+0.0004233im                0.0+0.0im
  -0.00133846-0.000632424im              0.0+0.0im
 -0.000324178-0.00294804im   …           0.0+0.0im
  0.000853275+0.000554057im              0.0+0.0im
  0.000503763-0.00194255im               0.0+0.0im
   0.00035071-0.000752746im              0.0+0.0im
   4.31725f-5-0.000380091im              0.0+0.0im
 -0.000595802-0.000639392im  …           0.0+0.0im
   4.03396f-5+0.000337018im     -0.000195414-0.0013991im

We first create a Julia Matrix so that the matrix-vector broadcasting .*= k is correctly applied across dimensions of A and then convert to a LowerTriangularMatrix.

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

map = transform(alms)

using CairoMakie
heatmap(map, title="k⁻²-distributed noise")

Random noise

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!

Precomputed polynomials and allocated memory

Reuse `SpectralTransform`s wherever possible

Depending on horizontal and vertical resolution of spectral and grid space, a SpectralTransform can be become very large in memory. Also the recomputation of the polynomials and the planning of the FFTs are expensive compared to the actual transform itself. Therefore reuse a SpectralTransform wherever possible.

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. It is therefore advised to reuse a precomputed SpectralTransform object wherever possible. This prevents transforms to allocate large memory which would otherwise be garbage collected again after the transform.

You get information about the size of that memory (both polynomials and required scratch memory) in the terminal "show" of a SpectralTransform object, e.g. at T127 resolution with 8 layers these are

spectral_grid = SpectralGrid(trunc=127, nlayers=8)
SpectralTransform(spectral_grid)
SpectralTransform{Float32, Array}:
├ Spectral:   T127, 129x128 LowerTriangularMatrix{Complex{Float32}}
├ Grid:       192-ring OctahedralGaussianArray{Float32}
├ Truncation: dealiasing = 2 (quadratic)
├ Legendre:   Polynomials 3.22 MB, shortcut: linear
└ Memory:     for 8 layers (2.50 MB)

Batched Transforms

SpeedyTransforms also supports batched transforms. With batched input data the transform is performed along the leading dimension, and all further dimensions are interpreted as batch dimensions. Take for example 

alms = randn(LowerTriangularMatrix{Complex{Float32}}, 32, 32, 5)
grids = transform(alms)
4608×5, 48-ring FullGaussianArray{Float32, 2, Matrix{Float32}}:
 -5.72033    -10.1227    4.39881    4.1058   -20.7772
 -4.67266    -10.0137    4.33609    3.70019  -20.6815
 -3.56206     -9.86772   4.22917    3.32544  -20.5448
 -2.39658     -9.68413   4.08117    2.98336  -20.3637
 -1.18628     -9.46202   3.8956     2.67527  -20.134
  0.0567691   -9.19991   3.67614    2.40201  -19.8509
  1.31865     -8.89588   3.42662    2.16386  -19.5096
  2.58382     -8.54761   3.15081    1.96055  -19.1048
  3.8354      -8.1525    2.85239    1.79119  -18.6315
  5.0556      -7.70788   2.53486    1.65433  -18.0849
  ⋮                                          
  3.23457     -1.31531   1.76262   -6.57055    1.18011
  3.27744     -0.56633   1.46828   -6.78818    1.02957
  3.30253      0.199651  1.20894   -6.98588    0.915881
  3.30962      0.979607  0.987442  -7.15571    0.838799
  3.29818      1.77067   0.806369  -7.29061    0.79669
  3.2674       2.57013   0.667953  -7.38461    0.786722
  3.21626      3.37547   0.574025  -7.43296    0.805047
  3.14367      4.18431   0.525958  -7.43225    0.847034
  3.04854      4.99437   0.524624  -7.38041    0.907517

In this case we first randomly generated five (32x32) LowerTriangularMatrix that hold the coefficients and then transformed all five matrices batched to the grid space with the transform command, yielding 5 RingGrids with each 48-rings.

Functions and type index

SpeedyWeather.SpeedyTransforms.AbstractLegendreShortcut — Type

Legendre shortcut is the truncation of the loop over the order m of the spherical harmonics in the Legendre transform. For reduced grids with as few as 4 longitudes around the poles (HEALPix grids) or 20 (octahedral grids) the higher wavenumbers in large orders m do not project (significantly) onto such few longitudes. For performance reasons the loop over m can therefore but truncated early. A Legendre shortcut <: AbstractLegendreShortcut is implemented as a functor that returns the 0-based maximum order m to retain per latitude ring, i.e. to be used for m in 0:mmax_truncation.

New shortcuts can be added by defining struct LegendreShortcutNew <: AbstractLegendreShortcut end and the functor method LegendreShortcutNew(nlon::Integer, lat::Real) = ..., with nlon the number of longitude points on that ring, and latd its latitude in degrees (-90˚ to 90˚N). Many implementations may not use the latitude latd but it is included for compatibility. If unused set to default value to 0. Also define short_name(::Type{<:LegendreShortcutNew}) = "new".

Implementions are LegendreShortcutLinear, LegendreShortcutQuadratic, LegendreShortcutCubic, LegendreShortcutLinQuadCosLat² and LegendreShortcutLinCubCoslat.

source
SpeedyWeather.SpeedyTransforms.AssociatedLegendrePolArray — Type
AssociatedLegendrePolArray{T, N, M, V} <: AbstractArray{T,N}

Type that wraps around a LowerTriangularArray{T,M,V} but is a subtype of AbstractArray{T,M+1}. This enables easier use with AssociatedLegendrePolynomials.jl which otherwise couldn't use the "matrix-style" (l, m) indexing of LowerTriangularArray. This type however doesn't support any other operations than indexing and is purerly intended for internal purposes.

  • data::LowerTriangularArray{T, M, V} where {T, M, V}
source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Type

SpectralTransform struct that contains all parameters and precomputed arrays to perform a spectral transform. Fields are

  • Grid::Type{<:AbstractGridArray}

  • nlat_half::Int64

  • nlayers::Int64

  • lmax::Int64

  • mmax::Int64

  • nfreq_max::Int64

  • LegendreShortcut::Type{<:SpeedyWeather.SpeedyTransforms.AbstractLegendreShortcut}

  • mmax_truncation::Any

  • nlon_max::Int64

  • nlons::Any

  • nlat::Int64

  • coslat::Any

  • coslat⁻¹::Any

  • lon_offsets::Any

  • norm_sphere::Any

  • rfft_plans::Vector{AbstractFFTs.Plan}

  • brfft_plans::Vector{AbstractFFTs.Plan}

  • rfft_plans_1D::Vector{AbstractFFTs.Plan}

  • brfft_plans_1D::Vector{AbstractFFTs.Plan}

  • legendre_polynomials::Any

  • scratch_memory_north::Any

  • scratch_memory_south::Any

  • scratch_memory_grid::Any

  • scratch_memory_spec::Any

  • scratch_memory_column_north::Any

  • scratch_memory_column_south::Any

  • solid_angles::Any

  • grad_y1::Any

  • grad_y2::Any

  • grad_y_vordiv1::Any

  • grad_y_vordiv2::Any

  • vordiv_to_uv_x::Any

  • vordiv_to_uv1::Any

  • vordiv_to_uv2::Any

  • eigenvalues::Any

  • eigenvalues⁻¹::Any

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
SpectralTransform(
    grids::AbstractGridArray{NF, N, ArrayType};
    trunc,
    dealiasing,
    one_more_degree,
    kwargs...
) -> SpectralTransform{NF, _A, _B, _C, _D, _E, _F, LowerTriangularArray{NF1, 1, _A1}, LowerTriangularArray{NF2, 2, _A2}} where {NF, _A, _B, _C, _D, _E, _F, NF1, _A1, NF2, _A2}

Generator function for a SpectralTransform struct based on the size and grid type of grids. Use keyword arugments trunc, dealiasing (ignored if trunc is used) or one_more_degree to define the spectral truncation.

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
SpectralTransform(
    grids::AbstractGridArray{NF1, N, ArrayType1},
    specs::LowerTriangularArray{NF2, N, ArrayType2};
    kwargs...
) -> SpectralTransform{NF, _A, _B, _C, _D, _E, _F, LowerTriangularArray{NF1, 1, _A1}, LowerTriangularArray{NF2, 2, _A2}} where {NF, _A, _B, _C, _D, _E, _F, NF1, _A1, NF2, _A2}

Generator function for a SpectralTransform struct to transform between grids and specs.

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
SpectralTransform(
    specs::LowerTriangularArray{NF, N, ArrayType};
    nlat_half,
    dealiasing,
    kwargs...
) -> SpectralTransform{NF, _A, _B, _C, _D, _E, _F, LowerTriangularArray{NF1, 1, _A1}, LowerTriangularArray{NF2, 2, _A2}} where {NF, _A, _B, _C, _D, _E, _F, NF1, _A1, NF2, _A2}

Generator function for a SpectralTransform struct based on the size of the spectral coefficients specs. Use keyword arguments nlat_half, Grid or deliasing (if nlat_half not provided) to define the grid.

source
SpeedyWeather.SpeedyTransforms.SpectralTransform — Method
SpectralTransform(
    ::Type{NF},
    lmax::Integer,
    mmax::Integer,
    nlat_half::Integer;
    Grid,
    ArrayType,
    nlayers,
    LegendreShortcut
) -> SpectralTransform{T, Array, Vector{T1}, Array{Complex{T2}, 1}, Vector{Int64}, Array{Complex{T3}, 2}, Array{Complex{T4}, 3}, LowerTriangularArray{T5, 1, Vector{T6}}, LowerTriangularArray{T7, 2, Matrix{T8}}} where {T, T1, T2, T3, T4, T5, T6, T7, T8}

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 and quadrature weights.

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::LowerTriangularArray,
    V::LowerTriangularArray,
    vor::LowerTriangularArray,
    S::SpectralTransform;
    radius
) -> Tuple{LowerTriangularArray, LowerTriangularArray}

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. Acts on the unit sphere, i.e. it omits any radius scaling as all inplace gradient operators, unless the radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.UV_from_vordiv! — Method
UV_from_vordiv!(
    U::LowerTriangularArray,
    V::LowerTriangularArray,
    vor::LowerTriangularArray,
    div::LowerTriangularArray,
    S::SpectralTransform;
    radius
) -> Tuple{LowerTriangularArray, LowerTriangularArray}

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. Acts on the unit sphere, i.e. it omits any radius scaling as all inplace gradient operators.

source
SpeedyWeather.SpeedyTransforms._divergence! — Method
_divergence!(
    kernel,
    div::LowerTriangularArray,
    u::LowerTriangularArray,
    v::LowerTriangularArray,
    S::SpectralTransform;
    radius
) -> LowerTriangularArray

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. Acts on the unit sphere, i.e. it omits 1/radius scaling as all gradient operators, unless the radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms._fourier_batched! — Method
_fourier_batched!(
    f_north::AbstractArray{<:Complex, 3},
    f_south::AbstractArray{<:Complex, 3},
    grids::AbstractGridArray,
    S::SpectralTransform
)

(Forward) Fast Fourier transform (grid to spectral) in zonal direction of grids, stored in scratch memories f_north, f_south to be passed on to the Legendre transform. Batched version that requires the number of vertical layers to be the same as precomputed in S. Not to be called directly, use transform! instead.

source
SpeedyWeather.SpeedyTransforms._fourier_batched! — Method
_fourier_batched!(
    grids::AbstractGridArray,
    g_north::AbstractArray{<:Complex, 3},
    g_south::AbstractArray{<:Complex, 3},
    S::SpectralTransform
)

Inverse fast Fourier transform (spectral to grid) of Legendre-transformed inputs g_north and g_south to be stored in grids. Not to be called directly, use transform! instead.

source
SpeedyWeather.SpeedyTransforms._fourier_serial! — Method
_fourier_serial!(
    f_north::AbstractArray{<:Complex, 3},
    f_south::AbstractArray{<:Complex, 3},
    grids::AbstractGridArray,
    S::SpectralTransform
)

(Forward) Fast Fourier transform (grid to spectral) in zonal direction of grids, stored in scratch memories f_north, f_south to be passed on to the Legendre transform. Serial version that does not require the number of vertical layers to be the same as precomputed in S. Not to be called directly, use transform! instead.

source
SpeedyWeather.SpeedyTransforms._fourier_serial! — Method
_fourier_serial!(
    grids::AbstractGridArray,
    g_north::AbstractArray{<:Complex, 3},
    g_south::AbstractArray{<:Complex, 3},
    S::SpectralTransform
)

(Inverse) Fast Fourier transform (spectral to grid) of Legendre-transformed inputs g_north and g_south to be stored in grids. Serial version that does not require the number of vertical layers to be the same as precomputed in S. Not to be called directly, use transform! instead.

source
SpeedyWeather.SpeedyTransforms._legendre! — Method
_legendre!(
    g_north::AbstractArray{<:Complex, 3},
    g_south::AbstractArray{<:Complex, 3},
    specs::LowerTriangularArray,
    S::SpectralTransform;
    unscale_coslat
)

Inverse Legendre transform, batched in the vertical. Not to be used directly, but called from transform!.

source
SpeedyWeather.SpeedyTransforms._legendre! — Method
_legendre!(
    specs::LowerTriangularArray,
    f_north::AbstractArray{<:Complex, 3},
    f_south::AbstractArray{<:Complex, 3},
    S::SpectralTransform
)

(Forward) Legendre transform, batched in the vertical. Not to be used directly, but called from transform!.

source
SpeedyWeather.SpeedyTransforms.curl! — Method
curl!(
    curl::LowerTriangularArray,
    u::LowerTriangularArray,
    v::LowerTriangularArray,
    S::SpectralTransform;
    flipsign,
    add,
    kwargs...
) -> LowerTriangularArray

Curl of a vector u, v written into curl, curl = ∇×(u, v). u, v are expected to have a 1/coslat-scaling included, otherwise curl is scaled. Acts on the unit sphere, i.e. it omits 1/radius scaling as all gradient operators unless the radius keyword argument is provided. 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::LowerTriangularArray,
    v::LowerTriangularArray;
    kwargs...
) -> 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). Acts on the unit sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided. An example usage is therefore

RingGrids.scale_coslat⁻¹!(u_grid)
RingGrids.scale_coslat⁻¹!(v_grid)
u = transform(u_grid)
v = transform(v_grid)
vor = curl(u, v, radius=6.371e6)
vor_grid = transform(div)
source
SpeedyWeather.SpeedyTransforms.curl — Method
curl(
    u::AbstractGridArray,
    v::AbstractGridArray;
    kwargs...
) -> Any

Curl (∇×) of two vector components u, v on a grid. Applies 1/coslat scaling, transforms to spectral space and returns the spectral curl. Acts on the unit sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided.

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

Divergence of a vector u, v written into div, div = ∇⋅(u, v). u, v are expected to have a 1/coslat-scaling included, otherwise div is scaled. Acts on the unit sphere, i.e. it omits 1/radius scaling as all gradient operators, unless the radius keyword argument is provided. 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::LowerTriangularArray,
    v::LowerTriangularArray;
    kwargs...
) -> LowerTriangularArray

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). Acts on the unit sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided. An example usage is therefore

RingGrids.scale_coslat⁻¹!(u_grid)
RingGrids.scale_coslat⁻¹!(v_grid)
u = transform(u_grid, one_more_degree=true)
v = transform(v_grid, one_more_degree=true)
div = divergence(u, v, radius = 6.371e6)
div_grid = transform(div)
source
SpeedyWeather.SpeedyTransforms.divergence — Method
divergence(
    u::AbstractGridArray,
    v::AbstractGridArray;
    kwargs...
) -> Any

Divergence (∇⋅) of two vector components u, v on a grid. Applies 1/coslat scaling, transforms to spectral space and returns the spectral divergence. Acts on the unit sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided.

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.ismatching — Method
ismatching(
    S::SpectralTransform,
    grid::AbstractGridArray
) -> Any

Spectral transform S and grid match if the resolution nlat_half and the type of the grid match and the number of vertical layers is equal or larger in the transform (constraints due to allocated scratch memory size).

source
SpeedyWeather.SpeedyTransforms.ismatching — Method
ismatching(
    S::SpectralTransform,
    L::LowerTriangularArray
) -> Any

Spectral transform S and lower triangular matrix L match if the spectral dimensions (lmax, mmax) match and the number of vertical layers is equal or larger in the transform (constraints due to allocated scratch memory size).

source
SpeedyWeather.SpeedyTransforms.power_spectrum — Method
power_spectrum(
    alms::LowerTriangularArray{Complex{NF}, 1, Array{Complex{NF}, 1}};
    normalize
) -> Any

Compute the power spectrum of the spherical harmonic coefficients alms (lower triangular matrix) of type Complex{NF}.

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_interpolation — Method
spectral_interpolation(
    _::Type{NF},
    alms::LowerTriangularArray{T, N, ArrayType},
    ltrunc::Integer,
    mtrunc::Integer
) -> Any

Returns a LowerTriangularArray that is interpolated from alms to the size (ltrunc+1) x (mtrunc+1), both inputs are 0-based, by padding zeros for higher wavenumbers. If ltrunc or mtrunc are smaller than the corresponding size ofalms than spectral_truncation is automatically called instead, returning a smaller LowerTriangularArray.

source
SpeedyWeather.SpeedyTransforms.spectral_smoothing! — Method
spectral_smoothing!(
    L::LowerTriangularArray,
    c::Real;
    power,
    truncation
)

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
spectral_smoothing(
    A::LowerTriangularArray,
    c::Real;
    power
) -> Any

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::LowerTriangularArray,
    ltrunc::Integer,
    mtrunc::Integer
) -> LowerTriangularArray

Triangular truncation to degree ltrunc and order mtrunc (both 0-based). 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.

source
SpeedyWeather.SpeedyTransforms.spectral_truncation — Method
spectral_truncation(
    _::Type{NF},
    alms::LowerTriangularArray{T, N, ArrayType},
    ltrunc::Integer,
    mtrunc::Integer
) -> Any

Returns a LowerTriangularArray that is truncated from alms to the size (ltrunc+1) x (mtrunc+1), both inputs are 0-based. If ltrunc or mtrunc is larger than the corresponding size ofalms than spectral_interpolation is automatically called instead, returning a LowerTriangularArray padded zero coefficients for higher wavenumbers.

source
SpeedyWeather.SpeedyTransforms.transform! — Method
transform!(
    grids::AbstractGridArray,
    specs::LowerTriangularArray,
    S::SpectralTransform;
    unscale_coslat
) -> AbstractGridArray

Spectral transform (spectral to grid space) from n-dimensional array specs of spherical harmonic coefficients to an n-dimensional array grids of ring grids. Uses FFT in the zonal direction, and a Legendre Transform in the meridional direction exploiting symmetries. The spectral transform is number format-flexible but grids and the spectral transform S have to have the same number format. Uses the precalculated arrays, FFT plans and other constants in the SpectralTransform struct S. The spectral transform is grid-flexible as long as the typeof(grids)<:AbstractGridArray and S.Grid matches.

source
SpeedyWeather.SpeedyTransforms.transform! — Method
transform!(
    specs::LowerTriangularArray,
    grids::AbstractGridArray,
    S::SpectralTransform
) -> LowerTriangularArray

Spectral transform (grid to spectral space) from n-dimensional array of grids to an n-dimensional array specs of spherical harmonic coefficients. Uses FFT in the zonal direction, and a Legendre Transform in the meridional direction exploiting symmetries. The spectral transform is number format-flexible but grids and the spectral transform S have to have the same number format. Uses the precalculated arrays, FFT plans and other constants in the SpectralTransform struct S. The spectral transform is grid-flexible as long as the typeof(grids)<:AbstractGridArray and S.Grid matches.

source
SpeedyWeather.SpeedyTransforms.transform — Method
transform(grids::AbstractGridArray; kwargs...) -> Any

Spectral transform (grid to spectral space) from grids to a newly allocated LowerTriangularArray. Based on the size of grids and the keyword dealiasing the spectral resolution trunc is retrieved. SpectralTransform struct S is allocated to execute transform(grids, S).

source
SpeedyWeather.SpeedyTransforms.transform — Method
transform(
    specs::LowerTriangularArray;
    unscale_coslat,
    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 transform(alms, S).

source
SpeedyWeather.SpeedyTransforms.transform — Method
transform(
    grids::AbstractGridArray,
    S::SpectralTransform{NF}
) -> Any

Spherical harmonic transform from grids to a newly allocated specs::LowerTriangularArray using the precomputed spectral transform S.

source
SpeedyWeather.SpeedyTransforms.transform — Method
transform(
    specs::LowerTriangularArray,
    S::SpectralTransform{NF};
    kwargs...
) -> Any

Spherical harmonic transform from specs to a newly allocated grids::AbstractGridArray using the precomputed spectral transform S.

source
SpeedyWeather.SpeedyTransforms.∇! — Method
∇!(
    dpdx::LowerTriangularArray,
    dpdy::LowerTriangularArray,
    p::LowerTriangularArray,
    S::SpectralTransform;
    radius
) -> Tuple{LowerTriangularArray, LowerTriangularArray}

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 1/radius scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇ — Method
∇(
    grid::AbstractGridArray,
    S::SpectralTransform;
    kwargs...
) -> Tuple{Any, Any}

The zonal and meridional gradient of grid. Transform to spectral space, takes the gradient and unscales the 1/coslat scaling in the gradient. Acts on the unit-sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided. Makes use of an existing spectral transform S.

source
SpeedyWeather.SpeedyTransforms.∇ — Method
∇(grid::AbstractGridArray; kwargs...) -> Tuple{Any, Any}

The zonal and meridional gradient of grid. Transform to spectral space, takes the gradient and unscales the 1/coslat scaling in the gradient. Acts on the unit-sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇ — Method
∇(
    p::LowerTriangularArray,
    S::SpectralTransform;
    kwargs...
) -> Tuple{Any, Any}

The zonal and meridional gradient of p using an existing SpectralTransform S. Acts on the unit sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇ — Method
∇(p::LowerTriangularArray; kwargs...) -> Tuple{Any, Any}

The zonal and meridional gradient of p. Precomputes a SpectralTransform S. Acts on the unit-sphere, i.e. it omits 1/radius scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇²! — Method
∇²!(
    ∇²alms::LowerTriangularArray,
    alms::LowerTriangularArray,
    S::SpectralTransform;
    add,
    flipsign,
    inverse,
    radius
) -> LowerTriangularArray

Laplace operator ∇² applied to the spectral coefficients alms in spherical coordinates. The eigenvalues which are precomputed in S. ∇²! is the in-place version which directly stores the output in the first argument ∇²alms. Acts on the unit sphere, i.e. it omits any radius scaling as all inplace gradient operators, unless the radius keyword argument is provided.

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::LowerTriangularArray,
    S::SpectralTransform;
    kwargs...
) -> Any

Laplace operator ∇² applied to input alms, using precomputed eigenvalues from S. Acts on the unit sphere, i.e. it omits 1/radius^2 scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇² — Method
∇²(alms::LowerTriangularArray; kwargs...) -> Any

Returns the Laplace operator ∇² applied to input alms. Acts on the unit sphere, i.e. it omits 1/radius^2 scaling unless radius keyword argument is provided.

source
SpeedyWeather.SpeedyTransforms.∇⁻²! — Method
∇⁻²!(
    ∇⁻²alms::LowerTriangularArray,
    alms::LowerTriangularArray,
    S::SpectralTransform;
    add,
    flipsign,
    kwargs...
) -> LowerTriangularArray

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

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

InverseLaplace operator ∇⁻² applied to input alms, using precomputed eigenvalues from S. Acts on the unit sphere, i.e. it omits radius^2 scaling unless radius keyword argument is provided.

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

Returns the inverse Laplace operator ∇⁻² applied to input alms. Acts on the unit sphere, i.e. it omits radius^2 scaling unless radius keyword argument is provided.

source