Homec4science
Diffusion BosonicCodes (master) Browse

R10106/steady_state_solver.py
22796502c783master

Last Change
Raw File

steady_state_solver.py
View Options

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Mar 6 10:56:10 2020
@author: david
"""
import qutip
def solve_steady_state(h, c_ops, verbose = True, **kwargs):
"""
Solve the steadystate given the Hamiltonian h and a list of collapse
operators c_ops.
Parameters
----------
h : Qutip.Qobj
Hamiltonian of the system
c_ops : list of Qutip.Qobj ([Qutip.Qobj])
list of collapse operators in the Dissipator.
NOTE: when specifying the rate with which the collapse operators act
on the system (e.g. single-photon loss rate, etc.), each of the collapse
operators should be multiplied by this rate.
verbose : bool (optional, default: True)
If True: verbose mode on (for logging and debugging)
**kwargs : (optional)
additional parameters passed to the Qutip Solver
Returns
-------
rho_ss : Qutip.Qobj
Returns the steadystate density matrix
"""
# Define verbose printing function
vprint = print if verbose else lambda *a, **k: None
nbytes = h.data.data.nbytes + h.data.indptr.nbytes + h.data.indices.nbytes
vprint(sizeof_fmt(nbytes), " of memory consumption for Hamiltonian.")
if nbytes < 1024*1e4: # 10MiB of data for hamiltonian
#perform direct LU decomposition to invert the matrix directly
vprint("using direct Matrix inversion...")
# start_time = time.time()
res = qutip.steadystate(h, c_ops, method = 'direct',
return_info = verbose, **kwargs)
# vprint("--- %s seconds ---" % (time.time() - start_time))
return res
else :
# perform iterative biconjugate gradient method
vprint("using iterative bicstab with preconditioner")
# start_time = time.time()
res = qutip.steadystate(h, c_ops, method = 'iterative-bicstab',
return_info = verbose, use_precond=True)
# vprint("--- %s seconds ---" % (time.time() - start_time))
return res
def sizeof_fmt(num, suffix='B'):
"helper function to convert bytes to human readable"
for unit in ['','Ki','Mi','Gi','Ti','Pi','Ei','Zi']:
if abs(num) < 1024.0:
return "%3.1f%s%s" % (num, unit, suffix)
num /= 1024.0
return "%.1f%s%s" % (num, 'Yi', suffix)
c4science · Help