Return sesolve when mesolve allows it (#455)

Author: albertomercurio

CHANGELOG.md

@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased](https://github.com/qutip/QuantumToolbox.jl/tree/main)
 
+- Return `sesolve` when `mesolve` allows it. ([#455])
+
 ## [v0.30.1]
 Release date: 2025-04-24
 
@@ -212,3 +214,4 @@ Release date: 2024-11-13
 [#448]: https://github.com/qutip/QuantumToolbox.jl/issues/448
 [#450]: https://github.com/qutip/QuantumToolbox.jl/issues/450
 [#453]: https://github.com/qutip/QuantumToolbox.jl/issues/453
+[#455]: https://github.com/qutip/QuantumToolbox.jl/issues/455

src/time_evolution/mesolve.jl

@@ -46,6 +46,7 @@ where
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
+- If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolveProblem`](@ref) instead.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
 
@@ -67,6 +68,17 @@ function mesolveProblem(
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject,OperatorKetQuantumObject},
 }
+    (isoper(H) && isket(ψ0) && isnothing(c_ops)) && return sesolveProblem(
+        H,
+        ψ0,
+        tlist;
+        e_ops = e_ops,
+        params = params,
+        progress_bar = progress_bar,
+        inplace = inplace,
+        kwargs...,
+    )
+
     haskey(kwargs, :save_idxs) &&
         throw(ArgumentError("The keyword argument \"save_idxs\" is not supported in QuantumToolbox."))
 
@@ -141,6 +153,7 @@ where
 
 - The states will be saved depend on the keyword argument `saveat` in `kwargs`.
 - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately.
+- If `H` is an [`Operator`](@ref), `ψ0` is a [`Ket`](@ref) and `c_ops` is `Nothing`, the function will call [`sesolve`](@ref) instead.
 - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`.
 - For more details about `alg` please refer to [`DifferentialEquations.jl` (ODE Solvers)](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
 - For more details about `kwargs` please refer to [`DifferentialEquations.jl` (Keyword Arguments)](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/)
@@ -164,6 +177,18 @@ function mesolve(
     HOpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     StateOpType<:Union{KetQuantumObject,OperatorQuantumObject,OperatorKetQuantumObject},
 }
+    (isoper(H) && isket(ψ0) && isnothing(c_ops)) && return sesolve(
+        H,
+        ψ0,
+        tlist;
+        alg = alg,
+        e_ops = e_ops,
+        params = params,
+        progress_bar = progress_bar,
+        inplace = inplace,
+        kwargs...,
+    )
+
     prob = mesolveProblem(
         H,
         ψ0,

test/core-test/time_evolution.jl

@@ -183,6 +183,9 @@
             rng = MersenneTwister(12),
         )
 
+        # Redirect to `sesolve`
+        sol_me5 = mesolve(H, ψ0, tlist, progress_bar = Val(false))
+
         ρt_mc = [ket2dm.(normalize.(states)) for states in sol_mc_states.states]
         expect_mc_states = mapreduce(states -> expect.(Ref(e_ops[1]), states), hcat, ρt_mc)
         expect_mc_states_mean = sum(expect_mc_states, dims = 2) / size(expect_mc_states, 2)
@@ -198,6 +201,7 @@
         sol_sme_string = sprint((t, s) -> show(t, "text/plain", s), sol_sme)
         @test prob_me.prob.f.f isa MatrixOperator
         @test prob_mc.prob.f.f isa MatrixOperator
+        @test isket(sol_me5.states[1])
         @test sum(abs, sol_mc.expect .- sol_me.expect) / length(tlist) < 0.1
         @test sum(abs, sol_mc2.expect .- sol_me.expect) / length(tlist) < 0.1
         @test sum(abs, vec(expect_mc_states_mean) .- vec(sol_me.expect[1, saveat_idxs])) / length(tlist) < 0.1