Skip to content

Commit

Permalink
More tweaks to tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@mtfishman
mtfishman committed Feb 16, 2024
1 parent a39321e commit daea7ab
Show file tree
Hide file tree
Showing 6 changed files with 268 additions and 208 deletions.
6 changes: 3 additions & 3 deletions examples/03_solvers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
265 changes: 134 additions & 131 deletions examples/03_tdvp_time_dependent.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
6 changes: 3 additions & 3 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
using Test
using ITensorTDVP
using Test: @testset
using ITensorTDVP: ITensorTDVP

test_path = joinpath(pkgdir(ITensorTDVP), "test")
test_files = filter(
Expand All @@ -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

Expand Down
77 changes: 77 additions & 0 deletions test/tdvp_time_dependent.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit daea7ab

Please sign in to comment.