Grid-flexible spectral transform

[milankl]

This PR is intended to fix #112 by allowing for a flexible grid<:AbstractGrid to be used in the spectral transform.

[milankl]

@hottad at the moment this PR lets one use any grid<:AbstractGrid without errors, but the transforms are not accurate yet for Octahedral/HEALPix. I assume I got some normalisation in the grid->spectral transform wrong. For a FullGaussianGrid we have

julia> alms = LowerTriangularMatrix(randn(ComplexF64,33,32));
julia> alms[:,1] .= real.(alms[:,1]);
julia> map = gridded(alms,grid=FullGaussianGrid);
julia> alms2 = spectral(map);
julia> imshow(log10.(abs.(alms-alms2)))

an absolute error that is O(1e-14) with Float64. That's good.

with FullClenshawGrid similar, that's great too, except for the last row. I assume this is because we use lmax+2,mmax+1 LowerTriangularMatrices for the vector quantities, but so that's technically T32 which violates the exactness of the Clenshaw-Curtis quadrature with 64 latitudes. I'm not worried about this as the last row is only included for consistent matrix sizes (and the issue disappears if one starts with the lmax+2 row being zeros).

However, for Octahedral/HEALPix the absolute errors are O(1e-2) and larger. I believe this is because I messed up the normalisation with the solid angle $\Delta \Omega = sin(\theta) \Delta \theta \Delta \phi$ that a grid point covers. For the full Gaussian grids I use (96x48)

julia> nlat_half = 24
julia> w = SpeedyWeather.get_quadrature_weights(FullGaussianGrid,nlat_half);
julia> w *= π/2nlat_half
24-element Vector{Float64}:
 0.0002063860165031381
 0.0004795872813548503
 0.0007511832466177125
 ...
 0.0041838316317469116
 0.00421930059577275
 0.004237072348271895

so quadrature weights *pi/nlat and the same for the Clenshaw-Curtis grid but with other weights obviously. I believe this has to be changed for Octahedral/HEALPix. I'll try to figure this out in the next days.

[milankl]

@hottad With the latest commit the transform errors on the HEALPix grid are much smaller

from

julia> alms = LowerTriangularMatrix(randn(ComplexF64,33,32))
julia> alms[end,:] .= 0
julia> alms[:,1] .= real.(alms[:,1])
julia> map = gridded(alms,grid=HEALPixGrid)
julia> alms2 = spectral(map)
julia> maximum(abs.(alms)-abs.(alms2))
0.018819099369166703

julia> imshow(abs.(alms-alms2))

but only when the lon offset rotation is deactivated

[milankl]

These are the best quadrature weights for the HEALPix grid I was sofar able to come up with.

[milankl]

However, the longitude offset rotation had to be deactivated. In general we do here

    _, lons = get_colatlons(grid,nresolution)
    lon1s = [lons[each_index_in_ring(grid,i,nresolution)[1]] for i in 1:nlat]
    lon_offsets = [cispi(m*lon1/π) for m in 0:mmax, lon1 in lon1s]

with lon1s the longitude in radians of the first grid point per ring.

[milankl]

@hottad While I'm still not sure about the accuracy of the HEALPix spectral transform, I'm tempted to merge this once all tests pass as all the other structure has been set up such that you can start contributing/experimenting and do whatever you'd like to test with new grids in SpeedyWeather.jl, I think we probably just ran the first HEALPix atmospheric simulation, see below? FullGaussianGrid will remain default for now and we can independently improve the transform and support more grids. For now I'd suggest to keep all grid-related stuff in src/grids.jl but eventually we could even go to src/Grids/grids.jl, healpix.jl, octahedral.jl or outsurce all of that into its won package. But I think that's an overkill for now.

[hottad]

@milankl Fantastic!
Thank you for the great job! So shallow water on FllClenshawGrid, OctahedralGaussianGrid and HEALPIxGrid now runs up to at least 10 days and the results are all visually very similar. Impressive!
I will fork this project and start experimenting probably next week. Congrats!

[milankl]

Yeah, the HEALPixGrid still doesn't have the correct longitude offset rotation implemented, that's something that I need to play around with. As you suggest, it's probably some implementation bug that I don't currently understand as the mathematics are rather clear I believe. However, it's still good enough to run even 1000 days or longer and also at different resolution.