fix QobjEvo in-place operator evaluation

Author: ytdHuang

src/qobj/quantum_object_evo.jl

@@ -513,7 +513,11 @@ function (A::QuantumObjectEvolution)(
     end
 
     # We put `nothing` in the place of `u` because the time-dependence doesn't depend on the state
-    A.data(ψout.data, ψin.data, nothing, p, t)
+    ## TODO: should use SciMLOperators v1 API: A.data(ψout.data, ψin.data, nothing, p, t)
+    ## But somehow that doesn't work for some specific cases in our runtests.
+    ## Therefore, we use the most basic way to do it for now
+    update_coefficients!(A.data, nothing, p, t)
+    mul!(ψout.data, A.data, ψin.data)
 
     return ψout
 end

src/qobj/superoperators.jl

@@ -14,32 +14,31 @@ _sprepost(A::AbstractMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), spa
 _sprepost(A::AbstractMatrix, B::AbstractSparseMatrix) = kron(transpose(B), sparse(A))
 _sprepost(A::AbstractSparseMatrix, B::AbstractMatrix) = kron(transpose(sparse(B)), A)
 _sprepost(A::AbstractSparseMatrix, B::AbstractSparseMatrix) = kron(transpose(B), A)
+function _sprepost(A, B) # for any other input types
+    Id_cache = I(size(A, 1))
+    return _spre(A, Id_cache) * _spost(B, Id_cache)
+end
 
 ## if input is AbstractSciMLOperator 
 ## some of them are optimized to speed things up
 ## the rest of the SciMLOperators will just use lazy tensor (and prompt a warning)
 _spre(A::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spre(A.A, Id))
 _spre(A::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(A.λ, _spre(A.L, Id))
 _spre(A::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _spre(op, Id), A.ops))
-function _spre(A::AbstractSciMLOperator, Id::Union{AbstractMatrix,AbstractSciMLOperator})
+function _spre(A::AbstractSciMLOperator, Id::AbstractMatrix)
     _lazy_tensor_warning(Id, A)
     return kron(Id, A)
 end
 
 _spost(B::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_spost(B.A, Id))
 _spost(B::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(B.λ, _spost(B.L, Id))
 _spost(B::AddedOperator, Id::AbstractMatrix) = AddedOperator(map(op -> _spost(op, Id), B.ops))
-function _spost(B::AbstractSciMLOperator, Id::Union{AbstractMatrix,AbstractSciMLOperator})
+function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     B_T = transpose(B)
     _lazy_tensor_warning(B_T, Id)
     return kron(B_T, Id)
 end
 
-function _sprepost(A::AbstractSciMLOperator, B::AbstractSciMLOperator)
-    Id_cache = IdentityOperator(size(A, 1))
-    return _spre(A, Id_cache) * _spost(B, Id_cache)
-end
-
 ## intrinsic liouvillian 
 _liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
     -1im * (_spre(H, Id) - _spost(H, Id))

test/core-test/quantum_objects_evo.jl

@@ -215,7 +215,7 @@
             coef2(p, t) * conj(coef3(p, t)) * (spre(a) * spost(X') - 0.5 * spre(X' * a) - 0.5 * spost(X' * a)) + # cross terms
             conj(coef2(p, t)) * coef3(p, t) * (spre(X) * spost(a') - 0.5 * spre(a' * X) - 0.5 * spost(a' * X))   # cross terms
         L_ti = liouvillian(H_ti) + D1_ti + D2_ti
-        L_td = @test_logs (:warn,) (:warn,) (:warn,) (:warn,) liouvillian(H_td, c_ops) # warnings from lazy tensor in `lindblad_dissipator(c_op2)`
+        L_td = @test_logs (:warn,) (:warn,) liouvillian(H_td, c_ops) # warnings from lazy tensor in `lindblad_dissipator(c_op2)`
         ρvec = mat2vec(rand_dm(N))
         @test L_td(p, t) ≈ L_ti
         @test iscached(L_td) == false