From 9c025b5c602d413ba5169a388b75dcb4ebc7f08b Mon Sep 17 00:00:00 2001 From: Ashley Milsted Date: Tue, 1 Aug 2023 13:46:29 -0700 Subject: [PATCH] Add const checking for non-dynamic evolution (#371) --- src/bloch_redfield_master.jl | 3 ++- src/master.jl | 10 ++++++++++ src/mcwf.jl | 8 ++++++++ src/schroedinger.jl | 1 + src/stochastic_base.jl | 12 ++++++++++++ src/stochastic_master.jl | 4 +++- src/stochastic_schroedinger.jl | 2 ++ src/timeevolution_base.jl | 11 +++++++++++ test/test_timeevolution_tdops.jl | 7 +++++++ 9 files changed, 56 insertions(+), 2 deletions(-) diff --git a/src/bloch_redfield_master.jl b/src/bloch_redfield_master.jl index cb33f9cb..4c68265c 100644 --- a/src/bloch_redfield_master.jl +++ b/src/bloch_redfield_master.jl @@ -16,7 +16,8 @@ See QuTiP's documentation (http://qutip.org/docs/latest/guide/dynamics/dynamics- This argument is only taken into account if use_secular=true. """ function bloch_redfield_tensor(H::AbstractOperator, a_ops; J=SparseOpType[], use_secular=true, secular_cutoff=0.1) - + _check_const(H) + _check_const.(J) # Use the energy eigenbasis H_evals, transf_mat = eigen(Array(H.data)) #Array call makes sure H is a dense array H_ekets = [Ket(H.basis_l, transf_mat[:, i]) for i in 1:length(H_evals)] diff --git a/src/master.jl b/src/master.jl index 3f6a6be5..f3fc3e8a 100644 --- a/src/master.jl +++ b/src/master.jl @@ -10,6 +10,9 @@ function master_h(tspan, rho0::Operator, H::AbstractOperator, J; Jdagger=dagger.(J), fout=nothing, kwargs...) + _check_const(H) + _check_const.(J) + _check_const.(Jdagger) check_master(rho0, H, J, Jdagger, rates) tmp = copy(rho0) dmaster_(t, rho, drho) = dmaster_h!(drho, H, J, Jdagger, rates, rho, tmp) @@ -33,6 +36,10 @@ function master_nh(tspan, rho0::Operator, Hnh::AbstractOperator, J; Jdagger=dagger.(J), fout=nothing, kwargs...) + _check_const(Hnh) + _check_const(Hnhdagger) + _check_const.(J) + _check_const.(Jdagger) check_master(rho0, Hnh, J, Jdagger, rates) tmp = copy(rho0) dmaster_(t, rho, drho) = dmaster_nh!(drho, Hnh, Hnhdagger, J, Jdagger, rates, rho, tmp) @@ -76,6 +83,9 @@ function master(tspan, rho0::Operator, H::AbstractOperator, J; Jdagger=dagger.(J), fout=nothing, kwargs...) + _check_const(H) + _check_const.(J) + _check_const.(Jdagger) isreducible = check_master(rho0, H, J, Jdagger, rates) if !isreducible tmp = copy(rho0) diff --git a/src/mcwf.jl b/src/mcwf.jl index 5e2b481b..1e1646e0 100644 --- a/src/mcwf.jl +++ b/src/mcwf.jl @@ -18,6 +18,9 @@ function mcwf_h(tspan, psi0::Ket, H::AbstractOperator, J; tmp=copy(psi0), display_beforeevent=false, display_afterevent=false, kwargs...) + _check_const(H) + _check_const.(J) + _check_const.(Jdagger) check_mcwf(psi0, H, J, Jdagger, rates) f(t, psi, dpsi) = dmcwf_h!(dpsi, H, J, Jdagger, rates, psi, tmp) j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, rates) @@ -42,6 +45,8 @@ function mcwf_nh(tspan, psi0::Ket, Hnh::AbstractOperator, J; seed=rand(UInt), fout=nothing, display_beforeevent=false, display_afterevent=false, kwargs...) + _check_const(Hnh) + _check_const.(J) check_mcwf(psi0, Hnh, J, J, nothing) f(t, psi, dpsi) = dschroedinger!(dpsi, Hnh, psi) j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new, nothing) @@ -96,6 +101,9 @@ function mcwf(tspan, psi0::Ket, H::AbstractOperator, J; fout=nothing, Jdagger=dagger.(J), display_beforeevent=false, display_afterevent=false, kwargs...) + _check_const(H) + _check_const.(J) + _check_const.(Jdagger) isreducible = check_mcwf(psi0, H, J, Jdagger, rates) if !isreducible tmp = copy(psi0) diff --git a/src/schroedinger.jl b/src/schroedinger.jl index c5aa840d..ede39534 100644 --- a/src/schroedinger.jl +++ b/src/schroedinger.jl @@ -15,6 +15,7 @@ Integrate Schroedinger equation to evolve states or compute propagators. function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing, kwargs...) where {B,Bo,T<:Union{AbstractOperator{B,Bo},StateVector{B}}} + _check_const(H) dschroedinger_(t, psi, dpsi) = dschroedinger!(dpsi, H, psi) tspan, psi0 = _promote_time_and_state(psi0, H, tspan) # promote only if ForwardDiff.Dual x0 = psi0.data diff --git a/src/stochastic_base.jl b/src/stochastic_base.jl index 339ac3fc..80f2fa8c 100644 --- a/src/stochastic_base.jl +++ b/src/stochastic_base.jl @@ -92,3 +92,15 @@ function integrate_stoch(tspan, df, dg, x0, end integrate_stoch(tspan, df, dg, x0, state, dstate, fout, n; kwargs...) end + +function _check_const(op) + if !QuantumOpticsBase.is_const(op) + throw( + ArgumentError("You are attempting to use a time-dependent dynamics generator " * + "(a Hamiltonian or Lindbladian) with a solver that assumes constant " * + "dynamics. To avoid errors, please use the _dynamic solvers instead, " * + "e.g. schroedinger_dynamic instead of schroedinger") + ) + end + nothing +end diff --git a/src/stochastic_master.jl b/src/stochastic_master.jl index 42d0e4a9..88baff30 100644 --- a/src/stochastic_master.jl +++ b/src/stochastic_master.jl @@ -35,7 +35,9 @@ function master(tspan, rho0::T, H::AbstractOperator{B,B}, Jdagger=dagger.(J), Cdagger=dagger.(C), fout=nothing, kwargs...) where {B,T<:Operator{B,B}} - + _check_const(H) + _check_const.(J) + _check_const.(Jdagger) tmp = copy(rho0) n = length(C) diff --git a/src/stochastic_schroedinger.jl b/src/stochastic_schroedinger.jl index 18fc7582..bd8066b2 100644 --- a/src/stochastic_schroedinger.jl +++ b/src/stochastic_schroedinger.jl @@ -24,6 +24,8 @@ function schroedinger(tspan, psi0::T, H::AbstractOperator{B,B}, Hs::Vector; normalize_state=false, calback=nothing, kwargs...) where {B,T<:Ket{B}} + _check_const(H) + _check_const.(Hs) tspan_ = convert(Vector{float(eltype(tspan))}, tspan) n = length(Hs) diff --git a/src/timeevolution_base.jl b/src/timeevolution_base.jl index 43ea6b77..3e8c270c 100644 --- a/src/timeevolution_base.jl +++ b/src/timeevolution_base.jl @@ -88,6 +88,17 @@ function (c::SteadyStateCondtion)(rho,t,integrator) drho/dt < c.tol end +function _check_const(op) + if !QuantumOpticsBase.is_const(op) + throw( + ArgumentError("You are attempting to use a time-dependent dynamics generator " * + "(a Hamiltonian or Lindbladian) with a solver that assumes constant " * + "dynamics. To avoid errors, please use the _dynamic solvers instead, " * + "e.g. schroedinger_dynamic instead of schroedinger") + ) + end + nothing +end const QO_CHECKS = Ref(true) """ diff --git a/test/test_timeevolution_tdops.jl b/test/test_timeevolution_tdops.jl index de76b6e4..1d45faad 100644 --- a/test/test_timeevolution_tdops.jl +++ b/test/test_timeevolution_tdops.jl @@ -19,6 +19,8 @@ _getf = (H0, Hd) -> (t,_) -> _h(t, H0, Hd) fman = _getf(H0, Hd) psi0 = basisstate(b, 1) + +@test_throws ArgumentError timeevolution.schroedinger(ts, psi0, H) ts_out, psis = timeevolution.schroedinger_dynamic(ts, psi0, H) # check this is not trivial @test !(psis[1].data ≈ psis[end].data) @@ -29,6 +31,7 @@ ts_out2, psis2 = timeevolution.schroedinger_dynamic(ts, psi0, fman) _getf = (H0, Hd, a) -> (t,_) -> (_h(t, H0, Hd), (), ()) fman = _getf(H0, Hd, a) +@test_throws ArgumentError timeevolution.master(ts, psi0, H, []) ts_out, rhos = timeevolution.master_dynamic(ts, psi0, H, []) ts_out2, rhos2 = timeevolution.master_dynamic(ts, psi0, fman) @test rhos[end].data ≈ rhos2[end].data @@ -41,10 +44,12 @@ _js(t, a) = (cos(t)*a, 0.01*a', 0.01*a'*a) _getf = (H0, Hd, a) -> (t,_) -> (_h(t, H0, Hd), _js(t, a), dagger.(_js(t, a))) fman = _getf(H0, Hd, a) +@test_throws ArgumentError timeevolution.mcwf(ts, psi0, H, Js; seed=0) ts_out, psis = timeevolution.mcwf_dynamic(ts, psi0, H, Js; seed=0) ts_out2, psis2 = timeevolution.mcwf_dynamic(ts, psi0, fman; seed=0) @test psis[end].data ≈ psis2[end].data +@test_throws ArgumentError timeevolution.master(ts, psi0, H, Js) ts_out, rhos = timeevolution.master_dynamic(ts, psi0, H, Js) ts_out2, rhos2 = timeevolution.master_dynamic(ts, psi0, fman) @test rhos[end].data ≈ rhos2[end].data @@ -58,6 +63,7 @@ _getf = (H0, Hd, a) -> (t,_) -> ( fman = _getf(H0, Hd, a) +@test_throws ArgumentError timeevolution.mcwf_nh(ts, psi0, Hnh, Js; seed=0) ts_out, psis = timeevolution.mcwf_nh_dynamic(ts, psi0, Hnh, Js; seed=0) ts_out2, psis2 = timeevolution.mcwf_nh_dynamic(ts, psi0, fman; seed=0) @test psis[end].data ≈ psis2[end].data @@ -70,6 +76,7 @@ _getf = (H0, Hd, a) -> (t,_) -> ( fman = _getf(H0, Hd, a) +@test_throws ArgumentError timeevolution.master_nh(ts, psi0, Hnh, Js) ts_out, rhos = timeevolution.master_nh_dynamic(ts, psi0, Hnh, Js) ts_out2, rhos2 = timeevolution.master_nh_dynamic(ts, psi0, fman) @test rhos[end].data ≈ rhos2[end].data