qojulia / QuantumOptics.jl

Library for the numerical simulation of closed as well as open quantum systems.
http://qojulia.org
Other
543 stars 110 forks source link

Cache input to master derivative functions #406

Open apkille opened 3 months ago

apkille commented 3 months ago

In several of the master functions, e.g., dmaster_h!, a cached copy of drho is required as input:

dmaster_h!(drho, H, J, Jdagger, rates, rho, drho_cache)

Presumably, this is to improve performance of solvers such as timeevolution.master. However, this feels excessive if we don't want to use the corresponding solver, and instead create an ODE problem directly with our density operators and dmaster_h!. For example, it would be great if we could define an in-place method without the cache input, say

f!(drho, rho, p, t) = timeevolution.dmaster_h!(drho, H, J, Jdagger, rates, rho)
prob  = ODEProblem(f!, rho0, (t0, t1))

This is a natural interface that I'm aiming for with all of the recent broadcasting work (#404, https://github.com/qojulia/QuantumOpticsBase.jl/pull/172). Do the QO maintainers have any comments or oppositions to defining such a method in QO.jl?

david-pl commented 3 months ago

It shouldn't create any ambiguities to add that method, so go ahead :+1:

Krastanov commented 2 months ago

this might be related https://github.com/qojulia/QuantumOptics.jl/issues/298

apkille commented 2 months ago

@krastanov to build off our discussion in person: a workaround is assigning the cached drho to the p parameter of the ODEProblem object, which prevents us from caching at each time step if we use OrdinaryDiffEq.jl directly. So we don't have to define any new methods. For example, using dmaster_h! which is already a part of the public API, we have

using QuantumOptics, OrdinaryDiffEq

b = FockBasis(5)
psi0 = fockstate(b, 4)
rho0 = dm(psi0)
tmp = copy(rho0)
H = number(b)
J = [destroy(b), create(b)]
Jdagger = dagger.(J)
rates = [0.7, 0.5]

f!(drho, rho, p, t) = timeevolution.dmaster_h!(drho, H, J, Jdagger, rates, rho, p)
prob = ODEProblem(f!, rho0, (0.0, 1.0), tmp)
sol = solve(prob, DP5(); save_everystep=false)