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
 
(5 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 35: Linha 35:
const int L2=4761;
const int L2=4761;
const double dh=0.014;
const double dh=0.014;
const double dt=0.01;
const double dt=0.02;


double tmax;
double tmax;
Linha 64: Linha 64:
   {
   {
     update(u, v, n_u, n_v);
     update(u, v, n_u, n_v);
 
   
/*Se usado -DGNU na compilação, executa o código abaixo*/   
#ifdef GNU
#ifdef GNU
     if(count % 50 == 0)
     if(count % 50 == 0)
Linha 78: Linha 77:
   }
   }


/*Se usado -DGNU na compilação, executa o código abaixo*/
#ifdef GNU
#ifdef GNU
   gnuplot(u);
   gnuplot(u);
   printf("pause 5\n");
   printf("pause 5\n");
#endif
#endif
/*Gera a saída final (por enquanto, apenas de u), no tempo ~tmax, no arquivo.*/


   sprintf(output,"gs_f_%.3f_k_%.3f_tmax_%.0f.dat",f, k, tmax);
   sprintf(output,"gs_f_%.3f_k_%.3f_tmax_%.0f.dat",f, k, tmax);
Linha 120: Linha 116:
     {
     {
       *(u+i*L+j) = (i*L+j == (L2-1)/2)? 0.0 : 1.0;
       *(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;
       *(v+i*L+j) = 1.0-*(u+i*L+j);
     }
     }


Linha 131: 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 contorno periódicas.
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
* Se i=0 (primeira linha da matriz), o vizinho da esquerda está na última
Linha 156: Linha 153:
     for(j=0; j<L; j++)
     for(j=0; j<L; j++)
     {
     {
      /*Condições de fronteira periódicas:*/
       uC = *(u+i*L+j);
       uC = *(u+i*L+j);
       uR = (i==L-1)? *(u+j) : *(u+(i+1)*L+j);
       uU = *(u+i*L+(j+1)%L);
       uU = (j==L-1)? *(u+i*L) : *(u+i*L+j+1);
       uR = *(u+((i+1)%L)*L+j);
       uL = (i==0)? *(u+L2-L+j) : *(u+(i-1)*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));
       uD = (j==0)? *(u+i*L+(L-1)) : *(u+i*L+(j-1));                


       vC = *(v+i*L+j);
       vC = *(v+i*L+j);
       vR = (i==L-1)? *(v+j) : *(v+(i+1)*L+j);
       vU = *(v+i*L+(j+1)%L);
       vU = (j==L-1)? *(v+i*L) : *(v+i*L+j+1);
       vR = *(v+((i+1)%L)*L+j);
       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));
 
     
       /*Laplacianos discretizados*/
       /*Laplacianos discretizados*/
       uLap = (uR + uL + uU + uD - 4*uC)/(dh*dh);
       uLap = (uR + uL + uU + uD - 4*uC)/(dh*dh);
Linha 185: Linha 184:
  }
  }


  /*Visualização no gnuplot. Autoexplicativo.*/
  /*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)
  void gnuplot(double *arr)
  {
  {
Linha 193: Linha 194:
   printf("set size square\n");
   printf("set size square\n");
   printf("set key off\n");
   printf("set key off\n");
   printf("set palette grey\n");
   printf("set palette gray\n");
   printf("plot \"-\" matrix w image\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 i=0; i<L; i++)
   {
   {
     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");
   
}