minor changes

Author: ytdHuang

src/qobj/superoperators.jl

@@ -19,15 +19,9 @@ function _sprepost(A, B) # for any other input types
     return _spre(A, Id_cache) * _spost(B, Id_cache)
 end
 
-_liouvillian(H::AbstractMatrix, Id::AbstractMatrix) = -1im * (_spre(H, Id) - _spost(H, Id))
-function _lindblad_dissipator(O::AbstractMatrix, Id::AbstractMatrix)
-    Od_O = O' * O
-    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
-end
-
-## if input is AbstractSciMLOperator
-_lazy_tensor_warning(func_name::String, data::AbstractSciMLOperator) =
-    @warn "The function `$func_name` uses lazy tensor (which can hurt performance) for data type: $(get_typename_wrapper(data))"
+## if input is AbstractSciMLOperator 
+## some of them are optimzed 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) = mapreduce(op -> _spre(op, Id), +, A.ops)
@@ -44,10 +38,18 @@ function _spost(B::AbstractSciMLOperator, Id::AbstractMatrix)
     return kron(transpose(B), Id)
 end
 
+## intrinsic liouvillian 
+_liouvillian(H::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}} =
+    -1im * (_spre(H, Id) - _spost(H, Id))
 _liouvillian(H::MatrixOperator, Id::AbstractMatrix) = MatrixOperator(_liouvillian(H.A, Id))
 _liouvillian(H::ScaledOperator, Id::AbstractMatrix) = ScaledOperator(H.λ, _liouvillian(H.L, Id))
 _liouvillian(H::AddedOperator, Id::AbstractMatrix) = mapreduce(op -> _liouvillian(op, Id), +, H.ops)
-_liouvillian(H::AbstractSciMLOperator, Id::AbstractMatrix) = -1im * (_spre(H, Id) - _spost(H, Id))
+
+# intrinsic lindblad_dissipator
+function _lindblad_dissipator(O::MT, Id::AbstractMatrix) where {MT<:Union{AbstractMatrix,AbstractSciMLOperator}}
+    Od_O = O' * O
+    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
+end
 function _lindblad_dissipator(O::MatrixOperator, Id::AbstractMatrix)
     _O = O.A
     Od_O = _O' * _O
@@ -57,10 +59,6 @@ function _lindblad_dissipator(O::ScaledOperator, Id::AbstractMatrix)
     λc_λ = conj(O.λ) * O.λ
     return ScaledOperator(λc_λ, _lindblad_dissipator(O.L, Id))
 end
-function _lindblad_dissipator(O::AbstractSciMLOperator, Id::AbstractMatrix)
-    Od_O = O' * O
-    return _sprepost(O, O') - (_spre(Od_O, Id) + _spost(Od_O, Id)) / 2
-end
 
 @doc raw"""
     spre(A::AbstractQuantumObject, Id_cache=I(size(A,1)))

src/utilities.jl

@@ -153,6 +153,9 @@ _non_static_array_warning(argname, arg::AbstractVector{T}) where {T} =
           join(arg, ", ") *
           ")` instead of `$argname = $arg`." maxlog = 1
 
+_lazy_tensor_warning(func_name::String, data::AbstractSciMLOperator) =
+    @warn "The function `$func_name` uses lazy tensor (which can hurt performance) for data type: $(get_typename_wrapper(data))"
+
 # functions for getting Float or Complex element type
 _FType(::AbstractArray{T}) where {T<:Number} = _FType(T)
 _FType(::Type{Int32}) = Float32

test/core-test/quantum_objects_evo.jl

@@ -181,13 +181,18 @@
         @test isoper(H_td) == true
 
         # SuperOperator
+        X = a * a'
         c_op1 = QobjEvo(((a', coef1),))
-        c_op2 = QobjEvo(((a, coef2), (a * a', coef3)))
+        c_op2 = QobjEvo(((a, coef2), (X, coef3)))
         c_ops = [c_op1, c_op2]
         D1_ti = abs2(coef1(p, t)) * lindblad_dissipator(a')
-        D2_ti = @test_logs (:warn,) (:warn,) lindblad_dissipator(c_op2)(p, t)
+        D2_ti =
+            abs2(coef2(p, t)) * lindblad_dissipator(a) + # normal dissipator for first  element in c_op2
+            abs2(coef3(p, t)) * lindblad_dissipator(X) + # normal dissipator for second element in c_op2
+            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,) liouvillian(H_td, c_ops)
+        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
         # TODO: L_td here is ComposedOperator and need to setup cache first for the following test
@@ -199,6 +204,7 @@
         @test isconstant(L_td) == false
         @test issuper(L_td) == true
 
+        @test_logs (:warn,) (:warn,) liouvillian(H_td * H_td) # warnings from lazy tensor
         @test_throws MethodError QobjEvo([[a, coef1], a' * a, [a', coef2]])
         @test_throws ArgumentError H_td(ρvec, p, t)
         @test_throws ArgumentError L_td(ψ, p, t)

test/core-test/time_evolution.jl

@@ -229,6 +229,7 @@
         # @test sol_sse.expect ≈ sol_sse_td2.expect atol = 1e-2 * length(tlist)
 
         @testset "Type Inference mesolve" begin
+            coef(p, t) = exp(-t)
             @inferred mesolveProblem(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolveProblem(H, ψ0, [0, 10], c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolveProblem(
@@ -242,7 +243,7 @@
             @inferred mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false))
             @inferred mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false))
             @inferred mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, saveat = tlist, progress_bar = Val(false))
-            @inferred mesolve(H, ψ0, tlist, (a, a'), e_ops = (a' * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple
+            @inferred mesolve(H, ψ0, tlist, (a, QobjEvo(((a', coef),))), e_ops = (a' * a, a'), progress_bar = Val(false)) # We test the type inference for Tuple
             @inferred mesolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
             @inferred mesolve(H_td2, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)
             @inferred mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p)