Cahn1ddiflog

De Física Computacional
Revisão de 00h52min de 1 de outubro de 2022 por Leomigotto (discussão | contribs) (Criou página com '<code lang = "python"> import numpy as np import matplotlib.pyplot as plt from scipy.fft import rfft, irfft, rfftfreq import os def cahnexp1d(cc, c1, c2): cc1 = np.copy(cc...')
(dif) ← Edição anterior | Revisão atual (dif) | Versão posterior → (dif)
Ir para navegação Ir para pesquisar

import numpy as np import matplotlib.pyplot as plt from scipy.fft import rfft, irfft, rfftfreq import os def cahnexp1d(cc, c1, c2):

   cc1 = np.copy(cc)*c1
   cc2 = np.copy(cc)*c2
   cc3 = np.copy(cc**3)*c1
   for i in range(l):
           n = i - 2
           cc[n] =   (   2*(cc1[n] - cc3[n]) + (cc3[n-1] + cc3[n+1] -1*(cc1[n-1] + cc1[n+1]))
                           -1*(cc2[n-2] + cc2[n+2] + 6*cc2[n])
                           +(4*cc2[n-1] + 4*cc2[n+1])
                       ) + cc[n]
   return cc

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

np.random.seed(4)

xmin = 0

xmax = 1

tmin = 0

tmax = 0.5

dt = 1/22000000

dx = 1/128

difd = 1

gamma = (3.4/128)**2

l = int((xmax/dx))

k = rfftfreq(l, dx/(2*np.pi))

c1 = difd*dt/(dx**2)

c2 = gamma*c1/(dx**2)

k2 = k**2

k4 = (k2**2)*gamma

ccexp = np.random.rand(l)*2 - 1

ccfourier = np.copy(ccexp)

t = 0

tamanho = l/2

plt.figure(figsize=(8, 8))

plt.yscale("log")

os.makedirs(f".\\1d", exist_ok = True)

dif1, dif2, dif3, dif4 = [], [], [], []

tt = np.arange(dt, tmax + dt, dt) u = 0 while (t < tmax):

   ccfourier = cahnfourier1d(ccfourier, k2, k4)
   ccexp = cahnexp1d(ccexp, c1, c2)
   if u%1000 == 0:
       diferenca = abs(ccexp-ccfourier)
       diferencas = [diferenca[0], diferenca[42], diferenca[75], diferenca[103]]
       dif1.append(diferencas[0])
       dif2.append(diferencas[1])
       dif3.append(diferencas[2])
       dif4.append(diferencas[3])
   t = round(t + dt, int(-np.log10(dt) + 5))
   u+=1

tt = np.linspace(dt, tmax, len(dif1)) plt.title(f"Tempo final: {t:.7f}") plt.plot(tt, dif1) plt.plot(tt, dif2) plt.plot(tt, dif3) plt.plot(tt, dif4) plt.yscale("log") plt.savefig(f".\\1d\\aaaaa{t:.6f}.png", dpi = tamanho)