-
Notifications
You must be signed in to change notification settings - Fork 24
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
A new grid for SpeedyWeather.jl? #112
Comments
Also useful https://arxiv.org/abs/astro-ph/0409513 / https://iopscience.iop.org/article/10.1086/427976/pdf
We can create a dependency on Healpix.jl and use their functions and the data structure Nside = 64
H = HealpixMap{Float32,RingOrder}(Nside) to hold the data in grid point space. We could use the following resolutions |
After playing around with creating different grids, I believe we'll always have the Fourier constraints that at the equator
Then the number of grid points for HEALPix is always 2/3 of the regular Gaussian grid. The resolution of a HEALPix grid could be estimated via
|
Following Hotta 2018 one should be able to swap out the Legendre weights we currently use with the Clenshaw-Curtis quadrature that uses the following weights with J being the number of latitudes so that's as simple as julia> nlat = 32
julia> θjs = [j/(nlat+1)*π for j in 1:nlat]
32-element Vector{Float64}:
0.09519977738150888
0.19039955476301776
0.28559933214452665
⋮
2.855993321445266
2.9511930988267756
3.046392876208284
julia> wj = [4sin(θj)/(nlat+1)*sum([sin(p*θj)/p for p in 1:2:nlat]) for θj in θjs]
32-element Vector{Float64}:
0.010663416953504472
0.01628798008667748
0.028546546759225664
⋮
0.02854654675922572
0.016287980086677478
0.010663416953504501
julia> sum(wj)
2.0000000000000004 |
Started the implementation in #122 |
#124 implements further grid abstractions towards a single spectral transform that can support all of the grids discussed above |
Daisuke Hotta @hottad was proposing HEALPix4x1 grid, that uses 4 base pixels, no equatorial zone, from Gorski's paper the pixel positions then become It has at most |
Summing this up, I think if we want cubic truncation with T31 then these would be good alternative candidates, F32 (128x64 grid points) the current default but higher resolution for cubic, O32 (64 latitude rings, 144 longitude points around the equator), H32 (63 latitude rings, 128 points around the equator) Cool thing is, H32 would have even fewer grid points than F24 what we currently use for T31. In terms of the FFT complexity, which is |
We could also, without any interpolation needed, store the data on a HEALPix4 grid directly ncview-able into a square matrix like so |
With #127 the spectral transform in T31 is timed with for grid in (FullGaussianGrid,FullClenshawGrid,OctahedralGaussianGrid,HEALPixGrid)
map = gridded(alms;grid)
S = SpeedyWeather.SpectralTransform(map,recompute_legendre=false)
@btime SpeedyWeather.gridded!($map,$alms,$S)
@btime SpeedyWeather.spectral!($alms,$map,$S)
end as (all Float64)
(EDIT: Now with shortened loops over the zonal wavenumber Note that the new grids are compared to F32 (so the current non-default cubic version of linear F24 we normally use with T31). HEALPix is only cubic in latitude, but linear in longitude. But sanity check: A single spectral transform can be fast for a range of grids ✅ |
Great! |
just a little overview about the characteristics of the latitudinal resolution of the different grids / quadratures with 1024 latitude rings (Gaussian) or 1023 (Clenshaw-Curtis, HEALPix).
|
Looping over rings was changed with #132 |
At the moment we use a regular Gaussian grid for grid-point space. This comes with advantages of easy storage in a matrix internally and for output but with the disadvantage that a lot of points are put near the poles that aren't actually needed. Alternatives are
This is an issue to discuss and collect ideas and reasons to implement either of these grids at some point in the future.
The text was updated successfully, but these errors were encountered: