thermal lbm

lbm.py

# paieska google videos 'thermal lattice bolzman'
# Robert Straka: https://www.youtube.com/watch?v=cvYfGpiJ5ec
# geras review, sharma 2020: https://sci-hub.se/10.1016/j.paerosci.2020.100616
# kuris cituoja:
# (yra apie D2Q5) Sharma 2018, Keerti Vardhan, Robert Straka, and Frederico Wanderley Tavares. "Natural convection heat transfer modeling by the cascaded thermal lattice Boltzmann method." International Journal of Thermal Sciences 134 (2018): 552-564.

#jei be standalone (async) serverio, paleisti reikia su:
#panel serve lbm.py --autoreload --dev

import numpy as np
from plotly import express as px
import holoviews as hv
import panel as pn

# Simulation parameters
Nx          = 100    # resolution x-dir
Ny          = 100    # resolution y-dir
# rho0        = 100    # average density
tau_f         = 0.6    # collision timescale
tau_g         = 0.99    # collision timescale
Tref = 1
beta = 0.9# linear isobaric thermal expansion coefficient 
g = np.array([0, -0.01]) # gravity

# Lattice speeds / weights
NL5 = 5
NL9 = 9
idxs5 = np.arange(NL5)
idxs9 = np.arange(NL9)
cxs5 = np.array([0, -1, 0, 1, 0])
cys5 = np.array([0, 0, -1, 0, 1])
cxs9 = np.array([0,-1,-1, -1, 0,  1, 1, 1, 0])
cys9 = np.array([0, 1, 0, -1,-1, -1, 0, 1, 1])
weights5 = np.array([1/3, 1/6, 1/6, 1/6, 1/6]) # sums to 1
weights9 = np.array([4/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9, 1/36, 1/9]) # sums to 1

X, Y = np.meshgrid(range(Nx), range(Ny))

# Initial Conditions - flow to the right with some perturbations
global G
global F
G = np.ones((Ny,Nx,NL5), dtype=np.float64) #+ 0.01*np.random.randn(Ny,Nx,NL)
F = np.ones((Ny,Nx,NL9), dtype=np.float64) #+ 0.01*np.random.randn(Ny,Nx,NL)

rho = np.sum(F,2)
T = np.sum(G,2)

#init cond
G[Ny//2,Nx//2,:] *= 1.1
# F[Ny//2,Nx//2,:] *= 1.1

# Cylinder boundary
cylinder = (X - Nx/4)**2 + (Y - Ny/2)**2 < (Ny/4)**2


def calc_iter(i):
  print(f'iter={i}')
  global G
  global F
  # Drift
  for i, cx, cy in zip(idxs5, cxs5, cys5):
    G[:,:,i] = np.roll(G[:,:,i], cx, axis=1)
    G[:,:,i] = np.roll(G[:,:,i], cy, axis=0)

  for i, cx, cy in zip(idxs9, cxs9, cys9):
    F[:,:,i] = np.roll(F[:,:,i], cx, axis=1)
    F[:,:,i] = np.roll(F[:,:,i], cy, axis=0)    

 
  # Calculate fluid variables
  T = np.sum(G,2)
  ux_g  = np.sum(G*cxs5,2) / T
  uy_g  = np.sum(G*cys5,2) / T
  rho = np.sum(F,2)
  ux_f  = np.sum(F*cxs9,2) / rho
  uy_f  = np.sum(F*cys9,2) / rho  
  
  # Apply Collision
  Geq = np.zeros(G.shape)
  Gi = np.zeros(G.shape)
  for i, cx, cy, w in zip(idxs5, cxs5, cys5, weights5):
    Geq[:,:,i] = T*w* (1 + (cx*ux_f+cy*uy_f))
    Gi[:, :, i] = cxs5[i]*g[0]*beta*(T-Tref) + cys5[i]*g[1]*beta*(T-Tref) # Boussinesq    
  
  Feq = np.zeros(F.shape)
  Fi = np.zeros(F.shape)
  for i, cx, cy, w in zip(idxs9, cxs9, cys9, weights9):
    Feq[:,:,i] = rho*w* (1 + 3*(cx*ux_f+cy*uy_f) + 9*(cx*ux_f+cy*uy_f)**2/2 - 3*(ux_f**2+uy_f**2)/2)
    Fi[:, :, i] = cxs9[i]*g[0]*beta*(T-Tref) + cys9[i]*g[1]*beta*(T-Tref) # Boussinesq

  G += -(1.0/tau_g) * (G - Geq)# + Gi
  F += -(1.0/tau_f) * (F - Feq) + Fi
  
pn.extension(sizing_mode="stretch_width")
hv.extension('plotly')
# hv.extension('bokeh')

button = pn.widgets.Button(name='Button')
i = 0
def update(e):
  global i
  # global F
  # print(Z.shape)
  calc_iter(i)
  # F[5,5,:] = i
  layout[0].object.data[0].z = F.sum(axis=2)
  layout[1].object.data[0].z = G.sum(axis=2)
  layout[0].object.layout.title = f'F, i={i}'
  layout[1].object.layout.title = 'T'
  i += 1

button.on_click(update)
layout = pn.Column(px.imshow(F.sum(axis=2)),px.imshow(G.sum(axis=2)), button)

layout.servable() # run command: panel serve lbm.py --autoreload --dev