Using Split ODE Solvers within Method of Steps

General Usage
1

I am interesting in simulating delayed reaction-diffusion equations of the form

u_t = \Delta u + f(u,u(t-\tau))

in a 2-D spatial domain of modest size (10,000 equations) and a non-stiff nonlinearity (for now). My current progress involves optimizing and benchmarking the ODE function and algorithm/settings choices for the undelayed problem, which splitting methods seem well-suited for at some tolerances and integration lengths. As I understand it, the extension to a DDEProblem more or less just involves rewriting your ODE function to incorporate the history function, then MethodOfSteps just wraps some adaptive time-stepping integrator and makes sure that the integration resolves the discontinuities well. Therefore, I expect at least some of the optimization for the non-delayed case to transfer to the delayed case.

Of course I can probably get away with just using a fully implicit method in MethodOfSteps, but I’d at least like to know if it’s possible to use a splitting method within MethodOfSteps. I will eventually be using a variant of this problem in some Bayesian inference, so trying eke out as much performance as possible is preferable. As I currently understand, to solve a DDE, you create a DDEProblem from a function of the state u and history function h, then the solver alg is just MethodOfSteps(ode_alg_of_choice()). For a SplitODEProblem, I have to give it two functions/a linear operator + function representing the stiff and non-stiff RHS terms then call a split algorithm in solve. It’s not clear to me from my own understanding and from browsing the many examples in the DE.jl documentation on how these two methods can interface (if at all). I’m not at the stage yet to provide a MWE of the delayed problem, but I wanted to know if there was some general framework to stitch together these two features for future use.

Understanding this potential code interface is important for future extensions of this problem such as

  • Introduction of a stiff, but easily separable, component to the nonlinearity
  • Galerkin projection onto a lower dimensional subspace (problem size decreases drastically but matrices/tensors become dense)
  • Increasing tolerances to allow for faster short-time integration
2

It’s in theory possible. I’m not sure if it currently dues the wrapping correctly to keep the special form, but in theory it should be possible. It probably needs an issue. This is something I personally haven’t tried myself.

3

Thanks for the quick response. I haven’t gotten to the stage of actually implementing it yet either, so within the next couple of days I will test some naive methods of wrapping the split solver and post here with the results. For the full-order system, this actually isn’t too big of a deal since the Jacobian of the complete system is so sparse that fully implicit methods still perform alright. The actual PDE I am solving contains multiple fields interacting so the bandwidth is too large for sparse direct solvers, but GMRES handles it pretty well.

For the eventual projected system though, I envision splitting algorithms will provide a marked increase over the fully-implict methods since the Jacobian will be dense and stiff, but one could prefactor the LU decomposition of the stiff part for use in a splitting method. I could still use this LU as a preconditioner to the Jacobian solve in a fully implicit method, but splitting would be ideal I think.