1 /*
2 * This is a small 4th-order Computational Fluid Dynamics solver. It
3 * illustrates multicomponent arrays (i.e. vector fields) and stencil
4 * objects. It is implemented in "C++Tran" style for simplicity.
5 *
6 * Incompressible flow is assumed.
7 *
8 * The simulation is of water in a closed cube. The water is initially
9 * quiescent. At a single layer near the bottom of the cube, a forcing
10 * function simulates the effect of a (highly idealized) propeller.
11 *
12 * Thanks to Stephen Smith (ACL/LANL) for help with the math.
13 */
14
15 /*
16 * This is a "work in progress". Things I am unhapppy about:
17 * - arrays, geometry, and boundary conditions are separate objects.
18 * This results in a certain amount of kludginess. For example,
19 * the stencil operators take geometry as a 2nd parameter.
20 * On the plus side, one geometry object is shared among many
21 * arrays.
22 * - The new stencil objects are nice, but declaring a new stencil
23 * object for every single equation seems like overkill. It would
24 * be nice if stencil operators could be "expression templatized"
25 * so that these equations could be moved to the appropriate place
26 * in the code. This probably wouldn't be too much work.
27 * - Stencil objects can't take scalar arguments, so all the scalars
28 * have to be globals. Big ugh. There is a fix for this,
29 * but it hasn't been implemented yet.
30 * - There are serious bugs, this simulation doesn't work very well.
31 */
32
33 #define NO_STIRRING
34 #undef NO_GRAVITY
35
36 #include <blitz/array-old.h>
37 #include <blitz/array/cgsolve.h>
38 #ifdef BZ_HAVE_STD
39 #include <fstream>
40 #else
41 #include <fstream.h>
42 #endif
43
44 using namespace blitz;
45
46 /*
47 * The current implementation of stencil objects forces these variables
48 * to be placed in global scope. Ugh. This restriction will be removed
49 * eventually.
50 */
51 double rho; // Density of fluid
52 double recip_rho; // 1/rho
53 double eta; // Kinematic viscosity
54 double time_now; // Elapsed seconds
55 double delta_t; // Time step
56 double volume; // Volume of a cell
57 double airPressure; // Air pressure (Pa)
58 double spatialStep; // Grid element size
59 double gravity; // Acceleration due to gravity
60 double gravityPressureGradient; // Pressure gradient due to gravity
61
62 /*
63 * The "geometry" object specifies how an array is mapped into real-world
64 * space. In this case, "UniformCubicGeometry" is used, which means that
65 * the real-world grid is orthogonal, regularly spaced, with the same spatial
66 * step in each dimension.
67 */
68
69 UniformCubicGeometry<3> geom; // Geometry
70
71 /*
72 * Some typedefs to make life easier.
73 */
74
75 typedef TinyVector<double,3> vector3d;
76 typedef Array<vector3d,3> vectorField;
77 typedef Array<double,3> scalarField;
78
79 /*
80 * Function prototypes
81 */
82
83 void record(vectorField& V);
84 void snapshot(scalarField&);
85 void snapshot(vectorField&);
86
87
88 /***********************************************************************
89 *
90 * Stencil objects involve these arrays:
91 *
92 * V Velocity (vector field)
93 * nextV Velocity field at next time step (vector field)
94 * advect Advection (vector field)
95 * P Pressure (scalar field)
96 * force Force (vector field)
97 * rhs Right-hand side of the pressure equation (scalar field)
98 *
99 * These stencil operators are used:
100 *
101 * grad3DVec4 4th-order gradient of a vector field (Jacobian matrix)
102 * Laplacian3DVec4 4th-order Laplacian of a vector field
103 * grad3D4 4th-order gradient of a scalar field
104 * div3DVec4 4th-order divergence of a vector field
105 *
106 ***********************************************************************/
107
108
109 /*********** Advection Calculation Stencil ************/
110
111 BZ_DECLARE_STENCIL2(calc_advect, advect, V)
112
113 advect = product(grad3DVec4(V,geom), *V);
114
115 BZ_END_STENCIL
116
117
118 /*********** Timestep the velocity field ************
119 * This is a 63-point stencil. For example, Laplacian3DVec4 turns into
120 * a 45-point stencil: each 2nd derivative is a 5-point stencil, and
121 * there are 9 of these derivatives to take the Laplacian of a 3D vector
122 * field.
123 */
124
125 BZ_DECLARE_STENCIL5(timestep, V, nextV, P, advect, force)
126
127 nextV = *V + delta_t * ( recip_rho * (
128 eta * Laplacian3DVec4(V,geom) - grad3D4(P, geom) + *force) - *advect);
129
130 BZ_END_STENCIL
131
132 /*********** Calculate the RHS of the pressure eqn ************/
133
134
135 BZ_DECLARE_STENCIL2(calc_P_rhs, rhs, advect)
136
137 rhs = - rho * div3DVec4(advect, geom);
138
139 BZ_END_STENCIL
140
141 /*********** Pressure equation update for the solver ************/
142
143 BZ_DECLARE_STENCIL2(P_solver_update, result, P)
144
145 result += Laplacian3D4(P,geom);
146
147 BZ_END_STENCIL
148
149
150 /*
151 * Allocate arrays and set their initial state
152 */
setup(const int N,vectorField & V,vectorField & nextV,scalarField & P,scalarField & P_rhs,vectorField & advect,vectorField & force)153 void setup(const int N, vectorField& V, vectorField& nextV, scalarField& P,
154 scalarField& P_rhs, vectorField& advect, vectorField& force)
155 {
156 // A 1m x 1m x 1m domain
157 spatialStep = 1.0 / (N - 1);
158 geom = UniformCubicGeometry<3>(spatialStep);
159
160 // Allocate arrays
161 allocateArrays(shape(N,N,N), advect, V, nextV, force); // vector fields
162 allocateArrays(shape(N,N,N), P, P_rhs); // scalar fields
163
164 // Since incompressibility is assumed, pressure only shows up as
165 // derivative terms in the equations. We choose airPressure = 0
166 // as an arbitrary datum.
167
168 airPressure = 0; // Pa
169 rho = 1000; // density of fluid, kg/m^3
170 recip_rho = 1.0 / rho; // inverse of density
171 eta = 1.0e-6; // kinematic viscosity of fluid, m^2/s
172 gravity = 9.81; // m/s^2
173 delta_t = 0.001; // initial time step, in seconds
174 volume = pow3(spatialStep); // cubic volume associated with grid point
175
176 // Kludge: Set eta high, so that the flow will spread faster.
177 // This means the cube is filled with molasses, rather than water.
178 eta *= 1000;
179
180 // Initial conditions: quiescent
181 V = 0.0;
182 P_rhs = 0.0;
183 advect = 0.0;
184 nextV = 0.0;
185
186 // Initial pressure condition: gravity causes a linear increase
187 // in pressure with depth. Note that tensor::k means the index
188 // associated with the z axis (they are labelled i, j, k).
189 #ifdef NO_GRAVITY
190 gravityPressureGradient = 0.0;
191 P = 0.0;
192 #else
193 gravityPressureGradient = spatialStep * gravity * rho;
194
195 P = airPressure + tensor::k * gravityPressureGradient;
196
197 #endif
198
199 // Set up the forcing function: gravity plus a stirring force
200 // at the bottom
201 double gravity_z = gravity * rho;
202
203 const int x = 0, y = 1, z = 2;
204 force[x] = 0.0;
205 force[y] = 0.0;
206 #ifdef NO_GRAVITY
207 force[z] = 0.0;
208 #else
209 force[z] = gravity_z;
210 #endif
211
212 #ifndef NO_STIRRING
213 // Centre of the stirring
214 int centrex = int(2 * N / 3.0);
215 int centrey = int(2 * N / 3.0);
216 int centrez = int(4 * N / 5.0);
217
218 const double stirRadius = 1.0 / 3.0;
219
220 vector3d zaxis(0,0,1);
221
222 // Loop through the 2D slice where the stirring occurs
223
224 for (int i=force.lbound(firstDim); i <= force.ubound(firstDim); ++i)
225 {
226 for (int j=force.lbound(secondDim); j <= force.ubound(secondDim); ++j)
227 {
228 // Vector from the centre of the stirring to the current
229 // coordinate
230
231 vector3d r((i-centrex) * spatialStep, (j-centrey) * spatialStep, 0.0);
232
233 if (norm(r) < stirRadius)
234 {
235 // The cross product of the z-axis and the vector3d to this
236 // coordinate yields the direction of the force. Multiply
237 // by gravity to get a reasonable magnitude (max force =
238 // 5 * gravity)
239 force(i,j,centrez) += cross(zaxis, r)
240 * (5 * gravity_z / stirRadius);
241 }
242 }
243 }
244 #endif // NO_STIRRING
245 }
246
247 // Calculate a simple check on a vector field
record(vectorField & V)248 void record(vectorField& V)
249 {
250 // Calculate the magnitude of a field
251 const int x=0, y=1, z=2;
252 double magx = sum(pow2(V[x])) / V.numElements();
253 double magy = sum(pow2(V[y])) / V.numElements();
254 double magz = sum(pow2(V[z])) / V.numElements();
255
256 cout << "norm = [" << magx
257 << " " << magy << " " << magz << " ]" << endl;
258 }
259
260 // Boundary conditions for the pressure field
261
262 class PressureBCs {
263 public:
PressureBCs(double gravityPressureGradient)264 PressureBCs(double gravityPressureGradient)
265 : gravityPressureGradient_(gravityPressureGradient)
266 { }
267
applyBCs(scalarField & P) const268 void applyBCs(scalarField& P) const
269 {
270 // Apply the Neumann boundary condition that grad(P) dot surface = 0
271 int xN = P.ubound(firstDim);
272 int yN = P.ubound(secondDim);
273 int zN = P.ubound(thirdDim);
274
275 Range all = Range::all();
276
277 // lower x
278 P(0,all,all) = P(2,all,all);
279 P(1,all,all) = P(2,all,all);
280
281 // upper x
282 P(xN,all,all) = P(xN-2,all,all);
283 P(xN-1,all,all) = P(xN-2,all,all);
284
285 // lower y
286 P(all,0,all) = P(all,2,all);
287 P(all,1,all) = P(all,2,all);
288
289 // upper y
290 P(all,yN,all) = P(all,yN-2,all);
291 P(all,yN-1,all) = P(all,yN-2,all);
292
293 // lower z
294 P(all,all,0) = P(all,all,2) - 2*gravityPressureGradient_;
295 P(all,all,1) = P(all,all,2) - gravityPressureGradient_;
296
297 // upper z
298 P(all,all,zN) = P(all,all,zN-2) + 2*gravityPressureGradient_;
299 P(all,all,zN-1) = P(all,all,zN-2) + gravityPressureGradient_;
300 }
301
302 private:
303 double gravityPressureGradient_;
304 };
305
306 /*
307 * One iteration: calculate advection, solve the pressure equation,
308 * then time step.
309 */
310
iterate(vectorField & V,vectorField & nextV,scalarField & P,scalarField & P_rhs,vectorField & advect,vectorField & force)311 void iterate(vectorField& V, vectorField& nextV, scalarField& P,
312 scalarField& P_rhs, vectorField& advect, vectorField& force)
313 {
314 // Calculate advection term
315 applyStencil(calc_advect(), advect, V);
316
317 // Solve to find the pressure
318 applyStencil(calc_P_rhs(), P_rhs, advect);
319 conjugateGradientSolver(P_solver_update(), P, P_rhs, 1e-8,
320 PressureBCs(gravityPressureGradient));
321
322 // Time step
323 applyStencil(timestep(), V, nextV, P, advect, force);
324 cycleArrays(V, nextV);
325
326 // Record some information about this time step
327 cout << "Velocity field: "; record(V);
328 }
329
330 /*
331 * Adjust the time step according to the CFL stability criterion
332 */
adjustTimeStep(vectorField & V)333 void adjustTimeStep(vectorField& V)
334 {
335 // Find maximum velocity magnitude
336 double maxV = 0.0;
337
338 // NEEDS_WORK: Blitz should provide a norm(vectorField) function.
339 // This is ugly.
340
341 for (int i=V.lbound(0); i <= V.ubound(0); ++i)
342 for (int j=V.lbound(1); j <= V.ubound(1); ++j)
343 for (int k=V.lbound(2); k <= V.ubound(2); ++k)
344 {
345 double normV = norm(V(i,j,k));
346 if (normV > maxV)
347 maxV = normV;
348 }
349
350 cout << "Maximum velocity is " << maxV << " m/s" << endl;
351
352 maxV += 1e-10; // Avoid divide-by-zero
353
354 // Steve K: need to have spatialStep^2 / diffusion constant
355 // diffusion constant = eta * recip_rho
356
357 delta_t = 0.3 * spatialStep / maxV;
358 const double maxTimeStep = 0.01;
359
360 if (delta_t > maxTimeStep)
361 delta_t = maxTimeStep;
362
363 cout << "Set time step to " << delta_t << " s" << endl;
364 }
365
366 /*
367 * Main program loop
368 */
369
main()370 int main()
371 {
372 vectorField V, nextV; // Velocity fields
373 scalarField P, P_rhs; // Pressure fields
374 vectorField advect; // Advection field
375 vectorField force; // Forcing function
376
377 const int N = 50; // Arrays are NxNxN
378
379 setup(N, V, nextV, P, P_rhs, advect, force);
380
381 const int nIters = 999;
382
383 // Stirring motion will stop at this time
384 const double forceTurnOff = 25000.0; // seconds
385
386 for (int i=0; i < nIters; ++i)
387 {
388 cout << "Iteration " << i << " Time = " << time_now << " s" << endl;
389
390 iterate(V, nextV, P, P_rhs, advect, force);
391
392 // Update the current time, turn off the force when we pass
393 // forceTurnOff
394
395 double oldtime_now = time_now;
396 time_now += delta_t;
397
398 if ((time_now > forceTurnOff) && (oldtime_now <= forceTurnOff))
399 force = 0.0;
400
401 // Dump some state
402 snapshot(V);
403 snapshot(P);
404
405 // Adjust the time step for the next iteration
406 adjustTimeStep(V);
407 }
408
409 return 0;
410 }
411
snapshot(scalarField & P)412 void snapshot(scalarField& P)
413 {
414 static int snapshotNum = 0;
415
416 ++snapshotNum;
417 char filename[128];
418 sprintf(filename, "pressure%03d.m", snapshotNum);
419
420 ofstream ofs(filename);
421 int N = P.length(firstDim);
422
423 int k = N/2;
424
425 ofs << "P" << snapshotNum << " = [ ";
426 for (int i=0; i < N; ++i)
427 {
428 for (int j=0; j < N; ++j)
429 {
430 float value = P(k,j,N-i-1);
431 ofs << value << " ";
432 }
433 if (i < N-1)
434 ofs << ";" << endl;
435 }
436 ofs << "];" << endl;
437 }
438
snapshot(vectorField & P)439 void snapshot(vectorField& P)
440 {
441 static int snapshotNum = 0;
442
443 ++snapshotNum;
444 char filename[128];
445 sprintf(filename, "velocity%03d.m", snapshotNum);
446
447 ofstream ofs(filename);
448 int N = P.length(firstDim);
449
450 int k = N/2;
451
452 ofs << "P" << snapshotNum << " = [ ";
453 for (int i=0; i < N; ++i)
454 {
455 for (int j=0; j < N; ++j)
456 {
457 float value = norm(P(k,j,N-i-1));
458 ofs << value << " ";
459 }
460 if (i < N-1)
461 ofs << ";" << endl;
462 }
463 ofs << "];" << endl;
464 }
465
466