Time-implicit discretisation of reaction–diffusion equations

    We want to discretise implicitly the reaction–diffusion equation \displaystyle \theta\partial_tf = \nabla\cdot(D\nabla f) + \beta f + r where \beta f + r is a reactive term, D is the diffusion coefficient and \theta can be a density term.

    Using a time-implicit backward Euler discretisation, this can be written \displaystyle \theta\frac{f^{n+1} - f^{n}}{dt} = \nabla\cdot(D\nabla f^{n+1}) + \beta f^{n+1} + r Rearranging the terms we get \displaystyle \nabla\cdot(D\nabla f^{n+1}) + (\beta - \frac{\theta}{dt}) f^{n+1} = - \frac{\theta}{dt}f^{n} - r This is a Poisson–Helmholtz problem which can be solved with a multigrid solver.

    #include "poisson.h"

    The parameters of the diffusion() function are a scalar field f, scalar fields r and \beta defining the reactive term, the timestep dt and a face vector field containing the diffusion coefficient D. If D or \theta are omitted they are set to one. If \beta is omitted it is set to zero. Both D and \beta may be constant fields.

    Note that the r, \beta and \theta fields will be modified by the solver.

    The function returns the statistics of the Poisson solver.

    struct Diffusion {
      // mandatory
      scalar f;
      double dt;
      // optional
      face vector D;  // default 1
      scalar r, beta; // default 0
      scalar theta;   // default 1
    mgstats diffusion (struct Diffusion p)

    If dt is zero we don’t do anything.

      if (p.dt == 0.) {
        mgstats s = {0};
        return s;

    We define f and r for convenience.

      scalar f = p.f, r = automatic (p.r);

    We define a (possibly constant) field equal to \theta/dt.

      const scalar idt[] = - 1./p.dt;
      (const) scalar theta_idt = p.theta.i ? p.theta : idt;
      if (p.theta.i) {
        scalar theta_idt = p.theta;
          theta_idt[] *= idt[];

    We use r to store the r.h.s. of the Poisson–Helmholtz solver.

      if (p.r.i)
          r[] = theta_idt[]*f[] - r[];
      else // r was not passed by the user
          r[] = theta_idt[]*f[];

    If \beta is provided, we use it to store the diagonal term \lambda.

      scalar lambda = theta_idt;
      if (p.beta.i) {
        scalar beta = p.beta;
          beta[] += theta_idt[];
        lambda = beta;
      boundary ({lambda});

    Finally we solve the system.

      return poisson (f, r, p.D, lambda);