Change default method to direct solver and add tests

LorenzoFioroni · May 30, 2024 · 120369b · 120369b
1 parent ef0090a
commit 120369b
Show file tree

Hide file tree

Showing 2 changed files with 19 additions and 4 deletions.
diff --git a/src/steadystate.jl b/src/steadystate.jl
@@ -13,7 +13,7 @@ end
 
 function steadystate(
     L::QuantumObject{<:AbstractArray{T},SuperOperatorQuantumObject};
-    solver::SteadyStateSolver = SteadyStateLinearSolver(),
+    solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
 ) where {T}
     return _steadystate(L, solver; kwargs...)
@@ -22,7 +22,7 @@ end
 function steadystate(
     H::QuantumObject{<:AbstractArray{T},OpType},
     c_ops::AbstractVector;
-    solver::SteadyStateSolver = SteadyStateLinearSolver(),
+    solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
 ) where {T,OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
     L = liouvillian(H, c_ops)

diff --git a/test/steady_state.jl b/test/steady_state.jl
@@ -8,8 +8,23 @@
     psi0 = fock(N, 3)
     t_l = LinRange(0, 200, 1000)
     sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = false)
-    ρ_ss = steadystate(H, c_ops)
-    @test abs(sol_me.expect[1, end] - expect(e_ops[1], ρ_ss)) < 1e-3
+    rho_me = sol_me.states[end]
+
+    solver = SteadyStateDirectSolver()
+    ρ_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
+
+    solver = SteadyStateLinearSolver(alg=KrylovJL_GMRES())
+    ρ_ss = steadystate(H, c_ops, solver = solver)
+    @test tracedist(rho_me, ρ_ss) < 1e-4
+
+    solver = SteadyStateEigenSolver()
+    ρ_ss = steadystate(H, c_ops, solver = solver)
+    @test tracedist(rho_me, ρ_ss) < 1e-4
 
     H = a_d * a
     H_t = 0.1 * (a + a_d)