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
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
|
import numpy
import time
import matplotlib
matplotlib.use('AGG')
import matplotlib.pyplot as plt
from simulation import Lattice, Geometry
from symbolic.generator import LBM
import symbolic.D3Q19 as D3Q19
from evtk.hl import imageToVTK
def MLUPS(cells, steps, time):
return cells * steps / time * 1e-6
def export_vtk(lattice, moments):
for i, m in enumerate(moments):
print("Export VTK file %d of %d." % (i+1, len(moments)))
imageToVTK("result/ldc_3d_%03d" % i, cellData = {
'velocity_x': m[1,:].reshape(lattice.geometry.size(), order = 'F'),
'velocity_y': m[2,:].reshape(lattice.geometry.size(), order = 'F'),
'velocity_z': m[3,:].reshape(lattice.geometry.size(), order = 'F')
})
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())))
# extract x-z-plane
y_slice = lattice.geometry.size_y//2
for z, x in numpy.ndindex(lattice.geometry.size_z-2, lattice.geometry.size_x-2):
gid = lattice.memory.gid(x+1,y_slice,z+1)
velocity[z,y_slice,x] = numpy.sqrt(m[1,gid]**2 + m[2,gid]**2 + m[3,gid]**2)
# extract y-z-plane
x_slice = lattice.geometry.size_x//2
for z, y in numpy.ndindex(lattice.geometry.size_z-2, lattice.geometry.size_y-2):
gid = lattice.memory.gid(x_slice,y+1,z+1)
velocity[z,y,x_slice] = numpy.sqrt(m[1,gid]**2 + m[2,gid]**2 + m[3,gid]**2)
plt.figure(figsize=(20, 10))
# plot x-z-plane
plt.subplot(1, 2, 1)
plt.imshow(velocity[:,y_slice,:], origin='lower', vmin=0.0, vmax=0.15, cmap=plt.get_cmap('seismic'))
# plot y-z-plane
plt.subplot(1, 2, 2)
plt.imshow(velocity[:,:,x_slice], origin='lower', vmin=0.0, vmax=0.15, cmap=plt.get_cmap('seismic'))
plt.savefig("result/ldc_3d_%02d.png" % i, bbox_inches='tight', pad_inches=0)
def get_cavity_material_map(geometry):
return [
(lambda x, y, z: x > 0 and x < geometry.size_x-1 and
y > 0 and y < geometry.size_y-1 and
z > 0 and z < geometry.size_z-1, 1), # bulk fluid
(lambda x, y, z: x == 1 or y == 1 or z == 1 or x == geometry.size_x-2 or y == geometry.size_y-2, 2), # walls
(lambda x, y, z: z == geometry.size_z-2, 3), # lid
(lambda x, y, z: x == 0 or x == geometry.size_x-1 or
y == 0 or y == geometry.size_y-1 or
z == 0 or z == geometry.size_z-1, 0) # ghost cells
]
boundary = """
if ( m == 2 ) {
u_0 = 0.0;
u_1 = 0.0;
u_2 = 0.0;
}
if ( m == 3 ) {
u_0 = 0.1;
u_1 = 0.0;
u_2 = 0.0;
}
"""
nUpdates = 10000
nStat = 1000
moments = []
print("Initializing simulation...\n")
lbm = LBM(D3Q19)
lattice = Lattice(
descriptor = D3Q19,
geometry = Geometry(64, 64, 64),
moments = lbm.moments(optimize = False),
collide = lbm.bgk(f_eq = lbm.equilibrium(), tau = 0.56),
boundary_src = boundary)
lattice.apply_material_map(
get_cavity_material_map(lattice.geometry))
lattice.sync_material()
print("Starting simulation using %d cells...\n" % lattice.geometry.volume)
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")
#export_vtk(lattice, moments)
generate_moment_plots(lattice, moments)
|
