1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
|
import numpy
import time
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)))
gid = lattice.memory.gid
velocity = numpy.reshape(
[ numpy.sqrt(m[gid(x,y)*3+1]**2 + m[gid(x,y)*3+2]**2) for x, y in lattice.geometry.inner_cells() ],
lattice.geometry.inner_size())
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 = 10000
geometry = Geometry(512, 512)
print("Generating kernel using boltzgen...\n")
functions = ['collide_and_stream', 'equilibrilize', 'collect_moments', 'momenta_boundary']
extras = ['cell_list_dispatch']
precision = 'single'
lbm = LBM(D2Q9)
generator = Generator(
descriptor = D2Q9,
moments = lbm.moments(),
collision = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.6),
target = 'cl',
precision = precision,
index = 'ZYX',
layout = 'SOA')
kernel_src = generator.kernel(geometry, functions, extras)
kernel_src += generator.custom(geometry, """
__kernel void equilibrilize(__global ${float_type}* f_next,
__global ${float_type}* f_prev)
{
const unsigned int gid = ${index.gid('get_global_id(0)', 'get_global_id(1)')};
equilibrilize_gid(f_next, f_prev, gid);
}
__kernel void collect_moments(__global ${float_type}* f,
__global ${float_type}* moments)
{
const unsigned int gid = ${index.gid('get_global_id(0)', 'get_global_id(1)')};
collect_moments_gid(f, moments, gid);
}
""")
print("Initializing simulation...\n")
lattice = Lattice(geometry, kernel_src, D2Q9, precision = precision)
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)
|
