Operator splitting schemes for reaction-diffusion systems

(Here, I describe our discussion into an issue for later reference.)

Currently, Morpheus solves Systems that contain reaction as well as diffusion by means of operator-splitting. More precisely, it uses a sequential splitting (SEQ) scheme where diffusion and reactions steps are solved independently for a time-step \tau: R,D. The step size is either the user-specified "System/time-step" or it is automatically reduced (halved) until is satisfied the CFL stability criterion.

However, this SEQ operator splitting scheme is known to have an error (wrt the analytical solution) of order 1. Alternative schemes include the so-called Marchuk-Strang (MS) splitting scheme in which first, diffusion is computed for half the time-step then reaction is computed for the full time-step, followed by a second half time-step of diffusion. Schemetically: 0.5D, R, 0.5D. A third option is the so-called weighted sequential splitting in which (in the symmetric case) one computes both sequential splittings and take the arithmetic mean of the results: (D,R + R,D)/2. Both alternatives have second order accuracy (when the reaction step is solved using Heun or Runge-Kutta solvers).

For more information, see this paper:

"Application of Operator Splitting to Solve Reaction-Diffusion Equations" Tamas Ladics, Electronic Journal of Qualitative Theory of Differential Equations, 2012. http://www.kurims.kyoto-u.ac.jp/EMIS/journals/EJQTDE/p1082.pdf

For Morpheus, this would imply that Systems that include both reactions and diffusion should be scheduled separately, in a "reaction-diffusion block". Within the block, we can implement the different operator-splitting schemes mentioned above. These could be selected by the user from XML (probably in Field/Diffusion).

Thanks to @alvarokohn for the reference.