diff --git a/src/RingGrids/RingGrids.jl b/src/RingGrids/RingGrids.jl index f63d9e6f3..fa3dcb340 100644 --- a/src/RingGrids/RingGrids.jl +++ b/src/RingGrids/RingGrids.jl @@ -2,7 +2,6 @@ module RingGrids import Parameters: @unpack import Statistics: mean - import Primes import FastGaussQuadrature # GRIDS diff --git a/src/RingGrids/full_grids.jl b/src/RingGrids/full_grids.jl index 28d3b1636..aadb476c8 100644 --- a/src/RingGrids/full_grids.jl +++ b/src/RingGrids/full_grids.jl @@ -102,9 +102,6 @@ nlat_half_clenshaw(npoints::Integer) = round(Int,1/4 + sqrt(1/16 + npoints/8)) FullClenshawGrid{T}(data::AbstractVector) where T = FullClenshawGrid{T}(data,nlat_half_clenshaw(length(data))) FullClenshawGrid(data::AbstractVector,n::Integer...) = FullClenshawGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:FullClenshawGrid}) = 3 # cubic -get_truncation(::Type{<:FullClenshawGrid},nlat_half::Integer) = floor(Int,(4nlat_half-1)/4) -get_resolution(::Type{<:FullClenshawGrid},trunc::Integer) = roundup_fft(ceil(Int,(4*trunc+1)/4)) nlat_odd(::Type{<:FullClenshawGrid}) = true get_npoints(::Type{<:FullClenshawGrid},nlat_half::Integer) = npoints_clenshaw(nlat_half) get_colat(::Type{<:FullClenshawGrid},nlat_half::Integer) = [j/(2nlat_half)*π for j in 1:2nlat_half-1] @@ -134,9 +131,6 @@ nlat_half_gaussian(npoints::Integer) = round(Int,sqrt(npoints/8)) FullGaussianGrid{T}(data::AbstractVector) where T = FullGaussianGrid{T}(data,nlat_half_gaussian(length(data))) FullGaussianGrid(data::AbstractVector,n::Integer...) = FullGaussianGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:FullGaussianGrid}) = 2 # quadratic -get_truncation(::Type{<:FullGaussianGrid},nlat_half::Integer) = floor(Int,(4nlat_half-1)/3) -get_resolution(::Type{<:FullGaussianGrid},trunc::Integer) = roundup_fft(ceil(Int,(3*trunc+1)/4)) nlat_odd(::Type{<:FullGaussianGrid}) = false get_npoints(::Type{<:FullGaussianGrid},nlat_half::Integer) = npoints_gaussian(nlat_half) get_colat(::Type{<:FullGaussianGrid},nlat_half::Integer) = @@ -166,9 +160,6 @@ nlat_half_fullhealpix(npoints::Integer) = round(Int,1/4 + sqrt(1/16 + npoints/8) # infer nlat_half from data vector length, infer parametric type from eltype of data FullHEALPixGrid{T}(data::AbstractVector) where T = FullHEALPixGrid{T}(data,nlat_half_fullhealpix(length(data))) FullHEALPixGrid(data::AbstractVector,n::Integer...) = FullHEALPixGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:FullHEALPixGrid}) = 2 # quadratic -get_truncation(::Type{<:FullHEALPixGrid},nlat_half::Integer) = get_truncation(HEALPixGrid,nlat_half) -get_resolution(::Type{<:FullHEALPixGrid},trunc::Integer) = get_resolution(HEALPixGrid,trunc) nlat_odd(::Type{<:FullHEALPixGrid}) = true get_npoints(::Type{<:FullHEALPixGrid},nlat_half::Integer) = npoints_fullhealpix(nlat_half) get_colat(::Type{<:FullHEALPixGrid},nlat_half::Integer) = get_colat(HEALPixGrid,nlat_half) @@ -197,9 +188,6 @@ nlat_half_fulloctahealpix(npoints::Integer) = round(Int,1/4 + sqrt(1/16 + npoint FullOctaHEALPixGrid{T}(data::AbstractVector) where T = FullOctaHEALPixGrid{T}(data,nlat_half_fulloctahealpix(length(data))) FullOctaHEALPixGrid(data::AbstractVector,n::Integer...) = FullOctaHEALPixGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:FullOctaHEALPixGrid}) = 3 # cubic -get_truncation(::Type{<:FullOctaHEALPixGrid},nlat_half::Integer) = floor(Int,(4nlat_half-1)/4) -get_resolution(::Type{<:FullOctaHEALPixGrid},trunc::Integer) = roundup_fft(ceil(Int,(4*trunc+1)/4)) nlat_odd(::Type{<:FullOctaHEALPixGrid}) = true get_npoints(::Type{<:FullOctaHEALPixGrid},nlat_half::Integer) = npoints_fulloctahealpix(nlat_half) get_colat(::Type{<:FullOctaHEALPixGrid},nlat_half::Integer) = get_colat(OctaHEALPixGrid,nlat_half) diff --git a/src/RingGrids/grids_general.jl b/src/RingGrids/grids_general.jl index 95b83d2a5..d2ec9eecd 100644 --- a/src/RingGrids/grids_general.jl +++ b/src/RingGrids/grids_general.jl @@ -105,10 +105,8 @@ Base.rand(::Type{Grid},nlat_half::Integer) where {Grid<:AbstractGrid{T}} where T Grid(rand(T,get_npoints(Grid,nlat_half)),nlat_half) # truncation is the spectral truncation corresponding to size of grid and lin/quad/cubic truncation -get_truncation(grid::Grid) where {Grid<:AbstractGrid} = get_truncation(Grid,grid.nlat_half) get_resolution(grid::AbstractGrid) = get_nlat_half(grid) get_nlat_half(grid::AbstractGrid) = grid.nlat_half -get_nlat_half(G::Type{<:AbstractGrid},trunc::Integer) = get_resolution(G,trunc) # does the grid have an odd number of latitudes? nlat_odd(grid::AbstractGrid) = nlat_odd(typeof(grid)) diff --git a/src/RingGrids/healpix.jl b/src/RingGrids/healpix.jl index eac7e84a2..3c78783e1 100644 --- a/src/RingGrids/healpix.jl +++ b/src/RingGrids/healpix.jl @@ -127,10 +127,6 @@ end HEALPixGrid{T}(data::AbstractVector) where T = HEALPixGrid{T}(data,nlat_half_healpix(length(data))) HEALPixGrid(data::AbstractVector,n::Integer...) = HEALPixGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:HEALPixGrid}) = 1 # linear (in longitude) -get_truncation(::Type{<:HEALPixGrid},nlat_half::Integer) = nlat_half-1 -get_resolution(::Type{<:HEALPixGrid},trunc::Integer) = trunc+1 - function get_colat(::Type{<:HEALPixGrid},nlat_half::Integer) nlat_half == 0 && return Float64[] diff --git a/src/RingGrids/octahealpix.jl b/src/RingGrids/octahealpix.jl index f80e16c8f..d7280bf5a 100644 --- a/src/RingGrids/octahealpix.jl +++ b/src/RingGrids/octahealpix.jl @@ -105,10 +105,6 @@ nlon_octahealpix(nlat_half::Integer,j::Integer) = min(4j,8nlat_half-4j) OctaHEALPixGrid{T}(data::AbstractVector) where T = OctaHEALPixGrid{T}(data,nlat_half_octahealpix(length(data))) OctaHEALPixGrid(data::AbstractVector,n::Integer...) = OctaHEALPixGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:OctaHEALPixGrid}) = 3 # cubic -get_truncation(::Type{<:OctaHEALPixGrid},nlat_half::Integer) = nlat_half-1 -get_resolution(::Type{<:OctaHEALPixGrid},trunc::Integer) = trunc+1 - function get_colat(::Type{<:OctaHEALPixGrid},nlat_half::Integer) nlat_half == 0 && return Float64[] colat = zeros(nlat_octahealpix(nlat_half)) diff --git a/src/RingGrids/octahedral.jl b/src/RingGrids/octahedral.jl index 4b8c092a0..1c5d3d609 100644 --- a/src/RingGrids/octahedral.jl +++ b/src/RingGrids/octahedral.jl @@ -109,9 +109,6 @@ nlon_octahedral(j::Integer) = 16+4j OctahedralGaussianGrid{T}(data::AbstractVector) where T = OctahedralGaussianGrid{T}(data,nlat_half_octahedral(length(data),false)) OctahedralGaussianGrid(data::AbstractVector,n::Integer...) = OctahedralGaussianGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:OctahedralGaussianGrid}) = 3 # cubic -get_truncation(::Type{<:OctahedralGaussianGrid},nlat_half::Integer) = nlat_half-1 -get_resolution(::Type{<:OctahedralGaussianGrid},trunc::Integer) = trunc+1 nlat_odd(::Type{<:OctahedralGaussianGrid}) = false get_npoints(::Type{<:OctahedralGaussianGrid},nlat_half::Integer) = npoints_octahedral(nlat_half,false) get_colat(::Type{<:OctahedralGaussianGrid},nlat_half::Integer) = get_colat(FullGaussianGrid,nlat_half) @@ -144,9 +141,6 @@ OctahedralClenshawGrid{T}(data::AbstractVector) where T = OctahedralClenshawGrid nlat_half_octahedral(length(data),true)) OctahedralClenshawGrid(data::AbstractVector,n::Integer...) = OctahedralClenshawGrid{eltype(data)}(data,n...) -truncation_order(::Type{<:OctahedralClenshawGrid}) = 3 # cubic -get_truncation(::Type{<:OctahedralClenshawGrid},nlat_half::Integer) = nlat_half-1 -get_nlat_half(::Type{<:OctahedralClenshawGrid},trunc::Integer) = roundup_fft(trunc+1) nlat_odd(::Type{<:OctahedralClenshawGrid}) = true get_npoints(::Type{<:OctahedralClenshawGrid},nlat_half::Integer) = npoints_octahedral(nlat_half,true) get_colat(::Type{<:OctahedralClenshawGrid},nlat_half::Integer) = get_colat(FullClenshawGrid,nlat_half) diff --git a/src/RingGrids/utility_functions.jl b/src/RingGrids/utility_functions.jl index 8f5a3e657..a7fc379aa 100644 --- a/src/RingGrids/utility_functions.jl +++ b/src/RingGrids/utility_functions.jl @@ -1,38 +1,3 @@ -""" - true/false = is_power_2(i::Integer) - -Checks whether an integer `i` is a power of 2, i.e. i = 2^k, with k = 0,1,2,3,....""" -is_power_2(i::Integer) = i != 0 ? i & (i-1) == 0 : false -is_power_2_or_0(i::Integer) = i & (i-1) == 0 - -""" - m = roundup_fft(n::Int; - small_primes::Vector{Int}=[2,3,5]) - -Returns an integer `m >= n` with only small prime factors 2, 3, 5 (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. -""" -function roundup_fft(n::Integer;small_primes::Vector{T}=[2,3,5]) where {T<:Integer} - factors_not_in_small_primes = true # starting condition for while loop - n += isodd(n) ? 1 : 0 # start with an even n - while factors_not_in_small_primes - - factors = Primes.factor(n) # prime factorization - all_factors_small = true # starting condition - - for i in 1:length(factors) # loop over factors and check they are small - factor = factors.pe[i].first # extract factor from factors - all_factors_small &= factor in small_primes - end - - factors_not_in_small_primes = ~all_factors_small # all factors small will abort while loop - n += 2 # test for next larger even n - - end - return n-2 # subtract unnecessary last += 2 addition -end - """ true/false = isincreasing(v::Vector) diff --git a/src/SpeedyTransforms/SpeedyTransforms.jl b/src/SpeedyTransforms/SpeedyTransforms.jl index 50f317253..aad88fcd1 100644 --- a/src/SpeedyTransforms/SpeedyTransforms.jl +++ b/src/SpeedyTransforms/SpeedyTransforms.jl @@ -8,6 +8,7 @@ module SpeedyTransforms import FFTW import GenericFFT import LinearAlgebra + import Primes # SPEEDYWEATHER MODULES using ..LowerTriangularMatrices @@ -22,6 +23,9 @@ module SpeedyTransforms spectral, spectral! + # ALIASING + export get_nlat_half + # GRADIENTS export curl!, divergence!, @@ -36,6 +40,7 @@ module SpeedyTransforms spectral_truncation!, spectral_interpolation + include("aliasing.jl") include("spectral_transform.jl") include("spectral_gradients.jl") include("spectral_truncation.jl") diff --git a/src/SpeedyTransforms/aliasing.jl b/src/SpeedyTransforms/aliasing.jl new file mode 100644 index 000000000..42cc33d46 --- /dev/null +++ b/src/SpeedyTransforms/aliasing.jl @@ -0,0 +1,48 @@ +const DEFAULT_DEALIASING = 2.0 + +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)/(dealiasing+1)) +end + +""" + m = roundup_fft(n::Int; + small_primes::Vector{Int}=[2,3]) + +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. +""" +function roundup_fft(n::Integer;small_primes::Vector{T}=[2,3,5]) where {T<:Integer} + factors_not_in_small_primes = true # starting condition for while loop + n += isodd(n) ? 1 : 0 # start with an even n + while factors_not_in_small_primes + + factors = Primes.factor(n) # prime factorization + all_factors_small = true # starting condition + + for i in 1:length(factors) # loop over factors and check they are small + factor = factors.pe[i].first # extract factor from factors + all_factors_small &= factor in small_primes + end + + factors_not_in_small_primes = ~all_factors_small # all factors small will abort while loop + n += 2 # test for next larger even n + + end + return n-2 # subtract unnecessary last += 2 addition +end + +""" + true/false = is_power_2(i::Integer) + +Checks whether an integer `i` is a power of 2, i.e. i = 2^k, with k = 0,1,2,3,....""" +is_power_2(i::Integer) = i != 0 ? i & (i-1) == 0 : false +is_power_2_or_0(i::Integer) = i & (i-1) == 0 \ No newline at end of file diff --git a/src/SpeedyTransforms/spectral_transform.jl b/src/SpeedyTransforms/spectral_transform.jl index ada4feee0..679e50093 100644 --- a/src/SpeedyTransforms/spectral_transform.jl +++ b/src/SpeedyTransforms/spectral_transform.jl @@ -8,7 +8,6 @@ struct SpectralTransform{NF<:AbstractFloat} # GRID Grid::Type{<:AbstractGrid} # grid type used nlat_half::Int # resolution parameter of grid (# of latitudes on one hemisphere, Eq incl) - order::Int # 1,2,3 for linear, quadratic, cubic grid # SPECTRAL RESOLUTION lmax::Int # Maximum degree l=[0,lmax] of spherical harmonics @@ -76,16 +75,16 @@ function SpectralTransform( ::Type{NF}, # Number format NF Grid::Type{<:AbstractGrid}, # type of spatial grid used trunc::Int, # Spectral truncation recompute_legendre::Bool; # re or precompute legendre polynomials? - legendre_shortcut::Symbol=:linear # shorten Legendre loop over order m + legendre_shortcut::Symbol=:linear, # shorten Legendre loop over order m + dealiasing::Real=DEFAULT_DEALIASING ) where NF # SPECTRAL RESOLUTION lmax = trunc # Maximum degree l=[0,lmax] of spherical harmonics mmax = trunc # Maximum order m=[0,l] of spherical harmonics - order = RingGrids.truncation_order(Grid) # RESOLUTION PARAMETERS - nlat_half = get_nlat_half(Grid,trunc) # resolution parameter nlat_half, + nlat_half = get_nlat_half(trunc,dealiasing) # resolution parameter nlat_half, # number of latitude rings on one hemisphere incl equator nlat = get_nlat(Grid,nlat_half) # 2nlat_half but one less if grids have odd # of lat rings nlon_max = get_nlon_max(Grid,nlat_half) # number of longitudes around the equator @@ -201,7 +200,7 @@ function SpectralTransform( ::Type{NF}, # Number format NF eigenvalues⁻¹[1] = 0 # set the integration constant to 0 # conversion to NF happens here - SpectralTransform{NF}( Grid,nlat_half,order, + SpectralTransform{NF}( Grid,nlat_half, lmax,mmax,nfreq_max,m_truncs, nlon_max,nlons,nlat, colat,cos_colat,sin_colat,lon_offsets, @@ -226,7 +225,7 @@ function spectral_transform_for_full_grid(S::SpectralTransform{NF}) where NF FullGrid = full_grid(S.Grid) # corresponding full grid # unpack everything that does not have to be recomputed for the full grid - @unpack nresolution, order, lmax, mmax, nfreq_max, nlon_max, nlat, nlat_half = S + @unpack nresolution, lmax, mmax, nfreq_max, nlon_max, nlat, nlat_half = S @unpack colat, cos_colat, sin_colat, norm_sphere, recompute_legendre, Λ, Λs = S @unpack ϵlms, grad_x, grad_y1, grad_y2, grad_y_vordiv1, grad_y_vordiv2 = S @unpack vordiv_to_uv_x, vordiv_to_uv1, vordiv_to_uv2, eigenvalues, eigenvalues⁻¹ = S @@ -244,7 +243,7 @@ function spectral_transform_for_full_grid(S::SpectralTransform{NF}) where NF solid_angles = get_solid_angles(FullGrid,nlat_half) - return SpectralTransform{NF}( FullGrid,nresolution,order, + return SpectralTransform{NF}( FullGrid,nresolution, lmax,mmax,nfreq_max,m_truncs, nlon_max,nlons,nlat,nlat_half, colat,cos_colat,sin_colat,lon_offsets, @@ -344,7 +343,7 @@ function gridded!( map::AbstractGrid{NF}, # gridded output unscale_coslat::Bool=false # unscale with cos(lat) on the fly? ) where {NF<:AbstractFloat} # number format NF - @unpack nlat, nlons, nlat_half, nfreq_max, order = S + @unpack nlat, nlons, nlat_half, nfreq_max = S @unpack cos_colat, sin_colat, lon_offsets = S @unpack recompute_legendre, Λ, Λs, m_truncs = S @unpack brfft_plans = S @@ -353,7 +352,7 @@ function gridded!( map::AbstractGrid{NF}, # gridded output recompute_legendre || @boundscheck size(alms) == size(Λs[1]) || throw(BoundsError) lmax, mmax = size(alms) .- 1 # maximum degree l, order m of spherical harmonics - @boundscheck mmax+1 <= nfreq_max || throw(BoundsError) + @boundscheck maximum(m_truncs) <= nfreq_max || throw(BoundsError) @boundscheck nlat == length(cos_colat) || throw(BoundsError) @boundscheck typeof(map) <: S.Grid || throw(BoundsError) @boundscheck get_nlat_half(map) == S.nlat_half || throw(BoundsError) @@ -452,7 +451,7 @@ function spectral!( alms::LowerTriangularMatrix{Complex{NF}}, # output: spectr S::SpectralTransform{NF} ) where {NF<:AbstractFloat} - @unpack nlat, nlat_half, nlons, nfreq_max, order, cos_colat = S + @unpack nlat, nlat_half, nlons, nfreq_max, cos_colat = S @unpack recompute_legendre, Λ, Λs, solid_angles = S @unpack rfft_plans, lon_offsets, m_truncs = S @@ -460,7 +459,7 @@ function spectral!( alms::LowerTriangularMatrix{Complex{NF}}, # output: spectr recompute_legendre || @boundscheck size(alms) == size(Λs[1]) || throw(BoundsError) lmax, mmax = size(alms) .- 1 # maximum degree l, order m of spherical harmonics - @boundscheck mmax+1 <= nfreq_max || throw(BoundsError) + @boundscheck maximum(m_truncs) <= nfreq_max || throw(BoundsError) @boundscheck nlat == length(cos_colat) || throw(BoundsError) @boundscheck typeof(map) <: S.Grid || throw(BoundsError) @boundscheck get_nlat_half(map) == S.nlat_half || throw(BoundsError) @@ -600,7 +599,7 @@ function spectral( map::AbstractGrid{NF}; # gridded field ) where NF # number format NF Grid = typeof(map) - trunc = get_truncation(map) + trunc = get_truncation(map.nlat_half) S = SpectralTransform(NF,Grid,trunc,recompute_legendre) return spectral(map,S) end diff --git a/src/default_parameters.jl b/src/default_parameters.jl index 858a94760..6a8e686a7 100644 --- a/src/default_parameters.jl +++ b/src/default_parameters.jl @@ -26,8 +26,11 @@ Base.@kwdef struct Parameters{Model<:ModelSetup} <: AbstractParameters{Model} "grid in use" Grid::Type{<:AbstractGrid} = FullGaussianGrid + "dealiasing factor, 1=linear, 2=quadratic, 3=cubic grid" + dealiasing::Float64 = 2 - # EARTH'S PROPERTIES + + # PLANET'S PROPERTIES "planet" planet::Planet = Earth() diff --git a/src/dynamics/geometry.jl b/src/dynamics/geometry.jl index 27d04edf0..b64be1d9b 100644 --- a/src/dynamics/geometry.jl +++ b/src/dynamics/geometry.jl @@ -67,7 +67,7 @@ Generator function to create the Geometry struct from parameters in `P`. """ function Geometry(P::Parameters,Grid::Type{<:AbstractGrid}) - @unpack trunc, nlev = P # grid type, spectral truncation, # of vertical levels + @unpack trunc, dealiasing, nlev = P # grid type, spectral truncation, # of vertical levels @unpack radius, rotation, gravity = P.planet # radius of planet, angular frequency, gravity @unpack R_dry, cₚ = P # gas constant for dry air, heat capacity @unpack σ_tropopause = P # number of vertical levels used for stratosphere @@ -75,8 +75,9 @@ function Geometry(P::Parameters,Grid::Type{<:AbstractGrid}) @unpack ΔT_stratosphere = P # used for stratospheric temperature increase # RESOLUTION PARAMETERS - nlat_half = get_nlat_half(Grid,trunc) # resolution parameter nlat_half - # = number of lat rings on one hemisphere (Equator incl) + # resolution parameter nlat_half (= number of lat rings on one hemisphere (Equator incl) + # from spectral resolution and dealiasing parameter (e.g. quadratic grid for T31) + nlat_half = SpeedyTransforms.get_nlat_half(trunc,dealiasing) nlat = get_nlat(Grid,nlat_half) # 2nlat_half but -1 if grids have odd # of lat rings nlon_max = get_nlon_max(Grid,nlat_half) # number of longitudes around the equator nlon = nlon_max # same (used for compatibility) @@ -194,6 +195,6 @@ end Generator function for a SpectralTransform struct pulling in parameters from a Parameters struct.""" function SpeedyTransforms.SpectralTransform(P::Parameters) - @unpack NF, Grid, trunc, recompute_legendre, legendre_shortcut = P - return SpectralTransform(NF,Grid,trunc,recompute_legendre;legendre_shortcut) + @unpack NF, Grid, trunc, dealiasing, recompute_legendre, legendre_shortcut = P + return SpectralTransform(NF,Grid,trunc,recompute_legendre;legendre_shortcut,dealiasing) end \ No newline at end of file diff --git a/test/spectral_transform.jl b/test/spectral_transform.jl index 193f73ec5..a95f843fb 100644 --- a/test/spectral_transform.jl +++ b/test/spectral_transform.jl @@ -1,19 +1,40 @@ @testset "FullGaussianGrid: Test grid and spectral resolution match" begin - Ts = (31,42,85,170,341,682) # spectral resolutions - nlons = (96,128,256,512,1024) # number of longitudes - nlats = (48,64,128,256,512) # number of latitudes - for (T,nlon,nlat) in zip(Ts,nlons,nlats) - - nlat_half = SpeedyWeather.get_resolution(FullGaussianGrid,T) - @test (nlon,nlat) == (4nlat_half,2nlat_half) + + p = 5:10 + + @testset for dealiasing in [2,3] + # powers of two minus 1, T31, T63, T127, etc + Ts = [2^i - 1 for i in p] # spectral resolutions + nlons = [(dealiasing+1)*2^i for i in p] # number of longitudes + nlats = nlons/2 # number of latitudes + for (T,nlon,nlat) in zip(Ts,nlons,nlats) + nlat_half = SpeedyTransforms.get_nlat_half(T,dealiasing) + @test (nlon,nlat) == (4nlat_half,2nlat_half) + + trunc = SpeedyTransforms.get_truncation(nlat_half,dealiasing) + @test T == trunc + end + end - trunc = SpeedyWeather.get_truncation(FullGaussianGrid,nlat÷2) - @test T == trunc + # for these resolutions just test idempotence as the roundup_fft may + # give various other options than just the 3*2^n-matching + @testset for dealiasing in [2,3] + # T42,T85,T170,T341,T682,T1365 etc + Ts = [floor(Int,2^(i+2)/3) for i in p] + nlons = [dealiasing*2^(i+1) for i in p] + nlats = nlons/2 # number of latitudes + for (T,nlon,nlat) in zip(Ts,nlons,nlats) + nlat_half = SpeedyTransforms.get_nlat_half(T,dealiasing) + trunc = SpeedyTransforms.get_truncation(nlat_half,dealiasing) + nlat_half2 = SpeedyTransforms.get_nlat_half(trunc,dealiasing) + @test nlat_half == nlat_half2 + end end end # for the following testsets test some spectral truncations # but not too large ones as they take so long + spectral_resolutions = (31,63,127) spectral_resolutions_inexact = (127,255) @@ -132,8 +153,11 @@ end FullClenshawGrid, OctahedralGaussianGrid, OctahedralClenshawGrid) - - P = Parameters{ShallowWaterModel}(;NF,Grid,trunc) + + # clenshaw-curtis grids are only exact for cubic truncation + dealiasing = Grid in (FullGaussianGrid,OctahedralGaussianGrid) ? 2 : 3 + + P = Parameters{ShallowWaterModel}(;NF,Grid,trunc,dealiasing) S = SpectralTransform(P) G = Geometry(P) B = Boundaries(P,S,G) diff --git a/test/utility_functions.jl b/test/utility_functions.jl index 60a614f10..61d2765f6 100644 --- a/test/utility_functions.jl +++ b/test/utility_functions.jl @@ -27,13 +27,13 @@ end @testset "roundup nlon for FFT" begin for i in 1:10 - @test 2^i == RingGrids.roundup_fft(2^i) - @test 2^i*3 == RingGrids.roundup_fft(2^i*3) - @test 2^i*5 == RingGrids.roundup_fft(2^i*5) + @test 2^i == SpeedyTransforms.roundup_fft(2^i) + @test 2^i*3 == SpeedyTransforms.roundup_fft(2^i*3) + @test 2^i*5 == SpeedyTransforms.roundup_fft(2^i*5) end for n in 1:10 i = rand(2:1000) - @test i <= RingGrids.roundup_fft(i) + @test i <= SpeedyTransforms.roundup_fft(i) end end