# Time-implicit discretisation of reaction–diffusion equations

We want to discretise implicitly the reaction–diffusion equation $\theta {\partial }_{t}f=\nabla \cdot \left(D\nabla f\right)+\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 $\theta \frac{{f}^{n+1}-{f}^{n}}{dt}=\nabla \cdot \left(D\nabla {f}^{n+1}\right)+\beta {f}^{n+1}+r$ Rearranging the terms we get $\nabla \cdot \left(D\nabla {f}^{n+1}\right)+\left(\beta -\frac{\theta }{dt}\right){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, β; // default 0
scalar θ;   // default 1
};

trace
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.θ.i ? p.θ : idt;

if (p.θ.i) {
scalar theta_idt = p.θ;
foreach()
theta_idt[] *= idt[];
}``````

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

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

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

``````  scalar λ = theta_idt;
if (p.β.i) {
scalar β = p.β;
foreach()
β[] += theta_idt[];
λ = β;
}
boundary ({λ});``````

Finally we solve the system.

``````  return poisson (f, r, p.D, λ);
}``````