From 56ad3a9e60cee785b105e8e96a80ad6653774e32 Mon Sep 17 00:00:00 2001
From: meudnaes <eliasudnes@hotmail.com>
Date: Thu, 13 Jun 2024 09:06:15 +0200
Subject: [PATCH] [DEV] testing interpolation and derivative routines against
 analytic

---
 test/stagger_operators.jl | 103 ++++++++++++++++++++++----------------
 1 file changed, 60 insertions(+), 43 deletions(-)

diff --git a/test/stagger_operators.jl b/test/stagger_operators.jl
index 37f32bf..e697548 100644
--- a/test/stagger_operators.jl
+++ b/test/stagger_operators.jl
@@ -1,48 +1,20 @@
+# Test 5th order interpolation with 5th order polynomial
+function p5(x::Real)
+    Float32(-0.001x^5 + 0.03x^4 - 0.3x^3 + 1.5x^2 - 2x + 5)
+end
+
+function dpdx(x::Real)
+    Float32(-0.005x^4 + 0.12x^3 - 0.9x^2 + 3x - 2)
+end
+
+# Second order polynomial to test boundaries
+function p2(x::Real)
+    Float32(-0.5x^2 + 2x - 3)
+end
+
 @testset "interpolation" begin
-    ####
-    # Not all combinations of slices/columns are going to pass these tests.
-    # If there is a rapid change in the variable, the 5th order Bifrost interpolation
-    # and cubic spline interpolation are going to give different results.  
-    ####
-    @testset "z direction" begin
-        # Load a column
-        bz_stagger = get_var(xp,xp.snaps,"bz",
-            destagger=false,slicex=[10],slicey=[10])
-
-        # Destagger by the 5th order Bifrost Interpolation
-        bz = zup(bz_stagger)
-        
-        # Use Intepolations cubic spline to interpolate
-        x = 1:xp.mesh.mz
-        x_new = x .+ 0.5
-        itp = cubic_spline_interpolation(x,bz_stagger[1,1,:],extrapolation_bc=Line())
-        bz_itp = itp(x_new)
-
-        @test isapprox(bz_itp,bz[1,1,:],atol=1e-3)
-    end
-    
-    @testset "y direction" begin
-        # load a strange narrow 3x48x3 cube
-        by_stagger = get_var(xp,xp.snaps,"by",
-            destagger=true,slicex=1:3,slicez=1:3)
-
-        # Destagger by the 5th order Bifrost Interpolation
-        by = yup(by_stagger,true)
-        
-        # Use Intepolations cubic spline to interpolate
-        x = 1:xp.mesh.my
-        x_new = x .+ 0.5
-
-        xx = 1:3
-        yy = 1:3
-
-        itp = cubic_spline_interpolation((xx,x,yy),by_stagger,extrapolation_bc=Line())
-        by_itp = itp[xx,x_new,yy]
-
-        @test isapprox(by_itp,by,atol=1e-3)
-    end
 
-    @testset "x direction" begin   
+    @testset "Extrapolation" begin   
         x_stagger = Float32.(1:10) .- 0.5f0
         x_stagger = reshape(x_stagger,(10,1,1))
 
@@ -53,4 +25,49 @@
         @test x[:,1,1] == Float32.(0:9)
     end
 
+    @testset "5th order exact interpolation" begin   
+        x = Float32.(1:10) .- 0.5f0
+        y = p5.(x)
+
+        # Reshape to 3D array, check against inner parts
+
+        # x direction
+        p_x = reshape(y,(10,1,1))
+        @test xup(p_x,false)[3:end-3,1,1] ≈ p5.(x .+ 0.5)[3:end-3]
+        @test xdn(p_x,false)[4:end-2,1,1] ≈ p5.(x .- 0.5)[4:end-2]
+
+        # y direction
+        p_y = reshape(y,(1,10,1))
+        @test yup(p_y,false)[1,3:end-3,1] ≈ p5.(x .+ 0.5)[3:end-3]
+        @test ydn(p_y,false)[1,4:end-2,1] ≈ p5.(x .- 0.5)[4:end-2]
+
+        # z direction
+        p_z = reshape(y,(1,1,10))
+        @test zup(p_z,false)[1,1,3:end-3] ≈ p5.(x .+ 0.5)[3:end-3]
+        @test zdn(p_z,false)[1,1,4:end-2] ≈ p5.(x .- 0.5)[4:end-2]
+    end
+
+    @testset "6th order exact derivative" begin   
+        x = Float32.(1:10) .- 0.5f0
+        dx = ones(Float32,10)
+        y = p5.(x)
+
+        # Reshape to 3D array, check against inner parts
+
+        # x direction
+        p_x = reshape(y,(10,1,1))
+        @test dxup(p_x,dx,false)[3:end-3,1,1] ≈ dpdx.(x .+ 0.5)[3:end-3]
+        @test dxdn(p_x,dx,false)[4:end-2,1,1] ≈ dpdx.(x .- 0.5)[4:end-2]
+
+        # y direction
+        p_y = reshape(y,(1,10,1))
+        @test dyup(p_y,dx,false)[1,3:end-3,1] ≈ dpdx.(x .+ 0.5)[3:end-3]
+        @test dydn(p_y,dx,false)[1,4:end-2,1] ≈ dpdx.(x .- 0.5)[4:end-2]
+
+        # z direction
+        p_z = reshape(y,(1,1,10))
+        @test dzup(p_z,dx,false)[1,1,3:end-3] ≈ dpdx.(x .+ 0.5)[3:end-3]
+        @test dzdn(p_z,dx,false)[1,1,4:end-2] ≈ dpdx.(x .- 0.5)[4:end-2]
+    end
+
 end
\ No newline at end of file