Time-implicit discretisation of reaction–diffusion equations

We want to discretise implicitly the reaction–diffusion equation θtf=(Df)+βf+r where βf+r is a reactive term, D is the diffusion coefficient and θ can be a density term.

Using a time-implicit backward Euler discretisation, this can be written θfn+1fndt=(Dfn+1)+βfn+1+r Rearranging the terms we get (Dfn+1)+(βθdt)fn+1=θdtfnr 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 β defining the reactive term, the timestep dt and a face vector field containing the diffusion coefficient D. If D or θ are omitted they are set to one. If β is omitted it is set to zero. Both D and β may be constant fields.

Note that the r, β and θ 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

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 θ/dt.

  const scalar idt[] = - 1./p.dt;
  (const) scalar theta_idt = p.θ.i ? p.θ : idt;
  if (p.θ.i) {
    scalar theta_idt = p.θ;
      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 β is provided, we use it to store the diagonal term λ.

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

Finally we solve the system.

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