-
Notifications
You must be signed in to change notification settings - Fork 6
/
FlowSimulation_Pseudospectral_VortexPair.py
64 lines (49 loc) · 1.27 KB
/
FlowSimulation_Pseudospectral_VortexPair.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Dec 2 16:45:29 2018
@author: jiahan
"""
import scipy as sc
import matplotlib.pyplot as plt
import IncompNST_Tools.Pseudospectral as nst
from scipy.fftpack import fft2, ifft2
# Physical constants
nu = 0.001
# Length of field
lx = 1
ly = 1
# Number of grid points
nx = 512
ny = 512
# Grid increments
dx = lx/nx
dy = ly/ny
# Initialize vortex
omega, p = nst.vortex_pair(nx, ny, dx, dy)
# Gradient operators in Fourier domain for x- and y-direction
Kx, Ky = nst.Spectral_Gradient(nx, ny, lx, ly)
# 2D Laplace operator and 2D inverse Laplace operator in Fourier domain
K2, K2inv = nst.Spectral_Laplace(nx, ny, Kx, Ky)
# Simulation time
T_simu = 10000
# Set discrete time step by choosing CFL number (condition: CFL <= 1)
CFL = 1
u = sc.real(ifft2(-Ky*K2inv*fft2(omega)))
v = sc.real(ifft2(Kx*K2inv*fft2(omega)))
u_max = sc.amax(sc.absolute(u))
v_max = sc.amax(sc.absolute(v))
t_step = (CFL*dx*dy)/(u_max*dy+v_max*dx)
# Start Simulation
t_sum = 0
i = 0
while t_sum <= T_simu:
# Runge-Kitta 4 time simulation
omega = nst.RK4(t_step, omega, Kx, Ky, K2, K2inv, nu)
# Plot every 100th frame
if 0 == i % 100:
plt.imshow(omega)
plt.pause(0.1)
i += 1
t_sum += t_step
print(i, t_sum)