Simulação do Modelo de Gray-Scott: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
(3 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 127: Linha 127:


/*Atualização das matrizes. Esse é o núcleo do programa. Integram-se as
/*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:
equações do modelo com FTCS. São usadas as seguintes condições de fronteira
 
periódicas:
* 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
* Se i=0 (primeira linha da matriz), o vizinho da esquerda está na última
Linha 154: Linha 141:
* Se i=L-1 (última linha da matriz, o vizinho da direita está na primeira
* Se i=L-1 (última linha da matriz, o vizinho da direita está na primeira
linha da matriz, com mesmo j).
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.
*/
*/


Linha 170: Linha 153:
     for(j=0; j<L; j++)
     for(j=0; j<L; j++)
     {
     {
       /*Condições totalmente periódicas:
       /*Condições de fronteira periódicas:*/


       uC = *(u+i*L+j);
       uC = *(u+i*L+j);
Linha 183: Linha 166:
       vL = (i==0)? *(v+L2-L+j) : *(v+(i-1)*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));
       vD = (j==0)? *(v+i*L+(L-1)) : *(v+i*L+(j-1));
      */
      /*Condições usadas*/
      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*/
       /*Laplacianos discretizados*/
Linha 228: Linha 196:
   printf("set palette gray\n");
   printf("set palette gray\n");
   printf("set pm3d map\n");
   printf("set pm3d map\n");
  printf("set cbrange [0:1]\n");
   printf("splot \"-\" matrix w image\n");
   printf("splot \"-\" matrix w image\n");


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

Edição atual tal como às 18h54min de 4 de março de 2022

/********************************************************
 * 		 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 as seguintes condições de fronteira
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++)
    {
      /*Condições de fronteira 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));
      
      /*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("set cbrange [0:1]\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+j*L));
     
     printf("\n");
   }

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