This is an interactive Python notebook. The grey cells contain runnable Python code which can be executed with Ctrl+Enter. Make sure to execute all of them in the right order, since they built upon each other. You can also execute them all in the beginning using "Cell->Run all". To get back to a clean state use "Kernel->Restart".
This first tutorial introduces waLBerla's basic data structures, by implementing a version of Conway's Game of Life cellular automaton. The Game of Life algorithm is formulated on a regular grid of cells. Each cell can be in one of two states: dead or alive. The time evolution of this game is a simple rule how to get to the next cell state by only using cell states of neighboring cells. For details see the Wikipedia page.
Before working with waLBerla we first start with a pure Python implementation of this algorithm using the popular numpy package. If you are not familiar with that package, now is a good time to read up on it.
In [ ]:
import numpy as np
def makeGrid(shape):
return np.zeros(shape, dtype=np.int8)
print(makeGrid( [5,5] ))
We have created an empty playing field consisting of 8-bit integer values, all initialized with zeros. Now we write a function iterating over all cells, applying the update rule on each one of them. In waLBerla these update functions are called sweeps, since they iterate (sweep) over the complete domain and update each cell separately. A crucial point here is that at each cell, we only have to access neighboring values, and, since we create a new temporary copy, all of these cell updates could in principle happen in parallel.
In [ ]:
ALIVE = 1
DEAD = 0
neighborhoodD2Q9 = [ (i,j) for i in [-1,0,1]
for j in [-1,0,1]
if i != 0 or j!=0 ]
def gameOfLifeSweep(grid):
temporaryGrid = np.copy(grid)
for i in range(1, grid.shape[0]-1):
for j in range(1, grid.shape[1]-1):
numberOfAliveNeighbors = 0
for neighborCellOffset in neighborhoodD2Q9:
ni = i + neighborCellOffset[0]
nj = j + neighborCellOffset[1]
if temporaryGrid[ni,nj] == ALIVE:
numberOfAliveNeighbors += 1
if numberOfAliveNeighbors < 2 or numberOfAliveNeighbors > 3:
grid[i,j] = DEAD
if numberOfAliveNeighbors == 3:
grid[i,j] = ALIVE
return grid
This code snippet takes a grid and returns the grid in the next timestep. For the GameOfLife we use a D2Q9 neighborhood, meaning 2 dimensional ( D2 ) with in total 9 cells (Q9). Since we leave out the center (0,0) we strictly speaking have only 8 cells.
In [ ]:
print(neighborhoodD2Q9)
Lets first initialize a so-called blinker and run a few timesteps
In [ ]:
grid = makeGrid( [5,5] )
grid[2,1:4] = ALIVE
print ( "Initial Setup:" )
print ( grid )
for t in range(2):
grid = gameOfLifeSweep(grid)
print("After timestep %d: " % (t+1,) )
print(grid)
Our implementation seems to work: The blinker is iterating ( blinking ) nicely between these two configurations.
Looking at these configurations in text representation is not a good idea for bigger grids, so lets display our grid using matplotlib. The setup code for matplotlib was put to a separate file called matplotlib_setup.py which also contains code for creating animations.
In [ ]:
from material.matplotlib_setup import * # import matplotlib and configures it to play nicely with iPython notebook
matplotlib.rcParams['image.cmap'] = 'Blues' # switch default colormap
im = plt.imshow(grid, interpolation='none')
Now lets create an animation of this blinker:
In [ ]:
ani = makeImshowAnimation(grid, gameOfLifeSweep, frames=6)
displayAsHtmlVideo(ani, fps=2)
Now lets load some more interesting starting configuration. Here we choose the 'Gosper Gliding Gun' scenario, taken from Wikipedia.
In [ ]:
from scipy.misc import imread
grid = imread('material/GosperGliderGun.png',flatten=True).astype(int)
grid[grid>0] = ALIVE # values are from 0 to 255 - set everything nonzero to ALIVE
ani = makeImshowAnimation(grid, gameOfLifeSweep, frames=6*15)
displayAsHtmlVideo(ani, fps=15)
waLBerla is parallelized using MPI (message passing interface). That means that multiple processes are started, possibly on different machines which all execute the same program and communicate by sending messages. In a normal Python environment we would start multiple Python interpreters all executing the same script using mpirun. The following command would run a script using 4 processes:
mpirun -np 4 python3 my_waLBerla_script.py
To make use of multiple processes in this IPython environment we use the ipyparallel package. Therefore you first have to start an ipyparallel cluster here in the IPython Cluster tab. Otherwise the following code will not work.
Now lets test our parallel environment. After importing and initializing ipyparallel we only have to prefix a code cell with the IPython magic px and the code cell will be executed by each process of the cluster.
In [ ]:
from ipyparallel import Client
rc = Client()
numberOfProcesses = len(rc.ids)
print("There are %d processes in your ipyparallel cluster" % ( numberOfProcesses,))
In [ ]:
%%px
print("Hello parallel world")
If above cells produce errors, go back to the IPCluster setup to make sure your cluster is running. If everything works you see the hello world message printed by each process in the cluster.
Now lets implement the Game of Life algorithm in parallel using waLBerla. waLBerla divides the complete domain into blocks. These blocks can then be distributed to the participating processes. We will set up here a simple example where each process gets one block. waLBerla can put multiple blocks per process. This makes sense if the computational load varies for each block, then a process receives either few expensive blocks or many cheap blocks.
While blocks are the basic unit of load balancing they also act as container for distributed data. In this Game of Life example we have to distribute our grid to the blocks e.g. each block internally stores only part of the complete domain. The local grids are extended by one ghost layer which are a shadow copy of the outermost layer of the neighboring block. All these details are handled by waLBerla.
We only have to specify the number of blocks (which equals the number of processes in our case) and how many cells we want to have on each block. Each of these size informations have to be 3 tuples since waLBerla inherently is a 3D framework. However we can mimic our 2D setup by choosing the size of the z coordinate as 1.
In [ ]:
from scipy.misc import imread
import os
# Read the initial scenario
initialConfig = np.rot90( imread('material/GosperGliderGun.png',flatten=True).astype(int), 3 )
initialConfig[initialConfig>0] = ALIVE # values are from 0 to 255 - set everything nonzero to ALIVE
# and send it to all execution engines together with other local variables
rc[:].push( {'initialConfig' : initialConfig,
'numberOfProcesses' : numberOfProcesses,
'ALIVE' : ALIVE,
'DEAD' : DEAD,
'neighborhoodD2Q9' : neighborhoodD2Q9,
'gameOfLifeSweep' : gameOfLifeSweep,
'cwd' : os.getcwd() } )
pass
In [ ]:
%%px
import sys
import waLBerla as wlb
import numpy as np
import os
os.chdir( cwd ) # change working directory of all engines to working directory of notebook kernel
domainSize = (initialConfig.shape[0], initialConfig.shape[1], 1)
# We can either specify the detailed domain partitioning ...
blocks = wlb.createUniformBlockGrid(cellsPerBlock=(domainSize[0]//numberOfProcesses, domainSize[1], domainSize[2]),
blocks=(numberOfProcesses,1,1),
periodic=(1,1,1))
# ... or let waLBerla do it automatically -
# if the domainSize is not divisible by the block count the domain is slightly extended
blocks = wlb.createUniformBlockGrid(cells = domainSize, periodic=(1,1,1))
# Now put one field (i.e. grid) on each block
wlb.field.addToStorage(blocks, name='PlayingField', type=np.int, ghostLayers=1)
# Iterate over local blocks - in our setup we have exactly one block per process - but lets be general
for block in blocks:
offsetInGlobalDomain = blocks.transformLocalToGlobal(block, (0,0,0))
myRank = wlb.mpi.rank()
print("Block on rank %d: with offset %s" % (myRank, offsetInGlobalDomain ))
Having a distributed domain, we can now run the stencil algorithm locally on each block. We can safely do this, since the simulation domain was extended with one ghost layer. After each sweep over the domain the ghost layers have to be synchronized again.
In [ ]:
%%px
import waLBerla.plot as wlbPlt
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from material.matplotlib_setup import *
matplotlib.rcParams['image.cmap'] = 'Blues' # switch default colormap
# Initialize our grid
for block in blocks:
grid = wlb.field.toArray( block['PlayingField'] ).squeeze()
offsetInGlobalDomain = blocks.transformLocalToGlobal(block, (0,0,0))
blockSize = grid.shape[:2]
xBegin, xEnd = offsetInGlobalDomain[0], offsetInGlobalDomain[0] + grid.shape[0]
yBegin, yEnd = offsetInGlobalDomain[1], offsetInGlobalDomain[1] + grid.shape[1]
grid[:,:] = initialConfig[xBegin:xEnd, yBegin:yEnd]
wlbPlt.scalarField( blocks, 'PlayingField', wlb.makeSlice[:,:,0] )
In [ ]:
%%px
communication = wlb.createUniformBufferedScheme( blocks, 'D2Q9')
communication.addDataToCommunicate( wlb.field.createPackInfo( blocks, 'PlayingField') )
def runTimestep():
communication()
for block in blocks:
grid = wlb.field.toArray( block['PlayingField'], withGhostLayers=True )[:,:,1,0]
gameOfLifeSweep( grid )
ani = wlbPlt.scalarFieldAnimation( blocks, 'PlayingField', wlb.makeSlice[:,:,0], runFunction=runTimestep, frames=100 )
displayAsHtmlVideo( ani, fps=30, show=(wlb.mpi.rank()==0) )