Código 1 - Forma Conservativa das Equações de Águas Rasas
Ir para navegação
Ir para pesquisar
#%% Bibliotecas
import numpy as np
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation
import matplotlib.patches as mpatches
from IPython.display import HTML, Image
#%% Parametros
L_xf = 10.0 # m
L_x0 = -L_xf
NX = 100
dx = (L_xf - L_x0) / NX
L_yf = 10.0 # m
L_y0 = -L_yf
NY = 100
dy = (L_yf - L_y0) / NY
N_INNER = (NX - 2) * (NY - 2)
g = 9.8 # m /s^2
# Tempo
dt = 0.002
Nt = 1500
time_interval = np.arange(0, Nt*dt, dt)
#%% Discretização do espaço x-y
x = np.linspace(L_x0, L_xf, NX-1)
y = np.linspace(L_y0, L_yf, NY-1)
X, Y = np.meshgrid(x, y)
#%% Condições iniciais, em distribuição gaussiana
sigma = 1.0
sigma_v = 1.0
h_2d = 6 * np.ones(shape=np.shape(X)) + np.exp(-(((X)**2 / 2*(sigma**2)) + ((Y)**2 / 2*(sigma**2))))
u_2d = 0.1 * np.exp(-(((X)**2 / 2*(sigma_v**2)) + ((Y)**2 / 2*(sigma_v**2))))
v_2d = 0.1 * np.exp(-(((X)**2 / 2*(sigma_v**2)) + ((Y)**2 / 2*(sigma_v**2))))
b = np.ones(shape=np.shape(X)) # Fundo da superficie constante
h_2d -= b # diminui a altura da superficie
#%% Vetores das variáveis anteriores e historico das variaveis
h_ant = np.copy(h_2d)
v_ant = np.copy(v_2d)
u_ant = np.copy(u_2d)
# Inicilização das listas para armazenar os valores
hist_h, hist_u, hist_v = [], [], []
Função para resolução das equações diferenciais com FTCS.
# Fator de multiplicacao
mx = (dt / (2*(dx)))
my = (dt / (2*(dy)))
def resolve_pdes(h, vx, vy):
"""
Função que retorna os valores de profunidade, velocidade em x e em y
no próximo tempo
Parametros
-----------
h : array
profundidade no tempo t
vx : aray
velocidade em x no tempo t
vy : array
velocidade em y no tempo t
Retorna
-----------
prox_h : array
profundidade no tempo t + dt
prox_u : array
velocidade em x no tempo t + dt
prox_v : array
velocidade em y no tempo t + dt
"""
# Inicializa os vetores para armazenarem os resultados calculados
prox_h = np.ones(shape = (np.shape(h)), dtype = np.float64)
prox_u = np.ones(shape = (np.shape(vx)), dtype = np.float64)
prox_v = np.ones(shape = (np.shape(vy)), dtype = np.float64)
# Loop nos pontos discretizados
for i in range(1, NX - 1):
for j in range(1, NY - 1):
# Alturas e velocidades conforme a posicao do ponto:
# _l : ponto a esquerda, _r: ponto a direita, _up: ponto acima, _d: abaixo
# Condicao a esquerda ------------------
if i == 1: # primeiro x interno
h_l = h[i, j]
u_l = -vx[i, j]
v_l = vy[i, j]
else:
h_l = h[i-1, j]
u_l = vx[i-1, j]
v_l = vy[i-1, j]
# --------------------------------------
# Condicao a direita -------------------
if i == NX - 2: # ultimo x interno
h_r = h[i, j]
u_r = -vx[i, j]
v_r = vy[i, j]
else:
h_r = h[i+1, j]
u_r = vx[i+1, j]
v_r = vy[i+1, j]
# --------------------------------------
# Condicao abaixo ----------------------
if j == 1: # primeiro y interno
h_d = h[i, j]
u_d = vx[i, j]
v_d = - vy[i, j]
else:
h_d = h[i, j - 1]
v_d = vy[i, j - 1]
u_d = vx[i, j - 1]
# --------------------------------------
# Condicao acima ----------------------
if j == NY - 2: # utlimo y interno
h_up = h[i, j]
u_up = vx[i, j]
v_up = - vy[i, j]
else:
h_up = h[i, j + 1]
v_up = vy[i, j + 1]
u_up = vx[i, j + 1]
# --------------------------------------
## Primeira Equação
h_ij = h[i, j] - \
(h_r * u_r - h_l * u_l) * fator_x - \
(h_up * v_up - h_d * v_d) * fator_y
prox_h[i, j] = h_ij
# ## Segunda equação
hu_ij = ((h[i, j]) * vx[i, j]) - \
mx * (
((h_r * (u_r ** 2)) + (1/2 * g * (h_r ** 2))) -\
((h_l * (u_l ** 2)) + (1/2 * g * (h_l ** 2)))
) - \
my * (
(h_up * u_up * v_up ) - (h_d * u_d * v_d)
) - \
mx * g * h[i, j] * (b[i, j+1] - b[i, j-1])
# # ## Terceira Equação
hv_ij = (h[i, j] * vy[i, j]) - \
m_x * (
(h_r * u_r * v_r) - (h_l * u_l * v_l)
) - \
my * (
((h_up * (v_up ** 2)) + (1/2 * g * (h_up ** 2))) - \
((h_d * (v_d ** 2)) + (1/2 * g * (h_d ** 2)))
) - \
my * g * h[i, j] * (b[i, j+1] - b[i, j-1])
prox_v[i, j] = hv_ij / (h_ij)
return prox_h, prox_u, prox_v
#%% Cálculo
# Resolve para cada passo de tempo
for t in time_interval:
h, u, v = resolve_pdes(h_ant, u_ant, v_ant) # valores das variaveis no tempo = t + dt
# faz isso pq tava dando um erro estranho nesses pontos?
h[0, :] = h[1, :]
h[:, 0] = h[:, 1]
v[0, :] = v[1, :]
v[:, 0] = v[:, 1]
u[0, :] = u[1, :]
u[:, 0] = u[:, 1]
# adicionar essas variáveis em listas pra conseguir plotar dps
hist_h.append(h); hist_u.append(u); hist_v.append(v)
# Coloca as variaveis atuais como anteriores pro proximo calculo
h_ant = np.copy(h)
u_ant = np.copy(u)
v_ant = np.copy(v)
#%% Gráfico
# Reorganiza os vetores para plotar
x_2d = X[0]
y_2d = Y[0]
for i in range(1,len(X)):
x_2d = np.append(x_2d, X[i])
y_2d = np.append(y_2d, Y[i])
def animate(i):
plt.clf() # limpa a figura, pra nao ficar sobrepondo figs
# titulos
plt.suptitle('Evolução da onda', fontsize=14)
plt.title(f'Tempo: {round(dt*i, 3)}', fontsize=12)
# labels
plt.ylabel('y', fontsize=8)
plt.xlabel('x', fontsize=8)
# plot
graph = plt.scatter(x_2d, y_2d, c= hist_h[0::8][i], vmin=4.5, vmax=5.5, marker='.')
plt.colorbar()
plt.savefig(f'caminho-frames/{i}.png')
return graph
# fig = plt.figure(figsize=(16,9))
# ax = fig.gca(projection='3d')
fig, ax = plt.subplots()
ani = animation.FuncAnimation(fig, animate, frames = len(hist_h[0::8]), repeat=False, interval=0.1)
#%% Salva o gif, com os frames salvos anteriormente
images = []
for i in range(len(hist_h[0::80])):
images.append(imageio.imread(f'caminho-frames/{i}.png'))
imageio.mimsave("gif2.gif", images,fps=8)