VoxelEngine/FluidSim/LatticeBoltzmann.py

353 lines
12 KiB
Python
Raw Normal View History

2021-08-03 10:33:07 +02:00
import matplotlib.pyplot as plt
import numpy as np
2021-10-15 19:21:17 +02:00
from FluidSim.FluidSimParameters import *
2021-08-03 10:33:07 +02:00
"""
Create Your Own Lattice Boltzmann Simulation (With Python)
Philip Mocz (2020) Princeton Univeristy, @PMocz
Simulate flow past cylinder
for an isothermal fluid
"""
def main():
""" Finite Volume simulation """
# Simulation parameters
2021-10-15 19:21:17 +02:00
epsilon = 0.000000001
2021-08-03 10:33:07 +02:00
Nx = 400 # resolution x-dir
Ny = 100 # resolution y-dir
2021-10-15 19:21:17 +02:00
rho0 = 1 # average density
2021-08-03 10:33:07 +02:00
tau = 0.6 # collision timescale
Nt = 80000 # number of timesteps
plotRealTime = True # switch on for plotting as the simulation goes along
2022-01-02 11:33:28 +01:00
render_frequency = 10
close_up_frequency = 10
friction = 0.0001
2021-08-03 10:33:07 +02:00
2021-10-15 19:21:17 +02:00
params = FluidSimParameter(Ny)
# params = WaterParameter(Ny)
# params = MagmaParameter(Ny)
2021-08-03 10:33:07 +02:00
# Lattice speeds / weights
NL = 9
idxs = np.arange(NL)
cxs = np.array([0, 0, 1, 1, 1, 0, -1, -1, -1])
cys = np.array([0, 1, 1, 0, -1, -1, -1, 0, 1])
weights = np.array([4 / 9, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36]) # sums to 1
2021-10-15 19:21:17 +02:00
xx, yy = np.meshgrid(range(Nx), range(Ny))
2021-08-03 10:33:07 +02:00
# Initial Conditions
2021-10-15 19:21:17 +02:00
N = np.ones((Ny, Nx, NL)) # * rho0 / NL
temperature = np.ones((Ny, Nx, NL), np.float) # * rho0 / NL
2021-12-20 17:21:46 +01:00
tracked_fluid = np.zeros((Ny, Nx, NL), np.float) # * rho0 / NL
tracked_fluid[50:, :, 0] = 1
2021-08-03 10:33:07 +02:00
has_fluid = np.ones((Ny, Nx), dtype=np.bool)
has_fluid[int(Ny/2):, :] = False
np.random.seed(42)
2021-10-15 19:21:17 +02:00
N += 0.01 * np.random.randn(Ny, Nx, NL)
2021-08-03 10:33:07 +02:00
X, Y = np.meshgrid(range(Nx), range(Ny))
2021-10-15 19:21:17 +02:00
N[:, :, 3] += 2 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4))
# N[:, :, 5] += 1
rho = np.sum(N, 2)
temperature_rho = np.sum(temperature, 2)
2021-08-03 10:33:07 +02:00
for i in idxs:
2021-10-15 19:21:17 +02:00
N[:, :, i] *= rho0 / rho
temperature[:, :, i] *= 1 / temperature_rho
# N[50:, :] = 0
temperature[:, :] = 0
# temperature += 0.01 * np.random.randn(Ny, Nx, NL)
2021-08-03 10:33:07 +02:00
# Cylinder boundary
X, Y = np.meshgrid(range(Nx), range(Ny))
2021-12-20 17:21:46 +01:00
2021-08-03 10:33:07 +02:00
cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4) ** 2
2021-12-20 17:21:46 +01:00
cylinder[:, :] = False
2021-10-15 19:21:17 +02:00
N[cylinder] = 0
N[0, :] = 0
N[Ny - 1, :] = 0
temperature[cylinder] = 0
2021-12-20 17:21:46 +01:00
tracked_fluid[cylinder] = 0
2021-10-15 19:21:17 +02:00
# N[int(Ny/2):, :] = 0
2021-08-03 10:33:07 +02:00
has_fluid[cylinder] = False
has_fluid[0, :] = False
has_fluid[Ny - 1, :] = False
# for i in idxs:
2021-10-15 19:21:17 +02:00
# N[:, :, i] *= has_fluid
2021-08-03 10:33:07 +02:00
# Prep figure
fig = plt.figure(figsize=(4, 2), dpi=80)
reflection_mapping = [0, 5, 6, 7, 8, 1, 2, 3, 4]
# Simulation Main Loop
for it in range(Nt):
print(it)
# Drift
new_has_fluid = np.zeros((Ny, Nx))
2021-10-15 19:21:17 +02:00
F_sum = np.sum(N, 2)
2021-08-03 10:33:07 +02:00
for i, cx, cy in zip(idxs, cxs, cys):
2021-10-15 19:21:17 +02:00
F_part = N[:, :, i] / F_sum
2021-08-03 10:33:07 +02:00
F_part[F_sum == 0] = 0
to_move = F_part * (has_fluid * 1.0)
to_move = (np.roll(to_move, cx, axis=1))
to_move = (np.roll(to_move, cy, axis=0))
new_has_fluid += to_move
2021-10-15 19:21:17 +02:00
N[:, :, i] = np.roll(N[:, :, i], cx, axis=1)
N[:, :, i] = np.roll(N[:, :, i], cy, axis=0)
temperature[:, :, i] = np.roll(temperature[:, :, i], cx, axis=1)
temperature[:, :, i] = np.roll(temperature[:, :, i], cy, axis=0)
2021-08-03 10:33:07 +02:00
2021-12-20 17:21:46 +01:00
tracked_fluid[:, :, i] = np.roll(tracked_fluid[:, :, i], cx, axis=1)
tracked_fluid[:, :, i] = np.roll(tracked_fluid[:, :, i], cy, axis=0)
2021-08-03 10:33:07 +02:00
# has_fluid = new_has_fluid > 0.5
# new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0
# new_has_fluid[(np.abs(F_sum) < 0.000000001)] = 0
2021-12-20 17:21:46 +01:00
#
# fluid_sum = np.sum(has_fluid * 1.0)
# has_fluid = (new_has_fluid / np.sum(new_has_fluid * 1.0)) * fluid_sum
#
# print('fluid_cells: %d' % np.sum(has_fluid * 1))
2021-08-03 10:33:07 +02:00
# for i in idxs:
2021-10-15 19:21:17 +02:00
# N[:, :, i] *= has_fluid
2021-08-03 10:33:07 +02:00
bndry = np.zeros((Ny, Nx), dtype=np.bool)
bndry[0, :] = True
bndry[Ny - 1, :] = True
# bndry[:, 0] = True
# bndry[:, Nx - 1] = True
bndry = np.logical_or(bndry, cylinder)
# bndry = np.logical_or(bndry, has_fluid < 0.5)
# Set reflective boundaries
2021-10-15 19:21:17 +02:00
bndryN = N[bndry, :]
bndryN = bndryN[:, reflection_mapping]
bndryTemp = temperature[bndry, :]
bndryTemp = bndryTemp[:, reflection_mapping]
2021-08-03 10:33:07 +02:00
2021-12-20 17:21:46 +01:00
bndryTracked = tracked_fluid[bndry, :]
bndryTracked = bndryTracked[:, reflection_mapping]
2021-10-15 19:21:17 +02:00
sum_f = np.sum(N)
print('Sum of Particles: %f' % sum_f)
print('Sum of Temperature: %f' % np.sum(temperature))
2021-12-20 17:21:46 +01:00
print('Sum of tacked particles: %f' % np.sum(tracked_fluid))
2021-08-03 10:33:07 +02:00
2021-10-15 19:21:17 +02:00
# sum_f_cyl = np.sum(N[cylinder])
2021-08-03 10:33:07 +02:00
# print('Sum of Forces in cylinder: %f' % sum_f_cyl)
2021-10-15 19:21:17 +02:00
# sum_f_inner_cyl = np.sum(N[inner_cylinder])
2021-08-03 10:33:07 +02:00
# print('Sum of Forces in inner cylinder: %f' % sum_f_inner_cyl)
# if sum_f > 4000000.000000:
# test = 1
2021-10-15 19:21:17 +02:00
# N[Ny - 1, :, 5] += 0.1
# N[0, :, 1] -= 0.1
# N[0, :, 5] += 0.1
# N[Ny - 1, :, 1] -= 0.1
2021-08-03 10:33:07 +02:00
# Calculate fluid variables
2021-10-15 19:21:17 +02:00
rho = np.sum(N, 2)
temp_rho = np.sum(temperature, 2)
2021-12-20 17:21:46 +01:00
tracked_rho = np.sum(tracked_fluid, 2)
2021-10-15 19:21:17 +02:00
ux = np.sum(N * cxs, 2) / rho
uy = np.sum(N * cys, 2) / rho
ux[(np.abs(rho) < epsilon)] = 0
uy[(np.abs(rho) < epsilon)] = 0
g = -params.g * (temp_rho - yy / Ny)
2021-12-20 17:21:46 +01:00
uy1 = (uy + g * params.t1 * (tracked_rho * 0.9 + 0.1))
uy2 = (uy + g * params.t2 * (tracked_rho * 0.9 + 0.1))
2021-10-15 19:21:17 +02:00
# uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0
2021-12-20 17:21:46 +01:00
# uy += g / 2.0
2021-10-15 19:21:17 +02:00
# u_length = np.maximum(np.abs(ux), np.abs(uy))
2022-01-02 11:33:28 +01:00
# safe guard against supersonic streams WIP
2021-12-20 17:21:46 +01:00
u_length1 = np.sqrt(np.square(ux) + np.square(uy1))
u_length2 = np.sqrt(np.square(ux) + np.square(uy2))
u_max_length1 = np.max(u_length1)
u_max_length2 = np.max(u_length2)
u_max_length = max(u_max_length1, u_max_length2)
if u_max_length > 2:
test = 1
2021-10-15 19:21:17 +02:00
if u_max_length > np.sqrt(2):
ux = (ux / u_max_length) * np.sqrt(2)
2021-12-20 17:21:46 +01:00
uy1 = (uy1 / u_max_length) * np.sqrt(2)
uy2 = (uy2 / u_max_length) * np.sqrt(2)
2021-10-15 19:21:17 +02:00
2022-01-02 11:33:28 +01:00
# apply friction
ux *= (1 - friction)
uy1 *= (1 - friction)
uy2 *= (1 - friction)
2021-10-15 19:21:17 +02:00
print('max vector part: %f' % u_max_length)
# ux /= u_max_length
# uy /= u_max_length
# scale = abs(np.max(np.maximum(np.abs(ux), np.abs(uy))) - 1.0) < epsilon
# if scale:
# g = 0.01 * (temp_rho - yy / Ny)
#
# # F = np.zeros((Ny, Nx), dtype=np.bool)
# # F = -0.1 * rho
#
# # uy[np.abs(rho) >= epsilon] += tau * F[np.abs(rho) >= epsilon] / rho[np.abs(rho) >= epsilon]
# uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0
#
# u_length = np.maximum(np.abs(ux), np.abs(uy))
# u_max_length = np.max(u_length)
#
# print('max vector part: %f' % u_max_length)
#
# ux /= u_max_length
# uy /= u_max_length
2021-08-03 10:33:07 +02:00
# print('minimum rho: %f' % np.min(np.abs(rho)))
2021-10-15 19:21:17 +02:00
# print('Maximum N: %f' % np.max(N))
# print('Minimum N: %f' % np.min(N))
2021-08-03 10:33:07 +02:00
# Apply Collision
2021-10-15 19:21:17 +02:00
temperature_eq = np.zeros(temperature.shape)
2021-12-20 17:21:46 +01:00
tracked_eq = np.zeros(temperature.shape)
2021-10-15 19:21:17 +02:00
Neq = np.zeros(N.shape)
2021-08-03 10:33:07 +02:00
for i, cx, cy, w in zip(idxs, cxs, cys, weights):
2021-10-15 19:21:17 +02:00
Neq[:, :, i] = rho * w * (
2021-12-20 17:21:46 +01:00
1 + 3 * (cx * ux + cy * uy1) + 9 * (cx * ux + cy * uy1) ** 2 / 2 - 3 * (ux ** 2 + uy1 ** 2) / 2)
2021-08-03 10:33:07 +02:00
2021-10-15 19:21:17 +02:00
temperature_eq[:, :, i] = temp_rho * w * (
2021-12-20 17:21:46 +01:00
1 + 3 * (cx * ux + cy * uy2) + 9 * (cx * ux + cy * uy2) ** 2 / 2 - 3 * (ux ** 2 + uy2 ** 2) / 2)
tracked_eq[:, :, i] = tracked_rho * w * (
1 + 3 * (cx * ux + cy * uy1) + 9 * (cx * ux + cy * uy1) ** 2 / 2 - 3 * (ux ** 2 + uy1 ** 2) / 2)
2021-10-15 19:21:17 +02:00
# test1 = np.sum(Neq)
test2 = np.sum(N-Neq)
if abs(test2) > 0.0001:
test = ''
print('Overall change: %f' % test2)
n_pre_sum = np.sum(N[np.logical_not(bndry)])
temperature_pre_sum = np.sum(temperature[np.logical_not(bndry)])
N += -(1.0 / params.t1) * (N - Neq)
temperature += -(1.0 / params.t2) * (temperature - temperature_eq)
2021-12-20 17:21:46 +01:00
tracked_fluid += -(1.0 / params.t1) * (tracked_fluid - tracked_eq)
2021-08-03 10:33:07 +02:00
# Apply boundary
2021-10-15 19:21:17 +02:00
N[bndry, :] = bndryN
temperature[bndry, :] = bndryTemp
2021-12-20 17:21:46 +01:00
tracked_fluid[bndry, :] = bndryTracked
2021-10-15 19:21:17 +02:00
# temperature[0, :, 0] = 0
# temperature[1, :, 0] = 0
temperature[0, :, 0] /= 2
temperature[1, :, 0] /= 2
2021-12-20 17:21:46 +01:00
if it <= 3000:
temperature[Ny - 1, :, 0] = 1
temperature[Ny - 2, :, 0] = 1
2021-10-15 19:21:17 +02:00
# n_sum = np.sum(N, 2)
# n_sum_min = np.min(n_sum)
# if n_sum_min < 0:
# N[np.logical_not(bndry)] += abs(n_sum_min)
# N[np.logical_not(bndry)] /= np.sum(N[np.logical_not(bndry)])
# N[np.logical_not(bndry)] *= n_pre_sum
# print('Sum of Forces: %f' % np.sum(N))
# temperature_sum = np.sum(temperature, 2)
# temperature_sum_min = np.min(temperature_sum)
# if temperature_sum_min < 0:
# temperature[np.logical_not(bndry)] += abs(temperature_sum_min)
# temperature[np.logical_not(bndry)] /= np.sum(temperature[np.logical_not(bndry)])
# temperature[np.logical_not(bndry)] *= temperature_pre_sum
# print('Sum of Temperature: %f' % np.sum(temperature))
no_cylinder_mask = np.sum(N, 2) != 0
print('min N: %f' % np.min(np.sum(N, 2)[no_cylinder_mask]))
print('max N: %f' % np.max(np.sum(N, 2)))
print('min Temp: %f' % np.min(np.sum(temperature, 2)[no_cylinder_mask]))
print('max Temp: %f' % np.max(np.sum(temperature, 2)))
2021-08-03 10:33:07 +02:00
2022-01-02 11:33:28 +01:00
if it > render_frequency:
render_frequency = close_up_frequency
2021-08-03 10:33:07 +02:00
# plot in real time - color 1/2 particles blue, other half red
2022-01-02 11:33:28 +01:00
if (plotRealTime and (it % render_frequency) == 0) or (it == Nt - 1):
2021-10-15 19:21:17 +02:00
fig.clear()
2021-08-03 10:33:07 +02:00
plt.cla()
ux[cylinder] = 0
uy[cylinder] = 0
vorticity = (np.roll(ux, -1, axis=0) - np.roll(ux, 1, axis=0)) - (
np.roll(uy, -1, axis=1) - np.roll(uy, 1, axis=1))
vorticity[cylinder] = np.nan
# vorticity *= has_fluid
cmap = plt.cm.bwr
cmap.set_bad('black')
2021-10-15 19:21:17 +02:00
plt.subplot(2, 2, 1)
plt.imshow(vorticity, cmap='bwr')
plt.clim(-.1, .1)
# plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr')
2021-08-03 10:33:07 +02:00
# plt.imshow(bndry * 2.0 - 1.0, cmap='bwr')
2021-10-15 19:21:17 +02:00
# plt.imshow(np.sum(N, 2) * 2.0 - 1.0, cmap='bwr')
# plt.imshow((np.sum(temperature, 2) / np.max(np.sum(temperature, 2))) * 2.0 - 1.0, cmap='bwr')
plt.subplot(2, 2, 2)
max_temp = np.max(np.sum(temperature, 2))
# plt.imshow(np.sum(temperature, 2) / max_temp * 2.0 - 1.0, cmap='bwr')
plt.imshow(np.sum(temperature, 2) * 2.0 - 1.0, cmap='bwr')
plt.clim(-.1, .1)
plt.subplot(2, 2, 3)
max_N = np.max(np.sum(N, 2))
plt.imshow(np.sum(N, 2) / max_N * 2.0 - 1.0, cmap='bwr')
2021-12-20 17:21:46 +01:00
plt.clim(-.1, .1)
2021-08-03 10:33:07 +02:00
2021-12-20 17:21:46 +01:00
plt.subplot(2, 2, 4)
plt.imshow(np.sum(tracked_fluid, 2) * 2.0 - 1.0, cmap='bwr')
2021-08-03 10:33:07 +02:00
plt.clim(-.1, .1)
2021-12-20 17:21:46 +01:00
2021-10-15 19:21:17 +02:00
# ax = plt.gca()
# ax.invert_yaxis()
# ax.get_xaxis().set_visible(False)
# ax.get_yaxis().set_visible(False)
# ax.set_aspect('equal')
2021-08-03 10:33:07 +02:00
plt.pause(0.001)
# Save figure
# plt.savefig('latticeboltzmann.png', dpi=240)
plt.show()
return 0
if __name__ == "__main__":
main()