[no ci] Differentiate tolerances between pure and mixed state dynamics

Author: albertomercurio

src/time_evolution/lr_mesolve.jl

@@ -450,7 +450,7 @@ function lr_mesolveProblem(
     p.Si .= pinv(Hermitian(p.S), atol = opt.atol_inv)
 
     # Initialization of ODEProblem's kwargs
-    default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., saveat = [t_l[end]])
+    default_values = (DEFAULT_ODE_ME_SOLVER_OPTIONS..., saveat = [t_l[end]])
     kwargs2 = merge(default_values, kwargs)
 
     # Initialization of Callbacks

src/time_evolution/mcsolve.jl

@@ -132,7 +132,7 @@ function mcsolveProblem(
     T = Base.promote_eltype(H_eff_evo, ψ0)
 
     # We disable the progress bar of the sesolveProblem because we use a global progress bar for all the trajectories
-    default_values = (DEFAULT_ODE_SOLVER_OPTIONS..., progress_bar = Val(false))
+    default_values = (DEFAULT_ODE_SE_MC_SOLVER_OPTIONS..., progress_bar = Val(false))
     kwargs2 = _merge_saveat(tlist, e_ops, default_values; kwargs...)
     kwargs3 = _generate_mcsolve_kwargs(ψ0, T, e_ops, tlist, c_ops, jump_callback, rng, kwargs2)
 

src/time_evolution/mesolve.jl

@@ -79,7 +79,7 @@ function mesolveProblem(
     ρ0 = to_dense(_CType(T), mat2vec(ket2dm(ψ0).data)) # Convert it to dense vector with complex element type
     L = L_evo.data
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_ME_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_se_me_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2, SaveFuncMESolve)
 
     tspan = (tlist[1], tlist[end])

src/time_evolution/sesolve.jl

@@ -69,7 +69,7 @@ function sesolveProblem(
     ψ0 = to_dense(_CType(T), get_data(ψ0)) # Convert it to dense vector with complex element type
     U = H_evo.data
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_ODE_SE_MC_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_se_me_kwargs(e_ops, makeVal(progress_bar), tlist, kwargs2, SaveFuncSESolve)
 
     tspan = (tlist[1], tlist[end])

src/time_evolution/smesolve.jl

@@ -118,7 +118,7 @@ function smesolveProblem(
     end
     D = DiffusionOperator(D_l)
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_ME_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_stochastic_kwargs(
         e_ops,
         sc_ops_list,

src/time_evolution/ssesolve.jl

@@ -118,7 +118,7 @@ function ssesolveProblem(
     D_l = map(op -> op + _ScalarOperator_e(op, -) * IdentityOperator(prod(dims)), sc_ops_evo_data)
     D = DiffusionOperator(D_l)
 
-    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SOLVER_OPTIONS; kwargs...)
+    kwargs2 = _merge_saveat(tlist, e_ops, DEFAULT_SDE_SE_SOLVER_OPTIONS; kwargs...)
     kwargs3 = _generate_stochastic_kwargs(
         e_ops,
         sc_ops_list,

src/time_evolution/time_evolution.jl

@@ -2,8 +2,11 @@ export TimeEvolutionSol, TimeEvolutionMCSol, TimeEvolutionStochasticSol
 
 export liouvillian_floquet, liouvillian_generalized
 
-const DEFAULT_ODE_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep = false, save_end = true)
-const DEFAULT_SDE_SOLVER_OPTIONS = (abstol = 1e-2, reltol = 1e-2, save_everystep = false, save_end = true)
+# Dynamics of density matrices generally support larger tolerances
+const DEFAULT_ODE_SE_MC_SOLVER_OPTIONS = (abstol = 1e-8, reltol = 1e-6, save_everystep = false, save_end = true)
+const DEFAULT_ODE_ME_SOLVER_OPTIONS = (abstol = 1e-7, reltol = 1e-5, save_everystep = false, save_end = true)
+const DEFAULT_SDE_SE_SOLVER_OPTIONS = (abstol = 5e-3, reltol = 1e-2, save_everystep = false, save_end = true)
+const DEFAULT_SDE_ME_SOLVER_OPTIONS = (abstol = 1e-2, reltol = 1e-2, save_everystep = false, save_end = true)
 const COL_TIMES_WHICH_INIT_SIZE = 200
 
 @doc raw"""