I’m wondering if there is already a mechanism for ugly coupling of differential equations.
In my application, I have one large-dimensional SDE and a reaction-diffusion PDE, which are coupled and have different time scales. Ideally, one might want to perform smaller time steps for the PDE, but depending on the parameters, it can swap, and the SDE requires smaller time steps due to dynamic stiffness.
I would like to solve both systems with an uncoupled adaptive time-stepping.
(In addition, the coupling also has some non-mathematical parts like updating neighbor lists periodically, which means one has callbacks that affect both systems in a sense.)
Below is a minimal example:
using OrdinaryDiffEq
ts = LinRange(0, 40, 10)
ode1 = ODEProblem((u,p,t) -> 1.0, 0.0, (0, 40), nothing, tstops = ts)
ode2 = ODEProblem((u,p,t) -> -u[1], 0.0, (0, 40), nothing, tstops = ts)
integrators = (init(ode1, Heun()), init(ode2, Heun()))
stepnext!(intgrs) = step!(argmin(i -> i.t, intgrs))
tmin(intgrs) = minimum((integr.t for integr in intgrs))
currentsol(intgrs) = (a = intgrs[1].u, b = intgrs[2].u, t1 = intgrs[1].t, t2 = intgrs[2].t)
# start solution
reinit!.(integrators)
sol = [currentsol(integrators)]
for t in ts
while tmin(integrators) < t
stepnext!(integrators)
end
# possibly do ugly stuff...
push!(sol, currentsol(integrators))
end
Is there a more direct way which allows to mix systems with such a staggered adaptive time-stepping?
[I am aware that this is an intermediate mathematical crime, as ideally one would prefer an error estimate for the complete system and the individual error estimates are blind to possible effects from the coupling, but I am quite certain that I can ignore these feedback effects.]