# 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.

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

## Initial conditions

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].

## 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) |