tdarcy.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
---
tdarcy.cpp (11192B)
---
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::initDarcyMem()
18 {
19 // Number of cells
20 darcy.nx = grid.num[0];
21 darcy.ny = grid.num[1];
22 darcy.nz = grid.num[2];
23 unsigned int ncells = darcyCells();
24 //unsigned int ncells_st = darcyCellsVelocity();
25
26 darcy.p = new Float[ncells]; // hydraulic pressure
27 darcy.v = new Float3[ncells]; // hydraulic velocity
28 darcy.phi = new Float[ncells]; // porosity
29 darcy.dphi = new Float[ncells]; // porosity change
30 darcy.norm = new Float[ncells]; // normalized residual of epsilon
31 darcy.f_p = new Float4[np]; // pressure force on particles
32 darcy.k = new Float[ncells]; // hydraulic pressure
33 darcy.p_constant = new int[ncells]; // constant pressure (0: no, 1: yes)
34 }
35
36 unsigned int DEM::darcyCells()
37 {
38 //return darcy.nx*darcy.ny*darcy.nz; // without ghost nodes
39 return (darcy.nx+2)*(darcy.ny+2)*(darcy.nz+2); // with ghost nodes
40 }
41
42 // Returns the number of velocity nodes in a congruent padded grid. There are
43 // velocity nodes between the boundary points and the pressure ghost nodes, but
44 // not on the outer side of the ghost nodes
45 unsigned int DEM::darcyCellsVelocity()
46 {
47 // Congruent padding for velocity grids. See "Cohen and Molemaker 'A fast
48 // double precision CFD code using CUDA'" for details
49 //return (darcy.nx+1)*(darcy.ny+1)*(darcy.nz+1); // without ghost nodes
50 return (darcy.nx+3)*(darcy.ny+3)*(darcy.nz+3); // with ghost nodes
51 }
52
53 // Free memory
54 void DEM::freeDarcyMem()
55 {
56 delete[] darcy.p;
57 delete[] darcy.v;
58 delete[] darcy.phi;
59 delete[] darcy.dphi;
60 delete[] darcy.norm;
61 delete[] darcy.f_p;
62 delete[] darcy.k;
63 delete[] darcy.p_constant;
64 }
65
66 // 3D index to 1D index
67 unsigned int DEM::d_idx(
68 const int x,
69 const int y,
70 const int z)
71 {
72 // without ghost nodes
73 //return x + d.nx*y + d.nx*d.ny*z;
74
75 // with ghost nodes
76 // the ghost nodes are placed at x,y,z = -1 and WIDTH
77 return (x+1) + (darcy.nx+2)*(y+1) + (darcy.nx+2)*(darcy.ny+2)*(z+1);
78 }
79
80 // 3D index to 1D index of cell-face velocity nodes. The cell-face velocities
81 // are placed at x = [0;nx], y = [0;ny], z = [0;nz].
82 // The coordinate x,y,z corresponds to the lowest corner of cell(x,y,z).
83 unsigned int DEM::d_vidx(
84 const int x,
85 const int y,
86 const int z)
87 {
88 //return x + (darcy.nx+1)*y + (darcy.nx+1)*(darcy.ny+1)*z; // without ghost nodes
89 return (x+1) + (darcy.nx+3)*(y+1) + (darcy.nx+3)*(darcy.ny+3)*(z+1); // with ghost nodes
90 }
91
92 Float DEM::largestDarcyPermeability()
93 {
94 Float k_max = 0.0;
95 for (unsigned int z=0; z<grid.num[2]; z++)
96 for (unsigned int y=0; y<grid.num[1]; y++)
97 for (unsigned int x=0; x<grid.num[0]; x++)
98 if (darcy.k[d_idx(x,y,z)] > k_max)
99 k_max = darcy.k[d_idx(x,y,z)];
100 return k_max;
101 }
102
103 Float DEM::smallestDarcyPorosity()
104 {
105 Float phi_min = 10.0;
106 for (unsigned int z=0; z<grid.num[2]; z++)
107 for (unsigned int y=0; y<grid.num[1]; y++)
108 for (unsigned int x=0; x<grid.num[0]; x++)
109 if (darcy.phi[d_idx(x,y,z)] < phi_min)
110 phi_min = darcy.phi[d_idx(x,y,z)];
111 return phi_min;
112 }
113
114 // Component-wise max absolute velocities
115 Float3 DEM::largestDarcyVelocities()
116 {
117 Float3 v_max_abs = MAKE_FLOAT3(0.0, 0.0, 0.0);
118 Float3 v;
119 for (unsigned int z=0; z<grid.num[2]; z++)
120 for (unsigned int y=0; y<grid.num[1]; y++)
121 for (unsigned int x=0; x<grid.num[0]; x++) {
122 v = darcy.v[d_idx(x,y,z)];
123 if (v.x > v_max_abs.x)
124 v_max_abs.x = fabs(v.x);
125 if (v.y > v_max_abs.y)
126 v_max_abs.y = fabs(v.y);
127 if (v.z > v_max_abs.z)
128 v_max_abs.z = fabs(v.z);
129 }
130 return v_max_abs;
131 }
132
133 // Determine if the FTCS (forward time, central space) solution of the Navier
134 // Stokes equations is unstable
135 void DEM::checkDarcyStability()
136 {
137 // Cell dimensions
138 const Float dx = grid.L[0]/grid.num[0];
139 const Float dy = grid.L[1]/grid.num[1];
140 const Float dz = grid.L[2]/grid.num[2];
141
142 /*const Float alpha_max = largestDarcyPermeability()
143 /(darcy.beta_f*smallestDarcyPorosity()*darcy.mu);
144
145 // von Neumann stability analysis
146 if (time.dt >= 1.0/(2.0*alpha_max) *
147 1.0/(1.0/(dx*dx) + 1.0/(dy*dy) + 1.0/(dz*dz))) {
148 std::cerr
149 << "\nError: The time step is too large to ensure stability in "
150 "the diffusive term of the fluid momentum equation.\n"
151 "Decrease the time step, increase the fluid viscosity, increase "
152 "the fluid compressibility and/or increase "
153 "the fluid grid cell size." << std::endl;
154 //exit(1);
155 }*/
156
157 // Courant-Friedrichs-Lewy criteria
158 Float3 v_max_abs = largestDarcyVelocities();
159 if (v_max_abs.x*time.dt > dx ||
160 v_max_abs.y*time.dt > dy ||
161 v_max_abs.z*time.dt > dz) {
162 std::cerr
163 << "\nError: The time step is too large to ensure stability due to "
164 "large fluid velocities.\n v_max_abs = "
165 << v_max_abs.x << ", "
166 << v_max_abs.y << ", "
167 << v_max_abs.z <<
168 " m/s.\nDecrease the time step "
169 "and/or increase the fluid grid cell size." << std::endl;
170 }
171 }
172
173 // Print array values to file stream (stdout, stderr, other file)
174 void DEM::printDarcyArray(FILE* stream, Float* arr)
175 {
176 int x, y, z;
177
178 // show ghost nodes
179 for (z = darcy.nz; z >= -1; z--) { // top to bottom
180
181 fprintf(stream, "z = %d\n", z);
182
183 for (y=-1; y<=darcy.ny; y++) {
184 for (x=-1; x<=darcy.nx; x++) {
185
186 if (x > -1 && x < darcy.nx &&
187 y > -1 && y < darcy.ny &&
188 z > -1 && z < darcy.nz) {
189 fprintf(stream, "%f\t", arr[d_idx(x,y,z)]);
190 } else { // ghost node
191 if (color_output) {
192 fprintf(stream, "\x1b[30;1m%f\x1b[0m\t",
193 arr[d_idx(x,y,z)]);
194 } else {
195 fprintf(stream, "%f\t", arr[d_idx(x,y,z)]);
196 }
197 }
198 }
199 fprintf(stream, "\n");
200 }
201 fprintf(stream, "\n");
202 }
203 }
204
205 // Overload printDarcyArray to add optional description
206 void DEM::printDarcyArray(FILE* stream, Float* arr, std::string desc)
207 {
208 std::cout << "\n" << desc << ":\n";
209 printDarcyArray(stream, arr);
210 }
211
212 // Print array values to file stream (stdout, stderr, other file)
213 void DEM::printDarcyArray(FILE* stream, Float3* arr)
214 {
215 int x, y, z;
216 for (z=0; z<darcy.nz; z++) {
217 for (y=0; y<darcy.ny; y++) {
218 for (x=0; x<darcy.nx; x++) {
219 fprintf(stream, "%f,%f,%f\t",
220 arr[d_idx(x,y,z)].x,
221 arr[d_idx(x,y,z)].y,
222 arr[d_idx(x,y,z)].z);
223 }
224 fprintf(stream, "\n");
225 }
226 fprintf(stream, "\n");
227 }
228 }
229
230 // Overload printDarcyArray to add optional description
231 void DEM::printDarcyArray(FILE* stream, Float3* arr, std::string desc)
232 {
233 std::cout << "\n" << desc << ":\n";
234 printDarcyArray(stream, arr);
235 }
236
237 // Returns the average value of the normalized residuals
238 double DEM::avgNormResDarcy()
239 {
240 double norm_res_sum, norm_res;
241
242 // do not consider the values of the ghost nodes
243 for (int z=0; z<grid.num[2]; ++z) {
244 for (int y=0; y<grid.num[1]; ++y) {
245 for (int x=0; x<grid.num[0]; ++x) {
246 norm_res = static_cast<double>(darcy.norm[d_idx(x,y,z)]);
247 if (norm_res != norm_res) {
248 std::cerr << "\nError: normalized residual is NaN ("
249 << norm_res << ") in cell "
250 << x << "," << y << "," << z << std::endl;
251 std::cerr << "\tt = " << time.current << ", iter = "
252 << int(time.current/time.dt) << std::endl;
253 std::cerr << "This often happens if the system has become "
254 "unstable." << std::endl;
255 exit(1);
256 }
257 norm_res_sum += norm_res;
258 }
259 }
260 }
261 return norm_res_sum/(grid.num[0]*grid.num[1]*grid.num[2]);
262 }
263
264
265 // Returns the average value of the normalized residuals
266 double DEM::maxNormResDarcy()
267 {
268 double max_norm_res = -1.0e9; // initialized to a small number
269 double norm_res;
270
271 // do not consider the values of the ghost nodes
272 for (int z=0; z<grid.num[2]; ++z) {
273 for (int y=0; y<grid.num[1]; ++y) {
274 for (int x=0; x<grid.num[0]; ++x) {
275 norm_res = fabs(static_cast<double>(darcy.norm[d_idx(x,y,z)]));
276 if (norm_res != norm_res) {
277 std::cerr << "\nError: normalized residual is NaN ("
278 << norm_res << ") in cell "
279 << x << "," << y << "," << z << std::endl;
280 std::cerr << "\tt = " << time.current << ", iter = "
281 << int(time.current/time.dt) << std::endl;
282 std::cerr << "This often happens if the system has become "
283 "unstable." << std::endl;
284 exit(1);
285 }
286 if (norm_res > max_norm_res)
287 max_norm_res = norm_res;
288 }
289 }
290 }
291 return max_norm_res;
292 }
293
294 // Initialize fluid parameters
295 void DEM::initDarcy()
296 {
297 // Cell size
298 darcy.dx = grid.L[0]/darcy.nx;
299 darcy.dy = grid.L[1]/darcy.ny;
300 darcy.dz = grid.L[2]/darcy.nz;
301
302 if (verbose == 1) {
303 std::cout << " - Fluid grid dimensions: "
304 << darcy.nx << "*"
305 << darcy.ny << "*"
306 << darcy.nz << std::endl;
307 std::cout << " - Fluid grid cell size: "
308 << darcy.dx << "*"
309 << darcy.dy << "*"
310 << darcy.dz << std::endl;
311 }
312 }
313
314 // Write values in scalar field to file
315 void DEM::writeDarcyArray(Float* array, const char* filename)
316 {
317 FILE* file;
318 if ((file = fopen(filename,"w"))) {
319 printDarcyArray(file, array);
320 fclose(file);
321 } else {
322 fprintf(stderr, "Error, could not open %s.\n", filename);
323 }
324 }
325
326 // Write values in vector field to file
327 void DEM::writeDarcyArray(Float3* array, const char* filename)
328 {
329 FILE* file;
330 if ((file = fopen(filename,"w"))) {
331 printDarcyArray(file, array);
332 fclose(file);
333 } else {
334 fprintf(stderr, "Error, could not open %s.\n", filename);
335 }
336 }
337
338
339 // Print final heads and free memory
340 void DEM::endDarcy()
341 {
342 // Write arrays to stdout/text files for debugging
343 //writeDarcyArray(darcy.phi, "ns_phi.txt");
344
345 freeDarcyMem();
346 }
347
348 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4