diff --git a/src/elementary_operators.jl b/src/elementary_operators.jl index 48d9f8a..b9461ae 100644 --- a/src/elementary_operators.jl +++ b/src/elementary_operators.jl @@ -45,8 +45,12 @@ function ITensors.op(::OpName"Dup", ::SiteType"Digit", s::Index) return ITensor(o, s, s') end -function plus_shift_ttn( - s::IndsNetwork, bit_map; dimension=default_dimension(), boundary=default_boundary(), n::Int64 = 0 +function plus_shift_opsum( + s::IndsNetwork, + bit_map; + dimension=default_dimension(), + boundary=default_boundary(), + n::Int64=0, ) @assert is_tree(s) @assert base(bit_map) == 2 @@ -56,7 +60,7 @@ function plus_shift_ttn( string_site = [("D+", vertex(bit_map, dimension, L - n))] add!(ttn_op, 1.0, "D+", vertex(bit_map, dimension, L - n)) - for i in (L-n):-1:2 + for i in (L - n):-1:2 pop!(string_site) push!(string_site, ("D-", vertex(bit_map, dimension, i))) push!(string_site, ("D+", vertex(bit_map, dimension, i - 1))) @@ -64,24 +68,34 @@ function plus_shift_ttn( end if boundary == "Neumann" - #TODO: Not convinced this is right.... - for i in 0:n - string_site = [j <= (L -i) ? ("Dup", vertex(bit_map, dimension, j)) : ("Ddn", vertex(bit_map, dimension, j)) for j in 1:L] - add!(ttn_op, 1.0, (string_site...)...) - end + string_site = [ + if j <= (L - n) + ("Dup", vertex(bit_map, dimension, j)) + else + ("I", vertex(bit_map, dimension, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) elseif boundary == "Periodic" - #TODO: Not convinced this is right.... - for i in 0:n - string_site = [j <= (L -i) ? ("D-", vertex(bit_map, dimension, j)) : ("D+", vertex(bit_map, dimension, j)) for j in 1:L] - add!(ttn_op, 1.0, (string_site...)...) - end + string_site = [ + if j <= (L - n) + ("D-", vertex(bit_map, dimension, j)) + else + ("I", vertex(bit_map, dimension, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) end - return ttn(ttn_op, s; algorithm="svd") + return ttn_op end -function minus_shift_ttn( - s::IndsNetwork, bit_map; dimension=default_dimension(), boundary=default_boundary(), n::Int64 = 0 +function minus_shift_opsum( + s::IndsNetwork, + bit_map; + dimension=default_dimension(), + boundary=default_boundary(), + n::Int64=0, ) @assert is_tree(s) @assert base(bit_map) == 2 @@ -91,7 +105,7 @@ function minus_shift_ttn( string_site = [("D-", vertex(bit_map, dimension, L - n))] add!(ttn_op, 1.0, "D-", vertex(bit_map, dimension, L - n)) - for i in (L-n):-1:2 + for i in (L - n):-1:2 pop!(string_site) push!(string_site, ("D+", vertex(bit_map, dimension, i))) push!(string_site, ("D-", vertex(bit_map, dimension, i - 1))) @@ -99,25 +113,33 @@ function minus_shift_ttn( end if boundary == "Neumann" - for i in 0:n - string_site = [j <= (L -i) ? ("Ddn", vertex(bit_map, dimension, j)) : ("Dup", vertex(bit_map, dimension, j)) for j in 1:L] - add!(ttn_op, 1.0, (string_site...)...) - end + string_site = [ + if j <= (L - n) + ("Ddn", vertex(bit_map, dimension, j)) + else + ("I", vertex(bit_map, dimension, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) elseif boundary == "Periodic" - for i in 0:n - string_site = [j <= (L -i) ? ("D+", vertex(bit_map, dimension, j)) : ("D-", vertex(bit_map, dimension, j)) for j in 1:L] - add!(ttn_op, 1.0, (string_site...)...) - end + string_site = [ + if j <= (L - n) + ("D+", vertex(bit_map, dimension, j)) + else + ("I", vertex(bit_map, dimension, j)) + end for j in 1:L + ] + add!(ttn_op, 1.0, (string_site...)...) end - return ttn(ttn_op, s; algorithm="svd") + return ttn_op end -function no_shift_ttn(s::IndsNetwork) +function no_shift_opsum(s::IndsNetwork) ttn_op = OpSum() string_site_full = [("I", v) for v in vertices(s)] add!(ttn_op, 1.0, (string_site_full...)...) - return ttn(ttn_op, s; algorithm="svd") + return ttn_op end function stencil( @@ -129,25 +151,28 @@ function stencil( left_boundary=default_boundary(), right_boundary=default_boundary(), scale=true, - truncate_kwargs..., + truncate_op=true, + kwargs..., ) @assert length(shifts) == 5 - stencil_op = shifts[3] * no_shift_ttn(s) - for i in [1,2] - n = i == 1 ? 1 : 0 - if !iszero(i) - stencil_op += shifts[i] * plus_shift_ttn(s, bit_map; dimension, boundary=right_boundary, n) + stencil_opsum = shifts[3] * no_shift_opsum(s) + for i in [1, 2] + n = i == 1 ? 1 : 0 + if !iszero(shifts[i]) + stencil_opsum += + shifts[i] * plus_shift_opsum(s, bit_map; dimension, boundary=right_boundary, n) end end - for i in [4,5] - n = i == 5 ? 1 : 0 - if !iszero(i) - stencil_op += shifts[i] * minus_shift_ttn(s, bit_map; dimension, boundary=left_boundary, n) + for i in [4, 5] + n = i == 5 ? 1 : 0 + if !iszero(shifts[i]) + stencil_opsum += + shifts[i] * minus_shift_opsum(s, bit_map; dimension, boundary=left_boundary, n) end end - stencil_op = truncate(stencil_op; truncate_kwargs...) + stencil_op = ttn(stencil_opsum, s; algorithm="svd", kwargs...) if scale for v in vertices(bit_map, dimension) @@ -167,7 +192,7 @@ function second_derivative_operator(s::IndsNetwork, bit_map; kwargs...) end function third_derivative_operator(s::IndsNetwork, bit_map; kwargs...) - return stencil(s, bit_map, [-0.5, 1.0, 0.0, -1.0, 0.5], 3; kwargs...) + return stencil(s, bit_map, [0.5, -1.0, 0.0, 1.0, -0.5], 3; kwargs...) end function fourth_derivative_operator(s::IndsNetwork, bit_map; kwargs...) @@ -227,12 +252,9 @@ end Base.:*(fs::ITensorNetworkFunction...) = multiply(fs...) -function operate( - operator::TreeTensorNetwork, ψ::ITensorNetworkFunction; truncate_kwargs=(;), kwargs... -) +function operate(operator::TreeTensorNetwork, ψ::ITensorNetworkFunction; kwargs...) ψ_tn = ttn(itensornetwork(ψ)) ψO_tn = noprime(contract(operator, ψ_tn; init=prime(copy(ψ_tn)), kwargs...)) - ψO_tn = truncate(ψO_tn; truncate_kwargs...) return ITensorNetworkFunction(ITensorNetwork(ψO_tn), bit_map(ψ)) end diff --git a/test/test_derivative_operators.jl b/test/test_derivative_operators.jl index d5a7f97..727d9e6 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 differentiation in 1D on MPS" begin - g = named_grid((12, 1)) + g = named_grid((9, 1)) L = nv(g) delta = (2.0)^(-Float64(L)) bit_map = BitMap(g) @@ -22,114 +22,118 @@ using Dictionaries: Dictionary f4 = fourth_derivative_operator(s, bit_map; cutoff=1e-12, left_boundary, right_boundary) ψ_fx = sin_itn(s, bit_map; k=2.0 * Float64(pi)) - ψ_f1x = operate(f1, ψ_fx; truncate_kwargs=(; cutoff=1e-12)) - ψ_f2x = operate(f2, ψ_fx; truncate_kwargs=(; cutoff=1e-12)) - ψ_f3x = operate(f3, ψ_fx; truncate_kwargs=(; cutoff=1e-12)) - ψ_f4x = operate(f4, ψ_fx; truncate_kwargs=(; cutoff=1e-12)) - xs = [0.0, delta, 0.25, 0.625, 0.875, 1.0 - delta] + ψ_f1x = operate(f1, ψ_fx; cutoff=1e-8) + ψ_f2x = operate(f2, ψ_fx; cutoff=1e-8) + ψ_f3x = operate(f3, ψ_fx; cutoff=1e-8) + ψ_f4x = operate(f4, ψ_fx; cutoff=1e-8) + + xs = [0.0, 0.25, 0.625, 0.875, 1.0 - delta] + for x in xs + @test 1.0 + calculate_fx(ψ_fx, x) ≈ 1.0 + sin(2.0 * pi * x) rtol = 1e-3 + @test 1.0 + calculate_fx(ψ_f1x, x) ≈ 1.0 + 2.0 * pi * cos(2.0 * pi * x) rtol = 1e-3 + @test 1.0 + calculate_fx(ψ_f2x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^2 * sin(2.0 * pi * x) rtol = + 1e-3 + @test 1.0 + calculate_fx(ψ_f3x, x) ≈ 1.0 + -1.0 * (2.0 * pi)^3 * cos(2.0 * pi * x) rtol = + 1e-3 + @test 1.0 + calculate_fx(ψ_f4x, x) ≈ 1.0 + 1.0 * (2.0 * pi)^4 * sin(2.0 * pi * x) rtol = + 1e-3 + end +end + +@testset "test derivative in 1D on tree" begin + g = named_comb_tree((4, 3)) + L = nv(g) + delta = 2.0^(-Float64(L)) + bit_map = BitMap(g) + s = siteinds(g, bit_map) + + ∂_∂x = first_derivative_operator(s, bit_map; cutoff=1e-10) + + ψ_fx = sin_itn(s, bit_map; k=Float64(pi)) + ∂x_ψ_fx = operate(∂_∂x, ψ_fx; cutoff=1e-12) + + 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-3 + end +end + +@testset "test multiplication_operator_in_1D" begin + g = named_comb_tree((4, 3)) + L = nv(g) + bit_map = BitMap(g) + s = siteinds(g, bit_map) + + ψ_gx = sin_itn(s, bit_map; k=0.5 * Float64(pi)) + ψ_fx = cos_itn(s, bit_map; k=0.25 * Float64(pi)) + + ψ_fxgx = ψ_gx * ψ_fx + xs = [0.025, 0.1, 0.25, 0.625, 0.875] for x in xs - @test calculate_fx(ψ_fx, x) ≈ sin(2.0 * pi * x) atol = 1e-3 - @test calculate_fx(ψ_f1x, x) ≈ 2.0 * pi * cos(2.0 * pi * x) atol = 1e-3 - @test calculate_fx(ψ_f2x, x) ≈ -(2.0 * pi)^(2) * sin(2.0 * pi * x) atol = 1e-3 - # @test calculate_fx(ψ_f3x, x) ≈ -(2.0 * pi)^(3) * cos(2.0 * pi * x) atol = 1e-3 - @test calculate_fx(ψ_f4x, x) ≈ (2.0 * pi)^(4) * sin(2.0 * pi * x) atol = 1e-3 + ψ_fxgx_x = real(calculate_fx(ψ_fxgx, x)) + @test ψ_fxgx_x ≈ sin(0.5 * pi * x) * cos(0.25 * pi * x) atol = 1e-3 end end -# @testset "test derivative in 1D on tree" begin -# g = named_comb_tree((4, 3)) -# L = nv(g) -# delta = 2.0^(-Float64(L)) -# bit_map = BitMap(g) -# s = siteinds(g, bit_map) - -# ∂_∂x = first_derivative_operator(s, bit_map; cutoff=1e-10) - -# ψ_fx = sin_itn(s, bit_map; k=Float64(pi)) -# ∂x_ψ_fx = operate(∂_∂x, ψ_fx; truncate_kwargs=(; cutoff=1e-12)) - -# 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-3 -# end -# end - -# @testset "test multiplication_operator_in_1D" begin -# g = named_comb_tree((4, 3)) -# L = nv(g) -# bit_map = BitMap(g) -# s = siteinds(g, bit_map) - -# ψ_gx = sin_itn(s, bit_map; k=0.5 * Float64(pi)) -# ψ_fx = cos_itn(s, bit_map; k=0.25 * Float64(pi)) - -# ψ_fxgx = ψ_gx * ψ_fx -# 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-3 -# end -# end - -# @testset "test multiplication_operator_in_2D" begin -# L = 8 -# g = NamedGraph(SimpleGraph(uniform_tree(L))) - -# bit_map = BitMap(g; map_dimension=2) -# s = siteinds(g, bit_map) - -# ψ_fx = cos_itn(s, bit_map; k=0.25 * Float64(pi), dimension=1) -# ψ_gy = sin_itn(s, bit_map; k=0.5 * Float64(pi), dimension=2) -# @assert dimension(ψ_fx) == dimension(ψ_gy) == 2 - -# ψ_fxgy = ψ_fx * ψ_gy - -# xs = [0.125, 0.25, 0.625, 0.875] -# ys = [0.125, 0.25, 0.625, 0.875] -# for x in xs -# for y in ys -# ψ_fx_x = real(calculate_fxyz(ψ_fx, [x, y])) -# ψ_gy_y = real(calculate_fxyz(ψ_gy, [x, y])) -# @test ψ_fx_x ≈ cos(0.25 * pi * x) -# @test ψ_gy_y ≈ sin(0.5 * pi * y) -# ψ_fxgy_xy = real(calculate_fxyz(ψ_fxgy, [x, y])) -# @test ψ_fxgy_xy ≈ cos(0.25 * pi * x) * sin(0.5 * pi * y) atol = 1e-3 -# end -# end -# end - -# @testset "test differentiation_operator_on_3D_function" begin -# 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 = sin_itn(s, bit_map; k=Float64(pi), dimension=3) -# @assert dimension(ψ_fx) == dimension(ψ_gy) == dimension(ψ_hz) == 3 - -# ψ_fxgyhz = ψ_fx * ψ_gy * ψ_hz - -# ∂_∂y = first_derivative_operator(s, bit_map; dimension=2, cutoff=1e-10) - -# ∂_∂y_ψ_fxgyhz = operate([∂_∂y], ψ_fxgyhz; truncate_kwargs=(; cutoff=1e-10)) - -# xs = [0.125, 0.25, 0.675] -# ys = [0.125, 0.25, 0.675] -# zs = [0.125, 0.25, 0.675] -# for x in xs -# 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) * 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) * sin(pi * z) atol = 1e-3 -# end -# end -# end -# end +@testset "test multiplication_operator_in_2D" begin + L = 8 + g = NamedGraph(SimpleGraph(uniform_tree(L))) + + bit_map = BitMap(g; map_dimension=2) + s = siteinds(g, bit_map) + + ψ_fx = cos_itn(s, bit_map; k=0.25 * Float64(pi), dimension=1) + ψ_gy = sin_itn(s, bit_map; k=0.5 * Float64(pi), dimension=2) + @assert dimension(ψ_fx) == dimension(ψ_gy) == 2 + + ψ_fxgy = ψ_fx * ψ_gy + + xs = [0.125, 0.25, 0.625, 0.875] + ys = [0.125, 0.25, 0.625, 0.875] + for x in xs + for y in ys + ψ_fx_x = real(calculate_fxyz(ψ_fx, [x, y])) + ψ_gy_y = real(calculate_fxyz(ψ_gy, [x, y])) + @test ψ_fx_x ≈ cos(0.25 * pi * x) + @test ψ_gy_y ≈ sin(0.5 * pi * y) + ψ_fxgy_xy = real(calculate_fxyz(ψ_fxgy, [x, y])) + @test ψ_fxgy_xy ≈ cos(0.25 * pi * x) * sin(0.5 * pi * y) atol = 1e-3 + end + end +end + +@testset "test differentiation_operator_on_3D_function" begin + 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 = sin_itn(s, bit_map; k=Float64(pi), dimension=3) + @assert dimension(ψ_fx) == dimension(ψ_gy) == dimension(ψ_hz) == 3 + + ψ_fxgyhz = ψ_fx * ψ_gy * ψ_hz + + ∂_∂y = first_derivative_operator(s, bit_map; dimension=2, cutoff=1e-10) + + ∂_∂y_ψ_fxgyhz = operate([∂_∂y], ψ_fxgyhz; cutoff=1e-10) + + xs = [0.125, 0.25, 0.675] + ys = [0.125, 0.25, 0.675] + zs = [0.125, 0.25, 0.675] + for x in xs + 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) * 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) * sin(pi * z) atol = 1e-3 + end + end + end +end diff --git a/test/test_shift_operators.jl b/test/test_shift_operators.jl index 6854ed3..07aa9fa 100644 --- a/test/test_shift_operators.jl +++ b/test/test_shift_operators.jl @@ -8,7 +8,7 @@ using NamedGraphs: named_grid, named_comb_tree, NamedGraph, nv, vertices using ITensorNumericalAnalysis: itensornetwork using Dictionaries: Dictionary -using ITensorNumericalAnalysis: plus_shift_ttn, minus_shift_ttn +using ITensorNumericalAnalysis: stencil @testset "test shift operators in 1D on MPS" begin g = named_grid((6, 1)) @@ -19,23 +19,30 @@ using ITensorNumericalAnalysis: plus_shift_ttn, minus_shift_ttn 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) + plus_shift_dirichlet = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Dirichlet" + ) + minus_shift_dirichlet = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Dirichlet" + ) + plus_shift_pbc = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Periodic" + ) + minus_shift_pbc = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Periodic" ) - ψ_fx_mshift_dirichlet = operate( - minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12) + plus_shift_neumann = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Neumann" ) - ψ_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)) + minus_shift_neumann = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Neumann" + ) + ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; cutoff=1e-12) for x in xs if x + delta < 1 @@ -72,23 +79,31 @@ end 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) + plus_shift_dirichlet = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Dirichlet" + ) + minus_shift_dirichlet = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Dirichlet" + ) + plus_shift_pbc = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Periodic" ) - ψ_fx_mshift_dirichlet = operate( - minus_shift_dirichlet, ψ_fx; truncate_kwargs=(; cutoff=1e-12) + minus_shift_pbc = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Periodic" + ) + plus_shift_neumann = stencil( + s, bit_map, [0.0, 1.0, 0.0, 0.0, 0.0], 0; scale=false, right_boundary="Neumann" + ) + minus_shift_neumann = stencil( + s, bit_map, [0.0, 0.0, 0.0, 1.0, 0.0], 0; scale=false, left_boundary="Neumann" ) - ψ_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)) + + ψ_fx_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_pbc = operate(plus_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_pbc = operate(minus_shift_pbc, ψ_fx; cutoff=1e-12) + ψ_fx_pshift_neumann = operate(plus_shift_neumann, ψ_fx; cutoff=1e-12) + ψ_fx_mshift_neumann = operate(minus_shift_neumann, ψ_fx; cutoff=1e-12) for x in xs if x + delta < 1 @@ -128,27 +143,67 @@ end ψ_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) + plus_shift_dirichlet = stencil( + s, + bit_map, + [0.0, 1.0, 0.0, 0.0, 0.0], + 0; + scale=false, + dimension=2, + right_boundary="Dirichlet", ) - ψ_fxy_mshift_dirichlet = operate( - minus_shift_dirichlet, ψ_fxy; truncate_kwargs=(; cutoff=1e-12) + minus_shift_dirichlet = stencil( + s, + bit_map, + [0.0, 0.0, 0.0, 1.0, 0.0], + 0; + scale=false, + dimension=2, + left_boundary="Dirichlet", ) - ψ_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) + plus_shift_pbc = stencil( + s, + bit_map, + [0.0, 1.0, 0.0, 0.0, 0.0], + 0; + scale=false, + dimension=2, + right_boundary="Periodic", ) - ψ_fxy_mshift_neumann = operate( - minus_shift_neumann, ψ_fxy; truncate_kwargs=(; cutoff=1e-12) + minus_shift_pbc = stencil( + s, + bit_map, + [0.0, 0.0, 0.0, 1.0, 0.0], + 0; + scale=false, + dimension=2, + left_boundary="Periodic", ) + plus_shift_neumann = stencil( + s, + bit_map, + [0.0, 1.0, 0.0, 0.0, 0.0], + 0; + scale=false, + dimension=2, + right_boundary="Neumann", + ) + minus_shift_neumann = stencil( + s, + bit_map, + [0.0, 0.0, 0.0, 1.0, 0.0], + 0; + scale=false, + dimension=2, + left_boundary="Neumann", + ) + + ψ_fxy_pshift_dirichlet = operate(plus_shift_dirichlet, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_dirichlet = operate(minus_shift_dirichlet, ψ_fxy; cutoff=1e-12) + ψ_fxy_pshift_pbc = operate(plus_shift_pbc, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_pbc = operate(minus_shift_pbc, ψ_fxy; cutoff=1e-12) + ψ_fxy_pshift_neumann = operate(plus_shift_neumann, ψ_fxy; cutoff=1e-12) + ψ_fxy_mshift_neumann = operate(minus_shift_neumann, ψ_fxy; cutoff=1e-12) for y in ys if y + delta < 1