Add more steadystate and Dynamical Shifted Fock benchmarks (#557)

Author: albertomercurio

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Fix `cite()` bibtex output. ([#552])
 - Add `qeye_like` and `qzero_like`, which are synonyms of `one` and `zero`. ([#555])
+- Add steadystate and DSF benchmarks. The `SteadyStateODESOlver` tolerances are lowered to `terminate_reltol=1e-4` and `terminate_abstol=1e-6` to improve speed at the cost of accuracy. ([#557])
 
 ## [v0.36.0]
 Release date: 2025-09-29
@@ -328,3 +329,4 @@ Release date: 2024-11-13
 [#546]: https://github.com/qutip/QuantumToolbox.jl/issues/546
 [#552]: https://github.com/qutip/QuantumToolbox.jl/issues/552
 [#555]: https://github.com/qutip/QuantumToolbox.jl/issues/555
+[#557]: https://github.com/qutip/QuantumToolbox.jl/issues/557

benchmarks/Project.toml

@@ -1,5 +1,4 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 QuantumToolbox = "6c2fb7c5-b903-41d2-bc5e-5a7c320b9fab"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

benchmarks/dynamical_shifted_fock.jl

@@ -0,0 +1,81 @@
+function benchmark_dsf!(SUITE)
+    F = 2
+    Δ = 0.25
+    κ = 1
+    U = 0.01
+    J = 0.5
+    tlist = LinRange(0, 15, 300)
+
+    N = 5
+    a1 = kron(destroy(N), qeye(N))
+    a2 = kron(qeye(N), destroy(N))
+    function H_dsf(op_list, p)
+        Δ = p.Δ
+        F = p.F
+        U = p.U
+        J = p.J
+        a1, a2 = op_list
+        return Δ * a1' * a1 +
+               Δ * a2' * a2 +
+               U * a1'^2 * a1^2 +
+               U * a2'^2 * a2^2 +
+               F * (a1 + a1') +
+               J * (a1' * a2 + a1 * a2')
+    end
+    function c_ops_dsf(op_list, p)
+        κ = p.κ
+        a1, a2 = op_list
+        return [√κ * a1, √κ * a2]
+    end
+    function e_ops_dsf(op_list, p)
+        a1, a2 = op_list
+        return [a1' * a1, a2' * a2]
+    end
+    op_list = [a1, a2]
+    ψ0 = kron(fock(N, 0), fock(N, 0))
+    α0 = 1.5
+    α0_l = [α0, α0]
+    dsf_params = (Δ = Δ, F = F, κ = κ, U = U, J = J)
+
+    SUITE["Time Evolution"]["Dynamical Shifted Fock"]["mesolve"] = @benchmarkable dsf_mesolve(
+        $H_dsf,
+        $ψ0,
+        $tlist,
+        $c_ops_dsf,
+        $op_list,
+        $α0_l,
+        $dsf_params,
+        e_ops = $e_ops_dsf,
+        progress_bar = Val(false),
+    )
+
+    SUITE["Time Evolution"]["Dynamical Shifted Fock"]["mcsolve"]["Serial"] = @benchmarkable dsf_mcsolve(
+        $H_dsf,
+        $ψ0,
+        $tlist,
+        $c_ops_dsf,
+        $op_list,
+        $α0_l,
+        $dsf_params,
+        ntraj = 100,
+        e_ops = $e_ops_dsf,
+        progress_bar = Val(false),
+        ensemblealg = EnsembleSerial(),
+    )
+
+    SUITE["Time Evolution"]["Dynamical Shifted Fock"]["mcsolve"]["Multithreaded"] = @benchmarkable dsf_mcsolve(
+        $H_dsf,
+        $ψ0,
+        $tlist,
+        $c_ops_dsf,
+        $op_list,
+        $α0_l,
+        $dsf_params,
+        ntraj = 100,
+        e_ops = $e_ops_dsf,
+        progress_bar = Val(false),
+        ensemblealg = EnsembleThreads(),
+    )
+
+    return nothing
+end

benchmarks/runbenchmarks.jl

@@ -2,21 +2,22 @@ using BenchmarkTools
 using LinearAlgebra
 using SparseArrays
 using QuantumToolbox
-using OrdinaryDiffEq
-using LinearSolve
+using SciMLBase: EnsembleSerial, EnsembleThreads
 
 BLAS.set_num_threads(1)
 
 const SUITE = BenchmarkGroup()
 
 include("correlations_and_spectrum.jl")
 include("dynamical_fock_dimension.jl")
+include("dynamical_shifted_fock.jl")
 include("eigenvalues.jl")
 include("steadystate.jl")
 include("timeevolution.jl")
 
 benchmark_correlations_and_spectrum!(SUITE)
 benchmark_dfd!(SUITE)
+benchmark_dsf!(SUITE)
 benchmark_eigenvalues!(SUITE)
 benchmark_steadystate!(SUITE)
 benchmark_timeevolution!(SUITE)

benchmarks/steadystate.jl

@@ -11,5 +11,9 @@ function benchmark_steadystate!(SUITE)
 
     SUITE["Steadystate"]["Direct"] = @benchmarkable steadystate($H, $c_ops)
 
+    SUITE["Steadystate"]["Iterative"] = @benchmarkable steadystate($H, $c_ops; solver = SteadyStateLinearSolver())
+
+    SUITE["Steadystate"]["ODE"] = @benchmarkable steadystate($H, $c_ops; solver = SteadyStateODESolver())
+
     return nothing
 end

src/steadystate.jl

@@ -47,8 +47,8 @@ end
         alg = Tsit5(),
         ψ0 = nothing,
         tmax = Inf,
-        terminate_reltol = 1e-5,
-        terminate_abstol = 1e-7
+        terminate_reltol = 1e-4,
+        terminate_abstol = 1e-6
     )
 
 An ordinary differential equation (ODE) solver for solving [`steadystate`](@ref). It solves the stationary state based on [`mesolve`](@ref) with a termination condition.
@@ -69,8 +69,8 @@ or
 - `alg::OrdinaryDiffEqAlgorithm=Tsit5()`: The algorithm to solve the ODE.
 - `ψ0::Union{Nothing,QuantumObject}=nothing`: The initial state of the system. If not specified, a random pure state will be generated.
 - `tmax::Real=Inf`: The final time step for the steady state problem.
-- `terminate_reltol` = The relative tolerance for stationary state terminate condition. Default to `1e-5`.
-- `terminate_abstol` = The absolute tolerance for stationary state terminate condition. Default to `1e-7`.
+- `terminate_reltol` = The relative tolerance for stationary state terminate condition. Default to `1e-4`.
+- `terminate_abstol` = The absolute tolerance for stationary state terminate condition. Default to `1e-6`.
 
 !!! warning "Tolerances for terminate condition"
     The terminate condition tolerances `terminate_reltol` and `terminate_abstol` should be larger than `reltol` and `abstol` of [`mesolve`](@ref), respectively.
@@ -87,8 +87,8 @@ Base.@kwdef struct SteadyStateODESolver{
     alg::MT = Tsit5()
     ψ0::ST = nothing
     tmax::TT = Inf
-    terminate_reltol::RT = 10 * DEFAULT_ODE_SOLVER_OPTIONS.reltol
-    terminate_abstol::AT = 10 * DEFAULT_ODE_SOLVER_OPTIONS.abstol
+    terminate_reltol::RT = 1e-4
+    terminate_abstol::AT = 1e-6
 end
 
 @doc raw"""