diff --git a/examples/03_solvers.jl b/examples/03_solvers.jl index ae4425d..7ece30a 100644 --- a/examples/03_solvers.jl +++ b/examples/03_solvers.jl @@ -17,13 +17,13 @@ function ode_solver( end time_span = (current_time, current_time + time_step) - u₀, ITensor_from_vec = to_vec(ψ₀) + u₀, to_itensor = to_vec(ψ₀) f(ψ::ITensor, p, t) = H(t)(ψ) - f(u::Vector, p, t) = to_vec(f(ITensor_from_vec(u), p, t))[1] + f(u::Vector, p, t) = to_vec(f(to_itensor(u), p, t))[1] prob = ODEProblem(f, u₀, time_span) sol = solve(prob, solver_alg; kwargs...) uₜ = sol.u[end] - return ITensor_from_vec(uₜ), nothing + return to_itensor(uₜ), nothing end function krylov_solver( diff --git a/examples/03_tdvp_time_dependent.jl b/examples/03_tdvp_time_dependent.jl index 3cd5b96..124caf3 100644 --- a/examples/03_tdvp_time_dependent.jl +++ b/examples/03_tdvp_time_dependent.jl @@ -3,149 +3,152 @@ using ITensorTDVP: tdvp using LinearAlgebra: norm using Random: Random -Random.seed!(1234) - -# Define the time-independent model include("03_models.jl") - -# Define the solvers needed for TDVP include("03_solvers.jl") -# Time dependent Hamiltonian is: -# H(t) = H₁(t) + H₂(t) + … -# = f₁(t) H₁(0) + f₂(t) H₂(0) + … -# = cos(ω₁t) H₁(0) + cos(ω₂t) H₂(0) + … - -# Number of sites -n = 6 - -# How much information to output from TDVP -# Set to 2 to get information about each bond/site -# evolution, and 3 to get information about the -# solver. -outputlevel = 1 - -# Frequency of time dependent terms -ω₁ = 0.1 -ω₂ = 0.2 - -# Nearest and next-nearest neighbor -# Heisenberg couplings. -J₁ = 1.0 -J₂ = 1.0 - -time_step = 0.1 -time_stop = 1.0 - -# nsite-update TDVP -nsite = 2 - -# Starting state bond/link dimension. -# A product state starting state can -# cause issues for TDVP without -# subspace expansion. -start_linkdim = 4 - -# TDVP truncation parameters -maxdim = 100 -cutoff = 1e-8 - -tol = 1e-15 - -# ODE solver parameters -ode_alg = Tsit5() -ode_kwargs = (; reltol=tol, abstol=tol) - -# Krylov solver parameters -krylov_kwargs = (; tol=tol, eager=true) - -@show n -@show ω₁, ω₂ -@show J₁, J₂ -@show maxdim, cutoff, nsite -@show start_linkdim -@show time_step, time_stop -@show ode_alg -@show ode_kwargs -@show krylov_kwargs - -ω⃗ = [ω₁, ω₂] -f⃗ = [t -> cos(ω * t) for ω in ω⃗] - -# H₀ = H(0) = H₁(0) + H₂(0) + … -ℋ₁₀ = heisenberg(n; J=J₁, J2=0.0) -ℋ₂₀ = heisenberg(n; J=0.0, J2=J₂) -ℋ⃗₀ = [ℋ₁₀, ℋ₂₀] - -s = siteinds("S=1/2", n) - -H⃗₀ = [MPO(ℋ₀, s) for ℋ₀ in ℋ⃗₀] - -# Initial state, ψ₀ = ψ(0) -# Initialize as complex since that is what OrdinaryDiffEq.jl/DifferentialEquations.jl -# expects. -ψ₀ = complex.(randomMPS(s, j -> isodd(j) ? "↑" : "↓"; linkdims=start_linkdim)) - -@show norm(ψ₀) - -println() -println("#"^100) -println("Running TDVP with ODE solver") -println("#"^100) -println() - -function ode_solver(H⃗₀, time_step, ψ₀; kwargs...) - return ode_solver( - -im * TimeDependentSum(f⃗, H⃗₀), - time_step, - ψ₀; - solver_alg=ode_alg, - ode_kwargs..., - kwargs..., +function main() + Random.seed!(1234) + + # Time dependent Hamiltonian is: + # H(t) = H₁(t) + H₂(t) + … + # = f₁(t) H₁(0) + f₂(t) H₂(0) + … + # = cos(ω₁t) H₁(0) + cos(ω₂t) H₂(0) + … + + # Number of sites + n = 6 + + # How much information to output from TDVP + # Set to 2 to get information about each bond/site + # evolution, and 3 to get information about the + # solver. + outputlevel = 1 + + # Frequency of time dependent terms + ω₁ = 0.1 + ω₂ = 0.2 + + # Nearest and next-nearest neighbor + # Heisenberg couplings. + J₁ = 1.0 + J₂ = 1.0 + + time_step = 0.1 + time_stop = 1.0 + + # nsite-update TDVP + nsite = 2 + + # Starting state bond/link dimension. + # A product state starting state can + # cause issues for TDVP without + # subspace expansion. + start_linkdim = 4 + + # TDVP truncation parameters + maxdim = 100 + cutoff = 1e-8 + + tol = 1e-15 + + # ODE solver parameters + ode_alg = Tsit5() + ode_kwargs = (; reltol=tol, abstol=tol) + + # Krylov solver parameters + krylov_kwargs = (; tol=tol, eager=true) + + @show n + @show ω₁, ω₂ + @show J₁, J₂ + @show maxdim, cutoff, nsite + @show start_linkdim + @show time_step, time_stop + @show ode_alg + @show ode_kwargs + @show krylov_kwargs + + ω⃗ = [ω₁, ω₂] + f⃗ = [t -> cos(ω * t) for ω in ω⃗] + + # H₀ = H(0) = H₁(0) + H₂(0) + … + ℋ₁₀ = heisenberg(n; J=J₁, J2=0.0) + ℋ₂₀ = heisenberg(n; J=0.0, J2=J₂) + ℋ⃗₀ = [ℋ₁₀, ℋ₂₀] + + s = siteinds("S=1/2", n) + + H⃗₀ = [MPO(ℋ₀, s) for ℋ₀ in ℋ⃗₀] + + # Initial state, ψ₀ = ψ(0) + # Initialize as complex since that is what OrdinaryDiffEq.jl/DifferentialEquations.jl + # expects. + ψ₀ = complex.(randomMPS(s, j -> isodd(j) ? "↑" : "↓"; linkdims=start_linkdim)) + + @show norm(ψ₀) + + println() + println("#"^100) + println("Running TDVP with ODE solver") + println("#"^100) + println() + + function ode_solver(H⃗₀, time_step, ψ₀; kwargs...) + return ode_solver( + -im * TimeDependentSum(f⃗, H⃗₀), + time_step, + ψ₀; + solver_alg=ode_alg, + ode_kwargs..., + kwargs..., + ) + end + + ψₜ_ode = tdvp( + ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel ) -end -ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite, outputlevel) + println() + println("Finished running TDVP with ODE solver") + println() -println() -println("Finished running TDVP with ODE solver") -println() + println() + println("#"^100) + println("Running TDVP with Krylov solver") + println("#"^100) + println() -println() -println("#"^100) -println("Running TDVP with Krylov solver") -println("#"^100) -println() + function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...) + return krylov_solver( + -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs... + ) + end -function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...) - return krylov_solver( - -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs... - ) -end + ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel) -ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite, outputlevel) + println() + println("Finished running TDVP with Krylov solver") + println() -println() -println("Finished running TDVP with Krylov solver") -println() + println() + println("#"^100) + println("Running full state evolution with ODE solver") + println("#"^100) + println() -println() -println("#"^100) -println("Running full state evolution with ODE solver") -println("#"^100) -println() + @disable_warn_order begin + ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel) + end -@disable_warn_order begin - ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀); outputlevel) -end + println() + println("Finished full state evolution with ODE solver") + println() -println() -println("Finished full state evolution with ODE solver") -println() + @show norm(ψₜ_ode) + @show norm(ψₜ_krylov) + @show norm(ψₜ_full) -@show norm(ψₜ_ode) -@show norm(ψₜ_krylov) -@show norm(ψₜ_full) + @show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full)) + @show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full)) +end -@show 1 - abs(inner(prod(ψₜ_ode), ψₜ_full)) -@show 1 - abs(inner(prod(ψₜ_krylov), ψₜ_full)) +main() diff --git a/test/runtests.jl b/test/runtests.jl index 3c5908e..b3423c9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ -using Test -using ITensorTDVP +using Test: @testset +using ITensorTDVP: ITensorTDVP test_path = joinpath(pkgdir(ITensorTDVP), "test") test_files = filter( @@ -8,7 +8,7 @@ test_files = filter( @testset "ITensorTDVP.jl" begin @testset "$filename" for filename in test_files println("Running $filename") - include(joinpath(test_path, filename)) + @time include(joinpath(test_path, filename)) end end diff --git a/test/tdvp_time_dependent.jl b/test/tdvp_time_dependent.jl new file mode 100644 index 0000000..b604760 --- /dev/null +++ b/test/tdvp_time_dependent.jl @@ -0,0 +1,77 @@ +# These tests take too long to compile, skip for now. +using ITensors: MPO, MPS, siteinds +using ITensorTDVP: ITensorTDVP, tdvp +using Test: @test, @testset + +include(joinpath(pkgdir(ITensorTDVP), "examples", "03_models.jl")) +include(joinpath(pkgdir(ITensorTDVP), "examples", "03_solvers.jl")) + +@testset "Time dependent Hamiltonian (eltype=$elt)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} +) + function ode_solver(H⃗₀, time_step, ψ₀; kwargs...) + ω₁ = typeof(time_step)(0.1) + ω₂ = typeof(time_step)(0.2) + ω⃗ = [ω₁, ω₂] + f⃗ = [t -> cos(ω * t) for ω in ω⃗] + ode_alg = Tsit5() + tol = √eps(real(time_step)) + ode_kwargs = (; reltol=tol, abstol=tol) + return ode_solver( + -im * TimeDependentSum(f⃗, H⃗₀), + time_step, + ψ₀; + solver_alg=ode_alg, + ode_kwargs..., + kwargs..., + ) + end + + function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...) + ω₁ = typeof(time_step)(0.1) + ω₂ = typeof(time_step)(0.2) + ω⃗ = [ω₁, ω₂] + f⃗ = [t -> cos(ω * t) for ω in ω⃗] + tol = √eps(real(time_step)) + krylov_kwargs = (; tol, eager=true) + return krylov_solver( + -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs... + ) + end + + n = 4 + J₁ = elt(1) + J₂ = elt(0.1) + time_step = real(elt)(0.1) + time_stop = real(elt)(1) + nsite = 2 + maxdim = 100 + cutoff = √(eps(real(elt))) + s = siteinds("S=1/2", n) + ℋ₁₀ = heisenberg(n; J=J₁, J2=elt(0)) + ℋ₂₀ = heisenberg(n; J=elt(0), J2=J₂) + ℋ⃗₀ = [ℋ₁₀, ℋ₂₀] + H⃗₀ = [MPO(elt, ℋ₀, s) for ℋ₀ in ℋ⃗₀] + ψ₀ = complex.(MPS(elt, s, j -> isodd(j) ? "↑" : "↓")) + ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite) + ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite) + ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀)) + + @test ITensors.scalartype(ψ₀) == complex(elt) + @test ITensors.scalartype(ψₜ_ode) == complex(elt) + @test ITensors.scalartype(ψₜ_krylov) == complex(elt) + @test ITensors.scalartype(ψₜ_full) == complex(elt) + @test norm(ψ₀) ≈ 1 + @test norm(ψₜ_ode) ≈ 1 + @test norm(ψₜ_krylov) ≈ 1 rtol = √(eps(real(elt))) + @test norm(ψₜ_full) ≈ 1 + + ode_err = norm(prod(ψₜ_ode) - ψₜ_full) + krylov_err = norm(prod(ψₜ_krylov) - ψₜ_full) + + @test krylov_err > ode_err + @test ode_err < √(eps(real(elt))) * 10^4 + @test krylov_err < √(eps(real(elt))) * 10^5 +end + +nothing diff --git a/test/test_linsolve.jl b/test/test_linsolve.jl index bf5bf65..32f3efe 100644 --- a/test/test_linsolve.jl +++ b/test/test_linsolve.jl @@ -40,7 +40,7 @@ using Random cutoff = 1e-5 maxdim = 20 solver_kwargs = (; tol=1e-4, maxiter=20, krylovdim=30, ishermitian=true) - x = @time linsolve(H, b, x0; nsweeps, cutoff, maxdim, solver_kwargs) + x = linsolve(H, b, x0; nsweeps, cutoff, maxdim, solver_kwargs) @test ITensors.scalartype(x) == eltype diff --git a/test/test_tdvp_time_dependent.jl b/test/test_tdvp_time_dependent.jl index 2df9b9a..efcd30f 100644 --- a/test/test_tdvp_time_dependent.jl +++ b/test/test_tdvp_time_dependent.jl @@ -1,76 +1,56 @@ -using ITensors: MPO, MPS, siteinds -using ITensorTDVP: ITensorTDVP, tdvp -using Test: @test, @testset - -include(joinpath(pkgdir(ITensorTDVP), "examples", "03_models.jl")) -include(joinpath(pkgdir(ITensorTDVP), "examples", "03_solvers.jl")) - -@testset "Time dependent Hamiltonian (eltype=$elt)" for elt in ( - Float32, Float64, Complex{Float32}, Complex{Float64} -) - function ode_solver(H⃗₀, time_step, ψ₀; kwargs...) - ω₁ = typeof(time_step)(0.1) - ω₂ = typeof(time_step)(0.2) - ω⃗ = [ω₁, ω₂] - f⃗ = [t -> cos(ω * t) for ω in ω⃗] - ode_alg = Tsit5() - tol = √eps(real(time_step)) - ode_kwargs = (; reltol=tol, abstol=tol) - return ode_solver( - -im * TimeDependentSum(f⃗, H⃗₀), - time_step, - ψ₀; - solver_alg=ode_alg, - ode_kwargs..., - kwargs..., - ) +using ITensors: + Index, MPO, ProjMPO, ProjMPOSum, QN, randomITensor, randomMPS, position!, siteinds +using ITensorTDVP: ITensorTDVP, TimeDependentSum, to_vec +using Test: @test, @test_skip, @testset + +@testset "TDVP with ODE local solver" begin + @testset "to_vec (eltype=$elt)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + space in (2, [QN(0) => 1, QN(1) => 1]) + + i, j = Index.((space, space)) + a = randomITensor(elt, i, j) + v, to_itensor = to_vec(a) + @test v isa Vector{elt} + @test vec(Array(a, i, j)) == v + @test to_itensor(v) == a end - - function krylov_solver(H⃗₀, time_step, ψ₀; kwargs...) - ω₁ = typeof(time_step)(0.1) - ω₂ = typeof(time_step)(0.2) - ω⃗ = [ω₁, ω₂] - f⃗ = [t -> cos(ω * t) for ω in ω⃗] - tol = √eps(real(time_step)) - krylov_kwargs = (; tol, eager=true) - return krylov_solver( - -im * TimeDependentSum(f⃗, H⃗₀), time_step, ψ₀; krylov_kwargs..., kwargs... + @testset "TimeDependentSum (eltype=$elt)" for elt in ( + Float32, Float64, Complex{Float32}, Complex{Float64} + ), + conserve_qns in [false, true] + + n = 4 + s = siteinds("S=1/2", 4; conserve_qns) + H = MPO(elt, s, "I") + H⃗ = [H, H] + region = 2:3 + ψ = randomMPS(elt, s, j -> isodd(j) ? "↑" : "↓"; linkdims=2) + H⃗ᵣ = ProjMPO.(H⃗) + map(Hᵣ -> position!(Hᵣ, ψ, first(region)), H⃗ᵣ) + ∑Hᵣ = ProjMPOSum(H⃗) + position!(∑Hᵣ, ψ, first(region)) + f⃗ₜ = [t -> sin(elt(0.1) * t), t -> cos(elt(0.2) * t)] + α = elt(0.5) + ∑Hₜ = α * TimeDependentSum(f⃗ₜ, ∑Hᵣ) + t₀ = elt(0.5) + ∑Hₜ₀ = ∑Hₜ(t₀) + ψᵣ = reduce(*, map(v -> ψ[v], region)) + Hψ = ∑Hₜ₀(ψᵣ) + @test eltype(Hψ) == elt + @test Hψ ≈ sum(i -> α * f⃗ₜ[i](t₀) * H⃗ᵣ[i](ψᵣ), eachindex(H⃗)) + end + @testset "Run example" begin + # These tests take too long to compile, skip for now. + @test_skip include( + joinpath(pkgdir(ITensorTDVP), "examples", "03_tdvp_time_dependent.jl") ) end - - n = 4 - J₁ = elt(1) - J₂ = elt(0.1) - time_step = real(elt)(0.1) - time_stop = real(elt)(1) - nsite = 2 - maxdim = 100 - cutoff = √(eps(real(elt))) - s = siteinds("S=1/2", n) - ℋ₁₀ = heisenberg(n; J=J₁, J2=elt(0)) - ℋ₂₀ = heisenberg(n; J=elt(0), J2=J₂) - ℋ⃗₀ = [ℋ₁₀, ℋ₂₀] - H⃗₀ = [MPO(elt, ℋ₀, s) for ℋ₀ in ℋ⃗₀] - ψ₀ = complex.(MPS(elt, s, j -> isodd(j) ? "↑" : "↓")) - ψₜ_ode = tdvp(ode_solver, H⃗₀, time_stop, ψ₀; time_step, maxdim, cutoff, nsite) - ψₜ_krylov = tdvp(krylov_solver, H⃗₀, time_stop, ψ₀; time_step, cutoff, nsite) - ψₜ_full, _ = ode_solver(prod.(H⃗₀), time_stop, prod(ψ₀)) - - @test ITensors.scalartype(ψ₀) == complex(elt) - @test ITensors.scalartype(ψₜ_ode) == complex(elt) - @test ITensors.scalartype(ψₜ_krylov) == complex(elt) - @test ITensors.scalartype(ψₜ_full) == complex(elt) - @test norm(ψ₀) ≈ 1 - @test norm(ψₜ_ode) ≈ 1 - @test norm(ψₜ_krylov) ≈ 1 rtol = √(eps(real(elt))) - @test norm(ψₜ_full) ≈ 1 - - ode_err = norm(prod(ψₜ_ode) - ψₜ_full) - krylov_err = norm(prod(ψₜ_krylov) - ψₜ_full) - - @test krylov_err > ode_err - @test ode_err < √(eps(real(elt))) * 10^4 - @test krylov_err < √(eps(real(elt))) * 10^5 + @testset "Time dependent TDVP" begin + # These tests take too long to compile, skip for now. + @test_skip include("tdvp_time_dependent.jl") + end end nothing