Simulação do Modelo de Gray-Scott

/********************************************************
 * 		 Modelo de Gray-Scott		        *
 *******************************************************/
/********************************************************
 *      Grid de 69 x 69 (L x L)                         *
 ********************************************************
 * 	Parâmetros (na ordem):   		        *
 *	tmax F k Du Dv				        *
 *							*
 *	Exemplo:					*
 *	./a.out 1000 0.015 0.055 0.00002 0.00001	*
 *	./a.out tmax   F     k     Du	    Dv		*
 *							*
 *	Saída para gnuplot:         			*
 ********************************************************
 *      compilar com:                                   *
 *      gcc -DGNU gray-scott.c -lm                      *
 ********************************************************
 *	Executar com pipe:		                *
 *	./a.out PARAMETROS | gnuplot			*
 *							*
 ********************************************************
 *      O programa sempre gera um arquivo com o         *
 *      resultado da simulação para tmax. Exemplo:      *
 *      gs_f_0.015_k_0.055_tmax_1000.dat                *
 ********************************************************/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

const int L=69;
const int L2=4761;
const double dh=0.014;
const double dt=0.01;

double tmax;
double Du, Dv;
double f, k;

void init(double *, double *, char *[]);
void update(double *, double *, double *, double *);
void gnuplot(double *);

int main(int argc, char *argv[ ])
{
  FILE *fp;
  char output[100];
  double t = 0;
  double *u, *v;
  double *n_u, *n_v;
  int count = 0;

  u = (double*)calloc(L2, sizeof(double));
  v = (double*)calloc(L2, sizeof(double));
  n_u = (double*)calloc(L2, sizeof(double));
  n_v = (double*)calloc(L2, sizeof(double));

  init(u, v, argv);

  while(t < tmax)
  {
    update(u, v, n_u, n_v);
    
#ifdef GNU
    if(count % 50 == 0)
    {
      printf("set title 'Tempo %.2lf'\n", t);
      gnuplot(u);
    }
#endif
    
    t += dt;
    count += 1;
  }

#ifdef GNU
  gnuplot(u);
  printf("pause 5\n");
#endif

  sprintf(output,"gs_f_%.3f_k_%.3f_tmax_%.0f.dat",f, k, tmax);
  fp = fopen(output, "w");

  for(int i=0; i<L; i++)
  {
    for(int j=0; j<L; j++)
      fprintf(fp, "%lf ", *(u+i*L+j));
  
    fprintf(fp, "\n");
  }

  fclose(fp);

  free(u);
  free(v);
  free(n_u);
  free(n_v);

  return 0;
}

/*INICIALIZAÇÃO:
Escolhendo o ponto do meio do grid para tirar o sistema do equilíbrio,
os demais pontos estão no equilíbrio trivial (u, v) = (1, 0).
Pressupõe-se que L seja ímpar. Se L for par, o ponto inicialmente fora
do equilíbrio não será central.
*/

void init(double *u, double *v, char *argv[ ])
{
  for(int i=0; i<L; i++)
    for(int j=0; j<L; j++)
    {
      *(u+i*L+j) = (i*L+j == (L2-1)/2)? 0.0 : 1.0;
      *(v+i*L+j) = (i*L+j == (L2-1)/2)? 1.0 : 0.0;
    }

  tmax = atof(argv[1]);
  f = atof(argv[2]);
  k = atof(argv[3]);
  Du = atof(argv[4]);
  Dv = atof(argv[5]);
}

/*Atualização das matrizes. Esse é o núcleo do programa. Integram-se as
equações do modelo com FTCS. São usadas condições de contorno periódicas.

* Se i=0 (primeira linha da matriz), o vizinho da esquerda está na última
linha da matriz, com mesmo j; 

* Se j=0 (primeira coluna da matriz), o vizinho de baixo está na última 
coluna da matriz, com mesmo i; 

* Se j=L-1 (última coluna da matriz), o vizinho de cima está na primeira
coluna da matriz, com mesmo i

* Se i=L-1 (última linha da matriz, o vizinho da direita está na primeira
linha da matriz, com mesmo j).
*/

void update(double *u, double *v, double *n_u, double *n_v)
{
  int i, j;
  double uLap, vLap;
  double uC, uR, uL, uU, uD;
  double vC, vR, vL, vU, vD;

  for(i=0; i<L; i++)
    for(j=0; j<L; j++)
    {
      uC = *(u+i*L+j);
      uR = (i==L-1)? *(u+j) : *(u+(i+1)*L+j);
      uU = (j==L-1)? *(u+i*L) : *(u+i*L+j+1);
      uL = (i==0)? *(u+L2-L+j) : *(u+(i-1)*L+j);
      uD = (j==0)? *(u+i*L+(L-1)) : *(u+i*L+(j-1));

      vC = *(v+i*L+j);
      vR = (i==L-1)? *(v+j) : *(v+(i+1)*L+j);
      vU = (j==L-1)? *(v+i*L) : *(v+i*L+j+1);
      vL = (i==0)? *(v+L2-L+j) : *(v+(i-1)*L+j);
      vD = (j==0)? *(v+i*L+(L-1)) : *(v+i*L+(j-1));

      //Laplacianos discretizados
      uLap = (uR + uL + uU + uD - 4*uC)/(dh*dh);
      vLap = (vR + vL + vU + vD - 4*vC)/(dh*dh);

      //Equações do modelo discretizadas
      *(n_u+i*L+j) = uC + (-uC*vC*vC + f*(1-uC) + Du*uLap)*dt;
      *(n_v+i*L+j) = vC + (uC*vC*vC - (f+k)*vC + Dv*vLap)*dt;
    }

  for(i=0; i<L; i++)
    for(j=0; j<L; j++)
    {
      *(u+i*L+j) = *(n_u+i*L+j);
      *(v+i*L+j) = *(n_v+i*L+j);
    }
 }

 /*Visualização no gnuplot. Autoexplicativo.*/
 void gnuplot(double *arr)
 {
   printf("set xrange [0:%d]\nset yrange [0:%d]\n", (L-1), (L-1));
   printf("unset xtics\nunset ytics\n");
   printf("unset colorbox\n");
   printf("set size square\n");
   printf("set key off\n");
   printf("set palette grey\n");
   printf("plot \"-\" matrix w image\n");

   for(int i=0; i<L; i++)
   {
     for(int j=0; j<L; j++)
       printf("%lf ", *(arr+i*L+j));
     
     printf("\n");
   }

   printf("\n\n\ne\n\n");
   
}