Equação de Swift-Hohenberg 2D: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
(Criou página com 'Trabalho desenvolvido no semestre de 2022/1 pelo aluno Artur Uhlik Fröhlich para a disciplina de Métodos Computacionais da Física C, ministrada pelo Professor Heitor Carpes...')
 
Sem resumo de edição
Linha 7: Linha 7:
Boa parte dos aspectos teóricos desse trabalho tem como base o livro ''Pattern Formation and Dynamics in Nonequilibrium Systems'' do Cross & Greenside <ref>Cross, M., & Greenside, H. (2009). Pattern Formation and Dynamics in Nonequilibrium Systems. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511627200</ref>.
Boa parte dos aspectos teóricos desse trabalho tem como base o livro ''Pattern Formation and Dynamics in Nonequilibrium Systems'' do Cross & Greenside <ref>Cross, M., & Greenside, H. (2009). Pattern Formation and Dynamics in Nonequilibrium Systems. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511627200</ref>.


Essa equação possui várias versões mas a escolhida para esse trabalho foi a seguinte:
Essa equação pode aparecer nos estudos sobre convexão, equações diferenciais não lineares, formação de padrões e sistemas de não-equilíbrio. Ela possui diferentes versões mas a escolhida para esse trabalho foi a seguinte:


<math> \partial_{t}u = (r-1)u - 2\nabla^{2}u - \nabla^{4}u -u^{3} </math>
<math> \partial_{t}u = (r-1)u - 2\nabla^{2}u - \nabla^{4}u -u^{3} </math>
Linha 13: Linha 13:
.
.


O contexto de derivação dessa  equação foi a convexão de Rayleigh-Bénard e essa formulação tem virtudes mais analíticas do que práticas e apresenta vários aspectos comuns a modelos de formação de padrões e não é considerado uma descrição precisa de nenhum sistema experimental em específico.
O contexto de derivação dessa  equação foi a convexão de [https://en.wikipedia.org/wiki/Rayleigh%E2%80%93B%C3%A9nard_convection Rayleigh-Bénard]. Essa formulação tem virtudes mais analíticas do que práticas e apresenta vários aspectos comuns a modelos de formação de padrões e não é considerado uma descrição precisa de nenhum sistema experimental em específico.


Modificações dessa equação na sua parte não linear e nos coeficientes dos operadores diferenciais geram uma família extensa de EDPs não lineares que então podem servir de base para modelos de sistemas físicos reais.
Modificações dessa equação na sua parte não linear e nos coeficientes dos operadores diferenciais geram uma família extensa de EDPs não lineares que então podem servir de base para modelos de sistemas físicos reais.
Linha 28: Linha 28:


<math> \mathcal{F}\{u(x,t)\}=\hat{u}(k,t) </math>
<math> \mathcal{F}\{u(x,t)\}=\hat{u}(k,t) </math>
<math> \mathcal{F}\left\{\frac{d^{n}u}{dx^{n}}\right\} = (ik)^{n}\hat{u} </math>
<math> \mathcal{F}\left\{\frac{d^{n}u}{dx^{n}}\right\} = (ik)^{n}\hat{u} </math>


Linha 40: Linha 41:
O código é análogo ao que será apresentado posteriormente em duas dimensões então somente os resultados serão mostrados aqui (todos os códigos se encontram no link para o GitHub <ref> Link para o GitHub com todos os códigos: [https://github.com/Artur-UF/MetComp/tree/main/MetCompC/trabalho1 Repositório no GitHub]  </ref>)
O código é análogo ao que será apresentado posteriormente em duas dimensões então somente os resultados serão mostrados aqui (todos os códigos se encontram no link para o GitHub <ref> Link para o GitHub com todos os códigos: [https://github.com/Artur-UF/MetComp/tree/main/MetCompC/trabalho1 Repositório no GitHub]  </ref>)


Para testar o comportamento da equação perante diferentes valores do parâmetro r, teste esse inspirado pela análise de instabilidade linear feita no livro do Cross & Greenside, no qual é encontrado que o valor crítico de r é 0, e que para r<0 o estado tem uma evolução diferente do que para r>0, exemplos disso são as integrações a seguir:
Para testar o comportamento da equação perante diferentes valores do parâmetro r, teste esse inspirado pela análise de instabilidade linear feita no capítulo 2 do livro do Cross & Greenside, na qual é encontrado que o valor crítico de r é 0, e que para r<0 o estado tem uma evolução diferente em comparação com r>0, exemplos disso são as integrações a seguir:


[[Arquivo:plot-1D-SH.png|1000px|thumb|center|Diferentes integrações da Equação de Swift-Hohenberg em uma dimensão usando r=-0.1 (azul) e r=0.1 (vermelho)]]  
[[Arquivo:plot-1D-SH.png|1000px|thumb|center|Diferentes integrações da Equação de Swift-Hohenberg em uma dimensão usando r=-0.1 (azul) e r=0.1 (vermelho)]]  
Linha 103: Linha 104:
     duhaty4 = kpy4 * uhat
     duhaty4 = kpy4 * uhat
     duhatxy = kpx2 * kpy2 * uhat
     duhatxy = kpx2 * kpy2 * uhat
     duhat = (r - 1)*uhat + 2*duhatx2 + 2*duhaty2 - duhatx4  
     duhat = (r - 1)*uhat + 2*duhatx2 + 2*duhaty2 - duhatx4 - duhaty4 - 2*duhatxy - uhat3
            - duhaty4 - 2*duhatxy - uhat3
     dut = ifft2(duhat)
     dut = ifft2(duhat)
     return dut.real
     return dut.real
Linha 142: Linha 142:
Parece uma imagem chata e sem graça mas o que ela nos mostra é que o estado evoluiu para uma solução estável em 0, assim como encontrado em 1 dimensão anteriormente.
Parece uma imagem chata e sem graça mas o que ela nos mostra é que o estado evoluiu para uma solução estável em 0, assim como encontrado em 1 dimensão anteriormente.


Realizando a integração com um valor para r>0 obtemos o seguinte:
Realizando a integração com um valor de r>0 obtemos o seguinte:


[[Arquivo:SH-OFplot.png|500px|thumb|center|Integração de Swift-Hohenberg em 2D utilizando r=0.25 e dt=0.0001]]
[[Arquivo:SH-OFplot.png|500px|thumb|center|Integração de Swift-Hohenberg em 2D utilizando r=0.25 e dt=0.0001]]

Edição das 21h42min de 21 de setembro de 2022

Trabalho desenvolvido no semestre de 2022/1 pelo aluno Artur Uhlik Fröhlich para a disciplina de Métodos Computacionais da Física C, ministrada pelo Professor Heitor Carpes Marques Fernandes.

Introdução

O problema escolhido foi a integração numérica da equação de Swift-Hohenberg em duas dimensões.

Boa parte dos aspectos teóricos desse trabalho tem como base o livro Pattern Formation and Dynamics in Nonequilibrium Systems do Cross & Greenside [1].

Essa equação pode aparecer nos estudos sobre convexão, equações diferenciais não lineares, formação de padrões e sistemas de não-equilíbrio. Ela possui diferentes versões mas a escolhida para esse trabalho foi a seguinte:

.

O contexto de derivação dessa equação foi a convexão de Rayleigh-Bénard. Essa formulação tem virtudes mais analíticas do que práticas e apresenta vários aspectos comuns a modelos de formação de padrões e não é considerado uma descrição precisa de nenhum sistema experimental em específico.

Modificações dessa equação na sua parte não linear e nos coeficientes dos operadores diferenciais geram uma família extensa de EDPs não lineares que então podem servir de base para modelos de sistemas físicos reais.

Integração em 1D

O objetivo final é chegar na Equação de Swift-Hohenberg em duas dimensões mas antes é pertinente verificar a integração em uma dimensão.

Para realizar a integração foram utilizadas as FFTs (Transformadas Rápidas de Fourier) e o método foi baseado nos vídeos e livro do Professor Steven Bruton [2].

O princípio do método é utilizar as transformadas para calcular as derivadas no espaço transformado (que são bem mais simples) e retornar com esses valores para o espaço cartesiano e então com as derivadas calculadas é realizada a evolução temporal utilizando o método de Euler explícito (no caso desse trabalho a integração em 1-D utilizou o integrador da biblioteca Scipy).

Então antes de adentrar no código é preciso fazer a transformada da equação em 1-D, dado que:

Então a equação fica:

Dessa maneira transformamos uma equação diferencial parcial em uma equação diferencial ordinária, pois só possui a derivada temporal.

Com isso é possível escrever um código capaz de fazer a transformada do estado, calcular as derivadas no espaço de Fourier e voltar para o espaço cartesiano e realizar a evolução temporal, tudo isso dentro de um passo da integração.

O código é análogo ao que será apresentado posteriormente em duas dimensões então somente os resultados serão mostrados aqui (todos os códigos se encontram no link para o GitHub [3])

Para testar o comportamento da equação perante diferentes valores do parâmetro r, teste esse inspirado pela análise de instabilidade linear feita no capítulo 2 do livro do Cross & Greenside, na qual é encontrado que o valor crítico de r é 0, e que para r<0 o estado tem uma evolução diferente em comparação com r>0, exemplos disso são as integrações a seguir:

Diferentes integrações da Equação de Swift-Hohenberg em uma dimensão usando r=-0.1 (azul) e r=0.1 (vermelho)

Pode-se notar que para r<0 a solução se torna estável, estabiliza no 0, e para r>0 o estado é instável e se mantém em constante movimento mas com as "ondas" de magnitude constante.

Com esses resultados podemos ir para a integração em duas dimensões.

Para começar a explicar o código devemos discretizar o espaço de Fourier com os seus respectivos números de onda kappa para cada dimensão (já que o campo escalar tem duas dimensões a transformada retorna um campo também de duas dimensões só que no espaço de Fourier).

Os parâmetros N e L determinam o tamanho do campo escalar e a sua definição, testes realizados com esse código retornam um limite inferior de aproximadamente de 0.39 no parâmetro dx.

# Todos os parâmetros estão aqui:
N = 256
L = 100
dx = L/N
r = 0.25
dt = 0.0001
tf = 100
checkpoint = 500
# *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
x = np.arange(-L/2, L/2, dx)
y = np.arange(-L/2, L/2, dx)
size = len(x)

# Os coeficientes
kx = 2 * np.pi * np.fft.fftfreq(N, d=dx)
ky = 2 * np.pi * np.fft.fftfreq(N, d=dx)
kappax, kappay = np.meshgrid(kx, ky)

Em seguida está a função que utiliza a FFT para calcular o lado direito da equação. Com os métodos da biblioteca Scipy o campo escalar inicial será transformado e com ele serão aplicadas as propriedades das transformadas de Fourier para as derivadas da seguinte maneira.

Dado que:

O lado direito da equação transformado então fica assim:

Todo esses termos ao lado direito foram chamados de RHS (Right Hand Side) e calculados na função seguinte.

def rhs(u, kappax, kappay, r):
    '''
    Calcula o lado direito da EDP
    '''
    kpx2 = np.power(kappax, 2)
    kpx4 = np.power(kappax, 4)
    kpy2 = np.power(kappay, 2)
    kpy4 = np.power(kappay, 4)
    uhat = fft2(u)
    uhat3 = fft2(u**3)
    duhatx2 = kpx2 * uhat
    duhaty2 = kpy2 * uhat
    duhatx4 = kpx4 * uhat
    duhaty4 = kpy4 * uhat
    duhatxy = kpx2 * kpy2 * uhat
    duhat = (r - 1)*uhat + 2*duhatx2 + 2*duhaty2 - duhatx4 - duhaty4 - 2*duhatxy - uhat3
    dut = ifft2(duhat)
    return dut.real

Como essa função retorna nosso vetor de estado no espaço cartesiano isso possibilita a realização da evolução temporal da equação utilizando o método de Euler explícito, a partir da criação de um estado inicial (um campo escalar que nesse caso foi inicializado utilizando uma distribuição uniforme de números aleatórios).


# Matriz do estado inicial
u0 = np.random.uniform(-1, 1, (size, size))

# Array da evolução temporal
track = [copy.deepcopy(u0)]

# A integração ocorre aqui
c = 0
t = np.arange(0, tf, dt)
for ti in t:
    c += 1
    u0 += dt * rhs(u0, kappax, kappay, r)
    if c % checkpoint == 0:
        track.append(copy.deepcopy(u0))
    if ti % 2 == 0:
        print(f't = {ti}')

Resultados

Foi realizada a integração com diferentes valores de r para testar as características vistas em uma dimensão, abaixo está o resultado de uma integração usando um r<0:

Integração de Swift-Hohenberg em 2D utilizando r=-0.25

Parece uma imagem chata e sem graça mas o que ela nos mostra é que o estado evoluiu para uma solução estável em 0, assim como encontrado em 1 dimensão anteriormente.

Realizando a integração com um valor de r>0 obtemos o seguinte:

Integração de Swift-Hohenberg em 2D utilizando r=0.25 e dt=0.0001

Podemos ver que para um r>0 o estado não atinge o equilíbrio e evolui constantemente, exemplificando mais uma vez o resultado obtido em 1D.

Outra característica interessante de se notar em 2D se chama "coarsening" que pode ser traduzida como "engrossamento". Esse é o efeito que vemos na evolução do padrão no qual ele vai aumentando as áreas de domínio das suas faixas.

Podemos ver a evolução desses padrões no vídeo da integração gerado por um dos códigos presentes no GitHub.

Outro exemplo desse efeito (coarsenning) também pode ser visto na literatura:

Simulação da equação de Swift-Hohenberg em uma grande geometria periódica de tamanho 256 × 256 com r = 0.25, começando com condições iniciais aleatórias . O aumento do tamanho médio do domínio no painel (a) no tempo 10 para o painel (b) no tempo 10000 é o chamado “coarsening”[4].)

Referências

  1. Cross, M., & Greenside, H. (2009). Pattern Formation and Dynamics in Nonequilibrium Systems. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511627200
  2. Brunton, S., & Kutz, J. (2019). Frontmatter. In Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control (pp. I-Iv). Cambridge: Cambridge University Press.
  3. Link para o GitHub com todos os códigos: Repositório no GitHub
  4. K. R. Elder, Jorge Viñals, and Martin Grant. Dynamic scaling and quasi-ordered states in the 2-dimensional Swift-Hohenberg equation. Phys. Rev. A, 46(12):7618–29, 1992.