Basic 2D LDC using boltzgen for kernel generation

Using cell lists as parameters for multiple non-branching kernels seems
to reduce performance by ~50 MLUPS (for single precision D2Q9).
This might be alleviated by padding the cell lists to enable thread
layout control or by improved kernel dispatching.

On the upside this OpenCL program runs not only on GPUs but is also vectorized on Intel
CPUs yielding about 180 MLUPS (single precision) and - anticlimactically - 85 MLUPS for
double precision on a i7-4790K.
However both these values compare well to the performance of established CPU LBM codes.

1 files changed, 94 insertions, 0 deletions

diff --git a/ldc_2d.py b/ldc_2d.py
new file mode 100644
index 0000000..14e8a5c
--- /dev/null
+++ b/ldc_2d.py
@@ -0,0 +1,94 @@
+import numpy
+import time
+from string import Template
+
+import matplotlib
+matplotlib.use('AGG')
+import matplotlib.pyplot as plt
+
+from boltzgen import LBM, Generator, Geometry
+from boltzgen.lbm.model import D2Q9
+
+from simulation import Lattice, CellList
+
+def MLUPS(cells, steps, time):
+    return cells * steps / time * 1e-6
+
+def generate_moment_plots(lattice, moments):
+    for i, m in enumerate(moments):
+        print("Generating plot %d of %d." % (i+1, len(moments)))
+
+        velocity = numpy.ndarray(shape=tuple(reversed(lattice.geometry.inner_size())))
+        for x, y in lattice.geometry.inner_cells():
+            velocity[y-1,x-1] = numpy.sqrt(m[1,lattice.memory.gid(x,y)]**2 + m[2,lattice.memory.gid(x,y)]**2)
+
+        plt.figure(figsize=(10, 10))
+        plt.imshow(velocity, origin='lower', cmap=plt.get_cmap('seismic'))
+        plt.savefig("result/ldc_2d_%02d.png" % i, bbox_inches='tight', pad_inches=0)
+
+nUpdates = 100000
+nStat    = 5000
+
+geometry = Geometry(512, 512)
+
+print("Generating kernel using boltzgen...\n")
+
+lbm = LBM(D2Q9)
+generator = Generator(
+    descriptor = D2Q9,
+    moments    = lbm.moments(),
+    collision  = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.6))
+
+functions = ['collide_and_stream', 'equilibrilize', 'collect_moments', 'momenta_boundary']
+extras    = ['cell_list_dispatch']
+
+kernel_src = generator.kernel('cl', 'single', 'SOA', geometry, functions, extras) + Template("""
+__kernel void equilibrilize(__global $float_type* f_next,
+                            __global $float_type* f_prev)
+{
+    const unsigned int gid = get_global_id(1)*$size_x + get_global_id(0);
+    equilibrilize_gid(f_next, f_prev, gid);
+}
+
+__kernel void collect_moments(__global $float_type* f,
+                              __global $float_type* moments)
+{
+    const unsigned int gid = get_global_id(1)*$size_x + get_global_id(0);
+    collect_moments_gid(f, moments, gid);
+}
+""").substitute(float_type = 'float', size_x = geometry.size_x)
+
+print("Initializing simulation...\n")
+
+lattice = Lattice(geometry, kernel_src)
+gid = lattice.memory.gid
+
+bulk_cells = CellList(lattice.context, lattice.queue, lattice.float_type,
+    [ gid(x,y) for x, y in geometry.inner_cells() if x > 1 and x < geometry.size_x-2 and y > 1 and y < geometry.size_y-2 ])
+wall_cells = CellList(lattice.context, lattice.queue, lattice.float_type,
+    [ gid(x,y) for x, y in geometry.inner_cells() if x == 1 or y == 1 or x == geometry.size_x-2 ])
+lid_cells = CellList(lattice.context, lattice.queue, lattice.float_type,
+    [ gid(x,y) for x, y in geometry.inner_cells() if y == geometry.size_y-2 ])
+
+lattice.schedule('collide_and_stream_cells', bulk_cells)
+lattice.schedule('velocity_momenta_boundary_cells', wall_cells, numpy.array([0.0, 0.0], dtype=lattice.float_type[0]))
+lattice.schedule('velocity_momenta_boundary_cells', lid_cells,  numpy.array([0.1, 0.0], dtype=lattice.float_type[0]))
+
+print("Starting simulation using %d cells...\n" % lattice.geometry.volume)
+
+moments = []
+
+lastStat = time.time()
+
+for i in range(1,nUpdates+1):
+    lattice.evolve()
+
+    if i % nStat == 0:
+        lattice.sync()
+        print("i = %4d; %3.0f MLUPS" % (i, MLUPS(lattice.geometry.volume, nStat, time.time() - lastStat)))
+        moments.append(lattice.get_moments())
+        lastStat = time.time()
+
+print("\nConcluded simulation.\n")
+
+generate_moment_plots(lattice, moments)