tns2dfd.py - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynamics
(HTM) git clone git://src.adamsgaard.dk/ns2dfd
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tns2dfd.py (14344B)
---
1 #!/usr/bin/env python
2
3 import numpy
4 import subprocess
5 import vtk
6 import matplotlib
7 matplotlib.use('Agg')
8 import matplotlib.pyplot as plt
9
10 class fluid:
11
12 def __init__(self, name = 'unnamed'):
13 '''
14 A Navier-Stokes two-dimensional fluid flow simulation object. Most
15 simulation values are assigned default values upon initialization.
16 :param name: Simulation identifier
17 :type name: str
18 '''
19 self.sim_id = name
20
21 self.init_grid()
22 self.current_time()
23 self.end_time()
24 self.file_time()
25 self.safety_factor()
26 self.max_iterations()
27 self.tolerance_criteria()
28 self.relaxation_parameter()
29 self.upwind_differencing_factor()
30 self.boundary_conditions()
31 self.reynolds_number()
32 self.gravity()
33
34 def init_grid(self, nx = 10, ny = 10, dx = 0.1, dy = 0.1):
35 '''
36 Initializes the numerical grid.
37 :param nx: Fluid grid width in number of cells
38 :type nx: int
39 :param ny: Fluid grid height in number of cells
40 :type ny: int
41 :param dx: Grid cell width (meters)
42 :type dx: float
43 :param dy: Grid cell height (meters)
44 :type dy: float
45 '''
46 self.nx = numpy.asarray(nx)
47 self.ny = numpy.asarray(ny)
48 self.dx = numpy.asarray(dx)
49 self.dy = numpy.asarray(dy)
50 self.P = numpy.zeros((nx+2, ny+2))
51 self.U = numpy.zeros((nx+2, ny+2))
52 self.V = numpy.zeros((nx+2, ny+2))
53
54 def current_time(self, t = 0.0):
55 '''
56 Set the current simulation time. Default value = 0.0.
57 :param t: The current time value.
58 :type t: float
59 '''
60 self.t = numpy.asarray(t)
61
62 def end_time(self, t_end = 1.0):
63 '''
64 Set the simulation end time.
65 :param t_end: The time when to stop the simulation.
66 :type t_end: float
67 '''
68 self.t_end = numpy.asarray(t_end)
69
70 def file_time(self, t_file = 0.1):
71 '''
72 Set the simulation output file interval.
73 :param t_file: The time when to stop the simulation.
74 :type t_file: float
75 '''
76 self.t_file = numpy.asarray(t_file)
77
78 def safety_factor(self, tau = 0.5):
79 '''
80 Define the safety factor for the time step size control. Default value =
81 0.5.
82 :param tau: Safety factor in ]0;1]
83 :type tau: float
84 '''
85 self.tau = numpy.asarray(tau)
86
87 def max_iterations(self, itermax = 5000):
88 '''
89 Set the maximal allowed iterations per time step. Default value = 5000.
90 :param itermax: Max. solution iterations in [1;inf[
91 :type itermax: int
92 '''
93 self.itermax = numpy.asarray(itermax)
94
95 def tolerance_criteria(self, epsilon = 1.0e-4):
96 '''
97 Set the tolerance criteria for the fluid solver. Default value = 1.0e-4.
98 :param epsilon: Criteria value
99 :type epsilon: float
100 '''
101 self.epsilon = numpy.asarray(epsilon)
102
103 def relaxation_parameter(self, omega = 1.7):
104 '''
105 Set the relaxation parameter for the successive overrelaxation (SOR)
106 solver. The solver is identical to the Gauss-Seidel method when omega =
107 1. Default value = 1.7.
108 :param omega: Relaxation parameter value, in ]0;2[
109 :type omega: float
110 '''
111 self.omega = numpy.asarray(omega)
112
113 def upwind_differencing_factor(self, gamma = 0.9):
114 '''
115 Set the upwind diffencing factor used in the finite difference
116 approximations. Default value = 0.9.
117 :param gamma: Upward differencing factor value, in ]0;1[
118 :type gamma: float
119 '''
120 self.gamma = numpy.asarray(gamma)
121
122 def boundary_conditions(self, left = 1, right = 1, top = 1, bottom = 1):
123 '''
124 Set the wall boundary conditions. The values correspond to the following
125 conditions: 1) free-slip, 2) no-slip, 3) outflow, 4) periodic
126 :param left, right, top, bottom: The wall to specify the BC for
127 :type left, right, top, bottom: int
128 '''
129 self.w_left = numpy.asarray(left)
130 self.w_right = numpy.asarray(right)
131 self.w_top = numpy.asarray(top)
132 self.w_bottom = numpy.asarray(bottom)
133
134 def reynolds_number(self, re = 100):
135 '''
136 Define the simulation Reynolds number.
137 :param re: Reynolds number in ]0;infty[
138 :type re: float
139 '''
140 self.re = numpy.asarray(re)
141
142 def gravity(self, gx = 0.0, gy = 0.0):
143 '''
144 Set the gravitational acceleration on the fluid.
145 :param gx: Horizontal gravitational acceleration.
146 :type gx: float
147 :param gy: Vertical gravitational acceleration. Negative values are
148 downward.
149 :type gy: float
150 '''
151 self.gx = numpy.asarray(gx)
152 self.gy = numpy.asarray(gy)
153
154 def read(self, path, verbose = True):
155 '''
156 Read data file from disk.
157 :param path: Path to data file
158 :type path: str
159 '''
160 fh = None
161 try:
162 targetfile = path
163 if verbose == True:
164 print('Input file: ' + targetfile)
165 fh = open(targetfile, 'rb')
166
167 self.t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
168 self.t_end = numpy.fromfile(fh, dtype=numpy.float64, count=1)
169 self.t_file = numpy.fromfile(fh, dtype=numpy.float64, count=1)
170 self.tau = numpy.fromfile(fh, dtype=numpy.float64, count=1)
171
172 self.itermax = numpy.fromfile(fh, dtype=numpy.int32, count=1)
173 self.epsilon = numpy.fromfile(fh, dtype=numpy.float64, count=1)
174 self.omega = numpy.fromfile(fh, dtype=numpy.float64, count=1)
175 self.gamma = numpy.fromfile(fh, dtype=numpy.float64, count=1)
176
177 self.gx = numpy.fromfile(fh, dtype=numpy.float64, count=1)
178 self.gy = numpy.fromfile(fh, dtype=numpy.float64, count=1)
179 self.re = numpy.fromfile(fh, dtype=numpy.float64, count=1)
180
181 self.w_left = numpy.fromfile(fh, dtype=numpy.int32, count=1)
182 self.w_right = numpy.fromfile(fh, dtype=numpy.int32, count=1)
183 self.w_top = numpy.fromfile(fh, dtype=numpy.int32, count=1)
184 self.w_bottom = numpy.fromfile(fh, dtype=numpy.int32, count=1)
185
186 self.dx = numpy.fromfile(fh, dtype=numpy.float64, count=1)
187 self.dy = numpy.fromfile(fh, dtype=numpy.float64, count=1)
188 self.nx = numpy.fromfile(fh, dtype=numpy.int32, count=1)
189 self.ny = numpy.fromfile(fh, dtype=numpy.int32, count=1)
190
191 self.init_grid(dx = self.dx, dy = self.dy,\
192 nx = self.nx, ny = self.ny)
193
194 for i in range(self.nx+2):
195 for j in range(self.ny+2):
196 self.P[i,j] = \
197 numpy.fromfile(fh, dtype=numpy.float64, count=1)
198
199 for i in range(self.nx+2):
200 for j in range(self.ny+2):
201 self.U[i,j] = \
202 numpy.fromfile(fh, dtype=numpy.float64, count=1)
203
204 for i in range(self.nx+2):
205 for j in range(self.ny+2):
206 self.V[i,j] = \
207 numpy.fromfile(fh, dtype=numpy.float64, count=1)
208
209 finally:
210 if fh is not None:
211 fh.close()
212
213 def write(self, verbose = True, folder = './'):
214 '''
215 Write the simulation parameters to disk so that the fluid flow solver
216 can read it.
217 '''
218 fh = None
219 try:
220 targetfile = folder + '/' + self.sim_id + '.dat'
221 if verbose == True:
222 print('Output file: ' + targetfile)
223 fh = open(targetfile, 'wb')
224
225 fh.write(self.t.astype(numpy.float64))
226 fh.write(self.t_end.astype(numpy.float64))
227 fh.write(self.t_file.astype(numpy.float64))
228 fh.write(self.tau.astype(numpy.float64))
229
230 fh.write(self.itermax.astype(numpy.int32))
231 fh.write(self.epsilon.astype(numpy.float64))
232 fh.write(self.omega.astype(numpy.float64))
233 fh.write(self.gamma.astype(numpy.float64))
234
235 fh.write(self.gx.astype(numpy.float64))
236 fh.write(self.gy.astype(numpy.float64))
237 fh.write(self.re.astype(numpy.float64))
238
239 fh.write(self.w_left.astype(numpy.int32))
240 fh.write(self.w_right.astype(numpy.int32))
241 fh.write(self.w_top.astype(numpy.int32))
242 fh.write(self.w_bottom.astype(numpy.int32))
243
244 fh.write(self.dx.astype(numpy.float64))
245 fh.write(self.dy.astype(numpy.float64))
246 fh.write(self.nx.astype(numpy.int32))
247 fh.write(self.ny.astype(numpy.int32))
248
249 for i in range(self.nx+2):
250 for j in range(self.ny+2):
251 fh.write(self.P[i,j].astype(numpy.float64))
252
253 for i in range(self.nx+2):
254 for j in range(self.ny+2):
255 fh.write(self.U[i,j].astype(numpy.float64))
256
257 for i in range(self.nx+2):
258 for j in range(self.ny+2):
259 fh.write(self.V[i,j].astype(numpy.float64))
260
261 finally:
262 if fh is not None:
263 fh.close()
264
265 def run(self):
266 '''
267 Run the simulation using the C program.
268 '''
269 self.write()
270 subprocess.call('./ns2dfd ' + self.sim_id + '.dat', shell=True)
271
272 def writeVTK(self, folder = './', verbose = True):
273 '''
274 Writes a VTK file for the fluid grid to the current folder by default.
275 The file name will be in the format ``<self.sid>.vti``. The vti files
276 can be used for visualizing the fluid in ParaView.
277
278 The fluid grid is visualized by opening the vti files, and pressing
279 "Apply" to import all fluid field properties. To visualize the scalar
280 fields, such as the pressure, the porosity, the porosity change or the
281 velocity magnitude, choose "Surface" or "Surface With Edges" as the
282 "Representation". Choose the desired property as the "Coloring" field.
283 It may be desirable to show the color bar by pressing the "Show" button,
284 and "Rescale" to fit the color range limits to the current file. The
285 coordinate system can be displayed by checking the "Show Axis" field.
286 All adjustments by default require the "Apply" button to be pressed
287 before regenerating the view.
288
289 The fluid vector fields (e.g. the fluid velocity) can be visualizing by
290 e.g. arrows. To do this, select the fluid data in the "Pipeline
291 Browser". Press "Glyph" from the "Common" toolbar, or go to the
292 "Filters" mennu, and press "Glyph" from the "Common" list. Make sure
293 that "Arrow" is selected as the "Glyph type", and "Velocity" as the
294 "Vectors" value. Adjust the "Maximum Number of Points" to be at least as
295 big as the number of fluid cells in the grid. Press "Apply" to visualize
296 the arrows.
297
298 If several data files are generated for the same simulation (e.g. using
299 the :func:`writeVTKall()` function), it is able to step the
300 visualization through time by using the ParaView controls.
301
302 :param folder: The folder where to place the output binary file (default
303 (default = './')
304 :type folder: str
305 :param verbose: Show diagnostic information (default = True)
306 :type verbose: bool
307 '''
308 filename = folder + '/' + self.sim_id + '.vti' # image grid
309
310 # initalize VTK data structure
311 grid = vtk.vtkImageData()
312 grid.SetOrigin([0.0, 0.0, 0.0])
313 grid.SetSpacing([self.dx, self.dy, 1])
314 grid.SetDimensions([self.nx+2, self.ny+2, 1])
315
316 # array of scalars: hydraulic pressures
317 pres = vtk.vtkDoubleArray()
318 pres.SetName("Pressure")
319 pres.SetNumberOfComponents(1)
320 pres.SetNumberOfTuples(grid.GetNumberOfPoints())
321
322 # array of vectors: hydraulic velocities
323 vel = vtk.vtkDoubleArray()
324 vel.SetName("Velocity")
325 vel.SetNumberOfComponents(2)
326 vel.SetNumberOfTuples(grid.GetNumberOfPoints())
327
328 # insert values
329 for y in range(self.ny+2):
330 for x in range(self.nx+2):
331 idx = x + (self.nx+2)*y
332 pres.SetValue(idx, self.P[x,y])
333 vel.SetTuple(idx, [self.U[x,y], self.V[x,y]])
334
335 # add pres array to grid
336 grid.GetPointData().AddArray(pres)
337 grid.GetPointData().AddArray(vel)
338
339 # write VTK XML image data file
340 writer = vtk.vtkXMLImageDataWriter()
341 writer.SetFileName(filename)
342 writer.SetInput(grid)
343 writer.Update()
344 if (verbose == True):
345 print('Output file: {0}'.format(filename))
346
347 def plot_PUV(self, format = 'png'):
348 plt.figure(figsize=[8,8])
349
350 #ax = plt.subplot(1, 3, 1)
351 plt.title("Pressure")
352 imgplt = plt.imshow(self.P.T, origin='lower')
353 imgplt.set_interpolation('nearest')
354 #imgplt.set_interpolation('bicubic')
355 #imgplt.set_cmap('hot')
356 plt.xlabel('$x$')
357 plt.ylabel('$y$')
358 plt.colorbar()
359
360 # show velocities as arrows
361 Q = plt.quiver(self.U, self.V)
362
363 # show velocities as stream lines
364 #plt.streamplot(numpy.arange(self.nx+2),numpy.arange(self.ny+2),\
365 #self.U, self.V)
366
367 '''
368 # show velocities as heat maps
369 ax = plt.subplot(1, 3, 2)
370 plt.title("U")
371 imgplt = plt.imshow(self.U.T, origin='lower')
372 imgplt.set_interpolation('nearest')
373 #imgplt.set_interpolation('bicubic')
374 #imgplt.set_cmap('hot')
375 plt.xlabel('$x$')
376 plt.ylabel('$y$')
377 plt.colorbar()
378
379 ax = plt.subplot(1, 3, 3)
380 plt.title("V")
381 imgplt = plt.imshow(self.V.T, origin='lower')
382 imgplt.set_interpolation('nearest')
383 #imgplt.set_interpolation('bicubic')
384 #imgplt.set_cmap('hot')
385 plt.xlabel('$x$')
386 plt.ylabel('$y$')
387 plt.colorbar()
388 '''
389
390 plt.savefig(self.sim_id + '-PUV.' + format, transparent=False)