From 3d355b66fe231837e0051cf289d8a0f72dec798c Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Thu, 13 Jun 2019 21:45:49 +0200 Subject: Further the separation between descriptor and lattice --- D2Q9.py | 93 ------------------------------------------------- implosion.py | 9 ++--- lbm.py | 98 ++++++++++++++++++++++++++++++++++++++++++++++++++++ lid_driven_cavity.py | 3 +- template/kernel.mako | 8 ++--- 5 files changed, 109 insertions(+), 102 deletions(-) delete mode 100644 D2Q9.py create mode 100644 lbm.py diff --git a/D2Q9.py b/D2Q9.py deleted file mode 100644 index 27ae375..0000000 --- a/D2Q9.py +++ /dev/null @@ -1,93 +0,0 @@ -import pyopencl as cl -mf = cl.mem_flags - -import numpy - -import sympy -import symbolic.D2Q9 as D2Q9 - -from mako.template import Template - -class Lattice: - def idx(self, x, y): - return y * self.nX + x; - - def __init__(self, nX, nY, geometry, moments, collide, pop_eq_src = '', boundary_src = ''): - self.nX = nX - self.nY = nY - self.nCells = nX * nY - - self.moments = moments - self.collide = collide - - self.pop_eq_src = pop_eq_src - self.boundary_src = boundary_src - - self.platform = cl.get_platforms()[0] - self.context = cl.Context(properties=[(cl.context_properties.PLATFORM, self.platform)]) - self.queue = cl.CommandQueue(self.context) - - self.np_material = numpy.ndarray(shape=(self.nCells, 1), dtype=numpy.int32) - self.setup_geometry(geometry) - - self.tick = True - self.cl_pop_a = cl.Buffer(self.context, mf.READ_WRITE, size=9*self.nCells*numpy.float32(0).nbytes) - self.cl_pop_b = cl.Buffer(self.context, mf.READ_WRITE, size=9*self.nCells*numpy.float32(0).nbytes) - - self.cl_moments = cl.Buffer(self.context, mf.WRITE_ONLY, size=3*self.nCells*numpy.float32(0).nbytes) - self.cl_material = cl.Buffer(self.context, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=self.np_material) - - self.build_kernel() - - self.program.equilibrilize(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_pop_b).wait() - - def setup_geometry(self, geometry): - for y in range(1,self.nY-1): - for x in range(1,self.nX-1): - self.np_material[self.idx(x,y)] = geometry(self.nX,self.nY,x,y) - - def build_kernel(self): - program_src = Template(filename = './template/kernel.mako').render( - nX = self.nX, - nY = self.nY, - nCells = self.nCells, - moments_helper = self.moments[0], - moments_assignment = self.moments[1], - collide_helper = self.collide[0], - collide_assignment = self.collide[1], - c = D2Q9.c, - w = D2Q9.w, - ccode = sympy.ccode, - - pop_eq_src = Template(self.pop_eq_src).render( - nX = self.nX, - nY = self.nY, - nCells = self.nCells, - c = D2Q9.c, - w = D2Q9.w - ), - boundary_src = Template(self.boundary_src).render( - d = D2Q9.d - ) - ) - self.program = cl.Program(self.context, program_src).build('-cl-single-precision-constant') - - def evolve(self): - if self.tick: - self.tick = False - self.program.collide_and_stream(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_pop_b, self.cl_material) - else: - self.tick = True - self.program.collide_and_stream(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_b, self.cl_pop_a, self.cl_material) - - def sync(self): - self.queue.finish() - - def get_moments(self): - moments = numpy.ndarray(shape=(3, self.nCells), dtype=numpy.float32) - if self.tick: - self.program.collect_moments(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_b, self.cl_moments) - else: - self.program.collect_moments(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_moments) - cl.enqueue_copy(self.queue, moments, self.cl_moments).wait(); - return moments diff --git a/implosion.py b/implosion.py index fb101b4..a3bf8a1 100644 --- a/implosion.py +++ b/implosion.py @@ -5,7 +5,7 @@ import matplotlib import matplotlib.pyplot as plt matplotlib.use('AGG') -from D2Q9 import Lattice +from lbm import Lattice import symbolic.D2Q9 as D2Q9 @@ -33,12 +33,12 @@ def box(nX, nY, x, y): pop_eq = """ if ( sqrt(pow(get_global_id(0) - ${nX//2}.f, 2.f) + pow(get_global_id(1) - ${nY//2}.f, 2.f)) < ${nX//10} ) { -% for i, w_i in enumerate(w): +% for i, w_i in enumerate(descriptor.w): preshifted_f_a[${i*nCells}] = 1./24.f; preshifted_f_b[${i*nCells}] = 1./24.f; % endfor } else { -% for i, w_i in enumerate(w): +% for i, w_i in enumerate(descriptor.w): preshifted_f_a[${i*nCells}] = ${w_i}.f; preshifted_f_b[${i*nCells}] = ${w_i}.f; % endfor @@ -59,10 +59,11 @@ moments = [] print("Initializing simulation...\n") lattice = Lattice( + descriptor = D2Q9, nX = 1024, nY = 1024, - geometry = box, moments = D2Q9.moments(optimize = False), collide = D2Q9.bgk(tau = 0.8), + geometry = box, pop_eq_src = pop_eq, boundary_src = boundary) diff --git a/lbm.py b/lbm.py new file mode 100644 index 0000000..5118f17 --- /dev/null +++ b/lbm.py @@ -0,0 +1,98 @@ +import pyopencl as cl +mf = cl.mem_flags + +import numpy +import sympy + +from mako.template import Template + +class Lattice: + def __init__(self, descriptor, nX, nY, moments, collide, geometry, pop_eq_src = '', boundary_src = ''): + self.descriptor = descriptor + + self.nX = nX + self.nY = nY + self.nCells = nX * nY + + self.moments = moments + self.collide = collide + + self.pop_eq_src = pop_eq_src + self.boundary_src = boundary_src + + self.platform = cl.get_platforms()[0] + self.context = cl.Context(properties=[(cl.context_properties.PLATFORM, self.platform)]) + self.queue = cl.CommandQueue(self.context) + + self.np_material = numpy.ndarray(shape=(self.nCells, 1), dtype=numpy.int32) + self.setup_geometry(geometry) + + self.tick = True + + self.pop_size = descriptor.q * self.nCells * numpy.float32(0).nbytes + self.moments_size = (descriptor.d+1) * self.nCells * numpy.float32(0).nbytes + + self.cl_pop_a = cl.Buffer(self.context, mf.READ_WRITE, size=self.pop_size) + self.cl_pop_b = cl.Buffer(self.context, mf.READ_WRITE, size=self.pop_size) + + self.cl_moments = cl.Buffer(self.context, mf.WRITE_ONLY, size=self.moments_size) + self.cl_material = cl.Buffer(self.context, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=self.np_material) + + self.build_kernel() + + self.program.equilibrilize(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_pop_b).wait() + + def idx(self, x, y): + return y * self.nX + x; + + def setup_geometry(self, geometry): + for y in range(1,self.nY-1): + for x in range(1,self.nX-1): + self.np_material[self.idx(x,y)] = geometry(self.nX,self.nY,x,y) + + def build_kernel(self): + program_src = Template(filename = './template/kernel.mako').render( + descriptor = self.descriptor, + + nX = self.nX, + nY = self.nY, + nCells = self.nCells, + + moments_helper = self.moments[0], + moments_assignment = self.moments[1], + collide_helper = self.collide[0], + collide_assignment = self.collide[1], + + pop_eq_src = Template(self.pop_eq_src).render( + descriptor = self.descriptor, + nX = self.nX, + nY = self.nY, + nCells = self.nCells + ), + boundary_src = Template(self.boundary_src).render( + descriptor = self.descriptor + ), + + ccode = sympy.ccode + ) + self.program = cl.Program(self.context, program_src).build('-cl-single-precision-constant -cl-fast-relaxed-math') + + def evolve(self): + if self.tick: + self.tick = False + self.program.collide_and_stream(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_pop_b, self.cl_material) + else: + self.tick = True + self.program.collide_and_stream(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_b, self.cl_pop_a, self.cl_material) + + def sync(self): + self.queue.finish() + + def get_moments(self): + moments = numpy.ndarray(shape=(self.descriptor.d+1, self.nCells), dtype=numpy.float32) + if self.tick: + self.program.collect_moments(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_b, self.cl_moments) + else: + self.program.collect_moments(self.queue, (self.nX,self.nY), (32,1), self.cl_pop_a, self.cl_moments) + cl.enqueue_copy(self.queue, moments, self.cl_moments).wait(); + return moments diff --git a/lid_driven_cavity.py b/lid_driven_cavity.py index a4de67d..2bfdb7a 100644 --- a/lid_driven_cavity.py +++ b/lid_driven_cavity.py @@ -5,7 +5,7 @@ import matplotlib import matplotlib.pyplot as plt matplotlib.use('AGG') -from D2Q9 import Lattice +from lbm import Lattice import symbolic.D2Q9 as D2Q9 @@ -52,6 +52,7 @@ moments = [] print("Initializing simulation...\n") lattice = Lattice( + descriptor = D2Q9, nX = 256, nY = 256, geometry = cavity, moments = D2Q9.moments(optimize = False), diff --git a/template/kernel.mako b/template/kernel.mako index 6f82b20..24a6f64 100644 --- a/template/kernel.mako +++ b/template/kernel.mako @@ -7,7 +7,7 @@ __kernel void equilibrilize(__global __write_only float* f_a, __global __write_only float* preshifted_f_b = f_b + gid; % if pop_eq_src == '': -% for i, w_i in enumerate(w): +% for i, w_i in enumerate(descriptor.w): preshifted_f_a[${i*nCells}] = ${w_i}.f; preshifted_f_b[${i*nCells}] = ${w_i}.f; % endfor @@ -42,7 +42,7 @@ __kernel void collide_and_stream(__global __write_only float* f_a, __global __write_only float* preshifted_f_a = f_a + gid; __global __read_only float* preshifted_f_b = f_b + gid; -% for i, c_i in enumerate(c): +% for i, c_i in enumerate(descriptor.c): const float f_curr_${i} = preshifted_f_b[${direction_index(c_i)*nCells + neighbor_offset(-c_i)}]; % endfor @@ -64,7 +64,7 @@ __kernel void collide_and_stream(__global __write_only float* f_a, const float ${ccode(expr)} % endfor -% for i in range(0,len(c)): +% for i in range(0,descriptor.q): preshifted_f_a[${i*nCells}] = f_next_${i}; % endfor } @@ -76,7 +76,7 @@ __kernel void collect_moments(__global __read_only float* f, __global __read_only float* preshifted_f = f + gid; -% for i in range(0,len(c)): +% for i in range(0,descriptor.q): const float f_curr_${i} = preshifted_f[${i*nCells}]; % endfor -- cgit v1.2.3