diff --git a/CHANGELOG.md b/CHANGELOG.md index 9273811f9..8eb59f8d1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Fix time evolution output when using `saveat` keyword argument. ([#398]) - Align some attributes of `mcsolve`, `ssesolve` and `smesolve` results with `QuTiP`. ([#402]) - Improve ensemble generation of `ssesolve` and change parameters handling on stochastic processes. ([#403]) +- Set default trajectories to 500 and rename the keyword argument `ensemble_method` to `ensemblealg`. ([#405]) ## [v0.26.0] Release date: 2025-02-09 @@ -132,3 +133,4 @@ Release date: 2024-11-13 [#398]: https://github.com/qutip/QuantumToolbox.jl/issues/398 [#402]: https://github.com/qutip/QuantumToolbox.jl/issues/402 [#403]: https://github.com/qutip/QuantumToolbox.jl/issues/403 +[#405]: https://github.com/qutip/QuantumToolbox.jl/issues/405 diff --git a/benchmarks/timeevolution.jl b/benchmarks/timeevolution.jl index b76816ed4..6c1d18f25 100644 --- a/benchmarks/timeevolution.jl +++ b/benchmarks/timeevolution.jl @@ -50,7 +50,7 @@ function benchmark_timeevolution!(SUITE) ntraj = 100, e_ops = $e_ops, progress_bar = Val(false), - ensemble_method = EnsembleSerial(), + ensemblealg = EnsembleSerial(), ) SUITE["Time Evolution"]["time-independent"]["mcsolve"]["Multithreaded"] = @benchmarkable mcsolve( $H, @@ -60,7 +60,7 @@ function benchmark_timeevolution!(SUITE) ntraj = 100, e_ops = $e_ops, progress_bar = Val(false), - ensemble_method = EnsembleThreads(), + ensemblealg = EnsembleThreads(), ) return nothing diff --git a/docs/make.jl b/docs/make.jl index 0aec9e286..e624cbee3 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,7 +34,7 @@ const PAGES = [ "Home" => "index.md", "Getting Started" => [ "Brief Example" => "getting_started/brief_example.md", - "Key differences from QuTiP" => "getting_started/qutip_differences.md", + # "Key differences from QuTiP" => "getting_started/qutip_differences.md", "The Importance of Type-Stability" => "getting_started/type_stability.md", "Example: Create QuantumToolbox.jl Logo" => "getting_started/logo.md", # "Cite QuantumToolbox.jl" => "getting_started/cite.md", @@ -51,7 +51,7 @@ const PAGES = [ "Time Evolution Solutions" => "users_guide/time_evolution/solution.md", "Schrödinger Equation Solver" => "users_guide/time_evolution/sesolve.md", "Lindblad Master Equation Solver" => "users_guide/time_evolution/mesolve.md", - "Monte-Carlo Solver" => "users_guide/time_evolution/mcsolve.md", + "Monte Carlo Solver" => "users_guide/time_evolution/mcsolve.md", "Stochastic Solver" => "users_guide/time_evolution/stochastic.md", "Solving Problems with Time-dependent Hamiltonians" => "users_guide/time_evolution/time_dependent.md", ], diff --git a/docs/src/users_guide/cluster.md b/docs/src/users_guide/cluster.md index f290ae1b4..d36259f03 100644 --- a/docs/src/users_guide/cluster.md +++ b/docs/src/users_guide/cluster.md @@ -112,7 +112,7 @@ e_ops = [Sx, Sy, Sz] tlist = range(0, 10, 100) -sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=5000, ensemble_method=EnsembleSplitThreads()) +sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops=e_ops, ntraj=5000, ensemblealg=EnsembleSplitThreads()) ## @@ -131,7 +131,7 @@ println("Hello! You have $(nworkers()) workers with $(remotecall_fetch(Threads.n command, we test that the distributed network is correctly initialized. The `remotecall_fetch(Threads.nthreads, 2)` command returns the number of threads of the worker with ID `2`. -We then write the main part of the script, where we define the lattice through the [`Lattice`](@ref) function. We set the parameters and define the Hamiltonian and collapse operators with the [`DissipativeIsing`](@ref) function. We also define the expectation operators `e_ops` and the initial state `ψ0`. Finally, we perform the Monte Carlo quantum trajectories with the [`mcsolve`](@ref) function. The `ensemble_method=EnsembleSplitThreads()` argument is used to parallelize the Monte Carlo quantum trajectories, by splitting the ensemble of trajectories among the workers. For a more detailed explanation of the different ensemble methods, you can check the [official documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) of the [**DifferentialEquations.jl**](https://github.com/SciML/DifferentialEquations.jl/) package. Finally, the `rmprocs(workers())` command is used to remove the workers after the computation is finished. +We then write the main part of the script, where we define the lattice through the [`Lattice`](@ref) function. We set the parameters and define the Hamiltonian and collapse operators with the [`DissipativeIsing`](@ref) function. We also define the expectation operators `e_ops` and the initial state `ψ0`. Finally, we perform the Monte Carlo quantum trajectories with the [`mcsolve`](@ref) function. The `ensemblealg=EnsembleSplitThreads()` argument is used to parallelize the Monte Carlo quantum trajectories, by splitting the ensemble of trajectories among the workers. For a more detailed explanation of the different ensemble methods, you can check the [official documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) of the [**DifferentialEquations.jl**](https://github.com/SciML/DifferentialEquations.jl/) package. Finally, the `rmprocs(workers())` command is used to remove the workers after the computation is finished. The output of the script will be printed in the `output.out` file, which contains an output similar to the following: diff --git a/docs/src/users_guide/time_evolution/intro.md b/docs/src/users_guide/time_evolution/intro.md index b8ad2ca0f..e1755183d 100644 --- a/docs/src/users_guide/time_evolution/intro.md +++ b/docs/src/users_guide/time_evolution/intro.md @@ -17,6 +17,9 @@ - [Example: Harmonic oscillator in thermal bath](@ref doc-TE:Example:Harmonic-oscillator-in-thermal-bath) - [Example: Two-level atom coupled to dissipative single-mode cavity](@ref doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity) - [Monte Carlo Solver](@ref doc-TE:Monte-Carlo-Solver) + - [Monte Carlo wave-function](@ref doc-TE:Monte-Carlo-wave-function) + - [Example: Two-level atom coupled to dissipative single-mode cavity (MC)](@ref doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity-(MC)) + - [Running trajectories in parallel](@ref doc-TE:Running-trajectories-in-parallel) - [Stochastic Solver](@ref doc-TE:Stochastic-Solver) - [Stochastic Schrödinger equation](@ref doc-TE:Stochastic-Schrödinger-equation) - [Stochastic master equation](@ref doc-TE:Stochastic-master-equation) diff --git a/docs/src/users_guide/time_evolution/mcsolve.md b/docs/src/users_guide/time_evolution/mcsolve.md index d62540830..e5a3e0309 100644 --- a/docs/src/users_guide/time_evolution/mcsolve.md +++ b/docs/src/users_guide/time_evolution/mcsolve.md @@ -1,3 +1,134 @@ # [Monte Carlo Solver](@id doc-TE:Monte-Carlo-Solver) -This page is still under construction, please visit [API](@ref doc-API) first. \ No newline at end of file +## [Monte Carlo wave-function](@id doc-TE:Monte-Carlo-wave-function) + +Where as the density matrix formalism describes the ensemble average over many identical realizations of a quantum system, the Monte Carlo (MC), or quantum-jump approach to wave function evolution, allows for simulating an individual realization of the system dynamics. Here, the environment is continuously monitored, resulting in a series of quantum jumps in the system wave function, conditioned on the increase in information gained about the state of the system via the environmental measurements. In general, this evolution is governed by the Schrödinger equation with a non-Hermitian effective Hamiltonian + +```math +\hat{H}_{\textrm{eff}} = \hat{H} - \frac{i}{2} \sum_{n=1}^N \hat{C}_n^\dagger \hat{C}_n. +``` + +where ``\hat{H}`` is the system Hamiltonian and ``\hat{C}_n`` are collapse (jump) operators (assume ``N`` is the total number of collapse operators). Each collapse operator corresponds to a separate irreversible process with rate ``\gamma_n``. Here, the strictly negative non-Hermitian portion of the above equation gives rise to a reduction in the norm of the wave function, that to first-order in a small time ``\delta t``, is given by + +```math +\langle \psi(t + \delta t) | \psi(t + \delta t) \rangle = 1 - \delta p, +``` + +where + +```math +\delta p = \delta t \sum_{n=1}^N \langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle, +``` + +and ``\delta t`` is such that ``\delta p \ll 1``. With a probability of remaining in the state ``| \psi(t + \delta t) \rangle`` given by ``1 - \delta p``, the corresponding quantum jump probability is thus ``\delta p``. If the environmental measurements register a quantum jump, say via the emission of a photon into the environment, or a change in the spin of a quantum dot, the wave function undergoes a jump into a state defined by projecting ``| \psi(t) \rangle`` using the collapse operator ``\hat{C}_n`` corresponding to the measurement + +```math +| \psi(t+\delta t) \rangle = \frac{\hat{C}_n |\psi(t)\rangle}{ \sqrt{\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle} }. +``` + +If more than a single collapse operator is present in ``\hat{H}_{\textrm{eff}}``, the probability of collapse due to the ``n``-th operator ``\hat{C}_n`` is given by + +```math +P_n(t) = \frac{1}{\delta p}\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle. +``` + +Note that the probability of all collapses should be normalized to unity for all time ``t``, namely + +```math +\sum_{n=1}^N P_n(t) = 1 ~~~\forall~~t. +``` + +Evaluating the MC evolution to first-order in time is quite tedious. Instead, `QuantumToolbox.jl` provides the function [`mcsolve`](@ref) which uses the following algorithm to simulate a single realization of a quantum system. Starting from a pure state ``| \psi(0) \rangle``: + +1. Choose two random numbers (``r_1`` and ``r_2``) between 0 and 1, where ``r_1`` represents the probability that a quantum jump occurs and ``r_2`` is used to select which collapse operator was responsible for the jump. +1. Integrate the Schrödinger equation with respect to the effective Hamiltonian ``\hat{H}_{\textrm{eff}}`` until a time ``\tau`` such that the norm of the wave function satisfies ``\langle \psi(\tau) | \psi(\tau) \rangle = r_1``, at which point a jump occurs +1. The resultant jump projects the system at time ``\tau`` into one of the renormalized states ``| \psi(\tau + \delta t) \rangle``. The corresponding collapse operator ``\hat{C}_n`` is chosen such that ``\tilde{n} \leq N`` is the smallest integer satisfying ``\sum_{n=1}^{\tilde{n}} P_n(\tau) \geq r_2``. +1. Using the renormalized state from previous step as the new initial condition at time ``\tau`` and repeat the above procedure until the final simulation time is reached. + +## [Example: Two-level atom coupled to dissipative single-mode cavity (MC)](@id doc-TE:Example:Two-level-atom-coupled-to-dissipative-single-mode-cavity-(MC)) + +In `QuantumToolbox.jl`, Monte Carlo evolution is implemented with the [`mcsolve`](@ref) function. It takes nearly the same arguments as the [`mesolve`](@ref) function for [Lindblad master equation evolution](@ref doc-TE:Lindblad-Master-Equation-Solver), except that the initial state must be a [`Ket`](@ref) vector, as oppose to a density matrix, and there is an optional keyword argument `ntraj` that defines the number of stochastic trajectories to be simulated. By default, `ntraj=500` indicating that `500` Monte Carlo trajectories will be performed. + +To illustrate the use of the Monte Carlo evolution of quantum systems in `QuantumToolbox.jl`, let’s again consider the case of a two-level atom coupled to a leaky cavity. The only differences to the master equation treatment is that in this case we invoke the [`mcsolve`](@ref) function instead of [`mesolve`](@ref) + +```@setup mcsolve +using QuantumToolbox + +using CairoMakie +CairoMakie.enable_only_mime!(MIME"image/svg+xml"()) +``` + +```@example mcsolve +times = LinRange(0.0, 10.0, 200) + +ψ0 = tensor(fock(2, 0), fock(10, 8)) +a = tensor(qeye(2), destroy(10)) +σm = tensor(destroy(2), qeye(10)) +H = 2 * π * a' * a + 2 * π * σm' * σm + 2 * π * 0.25 * (σm * a' + σm' * a) + +c_ops = [sqrt(0.1) * a] +e_ops = [a' * a, σm' * σm] + +sol_500 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops) + +# plot by CairoMakie.jl +fig = Figure(size = (500, 350)) +ax = Axis(fig[1, 1], + xlabel = "Time", + ylabel = "Expectation values", + title = "Monte Carlo time evolution (500 trajectories)", +) +lines!(ax, times, real(sol_500.expect[1,:]), label = "cavity photon number", linestyle = :solid) +lines!(ax, times, real(sol_500.expect[2,:]), label = "atom excitation probability", linestyle = :dash) + +axislegend(ax, position = :rt) + +fig +``` + +The advantage of the Monte Carlo method over the master equation approach is that only the state vector ([`Ket`](@ref)) is required to be kept in the computers memory, as opposed to the entire density matrix. For large quantum system this becomes a significant advantage, and the Monte Carlo solver is therefore generally recommended for such systems. However, for small systems, the added overhead of averaging a large number of stochastic trajectories to obtain the open system dynamics, as well as starting the multiprocessing functionality, outweighs the benefit of the minor (in this case) memory saving. Master equation methods are therefore generally more efficient when Hilbert space sizes are on the order of a couple of hundred states or smaller. + +We can also change the number of trajectories (`ntraj`). This can be used to explore the convergence of the Monte Carlo solver. For example, the following code plots the expectation values for `1`, `10` and `100` trajectories: + +```@example mcsolve +e_ops = [a' * a] + +sol_1 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 1) +sol_10 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 10) +sol_100 = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ntraj = 100) + +# plot by CairoMakie.jl +fig = Figure(size = (500, 350)) +ax = Axis(fig[1, 1], + xlabel = "Time", + ylabel = "Expectation values", + title = "Monte Carlo time evolution", +) +lines!(ax, times, real(sol_1.expect[1,:]), label = "1 trajectory", linestyle = :dashdot) +lines!(ax, times, real(sol_10.expect[1,:]), label = "10 trajectories", linestyle = :dash) +lines!(ax, times, real(sol_100.expect[1,:]), label = "100 trajectories", linestyle = :solid) + +axislegend(ax, position = :rt) + +fig +``` + +## [Running trajectories in parallel](@id doc-TE:Running-trajectories-in-parallel) + +Monte Carlo evolutions often need hundreds of trajectories to obtain sufficient statistics. Since all trajectories are independent of each other, they can be computed in parallel. The keyword argument `ensemblealg` can specify how the multiple trajectories are handled. The common ensemble methods are: + +- `EnsembleSerial()`: No parallelism +- `EnsembleThreads()`: **The default.** This uses multithreading. +- `EnsembleDistributed()`: This uses as many processors as you have Julia processes. +- `EnsembleSplitThreads()`: This uses multithreading on each process. + +!!! note "Other Ensemble Algorithms" + See the [documentation of `DifferentialEquations.jl`](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/) for more details. Also, see Julia's documentation for more details about multithreading and adding more processes. + +```julia +sol_serial = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleSerial()) +sol_parallel = mcsolve(H, ψ0, times, c_ops, e_ops=e_ops, ensemblealg=EnsembleThreads()); +``` + +!!! tip "Parallelization on a Cluster" + See the section [Intensive parallelization on a Cluster](@ref doc:Intensive-parallelization-on-a-Cluster) for more details. diff --git a/src/QuantumToolbox.jl b/src/QuantumToolbox.jl index 4776c0d43..7f512d685 100644 --- a/src/QuantumToolbox.jl +++ b/src/QuantumToolbox.jl @@ -28,6 +28,7 @@ import SciMLBase: ODEProblem, SDEProblem, EnsembleProblem, + EnsembleAlgorithm, EnsembleSerial, EnsembleThreads, EnsembleSplitThreads, diff --git a/src/time_evolution/mcsolve.jl b/src/time_evolution/mcsolve.jl index 066fa68c0..172e8e099 100644 --- a/src/time_evolution/mcsolve.jl +++ b/src/time_evolution/mcsolve.jl @@ -148,8 +148,8 @@ end e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), jump_callback::TJC = ContinuousLindbladJumpCallback(), progress_bar::Union{Val,Bool} = Val(true), prob_func::Union{Function, Nothing} = nothing, @@ -201,7 +201,7 @@ If the environmental measurements register a quantum jump, the wave function und - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used. - `rng`: Random number generator for reproducibility. - `ntraj`: Number of trajectories to use. -- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`. +- `ensemblealg`: Ensemble algorithm to use. Default to `EnsembleThreads()`. - `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise. - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities. - `prob_func`: Function to use for generating the ODEProblem. @@ -227,8 +227,8 @@ function mcsolveEnsembleProblem( e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), jump_callback::TJC = ContinuousLindbladJumpCallback(), progress_bar::Union{Val,Bool} = Val(true), prob_func::Union{Function,Nothing} = nothing, @@ -238,7 +238,7 @@ function mcsolveEnsembleProblem( _prob_func = isnothing(prob_func) ? _ensemble_dispatch_prob_func(rng, ntraj, tlist, _mcsolve_prob_func) : prob_func _output_func = output_func isa Nothing ? - _ensemble_dispatch_output_func(ensemble_method, progress_bar, ntraj, _mcsolve_output_func) : output_func + _ensemble_dispatch_output_func(ensemblealg, progress_bar, ntraj, _mcsolve_output_func) : output_func prob_mc = mcsolveProblem( H, @@ -272,8 +272,8 @@ end e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), jump_callback::TJC = ContinuousLindbladJumpCallback(), progress_bar::Union{Val,Bool} = Val(true), prob_func::Union{Function, Nothing} = nothing, @@ -327,7 +327,7 @@ If the environmental measurements register a quantum jump, the wave function und - `params`: Parameters to pass to the solver. This argument is usually expressed as a `NamedTuple` or `AbstractVector` of parameters. For more advanced usage, any custom struct can be used. - `rng`: Random number generator for reproducibility. - `ntraj`: Number of trajectories to use. -- `ensemble_method`: Ensemble method to use. Default to `EnsembleThreads()`. +- `ensemblealg`: Ensemble algorithm to use. Default to `EnsembleThreads()`. - `jump_callback`: The Jump Callback type: Discrete or Continuous. The default is `ContinuousLindbladJumpCallback()`, which is more precise. - `progress_bar`: Whether to show the progress bar. Using non-`Val` types might lead to type instabilities. - `prob_func`: Function to use for generating the ODEProblem. @@ -337,7 +337,7 @@ If the environmental measurements register a quantum jump, the wave function und # Notes -- `ensemble_method` can be one of `EnsembleThreads()`, `EnsembleSerial()`, `EnsembleDistributed()` +- `ensemblealg` can be one of `EnsembleThreads()`, `EnsembleSerial()`, `EnsembleDistributed()` - The states will be saved depend on the keyword argument `saveat` in `kwargs`. - If `e_ops` is empty, the default value of `saveat=tlist` (saving the states corresponding to `tlist`), otherwise, `saveat=[tlist[end]]` (only save the final state). You can also specify `e_ops` and `saveat` separately. - The default tolerances in `kwargs` are given as `reltol=1e-6` and `abstol=1e-8`. @@ -357,8 +357,8 @@ function mcsolve( e_ops::Union{Nothing,AbstractVector,Tuple} = nothing, params = NullParameters(), rng::AbstractRNG = default_rng(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), jump_callback::TJC = ContinuousLindbladJumpCallback(), progress_bar::Union{Val,Bool} = Val(true), prob_func::Union{Function,Nothing} = nothing, @@ -376,7 +376,7 @@ function mcsolve( params = params, rng = rng, ntraj = ntraj, - ensemble_method = ensemble_method, + ensemblealg = ensemblealg, jump_callback = jump_callback, progress_bar = progress_bar, prob_func = prob_func, @@ -384,17 +384,17 @@ function mcsolve( kwargs..., ) - return mcsolve(ens_prob_mc, alg, ntraj, ensemble_method, normalize_states) + return mcsolve(ens_prob_mc, alg, ntraj, ensemblealg, normalize_states) end function mcsolve( ens_prob_mc::TimeEvolutionProblem, alg::OrdinaryDiffEqAlgorithm = Tsit5(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), normalize_states = Val(true), ) - sol = _ensemble_dispatch_solve(ens_prob_mc, alg, ensemble_method, ntraj) + sol = _ensemble_dispatch_solve(ens_prob_mc, alg, ensemblealg, ntraj) dims = ens_prob_mc.dimensions _sol_1 = sol[:, 1] diff --git a/src/time_evolution/time_evolution.jl b/src/time_evolution/time_evolution.jl index 672414928..c32805524 100644 --- a/src/time_evolution/time_evolution.jl +++ b/src/time_evolution/time_evolution.jl @@ -308,7 +308,7 @@ end function _ensemble_dispatch_solve( ens_prob_mc::TimeEvolutionProblem, alg::Union{<:OrdinaryDiffEqAlgorithm,<:StochasticDiffEqAlgorithm}, - ensemble_method::ET, + ensemblealg::ET, ntraj::Int, ) where {ET<:Union{EnsembleSplitThreads,EnsembleDistributed}} sol = nothing @@ -319,7 +319,7 @@ function _ensemble_dispatch_solve( end @async begin - sol = solve(ens_prob_mc.prob, alg, ensemble_method, trajectories = ntraj) + sol = solve(ens_prob_mc.prob, alg, ensemblealg, trajectories = ntraj) put!(ens_prob_mc.kwargs.channel, false) end end @@ -329,10 +329,10 @@ end function _ensemble_dispatch_solve( ens_prob_mc::TimeEvolutionProblem, alg::Union{<:OrdinaryDiffEqAlgorithm,<:StochasticDiffEqAlgorithm}, - ensemble_method, + ensemblealg, ntraj::Int, ) - sol = solve(ens_prob_mc.prob, alg, ensemble_method, trajectories = ntraj) + sol = solve(ens_prob_mc.prob, alg, ensemblealg, trajectories = ntraj) return sol end diff --git a/src/time_evolution/time_evolution_dynamical.jl b/src/time_evolution/time_evolution_dynamical.jl index b4ec427a5..77969751e 100644 --- a/src/time_evolution/time_evolution_dynamical.jl +++ b/src/time_evolution/time_evolution_dynamical.jl @@ -607,8 +607,8 @@ function dsf_mcsolveEnsembleProblem( dsf_params::NamedTuple = NamedTuple(); e_ops::Function = (op_list, p) -> (), params::NamedTuple = NamedTuple(), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), δα_list::Vector{<:Real} = fill(0.2, length(op_list)), jump_callback::TJC = ContinuousLindbladJumpCallback(), krylov_dim::Int = min(5, cld(length(ψ0.data), 3)), @@ -660,7 +660,7 @@ function dsf_mcsolveEnsembleProblem( e_ops = e_ops₀, params = params2, ntraj = ntraj, - ensemble_method = ensemble_method, + ensemblealg = ensemblealg, jump_callback = jump_callback, prob_func = _dsf_mcsolve_prob_func, progress_bar = progress_bar, @@ -679,8 +679,8 @@ end e_ops::Function=(op_list,p) -> Vector{TOl}([]), params::NamedTuple=NamedTuple(), δα_list::Vector{<:Real}=fill(0.2, length(op_list)), - ntraj::Int=1, - ensemble_method=EnsembleThreads(), + ntraj::Int=500, + ensemblealg::EnsembleAlgorithm=EnsembleThreads(), jump_callback::LindbladJumpCallbackType=ContinuousLindbladJumpCallback(), krylov_dim::Int=max(6, min(10, cld(length(ket2dm(ψ0).data), 4))), progress_bar::Union{Bool,Val} = Val(true) @@ -704,8 +704,8 @@ function dsf_mcsolve( e_ops::Function = (op_list, p) -> (), params::NamedTuple = NamedTuple(), δα_list::Vector{<:Real} = fill(0.2, length(op_list)), - ntraj::Int = 1, - ensemble_method = EnsembleThreads(), + ntraj::Int = 500, + ensemblealg::EnsembleAlgorithm = EnsembleThreads(), jump_callback::TJC = ContinuousLindbladJumpCallback(), krylov_dim::Int = min(5, cld(length(ψ0.data), 3)), progress_bar::Union{Bool,Val} = Val(true), @@ -723,7 +723,7 @@ function dsf_mcsolve( e_ops = e_ops, params = params, ntraj = ntraj, - ensemble_method = ensemble_method, + ensemblealg = ensemblealg, δα_list = δα_list, jump_callback = jump_callback, krylov_dim = krylov_dim, @@ -731,5 +731,5 @@ function dsf_mcsolve( kwargs..., ) - return mcsolve(ens_prob_mc, alg, ntraj, ensemble_method) + return mcsolve(ens_prob_mc, alg, ntraj, ensemblealg) end diff --git a/test/core-test/dynamical-shifted-fock.jl b/test/core-test/dynamical-shifted-fock.jl index bb1bebc16..663ac6494 100644 --- a/test/core-test/dynamical-shifted-fock.jl +++ b/test/core-test/dynamical-shifted-fock.jl @@ -60,7 +60,6 @@ dsf_params, e_ops = e_ops_dsf, progress_bar = Val(false), - ntraj = 500, ) val_ss = abs2(sol0.expect[1, end]) @test sum(abs2.(sol0.expect[1, :] .- sol_dsf_me.expect[1, :])) / (val_ss * length(tlist)) < 0.1 @@ -140,7 +139,6 @@ dsf_params, e_ops = e_ops_dsf2, progress_bar = Val(false), - ntraj = 500, ) val_ss = abs2(sol0.expect[1, end]) diff --git a/test/core-test/time_evolution.jl b/test/core-test/time_evolution.jl index 163a02bf6..2b9d0afe7 100644 --- a/test/core-test/time_evolution.jl +++ b/test/core-test/time_evolution.jl @@ -110,24 +110,22 @@ sol_me2 = mesolve(H, ψ0, tlist, c_ops, progress_bar = Val(false)) sol_me3 = mesolve(H, ψ0, tlist, c_ops, e_ops = e_ops, saveat = saveat, progress_bar = Val(false)) prob_mc = mcsolveProblem(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false)) - sol_mc = mcsolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false)) + sol_mc = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false)) sol_mc2 = mcsolve( H, ψ0, tlist, c_ops, - ntraj = 500, e_ops = e_ops, progress_bar = Val(false), jump_callback = DiscreteLindbladJumpCallback(), ) - sol_mc_states = mcsolve(H, ψ0, tlist, c_ops, ntraj = 500, saveat = saveat, progress_bar = Val(false)) + sol_mc_states = mcsolve(H, ψ0, tlist, c_ops, saveat = saveat, progress_bar = Val(false)) sol_mc_states2 = mcsolve( H, ψ0, tlist, c_ops, - ntraj = 500, saveat = saveat, progress_bar = Val(false), jump_callback = DiscreteLindbladJumpCallback(), @@ -256,7 +254,7 @@ sol_se = sesolve(H_dr_fr, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false)) sol_me = mesolve(H_dr_fr, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false)) - sol_mc = mcsolve(H_dr_fr, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng) + sol_mc = mcsolve(H_dr_fr, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng) # sol_sse = ssesolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng) # Time Evolution in the lab frame @@ -269,17 +267,7 @@ sol_se_td = sesolve(H_td, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p) sol_me_td = mesolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p) - sol_mc_td = mcsolve( - H_td, - ψ0, - tlist, - c_ops, - ntraj = 500, - e_ops = e_ops, - progress_bar = Val(false), - params = p, - rng = rng, - ) + sol_mc_td = mcsolve(H_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng) # sol_sse_td = ssesolve(H_td, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng) @test sol_se.expect ≈ sol_se_td.expect atol = 1e-6 * length(tlist) @@ -292,17 +280,7 @@ sol_se_td2 = sesolve(H_td2, ψ0, tlist, e_ops = e_ops, progress_bar = Val(false), params = p) sol_me_td2 = mesolve(L_td, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p) - sol_mc_td2 = mcsolve( - H_td2, - ψ0, - tlist, - c_ops, - ntraj = 500, - e_ops = e_ops, - progress_bar = Val(false), - params = p, - rng = rng, - ) + sol_mc_td2 = mcsolve(H_td2, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng) # sol_sse_td2 = # ssesolve(H_td2, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), params = p, rng = rng) @@ -609,7 +587,7 @@ @testset "mcsolve, ssesolve and smesolve reproducibility" begin rng = MersenneTwister(1234) - sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng) + sol_mc1 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng) rng = MersenneTwister(1234) sol_sse1 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng) rng = MersenneTwister(1234) @@ -626,7 +604,7 @@ ) rng = MersenneTwister(1234) - sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, ntraj = 500, e_ops = e_ops, progress_bar = Val(false), rng = rng) + sol_mc2 = mcsolve(H, ψ0, tlist, c_ops, e_ops = e_ops, progress_bar = Val(false), rng = rng) rng = MersenneTwister(1234) sol_sse2 = ssesolve(H, ψ0, tlist, c_ops, ntraj = 50, e_ops = e_ops, progress_bar = Val(false), rng = rng) rng = MersenneTwister(1234) @@ -714,7 +692,7 @@ psi0 = kron(psi0_1, psi0_2) t_l = LinRange(0, 20 / γ1, 1000) sol_me = mesolve(H, psi0, t_l, c_ops, e_ops = [sp1 * sm1, sp2 * sm2], progress_bar = false) # Here we don't put Val(false) because we want to test the support for Bool type - sol_mc = mcsolve(H, psi0, t_l, c_ops, ntraj = 500, e_ops = [sp1 * sm1, sp2 * sm2], progress_bar = Val(false)) + sol_mc = mcsolve(H, psi0, t_l, c_ops, e_ops = [sp1 * sm1, sp2 * sm2], progress_bar = Val(false)) @test sum(abs.(sol_mc.expect[1:2, :] .- sol_me.expect[1:2, :])) / length(t_l) < 0.1 @test expect(sp1 * sm1, sol_me.states[end]) ≈ expect(sigmap() * sigmam(), ptrace(sol_me.states[end], 1)) end