brusselator_with_DP.c

References

This code calculates the diffusiophoretic motion of two colloidal species in response to the solute gradients generated by a Brusselator reaction diffusion system. It is based off of the Basilisk example, where we simply add advective diffusive tracers which advect porportionally to the gradient of the solutes.

There are two pattern-forming chemicals in the Brusselator model, which in this code are called C1 and C2. The number of colloidal species is up to the user, and is set to two in this example (N1 and N2).

This plot shows the spatial distribution of C1 at the final simulated time step.

This plot shows the spatial distribution of C2 at the final simulated time step.

In this plot, we show the spatial distribution of N1 at the final simulated time step.

In this plot, we show the spatial distribution of N2 at the final simulated time step.

#include "run.h"
#include "diffusion.h"
#include "advection.h"
#include "tracer.h"

C1 and C2 are the chemical concentrations, N1 and N2 are the diffusiophoretic colloids, are uf1 and uf2 are their respective diffusiophoretic velocities. Note that the Brusselator model involves two pattern-forming species, C1 and C2, but the number of colloidal species, which form their own patterns by advecting along the gradients of C1 and C2, is determined by the user. In this example, we have two colloidal species.

scalar C1[], C2[];
scalar N1[], N2[], * tracers = {N1, N2};
face vector uf1[];
face vector uf2[];

Set no-flux boundary conditions.

N1[bottom] = neumann(0);
N1[top] = neumann(0);
N1[left] = neumann(0);
N1[right] = neumann(0);
N2[bottom] = neumann(0);
N2[top] = neumann(0);
N2[left] = neumann(0);
N2[right] = neumann(0);

k is the Damköhler number, ka is the initial concentration of C1, D is the ratio of diffusivities of C1 and C2, mu in the supercriticality, setL0 is the domain length (it’s a square), M’s are the diffusiophoretic mobilities, and D1 and D2 are the colloid diffusivities. See Peña and Pérez-García, 2001 for the background on the Brusselator reaction-diffusion system parameters.

double k = 1., ka = 1.5, D = 4;
double mu = 0.05, kb;
double setL0 = 64;
double M1_one = 0.025, M2_one = 0.035;
double M1_two = -0.025, M2_two = -0.035;
double hold_D1 = 1e-3, hold_D2 = 1e-3;
const face vector D1[] = {hold_D1,hold_D1};
const face vector D2[] = {hold_D2,hold_D2};

double dt;
mgstats mgd1, mgd2;

int main()
{
  init_grid (128);
  size (setL0);
  TOLERANCE = 1e-4;
  run();
}

event init (i = 0)
{
  double nu = sqrt(1./D);
  double kbcrit = sq(1. + ka*nu);
  kb = kbcrit*(1. + mu);

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

event tracer_diffusion (i++)
{
  diffusion (N1, dt, D1);
  diffusion (N2, dt, D2);
}

event tracer_advection (i++)
{
  foreach_face()
    uf1.x[] = (M1_one*(C1[] - C1[-1,0]) + M2_one*(C2[] - C2[-1,0]))/Delta;

  foreach_face()
    uf2.x[] = (M1_two*(C1[] - C1[-1,0]) + M2_two*(C2[] - C2[-1,0]))/Delta;

  advection({N1}, uf1, dt);
  advection({N2}, uf2, dt);
}

// In this event, we update the movie and text file containing all data points.
/*
event movie (i = 0; i += 20)
{
  output_ppm (C1, linear = true, spread = 2, file = "C1.mp4", n = 600);
  output_ppm (C2, linear = true, spread = 2, file = "C2.mp4", n = 600);
  output_ppm (N1, linear = true, spread = 2, file = "N1.mp4", n = 600);
  output_ppm (N2, linear = true, spread = 2, file = "N2.mp4", n = 600);
  fprintf (stderr, "%d %g %g %d %d\n", i, t, dt, mgd1.i, mgd2.i);

  if (i==0)
  {
    fclose(fopen("dataAll.txt", "w"));
  }
  FILE *out_file = fopen("dataAll.txt", "a"); // write only
  foreach() {
    fprintf(out_file, "%f %f %f %f %f %f %f\n", t, x, y, C1[], C2[], N1[], N2[]);
  }
  fclose(out_file);
}
*/
// This event saves a separate text file with just the final time step data.
event final (t = 1200)
{
  output_ppm (C1, file = "C1_final.png", n = 300, linear = true, spread = 2);
  output_ppm (C2, file = "C2_final.png", n = 300, linear = true, spread = 2);
  output_ppm (N1, file = "N1_final.png", n = 300, linear = true, spread = 2);
  output_ppm (N2, file = "N2_final.png", n = 300, linear = true, spread = 2);

  /*
  fclose(fopen("dataFinal.txt", "w"));
  FILE *out_file = fopen("dataFinal.txt", "a"); // write only
  foreach() {
    fprintf(out_file, "%f %f %f %f %f %f %f\n", t, x, y, C1[], C2[], N1[], N2[]);
  }
  fclose(out_file);
  */
}

event integration (i++)
{
  dt = dtnext (1.);

  scalar r[], beta[];
  
  // This is the Brusselator reaction diffusion system.
  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);
}

event adapt (i++) {
  adapt_wavelet ({N1,N2,uf}, (double[]){2e-1,2e-1,0.1,0.1},
          9, 5);
}

References

[alessio2023diffusiophoresis]	Benjamin M Alessio and Ankur Gupta. Diffusiophoresis-enhanced turing patterns. Science Advances, 45, 2023. [ http ]
[pena2001stability]	B Pena and C Perez-Garcia. Stability of turing patterns in the brusselator model. Physical review E, 64(5):056213, 2001. [ http ]