Código 1 - Forma Conservativa das Equações de Águas Rasas

De Física Computacional
Ir para: navegação, pesquisa
#%% 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)