 |
tns2dfd.py - ns2dfd - 2D finite difference Navier Stokes solver for fluid dynam… |
 |
git clone git://src.adamsgaard.dk/ns2dfd |
 |
Log |
|
Files |
|
Refs |
|
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. Mo… |
|
15 simulation values are assigned default values upon initializatio… |
|
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… |
|
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 … |
|
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 =… |
|
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 (… |
|
106 solver. The solver is identical to the Gauss-Seidel method when … |
|
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 =… |
|
123 ''' |
|
124 Set the wall boundary conditions. The values correspond to the f… |
|
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 … |
|
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=… |
|
168 self.t_end = numpy.fromfile(fh, dtype=numpy.float64, count=… |
|
169 self.t_file = numpy.fromfile(fh, dtype=numpy.float64, count=… |
|
170 self.tau = numpy.fromfile(fh, dtype=numpy.float64, count=… |
|
171 |
|
172 self.itermax = numpy.fromfile(fh, dtype=numpy.int32, count=1) |
|
173 self.epsilon = numpy.fromfile(fh, dtype=numpy.float64, count… |
|
174 self.omega = numpy.fromfile(fh, dtype=numpy.float64, count… |
|
175 self.gamma = numpy.fromfile(fh, dtype=numpy.float64, count… |
|
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=… |
|
182 self.w_right = numpy.fromfile(fh, dtype=numpy.int32, count=… |
|
183 self.w_top = numpy.fromfile(fh, dtype=numpy.int32, count=… |
|
184 self.w_bottom = numpy.fromfile(fh, dtype=numpy.int32, count=… |
|
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, coun… |
|
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, coun… |
|
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, coun… |
|
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 s… |
|
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 de… |
|
275 The file name will be in the format ``<self.sid>.vti``. The vti … |
|
276 can be used for visualizing the fluid in ParaView. |
|
277 |
|
278 The fluid grid is visualized by opening the vti files, and press… |
|
279 "Apply" to import all fluid field properties. To visualize the s… |
|
280 fields, such as the pressure, the porosity, the porosity change … |
|
281 velocity magnitude, choose "Surface" or "Surface With Edges" as … |
|
282 "Representation". Choose the desired property as the "Coloring" … |
|
283 It may be desirable to show the color bar by pressing the "Show"… |
|
284 and "Rescale" to fit the color range limits to the current file.… |
|
285 coordinate system can be displayed by checking the "Show Axis" f… |
|
286 All adjustments by default require the "Apply" button to be pres… |
|
287 before regenerating the view. |
|
288 |
|
289 The fluid vector fields (e.g. the fluid velocity) can be visuali… |
|
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 … |
|
293 that "Arrow" is selected as the "Glyph type", and "Velocity" as … |
|
294 "Vectors" value. Adjust the "Maximum Number of Points" to be at … |
|
295 big as the number of fluid cells in the grid. Press "Apply" to v… |
|
296 the arrows. |
|
297 |
|
298 If several data files are generated for the same simulation (e.g… |
|
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 … |
|
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) |
