From b8eb51127d7107d369684ba1b669a1c41f51c1b9 Mon Sep 17 00:00:00 2001 From: Milan Date: Sun, 24 Sep 2023 23:02:33 -0400 Subject: [PATCH 1/3] Coriolis on grid function --- src/dynamics/constants.jl | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/src/dynamics/constants.jl b/src/dynamics/constants.jl index d63c42191..2142fe39a 100644 --- a/src/dynamics/constants.jl +++ b/src/dynamics/constants.jl @@ -117,3 +117,39 @@ function DynamicsConstants( spectral_grid::SpectralGrid, Δp_geopot_half, Δp_geopot_full, temp_ref_profile) end + +# default angular frequency of Earth's rotation [1/s] +const DEFAULT_ROTATION = 7.29e-5 + +""" +$(TYPEDSIGNATURES) +Return the Coriolis parameter `f` on the grid `Grid` of resolution `nlat_half` +on a planet of `ratation` [1/s]. Default rotation of Earth.""" +function coriolis( + ::Type{Grid}, + nlat_half::Integer; + rotation=DEFAULT_ROTATION +) where {Grid<:AbstractGrid} + + f = zeros(Grid,nlat_half) + lat = get_lat(Grid,nlat_half) + + for (j,ring) in enumerate(eachring(f)) + fⱼ = 2rotation*sin(lat[j]) + for ij in ring + f[ij] = fⱼ + end + end + return f +end + +""" +$(TYPEDSIGNATURES) +Return the Coriolis parameter `f` on a grid like `grid` +on a planet of `ratation` [1/s]. Default rotation of Earth.""" +function coriolis( + grid::Grid; + rotation=DEFAULT_ROTATION +) where {Grid<:AbstractGrid} + return coriolis(Grid,grid.nlat_half;rotation) +end \ No newline at end of file From 1acd2d4e25caa313a518e6e14e74110599fdee0d Mon Sep 17 00:00:00 2001 From: Milan Date: Sun, 24 Sep 2023 23:03:02 -0400 Subject: [PATCH 2/3] Convenience functions for div,curl,laplace --- src/SpeedyTransforms/spectral_gradients.jl | 88 +++++++++++++++++++++- 1 file changed, 86 insertions(+), 2 deletions(-) diff --git a/src/SpeedyTransforms/spectral_gradients.jl b/src/SpeedyTransforms/spectral_gradients.jl index defb6a292..a4c202496 100644 --- a/src/SpeedyTransforms/spectral_gradients.jl +++ b/src/SpeedyTransforms/spectral_gradients.jl @@ -96,10 +96,12 @@ This function requires both `u,v` to be transforms of fields that are scaled wit RingGrids.scale_coslat⁻¹!(u_grid) RingGrids.scale_coslat⁻¹!(v_grid) - u = spectral(u_grid) - v = spectral(v_grid) + u = spectral(u_grid,one_more_degree=true) + v = spectral(v_grid,one_more_degree=true) div = divergence(u,v) div_grid = gridded(div) + +Both `div` and `div_grid` are scaled with the radius. """ function divergence(u::LowerTriangularMatrix, v::LowerTriangularMatrix) @@ -112,6 +114,43 @@ function divergence(u::LowerTriangularMatrix, return div end +function _div_or_curl( kernel!, + u::Grid, + v::Grid) where {Grid<:AbstractGrid} + + @assert size(u) == size(v) "Size $(size(u)) and $(size(v)) incompatible." + + u_grid = copy(u) + v_grid = copy(v) + + RingGrids.scale_coslat⁻¹!(u_grid) + RingGrids.scale_coslat⁻¹!(v_grid) + + S = SpectralTransform(u_grid,one_more_degree=true) + us = spectral(u_grid,S) + vs = spectral(v_grid,S) + + div_or_vor = similar(us) + kernel!(div_or_vor,us,vs,S,add=false,flipsign=false) + return div_or_vor +end + +""" +$(TYPEDSIGNATURES) +Divergence (∇⋅) of two vector components `u,v` on a grid. +Applies 1/coslat scaling, transforms to spectral space and returns +the spectral divergence, which is scaled with the radius +of the sphere. Divide by radius for unscaling.""" +divergence(u::Grid,v::Grid) where {Grid<:AbstractGrid} = _div_or_curl(divergence!,u,v) + +""" +$(TYPEDSIGNATURES) +Curl (∇×) of two vector components `u,v` on a grid. +Applies 1/coslat scaling, transforms to spectral space and returns +the spectral curl, which is scaled with the radius +of the sphere. Divide by radius for unscaling.""" +curl(u::Grid,v::Grid) where {Grid<:AbstractGrid} = _div_or_curl(curl!,u,v) + """ $(TYPEDSIGNATURES) Curl (∇×) of two vector components `u,v` of size (n+1)xn, the last row @@ -337,6 +376,46 @@ function ∇²!( ∇²alms::LowerTriangularMatrix{Complex{NF}}, # Output: (inv return ∇²alms end +""" +$(TYPEDSIGNATURES) +Laplace operator ∇² applied to input `alms`, using precomputed +eigenvalues from `S`. The Laplace operator acts on the unit +sphere and therefore omits the 1/radius^2 scaling""" +function ∇²(alms::LowerTriangularMatrix, # Input: spectral coefficients + S::SpectralTransform) # precomputed eigenvalues + + ∇²alms = similar(alms) + ∇²!(∇²alms,alms,S,add=false,flipsign=false,inverse=false) + return ∇²alms +end + +""" +$(TYPEDSIGNATURES) +Returns the Laplace operator ∇² applied to input `alms`. +The Laplace operator acts on the unit +sphere and therefore omits the 1/radius^2 scaling""" +∇²(alms::LowerTriangularMatrix) = ∇²(alms,SpectralTransform(alms)) + +""" +$(TYPEDSIGNATURES) +InverseLaplace operator ∇⁻² applied to input `alms`, using precomputed +eigenvalues from `S`. The Laplace operator acts on the unit +sphere and therefore omits the radius^2 scaling""" +function ∇⁻²( ∇²alms::LowerTriangularMatrix, # Input: spectral coefficients + S::SpectralTransform) # precomputed eigenvalues + + alms = similar(∇²alms) + ∇⁻²!(alms,∇²alms,S,add=false,flipsign=false) + return alms +end + +""" +$(TYPEDSIGNATURES) +Returns the inverse Laplace operator ∇⁻² applied to input `alms`. +The Laplace operator acts on the unit +sphere and therefore omits the radius^2 scaling""" +∇⁻²(∇²alms::LowerTriangularMatrix) = ∇⁻²(∇²alms,SpectralTransform(∇²alms)) + """ ∇⁻²!( ∇⁻²alms::LowerTriangularMatrix, alms::LowerTriangularMatrix, @@ -356,6 +435,11 @@ function ∇⁻²!( ∇⁻²alms::LowerTriangularMatrix{Complex{NF}},# Output: return ∇²!(∇⁻²alms,alms,S;add,flipsign,inverse) end +""" +$(TYPEDSIGNATURES) +Applies the gradient operator ∇ applied to input `p` and stores the result in +`dpdx` (zonal derivative) and `dpdy` (meridional derivative). The gradient operator acts +on the unit sphere and therefore omits the radius scaling""" function ∇!(dpdx::LowerTriangularMatrix{Complex{NF}}, # Output: zonal gradient dpdy::LowerTriangularMatrix{Complex{NF}}, # Output: meridional gradient p::LowerTriangularMatrix{Complex{NF}}, # Input: spectral coefficients From 0d3a05a8e6cf297f5dec8eeb8bef4a390f501ffa Mon Sep 17 00:00:00 2001 From: Milan Date: Sun, 24 Sep 2023 23:03:24 -0400 Subject: [PATCH 3/3] SpectralTransform generator with one_more_degree --- src/SpeedyTransforms/spectral_transform.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/SpeedyTransforms/spectral_transform.jl b/src/SpeedyTransforms/spectral_transform.jl index 8e20bdef6..da8655ec4 100644 --- a/src/SpeedyTransforms/spectral_transform.jl +++ b/src/SpeedyTransforms/spectral_transform.jl @@ -233,11 +233,12 @@ Generator function for a `SpectralTransform` struct based on the size and grid t gridded field `map`. Recomputes the Legendre polynomials by default.""" function SpectralTransform( map::AbstractGrid{NF}; # gridded field recompute_legendre::Bool=true, # saves memory + one_more_degree::Bool=false, ) where NF # number format NF Grid = typeof(map) trunc = get_truncation(map) - return SpectralTransform(NF,Grid,trunc,trunc;recompute_legendre) + return SpectralTransform(NF,Grid,trunc+one_more_degree,trunc;recompute_legendre) end """