orography smoothing was disabled through old size call

SpeedyWeather · Jul 4, 2024 · acbe0c5 · acbe0c5
1 parent 575529c
commit acbe0c5
Show file tree

Hide file tree

Showing 6 changed files with 21 additions and 17 deletions.
diff --git a/ext/SpeedyWeatherJLArraysExt.jl b/ext/SpeedyWeatherJLArraysExt.jl
@@ -4,7 +4,7 @@ using SpeedyWeather, JLArrays
 
 # for RingGrids and LowerTriangularMatrices:
 # every Array needs this method to strip away the parameters
-SpeedyWeather.RingGrids.nonparametric_type(::Type{<:JLArray}) = JLArray
-SpeedyWeather.LowerTriangularMatrices.nonparametric_type(::Type{<:JLArray}) = JLArray
+RingGrids.nonparametric_type(::Type{<:JLArray}) = JLArray
+LowerTriangularMatrices.nonparametric_type(::Type{<:JLArray}) = JLArray
 
 end # module
diff --git a/src/SpeedyTransforms/legendrepolarray.jl b/src/SpeedyTransforms/legendrepolarray.jl
@@ -6,8 +6,8 @@ This enables easier use with AssociatedLegendrePolynomials.jl which otherwise co
 "matrix-style" (l, m) indexing of `LowerTriangularArray`. This type however doesn't support any
 other operations than indexing and is purerly intended for internal purposes. 
 """
-struct AssociatedLegendrePolArray{T,N,M,V} <: AbstractArray{T,N}
-    data::LowerTriangularArray{T,M,V}
+struct AssociatedLegendrePolArray{T, N, M, V} <: AbstractArray{T, N}
+    data::LowerTriangularArray{T, M, V}
 end 
 
 """2-dimensional `AssociatedLegendrePolArray` of type `T`` with its non-zero entries unravelled into a `Vector{T}`"""

diff --git a/src/SpeedyTransforms/spectral_truncation.jl b/src/SpeedyTransforms/spectral_truncation.jl
@@ -176,12 +176,12 @@ end
 
 Smooth the spectral field `A` following A *= (1-(1-c)*∇²ⁿ) with power n of a normalised Laplacian
 so that the highest degree lmax is dampened by multiplication with c. Anti-diffusion for c>1."""
-function spectral_smoothing!(   A::LowerTriangularMatrix,
+function spectral_smoothing!(   L::LowerTriangularMatrix,
                                 c::Real;
-                                power::Real=1,          # power of Laplacian used for smoothing
+                                power::Real=1,          # power of Laplacian used for smoothing
                                 truncation::Int=-1)     # smoothing wrt wavenumber (0 = largest)
 
-    lmax, mmax = size(A; as=Matrix)
+    lmax, mmax = size(L; as=Matrix)
     # normalize by largest eigenvalue by default, or wrt to given truncation
     eigenvalue_norm = truncation == -1 ? -mmax*(mmax+1) : -truncation*(truncation+1)
 
@@ -191,7 +191,7 @@ function spectral_smoothing!(   A::LowerTriangularMatrix,
             eigenvalue_normalised = -l*(l-1)/eigenvalue_norm
             # for eigenvalue_norm < largest eigenvalue the factor becomes negative
             # set to zero in that case
-            A[lm] *= max(1 - (1-c)*eigenvalue_normalised^power, 0)
+            L[lm] *= max(1 - (1-c)*eigenvalue_normalised^power, 0)
             lm += 1
         end
     end

diff --git a/src/dynamics/orography.jl b/src/dynamics/orography.jl
@@ -163,13 +163,12 @@ function initialize!(   orog::EarthOrography,
     orography .*= scale                         # scale orography (default 1)
     spectral!(geopot_surf, orography, S)        # no *gravity yet
 
-    if orog.smoothing                       # smooth orography in spectral space?
-        # translate smoothing_fraction to trunc to truncate beyond
-        trunc = (size(geopot_surf, 1) - 2)
-        truncation = round(Int, trunc * (1-orog.smoothing_fraction))
-        SpeedyTransforms.spectral_smoothing!(geopot_surf, orog.smoothing_strength,
-                                                            power=orog.smoothing_power;
-                                                            truncation)
+    if orog.smoothing                           # smooth orography in spectral space?
+        trunc = (size(geopot_surf, 1, as=Matrix) - 2)   # get trunc=lmax from size of geopot_surf
+        truncation = round(Int, trunc * (1-orog.smoothing_fraction))    # degree of harmonics to be truncated
+        c = orog.smoothing_strength
+        power = orog.smoothing_power
+        SpeedyTransforms.spectral_smoothing!(geopot_surf, c; power, truncation)
     end
 
     gridded!(orography, geopot_surf, S)     # to grid-point space

diff --git a/test/lower_triangular_matrix.jl b/test/lower_triangular_matrix.jl
@@ -213,6 +213,9 @@ end
     end
 end
 
+# only needed when the extension isn't loaded
+# LowerTriangularMatrices.nonparametric_type(::Type{<:JLArray}) = JLArray
+
 @testset "Zeros, ones, rand, and randn constructors" begin
     for f in (ones, zeros, rand, randn)
         s = (5, 5)

diff --git a/test/spectral_transform.jl b/test/spectral_transform.jl
@@ -136,8 +136,10 @@ end
                         map = gridded(alms, S)
                         alms2 = spectral(map, S)
 
+                        tol = 1e-3
+
                         for lm in SpeedyWeather.eachharmonic(alms, alms2)
-                            @test alms[lm] ≈ alms2[lm] atol=1e-3 rtol=1e-3
+                            @test alms[lm] ≈ alms2[lm] atol=tol rtol=tol
                         end
                     end
                 end
@@ -173,7 +175,7 @@ end
                 oro_grid2 = gridded(oro_spec1, S)
                 oro_spec2 = spectral(oro_grid2, S)
 
-                tol = 1e-1
+                tol = 1e-3
 
                 for lm in SpeedyWeather.eachharmonic(oro_spec1, oro_spec2)
                     @test oro_spec1[lm] ≈ oro_spec2[lm] atol=tol rtol=tol