High allocations in timevolution.schroedinger_dynamic and master_dynamic

[mgoloshchapov]

Hi! I've recently installed QuantumOptics.jl and tried benchmarking solution of qubit dynamics with time-dependent field.

$$H(t) = \frac{\Omega(t)}{2} (\left|g\right>\left< e \right| + \left| e \right> \left< g \right| )$$

I receive enormous amount of allocations when using timeevolution.schroedinger_dynamic compared to explicit implementation using DifferentialEquations.jl. Can you please tell me what I'm doing wrong?

#Packages
using QuantumOptics
using BenchmarkTools
using StaticArrays
using DifferentialEquations, OrdinaryDiffEq

#Basis and operators
const basis = NLevelBasis(2);
const g = nlevelstate(basis, 1);
const e = nlevelstate(basis, 2);
const σge = g ⊗ dagger(e);
const σeg = e ⊗ dagger(g);

Implementations in QuantumOptics.jl that I tried:

const H2(t, psi) = TimeDependentSum([1.0 .+ 0.1*sin.(t), 1.0 .+ 0.1*sin.(t)], [σge, σeg])

function H3(t, psi)
    return (1.0 + 0.1*sin(t)) * (σge + σeg)
end;

function H4(t, psi)
    return (1.0 + 0.1*sin(t)) * (σge + σeg)
end;

Results of benchmarking:

ψ0 = g;
tspan = [0.0:0.01:100.0;];

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H2; alg=Tsit5());

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H3; alg=Tsit5());

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H4; alg=Tsit5());

3.293 ms (42553 allocations: 3.53 MiB)

2.821 ms (32929 allocations: 3.09 MiB)

2.849 ms (32929 allocations: 3.09 MiB)

Implementations in DifferentialEquations.jl

function rabi_problem1(ψ, params, t)
    dcg = -1.0im * (1.0 .+ 0.1*sin.(t)) * ψ[2];
    dce = -1.0im * (1.0 .+ 0.1*sin.(t)) * ψ[1];
    [dcg, dce]
end;

function rabi_problem1!(dψ, ψ, params, t)
    dψ[1] = -1.0im * (1.0 .+ 0.1*sin.(t)) * ψ[2];
    dψ[2] = -1.0im * (1.0 .+ 0.1*sin.(t)) * ψ[1];
    nothing
end;

function rabi_problem2(ψ, params, t)
    dcg = -1.0im * (ComplexF64(1.0) .+ ComplexF64(0.1)*sin.(t)) * ψ[2]
    dce = -1.0im * (ComplexF64(1.0) .+ ComplexF64(0.1)*sin.(t)) * ψ[1]
    SA[dcg, dce]
end;

Results for implementation in DifferentialEquations.jl:

tspan = (0.0, 100.0);
initial = [ComplexF64(1.0), ComplexF64(0.0)];
initial_static = SA[ComplexF64(1.0), ComplexF64(0.0)];

problem1 = ODEProblem(rabi_problem1, initial, tspan);

@btime solve(problem1, Tsit5());

452.545 μs (9822 allocations: 927.69 KiB)

problem2 = ODEProblem(rabi_problem1, initial, tspan; save_everystep = false)

@btime solve(problem2, Tsit5());

412.397 μs (8677 allocations: 813.20 KiB)

problem3 = ODEProblem(rabi_problem1!, initial, tspan; save_everystep = false)

@btime solve(problem3);

81.041 μs (153 allocations: 12.94 KiB)

problem4 = ODEProblem(rabi_problem2, initial_static, tspan)

@btime solve(problem4, Tsit5());

51.547 μs (152 allocations: 48.84 KiB)

problem5 = ODEProblem(rabi_problem2, initial_static, tspan; save_everystep = false)

@btime solve(problem5, Tsit5());

39.312 μs (22 allocations: 2.36 KiB)

I run my code in Jupyter(IJulia v1.24.2) and use QuantumOptics v1.0.14. I decided to do the benchmarking because simulations for 5-level system in my current project allocate up to 8GiB(

[Krastanov]

I suspect your use of TimeDependentSum is incorrect. Currently in H2 you create a new TimeDependentSum object (which allocates new operators, etc). Rather, you want to create a single TimeDependentSum object and just update the scalar weight in it.

Check out https://docs.qojulia.org/timeevolution/timedependent-problems/

In particular, it seems the constructor should be

H2 = TimeDependentSum([t -> 1.0 .+ 0.1*sin.(t), t -> 1.0 .+ 0.1*sin.(t)], [σge, σeg])

Notice that H2 is NOT a function anymore, so it will not be creating a new TimeDependentSum instance. Rather the solver will be appropriately updating (in-place modifying) the single instance you give it.

[mgoloshchapov]

@Krastanov, thank you for your answer. I tried examples from https://docs.qojulia.org/timeevolution/timedependent-problems/ and different constructor for H2 that you suggested. They work slightly better than my initial versions, but still far from realisations with DifferentialEquations.jl(

H5 = TimeDependentSum([t -> 1.0 .+ 0.1*sin.(t), t -> 1.0 .+ 0.1*sin.(t)], [σge, σeg]);

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H5; alg=Tsit5());

2.544 ms (20105 allocations: 2.02 MiB)

const H6 = LazySum(ComplexF64[0.0, 0.0],[σge, σeg]);
function H_pump(t, psi)
  H6.factors[1] = 1.0 + 0.1*sin(t);
  H6.factors[2] = 1.0 + 0.1*sin(t);
  return H6
end;

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H_pump; alg=Tsit5());

2.714 ms (26511 allocations: 2.16 MiB)

const coeff_funcs = [t->1.0 + 0.1*sin(t),t->1.0 + 0.1*sin(t)];
const H7 = LazySum([c(tspan[1]) for c∈coeff_funcs],[σge,σeg])

# Dynamic function
function Ht(t,psi)
    for i=1:length(H7.factors)
        H7.factors[i] = coeff_funcs[i](t)
    end
    return H7
end

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, Ht; alg=Tsit5());

2.937 ms (32929 allocations: 2.21 MiB)

I also tried example from tutorial with larger tspan:

# Generic Gaussian pulse
pulse(t,t0,Ω) = @. Ω*exp(-(t-t0)^2)

# Operators
b1 = SpinBasis(1//2)
sx1 = tensor(sigmax(b1), one(b1))
sx2 = tensor(one(b1), sigmax(b1))

# Define coefficients and Hamiltonian
tspan = [0.0:0.01:100.0;]
const coeff_funcs = [t->pulse(t,1,0.5),t->(pulse(t,5,1))]
const H = LazySum([c(tspan[1]) for c∈coeff_funcs],[sx1,sx2])

# Dynamic function
function Ht(t,psi)
    for i=1:length(H.factors)
        H.factors[i] = coeff_funcs[i](t)
    end
    return H
end

psi0 = tensor(spindown(b1), spindown(b1));

@btime timeevolution.schroedinger_dynamic(tspan, psi0, Ht);

1.913 ms (21072 allocations: 2.80 MiB)

Do you have the same performance on your computer?

[Krastanov]

Indeed, I confirm that I see the same large number of allocations on QuantumOptics v1.0.14 and julia 1.9.3

@amilsted , I think you are most familiar with this portion of the code. Any ideas?

[amilsted]

There are a couple of (related) things going on:

1. The QO evolution functions will store the state at every element of tspan, so allocations are needed for that.
2. The integrator will stop at every element of tspan, so has to do at least that many timesteps.
3. You are using a two-element tspan for DiffEq and a 10000-element one for QO.

julia> tspan = [0.0:0.01:100.0;];

julia> @time timeevolution.schroedinger_dynamic(tspan, ψ0, H5; alg=Tsit5());
  0.008411 seconds (20.11 k allocations: 2.017 MiB)

julia> fout = (x...)->nothing
#23 (generic function with 1 method)

julia> @time timeevolution.schroedinger_dynamic(tspan, ψ0, H5; alg=Tsit5(), fout=fout);
  0.001917 seconds (111 allocations: 489.125 KiB)

julia> tspan = (0.0, 100.0);

julia> @time timeevolution.schroedinger_dynamic(tspan, ψ0, H5; alg=Tsit5(), fout=fout);
  0.000661 seconds (100 allocations: 6.508 KiB)

Now we're much closer to DiffEq. The other thing is that you're effectively hardcoding a sparse representation of the operators in the DiffEq case, this probably accounts for the rest.

Btw, there's no need for . broadcasting syntax in defining the time-dependent operator. This H5_ = TimeDependentSum([t -> 1.0 + 0.1*sin(t), t -> 1.0 + 0.1*sin(t)], (σge, σeg)) is marginally faster.

[amilsted]

By the way, you should find you can use static arrays in the QO case too. Operator(my_basis, some_static_array) and Ket(my_basis, static_vector) should work.

[mgoloshchapov]

Thank you so much! I changed tspan and operator definition both in example and my project, everything works fast now.

It seems like static arrays don't help here:

Without SA:

tspan = (0.0, 100.0);
ψ0 = nlevelstate(basis, 1);

H5 = TimeDependentSum([t -> 1.0 + 0.1*sin(t), t -> 1.0 + 0.1*sin(t)], [σge, σeg]);

@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H5; alg=Tsit5());

641.848 μs (96 allocations: 6.53 KiB)

With SA:

σge_static = Operator(basis, 
                SA[0.0+0.0im  1.0+0.0im
                 0.0+0.0im  0.0+0.0im]);
σeg_static = Operator(basis, 
                SA[0.0+0.0im  0.0+0.0im
                 1.0+0.0im  0.0+0.0im]);

ψ0_static = Ket(basis, SA[1.0+0.0im,0.0+0.0im]);

H8 = TimeDependentSum([t -> 1.0 + 0.1*sin(t), t -> 1.0 + 0.1*sin(t)], [σge_static, σeg_static]);
@btime timeevolution.schroedinger_dynamic(tspan, ψ0, H8; alg=Tsit5());

611.492 μs (12924 allocations: 758.17 KiB)

One thing I don't really understand is why there is an error when I pass ψ0_static with H8:

@btime timeevolution.schroedinger_dynamic(tspan, ψ0_static, H8; alg=Tsit5());

Initial condition incompatible with functional form.
Detected an in-place function with an initial condition of type Number or SArray.
This is incompatible because Numbers cannot be mutated, i.e.
`x = 2.0; y = 2.0; x .= y` will error.

If using a immutable initial condition type, please use the out-of-place form.
I.e. define the function `du=f(u,p,t)` instead of attempting to "mutate" the immutable `du`.

If your differential equation function was defined with multiple dispatches and one is
in-place, then the automatic detection will choose in-place. In this case, override the
choice in the problem constructor, i.e. `ODEProblem{false}(f,u0,tspan,p,kwargs...)`.

For a longer discussion on mutability vs immutability and in-place vs out-of-place, see:
https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Stacktrace:
  [1] get_concrete_u0(prob::ODEProblem{SVector{2, ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, QuantumOptics.timeevolution.var"#df_#3"{QuantumOptics.timeevolution.var"#dschroedinger_#52"{QuantumOptics.timeevolution.var"#_tdop_schroedinger_wrapper#9"{TimeDependentSum{NLevelBasis{Int64}, NLevelBasis{Int64}, Tuple{var"#15#17", var"#16#18"}, LazySum{NLevelBasis{Int64}, NLevelBasis{Int64}, Vector{Float64}, Tuple{Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}, Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}}}, Float64}}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool, t0::Float64, kwargs::Base.Pairs{Symbol, Any, NTuple{8, Symbol}, NamedTuple{(:u0, :p, :reltol, :abstol, :save_everystep, :save_start, :save_end, :callback), Tuple{SVector{2, ComplexF64}, SciMLBase.NullParameters, Float64, Float64, Bool, Bool, Bool, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#30#31", DiffEqCallbacks.SavingAffect{QuantumOptics.timeevolution.var"#fout_#4"{Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, QuantumOptics.timeevolution.var"#fout#7"}, Float64, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, DataStructures.BinaryMinHeap{Float64}, Vector{Float64}}, typeof(DiffEqCallbacks.saving_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/MFgVe/src/solve.jl:1237
  [2] get_concrete_problem(prob::ODEProblem{SVector{2, ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, QuantumOptics.timeevolution.var"#df_#3"{QuantumOptics.timeevolution.var"#dschroedinger_#52"{QuantumOptics.timeevolution.var"#_tdop_schroedinger_wrapper#9"{TimeDependentSum{NLevelBasis{Int64}, NLevelBasis{Int64}, Tuple{var"#15#17", var"#16#18"}, LazySum{NLevelBasis{Int64}, NLevelBasis{Int64}, Vector{Float64}, Tuple{Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}, Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}}}, Float64}}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, isadapt::Bool; kwargs::Base.Pairs{Symbol, Any, NTuple{8, Symbol}, NamedTuple{(:u0, :p, :reltol, :abstol, :save_everystep, :save_start, :save_end, :callback), Tuple{SVector{2, ComplexF64}, SciMLBase.NullParameters, Float64, Float64, Bool, Bool, Bool, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#30#31", DiffEqCallbacks.SavingAffect{QuantumOptics.timeevolution.var"#fout_#4"{Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, QuantumOptics.timeevolution.var"#fout#7"}, Float64, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, DataStructures.BinaryMinHeap{Float64}, Vector{Float64}}, typeof(DiffEqCallbacks.saving_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/MFgVe/src/solve.jl:1093
  [3] solve_up(prob::ODEProblem{SVector{2, ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, QuantumOptics.timeevolution.var"#df_#3"{QuantumOptics.timeevolution.var"#dschroedinger_#52"{QuantumOptics.timeevolution.var"#_tdop_schroedinger_wrapper#9"{TimeDependentSum{NLevelBasis{Int64}, NLevelBasis{Int64}, Tuple{var"#15#17", var"#16#18"}, LazySum{NLevelBasis{Int64}, NLevelBasis{Int64}, Vector{Float64}, Tuple{Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}, Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}}}, Float64}}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::SVector{2, ComplexF64}, p::SciMLBase.NullParameters, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:reltol, :abstol, :save_everystep, :save_start, :save_end, :callback), Tuple{Float64, Float64, Bool, Bool, Bool, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#30#31", DiffEqCallbacks.SavingAffect{QuantumOptics.timeevolution.var"#fout_#4"{Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, QuantumOptics.timeevolution.var"#fout#7"}, Float64, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, DataStructures.BinaryMinHeap{Float64}, Vector{Float64}}, typeof(DiffEqCallbacks.saving_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/MFgVe/src/solve.jl:1000
  [4] solve(prob::ODEProblem{SVector{2, ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, QuantumOptics.timeevolution.var"#df_#3"{QuantumOptics.timeevolution.var"#dschroedinger_#52"{QuantumOptics.timeevolution.var"#_tdop_schroedinger_wrapper#9"{TimeDependentSum{NLevelBasis{Int64}, NLevelBasis{Int64}, Tuple{var"#15#17", var"#16#18"}, LazySum{NLevelBasis{Int64}, NLevelBasis{Int64}, Vector{Float64}, Tuple{Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}, Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}}}, Float64}}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Any, NTuple{6, Symbol}, NamedTuple{(:reltol, :abstol, :save_everystep, :save_start, :save_end, :callback), Tuple{Float64, Float64, Bool, Bool, Bool, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#30#31", DiffEqCallbacks.SavingAffect{QuantumOptics.timeevolution.var"#fout_#4"{Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, QuantumOptics.timeevolution.var"#fout#7"}, Float64, Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, DataStructures.BinaryMinHeap{Float64}, Vector{Float64}}, typeof(DiffEqCallbacks.saving_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}}}}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/MFgVe/src/solve.jl:929
  [5] integrate(tspan::Tuple{Float64, Float64}, df::QuantumOptics.timeevolution.var"#dschroedinger_#52"{QuantumOptics.timeevolution.var"#_tdop_schroedinger_wrapper#9"{TimeDependentSum{NLevelBasis{Int64}, NLevelBasis{Int64}, Tuple{var"#15#17", var"#16#18"}, LazySum{NLevelBasis{Int64}, NLevelBasis{Int64}, Vector{Float64}, Tuple{Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}, Operator{NLevelBasis{Int64}, NLevelBasis{Int64}, SMatrix{2, 2, ComplexF64, 4}}}}, Float64}}}, x0::SVector{2, ComplexF64}, state::Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, dstate::Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, fout::QuantumOptics.timeevolution.var"#fout#7"; alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, steady_state::Bool, tol::Float64, save_everystep::Bool, saveat::Tuple{Float64, Float64}, callback::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ QuantumOptics.timeevolution ~/.julia/packages/QuantumOptics/6utec/src/timeevolution_base.jl:59
  [6] #integrate#6
    @ ~/.julia/packages/QuantumOptics/6utec/src/timeevolution_base.jl:75 [inlined]
  [7] schroedinger_dynamic(tspan::Tuple{Float64, Float64}, psi0::Ket{NLevelBasis{Int64}, SVector{2, ComplexF64}}, f::Function; fout::Nothing, kwargs::Base.Pairs{Symbol, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, Tuple{Symbol}, NamedTuple{(:alg,), Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}}})
    @ QuantumOptics.timeevolution ~/.julia/packages/QuantumOptics/6utec/src/schroedinger.jl:54
  [8] schroedinger_dynamic
    @ ~/.julia/packages/QuantumOptics/6utec/src/schroedinger.jl:46 [inlined]
  [9] #schroedinger_dynamic#53
    @ ~/.julia/packages/QuantumOptics/6utec/src/schroedinger.jl:59 [inlined]
 [10] var"##core#1153"()
    @ Main ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:489
 [11] var"##sample#1154"(::Tuple{}, __params::BenchmarkTools.Parameters)
    @ Main ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:495
 [12] _run(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; verbose::Bool, pad::String, kwargs::Base.Pairs{Symbol, Integer, NTuple{4, Symbol}, NamedTuple{(:samples, :evals, :gctrial, :gcsample), Tuple{Int64, Int64, Bool, Bool}}})
    @ BenchmarkTools ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:99
 [13] #invokelatest#2
    @ ./essentials.jl:821 [inlined]
 [14] invokelatest
    @ ./essentials.jl:816 [inlined]
 [15] #run_result#45
    @ ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:34 [inlined]
 [16] run_result
    @ ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:34 [inlined]
 [17] run(b::BenchmarkTools.Benchmark, p::BenchmarkTools.Parameters; progressid::Nothing, nleaves::Float64, ndone::Float64, kwargs::Base.Pairs{Symbol, Integer, NTuple{5, Symbol}, NamedTuple{(:verbose, :samples, :evals, :gctrial, :gcsample), Tuple{Bool, Int64, Int64, Bool, Bool}}})
    @ BenchmarkTools ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:117
 [18] run (repeats 2 times)
    @ ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:117 [inlined]
 [19] #warmup#54
    @ ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:169 [inlined]
 [20] warmup(item::BenchmarkTools.Benchmark)
    @ BenchmarkTools ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:168
 [21] top-level scope
    @ ~/.julia/packages/BenchmarkTools/0owsb/src/execution.jl:575

[amilsted]

Ah, yeah - that's a point. We try to do in-place updates of the statevector and that won't work with static arrays.