tnavierstokes.cpp - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tnavierstokes.cpp (11524B)
---
1 #include <iostream>
2 #include <cstdio>
3 #include <cstdlib>
4 #include <string>
5 #include <vector>
6
7 #include "typedefs.h"
8 #include "datatypes.h"
9 #include "constants.h"
10 #include "sphere.h"
11 #include "utility.h"
12
13 // 1: Enable color output in array printing functions, 0: Disable
14 const int color_output = 0;
15
16 // Initialize memory
17 void DEM::initNSmem()
18 {
19 // Number of cells
20 ns.nx = grid.num[0];
21 ns.ny = grid.num[1];
22 ns.nz = grid.num[2];
23 unsigned int ncells = NScells();
24 unsigned int ncells_st = NScellsVelocity();
25
26 ns.p = new Float[ncells]; // hydraulic pressure
27 ns.v = new Float3[ncells]; // hydraulic velocity
28 ns.v_x = new Float[ncells_st]; // hydraulic velocity in staggered grid
29 ns.v_y = new Float[ncells_st]; // hydraulic velocity in staggered grid
30 ns.v_z = new Float[ncells_st]; // hydraulic velocity in staggered grid
31 //ns.v_p = new Float3[ncells]; // predicted hydraulic velocity
32 //ns.v_p_x = new Float[ncells_st]; // pred. hydraulic velocity in st. grid
33 //ns.v_p_y = new Float[ncells_st]; // pred. hydraulic velocity in st. grid
34 //ns.v_p_z = new Float[ncells_st]; // pred. hydraulic velocity in st. grid
35 ns.phi = new Float[ncells]; // porosity
36 ns.dphi = new Float[ncells]; // porosity change
37 ns.norm = new Float[ncells]; // normalized residual of epsilon
38 ns.epsilon = new Float[ncells]; // normalized residual of epsilon
39 ns.epsilon_new = new Float[ncells]; // normalized residual of epsilon
40 ns.f_d = new Float4[np]; // drag force on particles
41 ns.f_p = new Float4[np]; // pressure force on particles
42 ns.f_v = new Float4[np]; // viscous force on particles
43 ns.f_sum = new Float4[np]; // sum of fluid forces on particles
44 ns.p_constant = new int[ncells]; // unused
45 }
46
47 unsigned int DEM::NScells()
48 {
49 //return ns.nx*ns.ny*ns.nz; // without ghost nodes
50 return (ns.nx+2)*(ns.ny+2)*(ns.nz+2); // with ghost nodes
51 }
52
53 // Returns the number of velocity nodes in a congruent padded grid. There are
54 // velocity nodes between the boundary points and the pressure ghost nodes, but
55 // not on the outer side of the ghost nodes
56 unsigned int DEM::NScellsVelocity()
57 {
58 // Congruent padding for velocity grids. See "Cohen and Molemaker 'A fast
59 // double precision CFD code using CUDA'" for details
60 //return (ns.nx+1)*(ns.ny+1)*(ns.nz+1); // without ghost nodes
61 return (ns.nx+3)*(ns.ny+3)*(ns.nz+3); // with ghost nodes
62 }
63
64 // Free memory
65 void DEM::freeNSmem()
66 {
67 delete[] ns.p;
68 delete[] ns.v;
69 delete[] ns.v_x;
70 delete[] ns.v_y;
71 delete[] ns.v_z;
72 //delete[] ns.v_p;
73 //delete[] ns.v_p_x;
74 //delete[] ns.v_p_y;
75 //delete[] ns.v_p_z;
76 delete[] ns.phi;
77 delete[] ns.dphi;
78 delete[] ns.norm;
79 delete[] ns.epsilon;
80 delete[] ns.epsilon_new;
81 delete[] ns.f_d;
82 delete[] ns.f_p;
83 delete[] ns.f_v;
84 delete[] ns.f_sum;
85 delete[] ns.p_constant;
86 }
87
88 // 3D index to 1D index
89 unsigned int DEM::idx(
90 const int x,
91 const int y,
92 const int z)
93 {
94 // without ghost nodes
95 //return x + d.nx*y + d.nx*d.ny*z;
96
97 // with ghost nodes
98 // the ghost nodes are placed at x,y,z = -1 and WIDTH
99 return (x+1) + (ns.nx+2)*(y+1) + (ns.nx+2)*(ns.ny+2)*(z+1);
100 }
101
102 // 3D index to 1D index of cell-face velocity nodes. The cell-face velocities
103 // are placed at x = [0;nx], y = [0;ny], z = [0;nz].
104 // The coordinate x,y,z corresponds to the lowest corner of cell(x,y,z).
105 unsigned int DEM::vidx(
106 const int x,
107 const int y,
108 const int z)
109 {
110 //return x + (ns.nx+1)*y + (ns.nx+1)*(ns.ny+1)*z; // without ghost nodes
111 return (x+1) + (ns.nx+3)*(y+1) + (ns.nx+3)*(ns.ny+3)*(z+1); // with ghost nodes
112 }
113
114 // Determine if the FTCS (forward time, central space) solution of the Navier
115 // Stokes equations is unstable
116 void DEM::checkNSstability()
117 {
118 // Cell dimensions
119 const Float dx = grid.L[0]/grid.num[0];
120 const Float dy = grid.L[1]/grid.num[1];
121 const Float dz = grid.L[2]/grid.num[2];
122
123 // The smallest grid spacing
124 const Float dmin = fmin(dx, fmin(dy, dz));
125
126 // Check the diffusion term using von Neumann stability analysis
127 if (ns.mu*time.dt/(dmin*dmin) > 0.5) {
128 std::cerr << "Error: The time step is too large to ensure stability in "
129 "the diffusive term of the fluid momentum equation.\n"
130 "Decrease the viscosity, decrease the time step, and/or increase "
131 "the fluid grid cell size." << std::endl;
132 exit(1);
133 }
134
135 int x,y,z;
136 Float3 v;
137 for (x=0; x<ns.nx; ++x) {
138 for (y=0; y<ns.ny; ++y) {
139 for (z=0; z<ns.nz; ++z) {
140
141 v = ns.v[idx(x,y,z)];
142
143 // Check the advection term using the Courant-Friedrichs-Lewy
144 // condition
145 if (v.x*time.dt/dx + v.y*time.dt/dy + v.z*time.dt/dz > 1.0) {
146 std::cerr << "Error: The time step is too large to ensure "
147 "stability in the advective term of the fluid momentum "
148 "equation.\n"
149 "This is caused by too high fluid velocities. "
150 "You can try to decrease the time step, and/or "
151 "increase the fluid grid cell size.\n"
152 "v(" << x << ',' << y << ',' << z << ") = ["
153 << v.x << ',' << v.y << ',' << v.z << "] m/s"
154 << std::endl;
155 exit(1);
156 }
157 }
158 }
159 }
160
161
162
163 }
164
165 // Print array values to file stream (stdout, stderr, other file)
166 void DEM::printNSarray(FILE* stream, Float* arr)
167 {
168 int x, y, z;
169
170 // show ghost nodes
171 //for (z=-1; z<=ns.nz; z++) { // bottom to top
172 //for (z = ns.nz-1; z >= -1; z--) { // top to bottom
173 for (z = ns.nz; z >= -1; z--) { // top to bottom
174
175 fprintf(stream, "z = %d\n", z);
176
177 for (y=-1; y<=ns.ny; y++) {
178 for (x=-1; x<=ns.nx; x++) {
179
180 // hide ghost nodes
181 /*for (z=0; z<ns.nz; z++) {
182 for (y=0; y<ns.ny; y++) {
183 for (x=0; x<ns.nx; x++) {*/
184
185 if (x > -1 && x < ns.nx &&
186 y > -1 && y < ns.ny &&
187 z > -1 && z < ns.nz) {
188 fprintf(stream, "%f\t", arr[idx(x,y,z)]);
189 } else { // ghost node
190 if (color_output) {
191 fprintf(stream, "\x1b[30;1m%f\x1b[0m\t",
192 arr[idx(x,y,z)]);
193 } else {
194 fprintf(stream, "%f\t", arr[idx(x,y,z)]);
195 }
196 }
197 }
198 fprintf(stream, "\n");
199 }
200 fprintf(stream, "\n");
201 }
202 }
203
204 // Overload printNSarray to add optional description
205 void DEM::printNSarray(FILE* stream, Float* arr, std::string desc)
206 {
207 std::cout << "\n" << desc << ":\n";
208 printNSarray(stream, arr);
209 }
210
211 // Print array values to file stream (stdout, stderr, other file)
212 void DEM::printNSarray(FILE* stream, Float3* arr)
213 {
214 int x, y, z;
215 for (z=0; z<ns.nz; z++) {
216 for (y=0; y<ns.ny; y++) {
217 for (x=0; x<ns.nx; x++) {
218 fprintf(stream, "%f,%f,%f\t",
219 arr[idx(x,y,z)].x,
220 arr[idx(x,y,z)].y,
221 arr[idx(x,y,z)].z);
222 }
223 fprintf(stream, "\n");
224 }
225 fprintf(stream, "\n");
226 }
227 }
228
229 // Overload printNSarray to add optional description
230 void DEM::printNSarray(FILE* stream, Float3* arr, std::string desc)
231 {
232 std::cout << "\n" << desc << ":\n";
233 printNSarray(stream, arr);
234 }
235
236 // Returns the mean particle radius
237 Float DEM::meanRadius()
238 {
239 unsigned int i;
240 Float r_sum;
241 for (i=0; i<np; ++i)
242 r_sum += k.x[i].w;
243 return r_sum/((Float)np);
244 }
245
246 // Returns the average value of the normalized residuals
247 double DEM::avgNormResNS()
248 {
249 double norm_res_sum, norm_res;
250
251 // do not consider the values of the ghost nodes
252 for (int z=0; z<grid.num[2]; ++z) {
253 for (int y=0; y<grid.num[1]; ++y) {
254 for (int x=0; x<grid.num[0]; ++x) {
255 norm_res = static_cast<double>(ns.norm[idx(x,y,z)]);
256 if (norm_res != norm_res) {
257 std::cerr << "\nError: normalized residual is NaN ("
258 << norm_res << ") in cell "
259 << x << "," << y << "," << z << std::endl;
260 std::cerr << "\tt = " << time.current << ", iter = "
261 << int(time.current/time.dt) << std::endl;
262 std::cerr << "This often happens if the system has become "
263 "unstable." << std::endl;
264 exit(1);
265 }
266 norm_res_sum += norm_res;
267 }
268 }
269 }
270 return norm_res_sum/(grid.num[0]*grid.num[1]*grid.num[2]);
271 }
272
273
274 // Returns the average value of the normalized residuals
275 double DEM::maxNormResNS()
276 {
277 double max_norm_res = -1.0e9; // initialized to a small number
278 double norm_res;
279
280 // do not consider the values of the ghost nodes
281 for (int z=0; z<grid.num[2]; ++z) {
282 for (int y=0; y<grid.num[1]; ++y) {
283 for (int x=0; x<grid.num[0]; ++x) {
284 norm_res = fabs(static_cast<double>(ns.norm[idx(x,y,z)]));
285 if (norm_res != norm_res) {
286 std::cerr << "\nError: normalized residual is NaN ("
287 << norm_res << ") in cell "
288 << x << "," << y << "," << z << std::endl;
289 std::cerr << "\tt = " << time.current << ", iter = "
290 << int(time.current/time.dt) << std::endl;
291 std::cerr << "This often happens if the system has become "
292 "unstable." << std::endl;
293 exit(1);
294 }
295 if (norm_res > max_norm_res)
296 max_norm_res = norm_res;
297 }
298 }
299 }
300 return max_norm_res;
301 }
302
303 // Initialize fluid parameters
304 void DEM::initNS()
305 {
306 // Cell size
307 ns.dx = grid.L[0]/ns.nx;
308 ns.dy = grid.L[1]/ns.ny;
309 ns.dz = grid.L[2]/ns.nz;
310
311 if (verbose == 1) {
312 std::cout << " - Fluid grid dimensions: "
313 << ns.nx << "*"
314 << ns.ny << "*"
315 << ns.nz << std::endl;
316 std::cout << " - Fluid grid cell size: "
317 << ns.dx << "*"
318 << ns.dy << "*"
319 << ns.dz << std::endl;
320 }
321 }
322
323 // Write values in scalar field to file
324 void DEM::writeNSarray(Float* array, const char* filename)
325 {
326 FILE* file;
327 if ((file = fopen(filename,"w"))) {
328 printNSarray(file, array);
329 fclose(file);
330 } else {
331 fprintf(stderr, "Error, could not open %s.\n", filename);
332 }
333 }
334
335 // Write values in vector field to file
336 void DEM::writeNSarray(Float3* array, const char* filename)
337 {
338 FILE* file;
339 if ((file = fopen(filename,"w"))) {
340 printNSarray(file, array);
341 fclose(file);
342 } else {
343 fprintf(stderr, "Error, could not open %s.\n", filename);
344 }
345 }
346
347
348 // Print final heads and free memory
349 void DEM::endNS()
350 {
351 // Write arrays to stdout/text files for debugging
352 //writeNSarray(ns.phi, "ns_phi.txt");
353
354 //printNSarray(stdout, ns.K, "ns.K");
355 //printNSarray(stdout, ns.H, "ns.H");
356 //printNSarray(stdout, ns.H_new, "ns.H_new");
357 //printNSarray(stdout, ns.V, "ns.V");
358
359 freeNSmem();
360 }
361
362 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4