Simulação do Modelo de Gray-Scott

<source lang = "c">

/********************************************************

* 		 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                *
********************************************************/

1. include<stdio.h>
2. include<stdlib.h>
3. 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);
   

1. ifdef GNU

   if(count % 50 == 0)
   {
     printf("set title 'Tempo %.2lf'\n", t);
     gnuplot(u);
   }

1. endif

   t += dt;
   count += 1;
 }

1. ifdef GNU

 gnuplot(u);
 printf("pause 5\n");

1. 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");
  

}