Exibir código-fonte para Equação de Cahn-Hilliard em 2D

'''Leonardo Dasso Migottto'''

O objetivo deste trabalho é resolver computacionalmente a equação de Cahn-Hilliard, utilizando a Transformada Rápida de Fourier <ref>https://fiscomp.if.ufrgs.br/index.php/FFT</ref> em uma e (principalmente) em duas dimensões. Será explorado as variações em concentração inicial e seus respectivos padrões formados, dados coeficientes de difusão e dimensões da superfície fixos.

[[Arquivo:NOPBC03.gif|thumb|Cahn-Hilliard bidimensional com condições de contorno periódicas.]]

Esta equação já foi tratada em detalhes por colegas anteriores a mim <ref>https://fiscomp.if.ufrgs.br/index.php/Equa%C3%A7%C3%A3o_de_Cahn-Hilliard</ref>, e a leitura do trabalho por eles desenvolvido é recomendada para melhor entendimento da equação. O foco deste trabalho é explorar a solução numérica para a equação quando tratada em duas dimensões, onde a formação de padrões apresenta resultados mais interessantes. No entanto, a fim de facilitar a implementação e entendimento em duas dimensões, também será exibida uma implementação em uma dimensão.

'''Todos códigos desenvolvidos para duas dimensões foram feitos a fim de que qualquer pessoa, com Python e as bibliotecas necessárias instaladas, possa executá-los sem nenhum conhecimento prévio. Estes códigos podem ser encontrados no final desta página.'''

== Equação de Cahn-Hiliiard Utilizando Transformada de Fourier ==
Para encontrar a equação que implementaremos com o uso da Transformada Rápida de Fourier, precisamos encontrar a nossa equação representada no espaço de Fourier. Seguirei a literatura de S. Bulent Biner <ref>BINER, S. Bulent Biner. Programming Phase-Field Modeling. Idaho National Laboratory, Idaho Falls, ID, USA: Springer, 2017. 400 p. ISBN 978-3-319-41196-5.</ref>, onde há um capítulo dedicado a resolver equações de difusão com métodos desse tipo. Primeiro, resolveremos a equação em uma dimensão, que segue abaixo:

:<math>
\frac{\partial c}{\partial t} = D {\nabla}^2 (c^3 - c - \gamma {\nabla}^2 c)
</math>

Lembrando que, conforme explicado no trabalho anterior a este, <math>c</math> é a soma da concentração de duas fases, e seu valor varia de -1 (presença de só uma das fases) a 1 (presença somente da outra fase). Em uma dimensão, os laplacianos podem ser substituídos pela derivada segunda em relação a <math>x</math>, resultando na seguinte equação:

:<math>
\frac{\partial c}{\partial t} = D \Bigl( \frac{\partial^2 (c^3 - c)}{\partial x^2} - \gamma \frac{\partial^4 c}{\partial x^4} \Bigr)
</math>

Para solucioná-la numericamente, aplicaremos a Transformada de Fourier à frente em ambos os lados, da maneira descrita abaixo, onde k é o respectivo coeficiente de Fourier):

:<math>
\frac{\partial \{c\}_k}{\partial t} = D \Bigl( \Bigl\{ \frac{ \partial^2 c^3 - c}{\partial x^2} \Bigr\}_k - \gamma \Bigl\{ \frac{\partial^4 c}{\partial x^4} \Bigr\}_k \Bigr)
</math>

Em seguida, substituimos as derivadas espaciais pela sua equivalente no espaço de Fourier:

:<math>
\Bigl\{ \frac{\partial^n c}{\partial x^n} \Bigr\}_k \to (ik)^n \{c\}_k
</math>

Assim, obtemos a seguinte equação:

:<math>
\frac{\partial \{c\}_k}{\partial t} = D \bigl(-k^2 (\{c^3\}_k - \{c\}_k) - \gamma k^4 \{c\}_k \bigr)
</math>

O próximo passo é fazer a derivada à direita quanto ao tempo, da seguinte maneira:

:<math>
\frac{\partial \{c\}_k}{\partial t}\to \frac{\{c\}_k^{n+1}-\{c\}_k^n}{\Delta t}
</math>

Substituindo na equação e reescrevendo-a a fim de isolar <math>\{c\}_k^{n+1}</math>, obtemos a equação final:

:<math>
\{c\}_k^{n+1} = \{c\}_k^n + D \Delta t \bigl(-k^2 (\{c^3\}_k^n - \{c\}_k^n) - \gamma k^4 \{c\}_k^n \bigr)
</math>

A equação que descobrimos para uma dimensão possui a mesma forma quando em duas dimensões. No entanto, <math>k</math> representará um vetor com coordenada <math>(k_{1}, k_{2})</math> com módulo <math>\sqrt{k_1^2 + k_2^2}</math>, onde <math>k_1</math> e <math>k_2</math> são os coeficientes em <math>x</math> e <math>y</math> respectivamente.

== Código em Uma Dimensão ==
O código completo está disponível no final desta página. Abaixo há o excerto da funcão, em Python, que calcula o instante seguinte do nosso sistema, utilizando as funções RFFT e IRFFT do pacote Scipy.
<source lang="python">
def cahnfourier1d(cc, k2, k4):
    cct = rfft(cc)
    cct3 = rfft(cc**3 - cc)
    cct = cct + difd*dt*(-k2*(cct3) - k4*cct)
    cc = irfft(cct)
    return cc
</source>
No código, k2 e k4 são os coeficientes elevados às respectivas potências, além de que k4 está multiplicado pelo valor de <math>\gamma</math>. Calcular o quadrado de k2 e multiplica-lo por gamma fora da função evita que essas multiplicações sejam feito toda vez que a função é chamada. Também é importante destacar que as funções RFFT e IRFFT funcionam utilizando o vetor completo de pontos para fazer os cálculos, eliminando qualquer necessidade de iteração sobre o vetor de valores (uma diferença notável do método FTCS, onde é necessário iterar sobre o vetor).

== Resultados em Uma Dimensão e Discussão==
Como já há um trabalho que trata em detalhes a implementação unidimensional e seus resultados, irei comparar aqui ambas implementações. Abaixo, vemos alguns instantes comparando ambos métodos a partir de uma condição inicial aleatória, utillizando condições de contorno periódicas, com o maior valor de erro destacado no topo. O valor das constantes relevantes são:

<math>\Delta t = \frac{1}{22 \cdot 10^{6}}, \Delta x = \frac{1}{128}, D = 1, \gamma = \Bigl(\frac{3.4}{128}\Bigr)^2, x_{max} = 1</math>.

[[Arquivo:1d100passos.png]] [[Arquivo:1d1000passos.png]] [[Arquivo:1d10000.png]]

Como podemos ver, a diferença dos valores entre os resultados obtidos pelo método FTCS e o método das transformadas é pequena (após poucos instantes o maior módulo da diferença aproxima-se de 0,01). É interessante citar que a escolha de <math>dx</math> como um valor tão específico não é em vão: a função RFFT utilizada no código apresenta maior rapidez de execução quando o vetor utilizado é de tamanho <math>2^n , n = 1, 2, 3, ...</math>. Com este valor de <math>dx</math>, nosso vetor é composto de 128 elementos, ou <math>2^7</math> elementos, fazendo uso dessa vantagem. Dado a escala utilizada no sistema acima, os tempos de processamento foram pequenos, no entanto, para sistemas maiores (como o bidimensional), onde o número de operações necessárias aumenta drasticamente, o método das transformadas será distintamente ágil.

== Código em Duas Dimensões ==
Todos códigos estão disponíveis na íntegra no final desta página. A fim de poupar tempo, quatro programas foram criados:
:- Um programa que gera, a partir de valores aleatórios parametrizados, dados até um certo instante de tempo, salvando-os em um arquivo.
:- Um programa que, lê um arquivo de valores já gerados e gera mais valores até um certo instante de tempo, salvando-os neste mesmo arquivo.
:- Um programa que lê um arquivo gerado previamente, criando frames dos valores salvos no arquivo.
:- Um programa que lê um arquivo gerado previamente, criando uma animação dos valores salvos no arquivo.

Dois arquivos de formato .npy são criados, para armazenar os arrays e as informações necessárias para gerar mais valores ou plotar os gráficos, além de um .txt que serve como uma consulta local dos valores pelo usuário (arquivos .npy não podem ser lidos no bloco de notas). Todos arquivos são salvos em pastas individuias cujo nome é um número aleatório único, sendo este número o nome do arquivo de texto, o do arquivo .npy e o arquivo com os valores calculados.
Abaixo há o excerto da funcão, em Python, que calcula o instante seguinte do nosso sistema, utilizando as funções RFFT2 e IRFFT2 do pacote Scipy.
<source lang="python">
def cahnfourier2d(aa, kk2, kk4):
    cct = rfft2(aa)
    cct3 = rfft2(aa**3)
    cct = cct + difd*dt*(-kk2*(cct3 - cct) - kk4*cct)
    ccn = irfft2(cct)
    return ccn
</source>
No código, k2 e k4 são os coeficientes elevados às respectivas potências, além de que k4 está multiplicado pelo valor de <math>\gamma</math>. Vale destacar a semelhança entre este código e o unidimensional, fruto da discretização quase idêntica.

== Resultado em Duas Dimensões: Concentração Média 0 ==
Para todos resultados abaixo, foram utilizados as seguintes constantes relevantes:

<math>\Delta t = \frac{1.8}{10^7}, \Delta x = \Delta y = \frac{1}{128}, D = 1, \gamma = 0.01^2, x_{max} = 1</math>.

Também foram utilizadas condições de contorno periódicas. Como observado anteriormente, o valor de <math>dx</math> escolhido serve para que nosso array de valores possua <math>10^{14}</math> elementos, também fazendo uso da eficiência completa da função RFFT2 (equivalente bidimensional da função RFFT).

O motivo de não explorarmos valores diferentes de <math>\gamma</math> nem de <math>D</math> está na sensibilidade da estabilidade da equação: pequenas mudanças na ordem de grandeza de qualquer uma das constantes pode resultar num sistema instável.

Segue abaixo mídias da difusão sob estes parâmetros com concentração média 0. As cores variam de vermelho (-1) a violeta (1), de modo que as cores variam de acordo com o espectro visível, onde a cor amarela corresponde a 0.

A "baixa qualidade" das imagens se dá pelo tamanho de array utilizado (128x128). É possível executar o programa com um <math>\Delta x</math> menor, logo um array maior, porém o tamanho escolhido já permite analisarmos os dados, além de ter um tempo de processamento menor (lembre-se que estamos utilizando um <math>\Delta t</math> minúsculo, além de que, utilizar metade do <math>\Delta x</math> resultará em um array com 4 vezes o número de células).

[[Arquivo: NOPBCINICIAL.png|thumb|none|256px|Instante inicial com concentrações aleatórias.]]

[[Arquivo: NOPBCCURTO.gif|frame|none|256px|Tempo máximo da animação: 0.012. Repare na breve homogenização do sistema, seguido da rápida formação das fases.]]

[[Arquivo: PBC BEM CEDO1.png|thumb|none|384px|Instante expandido a fim de demonstrar as condições de contorno periódicas. Esta imagem é composta por uma grade 3x3 das de tamanho menor.]]

[[Arquivo: NOPBC03.gif|frame|none|256px|Tempo máximo da animação: 0.3. Este vídeo está acelerado em relação ao anterior. Repare na difusão das "bolhas" que não entraram em contato com outras concentrações.]]

[[Arquivo: NOPBC3NOVO.gif|frame|none|256px|Tempo máximo da animação: 3. Este vídeo está acelerado em relação ao anterior. Repare no arredondamento da concentração se formou.]]

[[Arquivo: MOSTRANDO ESTABILIDADE.gif|frame|none|512px|Tempo máximo da animação: 6. Este vídeo está acelerado em relação ao anterior. Utilizando da visualização expandida, podemos ver que, ao atingir a estabilidade, uma das fases formou uma concentração circular, enquanto a outra ocupou o espaço remanescente.]]

[[Arquivo: PBCFINAL5050.png|thumb|none|384px|Último instante da simulação, no qual as fases atingiram a estabilidade.]]

== Resultado em Duas Dimensões: Concentração Média 0.5 ==
Abaixo seguem mídias da simulação de condições iniciais idênticas às anteriores, porém com concentração média inicial 0.5 (equivalente a 75% em uma fase e 25% em outra).

[[Arquivo: FRAMEINICIAL.png|thumb|none|256px|Instante inicial com concentrações aleatórias. Repare nas cores serem majoritariamente verdes ou azuis, relativas à concentrações de valor positivo.]]

[[Arquivo: TEMPOCURTO.gif|frame|none|256px|Tempo máximo da animação: 0.012. Repare na breve homogenização do sistema, seguido da rápida formação de concentrações pequenas da fase em minoria.]]

[[Arquivo: DESTAQUEPBC.png|thumb|none|384px|Instante expandido a fim de demonstrar as condições de contorno periódicas. Esta imagem é composta por uma grade 3x3 das de tamanho menor.]]

[[Arquivo: TMAX03.gif|frame|none|256px|Tempo máximo da animação: 0.3. Este vídeo está acelerado em relação ao anterior. Repare na difusão das "bolhas" que não entraram em contato com outras concentrações.]]

[[Arquivo: TMAX3.gif|frame|none|256px|Tempo máximo da animação: 3. Este vídeo está acelerado em relação ao anterior. Repare na formação de 3 concentrações já circulares de forma, além do início da difusão de uma delas.]]

[[Arquivo:GIFINAL GRANDE.gif|frame|none|512px|Tempo máximo da animação: 7. Este vídeo está acelerado em relação ao anterior. Utilizando da visualização expandida, podemos ver novamente que, ao atingir a estabilidade, uma das fases formou uma concentração circular, enquanto a outra ocupou o espaço remanescente.]]

[[Arquivo: FINALPEQUENOCONC.png|thumb|none|384px|Último instante da simulação, no qual as fases atingiram a estabilidade.]]

== Discussão dos resultados em duas dimensões==
Para tempos curtos, vemos o fenômeno da homogenização do sistema, seguida da formação das regiões das diferentes fases, conforme as concentrações aleatórias difundem. Essa breve homogenização se dá por nosso sistema de concentrações aleatórias ser equivalente a uma mistura homogênea com perturbações.

Após poucos instantes, já percebemos a separação completa das fases, com exceção das fronteiras, que apresentam uma mistura de ambas. É importante destacar que, conforme maiores as concentrações de fases se tornam, menor a velocidade que elas se "movem". Ao lembrar que a equação de Cahn-Hilliard descreve a separação de fases em um sistema homogêneo termodiamicamente instável, que reduz sua energia conforme a separação delas, fica evidente que a desaçeleração da separação ocorre pela aproximação da sua energia atual à energia mínima do sistema.

Percebemos que o estado de menor energia é uma concentração circular, enquanto o restante do meio é ocupado pela outra. Mesmo sem saber que este é o estado estacionário, poderíamos intuir que seria ele ao ver os instantes iniciais da simulação: o círculo é a forma que possui maior área por perímetro, e vimos que as fases diferentes morfam de modo a diminuir as fronteiras entre elas. No caso de uma quantidade diferente entre as fases, a de menor volume (ou neste caso, área) total, forma o círculo, caso contrário, a concentração que forma é escolhida "aleatoriamente" (depende dos locais das perturbações, que, no nosso caso, são aleatórias).

== Códigos ==
Todos códigos referentes a duas dimensões possuem, nas primeiras linhas, comentários que destacam quais valores o usuário é recomendado a alterar, além da semente aleatória utilizada para gerar as mídias vistas na página. Caso você não possua experiência com programação, é recomendado não alterar as linhas de código abaixo dos locais indicados.

[[cahn2dfourier]] Programa utilizado para gerar valores sem partir de nenhum arquivo anterior. Os dados são identificados por um código numérico único, e serão salvos em uma pasta com nome de tal código, onde ali se encontrarão:

:- O arquivo código.npy principal, que consiste nos arrays.

:- O arquivo valcódigo.npy, que consiste em um arquivo que é lido pelos outros programas contendo informações cruciais.

:- O arquivo código.txt, que contém tais informações cruciais em um formato legível ao usuário.

Utilizado em conjunto com os programas "framer" e "animador", você pode gerar imagens como as vistas neste trabalho.

[[cahn2dfouriercontinue]] Programa utilizado para gerar valores partindo de outros já antes calculados. O arquivo com os dados novos (concatenados aos antigos) se encontrará na pasta original, substituindo o arquivo antigo. Recomendado para gerar valores por tempos longos sem precisar deixar o programa executando por longas horas sem intervalo.

[[cahn2dframer]] Programa utilizado para gerar as imagens estáticas vistas neste trabalho. Ele salvará os frames em uma pasta dentro da pasta com nome do código do arquivo, e serão numerados de acordo com qual frame da simulação ela se refere.

[[cahn2danimador]] Programa utilizado para gerar as animações vistas neste trabalho. Os gifs são salvos em uma pasta dentro da pasta com nome do código do arquivo, com o nome dele identificando algumas características dele. Atente que, para computadores sem placa de vídeo dedicada, a renderização dos gifs pode levar muito tempo, em casos como este recomendo gerar as animações em menor resolução ou com intervalos maiores entre as imagens.

[[cahn1dcompare]] Código utilizado para gerar o frame de comparação em uma dimensão. Não recomendo acessá-lo, por não possuir um resultado interessante, além de não estar otimizado nem comentado.

'''As bibliotecas utilizadas que requerem instalação são: Numpy, Scipy e Matplotlib. O renderizador de vídeo FFMPEG foi utilizado, porém na maioria dos computadores modernos executando Windows, na ausência deste, outro renderizador será utilizado sem necessidade de alterar o código. Porém, caso seja necessário instalar o renderizador, este link possui um tutorial simples em inglês para Windows: [https://stackoverflow.com/posts/65860115/revisions]'''

== Referências ==
<references/>