Dealiasing parameter for quad,cub truncation

src/RingGrids/RingGrids.jl

@@ -2,7 +2,6 @@ module RingGrids
 
     import Parameters: @unpack
     import Statistics: mean
-    import Primes
     import FastGaussQuadrature
 
     # GRIDS

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)

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))

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[]
 

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))

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)

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)
 

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")

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

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,15 +451,15 @@ 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
 
     recompute_legendre && @boundscheck size(alms) == size(Λ) || throw(BoundsError)
     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

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()

src/dynamics/geometry.jl

@@ -67,16 +67,17 @@ 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
     @unpack temp_ref, temp_top, lapse_rate = P      # for reference atmosphere
     @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