Reduce computational complexity of tests

JoeyT1994 · Apr 10, 2024 · 098c8f6 · 098c8f6
1 parent 22b870c
commit 098c8f6
Show file tree

Hide file tree

Showing 3 changed files with 159 additions and 122 deletions.
diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl
@@ -74,7 +74,9 @@ function plus_shift_ttn(
   return ttn(ttn_op, s; algorithm="svd")
 end
 
-function minus_shift_ttn(s::IndsNetwork, bit_map; dimension=default_dimension(), boundary=default_boundary())
+function minus_shift_ttn(
+  s::IndsNetwork, bit_map; dimension=default_dimension(), boundary=default_boundary()
+)
   @assert is_tree(s)
   @assert base(bit_map) == 2
   ttn_op = OpSum()
@@ -114,14 +116,16 @@ function stencil(
   shifts::Vector{Float64},
   delta_power::Int64;
   dimension=default_dimension(),
-  left_boundary = default_boundary(),
-  right_boundary = default_boundary(),
+  left_boundary=default_boundary(),
+  right_boundary=default_boundary(),
   scale=true,
   truncate_kwargs...,
 )
   @assert length(shifts) == 3
-  plus_shift = first(shifts) * plus_shift_ttn(s, bit_map; dimension, boundary = right_boundary)
-  minus_shift = last(shifts) * minus_shift_ttn(s, bit_map; dimension, boundary = left_boundary)
+  plus_shift =
+    first(shifts) * plus_shift_ttn(s, bit_map; dimension, boundary=right_boundary)
+  minus_shift =
+    last(shifts) * minus_shift_ttn(s, bit_map; dimension, boundary=left_boundary)
   no_shift = shifts[2] * no_shift_ttn(s)
 
   stencil_op = plus_shift + minus_shift + no_shift

diff --git a/test/test_derivative_operators.jl b/test/test_derivative_operators.jl
@@ -9,7 +9,7 @@ using ITensorNumericalAnalysis: itensornetwork
 using Dictionaries: Dictionary
 
 @testset "test laplacian in 1D on MPS" begin
-  g = named_grid((16, 1))
+  g = named_grid((12, 1))
   L = nv(g)
   delta = (2.0)^(-Float64(L))
   bit_map = BitMap(g)
@@ -26,11 +26,26 @@ using Dictionaries: Dictionary
     ∂2x_ψ_fx_x = real(calculate_fx(∂2x_ψ_fx, x))
     @test ∂2x_ψ_fx_x ≈ -pi * pi * sin(pi * x) atol = 1e-3
   end
+
+  ∇sq = laplacian_operator(
+    s, bit_map; cutoff=1e-12, left_boundary="Periodic", right_boundary="Periodic"
+  )
+
+  ψ_fx = sin_itn(s, bit_map; k=2.0 * Float64(pi))
+  ∂2x_ψ_fx = operate(∇sq, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+
+  #Can include 0.0 now because the boundary conditions fix this
+  xs = [0.0, delta, 0.25, 0.625, 0.875]
+  for x in xs
+    ∂2x_ψ_fx_x = real(calculate_fx(∂2x_ψ_fx, x))
+    @test ∂2x_ψ_fx_x ≈ -4.0 * pi * pi * sin(2.0 * pi * x) atol = 1e-3
+  end
 end
 
 @testset "test derivative in 1D on tree" begin
-  g = named_comb_tree((4, 4))
+  g = named_comb_tree((4, 3))
   L = nv(g)
+  delta = 2.0^(-Float64(L))
   bit_map = BitMap(g)
   s = siteinds(g, bit_map)
 
@@ -39,15 +54,15 @@ end
   ψ_fx = sin_itn(s, bit_map; k=Float64(pi))
   ∂x_ψ_fx = operate(∂_∂x, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
 
-  xs = [0.025, 0.1, 0.25, 0.625, 0.875]
+  xs = [delta, 0.125, 0.25, 0.625, 0.875]
   for x in xs
     ∂x_ψ_fx_x = real(calculate_fx(∂x_ψ_fx, x))
-    @test ∂x_ψ_fx_x ≈ pi * cos(pi * x) atol = 1e-4
+    @test ∂x_ψ_fx_x ≈ pi * cos(pi * x) atol = 1e-3
   end
 end
 
 @testset "test multiplication_operator_in_1D" begin
-  g = named_comb_tree((4, 4))
+  g = named_comb_tree((4, 3))
   L = nv(g)
   bit_map = BitMap(g)
   s = siteinds(g, bit_map)
@@ -59,12 +74,12 @@ end
   xs = [0.025, 0.1, 0.25, 0.625, 0.875]
   for x in xs
     ψ_fxgx_x = real(calculate_fx(ψ_fxgx, x))
-    @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-4
+    @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3
   end
 end
 
 @testset "test multiplication_operator_in_2D" begin
-  L = 10
+  L = 8
   g = NamedGraph(SimpleGraph(uniform_tree(L)))
 
   bit_map = BitMap(g; map_dimension=2)
@@ -91,15 +106,15 @@ end
 end
 
 @testset "test differentiation_operator_on_3D_function" begin
-  L = 60
+  L = 45
   g = named_grid((L, 1))
 
   bit_map = BitMap(g; map_dimension=3)
   s = siteinds(g, bit_map)
 
   ψ_fx = poly_itn(s, bit_map, [0.0, -1.0, 1.0]; dimension=1)
   ψ_gy = sin_itn(s, bit_map; k=Float64(pi), dimension=2)
-  ψ_hz = sinh_itn(s, bit_map; k=Float64(pi), dimension=3)
+  ψ_hz = sin_itn(s, bit_map; k=Float64(pi), dimension=3)
   @assert dimension(ψ_fx) == dimension(ψ_gy) == dimension(ψ_hz) == 3
 
   ψ_fxgyhz = ψ_fx * ψ_gy * ψ_hz
@@ -115,10 +130,10 @@ end
     for y in ys
       for z in zs
         ψ_fxgyhz_xyz = real(calculate_fxyz(ψ_fxgyhz, [x, y, z]))
-        @test ψ_fxgyhz_xyz ≈ (x^2 - x) * sin(pi * y) * sinh(pi * z) atol = 1e-3
+        @test ψ_fxgyhz_xyz ≈ (x^2 - x) * sin(pi * y) * sin(pi * z) atol = 1e-3
 
         ∂_∂y_ψ_fxgyhz_xyz = real(calculate_fxyz(∂_∂y_ψ_fxgyhz, [x, y, z]))
-        @test ∂_∂y_ψ_fxgyhz_xyz ≈ pi * (x^2 - x) * cos(pi * y) * sinh(pi * z) atol = 1e-3
+        @test ∂_∂y_ψ_fxgyhz_xyz ≈ pi * (x^2 - x) * cos(pi * y) * sin(pi * z) atol = 1e-3
       end
     end
   end

diff --git a/test/test_shift_operators.jl b/test/test_shift_operators.jl
@@ -11,10 +11,9 @@ using Dictionaries: Dictionary
 using ITensorNumericalAnalysis: plus_shift_ttn, minus_shift_ttn
 
 @testset "test shift operators in 1D on MPS" begin
-
-  g = named_grid((8, 1))
+  g = named_grid((6, 1))
   L = nv(g)
-  delta = 2.0^(-1.0*L)
+  delta = 2.0^(-1.0 * L)
   bit_map = BitMap(g)
   s = siteinds(g, bit_map)
   xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta]
@@ -27,8 +26,12 @@ using ITensorNumericalAnalysis: plus_shift_ttn, minus_shift_ttn
   plus_shift_neumann = plus_shift_ttn(s, bit_map; boundary="Neumann")
   minus_shift_neumann = minus_shift_ttn(s, bit_map; boundary="Neumann")
 
-  ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-  ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fx_pshift_dirichlet = operate(
+    plus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fx_mshift_dirichlet = operate(
+    minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12)
+  )
   ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
   ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
   ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
@@ -43,120 +46,135 @@ using ITensorNumericalAnalysis: plus_shift_ttn, minus_shift_ttn
     elseif x == 1.0 - delta
       @test calculate_fx(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8
       @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
-      @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
+      @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol =
+        1e-8
     end
 
     if x - delta >= 0.0
-        fx_xminus = calculate_fx(ψ_fx, x - delta)
-        @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8
-        @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8
-        @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8
+      fx_xminus = calculate_fx(ψ_fx, x - delta)
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8
     elseif x == 0.0
-        @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8
-        @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
-        @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
+      @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8
+      @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
+      @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
     end
   end
 end
 
-
 @testset "test shift operators in 1D on Tree" begin
+  g = named_comb_tree((2, 3))
+  L = nv(g)
+  delta = 2.0^(-1.0 * L)
+  bit_map = BitMap(g)
+  s = siteinds(g, bit_map)
+  xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta]
+  ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25])
+
+  plus_shift_dirichlet = plus_shift_ttn(s, bit_map; boundary="Dirichlet")
+  minus_shift_dirichlet = minus_shift_ttn(s, bit_map; boundary="Dirichlet")
+  plus_shift_pbc = plus_shift_ttn(s, bit_map; boundary="Periodic")
+  minus_shift_pbc = minus_shift_ttn(s, bit_map; boundary="Periodic")
+  plus_shift_neumann = plus_shift_ttn(s, bit_map; boundary="Neumann")
+  minus_shift_neumann = minus_shift_ttn(s, bit_map; boundary="Neumann")
+
+  ψ_fx_pshift_dirichlet = operate(
+    plus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fx_mshift_dirichlet = operate(
+    minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
+
+  for x in xs
+    if x + delta < 1
+      fx_xplus = calculate_fx(ψ_fx, x + delta)
+      @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_dirichlet, x) atol = 1e-8
+      @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_pbc, x) atol = 1e-8
+      @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_neumann, x) atol = 1e-8
+    elseif x == 1.0 - delta
+      @test calculate_fx(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8
+      @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
+      @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol =
+        1e-8
+    end
 
-    g = named_comb_tree((3,3))
-    L = nv(g)
-    delta = 2.0^(-1.0*L)
-    bit_map = BitMap(g)
-    s = siteinds(g, bit_map)
-    xs = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta]
-    ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25])
-
-    plus_shift_dirichlet = plus_shift_ttn(s, bit_map; boundary="Dirichlet")
-    minus_shift_dirichlet = minus_shift_ttn(s, bit_map; boundary="Dirichlet")
-    plus_shift_pbc = plus_shift_ttn(s, bit_map; boundary="Periodic")
-    minus_shift_pbc = minus_shift_ttn(s, bit_map; boundary="Periodic")
-    plus_shift_neumann = plus_shift_ttn(s, bit_map; boundary="Neumann")
-    minus_shift_neumann = minus_shift_ttn(s, bit_map; boundary="Neumann")
-
-    ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; truncate_kwargs=(; cutoff=1e-12))
-
-    for x in xs
-      if x + delta < 1
-        fx_xplus = calculate_fx(ψ_fx, x + delta)
-        @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_dirichlet, x) atol = 1e-8
-        @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_pbc, x) atol = 1e-8
-        @test fx_xplus ≈ calculate_fx(ψ_fx_pshift_neumann, x) atol = 1e-8
-      elseif x == 1.0 - delta
-        @test calculate_fx(ψ_fx_pshift_dirichlet, x) ≈ 0.0 atol = 1e-8
-        @test calculate_fx(ψ_fx_pshift_pbc, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
-        @test calculate_fx(ψ_fx_pshift_neumann, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
-      end
-
-      if x - delta >= 0.0
-          fx_xminus = calculate_fx(ψ_fx, x - delta)
-          @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8
-          @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8
-          @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8
-        elseif x == 0.0
-          @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8
-          @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
-          @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
-        end
+    if x - delta >= 0.0
+      fx_xminus = calculate_fx(ψ_fx, x - delta)
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8
+      @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8
+    elseif x == 0.0
+      @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8
+      @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
+      @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
     end
+  end
 end
 
 @testset "test shift operators in 2D on Tree" begin
+  g = named_comb_tree((3, 3))
+  bit_map = BitMap(g; map_dimension=2)
+  L = length(vertices(bit_map, 2))
+  delta = 2.0^(-1.0 * L)
+  s = siteinds(g, bit_map)
+  x = 0.5
+  ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta]
+  ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25]; dimension=1)
+  ψ_fy = cos_itn(s, bit_map; dimension=2)
+  ψ_fxy = ψ_fx + ψ_fx
+
+  plus_shift_dirichlet = plus_shift_ttn(s, bit_map; dimension=2, boundary="Dirichlet")
+  minus_shift_dirichlet = minus_shift_ttn(s, bit_map; dimension=2, boundary="Dirichlet")
+  plus_shift_pbc = plus_shift_ttn(s, bit_map; dimension=2, boundary="Periodic")
+  minus_shift_pbc = minus_shift_ttn(s, bit_map; dimension=2, boundary="Periodic")
+  plus_shift_neumann = plus_shift_ttn(s, bit_map; dimension=2, boundary="Neumann")
+  minus_shift_neumann = minus_shift_ttn(s, bit_map; dimension=2, boundary="Neumann")
 
-    g = named_comb_tree((3,3))
-    bit_map = BitMap(g; map_dimension = 2)
-    L = length(vertices(bit_map, 2))
-    delta = 2.0^(-1.0*L)
-    s = siteinds(g, bit_map)
-    x = 0.5
-    ys = [0.0, delta, 0.25, 0.5, 0.625, 0.875, 1.0 - delta]
-    ψ_fx = poly_itn(s, bit_map, [1.0, 0.5, 0.25]; dimension = 1) 
-    ψ_fy = cos_itn(s, bit_map; dimension = 2)
-    ψ_fxy = ψ_fx + ψ_fx
-
-    plus_shift_dirichlet = plus_shift_ttn(s, bit_map; dimension = 2, boundary="Dirichlet")
-    minus_shift_dirichlet = minus_shift_ttn(s, bit_map; dimension = 2, boundary="Dirichlet")
-    plus_shift_pbc = plus_shift_ttn(s, bit_map; dimension = 2, boundary="Periodic")
-    minus_shift_pbc = minus_shift_ttn(s, bit_map; dimension = 2, boundary="Periodic")
-    plus_shift_neumann = plus_shift_ttn(s, bit_map; dimension = 2, boundary="Neumann")
-    minus_shift_neumann = minus_shift_ttn(s, bit_map; dimension = 2, boundary="Neumann")
-
-    ψ_fxy_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fxy_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fxy_pshift_pbc = operate(plus_shift_pbc, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fxy_mshift_pbc = operate(minus_shift_pbc, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fxy_pshift_neumann = operate(plus_shift_neumann, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-    ψ_fxy_mshift_neumann = operate(minus_shift_neumann, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
-
-    for y in ys
-      if y + delta < 1
-        fxy_xyplus = calculate_fxy(ψ_fxy, [x, y + delta])
-        @test fxy_xyplus ≈ calculate_fxy(ψ_fxy_pshift_dirichlet, [x,y]) atol = 1e-8
-        @test fxy_xyplus ≈ calculate_fxy(ψ_fxy_pshift_pbc, [x,y]) atol = 1e-8
-        @test fxy_xyplus ≈ calculate_fxy(ψ_fxy_pshift_neumann, [x,y]) atol = 1e-8
-      elseif y == 1.0 - delta
-        @test calculate_fxy(ψ_fxy_pshift_dirichlet, [x,y]) ≈ 0.0 atol = 1e-8
-        @test calculate_fx(ψ_fxy_pshift_pbc, [x,y]) ≈ calculate_fxy(ψ_fxy, [x, 0.0]) atol = 1e-8
-        @test calculate_fx(ψ_fxy_pshift_neumann, [x,y]) ≈ calculate_fxy(ψ_fxy, [x, 1.0 - delta]) atol = 1e-8
-      end
-
-    #   if x - delta >= 0.0
-    #       fx_xminus = calculate_fx(ψ_fx, x - delta)
-    #       @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_dirichlet, x) atol = 1e-8
-    #       @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_pbc, x) atol = 1e-8
-    #       @test fx_xminus ≈ calculate_fx(ψ_fx_mshift_neumann, x) atol = 1e-8
-    #   elseif x == 0.0
-    #       @test calculate_fx(ψ_fx_mshift_dirichlet, x) ≈ 0.0 atol = 1e-8
-    #       @test calculate_fx(ψ_fx_mshift_pbc, x) ≈ calculate_fx(ψ_fx, 1.0 - delta) atol = 1e-8
-    #       @test calculate_fx(ψ_fx_mshift_neumann, x) ≈ calculate_fx(ψ_fx, 0.0) atol = 1e-8
-    #   end
+  ψ_fxy_pshift_dirichlet = operate(
+    plus_shift_dirichlet, ψ_fxy; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fxy_mshift_dirichlet = operate(
+    minus_shift_dirichlet, ψ_fxy; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fxy_pshift_pbc = operate(plus_shift_pbc, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fxy_mshift_pbc = operate(minus_shift_pbc, ψ_fxy; truncate_kwargs=(; cutoff=1e-12))
+  ψ_fxy_pshift_neumann = operate(
+    plus_shift_neumann, ψ_fxy; truncate_kwargs=(; cutoff=1e-12)
+  )
+  ψ_fxy_mshift_neumann = operate(
+    minus_shift_neumann, ψ_fxy; truncate_kwargs=(; cutoff=1e-12)
+  )
+
+  for y in ys
+    if y + delta < 1
+      fxy_xyplus = calculate_fxyz(ψ_fxy, [x, y + delta])
+      @test fxy_xyplus ≈ calculate_fxyz(ψ_fxy_pshift_dirichlet, [x, y]) atol = 1e-8
+      @test fxy_xyplus ≈ calculate_fxyz(ψ_fxy_pshift_pbc, [x, y]) atol = 1e-8
+      @test fxy_xyplus ≈ calculate_fxyz(ψ_fxy_pshift_neumann, [x, y]) atol = 1e-8
+    elseif y == 1.0 - delta
+      @test calculate_fxyz(ψ_fxy_pshift_dirichlet, [x, y]) ≈ 0.0 atol = 1e-8
+      @test calculate_fxyz(ψ_fxy_pshift_pbc, [x, y]) ≈ calculate_fxyz(ψ_fxy, [x, 0.0]) atol =
+        1e-8
+      @test calculate_fxyz(ψ_fxy_pshift_neumann, [x, y]) ≈
+        calculate_fxyz(ψ_fxy, [x, 1.0 - delta]) atol = 1e-8
+    end
+
+    if y - delta >= 0.0
+      fxy_xyminus = calculate_fxyz(ψ_fxy, [x, y - delta])
+      @test fxy_xyminus ≈ calculate_fxyz(ψ_fxy_mshift_dirichlet, [x, y]) atol = 1e-8
+      @test fxy_xyminus ≈ calculate_fxyz(ψ_fxy_mshift_pbc, [x, y]) atol = 1e-8
+      @test fxy_xyminus ≈ calculate_fxyz(ψ_fxy_mshift_neumann, [x, y]) atol = 1e-8
+    elseif y == 0.0
+      @test calculate_fxyz(ψ_fxy_mshift_dirichlet, [x, y]) ≈ 0.0 atol = 1e-8
+      @test calculate_fxyz(ψ_fxy_mshift_pbc, [x, y]) ≈
+        calculate_fxyz(ψ_fxy, [x, 1.0 - delta]) atol = 1e-8
+      @test calculate_fxyz(ψ_fxy_mshift_neumann, [x, y]) ≈ calculate_fxyz(ψ_fxy, [x, 0.0]) atol =
+        1e-8
     end
-end
+  end
+end