#!/usr/bin/env python3 # LatticeBoltzmannDemo.py: a two-dimensional lattice-Boltzmann "wind tunnel" simulation # Uses numpy to speed up all array handling. # Uses matplotlib to plot and animate the curl of the macroscopic velocity field. # Copyright 2013, Daniel V. Schroeder (Weber State University) 2013 # Permission is hereby granted, free of charge, to any person obtaining a copy of # this software and associated data and documentation (the "Software"), to deal in # the Software without restriction, including without limitation the rights to # use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies # of the Software, and to permit persons to whom the Software is furnished to do # so, subject to the following conditions: # The above copyright notice and this permission notice shall be included in all # copies or substantial portions of the Software. # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, # INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A # PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR # ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR # OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. # Except as contained in this notice, the name of the author shall not be used in # advertising or otherwise to promote the sale, use or other dealings in this # Software without prior written authorization. # Credits: # The "wind tunnel" entry/exit conditions are inspired by Graham Pullan's code # (http://www.many-core.group.cam.ac.uk/projects/LBdemo.shtml). Additional inspiration from # Thomas Pohl's applet (http://thomas-pohl.info/work/lba.html). Other portions of code are based # on Wagner (http://www.ndsu.edu/physics/people/faculty/wagner/lattice_boltzmann_codes/) and # Gonsalves (http://www.physics.buffalo.edu/phy411-506-2004/index.html; code adapted from Succi, # http://global.oup.com/academic/product/the-lattice-boltzmann-equation-9780199679249). # For related materials see http://physics.weber.edu/schroeder/fluids import sys import numpy, time, matplotlib.pyplot, matplotlib.animation from random import Random from math import sqrt from PIL import Image # Define constants: height = 80 # lattice dimensions width = 200 # shown width will be less by 4 viscosity = 0.01 # fluid viscosity viscosity = 0.015 omega = 1 / (3*viscosity + 0.5) # "relaxation" parameter u0x = 0.1 # initial and in-flow speed, x component. u0y = 0.0 # initial and in-flow speed, y component. four9ths = 4.0/9.0 # abbreviations for lattice-Boltzmann weight factors one9th = 1.0/9.0 one36th = 1.0/36.0 performanceData = True # set to True if performance data is desired brownian = 1.0e-2 #brownian = 0.0 steps = 20 # steps / rendered frame capture_frame = 0 # which frame to capture if > 0 rng_seed = 12345678 args = [float(x) for x in sys.argv[1:]] if len(args) >= 2: u0x = args[0] u0y = args[1] if len(args) == 3: capture_frame = args[2] rng = Random() rng.seed(rng_seed) # Initialize all the arrays to steady rightward flow: # particle densities along 9 directions n0 = four9ths * (numpy.ones((height,width)) - 1.5*u0x**2 - 1.5*u0y**2) nN = one9th * numpy.ones((height,width)) nS = one9th * numpy.ones((height,width)) nE = one9th * numpy.ones((height,width)) nW = one9th * numpy.ones((height,width)) nNE = one36th * numpy.ones((height,width)) nSE = one36th * numpy.ones((height,width)) nNW = one36th * numpy.ones((height,width)) nSW = one36th * numpy.ones((height,width)) nE += one9th * (3*u0x + 3*u0x**2 - 1.5*u0y**2) nW += one9th * (-3*u0x + 3*u0x**2 - 1.5*u0y**2) nN += one9th * (3*u0y + 3*u0y**2 - 1.5*u0x**2) nS += one9th * (-3*u0y + 3*u0y**2 - 1.5*u0x**2) nNE += one36th * (3*u0x + 3*u0x**2 + 3*u0y + 3*u0y**2) nSE += one36th * (3*u0x + 3*u0x**2 - 3*u0y + 3*u0y**2) nNW += one36th * (-3*u0x + 3*u0x**2 + 3*u0y + 3*u0y**2) nSW += one36th * (-3*u0x + 3*u0x**2 - 3*u0y + 3*u0y**2) # Add pseudo-Brownian perturbation for y in range(height): for x in range(width): bx = 0.1 * rng.randint(-9, 9) by = sqrt(1.0 - bx**2) * 2 * (rng.randint(0, 1) - 0.5) br = 0.1 * rng.randint(0, 10) * brownian bx *= br by *= br nE[y, x] += one9th * (3*bx + 3*bx**2 - 1.5*by**2) nW[y, x] += one9th * (-3*bx + 3*bx**2 - 1.5*by**2) nN[y, x] += one9th * (3*by + 3*by**2 - 1.5*bx**2) nS[y, x] += one9th * (-3*by + 3*by**2 - 1.5*bx**2) nNE[y, x] += one36th * (3*bx + 3*bx**2 + 3*by + 3*by**2) nSE[y, x] += one36th * (3*bx + 3*bx**2 - 3*by + 3*by**2) nNW[y, x] += one36th * (-3*bx + 3*bx**2 + 3*by + 3*by**2) nSW[y, x] += one36th * (-3*bx + 3*bx**2 - 3*by + 3*by**2) rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW # macroscopic density ux = (nE + nNE + nSE - nW - nNW - nSW) / rho # macroscopic x velocity uy = (nN + nNE + nNW - nS - nSE - nSW) / rho # macroscopic y velocity # Initialize barriers: barrier = numpy.zeros((height,width), bool) # True wherever there's a barrier barrier[(height//2)-8:(height//2)+8, height//2] = True # simple linear barrier diagBarrier = numpy.zeros((height,width), bool) for i in range(-8, 9): diagBarrier[(height // 2 + i), (height // 2 + i)] = True diagBarrier[(height // 2 + i), (height // 2 + i + 1)] = True # barrier = diagBarrier lBarrier = numpy.zeros((height, width), bool) for i in range(-8, 9): lBarrier[height // 2 + i, height // 2 - 8] = True lBarrier[height // 2 - 8, height // 2 + i] = True # barrier = lBarrier im = Image.open('wing2d.png') pix = im.load() (imWidth, imHeight) = im.size imBarrier = numpy.zeros((height, width), bool) for x in range(64): for y in range(64): (R, G, B, A) = pix[imWidth * x // 64, imHeight * y // 64] if (A > 127): imBarrier[79 - y, 24 + x] = True im.close() barrier = imBarrier barrierN = numpy.roll(barrier, 1, axis=0) # sites just north of barriers barrierS = numpy.roll(barrier, -1, axis=0) # sites just south of barriers barrierE = numpy.roll(barrier, 1, axis=1) # etc. barrierW = numpy.roll(barrier, -1, axis=1) barrierNE = numpy.roll(barrierN, 1, axis=1) barrierNW = numpy.roll(barrierN, -1, axis=1) barrierSE = numpy.roll(barrierS, 1, axis=1) barrierSW = numpy.roll(barrierS, -1, axis=1) # Move all particles by one step along their directions of motion (pbc): def stream(): global nN, nS, nE, nW, nNE, nNW, nSE, nSW nN = numpy.roll(nN, 1, axis=0) # axis 0 is north-south; + direction is north nNE = numpy.roll(nNE, 1, axis=0) nNW = numpy.roll(nNW, 1, axis=0) nS = numpy.roll(nS, -1, axis=0) nSE = numpy.roll(nSE, -1, axis=0) nSW = numpy.roll(nSW, -1, axis=0) nE = numpy.roll(nE, 1, axis=1) # axis 1 is east-west; + direction is east nNE = numpy.roll(nNE, 1, axis=1) nSE = numpy.roll(nSE, 1, axis=1) nW = numpy.roll(nW, -1, axis=1) nNW = numpy.roll(nNW, -1, axis=1) nSW = numpy.roll(nSW, -1, axis=1) """ # Use tricky boolean arrays to handle barrier collisions (bounce-back): nN[barrierN] = nS[barrier] nS[barrierS] = nN[barrier] nE[barrierE] = nW[barrier] nW[barrierW] = nE[barrier] nNE[barrierNE] = nSW[barrier] nNW[barrierNW] = nSE[barrier] nSE[barrierSE] = nNW[barrier] nSW[barrierSW] = nNE[barrier] """ # Use boolean indexing to handle barrier collisions (bounce-back): nN[barrierN] += nS[barrier] nS[barrierS] += nN[barrier] nE[barrierE] += nW[barrier] nW[barrierW] += nE[barrier] nNE[barrierNE] += nSW[barrier] nNW[barrierNW] += nSE[barrier] nSE[barrierSE] += nNW[barrier] nSW[barrierSW] += nNE[barrier] nS[barrier] = 0.0 nN[barrier] = 0.0 nW[barrier] = 0.0 nE[barrier] = 0.0 nSW[barrier] = 0.0 nSE[barrier] = 0.0 nNW[barrier] = 0.0 nNE[barrier] = 0.0 # Collide particles within each cell to redistribute velocities (could be optimized a little more): def collide(): global rho, ux, uy, n0, nN, nS, nE, nW, nNE, nNW, nSE, nSW rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW rho += 1e-9 ux = (nE + nNE + nSE - nW - nNW - nSW) / rho uy = (nN + nNE + nNW - nS - nSE - nSW) / rho ux2 = ux * ux # pre-compute terms used repeatedly... uy2 = uy * uy u2 = ux2 + uy2 omu215 = 1 - 1.5*u2 # "one minus u2 times 1.5" uxuy = ux * uy n0 = (1-omega)*n0 + omega * four9ths * rho * omu215 nN = (1-omega)*nN + omega * one9th * rho * (omu215 + 3*uy + 4.5*uy2) nS = (1-omega)*nS + omega * one9th * rho * (omu215 - 3*uy + 4.5*uy2) nE = (1-omega)*nE + omega * one9th * rho * (omu215 + 3*ux + 4.5*ux2) nW = (1-omega)*nW + omega * one9th * rho * (omu215 - 3*ux + 4.5*ux2) nNE = (1-omega)*nNE + omega * one36th * rho * (omu215 + 3*(ux+uy) + 4.5*(u2+2*uxuy)) nNW = (1-omega)*nNW + omega * one36th * rho * (omu215 + 3*(-ux+uy) + 4.5*(u2-2*uxuy)) nSE = (1-omega)*nSE + omega * one36th * rho * (omu215 + 3*(ux-uy) + 4.5*(u2-2*uxuy)) nSW = (1-omega)*nSW + omega * one36th * rho * (omu215 + 3*(-ux-uy) + 4.5*(u2+2*uxuy)) # Force steady flow at ends: n0[:, 0] = four9ths * (1 - 1.5*u0x**2 - 1.5*u0y**2) nN[:, 0] = one9th nS[:, 0] = one9th nE[:, 0] = one9th nW[:, 0] = one9th nNE[:, 0] = one36th nSE[:, 0] = one36th nNW[:, 0] = one36th nSW[:, 0] = one36th nE[:, 0] += one9th * (3*u0x + 3*u0x**2 - 1.5*u0y**2) nW[:, 0] += one9th * (-3*u0x + 3*u0x**2 - 1.5*u0y**2) nN[:, 0] += one9th * (3*u0y + 3*u0y**2 - 1.5*u0x**2) nS[:, 0] += one9th * (-3*u0y + 3*u0y**2 - 1.5*u0x**2) nNE[:, 0] += one36th * (3*u0x + 3*u0x**2 + 3*u0y + 3*u0y**2) nSE[:, 0] += one36th * (3*u0x + 3*u0x**2 - 3*u0y + 3*u0y**2) nNW[:, 0] += one36th * (-3*u0x + 3*u0x**2 + 3*u0y + 3*u0y**2) nSW[:, 0] += one36th * (-3*u0x + 3*u0x**2 - 3*u0y + 3*u0y**2) # Compute curl of the macroscopic velocity field: def curl(ux, uy): temp = (numpy.roll(uy, -1, axis=1) - numpy.roll(uy, 1, axis=1) - numpy.roll(ux, -1, axis=0) + numpy.roll(ux, 1, axis=0)) temp[:, 0:4] = temp[:, 4:8] return temp[:, :(width - 4)] # Here comes the graphics and animation... theFig = matplotlib.pyplot.figure(figsize=(8,3)) perfText = matplotlib.pyplot.text(width - 50, height - 10, '', fontsize=8) frameText = matplotlib.pyplot.text(width - 50, height - 20, '', fontsize=8) fluidImage = matplotlib.pyplot.imshow(curl(ux, uy), origin='lower', norm=matplotlib.pyplot.Normalize(-.1,.1), cmap=matplotlib.pyplot.get_cmap('jet'), interpolation='none') # See http://www.loria.fr/~rougier/teaching/matplotlib/#colormaps for other cmap options bImageArray = numpy.zeros((height, width - 4, 4), numpy.uint8) # an RGBA image bImageArray[barrier[:, :(width - 4)],3] = 255 # set alpha=255 only at barrier sites barrierImage = matplotlib.pyplot.imshow(bImageArray, origin='lower', interpolation='none') startTime = time.process_time() #frameList = open('frameList.txt','w') # file containing list of images (to make movie) # Function called for each successive animation frame: def nextFrame(arg): # (arg is the frame number, which we don't need) global startTime, perfText, fluidImage, steps, frameText, rng_seed, theFig if performanceData and (arg%25 == 0) and (arg>0): endTime = time.process_time() perfString = "{:1.1f} FPS, {:d} steps /F Velocity: ({:1.2f}, {:1.2f})".format(25 / (endTime - startTime), steps, u0x, u0y) for txt in theFig.texts: txt.set_visible(False) perfText.remove() perfText = matplotlib.pyplot.text(width - 100, height - 10, perfString, fontsize=8) startTime = endTime if (capture_frame > 0) and (arg == capture_frame): framename = "boltz2D_frame{:04d}S{:d}.png".format(arg, rng_seed) frameString = "Frame {:d}".format(arg) frameText = matplotlib.pyplot.text(width - 50, height - 20, frameString, fontsize=8) matplotlib.pyplot.savefig(framename) for step in range(steps): # adjust number of steps for smooth animation stream() collide() fluidImage.set_array(curl(ux, uy)) return (fluidImage, barrierImage, perfText) # return the figure elements to redraw animate = matplotlib.animation.FuncAnimation(theFig, nextFrame, interval=1, blit=True, save_count=10) matplotlib.pyplot.show()