From 42f4ae5f67f17ff37b3e95cab3c905668816eee8 Mon Sep 17 00:00:00 2001
From: Adrian Kummerlaender
Date: Mon, 28 Oct 2019 23:07:44 +0100
Subject: 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.
---
 ldc_2d.py | 94 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 94 insertions(+)
 create mode 100644 ldc_2d.py

(limited to 'ldc_2d.py')

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)
-- 
cgit v1.2.3