lava lamp lattice boltzmann
This commit is contained in:
zomseffen
parent
95bd3de45a
commit
8a5de47da3
2 changed files with 188 additions and 45 deletions
24
FluidSim/FluidSimParameters.py
Normal file
24
FluidSim/FluidSimParameters.py
Normal file
|
|
@ -0,0 +1,24 @@
|
||||||
|
class FluidSimParameter:
|
||||||
|
viscosity = 0.1 / 3.0
|
||||||
|
# Pr = 1.0
|
||||||
|
Pr = 100.0
|
||||||
|
# vc = 1.0
|
||||||
|
vc = 0.5
|
||||||
|
|
||||||
|
def __init__(self, height: int):
|
||||||
|
self.t1 = 3 * self.viscosity + 0.5
|
||||||
|
self.t2 = (2 * self.t1 - 1) / (2 * self.Pr) + 0.5
|
||||||
|
self.g = (self.vc ** 2) / height
|
||||||
|
|
||||||
|
self.R = self.Pr * self.g * (height ** 3) / (self.viscosity ** 2)
|
||||||
|
|
||||||
|
|
||||||
|
class MagmaParameter(FluidSimParameter):
|
||||||
|
viscosity = 10 ** 19
|
||||||
|
Pr = 10 ** 25
|
||||||
|
|
||||||
|
|
||||||
|
class WaterParameter(FluidSimParameter):
|
||||||
|
viscosity = 8.9 * 10 ** -4
|
||||||
|
Pr = 7.56
|
||||||
|
vc = 0.05
|
||||||
|
|
@ -1,5 +1,6 @@
|
||||||
import matplotlib.pyplot as plt
|
import matplotlib.pyplot as plt
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
from FluidSim.FluidSimParameters import *
|
||||||
|
|
||||||
"""
|
"""
|
||||||
Create Your Own Lattice Boltzmann Simulation (With Python)
|
Create Your Own Lattice Boltzmann Simulation (With Python)
|
||||||
|
|
@ -13,48 +14,63 @@ def main():
|
||||||
""" Finite Volume simulation """
|
""" Finite Volume simulation """
|
||||||
|
|
||||||
# Simulation parameters
|
# Simulation parameters
|
||||||
|
epsilon = 0.000000001
|
||||||
Nx = 400 # resolution x-dir
|
Nx = 400 # resolution x-dir
|
||||||
Ny = 100 # resolution y-dir
|
Ny = 100 # resolution y-dir
|
||||||
rho0 = 100 # average density
|
rho0 = 1 # average density
|
||||||
tau = 0.6 # collision timescale
|
tau = 0.6 # collision timescale
|
||||||
Nt = 80000 # number of timesteps
|
Nt = 80000 # number of timesteps
|
||||||
plotRealTime = True # switch on for plotting as the simulation goes along
|
plotRealTime = True # switch on for plotting as the simulation goes along
|
||||||
|
|
||||||
|
params = FluidSimParameter(Ny)
|
||||||
|
# params = WaterParameter(Ny)
|
||||||
|
# params = MagmaParameter(Ny)
|
||||||
|
|
||||||
# Lattice speeds / weights
|
# Lattice speeds / weights
|
||||||
NL = 9
|
NL = 9
|
||||||
idxs = np.arange(NL)
|
idxs = np.arange(NL)
|
||||||
cxs = np.array([0, 0, 1, 1, 1, 0, -1, -1, -1])
|
cxs = np.array([0, 0, 1, 1, 1, 0, -1, -1, -1])
|
||||||
cys = np.array([0, 1, 1, 0, -1, -1, -1, 0, 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
|
weights = np.array([4 / 9, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36, 1 / 9, 1 / 36]) # sums to 1
|
||||||
|
xx, yy = np.meshgrid(range(Nx), range(Ny))
|
||||||
|
|
||||||
# Initial Conditions
|
# Initial Conditions
|
||||||
F = np.ones((Ny, Nx, NL)) # * rho0 / NL
|
N = np.ones((Ny, Nx, NL)) # * rho0 / NL
|
||||||
|
temperature = np.ones((Ny, Nx, NL), np.float) # * rho0 / NL
|
||||||
has_fluid = np.ones((Ny, Nx), dtype=np.bool)
|
has_fluid = np.ones((Ny, Nx), dtype=np.bool)
|
||||||
has_fluid[int(Ny/2):, :] = False
|
has_fluid[int(Ny/2):, :] = False
|
||||||
np.random.seed(42)
|
np.random.seed(42)
|
||||||
F += 0.01 * np.random.randn(Ny, Nx, NL)
|
N += 0.01 * np.random.randn(Ny, Nx, NL)
|
||||||
X, Y = np.meshgrid(range(Nx), range(Ny))
|
X, Y = np.meshgrid(range(Nx), range(Ny))
|
||||||
F[:, :, 3] += 2 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4))
|
N[:, :, 3] += 2 * (1 + 0.2 * np.cos(2 * np.pi * X / Nx * 4))
|
||||||
# F[:, :, 5] += 1
|
# N[:, :, 5] += 1
|
||||||
rho = np.sum(F, 2)
|
rho = np.sum(N, 2)
|
||||||
|
temperature_rho = np.sum(temperature, 2)
|
||||||
for i in idxs:
|
for i in idxs:
|
||||||
F[:, :, i] *= rho0 / rho
|
N[:, :, i] *= rho0 / rho
|
||||||
|
temperature[:, :, i] *= 1 / temperature_rho
|
||||||
|
# N[50:, :] = 0
|
||||||
|
temperature[:, :] = 0
|
||||||
|
# temperature += 0.01 * np.random.randn(Ny, Nx, NL)
|
||||||
|
|
||||||
|
|
||||||
# Cylinder boundary
|
# Cylinder boundary
|
||||||
X, Y = np.meshgrid(range(Nx), range(Ny))
|
X, Y = np.meshgrid(range(Nx), range(Ny))
|
||||||
cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4) ** 2
|
cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4) ** 2
|
||||||
inner_cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4 - 2) ** 2
|
inner_cylinder = (X - Nx / 4) ** 2 + (Y - Ny / 2) ** 2 < (Ny / 4 - 2) ** 2
|
||||||
F[cylinder] = 0
|
N[cylinder] = 0
|
||||||
F[0, :] = 0
|
N[0, :] = 0
|
||||||
F[Ny - 1, :] = 0
|
N[Ny - 1, :] = 0
|
||||||
# F[int(Ny/2):, :] = 0
|
|
||||||
|
temperature[cylinder] = 0
|
||||||
|
# N[int(Ny/2):, :] = 0
|
||||||
|
|
||||||
has_fluid[cylinder] = False
|
has_fluid[cylinder] = False
|
||||||
has_fluid[0, :] = False
|
has_fluid[0, :] = False
|
||||||
has_fluid[Ny - 1, :] = False
|
has_fluid[Ny - 1, :] = False
|
||||||
|
|
||||||
# for i in idxs:
|
# for i in idxs:
|
||||||
# F[:, :, i] *= has_fluid
|
# N[:, :, i] *= has_fluid
|
||||||
|
|
||||||
# Prep figure
|
# Prep figure
|
||||||
fig = plt.figure(figsize=(4, 2), dpi=80)
|
fig = plt.figure(figsize=(4, 2), dpi=80)
|
||||||
|
|
@ -66,9 +82,9 @@ def main():
|
||||||
|
|
||||||
# Drift
|
# Drift
|
||||||
new_has_fluid = np.zeros((Ny, Nx))
|
new_has_fluid = np.zeros((Ny, Nx))
|
||||||
F_sum = np.sum(F, 2)
|
F_sum = np.sum(N, 2)
|
||||||
for i, cx, cy in zip(idxs, cxs, cys):
|
for i, cx, cy in zip(idxs, cxs, cys):
|
||||||
F_part = F[:, :, i] / F_sum
|
F_part = N[:, :, i] / F_sum
|
||||||
F_part[F_sum == 0] = 0
|
F_part[F_sum == 0] = 0
|
||||||
to_move = F_part * (has_fluid * 1.0)
|
to_move = F_part * (has_fluid * 1.0)
|
||||||
to_move = (np.roll(to_move, cx, axis=1))
|
to_move = (np.roll(to_move, cx, axis=1))
|
||||||
|
|
@ -76,8 +92,11 @@ def main():
|
||||||
|
|
||||||
new_has_fluid += to_move
|
new_has_fluid += to_move
|
||||||
|
|
||||||
F[:, :, i] = np.roll(F[:, :, i], cx, axis=1)
|
N[:, :, i] = np.roll(N[:, :, i], cx, axis=1)
|
||||||
F[:, :, i] = np.roll(F[:, :, i], cy, axis=0)
|
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)
|
||||||
|
|
||||||
# has_fluid = new_has_fluid > 0.5
|
# has_fluid = new_has_fluid > 0.5
|
||||||
# new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0
|
# new_has_fluid[F_sum == 0] += has_fluid[F_sum == 0] * 1.0
|
||||||
|
|
@ -89,7 +108,7 @@ def main():
|
||||||
print('fluid_cells: %d' % np.sum(has_fluid * 1))
|
print('fluid_cells: %d' % np.sum(has_fluid * 1))
|
||||||
|
|
||||||
# for i in idxs:
|
# for i in idxs:
|
||||||
# F[:, :, i] *= has_fluid
|
# N[:, :, i] *= has_fluid
|
||||||
|
|
||||||
bndry = np.zeros((Ny, Nx), dtype=np.bool)
|
bndry = np.zeros((Ny, Nx), dtype=np.bool)
|
||||||
bndry[0, :] = True
|
bndry[0, :] = True
|
||||||
|
|
@ -101,51 +120,138 @@ def main():
|
||||||
# bndry = np.logical_or(bndry, has_fluid < 0.5)
|
# bndry = np.logical_or(bndry, has_fluid < 0.5)
|
||||||
|
|
||||||
# Set reflective boundaries
|
# Set reflective boundaries
|
||||||
bndryF = F[bndry, :]
|
bndryN = N[bndry, :]
|
||||||
bndryF = bndryF[:, reflection_mapping]
|
bndryN = bndryN[:, reflection_mapping]
|
||||||
|
|
||||||
sum_f = np.sum(F)
|
bndryTemp = temperature[bndry, :]
|
||||||
print('Sum of Forces: %f' % sum_f)
|
bndryTemp = bndryTemp[:, reflection_mapping]
|
||||||
|
|
||||||
# sum_f_cyl = np.sum(F[cylinder])
|
sum_f = np.sum(N)
|
||||||
|
print('Sum of Particles: %f' % sum_f)
|
||||||
|
print('Sum of Temperature: %f' % np.sum(temperature))
|
||||||
|
|
||||||
|
# sum_f_cyl = np.sum(N[cylinder])
|
||||||
# print('Sum of Forces in cylinder: %f' % sum_f_cyl)
|
# print('Sum of Forces in cylinder: %f' % sum_f_cyl)
|
||||||
|
|
||||||
# sum_f_inner_cyl = np.sum(F[inner_cylinder])
|
# sum_f_inner_cyl = np.sum(N[inner_cylinder])
|
||||||
# print('Sum of Forces in inner cylinder: %f' % sum_f_inner_cyl)
|
# print('Sum of Forces in inner cylinder: %f' % sum_f_inner_cyl)
|
||||||
|
|
||||||
# if sum_f > 4000000.000000:
|
# if sum_f > 4000000.000000:
|
||||||
# test = 1
|
# test = 1
|
||||||
|
|
||||||
# F[Ny - 1, :, 5] += 0.1
|
# N[Ny - 1, :, 5] += 0.1
|
||||||
# F[0, :, 1] -= 0.1
|
# N[0, :, 1] -= 0.1
|
||||||
# F[0, :, 5] += 0.1
|
# N[0, :, 5] += 0.1
|
||||||
# F[Ny - 1, :, 1] -= 0.1
|
# N[Ny - 1, :, 1] -= 0.1
|
||||||
|
|
||||||
# Calculate fluid variables
|
# Calculate fluid variables
|
||||||
rho = np.sum(F, 2)
|
rho = np.sum(N, 2)
|
||||||
ux = np.sum(F * cxs, 2) / rho
|
temp_rho = np.sum(temperature, 2)
|
||||||
uy = np.sum(F * cys, 2) / rho
|
ux = np.sum(N * cxs, 2) / rho
|
||||||
|
uy = np.sum(N * cys, 2) / rho
|
||||||
|
|
||||||
ux[(np.abs(rho) < 0.000000001)] = 0
|
ux[(np.abs(rho) < epsilon)] = 0
|
||||||
uy[(np.abs(rho) < 0.000000001)] = 0
|
uy[(np.abs(rho) < epsilon)] = 0
|
||||||
|
|
||||||
|
g = -params.g * (temp_rho - yy / Ny)
|
||||||
|
# uy[np.abs(rho) >= epsilon] += g[np.abs(rho) >= epsilon] / 2.0
|
||||||
|
uy += g / 2.0
|
||||||
|
|
||||||
|
# u_length = np.maximum(np.abs(ux), np.abs(uy))
|
||||||
|
u_length = np.sqrt(np.square(ux) + np.square(uy))
|
||||||
|
|
||||||
|
u_max_length = np.max(u_length)
|
||||||
|
if u_max_length > np.sqrt(2):
|
||||||
|
ux = (ux / u_max_length) * np.sqrt(2)
|
||||||
|
uy = (uy / u_max_length) * np.sqrt(2)
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
# print('minimum rho: %f' % np.min(np.abs(rho)))
|
# print('minimum rho: %f' % np.min(np.abs(rho)))
|
||||||
# print('Maximum F: %f' % np.max(F))
|
# print('Maximum N: %f' % np.max(N))
|
||||||
# print('Minimum F: %f' % np.min(F))
|
# print('Minimum N: %f' % np.min(N))
|
||||||
|
|
||||||
# Apply Collision
|
# Apply Collision
|
||||||
Feq = np.zeros(F.shape)
|
temperature_eq = np.zeros(temperature.shape)
|
||||||
|
Neq = np.zeros(N.shape)
|
||||||
for i, cx, cy, w in zip(idxs, cxs, cys, weights):
|
for i, cx, cy, w in zip(idxs, cxs, cys, weights):
|
||||||
Feq[:, :, i] = rho * w * (
|
Neq[:, :, i] = rho * w * (
|
||||||
1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2)
|
1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2)
|
||||||
|
|
||||||
F += -(1.0 / tau) * (F - Feq)
|
temperature_eq[:, :, i] = temp_rho * w * (
|
||||||
|
1 + 3 * (cx * ux + cy * uy) + 9 * (cx * ux + cy * uy) ** 2 / 2 - 3 * (ux ** 2 + uy ** 2) / 2)
|
||||||
|
# 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)
|
||||||
|
|
||||||
# Apply boundary
|
# Apply boundary
|
||||||
F[bndry, :] = bndryF
|
N[bndry, :] = bndryN
|
||||||
|
temperature[bndry, :] = bndryTemp
|
||||||
|
|
||||||
|
# temperature[0, :, 0] = 0
|
||||||
|
# temperature[1, :, 0] = 0
|
||||||
|
|
||||||
|
temperature[0, :, 0] /= 2
|
||||||
|
temperature[1, :, 0] /= 2
|
||||||
|
|
||||||
|
temperature[Ny - 1, :, 0] = 1
|
||||||
|
temperature[Ny - 2, :, 0] = 1
|
||||||
|
|
||||||
|
# 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)))
|
||||||
|
|
||||||
# plot in real time - color 1/2 particles blue, other half red
|
# plot in real time - color 1/2 particles blue, other half red
|
||||||
if (plotRealTime and (it % 10) == 0) or (it == Nt - 1):
|
if (plotRealTime and (it % 10) == 0) or (it == Nt - 1):
|
||||||
|
fig.clear()
|
||||||
plt.cla()
|
plt.cla()
|
||||||
ux[cylinder] = 0
|
ux[cylinder] = 0
|
||||||
uy[cylinder] = 0
|
uy[cylinder] = 0
|
||||||
|
|
@ -158,16 +264,29 @@ def main():
|
||||||
cmap = plt.cm.bwr
|
cmap = plt.cm.bwr
|
||||||
cmap.set_bad('black')
|
cmap.set_bad('black')
|
||||||
|
|
||||||
# plt.imshow(vorticity, cmap='bwr')
|
plt.subplot(2, 2, 1)
|
||||||
plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr')
|
plt.imshow(vorticity, cmap='bwr')
|
||||||
|
plt.clim(-.1, .1)
|
||||||
|
# plt.imshow(has_fluid * 2.0 - 1.0, cmap='bwr')
|
||||||
# plt.imshow(bndry * 2.0 - 1.0, cmap='bwr')
|
# plt.imshow(bndry * 2.0 - 1.0, cmap='bwr')
|
||||||
|
# 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')
|
||||||
|
|
||||||
plt.clim(-.1, .1)
|
plt.clim(-.1, .1)
|
||||||
ax = plt.gca()
|
# ax = plt.gca()
|
||||||
ax.invert_yaxis()
|
# ax.invert_yaxis()
|
||||||
ax.get_xaxis().set_visible(False)
|
# ax.get_xaxis().set_visible(False)
|
||||||
ax.get_yaxis().set_visible(False)
|
# ax.get_yaxis().set_visible(False)
|
||||||
ax.set_aspect('equal')
|
# ax.set_aspect('equal')
|
||||||
plt.pause(0.001)
|
plt.pause(0.001)
|
||||||
|
|
||||||
# Save figure
|
# Save figure
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue