Add const checking for non-dynamic evolution (#371)

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)]

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)

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)

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

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

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)

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)

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)
 """

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