From 098c8f6721d2298a19a2c935c709d94abccce3d3 Mon Sep 17 00:00:00 2001 From: Joey Date: Wed, 10 Apr 2024 16:44:13 -0400 Subject: [PATCH] Reduce computational complexity of tests --- src/elementary_operators.jl | 14 +- test/test_derivative_operators.jl | 37 +++-- test/test_shift_operators.jl | 230 ++++++++++++++++-------------- 3 files changed, 159 insertions(+), 122 deletions(-) diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 7b7dc12..5193b67 100644 --- 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 index aa8e73b..23e8edb 100644 --- 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,7 +106,7 @@ 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) @@ -99,7 +114,7 @@ end ψ_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 index 84d4ebd..6854ed3 100644 --- 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 \ No newline at end of file + end +end