Skip to content
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

Dealiasing parameter for quad,cub truncation #290

Merged
merged 8 commits into from
Apr 7, 2023
Merged

Dealiasing parameter for quad,cub truncation #290

merged 8 commits into from
Apr 7, 2023

Conversation

milankl
Copy link
Member

@milankl milankl commented Apr 6, 2023

@hottad this introduces dealiasing as a parameter (1=linear, 2=quadratic, 3=cubic) but you can also chose anything in between and it chooses the right grid size given the truncation trunc (and dealiasing). I therefore removed anything that related to the spectral resolution from RingGrids and SpeedyTransforms now has the following methods

const DEFAULT_DEALIASING = 2.0    # quadratic truncation default

function RingGrids.get_nlat_half(   trunc::Integer,
                                    dealiasing::Real=DEFAULT_DEALIASING)
                    
    return roundup_fft(ceil(Int,((1+dealiasing)*trunc+1)/4))
end

function get_truncation(nlat_half::Integer,
                        dealiasing::Real=DEFAULT_DEALIASING)

    return floor(Int,(4nlat_half-1)/(1+dealiasing))
end

to match spectral and grid-point resolution as desired.

@milankl milankl added performance 🚀 Faster faster! transform ⬅️ ➡️ Our spherical harmonic transform, grid to spectral, spectral to grid labels Apr 6, 2023
@milankl
Copy link
Member Author

milankl commented Apr 6, 2023

OctaHEALPixGrid is 1.6x faster than FullGaussianGrid when both are used at quadratic truncation:

image

There are visible differences, does it matter? I don't know, we will see. The OctahedralGaussianGrid is still exact at quadratic so no visible differences but considerably slower than OctaHEALPix

image

@milankl
Copy link
Member Author

milankl commented Apr 6, 2023

For some reason this still throws an error in the first spectral transform

julia> run_speedy(PrimitiveDryCore,Grid=HEALPixGrid,dealiasing=2)
ERROR: BoundsError
Stacktrace:
 [1] gridded!(map::HEALPixGrid{Float64}, alms::LowerTriangularMatrix{ComplexF64}, S::SpectralTransform{Float64}; unscale_coslat::Bool)
   @ Main.SpeedyWeather.SpeedyTransforms ~/git/SpeedyWeather.jl/src/SpeedyTransforms/spectral_transform.jl:355
 [2] gridded!
   @ ~/git/SpeedyWeather.jl/src/SpeedyTransforms/spectral_transform.jl:340 [inlined]
 [3] initialize_orography!(orog::Main.SpeedyWeather.Orography{Float64, HEALPixGrid{Float64}}, #unused#::EarthOrography, P::Parameters{PrimitiveDryCore}, S::SpectralTransform{Float64}, G::Geometry{Float64})

@milankl milankl added this to the v0.5 milestone Apr 6, 2023
@milankl
Copy link
Member Author

milankl commented Apr 7, 2023

The last commit resolves the issue that HEALPix errored with quadratic truncation. However, there are clearly issues with the HEALPixGrid when using quadratic truncation: This is T63 with

  • (top) O48, OctahedralGaussian, 96 rings, 208 points at Equator
  • (middle) P48, OctaHEALPix, 95 rings, 192 points at Equator
  • (bottom) H48, HEALPix, 95 rings, 96 points at Equator

image

Note all the aliasing with HEALPix, which however goes basically away with OctaHEALPix! 🎉

@milankl milankl merged commit 1409af8 into main Apr 7, 2023
@hottad
Copy link
Contributor

hottad commented Apr 9, 2023

Great work! Thank you, Milan!

Meridional quadrature in HEALPixGrid and OctaHEALPixGrid are inexact by design, so it is natural that you see visibly different results with HEALPixGrid at a relatively low resolution of T63.
It will be interesting to compare at a much higher resolution since the inexact quadrature in HEALPixGrid and OctaHEALPixGrid becomes more and more accurate with increased resolution. From experience with ectrans (IFS's spectral transform code), the "inexactness" of the quadrature becomes irrelevant at around T479 resolution, in the sense that single-precision transforms on these grids become as accurate as single-precision transforms on the tried-and-truth Gaussian or Clenshaw-Curtis grids.

I found it surprising that T63 on H48 did not blow up (albeit with easily recognizable aliasing noises)!
With T63, purely zonal wave can have zonal wavenumber of 63, which requires (63+1)*2=128 zonal grid points along the Equator even to avoid linear aliasing, while H48 only has 96 points along the tropical belt.

It's great that this problem goes away with OctaHEALPix. It is exactly for this reason that I suggested OctaHEALPix, so it is just consistent with the theory, but it is always nice to see that an implementation actually works as expected from theory.

Thank you! and Congratulations!

@milankl milankl deleted the mk/truncation branch April 21, 2023 20:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
performance 🚀 Faster faster! transform ⬅️ ➡️ Our spherical harmonic transform, grid to spectral, spectral to grid
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants