190 lines
7.5 KiB
Python
190 lines
7.5 KiB
Python
# this code is based on the wave equation example from this article
|
|
# https://beltoforion.de/en/recreational_mathematics/2d-wave-equation.php
|
|
|
|
import pygame
|
|
import numpy as np
|
|
import time
|
|
|
|
hs = 0.66 # spatial step width
|
|
ts = 0.000001 # time step width in s
|
|
dimx = 300 # width of the simulation domain
|
|
dimy = 300 # height of the simulation domain
|
|
cellsize = 2 # display size of a cell in pixel
|
|
|
|
speed_of_sound = 343000 # mm/s Wave velocity of sound in air
|
|
f = 40000 # frequency
|
|
l = speed_of_sound/f # wavelength
|
|
dist = 10 # disctance between point sources
|
|
count = 10 # number of point sources
|
|
|
|
def create_arrays():
|
|
global velocity
|
|
global tau
|
|
global kappa
|
|
global u
|
|
|
|
# The three dimensional simulation grid
|
|
u = np.zeros((3, dimx, dimy))
|
|
|
|
# A field containing the velocity for each cell
|
|
velocity = np.zeros((dimx, dimy))
|
|
|
|
# A field containing the factor for the Laplace Operator that combines Velocity and Grid Constants for the Wave Equation
|
|
tau = np.zeros((dimx, dimy))
|
|
|
|
# A field containing the factor for the Laplace Operator that combines Velocity and Grid Constants for the Boundary Condition
|
|
kappa = np.zeros((dimx, dimy))
|
|
|
|
|
|
def set_initial_conditions(u):
|
|
global velocity
|
|
global tau
|
|
global kappa
|
|
|
|
velocity[0:dimx,0:dimy] = speed_of_sound
|
|
|
|
# compute tau and kappa from the velocity field
|
|
tau = ( (velocity*ts) / hs )**2
|
|
kappa = ts * velocity / hs
|
|
|
|
|
|
def update(u : any):
|
|
u[2] = u[1]
|
|
u[1] = u[0]
|
|
|
|
# ok, (8th) https://www.ams.org/journals/mcom/1988-51-184/S0025-5718-1988-0935077-0/S0025-5718-1988-0935077-0.pdf; Page 702
|
|
boundary_size = 4
|
|
u[0, 4:dimx-4, 4:dimy-4] = tau[4:dimx-4, 4:dimy-4]\
|
|
* ( - 1/560 * u[1, 4:dimx-4, 0:dimy-8] # c, r-4
|
|
+ 8/315 * u[1, 4:dimx-4, 1:dimy-7] # c, r-3
|
|
- 1/5 * u[1, 4:dimx-4, 2:dimy-6] # c, r-2
|
|
+ 8/5 * u[1, 4:dimx-4, 3:dimy-5] # c, r-1
|
|
|
|
- 1/560 * u[1, 0:dimx-8, 4:dimy-4] # c - 4, r
|
|
+ 8/315 * u[1, 1:dimx-7, 4:dimy-4] # c - 3, r
|
|
- 1/5 * u[1, 2:dimx-6, 4:dimy-4] # c - 2, r
|
|
+ 8/5 * u[1, 3:dimx-5, 4:dimy-4] # c - 1, r
|
|
- 410/72 * u[1, 4:dimx-4, 4:dimy-4] # c , r
|
|
+ 8/5 * u[1, 5:dimx-3, 4:dimy-4] # c + 1, r
|
|
- 1/5 * u[1, 6:dimx-2, 4:dimy-4] # c + 2, r
|
|
+ 8/315 * u[1, 7:dimx-1, 4:dimy-4] # c + 3, r
|
|
- 1/560 * u[1, 8:dimx , 4:dimy-4] # c + 4, r
|
|
|
|
+ 8/5 * u[1, 4:dimx-4, 5:dimy-3] # c , r+1
|
|
- 1/5 * u[1, 4:dimx-4, 6:dimy-2] # c , r+2
|
|
+ 8/315 * u[1, 4:dimx-4, 7:dimy-1] # c , r+3
|
|
- 1/560 * u[1, 4:dimx-4, 8:dimy ] # c , r+4
|
|
) \
|
|
+ 2*u[1, 4:dimx-4, 4:dimy-4] \
|
|
- u[2, 4:dimx-4, 4:dimy-4]
|
|
|
|
# Absorbing Boundary Conditions:
|
|
mur = True
|
|
if mur==True:
|
|
update_boundary(u, boundary_size)
|
|
|
|
|
|
def update_boundary(u, sz) -> None:
|
|
"""Update the boundary cells.
|
|
|
|
Implement MUR boundary conditions. This represents an open boundary were waves can leave the
|
|
simulation domain with little remaining reflection artifacts. Although this is of a low error
|
|
order it is good enough for this simulation.
|
|
"""
|
|
c = dimx-1
|
|
u[0, dimx-sz-1:c, 1:dimy-1] = u[1, dimx-sz-2:c-1, 1:dimy-1] + (kappa[dimx-sz-1:c, 1:dimy-1]-1)/(kappa[ dimx-sz-1:c, 1:dimy-1]+1) * (u[0, dimx-sz-2:c-1,1:dimy-1] - u[1, dimx-sz-1:c,1:dimy-1])
|
|
|
|
c = 0
|
|
u[0, c:sz, 1:dimy-1] = u[1, c+1:sz+1, 1:dimy-1] + (kappa[c:sz, 1:dimy-1]-1)/(kappa[c:sz, 1:dimy-1]+1) * (u[0, c+1:sz+1,1:dimy-1] - u[1,c:sz,1:dimy-1])
|
|
|
|
r = dimy-1
|
|
u[0, 1:dimx-1, dimy-1-sz:r] = u[1, 1:dimx-1, dimy-2-sz:r-1] + (kappa[1:dimx-1, dimy-1-sz:r]-1)/(kappa[1:dimx-1, dimy-1-sz:r]+1) * (u[0, 1:dimx-1, dimy-2-sz:r-1] - u[1, 1:dimx-1, dimy-1-sz:r])
|
|
|
|
r = 0
|
|
u[0, 1:dimx-1, r:sz] = u[1, 1:dimx-1, r+1:sz+1] + (kappa[1:dimx-1, r:sz]-1)/(kappa[1:dimx-1, r:sz]+1) * (u[0, 1:dimx-1, r+1:sz+1] - u[1, 1:dimx-1, r:sz])
|
|
|
|
|
|
def draw_waves(display, u, data, offset):
|
|
global velocity
|
|
global font
|
|
|
|
#data[1:dimx, 1:dimy, 0] = np.clip(u[0, 1:dimx, 1:dimy]*2 + 128, 0, 255)
|
|
#data[1:dimx, 1:dimy, 1] = np.clip(u[1, 1:dimx, 1:dimy]*2 + 128, 0, 255)
|
|
#data[1:dimx, 1:dimy, 2] = np.clip(u[2, 1:dimx, 1:dimy]*2 + 128, 0, 255)
|
|
|
|
data[1:dimx, 1:dimy, 0] = 255 - np.clip((u[0, 1:dimx, 1:dimy]>0) * 10 * u[0, 1:dimx, 1:dimy]+u[1, 1:dimx, 1:dimy]+u[2, 1:dimx, 1:dimy], 0, 255)
|
|
data[1:dimx, 1:dimy, 1] = 255 - np.clip(np.abs(u[0, 1:dimx, 1:dimy]) * 10, 0, 255)
|
|
data[1:dimx, 1:dimy, 2] = 255 - np.clip((u[0, 1:dimx, 1:dimy]<=0) * -10 * u[0, 1:dimx, 1:dimy] + u[1, 1:dimx, 1:dimy] + u[2, 1:dimx, 1:dimy], 0, 255)
|
|
|
|
surf = pygame.surfarray.make_surface(data)
|
|
display.blit(pygame.transform.scale(surf, (dimx * cellsize, dimy * cellsize)), offset)
|
|
|
|
|
|
def draw_text(display, fps, tick):
|
|
global font
|
|
|
|
text_surface = font.render('2D Wave Equation - Explicit Euler (Radiating Boundary Conditions)', True, (0, 0, 40))
|
|
display.blit(text_surface, (5,5))
|
|
|
|
text_surface = font.render(f'FPS: {fps:.1f}; t={tick:.6f} s; area={dimx*hs}x{dimy*hs} mm', True, (0, 0, 40))
|
|
display.blit(text_surface, (5, dimy*cellsize - 20))
|
|
|
|
def calculate_phase_for_focus(source, point):
|
|
return np.sqrt((point[0] + source[0])**2 + (point[1] + source[1])**2) * ((np.pi*2)/l)
|
|
|
|
|
|
def main():
|
|
global font
|
|
|
|
pygame.init()
|
|
pygame.font.init()
|
|
|
|
font = pygame.font.SysFont('Consolas', 15)
|
|
display = pygame.display.set_mode((dimx*cellsize, dimy*cellsize))
|
|
pygame.display.set_caption("Solving the 2d Wave Equation")
|
|
|
|
create_arrays()
|
|
set_initial_conditions(u)
|
|
|
|
image1data = np.zeros((dimx, dimy, 3), dtype = np.uint8 )
|
|
|
|
tick = 0
|
|
last_tick = 0
|
|
fps = 0
|
|
start_time = time.time()
|
|
|
|
while True:
|
|
for event in pygame.event.get():
|
|
if event.type == pygame.QUIT:
|
|
pygame.quit()
|
|
return
|
|
|
|
tick = tick + ts/2
|
|
|
|
current_time = time.time()
|
|
if current_time - start_time > 0.5:
|
|
fps = (tick*1000000-last_tick) / (current_time - start_time)
|
|
start_time = time.time()
|
|
last_tick = tick*1000000
|
|
|
|
x, y = pygame.mouse.get_pos()
|
|
w = (dimx/(2) - (x)/cellsize)*hs
|
|
h = (dimy*2 - ((y*2)/cellsize))*hs
|
|
|
|
x = np.linspace(-(count - 1)*dist*0.5, (count - 1)*dist*0.5, count)
|
|
|
|
for idx, i in enumerate(x):
|
|
phase = calculate_phase_for_focus((x[idx], 0), (w,h))
|
|
u[0:2, int(dimx/(2)+i/hs), int(dimy-10)] = np.sin(tick*2*np.pi*f + phase)*40
|
|
|
|
update(u)
|
|
draw_waves(display, u, image1data, (0,0))
|
|
draw_text(display, fps, tick)
|
|
|
|
#surf1 = pygame.Surface((5,5))
|
|
#display.blit(surf1, (int(dimx) - w*cellsize, int(dimy) - h))
|
|
|
|
pygame.display.update()
|
|
|
|
if __name__ == "__main__":
|
|
main() |