Add more generic solver for steadystate_floquet (#396)

Author: albertomercurio

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Fix erroneous definition of the stochastic term in `smesolve`. ([#393])
 - Change name of `MultiSiteOperator` to `multisite_operator`. ([#394])
 - Fix `smesolve` for specifying initial state as density matrix. ([#395])
+- Add more generic solver for `steadystate_floquet` to allow more linear solvers. ([#396])
 - Fix time evolution output when using `saveat` keyword argument. ([#398])
 
 ## [v0.26.0]
@@ -125,4 +126,5 @@ Release date: 2024-11-13
 [#393]: https://github.com/qutip/QuantumToolbox.jl/issues/393
 [#394]: https://github.com/qutip/QuantumToolbox.jl/issues/394
 [#395]: https://github.com/qutip/QuantumToolbox.jl/issues/395
+[#396]: https://github.com/qutip/QuantumToolbox.jl/issues/396
 [#398]: https://github.com/qutip/QuantumToolbox.jl/issues/398

docs/src/resources/api.md

@@ -212,11 +212,12 @@ liouvillian_generalized
 
 ```@docs
 steadystate
-steadystate_floquet
+steadystate_fourier
 SteadyStateDirectSolver
 SteadyStateEigenSolver
 SteadyStateLinearSolver
 SteadyStateODESolver
+SSFloquetEffectiveLiouvillian
 ```
 
 ### [Dynamical Shifted Fock method](@id doc-API:Dynamical-Shifted-Fock-method)

docs/src/users_guide/steadystate.md

@@ -121,4 +121,4 @@ fig
 
 ## Calculate steady state for periodically driven systems
 
-See the docstring of [`steadystate_floquet`](@ref) for more details.
+See the docstring of [`steadystate_fourier`](@ref) for more details.

src/steadystate.jl

@@ -1,15 +1,12 @@
-export steadystate, steadystate_floquet
+export steadystate, steadystate_fourier, steadystate_floquet
 export SteadyStateSolver,
     SteadyStateDirectSolver,
     SteadyStateEigenSolver,
     SteadyStateLinearSolver,
     SteadyStateODESolver,
-    SteadyStateFloquetSolver,
-    SSFloquetLinearSystem,
     SSFloquetEffectiveLiouvillian
 
 abstract type SteadyStateSolver end
-abstract type SteadyStateFloquetSolver end
 
 @doc raw"""
     SteadyStateDirectSolver()
@@ -58,8 +55,18 @@ Base.@kwdef struct SteadyStateODESolver{MT<:OrdinaryDiffEqAlgorithm} <: SteadySt
     alg::MT = Tsit5()
 end
 
-struct SSFloquetLinearSystem <: SteadyStateFloquetSolver end
-Base.@kwdef struct SSFloquetEffectiveLiouvillian{SSST<:SteadyStateSolver} <: SteadyStateFloquetSolver
+@doc raw"""
+    SSFloquetEffectiveLiouvillian(steadystate_solver = SteadyStateDirectSolver())
+
+A solver which solves [`steadystate_fourier`](@ref) by first extracting an effective time-independent Liouvillian and then using the `steadystate_solver` to extract the steadystate..
+
+# Parameters
+- `steadystate_solver::SteadyStateSolver=SteadyStateDirectSolver()`: The solver to use for the effective Liouvillian.
+
+!!! note
+    This solver is only available for [`steadystate_fourier`](@ref).
+"""
+Base.@kwdef struct SSFloquetEffectiveLiouvillian{SSST<:SteadyStateSolver} <: SteadyStateSolver
     steadystate_solver::SSST = SteadyStateDirectSolver()
 end
 
@@ -85,6 +92,12 @@ function steadystate(
     solver::SteadyStateSolver = SteadyStateDirectSolver(),
     kwargs...,
 ) where {OpType<:Union{OperatorQuantumObject,SuperOperatorQuantumObject}}
+    solver isa SSFloquetEffectiveLiouvillian && throw(
+        ArgumentError(
+            "The solver `SSFloquetEffectiveLiouvillian` is only available for the `steadystate_fourier` function.",
+        ),
+    )
+
     L = liouvillian(H, c_ops)
 
     return _steadystate(L, solver; kwargs...)
@@ -247,15 +260,15 @@ function _steadystate_ode_condition(integrator, abstol, reltol, min_t)
 end
 
 @doc raw"""
-    steadystate_floquet(
-        H_0::QuantumObject{OpType1},
-        H_p::QuantumObject{OpType2},
-        H_m::QuantumObject{OpType3},
+    steadystate_fourier(
+        H_0::QuantumObject,
+        H_p::QuantumObject,
+        H_m::QuantumObject,
         ωd::Number,
         c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
         n_max::Integer = 2,
         tol::R = 1e-8,
-        solver::FSolver = SSFloquetLinearSystem,
+        solver::FSolver = SteadyStateLinearSolver(),
         kwargs...,
     )
 
@@ -265,11 +278,11 @@ Considering a monochromatic drive at frequency ``\omega_d``, we divide it into t
 `H_p` and `H_m`, where `H_p` oscillates
 as ``e^{i \omega t}`` and `H_m` oscillates as ``e^{-i \omega t}``.
 There are two solvers available for this function:
-- `SSFloquetLinearSystem`: Solves the linear system of equations.
+- `SteadyStateLinearSolver`: Solves the linear system of equations.
 - `SSFloquetEffectiveLiouvillian`: Solves the effective Liouvillian.
 For both cases, `n_max` is the number of Fourier components to consider, and `tol` is the tolerance for the solver.
 
-In the case of `SSFloquetLinearSystem`, the full linear system is solved at once:
+In the case of `SteadyStateLinearSolver`, the full linear system is solved at once:
 
 ```math
 ( \mathcal{L}_0 - i n \omega_d ) \hat{\rho}_n + \mathcal{L}_1 \hat{\rho}_{n-1} + \mathcal{L}_{-1} \hat{\rho}_{n+1} = 0
@@ -312,7 +325,10 @@ This will allow to simultaneously obtain all the ``\hat{\rho}_n``.
 In the case of `SSFloquetEffectiveLiouvillian`, instead, the effective Liouvillian is calculated using the matrix continued fraction method.
 
 !!! note "different return"
-    The two solvers returns different objects. The `SSFloquetLinearSystem` returns a list of [`QuantumObject`](@ref), containing the density matrices for each Fourier component (``\hat{\rho}_{-n}``, with ``n`` from ``0`` to ``n_\textrm{max}``), while the `SSFloquetEffectiveLiouvillian` returns only ``\hat{\rho}_0``. 
+    The two solvers returns different objects. The `SteadyStateLinearSolver` returns a list of [`QuantumObject`](@ref), containing the density matrices for each Fourier component (``\hat{\rho}_{-n}``, with ``n`` from ``0`` to ``n_\textrm{max}``), while the `SSFloquetEffectiveLiouvillian` returns only ``\hat{\rho}_0``. 
+
+!!! note
+    `steadystate_floquet` is a synonym of `steadystate_fourier`.
 
 ## Arguments
 - `H_0::QuantumObject`: The Hamiltonian or the Liouvillian of the undriven system.
@@ -322,38 +338,38 @@ In the case of `SSFloquetEffectiveLiouvillian`, instead, the effective Liouvilli
 - `c_ops::Union{Nothing,AbstractVector} = nothing`: The optional collapse operators.
 - `n_max::Integer = 2`: The number of Fourier components to consider.
 - `tol::R = 1e-8`: The tolerance for the solver.
-- `solver::FSolver = SSFloquetLinearSystem`: The solver to use.
+- `solver::FSolver = SteadyStateLinearSolver`: The solver to use.
 - `kwargs...`: Additional keyword arguments to be passed to the solver.
 """
-function steadystate_floquet(
+function steadystate_fourier(
     H_0::QuantumObject{OpType1},
     H_p::QuantumObject{OpType2},
     H_m::QuantumObject{OpType3},
     ωd::Number,
     c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
     n_max::Integer = 2,
     tol::R = 1e-8,
-    solver::FSolver = SSFloquetLinearSystem(),
+    solver::FSolver = SteadyStateLinearSolver(),
     kwargs...,
 ) where {
     OpType1<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     OpType2<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     OpType3<:Union{OperatorQuantumObject,SuperOperatorQuantumObject},
     R<:Real,
-    FSolver<:SteadyStateFloquetSolver,
+    FSolver<:SteadyStateSolver,
 }
     L_0 = liouvillian(H_0, c_ops)
     L_p = liouvillian(H_p)
     L_m = liouvillian(H_m)
-    return _steadystate_floquet(L_0, L_p, L_m, ωd, solver; n_max = n_max, tol = tol, kwargs...)
+    return _steadystate_fourier(L_0, L_p, L_m, ωd, solver; n_max = n_max, tol = tol, kwargs...)
 end
 
-function _steadystate_floquet(
+function _steadystate_fourier(
     L_0::QuantumObject{SuperOperatorQuantumObject},
     L_p::QuantumObject{SuperOperatorQuantumObject},
     L_m::QuantumObject{SuperOperatorQuantumObject},
     ωd::Number,
-    solver::SSFloquetLinearSystem;
+    solver::SteadyStateLinearSolver;
     n_max::Integer = 1,
     tol::R = 1e-8,
     kwargs...,
@@ -387,9 +403,18 @@ function _steadystate_floquet(
     v0 = zeros(T, n_fourier * N)
     v0[n_max*N+1] = weight
 
-    Pl = ilu(M, τ = 0.01)
+    (haskey(kwargs, :Pl) || haskey(kwargs, :Pr)) && error("The use of preconditioners must be defined in the solver.")
+    if !isnothing(solver.Pl)
+        kwargs = merge((; kwargs...), (Pl = solver.Pl(M),))
+    elseif isa(M, SparseMatrixCSC)
+        kwargs = merge((; kwargs...), (Pl = ilu(M, τ = 0.01),))
+    end
+    !isnothing(solver.Pr) && (kwargs = merge((; kwargs...), (Pr = solver.Pr(M),)))
+    !haskey(kwargs, :abstol) && (kwargs = merge((; kwargs...), (abstol = tol,)))
+    !haskey(kwargs, :reltol) && (kwargs = merge((; kwargs...), (reltol = tol,)))
+
     prob = LinearProblem(M, v0)
-    ρtot = solve(prob, KrylovJL_GMRES(), Pl = Pl, abstol = tol, reltol = tol).u
+    ρtot = solve(prob, solver.alg; kwargs...).u
 
     offset1 = n_max * N
     offset2 = (n_max + 1) * N
@@ -409,7 +434,7 @@ function _steadystate_floquet(
     return ρ_list
 end
 
-function _steadystate_floquet(
+function _steadystate_fourier(
     L_0::QuantumObject{SuperOperatorQuantumObject},
     L_p::QuantumObject{SuperOperatorQuantumObject},
     L_m::QuantumObject{SuperOperatorQuantumObject},
@@ -425,3 +450,6 @@ function _steadystate_floquet(
 
     return steadystate(L_eff; solver = solver.steadystate_solver, kwargs...)
 end
+
+# TODO: Synonym to align with QuTiP. Track https://github.com/qutip/qutip/issues/2632 when this can be removed.
+const steadystate_floquet = steadystate_fourier

test/core-test/steady_state.jl

@@ -65,16 +65,23 @@
     H_td = (H, (H_t, coeff))
 
     sol_me = mesolve(H_td, psi0, t_l, c_ops, e_ops = e_ops, progress_bar = Val(false))
-    ρ_ss1 = steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetLinearSystem())[1]
+    ρ_ss1 = steadystate_fourier(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SteadyStateLinearSolver())[1]
     ρ_ss2 =
-        steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian())
+        steadystate_fourier(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetEffectiveLiouvillian())
 
     @test abs(sum(sol_me.expect[1, end-100:end]) / 101 - expect(e_ops[1], ρ_ss1)) < 1e-3
     @test abs(sum(sol_me.expect[1, end-100:end]) / 101 - expect(e_ops[1], ρ_ss2)) < 1e-3
 
-    @testset "Type Inference (steadystate_floquet)" begin
-        @inferred steadystate_floquet(H, -1im * 0.5 * H_t, 1im * 0.5 * H_t, 1, c_ops, solver = SSFloquetLinearSystem())
-        @inferred steadystate_floquet(
+    @testset "Type Inference (steadystate_fourier)" begin
+        @inferred steadystate_fourier(
+            H,
+            -1im * 0.5 * H_t,
+            1im * 0.5 * H_t,
+            1,
+            c_ops,
+            solver = SteadyStateLinearSolver(),
+        )
+        @inferred steadystate_fourier(
             H,
             -1im * 0.5 * H_t,
             1im * 0.5 * H_t,