Merge pull request #363 from SpeedyWeather/mk/docs

SpeedyTransforms documentation

Project.toml

@@ -44,10 +44,9 @@ JLD2 = "0.4"
 KernelAbstractions = "0.9"
 NetCDF = "0.10, 0.11"
 Primes = "0.5"
-ProgressMeter = "1.7"
-UnicodePlots = "3.3.2"
+ProgressMeter = "^1.7"
+UnicodePlots = "^3.3"
 julia = "1.8"
-
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 
 [compat]
 Documenter = "0.26, 0.27"
+UnicodePlots = "^3.3"

docs/src/grids.md

@@ -62,7 +62,7 @@ Currently the following full grids `<: AbstractFullGrid` are implemented
 - `FullClenshawGrid`, a full grid with equi-angle latitudes
 
 and additionally we have `FullHEALPixGrid` and `FullOctaHEALPixGrid` which are the full grid equivalents to the
-[HEALPix grid](@ref) and the [OctaHEALPix grid](@ref) discussed below. Full grid equivalent means that they have
+[HEALPix grid](@ref HEALPixGrid) and the [OctaHEALPix grid](@ref OctaHEALPixGrid) discussed below. Full grid equivalent means that they have
 the same latitude rings, but no reduction in the number of points per ring towards the poles and no longitude offset.
 Other implemented reduced grids are
 
@@ -144,11 +144,11 @@ You will obtain this information every time you create a `SpectralGrid(;Grid,tru
 
 ...
 
-## Octahedral Gaussian grid
+## [Octahedral Gaussian grid](@id OctahedralGaussianGrid)
 
 ...
 
-## HEALPix grid
+## [HEALPix grid](@id HEALPixGrid)
 
 Technically, HEALPix grids are a class of grids that tessalate the sphere into faces that are often
 called basepixels. For each member of this class there are ``N_\varphi`` basepixels in zonal direction
@@ -197,7 +197,7 @@ z = \frac{2}{3}-\frac{4k}{3N_{side}} \pm \frac{8\phi}{3\pi}
 ```
 
 
-## OctaHEALPix grid
+## [OctaHEALPix grid](@id OctaHEALPixGrid)
 
 While the classic HEALPix grid is based on a [dodecahedron](https://en.wikipedia.org/wiki/Rhombic_dodecahedron),
 other choices for ``N_\varphi`` and ``N_\theta`` in the class of HEALPix grids will change the number of faces

docs/src/speedytransforms.md

@@ -2,33 +2,253 @@
 
 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.
-
-### Example transforms
-
-```julia
-julia> using SpeedyWeather
-julia> alms = zeros(ComplexF64,3,3)    # spectral coefficients
-julia> alms[2,2] = 1                   # only l=1,m=1 harmonic
-julia> map = gridded(alms)             # convert to grid space
-8×4 Matrix{Float64}:
- -0.324541  -0.600363  -0.600363  -0.324541
- -0.134429  -0.248678  -0.248678  -0.134429
-  0.134429   0.248678   0.248678   0.134429
-  0.324541   0.600363   0.600363   0.324541
-  0.324541   0.600363   0.600363   0.324541
-  0.134429   0.248678   0.248678   0.134429
- -0.134429  -0.248678  -0.248678  -0.134429
- -0.324541  -0.600363  -0.600363  -0.324541
-
-julia> spectral(map)                   # back to spectral space
-3×3 Matrix{ComplexF64}:
- 0.0+0.0im  0.0+0.0im          0.0+0.0im
- 0.0+0.0im  1.0+3.60727e-17im  0.0+0.0im
- 0.0+0.0im  0.0+0.0im          0.0+0.0im
-```
-
-and we have successfully reobtained the ``l=m=1`` spherical harmonic.
+It is just not put into its own respective repository for now.
+
+The SpeedyTransforms are based on [RingGrids](@ref) and
+[LowerTriangularMatrices](@ref 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](@ref) here, but will focus on the usage of the
+SpeedyTransforms module.
+
+## 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
+
+```@example speedytransforms
+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.
+
+```@example speedytransforms
+alms = zeros(LowerTriangularMatrix{ComplexF64},6,6)     # spectral coefficients
+alms[2,2] = 1                                           # only l=1,m=1 harmonic
+alms
+```
+Now `gridded` is the function that takes spectral coefficients `alms` and converts
+them a grid-point space `map`
+
+```@example speedytransforms
+map = gridded(alms)
+```
+By default, the `gridded` transforms onto a [`FullGaussianGrid`](@ref FullGaussianGrid) unravelled here
+into a vector west to east, starting at the prime meridian, then north to south, see [RingGrids](@ref).
+We can visualize `map` quickly with a unicodeplot via `plot` (see [Visualising RingGrid data](@ref))
+```@example speedytransforms
+plot(map)
+```
+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`
+```@example speedytransforms
+alms2 = spectral(map)
+```
+Comparing with `alms` from above you can see that the transform is exact up to a typical rounding error
+from `Float64`. 
+```@example speedytransforms
+alms ≈ alms2
+```
+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`](@ref 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](@ref) is the [`FullGaussianGrid`](@ref FullGaussianGrid)
+we can transform onto other grids by specifying `Grid` too
+```@example speedytransforms
+map = gridded(alms,Grid=HEALPixGrid)
+plot(map)
+```
+which, if transformed back, however, can yield a larger transform error as discussed above
+```@example speedytransforms
+spectral(map)
+```
+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](@ref)) 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
+```@example speedytransforms
+spectral_truncation(alms,2)
+```
+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](@ref)
+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.
+```@example speedytransforms
+using SpeedyWeather
+spectral_grid = SpectralGrid(Float32,trunc=5,Grid=OctahedralGaussianGrid,dealiasing=3)
+```
+(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
+```@example speedytransforms
+S = SpectralTransform(spectral_grid)
+```
+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 precomputation of the Legendre polynomials, fore more information see
+[(P)recompute Legendre polynomials](@ref).
+
+Passing on `S` the `SpectralTransform` now allows us to transform directly on the grid
+defined therein.
+```@example speedytransforms
+map = gridded(alms,S)
+plot(map)
+```
+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
+```@example speedytransforms
+alms2 = spectral(map,S)
+```
+As you can see the rounding error is now more like ``10^{-8}`` as we are using Float32
+(the [OctahedralGaussianGrid](@ref 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
+`LowerTriangularMatrix`s 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 `LowerTriangularMatrix`s in SpeedyWeather.jl carry this additional degree
+but only the vector quantities explicitly make use of it. 
+See [Meridional derivative](@ref) 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)
+```@example speedytransforms
+spectral_truncation(alms2,5,5)
+```
+`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
+
+```@example speedytransforms
+S = SpectralTransform(spectral_grid,one_more_degree=false)
+```
+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
+```@example speedytransforms
+alms = randn(LowerTriangularMatrix{Complex{Float32}},32,32) # hide
+spectral_truncation!(alms,10) # hide
+map = gridded(alms,Grid=FullClenshawGrid) # hide
+m = Matrix(map) # hide
+m
+```
+You hopefully know which grid this data comes on, let us assume it is a regular
+latitde-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
+```@example speedytransforms
+map = FullClenshawGrid(m)
+plot(map)
+```
+Now we transform into spectral space and call `power_spectrum(::LowerTriangularMatrix)`
+```@example speedytransforms
+alms = spectral(map)
+power = SpeedyTransforms.power_spectrum(alms)
+nothing # hide
+```
+
+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
+```@example speedytransforms
+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)
+```
+
+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}``
+
+```@example speedytransforms
+k = 1:32
+alms = randn(LowerTriangularMatrix{Complex{Float32}},32,32)
+alms .*= k.^-2
+```
+Awesome. For higher degrees and orders the amplitude clearly decreases! Now
+to grid-point space and let us visualize the result
+```@example speedytransforms
+map = gridded(alms)
+plot(map)
+```
+
+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
+
+```@example speedytransforms
+spectral_grid = SpectralGrid(trunc=127)
+SpectralTransform(spectral_grid,recompute_legendre=false)
+```
+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
+```@example speedytransforms
+SpectralTransform(spectral_grid,recompute_legendre=true)
+```
+
 
 ## Functions and type index
 

src/RingGrids/RingGrids.jl

@@ -70,4 +70,5 @@ include("octahealpix.jl")
 include("quadrature_weights.jl")
 include("interpolation.jl")
 include("show.jl")
+include("similar.jl")
 end

src/RingGrids/similar.jl

@@ -0,0 +1,24 @@
+# todo this should be replaced by general broadcasting for <: AbstractGrid
+for Grid in (:FullGaussianGrid, :FullClenshawGrid, :FullHEALPixGrid, :FullOctaHEALPixGrid,
+            :OctahedralGaussianGrid, :OctahedralClenshawGrid, :HEALPixGrid,
+            :OctaHEALPixGrid)
+    @eval begin
+
+        function Base.similar(G::$Grid)
+            return $Grid{eltype(G)}(undef,G.nlat_half)
+        end
+
+        function Base.similar(G::$Grid,::Type{T}) where T
+            return $Grid{T}(undef,G.nlat_half)
+        end
+
+        function Base.similar(G::$Grid,nlat_half::Integer)
+            return $Grid{eltype(G)}(undef,nlat_half)
+        end
+
+        function Base.similar(G::$Grid,::Type{T},nlat_half::Integer) where T
+            return $Grid{T}(undef,nlat_half)
+        end
+
+    end
+end