sandbox/benalessio/giererMeinhardt_with_DP.c

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

scalar C1[], C2[];
scalar Ndp[], * tracers = {Ndp};
face vector uf[];

// See Yochelis et al for explanation of parameters.
double D = 0.005, G = 1.0, E = 2.0, S = 0.52;
double setL0 = 16;

double M1 = 0.01, M2 = 0.01, hold_Dp = 1e-3;
const face vector Dp[] = {hold_Dp,hold_Dp};

double dt;
mgstats mgd1, mgd2;

int main()
{
  periodic(right);
  periodic(top);
  Ndp[bottom] = neumann(0);
  Ndp[top] = neumann(0);
  Ndp[left] = neumann(0);
  Ndp[right] = neumann(0);

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

event init (i = 0)
{
  double alpha = 108 - 72 * S * (S + 1) + 2 * pow(S + 1, 3);
  double beta = 1 + 14 * S + pow(S, 2);
  double rho = pow(alpha / 2 + sqrt(pow(alpha, 2) / 4 - pow(beta, 3)), 1.0 / 3);
  double thet = beta / (3 * rho) + rho / 3 - 2 * (S + 1) / 3;
  double phi = -thet - 2 * (S + 1);
  double eta = 4 / sqrt(thet) + phi;
  double C1plus = (sqrt(thet) + sqrt(eta)) / 2;
  double C2plus = (sqrt(thet) * (sqrt(thet) + sqrt(eta))) / \
    (2 - (S - 1) * sqrt(thet) + thet * sqrt(eta));

  foreach() {
    C1[] = C1plus + 0.01*noise(); 
    C2[] = C2plus;// + 0.01*noise();
    // uncomment this portion in order to implement planar wave initial conditions.
    //double x, y;
    //C2[] = 1.0 + 0.1*cos(x) + 0.1*cos(y);
    Ndp[] = 1.0;
  }

  foreach_face()
    uf.x[] = (M1*(C1[] - C1[-1,0]) + M2*(C2[] - C2[-1,0]))/Delta;

}

event tracer_diffusion (i++)
{
  diffusion (Ndp, dt, Dp);
}

// In this event, we update the movie and text file containing all data points.
/*
event movie (i = 0; i += 3)
{
  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 (Ndp, linear = true, spread = 2, file = "Ndp.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\n", t, x, y, C1[], C2[], Ndp[]);

  fclose(out_file);
}
*/
// This event saves a separate text file with just the final time step data.
event final (t = 150)
{
  char name[80];
  //sprintf (name, "mu-%g.png", mu);
  sprintf (name, "C1_final.png");
  output_ppm (C1, file = name, n = 300, linear = true, spread = 2);
  sprintf (name, "C2_final.png");
  output_ppm (C2, file = name, n = 300, linear = true, spread = 2);
  sprintf (name, "Ndp_final.png");
  output_ppm (Ndp, file = name, n = 300, linear = true, spread = 2);

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

}

event integration (i++)
{

  dt = dtnext (1.);

  scalar r[], beta[];
  // This is the Gierer-Meinhardt model.
  foreach() {
    r[] = 0.0;
    beta[] = -1 + C1[]/(C2[]*(1+sq(C1[])));
  }
  const face vector du[] = {D, D};
  mgd1 = diffusion (C1, dt, du, r = r, beta = beta);

  foreach() {
    r[] = G*sq(C1[]) + S;
    beta[] = -E;
  }
  const face vector dv[] = {1.0, 1.0};
  mgd2 = diffusion (C2, dt, dv, r = r, beta = beta);
  foreach_face()
    uf.x[] = (M1*(C1[] - C1[-1,0]) + M2*(C2[] - C2[-1,0]))/Delta;

}

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