sandbox/benalessio/cellCellInteraction_with_DP.c

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

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

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);
C3[bottom] = neumann(0);
C3[top] = neumann(0);
C3[left] = neumann(0);
C3[right] = neumann(0);

// See [Nakamasu et al.](https://www.pnas.org/doi/epdf/10.1073/pnas.0808622106) for explanation of parameters.
double c1 = -0.04, c2 = -0.055, c3 = 0.37, c4 = -0.05, c5 = 0.0, c6 = 0.25;
double c7 = 0.016, c8 = -0.03, c9 = 0.24, cu = 0.02, cv = 0.025, cw = 0.06;
double Du = 0.02, Dv = 0.02, Dw = 0.2, U = 0.5, V = 0.5, W = 0.5;
double setL0 = 64;
// M1 and M2 are the diffusiophoretic mobilities.
double M1 = 0.05, M2 = -0.05;

double dt;

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

event init (i = 0)
{
  foreach() {
    N1[] = 1.0 + 0.1*noise();
    N2[] = 1.0 + 0.1*noise();
    double spot1 = 0.4*setL0;
    double spot2 = 0.6*setL0;
    // This initial condition imprints a pattern.
    C3[] = 1.0 + 0.1*cos(2*3.14*y/setL0*10)\
      *((x <= spot1) + (x >= spot2))*((y <= spot1) + (y >= spot2))\
      + ((x > spot1) * (x < spot2))*((y > spot1) * (y < spot2))\
      *(1.0 + 0.1*noise());
    // Uncomment either of these lines to have uniform noise or sine wave instead.
    //C3[] = 1.0 + 0.1*noise();
    //C3[] = 1.0 - 0.1*cos(2*3.14*y/setL0*2);
  }

}

event tracer_diffusion (i++)
{
  scalar r1[], beta1[];
  foreach() {
    beta1[] = -cu;
    r1[] = (c1*N2[] + c2*C3[] + c3)*(
      (c1*N2[] + c2*C3[] + c3) >= 0)*(U >= (c1*N2[] + c2*C3[] + c3)) + (
      U < (c1*N2[] + c2*C3[] + c3))*U;
  }
  const face vector D1[] = {Du,Du};
  diffusion (N1, dt, D1, r = r1, beta = beta1);

  scalar r2[], beta2[];
  foreach() {
    beta2[] = -cv;
    r2[] = (c4*N1[] + c5*C3[] + c6)*(
      (c4*N1[] + c5*C3[] + c6) >= 0)*(V >= (c4*N1[] + c5*C3[] + c6)) + (
      V < (c4*N1[] + c5*C3[] + c6))*V;
  }
  const face vector D2[] = {Dv,Dv};
  diffusion (N2, dt, D2, r = r2, beta = beta2);
}

event tracer_advection (i++)
{
  foreach_face()
    uf1.x[] = M1 * (C3[] - C3[-1,0]) / Delta;

  foreach_face()
    uf2.x[] = M2 * (C3[] - C3[-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 += 30)
{
  output_ppm (N1, linear = true, spread = 2, file = "N1.mp4", n = 300);
  output_ppm (N2, linear = true, spread = 2, file = "N2.mp4", n = 300);
  output_ppm (C3, linear = true, spread = 2, file = "C3.mp4", n = 300);

  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, N1[], N2[], C3[]);
  }
  fclose(out_file);
}
*/
// This event saves a separate text file with just the final time step data.
event final (t = 1000)
{
  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);
  output_ppm (C3, file = "C3_final.png", 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, N1[], N2[], C3[]);
  }
  fclose(out_file);
  */
}

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

  scalar r[], beta[];
  
  foreach() {
    r[] = (c7*N1[] + c8*N2[] + c9)*(
      (c7*N1[] + c8*N2[] + c9) >= 0)*(W >= (c7*N1[] + c8*N2[] + c9)) + (
      W < (c7*N1[] + c8*N2[] + c9))*W;;
    beta[] = -cw;
  }
  const face vector D3[] = {Dw,Dw};
  diffusion (C3, dt, D3, r = r, beta = beta);

}

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