From a61bfa2e9024fdbeae98f54b4bd6f7c6c53f2d93 Mon Sep 17 00:00:00 2001 From: cailixun Date: Thu, 28 Aug 2025 14:27:53 +0800 Subject: [PATCH 1/5] Add steaystate ode kwargs passing --- src/steadystate.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/steadystate.jl b/src/steadystate.jl index a81fe3b3f..6baf658ab 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -217,6 +217,7 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESol save_everystep = false, saveat = ftype[], callback = cb, + kwargs..., ) ρss = sol.states[end] From 14717268a2084e1296c05a9690c70e03d8ee034a Mon Sep 17 00:00:00 2001 From: cailixun Date: Thu, 28 Aug 2025 15:00:22 +0800 Subject: [PATCH 2/5] add support of `QobjEvo` for `SteadyStateODESolver` --- CHANGELOG.md | 3 +++ src/steadystate.jl | 9 ++++++--- test/core-test/steady_state.jl | 8 ++++---- test/ext-test/cpu/autodiff/autodiff.jl | 15 +++++++++++++++ 4 files changed, 28 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a1e7269f1..f71b63114 100644 --- a/CHANGELOG.md +++ b/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) +- Add support of `QobjEvo` for `steadystate` (ODE solver only). ([#536]) + ## [v0.34.1] Release date: 2025-08-23 @@ -301,3 +303,4 @@ Release date: 2024-11-13 [#517]: https://github.com/qutip/QuantumToolbox.jl/issues/517 [#520]: https://github.com/qutip/QuantumToolbox.jl/issues/520 [#531]: https://github.com/qutip/QuantumToolbox.jl/issues/531 +[#536]: https://github.com/qutip/QuantumToolbox.jl/issues/536 diff --git a/src/steadystate.jl b/src/steadystate.jl index 6baf658ab..5ff7cd9cd 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -96,7 +96,7 @@ end @doc raw""" steadystate( - H::QuantumObject{OpType}, + H::AbstractQuantumObject{OpType}, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing; solver::SteadyStateSolver = SteadyStateDirectSolver(), kwargs..., @@ -111,7 +111,7 @@ Solve the stationary state based on different solvers. - `kwargs`: The keyword arguments for the solver. """ function steadystate( - H::QuantumObject{OpType}, + H::AbstractQuantumObject{OpType}, c_ops::Union{Nothing,AbstractVector,Tuple} = nothing; solver::SteadyStateSolver = SteadyStateDirectSolver(), kwargs..., @@ -199,7 +199,7 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateDirect return QuantumObject(ρss, Operator(), L.dimensions) end -function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...) +function _steadystate(L::AbstractQuantumObject{SuperOperator}, solver::SteadyStateODESolver; kwargs...) tmax = solver.tmax ψ0 = isnothing(solver.ψ0) ? rand_ket(L.dimensions) : solver.ψ0 @@ -224,6 +224,9 @@ function _steadystate(L::QuantumObject{SuperOperator}, solver::SteadyStateODESol return ρss end +_steadystate(L::QuantumObjectEvolution{SuperOperator}, solver::SteadyStateSolver; kwargs...) = + throw(ArgumentError("$(get_typename_wrapper(solver)) does not support QobjEvo.")) + struct SteadyStateODECondition{CT<:AbstractArray} cache::CT end diff --git a/test/core-test/steady_state.jl b/test/core-test/steady_state.jl index d78c93c33..6151af349 100644 --- a/test/core-test/steady_state.jl +++ b/test/core-test/steady_state.jl @@ -3,6 +3,7 @@ a = destroy(N) a_d = a' H = a_d * a + 0.1 * (a + a_d) + Ht = QobjEvo(H, (p, t) -> 1) # for test throw c_ops = [sqrt(0.1) * a] e_ops = [a_d * a] psi0 = fock(N, 3) @@ -17,18 +18,17 @@ solver = SteadyStateDirectSolver() ρ_ss = steadystate(H, c_ops, solver = solver) @test tracedist(rho_me, ρ_ss) < 1e-4 + @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver) solver = SteadyStateLinearSolver() ρ_ss = steadystate(H, c_ops, solver = solver) @test tracedist(rho_me, ρ_ss) < 1e-4 - - solver = SteadyStateLinearSolver() - ρ_ss = steadystate(H, c_ops, solver = solver) - @test tracedist(rho_me, ρ_ss) < 1e-4 + @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver) solver = SteadyStateEigenSolver() ρ_ss = steadystate(H, c_ops, solver = solver) @test tracedist(rho_me, ρ_ss) < 1e-4 + @test_throws ArgumentError steadystate(Ht, c_ops, solver = solver) @testset "Type Inference (steadystate)" begin L = liouvillian(H, c_ops) diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl index 3f8e601a4..4ba0a82b6 100644 --- a/test/ext-test/cpu/autodiff/autodiff.jl +++ b/test/ext-test/cpu/autodiff/autodiff.jl @@ -67,6 +67,17 @@ function my_f_mesolve(p) return real(expect(a' * a, sol.states[end])) end +function my_f_steadystate(p) + ρss = steadystate( + L, + SteadyStateODESolver(ψ0 = ψ0_mesolve, tmax = tlist_mesolve[end]); + params = p, + sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()), + ) + + return real(expect(a' * a, ρss)) +end + # Analytical solution n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2)) @@ -113,8 +124,12 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2)) my_f_mesolve_direct(params) my_f_mesolve(params) + my_f_steadystate(params) + # calculate exact solution and check if steadystate works grad_exact = Zygote.gradient((p) -> n_ss(p[1], p[2], p[3]), params)[1] + grad_ss = Zygote.gradient(my_f_steadystate, params)[1] + @test grad_ss ≈ grad_exact atol=1e-6 @testset "ForwardDiff.jl" begin grad_qt = ForwardDiff.gradient(my_f_mesolve_direct, params) From 7b66757f32afed469dfcba1445bb1dbec665d560 Mon Sep 17 00:00:00 2001 From: cailixun Date: Thu, 28 Aug 2025 15:31:37 +0800 Subject: [PATCH 3/5] fix typo --- test/ext-test/cpu/autodiff/autodiff.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl index 4ba0a82b6..8c490ed91 100644 --- a/test/ext-test/cpu/autodiff/autodiff.jl +++ b/test/ext-test/cpu/autodiff/autodiff.jl @@ -69,8 +69,8 @@ end function my_f_steadystate(p) ρss = steadystate( - L, - SteadyStateODESolver(ψ0 = ψ0_mesolve, tmax = tlist_mesolve[end]); + L; + solver = SteadyStateODESolver(ψ0 = ψ0_mesolve, tmax = tlist_mesolve[end]), params = p, sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()), ) From 34f74c3d31ce4507d217b497ca0743d65e750d89 Mon Sep 17 00:00:00 2001 From: cailixun Date: Thu, 28 Aug 2025 16:26:45 +0800 Subject: [PATCH 4/5] fix methods --- src/steadystate.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/steadystate.jl b/src/steadystate.jl index 5ff7cd9cd..7e55911ac 100644 --- a/src/steadystate.jl +++ b/src/steadystate.jl @@ -212,7 +212,8 @@ function _steadystate(L::AbstractQuantumObject{SuperOperator}, solver::SteadySta sol = mesolve( L, ψ0, - [ftype(0), ftype(tmax)], + [ftype(0), ftype(tmax)]; + alg = solver.alg, progress_bar = Val(false), save_everystep = false, saveat = ftype[], @@ -224,7 +225,11 @@ function _steadystate(L::AbstractQuantumObject{SuperOperator}, solver::SteadySta return ρss end -_steadystate(L::QuantumObjectEvolution{SuperOperator}, solver::SteadyStateSolver; kwargs...) = +_steadystate( + L::QuantumObjectEvolution{SuperOperator}, + solver::T; + kwargs..., +) where {T<:Union{SteadyStateDirectSolver,SteadyStateEigenSolver,SteadyStateLinearSolver}} = throw(ArgumentError("$(get_typename_wrapper(solver)) does not support QobjEvo.")) struct SteadyStateODECondition{CT<:AbstractArray} From ed7e22ba595ac187b077af5b8245a1d2a2116ba1 Mon Sep 17 00:00:00 2001 From: cailixun Date: Thu, 28 Aug 2025 20:27:18 +0800 Subject: [PATCH 5/5] remove autodiff tests for steadystate --- test/ext-test/cpu/autodiff/autodiff.jl | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/test/ext-test/cpu/autodiff/autodiff.jl b/test/ext-test/cpu/autodiff/autodiff.jl index 8c490ed91..3f8e601a4 100644 --- a/test/ext-test/cpu/autodiff/autodiff.jl +++ b/test/ext-test/cpu/autodiff/autodiff.jl @@ -67,17 +67,6 @@ function my_f_mesolve(p) return real(expect(a' * a, sol.states[end])) end -function my_f_steadystate(p) - ρss = steadystate( - L; - solver = SteadyStateODESolver(ψ0 = ψ0_mesolve, tmax = tlist_mesolve[end]), - params = p, - sensealg = BacksolveAdjoint(autojacvec = EnzymeVJP()), - ) - - return real(expect(a' * a, ρss)) -end - # Analytical solution n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2)) @@ -124,12 +113,8 @@ n_ss(Δ, F, γ) = abs2(F / (Δ + 1im * γ / 2)) my_f_mesolve_direct(params) my_f_mesolve(params) - my_f_steadystate(params) - # calculate exact solution and check if steadystate works grad_exact = Zygote.gradient((p) -> n_ss(p[1], p[2], p[3]), params)[1] - grad_ss = Zygote.gradient(my_f_steadystate, params)[1] - @test grad_ss ≈ grad_exact atol=1e-6 @testset "ForwardDiff.jl" begin grad_qt = ForwardDiff.gradient(my_f_mesolve_direct, params)