diff options
| -rw-r--r-- | implosion.py | 28 |
1 files changed, 20 insertions, 8 deletions
diff --git a/implosion.py b/implosion.py index 25d53c2..8a3e530 100644 --- a/implosion.py +++ b/implosion.py @@ -68,7 +68,8 @@ float equilibrium(float d, float2 v, int i, int j) { } __kernel void collide_and_stream(__global float* pop_a, - __global float* pop_b) + __global float* pop_b, + __global int* material) { const unsigned int gid = get_global_id(0); const uint2 cell = cellAtGid(gid); @@ -76,7 +77,7 @@ __kernel void collide_and_stream(__global float* pop_a, float d = density(&pop_b[gid*9]); float2 v = velocity(&pop_b[gid*9],d); - if ( cell.x >= 2 && cell.x <= $nX-3 && cell.y >= 2 && cell.y <= $nY-3 ) { + if ( material[gid] == 1 ) { for ( int i = -1; i <= 1; ++i ) { for ( int j = -1; j <= 1; ++j ) { pop_a[gidOfCell(cell.x+i, cell.y+j)*9 + indexOfDirection(i,j)] = @@ -84,9 +85,7 @@ __kernel void collide_and_stream(__global float* pop_a, } } } - else if ( ((cell.y == 1 || cell.y == $nY-2) && (cell.x > 0 && cell.x < $nX-1)) || - ((cell.x == 1 || cell.x == $nX-2) && (cell.y > 0 && cell.y < $nY-1)) ) - { + else if ( material[gid] == 2 ) { for ( int i = -1; i <= 1; ++i ) { for ( int j = -1; j <= 1; ++j ) { pop_a[gidOfCell(cell.x-i, cell.y-j)*9 + indexOfDirection(-i,-j)] = @@ -112,15 +111,28 @@ class D2Q9_BGK_Lattice: self.np_pop_a = numpy.ndarray(shape=(self.nCells, 9), dtype=numpy.float32) self.np_pop_b = numpy.ndarray(shape=(self.nCells, 9), dtype=numpy.float32) + self.np_material = numpy.ndarray(shape=(self.nCells, 1), dtype=numpy.int32) + + self.setup_geometry() self.equilibrilize() self.setup_anomaly() self.cl_pop_a = cl.Buffer(self.context, mf.READ_WRITE | mf.USE_HOST_PTR, hostbuf=self.np_pop_a) self.cl_pop_b = cl.Buffer(self.context, mf.READ_WRITE | mf.USE_HOST_PTR, hostbuf=self.np_pop_b) + self.cl_material = cl.Buffer(self.context, mf.READ_ONLY | mf.USE_HOST_PTR, hostbuf=self.np_material) self.build_kernel() + def setup_geometry(self): + self.np_material[:] = 0 + for x in range(1,self.nX-1): + for y in range(1,self.nY-1): + if x == 1 or y == 1 or x == self.nX-2 or y == self.nY-2: + self.np_material[self.idx(x,y)] = 2 + else: + self.np_material[self.idx(x,y)] = 1 + def equilibrilize(self): self.np_pop_a[:,:] = [ 1./36., 1./9., 1./36., 1./9. , 4./9., 1./9. , 1./36 , 1./9., 1./36. ] self.np_pop_b[:,:] = [ 1./36., 1./9., 1./36., 1./9. , 4./9., 1./9. , 1./36 , 1./9., 1./36. ] @@ -148,11 +160,11 @@ class D2Q9_BGK_Lattice: def evolve(self): if self.tick: self.tick = False - self.program.collide_and_stream(self.queue, (self.nCells,), None, self.cl_pop_a, self.cl_pop_b) + self.program.collide_and_stream(self.queue, (self.nCells,), None, self.cl_pop_a, self.cl_pop_b, self.cl_material) self.queue.finish() else: self.tick = True - self.program.collide_and_stream(self.queue, (self.nCells,), None, self.cl_pop_b, self.cl_pop_a) + self.program.collide_and_stream(self.queue, (self.nCells,), None, self.cl_pop_b, self.cl_pop_a, self.cl_material) self.queue.finish() def show(self, i): @@ -180,7 +192,7 @@ nUpdates = 100 start = timer() -for i in range(1,nUpdates): +for i in range(0,nUpdates): LBM.evolve() end = timer() |