I have some code for simulating complex differential equations of a certain model that used to work fine, but fails after going from Julia 1.7 to 1.8 and updating all the packages.
using DifferentialEquations, DiffEqCallbacks, SparseArrays function du!(du,u,p,t) # Define differential equations J,Jm,U,γv,D=p mul!(D,Jm,u) @. du = im*J*D - (im*U*abs2(u)+γv)*u end save_func(u,t,integrator)=integrator.p saved_values=SavedValues(Float64,Float64) cb_save=SavingCallback(save_func,saved_values,saveat=tObs) M=size(ψ0,1); mat=build(jp); # Set parameters param=[jp.J, mat.Jm, jp.U, mat.γv, zeros(ComplexF64,M)] #--------------------------------------- prob=ODEProblem(du!, ψ0, (tObs, tObs[end]), param, callback=cb_save) # Define ODE problem & solve sol=solve(prob,Vern8(),dt=dt0,adaptive=true,saveat=tObs)
Without the SavingCallback there is no problem and the
Vern8() solver is more than suffucient with the default tolerances. When the Savingcallback is added, the exact same setup gives incorrect results. In the figure I plot the absolute value of one element of the vector
u that in this model is supposed to stabilize.
I find it counter-intuitive that simply storing a parameter changes the solving algorithm. Am I missing something trivial here?