src/examples/brusselator.c

Coupled reaction–diffusion equations

The Brusselator is a theoretical model for a type of autocatalytic reaction. The Brusselator model was proposed by Ilya Prigogine and his collaborators at the Free University of Brussels.

Two chemical compounds with concentrations C_1 and C_2 interact according to the coupled reaction–diffusion equations: \displaystyle \partial_t C_1 = \nabla^2 C_1 + k(ka - (kb + 1)C_1 + C_1^2 C_2) \displaystyle \partial_t C_2 = D \nabla^2 C_2 + k(kb C_1 - C_1^2 C_2)

We will use a Cartesian (multi)grid, the generic time loop and the time-implicit diffusion solver.

#include "grid/multigrid.h"
#include "run.h"
#include "diffusion.h"

We need scalar fields for the concentrations.

scalar C1[], C2[];

We use the same parameters as Pena and Perez-Garcia, 2001

double k = 1., ka = 4.5, D = 8.;
double mu, kb;

The generic time loop needs a timestep. We will store the statistics on the diffusion solvers in mgd1 and mgd2.

double dt;
mgstats mgd1, mgd2;

Parameters

We change the size of the domain L0 and set the tolerance of the implicit diffusion solver.

int main()
{
  init_grid (128);
  size (64);
  TOLERANCE = 1e-4;

Here \mu is the control parameter. For \mu > 0 the system is supercritical (Hopf bifurcation). We test several values of \mu.

  mu = 0.04; run();
  mu = 0.1;  run();
  mu = 0.98; run();
}

Initial conditions

event init (i = 0)
{

The marginal stability is obtained for kb = kbcrit.

  double nu = sqrt(1./D);
  double kbcrit = sq(1. + ka*nu);
  kb = kbcrit*(1. + mu);

The (unstable) stationary solution is C_1 = ka and C_2 = kb/ka. It is perturbed by a random noise in [-0.01:0.01].

  foreach() {
    C1[] = ka ; 
    C2[] = kb/ka + 0.01*noise();
  }
}

Outputs

Here we create an mpeg animation of the C_1 concentration. The spread parameter sets the color scale to \pm twice the standard deviation.

event movie (i = 1; i += 10)
{
  output_ppm (C1, linear = true, spread = 2, file = "f.mp4", n = 200);
  fprintf (stderr, "%d %g %g %d %d\n", i, t, dt, mgd1.i, mgd2.i);
}

We make a PNG image of the final “pseudo-stationary” solution.

event final (t = 3000)
{
  char name[80];
  sprintf (name, "mu-%g.png", mu);
  output_ppm (C1, file = name, n = 200, linear = true, spread = 2);
}

Time integration

event integration (i++)
{

We first set the timestep according to the timing of upcoming events. We choose a maximum timestep of 1 which ensures the stability of the reactive terms for this example.

  dt = dtnext (1.);

We can rewrite the evolution equations as \displaystyle \partial_t C_1 = \nabla^2 C_1 + k k_a + k (C_1 C_2 - k_b - 1) C_1 \displaystyle \partial_t C_2 = D \nabla^2 C_2 + k k_b C_1 - k C_1^2 C_2 And use the diffusion solver to advance the system from t to t+dt.

  scalar r[], beta[];
  
  foreach() {
    r[] = k*ka;
    beta[] = k*(C1[]*C2[] - kb - 1.);
  }
  mgd1 = diffusion (C1, dt, r = r, beta = beta);
  foreach() {
    r[] = k*kb*C1[];
    beta[] = - k*sq(C1[]);
  }
  const face vector c[] = {D, D};
  mgd2 = diffusion (C2, dt, c, r, beta);
}

Results

We get the following stable Turing patterns.

\mu=0.04	\mu=0.1 (stripes)	\mu=0.98 (hexagons)