Merge pull request #208 from SpeedyWeather/mk/temperature

test ∇×∇=0 and ∇⋅∇=∇²
SpeedyWeather · Dec 20, 2022 · bf4cdb4 · bf4cdb4
2 parents 73a44e8 + 9bba706
commit bf4cdb4
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 11 deletions.
diff --git a/src/default_parameters.jl b/src/default_parameters.jl
@@ -112,7 +112,6 @@ The default values of the keywords define the default model setup.
     orography::Union{Bool,Real} = true  # switch on/off orography or scale it by a factor
     orography_path::String = boundary_path
     orography_file::String = "orography_F512.nc"
-
 
     # INITIAL CONDITIONS
     seed::Int = abs(rand(Int))          # random seed for the global random number generator

diff --git a/test/spectral_gradients.jl b/test/spectral_gradients.jl
@@ -285,25 +285,45 @@ end
     end
 end
 
-@testset "∇ no errors" begin
+@testset "∇×∇=0 and ∇⋅∇=∇²" begin
     for NF in (Float32,Float64)
-        NF = Float64
         p,d,m = initialize_speedy(NF)
 
         a = rand(LowerTriangularMatrix{Complex{NF}},33,32)
+        SpeedyWeather.spectral_truncation!(a,31)
+        a[:,1] .= real.(a[:,1])
+
+        println(rand(NF))
+
         dadx = zero(a)
         dady = zero(a)
         SpeedyWeather.∇!(dadx,dady,a,m.spectral_transform)
+
+        dadx_grid = gridded(dadx,m.spectral_transform)
+        dady_grid = gridded(dady,m.spectral_transform)
 
-        # coslat unscaling missing in the following
-        # ∇²a_1 = zero(a)
-        # ∇²a_2 = zero(a)
+        SpeedyWeather.scale_coslat⁻²!(dadx_grid,m.geometry)
+        SpeedyWeather.scale_coslat⁻²!(dady_grid,m.geometry)
+
+        SpeedyWeather.spectral!(dadx,dadx_grid,m.spectral_transform)
+        SpeedyWeather.spectral!(dady,dady_grid,m.spectral_transform)
 
-        # SpeedyWeather.divergence!(∇²a_1,dadx,dady,m.spectral_transform)
-        # SpeedyWeather.∇²!(∇²a_2,a,m.spectral_transform)
+        # CURL(GRAD(A)) = 0
+        ∇x∇a = zero(a)
+        SpeedyWeather.curl!(∇x∇a,dadx,dady,m.spectral_transform)
 
-        # for lm in SpeedyWeather.eachharmonic(∇²a_1,∇²a_2)
-        #     @test ∇²a_1[lm] ≈ ∇²a_2[lm]
-        # end
+        for lm in SpeedyWeather.eachharmonic(∇x∇a)
+            @test ∇x∇a[lm] ≈ 0 atol=3*sqrt(eps(NF))
+        end
+
+        # DIV(GRAD(A)) = LAPLACE(A)
+        div_∇a = zero(a)
+        SpeedyWeather.divergence!(div_∇a,dadx,dady,m.spectral_transform)
+        ∇²a = zero(a)
+        SpeedyWeather.∇²!(∇²a,a,m.spectral_transform)
+
+        for lm in SpeedyWeather.eachharmonic(div_∇a,∇²a)
+            @test div_∇a[lm] ≈ ∇²a[lm] atol=sqrt(eps(NF)) rtol=sqrt(eps(NF))
+        end
     end
 end