Simulação do Modelo de Gray-Scott
Ir para navegação
Ir para pesquisar
/********************************************************
* 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.02;
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) = 1.0-*(u+i*L+j);
}
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 fronteira conforme o Sayama:
* Se i=0 (primeira linha da matriz), o vizinho da esquerda é o vizinho da direita;
* Se j=0 (primeira coluna da matriz), o vizinho de baixo é o vizinho de cima;
* 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).
*********************************************************************************
Também testamos várias simulações com condições completamente 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).
Os comportamentos gerais da formação de padrões são os mesmos, mas as figuras
diferem, principalmente para tempos longos, em que a difusão atinge as
fronteiras do grid.
*/
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++)
{
/*Condições totalmente periódicas:
uC = *(u+i*L+j);
uU = *(u+i*L+(j+1)%L);
uR = *(u+((i+1)%L)*L+j);
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);
vU = *(v+i*L+(j+1)%L);
vR = *(v+((i+1)%L)*L+j);
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));
*/
/*Condições do Sayama*/
uC = *(u+i*L+j);
uU = *(u+i*L+(j+1)%L);
uR = *(u+((i+1)%L)*L+j);
uL = *(u+abs((i-1)%L)*L+j);
uD = *(u+i*L+abs((j-1)%L));
vC = *(v+i*L+j);
vU = *(v+i*L+(j+1)%L);
vR = *(v+((i+1)%L)*L+j);
vL = *(v+abs((i-1)%L)*L+j);
vD = *(v+i*L+abs((j-1)%L));
/*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. Este módulo é praticamente
igual ao criado pelo grupo do modelo de Turing. Apenas foi acrescentado
o comando set palette gray*/
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 gray\n");
printf("set pm3d map\n");
printf("splot \"-\" 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");
}