Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add const checking for non-dynamic evolution #371

Merged
merged 5 commits into from
Aug 1, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion src/bloch_redfield_master.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)]
Expand Down
10 changes: 10 additions & 0 deletions src/master.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
8 changes: 8 additions & 0 deletions src/mcwf.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
1 change: 1 addition & 0 deletions src/schroedinger.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
12 changes: 12 additions & 0 deletions src/stochastic_base.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,3 +92,15 @@
end
integrate_stoch(tspan, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function _check_const(op)
if !QuantumOpticsBase.is_const(op)
throw(

Check warning on line 98 in src/stochastic_base.jl

View check run for this annotation

Codecov / codecov/patch

src/stochastic_base.jl#L98

Added line #L98 was not covered by tests
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
4 changes: 3 additions & 1 deletion src/stochastic_master.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 2 additions & 0 deletions src/stochastic_schroedinger.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
11 changes: 11 additions & 0 deletions src/timeevolution_base.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
"""
Expand Down
7 changes: 7 additions & 0 deletions test/test_timeevolution_tdops.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down
Loading