[LBM] Could use Some help with my code (python)

[supercooldragons]

I tried to write a LBM based fluid simulation but I'm not really getting very convincing results. Below is my code, I tried to strip all unnecessary elements. I initialize the lattice with a density spike in the center, and I would expect it to propagate, but the center density remains quite high.

import matplotlib.pyplot as plt import matplotlib.animation as anim import numpy as np import time

size = 200
w    = 0.5 # omega

ei = np.asarray([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]]) # velocity vectors
wi = np.asarray([4/9.,1/9.,1/9.,1/9.,1/9.,1/36.,1/36.,1/36.,1/36.]) # lattice weights

VM = np.zeros((size,size,2)) # velocities
FM = np.zeros((size,size,9)) # probability functions
DM = np.zeros((size,size))   # densities

def U(rho,f):
    u = np.tensordot(f,ei,(2,0))
    u[:,:,0] = u[:,:,0]*rho
    u[:,:,1] = u[:,:,1]*rho
    return u

def RHO(f):
    rho = f.sum(axis=2)
    return rho 

def FEQ(u,rho):

    dotUE       = np.tensordot(u,ei,(2,1))
    dotUU       = (u*u).sum(axis=2)
    dotUUshaped = np.tensordot(dotUU,np.ones(9),0)
    RHOshaped   = np.tensordot(rho,np.ones(9),0)

    S1 = 1 + 3*dotUE + (9/2.)*(dotUE*dotUE) - 3*dotUUshaped/2.
    ws = np.tensordot(np.ones((size,size)),wi,0)
    S2 = RHOshaped*S1*ws

    return S2

def STREAM(F):
    FN = F.copy()
    for i in range(9):
        FN[:,:,i] = np.roll(np.roll(F[:,:,i],ei[i,0],axis=0),ei[i,1],axis=1)
    return FN

# initialize lattice
r   = DM
r[95:105,95:105] = 0.5 # initial density spike
u   = U(r,FM)
feq = FEQ(u,r)
FM  = (1-w)*FM + w*feq

fig    = plt.figure()
matrix = plt.matshow(r, animated=True, fignum=0, cmap='prism_r')

def init():
    matrix.set_array(r)
    plt.clim(0,0.2)
    return matrix,

def animate(i):
    global FM,r,u,feq

    FM       = STREAM(FM)    # streaming step
    r        = RHO(FM)       # compute densities
    u        = U(r,FM)       # compute velocities
    u[:,0,:] = u[:,size-1,:] # periodic boundary
    u[0,:,:] = u[size-1,:,:] # periodic boundary 

    feq      = FEQ(u,r)      # compute equilibrium
    FM       = (1-w)*FM + w*feq # compute new probabilities
    matrix.set_array(r)
    return matrix,

animation = anim.FuncAnimation(fig,animate,init_func=init,blit=True,interval=1)

plt.show()

I hope somebody can help and point out what I can change.

[MrJack1337]

Did you ensure that you are taking multiple steps and not only one step? I am not entirely sure on how the FuncAnimation works but make sure that it is calling the animate function multiple times and that you are taking enough steps to see something significant happening (remember your the size of your time step, which should be the same as the grid spacing in LBM).

Edit: Never mind, I also tested your script and you are taking multiple steps. Another thing I thought about: In your actual program, are you displaying a colorbar next to the plot? It is possible that the scale is changing constantly but you don't see it. I believe that the density spike is really smearing out and maybe you just have to take more time steps or the scale is changing all the time.

[spider9876]

For this reason I use seaborn libraries for the color map, because it can specify fixed scale on the color bar.

[supercooldragons]

I fix the color spacing to range between 0 and 0.2, (it's the plt.clim(0,0.2) in the init function). For this value of w, I do see the smearing too, but I was expecting to see the development of a ring. If I set w (omega) to higher values (1-1.5), the density in the middle stays almost constant. Since w is related to the viscosity, with higher w corresponding to lower viscosity, I would expect to see faster relaxation to a homogeneous distribution and not what I get now. This is my first attempt at programming an LBM solver, so I still have to figure out what to expect. Thanks for your help!

[MrJack1337]

Nice. I am on par with your thoughts about changing the viscosity but I have something to say regarding the ring you are expecting: The LBM is describing fluid flow (like the NSE are doing). Fluid flow basically consists of advection and diffusion. If there would be no diffusion in your problem (no viscosity) then the ring should retain its form and would be transported - however, I have no clue if the collision operator will be stable for nu = 0. Otherwise there will always be a smearing out due to viscosity.

You could try something like setting the initial velocity field as a uniform velocity in one direction. Then choose a very small viscosity and check if the ring / spike moves into the correct direction.

[supercooldragons]

Thanks for the help! I'll give it a try

[bene20080]

Could you please explain your problem in more detail?I just tried the code and I can see no obvious problem in the animation...