Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Direct steadystate solver #338

Open
labay11 opened this issue Jul 26, 2022 · 1 comment
Open

Direct steadystate solver #338

labay11 opened this issue Jul 26, 2022 · 1 comment

Comments

@labay11
Copy link

labay11 commented Jul 26, 2022

I think it would be useful to include a direct steady-state solver which computes the solution of the system A x = 0 in an exact way, just like QuTip steady-state method.

For instance, using Pardiso.jl, the steady-state can be found as

function steadystate_mkl(L::SparseSuperOpType)
    n, m = length(L.basis_r[1]), length(L.basis_r[2])

    ps = MKLPardisoSolver()
    b = zeros(ComplexF64, n*m)

    # unity trace condition for the linear solvers
    w = mean(abs.(L.data))
    b[1] = w
    A = L.data + sparse(ones(n), 1:(n+1):n^2, fill(w, n), n*m, n*m)
    x = solve(ps, A, b)

    data = reshape(x, n, m)
    DenseOperator(L.basis_r[1], L.basis_r[2], data)
end

which is just slightly slower than the iterative method but requires fewer allocations and is more exact.

Here you can find a working example for the steady-state of a "pumped cavity": https://gist.github.com/labay11/66934e15e648d06d8511fc557c259e59

@PhilipVinc
Copy link
Contributor

PhilipVinc commented Oct 24, 2022

https://docs.qojulia.org/api/#QuantumOptics.steadystate.iterative

probably you can pass a function as an iterative method even if that solver is not iterative

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants