device.cu - 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
---
device.cu (121387B)
---
1 // device.cu -- GPU specific operations utilizing the CUDA API.
2 #include <iostream>
3 #include <fstream>
4 #include <string>
5 #include <cstdio>
6 #include <iomanip>
7 #include <time.h>
8
9 #ifdef SPHERE_GPU
10 #include <cuda.h>
11 #include <helper_math.h>
12 #include "thrust/device_ptr.h"
13 #include "thrust/sort.h"
14 #else
15 #include <algorithm>
16 #include <utility>
17 #include <vector>
18 #include "cpu_backend.h"
19 #endif
20
21 #include "vector_arithmetic.h" // for arbitrary prec. vectors
22 //#include <vector_functions.h> // for single prec. vectors
23 #include "launch.h"
24
25 #include "sphere.h"
26 #include "datatypes.h"
27 #include "utility.h"
28 #include "constants.cuh"
29 #include "debug.h"
30 #include "version.h"
31
32 #include "sorting.cuh"
33 #include "contactmodels.cuh"
34 #include "cohesion.cuh"
35 #include "contactsearch.cuh"
36 #include "integration.cuh"
37 #include "raytracer.cuh"
38 #include "navierstokes.cuh"
39 #include "darcy.cuh"
40
41 #ifdef SPHERE_GPU
42 // Returns the number of cores per streaming multiprocessor, which is
43 // a function of the device compute capability
44 // https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities
45 int cudaCoresPerSM(int major, int minor)
46 {
47 if (major == 1)
48 return 8;
49 else if (major == 2 && minor == 0)
50 return 32;
51 else if (major == 2 && minor == 1)
52 return 48;
53 else if (major == 3)
54 return 192;
55 else if (major == 4)
56 return 128;
57 else if (major == 5)
58 return 128;
59 else if (major == 6 && minor == 0)
60 return 64;
61 else if (major == 6 && minor == 1)
62 return 128;
63 else if (major == 6 && minor == 2)
64 return 128;
65 else if (major == 7)
66 return 32;
67 else if (major == 8)
68 return 64;
69 else
70 printf("Error in cudaCoresPerSM Device compute capability value "
71 "(%d.%d) not recognized.", major, minor);
72 return -1;
73 }
74
75 // Wrapper function for initializing the CUDA components.
76 // Called from main.cpp
77 void DEM::initializeGPU(void)
78 {
79 using std::cout; // stdout
80
81 // Specify target device
82 int cudadevice = 0;
83
84 // Variables containing device properties
85 cudaDeviceProp prop;
86 int deviceCount;
87 int cudaDriverVersion;
88 int cudaRuntimeVersion;
89
90 checkForCudaErrors("Before initializing CUDA device");
91
92 // Register number of devices
93 cudaGetDeviceCount(&deviceCount);
94 ndevices = deviceCount; // store in DEM class
95
96 if (deviceCount == 0) {
97 std::cerr << "\nERROR: No CUDA-enabled devices availible. Bye."
98 << std::endl;
99 exit(EXIT_FAILURE);
100 } else if (deviceCount == 1) {
101 if (verbose == 1)
102 cout << " System contains 1 CUDA compatible device.\n";
103 } else {
104 if (verbose == 1)
105 cout << " System contains " << deviceCount
106 << " CUDA compatible devices.\n";
107 }
108
109 // Loop through GPU's and choose the one with the most CUDA cores
110 if (device == -1) {
111 int ncudacores;
112 int max_ncudacores = 0;
113 for (int d=0; d<ndevices; d++) {
114 cudaGetDeviceProperties(&prop, d);
115 cudaDriverGetVersion(&cudaDriverVersion);
116 cudaRuntimeGetVersion(&cudaRuntimeVersion);
117
118 ncudacores = prop.multiProcessorCount
119 *cudaCoresPerSM(prop.major, prop.minor);
120 if (ncudacores > max_ncudacores) {
121 max_ncudacores = ncudacores;
122 cudadevice = d;
123 }
124
125 if (verbose == 1) {
126 cout << " CUDA device ID: " << d << "\n";
127 cout << " - Name: " << prop.name << ", compute capability: "
128 << prop.major << "." << prop.minor << ".\n";
129 cout << " - CUDA Driver version: " << cudaDriverVersion/1000
130 << "." << cudaDriverVersion%100
131 << ", runtime version " << cudaRuntimeVersion/1000 << "."
132 << cudaRuntimeVersion%100 << std::endl;
133 }
134 }
135
136 device = cudadevice; // store in DEM class
137 if (verbose == 1) {
138 cout << " Using CUDA device ID " << device << " with "
139 << max_ncudacores << " cores." << std::endl;
140 }
141
142 } else {
143
144 cudaGetDeviceProperties(&prop, device);
145 cudaDriverGetVersion(&cudaDriverVersion);
146 cudaRuntimeGetVersion(&cudaRuntimeVersion);
147
148 int ncudacores = prop.multiProcessorCount
149 *cudaCoresPerSM(prop.major, prop.minor);
150
151 if (verbose == 1) {
152 cout << " CUDA device ID: " << device << "\n";
153 cout << " - Name: " << prop.name << ", compute capability: "
154 << prop.major << "." << prop.minor << ".\n";
155 cout << " - CUDA Driver version: " << cudaDriverVersion/1000
156 << "." << cudaDriverVersion%100
157 << ", runtime version " << cudaRuntimeVersion/1000 << "."
158 << cudaRuntimeVersion%100
159 << "\n - " << ncudacores << " CUDA cores" << std::endl;
160 }
161 }
162
163 // The value of device is now 0 or larger
164 cudaSetDevice(device);
165
166 checkForCudaErrors("While initializing CUDA device");
167 }
168 #endif // SPHERE_GPU
169
170 // Start timer for kernel profiling
171 void startTimer(cudaEvent_t* kernel_tic)
172 {
173 cudaEventRecord(*kernel_tic);
174 }
175
176 // Stop timer for kernel profiling and time to function sum
177 void stopTimer(cudaEvent_t *kernel_tic,
178 cudaEvent_t *kernel_toc,
179 float *kernel_elapsed,
180 double* sum)
181 {
182 cudaEventRecord(*kernel_toc, 0);
183 cudaEventSynchronize(*kernel_toc);
184 cudaEventElapsedTime(kernel_elapsed, *kernel_tic, *kernel_toc);
185 *sum += *kernel_elapsed;
186 }
187
188 // Check values of parameters in constant memory
189 __global__ void checkConstantValues(int* dev_equal,
190 Grid* dev_grid,
191 Params* dev_params)
192 {
193 // Values ok (0)
194 *dev_equal = 0;
195
196 // Compare values between global- and constant
197 // memory structures
198 if (dev_grid->origo[0] != devC_grid.origo[0])
199 *dev_equal = 1;
200 if (dev_grid->origo[1] != devC_grid.origo[1])
201 *dev_equal = 2; // Not ok
202 if (dev_grid->origo[2] != devC_grid.origo[2])
203 *dev_equal = 3; // Not ok
204 if (dev_grid->L[0] != devC_grid.L[0])
205 *dev_equal = 4; // Not ok
206 if (dev_grid->L[1] != devC_grid.L[1])
207 *dev_equal = 5; // Not ok
208 if (dev_grid->L[2] != devC_grid.L[2])
209 *dev_equal = 6; // Not ok
210 if (dev_grid->num[0] != devC_grid.num[0])
211 *dev_equal = 7; // Not ok
212 if (dev_grid->num[1] != devC_grid.num[1])
213 *dev_equal = 8; // Not ok
214 if (dev_grid->num[2] != devC_grid.num[2])
215 *dev_equal = 9; // Not ok
216 if (dev_grid->periodic != devC_grid.periodic)
217 *dev_equal = 10; // Not ok
218
219 if (dev_params->g[0] != devC_params.g[0])
220 *dev_equal = 11; // Not ok
221 if (dev_params->g[1] != devC_params.g[1])
222 *dev_equal = 12; // Not ok
223 if (dev_params->g[2] != devC_params.g[2])
224 *dev_equal = 13; // Not ok
225 if (dev_params->k_n != devC_params.k_n)
226 *dev_equal = 14; // Not ok
227 if (dev_params->k_t != devC_params.k_t)
228 *dev_equal = 15; // Not ok
229 if (dev_params->k_r != devC_params.k_r)
230 *dev_equal = 16; // Not ok
231 if (dev_params->gamma_n != devC_params.gamma_n)
232 *dev_equal = 17; // Not ok
233 if (dev_params->gamma_t != devC_params.gamma_t)
234 *dev_equal = 18; // Not ok
235 if (dev_params->gamma_r != devC_params.gamma_r)
236 *dev_equal = 19; // Not ok
237 if (dev_params->mu_s != devC_params.mu_s)
238 *dev_equal = 20; // Not ok
239 if (dev_params->mu_d != devC_params.mu_d)
240 *dev_equal = 21; // Not ok
241 if (dev_params->mu_r != devC_params.mu_r)
242 *dev_equal = 22; // Not ok
243 if (dev_params->rho != devC_params.rho)
244 *dev_equal = 23; // Not ok
245 if (dev_params->contactmodel != devC_params.contactmodel)
246 *dev_equal = 24; // Not ok
247 if (dev_params->kappa != devC_params.kappa)
248 *dev_equal = 25; // Not ok
249 if (dev_params->db != devC_params.db)
250 *dev_equal = 26; // Not ok
251 if (dev_params->V_b != devC_params.V_b)
252 *dev_equal = 27; // Not ok
253 if (dev_params->lambda_bar != devC_params.lambda_bar)
254 *dev_equal = 28; // Not ok
255 if (dev_params->nb0 != devC_params.nb0)
256 *dev_equal = 29; // Not ok
257 if (dev_params->E != devC_params.E)
258 *dev_equal = 30; // Not ok
259 }
260
261 __global__ void checkParticlePositions(
262 const Float4* __restrict__ dev_x)
263 {
264 unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x; // Thread id
265
266 if (idx < devC_np) { // Condition prevents block size error
267 Float4 x = dev_x[idx];
268
269 // make sure grain doesn't have NaN or Inf position
270 if (!isfinite(x.x) || !isfinite(x.y) || !isfinite(x.z)) {
271 __syncthreads();
272 printf("\nParticle %d has non-finite position: x = %f %f %f",
273 idx, x.x, x.y, x.z);
274 }
275
276 /*__syncthreads();
277 printf("\nParticle %d: x = %f %f %f",
278 idx, x.x, x.y, x.z);*/
279
280 // check that the particle is inside of the simulation domain
281 if (x.x < devC_grid.origo[0] ||
282 x.y < devC_grid.origo[1] ||
283 x.z < devC_grid.origo[2] ||
284 x.x > devC_grid.L[0] ||
285 x.y > devC_grid.L[1] ||
286 x.z > devC_grid.L[2]) {
287 __syncthreads();
288 printf("\nParticle %d is outside the computational domain "
289 "(%f %f %f to %f %f %f): x = %f %f %f",
290 idx,
291 devC_grid.origo[0], devC_grid.origo[1], devC_grid.origo[2],
292 devC_grid.L[0], devC_grid.L[1], devC_grid.L[2],
293 x.x, x.y, x.z);
294 }
295 }
296 }
297
298
299 // Copy the constant data components to device memory,
300 // and check whether the values correspond to the
301 // values in constant memory.
302 void DEM::checkConstantMemory()
303 {
304 // Allocate space in global device memory
305 Grid* dev_grid;
306 Params* dev_params;
307 cudaMalloc((void**)&dev_grid, sizeof(Grid));
308 cudaMalloc((void**)&dev_params, sizeof(Params));
309
310 // Copy structure data from host to global device memory
311 cudaMemcpy(dev_grid, &grid, sizeof(Grid), cudaMemcpyHostToDevice);
312 cudaMemcpy(dev_params, ¶ms, sizeof(Params), cudaMemcpyHostToDevice);
313
314 // Compare values between global and constant memory
315 // structures on the device.
316 int* equal = new int; // The values are equal = 0, if not = 1
317 *equal = 0;
318 int* dev_equal;
319 cudaMalloc((void**)&dev_equal, sizeof(int));
320 KERNEL_LAUNCH(checkConstantValues, 1, 1, 0,
321 (dev_equal, dev_grid, dev_params));
322 checkForCudaErrors("After constant memory check");
323
324 // Copy result to host
325 cudaMemcpy(equal, dev_equal, sizeof(int), cudaMemcpyDeviceToHost);
326
327 // Free global device memory
328 cudaFree(dev_grid);
329 cudaFree(dev_params);
330 cudaFree(dev_equal);
331
332 // Are the values equal?
333 if (*equal != 0) {
334 std::cerr << "Error! The values in constant memory do not "
335 << "seem to be correct (" << *equal << ")." << std::endl;
336 exit(1);
337 } else {
338 if (verbose == 1)
339 std::cout << " Constant values ok (" << *equal << ")."
340 << std::endl;
341 }
342 }
343
344 // Copy selected constant components to constant device memory.
345 void DEM::transferToConstantDeviceMemory()
346 {
347 using std::cout;
348
349 if (verbose == 1)
350 cout << " Transfering data to constant device memory: ";
351
352 /*for (int d=0; d<ndevices; d++) {
353 cudaSetDevice(d);*/
354 cudaMemcpyToSymbol(devC_nd, &nd, sizeof(nd));
355 cudaMemcpyToSymbol(devC_np, &np, sizeof(np));
356 cudaMemcpyToSymbol(devC_nw, &walls.nw, sizeof(unsigned int));
357 cudaMemcpyToSymbol(devC_nc, &NC, sizeof(int));
358 cudaMemcpyToSymbol(devC_dt, &time.dt, sizeof(Float));
359 cudaMemcpyToSymbol(devC_grid, &grid, sizeof(Grid));
360 cudaMemcpyToSymbol(devC_params, ¶ms, sizeof(Params));
361 /*}
362 cudaSetDevice(device);*/
363
364 checkForCudaErrors("After transferring to device constant memory");
365
366 if (verbose == 1)
367 cout << "Done\n";
368
369 // only for device with most CUDA cores
370 checkConstantMemory();
371 }
372
373 __global__ void printWorldSize(Float4* dev_walls_nx)
374 {
375 printf("\nL = %f, %f, %f\n",
376 devC_grid.L[0], devC_grid.L[1], devC_grid.L[2]);
377 printf("\ndev_walls_nx[0] = %f, %f, %f, %f\n",
378 dev_walls_nx[0].x,
379 dev_walls_nx[0].y,
380 dev_walls_nx[0].z,
381 dev_walls_nx[0].w);
382 }
383
384 void DEM::updateGridSize()
385 {
386 //printf("\nDEM::updateGridSize() start\n");
387 Float* Lz = new Float;
388
389 // Get top wall position from dev_walls_nx[0].z
390 cudaMemcpy(Lz, &dev_walls_nx[0].w, sizeof(Float), cudaMemcpyDeviceToHost);
391 checkForCudaErrors("DEM::updateGridSize(): copying wall position");
392
393 //printWorldSize<<<1,1>>>(dev_walls_nx);
394 //cudaDeviceSynchronize();
395 //checkForCudaErrors("DEM::updateGridSize(): first printWorldSize");
396
397 //printf("\nLz = %f\n", *Lz);
398
399 // Write value to grid.L[2]
400 grid.L[2] = *Lz;
401
402 // Write value to devC_grid.L[2]
403 //cudaMemcpyToSymbol(devC_grid.L[2], &Lz, sizeof(Float));
404 cudaMemcpyToSymbol(devC_grid, &grid, sizeof(Grid));
405
406 checkForCudaErrors("DEM::updateGridSize(): write to devC_grid.L[2]");
407
408 //printWorldSize<<<1,1>>>(dev_walls_nx);
409 //cudaDeviceSynchronize();
410 //checkForCudaErrors("DEM::updateGridSize(): second printWorldSize");
411
412 // check value only during debugging
413 //checkConstantMemory();
414 }
415
416
417 // Allocate device memory for particle variables,
418 // tied to previously declared pointers in structures
419 void DEM::allocateGlobalDeviceMemory(void)
420 {
421 // Particle memory size
422 unsigned int memSizeF = sizeof(Float) * np;
423 unsigned int memSizeF4 = sizeof(Float4) * np;
424
425 if (verbose == 1)
426 std::cout << " Allocating global device memory: ";
427
428 k.acc = new Float4[np];
429 k.angacc = new Float4[np];
430 #pragma omp parallel for if(np>100)
431 for (unsigned int i = 0; i<np; ++i) {
432 k.acc[i] = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
433 k.angacc[i] = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
434 }
435
436 // Kinematics arrays
437 cudaMalloc((void**)&dev_x, memSizeF4);
438 cudaMalloc((void**)&dev_xyzsum, memSizeF4);
439 cudaMalloc((void**)&dev_vel, memSizeF4);
440 cudaMalloc((void**)&dev_vel0, memSizeF4);
441 cudaMalloc((void**)&dev_acc, memSizeF4);
442 cudaMalloc((void**)&dev_force, memSizeF4);
443 cudaMalloc((void**)&dev_angpos, memSizeF4);
444 cudaMalloc((void**)&dev_angvel, memSizeF4);
445 cudaMalloc((void**)&dev_angvel0, memSizeF4);
446 cudaMalloc((void**)&dev_angacc, memSizeF4);
447 cudaMalloc((void**)&dev_torque, memSizeF4);
448
449 // Particle contact bookkeeping arrays
450 cudaMalloc((void**)&dev_contacts,
451 sizeof(unsigned int)*np*NC);
452 cudaMalloc((void**)&dev_distmod, memSizeF4*NC);
453 cudaMalloc((void**)&dev_delta_t, memSizeF4*NC);
454 cudaMalloc((void**)&dev_bonds, sizeof(uint2)*params.nb0);
455 cudaMalloc((void**)&dev_bonds_delta, sizeof(Float4)*params.nb0);
456 cudaMalloc((void**)&dev_bonds_omega, sizeof(Float4)*params.nb0);
457
458 // Sorted arrays
459 cudaMalloc((void**)&dev_x_sorted, memSizeF4);
460 cudaMalloc((void**)&dev_vel_sorted, memSizeF4);
461 cudaMalloc((void**)&dev_angvel_sorted, memSizeF4);
462
463 // Energy arrays
464 cudaMalloc((void**)&dev_es_dot, memSizeF);
465 cudaMalloc((void**)&dev_ev_dot, memSizeF);
466 cudaMalloc((void**)&dev_es, memSizeF);
467 cudaMalloc((void**)&dev_ev, memSizeF);
468 cudaMalloc((void**)&dev_p, memSizeF);
469
470 // Cell-related arrays
471 cudaMalloc((void**)&dev_gridParticleCellID, sizeof(unsigned int)*np);
472 cudaMalloc((void**)&dev_gridParticleIndex, sizeof(unsigned int)*np);
473 cudaMalloc((void**)&dev_cellStart, sizeof(unsigned int)
474 *grid.num[0]*grid.num[1]*grid.num[2]);
475 cudaMalloc((void**)&dev_cellEnd, sizeof(unsigned int)
476 *grid.num[0]*grid.num[1]*grid.num[2]);
477
478 // Host contact bookkeeping arrays
479 k.contacts = new unsigned int[np*NC];
480 // Initialize contacts lists to np
481 #pragma omp parallel for if(np>100)
482 for (unsigned int i=0; i<(np*NC); ++i)
483 k.contacts[i] = np;
484 k.distmod = new Float4[np*NC];
485 k.delta_t = new Float4[np*NC];
486
487 // Wall arrays
488 cudaMalloc((void**)&dev_walls_wmode, sizeof(int)*walls.nw);
489 cudaMalloc((void**)&dev_walls_nx, sizeof(Float4)*walls.nw);
490 cudaMalloc((void**)&dev_walls_mvfd, sizeof(Float4)*walls.nw);
491 cudaMalloc((void**)&dev_walls_tau_x, sizeof(Float)*walls.nw);
492 cudaMalloc((void**)&dev_walls_tau_eff_x_pp, sizeof(Float)*walls.nw*np);
493 cudaMalloc((void**)&dev_walls_force_pp, sizeof(Float)*walls.nw*np);
494 cudaMalloc((void**)&dev_walls_acc, sizeof(Float)*walls.nw);
495 // dev_walls_force_partial allocated later
496 // dev_walls_tau_eff_x_partial allocated later
497
498 checkForCudaErrors("End of allocateGlobalDeviceMemory");
499 if (verbose == 1)
500 std::cout << "Done" << std::endl;
501 }
502
503 // Allocate global memory on other devices required for "interact" function.
504 // The values of domain_size[ndevices] must be set beforehand.
505 void DEM::allocateHelperDeviceMemory(void)
506 {
507 // Particle memory size
508 unsigned int memSizeF4 = sizeof(Float4) * np;
509
510 // Initialize pointers to per-GPU arrays
511 hdev_gridParticleIndex = (unsigned**)malloc(ndevices*sizeof(unsigned*));
512 hdev_gridCellStart = (unsigned**)malloc(ndevices*sizeof(unsigned*));
513 hdev_gridCellEnd = (unsigned**)malloc(ndevices*sizeof(unsigned*));
514 hdev_x = (Float4**)malloc(ndevices*sizeof(Float4*));
515 hdev_x_sorted = (Float4**)malloc(ndevices*sizeof(Float4*));
516 hdev_vel = (Float4**)malloc(ndevices*sizeof(Float4*));
517 hdev_vel_sorted = (Float4**)malloc(ndevices*sizeof(Float4*));
518 hdev_angvel = (Float4**)malloc(ndevices*sizeof(Float4*));
519 hdev_angvel_sorted = (Float4**)malloc(ndevices*sizeof(Float4*));
520 hdev_walls_nx = (Float4**)malloc(ndevices*sizeof(Float4*));
521 hdev_walls_mvfd = (Float4**)malloc(ndevices*sizeof(Float4*));
522 hdev_distmod = (Float4**)malloc(ndevices*sizeof(Float4*));
523
524 hdev_force = (Float4**)malloc(ndevices*sizeof(Float4*));
525 hdev_torque = (Float4**)malloc(ndevices*sizeof(Float4*));
526 hdev_delta_t = (Float4**)malloc(ndevices*sizeof(Float4*));
527 hdev_es_dot = (Float**)malloc(ndevices*sizeof(Float*));
528 hdev_es = (Float**)malloc(ndevices*sizeof(Float*));
529 hdev_ev_dot = (Float**)malloc(ndevices*sizeof(Float*));
530 hdev_ev = (Float**)malloc(ndevices*sizeof(Float*));
531 hdev_p = (Float**)malloc(ndevices*sizeof(Float*));
532 hdev_walls_force_pp = (Float**)malloc(ndevices*sizeof(Float*));
533 hdev_contacts = (unsigned**)malloc(ndevices*sizeof(unsigned*));
534
535 for (int d=0; d<ndevices; d++) {
536
537 // do not allocate memory on primary GPU
538 if (d == device)
539 continue;
540
541 cudaSetDevice(d);
542
543 // allocate space for full input arrays for interact()
544 cudaMalloc((void**)&hdev_gridParticleIndex[d], sizeof(unsigned int)*np);
545 cudaMalloc((void**)&hdev_gridCellStart[d], sizeof(unsigned int)
546 *grid.num[0]*grid.num[1]*grid.num[2]);
547 cudaMalloc((void**)&hdev_gridCellEnd[d], sizeof(unsigned int)
548 *grid.num[0]*grid.num[1]*grid.num[2]);
549 cudaMalloc((void**)&hdev_x[d], memSizeF4);
550 cudaMalloc((void**)&hdev_x_sorted[d], memSizeF4);
551 cudaMalloc((void**)&hdev_vel[d], memSizeF4);
552 cudaMalloc((void**)&hdev_vel_sorted[d], memSizeF4);
553 cudaMalloc((void**)&hdev_angvel[d], memSizeF4);
554 cudaMalloc((void**)&hdev_angvel_sorted[d], memSizeF4);
555 cudaMalloc((void**)&hdev_walls_nx[d], sizeof(Float4)*walls.nw);
556 cudaMalloc((void**)&hdev_walls_mvfd[d], sizeof(Float4)*walls.nw);
557 cudaMalloc((void**)&hdev_distmod[d], memSizeF4*NC);
558
559 // allocate space for partial output arrays for interact()
560 cudaMalloc((void**)&hdev_force[d], sizeof(Float4)*domain_size[d]);
561 cudaMalloc((void**)&hdev_torque[d], sizeof(Float4)*domain_size[d]);
562 cudaMalloc((void**)&hdev_es_dot[d], sizeof(Float)*domain_size[d]);
563 cudaMalloc((void**)&hdev_ev_dot[d], sizeof(Float)*domain_size[d]);
564 cudaMalloc((void**)&hdev_es[d], sizeof(Float)*domain_size[d]);
565 cudaMalloc((void**)&hdev_ev[d], sizeof(Float)*domain_size[d]);
566 cudaMalloc((void**)&hdev_p[d], sizeof(Float)*domain_size[d]);
567 cudaMalloc((void**)&hdev_walls_force_pp[d],
568 sizeof(Float)*domain_size[d]*walls.nw);
569 cudaMalloc((void**)&hdev_contacts[d],
570 sizeof(unsigned)*domain_size[d]*NC);
571 cudaMalloc((void**)&hdev_delta_t[d], sizeof(Float4)*domain_size[d]*NC);
572
573 checkForCudaErrors("During allocateGlobalDeviceMemoryOtherDevices");
574 }
575 cudaSetDevice(device); // select main GPU
576 }
577
578 void DEM::freeHelperDeviceMemory()
579 {
580 for (int d=0; d<ndevices; d++) {
581
582 // do not allocate memory on primary GPU
583 if (d == device)
584 continue;
585
586 cudaSetDevice(d);
587
588 cudaFree(hdev_gridParticleIndex[d]);
589 cudaFree(hdev_gridCellStart[d]);
590 cudaFree(hdev_gridCellEnd[d]);
591 cudaFree(hdev_x[d]);
592 cudaFree(hdev_vel[d]);
593 cudaFree(hdev_vel_sorted[d]);
594 cudaFree(hdev_angvel[d]);
595 cudaFree(hdev_angvel_sorted[d]);
596 cudaFree(hdev_walls_nx[d]);
597 cudaFree(hdev_walls_mvfd[d]);
598 cudaFree(hdev_distmod[d]);
599
600 cudaFree(hdev_force[d]);
601 cudaFree(hdev_torque[d]);
602 cudaFree(hdev_es_dot[d]);
603 cudaFree(hdev_ev_dot[d]);
604 cudaFree(hdev_es[d]);
605 cudaFree(hdev_ev[d]);
606 cudaFree(hdev_p[d]);
607 cudaFree(hdev_walls_force_pp[d]);
608 cudaFree(hdev_contacts[d]);
609 cudaFree(hdev_delta_t[d]);
610
611 checkForCudaErrors("During helper device cudaFree calls");
612 }
613 cudaSetDevice(device); // select primary GPU
614 }
615
616 void DEM::freeGlobalDeviceMemory()
617 {
618 if (verbose == 1)
619 printf("\nFreeing device memory: ");
620
621 // Particle arrays
622 cudaFree(dev_x);
623 cudaFree(dev_xyzsum);
624 cudaFree(dev_vel);
625 cudaFree(dev_vel0);
626 cudaFree(dev_acc);
627 cudaFree(dev_force);
628 cudaFree(dev_angpos);
629 cudaFree(dev_angvel);
630 cudaFree(dev_angvel0);
631 cudaFree(dev_angacc);
632 cudaFree(dev_torque);
633
634 cudaFree(dev_contacts);
635 cudaFree(dev_distmod);
636 cudaFree(dev_delta_t);
637 cudaFree(dev_bonds);
638 cudaFree(dev_bonds_delta);
639 cudaFree(dev_bonds_omega);
640
641 cudaFree(dev_es_dot);
642 cudaFree(dev_es);
643 cudaFree(dev_ev_dot);
644 cudaFree(dev_ev);
645 cudaFree(dev_p);
646
647 cudaFree(dev_x_sorted);
648 cudaFree(dev_vel_sorted);
649 cudaFree(dev_angvel_sorted);
650
651 // Cell-related arrays
652 cudaFree(dev_gridParticleIndex);
653 cudaFree(dev_cellStart);
654 cudaFree(dev_cellEnd);
655
656 // Wall arrays
657 cudaFree(dev_walls_nx);
658 cudaFree(dev_walls_mvfd);
659 cudaFree(dev_walls_tau_x);
660 cudaFree(dev_walls_force_partial);
661 cudaFree(dev_walls_force_pp);
662 cudaFree(dev_walls_acc);
663 cudaFree(dev_walls_tau_eff_x_pp);
664 cudaFree(dev_walls_tau_eff_x_partial);
665
666 // Fluid arrays
667 if (fluid == 1 && cfd_solver == 0) {
668 freeNSmemDev();
669 }
670 if (fluid == 1 && cfd_solver == 1) {
671 freeDarcyMemDev();
672 }
673
674 //checkForCudaErrors("During cudaFree calls");
675
676 if (verbose == 1)
677 std::cout << "Done" << std::endl;
678 }
679
680
681 void DEM::transferToGlobalDeviceMemory(int statusmsg)
682 {
683 if (verbose == 1 && statusmsg == 1)
684 std::cout << " Transfering data to the device: ";
685
686 // Commonly-used memory sizes
687 unsigned int memSizeF = sizeof(Float) * np;
688 unsigned int memSizeF4 = sizeof(Float4) * np;
689
690 // Copy static-size structure data from host to global device memory
691 //cudaMemcpy(dev_time, &time, sizeof(Time), cudaMemcpyHostToDevice);
692
693 // Kinematic particle values
694 cudaMemcpy( dev_x, k.x,
695 memSizeF4, cudaMemcpyHostToDevice);
696 cudaMemcpy( dev_xyzsum, k.xyzsum,
697 memSizeF4, cudaMemcpyHostToDevice);
698 cudaMemcpy( dev_vel, k.vel,
699 memSizeF4, cudaMemcpyHostToDevice);
700 cudaMemcpy( dev_vel0, k.vel,
701 memSizeF4, cudaMemcpyHostToDevice);
702 cudaMemcpy( dev_acc, k.acc,
703 memSizeF4, cudaMemcpyHostToDevice);
704 cudaMemcpy( dev_force, k.force,
705 memSizeF4, cudaMemcpyHostToDevice);
706 cudaMemcpy( dev_angpos, k.angpos,
707 memSizeF4, cudaMemcpyHostToDevice);
708 cudaMemcpy( dev_angvel, k.angvel,
709 memSizeF4, cudaMemcpyHostToDevice);
710 cudaMemcpy( dev_angvel0, k.angvel,
711 memSizeF4, cudaMemcpyHostToDevice);
712 cudaMemcpy( dev_angacc, k.angacc,
713 memSizeF4, cudaMemcpyHostToDevice);
714 cudaMemcpy( dev_torque, k.torque,
715 memSizeF4, cudaMemcpyHostToDevice);
716 cudaMemcpy( dev_contacts, k.contacts,
717 sizeof(unsigned int)*np*NC, cudaMemcpyHostToDevice);
718 cudaMemcpy( dev_distmod, k.distmod,
719 memSizeF4*NC, cudaMemcpyHostToDevice);
720 cudaMemcpy( dev_delta_t, k.delta_t,
721 memSizeF4*NC, cudaMemcpyHostToDevice);
722 cudaMemcpy( dev_bonds, k.bonds,
723 sizeof(uint2)*params.nb0, cudaMemcpyHostToDevice);
724 cudaMemcpy( dev_bonds_delta, k.bonds_delta,
725 sizeof(Float4)*params.nb0, cudaMemcpyHostToDevice);
726 cudaMemcpy( dev_bonds_omega, k.bonds_omega,
727 sizeof(Float4)*params.nb0, cudaMemcpyHostToDevice);
728
729 // Individual particle energy values
730 cudaMemcpy( dev_es_dot, e.es_dot,
731 memSizeF, cudaMemcpyHostToDevice);
732 cudaMemcpy( dev_es, e.es,
733 memSizeF, cudaMemcpyHostToDevice);
734 cudaMemcpy( dev_ev_dot, e.ev_dot,
735 memSizeF, cudaMemcpyHostToDevice);
736 cudaMemcpy( dev_ev, e.ev,
737 memSizeF, cudaMemcpyHostToDevice);
738 cudaMemcpy( dev_p, e.p,
739 memSizeF, cudaMemcpyHostToDevice);
740
741 // Wall parameters
742 cudaMemcpy( dev_walls_wmode, walls.wmode,
743 sizeof(int)*walls.nw, cudaMemcpyHostToDevice);
744 cudaMemcpy( dev_walls_nx, walls.nx,
745 sizeof(Float4)*walls.nw, cudaMemcpyHostToDevice);
746 cudaMemcpy( dev_walls_mvfd, walls.mvfd,
747 sizeof(Float4)*walls.nw, cudaMemcpyHostToDevice);
748 cudaMemcpy( dev_walls_tau_x, walls.tau_x,
749 sizeof(Float)*walls.nw, cudaMemcpyHostToDevice);
750
751 // Fluid arrays
752 if (fluid == 1) {
753 if (cfd_solver == 0) {
754 transferNStoGlobalDeviceMemory(1);
755 } else if (cfd_solver == 1) {
756 transferDarcyToGlobalDeviceMemory(1);
757 } else {
758 std::cerr << "Error: cfd_solver value not understood ("
759 << cfd_solver << ")" << std::endl;
760 }
761 }
762
763 checkForCudaErrors("End of transferToGlobalDeviceMemory");
764 if (verbose == 1 && statusmsg == 1)
765 std::cout << "Done" << std::endl;
766 }
767
768 void DEM::transferFromGlobalDeviceMemory()
769 {
770 //std::cout << " Transfering data from the device: ";
771
772 // Commonly-used memory sizes
773 unsigned int memSizeF = sizeof(Float) * np;
774 unsigned int memSizeF4 = sizeof(Float4) * np;
775
776 // Copy static-size structure data from host to global device memory
777 //cudaMemcpy(&time, dev_time, sizeof(Time), cudaMemcpyDeviceToHost);
778
779 // Kinematic particle values
780 cudaMemcpy( k.x, dev_x,
781 memSizeF4, cudaMemcpyDeviceToHost);
782 cudaMemcpy( k.xyzsum, dev_xyzsum,
783 memSizeF4, cudaMemcpyDeviceToHost);
784 cudaMemcpy( k.vel, dev_vel,
785 memSizeF4, cudaMemcpyDeviceToHost);
786 cudaMemcpy( k.acc, dev_acc,
787 memSizeF4, cudaMemcpyDeviceToHost);
788 cudaMemcpy( k.force, dev_force,
789 memSizeF4, cudaMemcpyDeviceToHost);
790 cudaMemcpy( k.angpos, dev_angpos,
791 memSizeF4, cudaMemcpyDeviceToHost);
792 cudaMemcpy( k.angvel, dev_angvel,
793 memSizeF4, cudaMemcpyDeviceToHost);
794 cudaMemcpy( k.angacc, dev_angacc,
795 memSizeF4, cudaMemcpyDeviceToHost);
796 cudaMemcpy( k.torque, dev_torque,
797 memSizeF4, cudaMemcpyDeviceToHost);
798 cudaMemcpy( k.contacts, dev_contacts,
799 sizeof(unsigned int)*np*NC, cudaMemcpyDeviceToHost);
800 cudaMemcpy( k.distmod, dev_distmod,
801 memSizeF4*NC, cudaMemcpyDeviceToHost);
802 cudaMemcpy( k.delta_t, dev_delta_t,
803 memSizeF4*NC, cudaMemcpyDeviceToHost);
804 cudaMemcpy( k.bonds, dev_bonds,
805 sizeof(uint2)*params.nb0, cudaMemcpyDeviceToHost);
806 cudaMemcpy( k.bonds_delta, dev_bonds_delta,
807 sizeof(Float4)*params.nb0, cudaMemcpyDeviceToHost);
808 cudaMemcpy( k.bonds_omega, dev_bonds_omega,
809 sizeof(Float4)*params.nb0, cudaMemcpyDeviceToHost);
810
811 // Individual particle energy values
812 cudaMemcpy( e.es_dot, dev_es_dot,
813 memSizeF, cudaMemcpyDeviceToHost);
814 cudaMemcpy( e.es, dev_es,
815 memSizeF, cudaMemcpyDeviceToHost);
816 cudaMemcpy( e.ev_dot, dev_ev_dot,
817 memSizeF, cudaMemcpyDeviceToHost);
818 cudaMemcpy( e.ev, dev_ev,
819 memSizeF, cudaMemcpyDeviceToHost);
820 cudaMemcpy( e.p, dev_p,
821 memSizeF, cudaMemcpyDeviceToHost);
822
823 // Wall parameters
824 cudaMemcpy( walls.wmode, dev_walls_wmode,
825 sizeof(int)*walls.nw, cudaMemcpyDeviceToHost);
826 cudaMemcpy( walls.nx, dev_walls_nx,
827 sizeof(Float4)*walls.nw, cudaMemcpyDeviceToHost);
828 cudaMemcpy( walls.mvfd, dev_walls_mvfd,
829 sizeof(Float4)*walls.nw, cudaMemcpyDeviceToHost);
830 cudaMemcpy( walls.tau_x, dev_walls_tau_x,
831 sizeof(Float)*walls.nw, cudaMemcpyDeviceToHost);
832
833 // Fluid arrays
834 if (fluid == 1 && cfd_solver == 0) {
835 transferNSfromGlobalDeviceMemory(0);
836 }
837 else if (fluid == 1 && cfd_solver == 1) {
838 transferDarcyFromGlobalDeviceMemory(0);
839 checkDarcyStability();
840 }
841
842 //checkForCudaErrors("End of transferFromGlobalDeviceMemory");
843 }
844
845
846 // Iterate through time by explicit time integration
847 void DEM::startTime()
848 {
849 using std::cout;
850 using std::cerr;
851 using std::endl;
852
853 std::string outfile;
854 char file[200];
855 FILE *fp;
856
857 // Synchronization point
858 cudaDeviceSynchronize();
859 checkForCudaErrors("Start of startTime()");
860
861 // Write initial data to output/<sid>.output00000.bin
862 if (time.step_count == 0)
863 writebin(("output/" + sid + ".output00000.bin").c_str());
864
865 // Time variables
866 clock_t tic, toc;
867 double filetimeclock, time_spent;
868 float dev_time_spent;
869
870 // Start CPU clock
871 tic = clock();
872
873 //// GPU workload configuration
874 unsigned int threadsPerBlock = 256;
875 //unsigned int threadsPerBlock = 512;
876
877 // Create enough blocks to accomodate the particles
878 unsigned int blocksPerGrid = iDivUp(np, threadsPerBlock);
879 dim3 dimGrid(blocksPerGrid, 1, 1); // Blocks arranged in 1D grid
880 dim3 dimBlock(threadsPerBlock, 1, 1); // Threads arranged in 1D block
881
882 unsigned int blocksPerGridBonds = iDivUp(params.nb0, threadsPerBlock);
883 dim3 dimGridBonds(blocksPerGridBonds, 1, 1); // Blocks arranged in 1D grid
884
885 // Use 3D block and grid layout for cell-centered fluid calculations
886 dim3 dimBlockFluid(8, 8, 8); // 512 threads per block
887 dim3 dimGridFluid(
888 iDivUp(grid.num[0], dimBlockFluid.x),
889 iDivUp(grid.num[1], dimBlockFluid.y),
890 iDivUp(grid.num[2], dimBlockFluid.z));
891 if (dimGridFluid.z > 64 && fluid == 1) {
892 cerr << "Error: dimGridFluid.z > 64" << endl;
893 exit(1);
894 }
895
896 // Use 3D block and grid layout for cell-face fluid calculations
897 dim3 dimBlockFluidFace(8, 8, 8); // 512 threads per block
898 dim3 dimGridFluidFace(
899 iDivUp(grid.num[0]+1, dimBlockFluidFace.x),
900 iDivUp(grid.num[1]+1, dimBlockFluidFace.y),
901 iDivUp(grid.num[2]+1, dimBlockFluidFace.z));
902 if (dimGridFluidFace.z > 64 && fluid == 1) {
903 cerr << "Error: dimGridFluidFace.z > 64" << endl;
904 exit(1);
905 }
906
907
908 // Shared memory per block
909 unsigned int smemSize = sizeof(unsigned int)*(threadsPerBlock+1);
910
911 // Pre-sum of force per wall
912 cudaMalloc((void**)&dev_walls_force_partial,
913 sizeof(Float)*dimGrid.x*walls.nw);
914
915 // Pre-sum of shear stress per wall
916 cudaMalloc((void**)&dev_walls_tau_eff_x_partial,
917 sizeof(Float)*dimGrid.x*walls.nw);
918
919 // Report to stdout
920 if (verbose == 1) {
921 cout << "\n Device memory allocation and transfer complete.\n"
922 << " - Blocks per grid: "
923 << dimGrid.x << "*" << dimGrid.y << "*" << dimGrid.z << "\n"
924 << " - Threads per block: "
925 << dimBlock.x << "*" << dimBlock.y << "*" << dimBlock.z << "\n"
926 << " - Shared memory required per block: " << smemSize << " bytes"
927 << endl;
928 if (fluid == 1) {
929 cout << " - Blocks per fluid grid: "
930 << dimGridFluid.x << "*" << dimGridFluid.y << "*" <<
931 dimGridFluid.z << "\n"
932 << " - Threads per fluid block: "
933 << dimBlockFluid.x << "*" << dimBlockFluid.y << "*" <<
934 dimBlockFluid.z << endl;
935 }
936 }
937
938 // Initialize counter variable values
939 filetimeclock = 0.0;
940 long iter = 0;
941 const int stdout_report = 10; // no of steps between reporting to stdout
942
943 // Create first status.dat
944 //sprintf(file,"output/%s.status.dat", sid);
945 outfile = "output/" + sid + ".status.dat";
946 fp = fopen(outfile.c_str(), "w");
947 fprintf(fp,"%2.4e %2.4e %d\n",
948 time.current,
949 100.0*time.current/time.total,
950 time.step_count);
951 fclose(fp);
952
953 if (verbose == 1) {
954 cout << "\n Entering the main calculation time loop...\n\n"
955 << " IMPORTANT: Do not close this terminal, doing so will \n"
956 << " terminate this SPHERE process. Follow the \n"
957 << " progress by executing:\n"
958 << " $ ./sphere_status " << sid << endl << endl;
959 }
960
961
962 // Start GPU clock
963 cudaEvent_t dev_tic, dev_toc;
964 cudaEventCreate(&dev_tic);
965 cudaEventCreate(&dev_toc);
966 cudaEventRecord(dev_tic, 0);
967
968 // If profiling is enabled, initialize timers for each kernel
969 cudaEvent_t kernel_tic, kernel_toc;
970 float kernel_elapsed;
971 double t_calcParticleCellID = 0.0;
972 double t_thrustsort = 0.0;
973 double t_reorderArrays = 0.0;
974 double t_topology = 0.0;
975 double t_interact = 0.0;
976 double t_bondsLinear = 0.0;
977 double t_latticeBoltzmannD3Q19 = 0.0;
978 double t_integrate = 0.0;
979 double t_summation = 0.0;
980 double t_integrateWalls = 0.0;
981
982 double t_findPorositiesDev = 0.0;
983 double t_findNSstressTensor = 0.0;
984 double t_findNSdivphiviv = 0.0;
985 double t_findNSdivtau = 0.0;
986 double t_findPredNSvelocities = 0.0;
987 double t_setNSepsilon = 0.0;
988 double t_setNSdirichlet = 0.0;
989 double t_setNSghostNodesDev = 0.0;
990 double t_findNSforcing = 0.0;
991 double t_jacobiIterationNS = 0.0;
992 double t_updateNSvelocityPressure = 0.0;
993
994 double t_findDarcyPorosities = 0.0;
995 double t_setDarcyGhostNodes = 0.0;
996 double t_findDarcyPressureForce = 0.0;
997 double t_setDarcyTopPressure = 0.0;
998 double t_findDarcyPermeabilities = 0.0;
999 double t_findDarcyPermeabilityGradients = 0.0;
1000 //double t_findDarcyPressureChange = 0.0;
1001 double t_updateDarcySolution = 0.0;
1002 double t_copyValues = 0.0;
1003 double t_findDarcyVelocities = 0.0;
1004
1005 if (PROFILING == 1) {
1006 cudaEventCreate(&kernel_tic);
1007 cudaEventCreate(&kernel_toc);
1008 }
1009
1010 // The model start time is saved for profiling performance
1011 double t_start = time.current;
1012 double t_ratio; // ration between time flow in model vs. reality
1013
1014 // Hard-coded parameters for stepwise velocity change (rate-state exp)
1015 int velocity_state = 1; // 1: v1, 2: v2
1016 int change_velocity_state = 0; // 1: increase velocity, 2: decrease vel.
1017 const Float velocity_factor = 10.0; // v2 = v1*velocity_factor
1018 const Float v2_start = 10.0; // seconds
1019 const Float v2_end = 15.0; // seconds
1020
1021 // Index of dynamic top wall (if it exists)
1022 unsigned int wall0_iz = 10000000;
1023 // weight of fluid between two cells in z direction
1024 Float dp_dz;
1025 if (fluid == 1) {
1026 if (cfd_solver == 0)
1027 dp_dz = fabs(ns.rho_f*params.g[2]*grid.L[2]/grid.num[2]);
1028 else if (cfd_solver == 1) {
1029 dp_dz = fabs(darcy.rho_f*params.g[2]*grid.L[2]/grid.num[2]);
1030
1031 // determine pressure at top wall at t=0
1032 darcy.p_top_orig = darcy.p[d_idx(0,0,darcy.nz-1)]
1033 - darcy.p_mod_A
1034 *sin(2.0*M_PI*darcy.p_mod_f*time.current
1035 + darcy.p_mod_phi);
1036 }
1037 }
1038 //std::cout << "dp_dz = " << dp_dz << std::endl;
1039
1040 // Write a log file of the number of iterations it took before
1041 // convergence in the fluid solver
1042 std::ofstream convlog;
1043 if (write_conv_log == 1) {
1044 std::string f = "output/" + sid + "-conv.log";
1045 convlog.open(f.c_str());
1046 }
1047
1048 if (verbose == 1)
1049 cout << " Current simulation time: " << time.current << " s.";
1050
1051 // MAIN CALCULATION TIME LOOP
1052 while (time.current <= time.total) {
1053
1054 // Print current step number to terminal
1055 //printf("\n\n@@@ DEM time step: %ld\n", iter);
1056
1057 // Routine check for errors
1058 checkForCudaErrors("Start of main while loop");
1059
1060 if (np > 0) {
1061
1062 // check if particle positions have finite values
1063 #ifdef CHECK_PARTICLES_FINITE
1064 KERNEL_LAUNCH(checkParticlePositions, dimGrid, dimBlock, 0,
1065 (dev_x));
1066 cudaDeviceSynchronize();
1067 checkForCudaErrorsIter("Post checkParticlePositions", iter);
1068 #endif
1069
1070 // If the grid is adaptive, readjust the grid height to equal the
1071 // positions of the dynamic walls
1072 if (grid.adaptive == 1 && walls.nw > 0) {
1073 updateGridSize();
1074 }
1075
1076 // For each particle:
1077 // Compute hash key (cell index) from position
1078 // in the fine, uniform and homogenous grid.
1079 if (PROFILING == 1)
1080 startTimer(&kernel_tic);
1081 KERNEL_LAUNCH(calcParticleCellID, dimGrid, dimBlock, 0,
1082 (dev_gridParticleCellID,
1083 dev_gridParticleIndex,
1084 dev_x));
1085
1086 // Synchronization point
1087 cudaDeviceSynchronize();
1088 if (PROFILING == 1)
1089 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1090 &t_calcParticleCellID);
1091 checkForCudaErrorsIter("Post calcParticleCellID", iter);
1092
1093
1094 // Sort particle (key, particle ID) pairs by hash key with Thrust
1095 // radix sort
1096 if (PROFILING == 1)
1097 startTimer(&kernel_tic);
1098 #ifdef SPHERE_GPU
1099 thrust::sort_by_key(
1100 thrust::device_ptr<uint>(dev_gridParticleCellID),
1101 thrust::device_ptr<uint>(dev_gridParticleCellID + np),
1102 thrust::device_ptr<uint>(dev_gridParticleIndex));
1103 #else
1104 { // stable sort (key, index) pairs to match thrust radix order
1105 std::vector<std::pair<unsigned int, unsigned int> > kv(np);
1106 for (unsigned int i = 0; i < np; ++i)
1107 kv[i] = std::make_pair(dev_gridParticleCellID[i],
1108 dev_gridParticleIndex[i]);
1109 std::stable_sort(kv.begin(), kv.end(),
1110 [](const std::pair<unsigned int, unsigned int>& a,
1111 const std::pair<unsigned int, unsigned int>& b)
1112 { return a.first < b.first; });
1113 for (unsigned int i = 0; i < np; ++i) {
1114 dev_gridParticleCellID[i] = kv[i].first;
1115 dev_gridParticleIndex[i] = kv[i].second;
1116 }
1117 }
1118 #endif
1119 cudaDeviceSynchronize(); // Maybe Thrust synchronizes implicitly?
1120 if (PROFILING == 1)
1121 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1122 &t_thrustsort);
1123 checkForCudaErrorsIter("Post thrust::sort_by_key", iter);
1124
1125
1126 // Zero cell array values by setting cellStart to its highest
1127 // possible value, specified with pointer value 0xffffffff, which
1128 // for a 32 bit unsigned int is 4294967295.
1129 cudaMemset(dev_cellStart, 0xffffffff,
1130 grid.num[0]*grid.num[1]*grid.num[2]*sizeof(unsigned int));
1131 cudaDeviceSynchronize();
1132 checkForCudaErrorsIter("Post cudaMemset", iter);
1133
1134 // Use sorted order to reorder particle arrays (position,
1135 // velocities, radii) to ensure coherent memory access. Save ordered
1136 // configurations in new arrays (*_sorted).
1137 if (PROFILING == 1)
1138 startTimer(&kernel_tic);
1139 KERNEL_LAUNCH(reorderArrays, dimGrid, dimBlock, smemSize,
1140 (dev_cellStart,
1141 dev_cellEnd,
1142 dev_gridParticleCellID,
1143 dev_gridParticleIndex,
1144 dev_x, dev_vel,
1145 dev_angvel,
1146 dev_x_sorted,
1147 dev_vel_sorted,
1148 dev_angvel_sorted));
1149
1150 // Synchronization point
1151 cudaDeviceSynchronize();
1152 if (PROFILING == 1)
1153 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1154 &t_reorderArrays);
1155 checkForCudaErrorsIter("Post reorderArrays", iter);
1156
1157 // The contact search in topology() is only necessary for
1158 // determining the accumulated shear distance needed in the linear
1159 // elastic and nonlinear contact force model
1160 if (params.contactmodel == 2 || params.contactmodel == 3) {
1161 // For each particle: Search contacts in neighbor cells
1162 if (PROFILING == 1)
1163 startTimer(&kernel_tic);
1164 KERNEL_LAUNCH(topology, dimGrid, dimBlock, 0,
1165 (dev_cellStart,
1166 dev_cellEnd,
1167 dev_gridParticleIndex,
1168 dev_x_sorted,
1169 dev_contacts,
1170 dev_distmod));
1171
1172 // Synchronization point
1173 cudaDeviceSynchronize();
1174 if (PROFILING == 1)
1175 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1176 &t_topology);
1177 checkForCudaErrorsIter(
1178 "Post topology: One or more particles moved "
1179 "outside the grid.\nThis could possibly be caused by a "
1180 "numerical instability.\nIs the computational time step"
1181 " too large?", iter);
1182 }
1183
1184 // For each particle process collisions and compute resulting forces
1185 //cudaPrintfInit();
1186 if (PROFILING == 1)
1187 startTimer(&kernel_tic);
1188 KERNEL_LAUNCH(interact, dimGrid, dimBlock, 0,
1189 (dev_gridParticleIndex,
1190 dev_cellStart,
1191 dev_cellEnd,
1192 dev_x,
1193 dev_x_sorted,
1194 dev_vel_sorted,
1195 dev_angvel_sorted,
1196 dev_vel,
1197 dev_angvel,
1198 dev_force,
1199 dev_torque,
1200 dev_es_dot,
1201 dev_ev_dot,
1202 dev_es,
1203 dev_ev,
1204 dev_p,
1205 dev_walls_nx,
1206 dev_walls_mvfd,
1207 dev_walls_force_pp,
1208 dev_contacts,
1209 dev_distmod,
1210 dev_delta_t));
1211
1212 // Synchronization point
1213 cudaDeviceSynchronize();
1214 //cudaPrintfDisplay(stdout, true);
1215 if (PROFILING == 1)
1216 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1217 &t_interact);
1218 checkForCudaErrorsIter(
1219 "Post interact - often caused if particles move "
1220 "outside the grid", iter);
1221
1222 // Process particle pairs
1223 if (params.nb0 > 0) {
1224 if (PROFILING == 1)
1225 startTimer(&kernel_tic);
1226 KERNEL_LAUNCH(bondsLinear, dimGridBonds, dimBlock, 0,
1227 (dev_bonds,
1228 dev_bonds_delta,
1229 dev_bonds_omega,
1230 dev_x,
1231 dev_vel,
1232 dev_angvel,
1233 dev_force,
1234 dev_torque));
1235 // Synchronization point
1236 cudaDeviceSynchronize();
1237 //cudaPrintfDisplay(stdout, true);
1238 if (PROFILING == 1)
1239 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1240 &t_bondsLinear);
1241 checkForCudaErrorsIter("Post bondsLinear", iter);
1242 }
1243 }
1244
1245 // Solve fluid flow through the grid
1246 if (fluid == 1) {
1247
1248 // Navier-Stokes solution
1249 if (cfd_solver == 0) {
1250
1251 checkForCudaErrorsIter("Before findPorositiesDev", iter);
1252 // Find cell porosities, average particle velocities, and
1253 // average particle diameters. These are needed for predicting
1254 // the fluid velocities
1255 if (PROFILING == 1)
1256 startTimer(&kernel_tic);
1257 //findPorositiesVelocitiesDiametersSphericalGradient
1258 KERNEL_LAUNCH(findPorositiesVelocitiesDiametersSpherical, dimGridFluid, dimBlockFluid, 0,
1259 (dev_cellStart,
1260 dev_cellEnd,
1261 dev_x_sorted,
1262 dev_vel_sorted,
1263 dev_ns_phi,
1264 dev_ns_dphi,
1265 dev_ns_vp_avg,
1266 dev_ns_d_avg,
1267 iter,
1268 np,
1269 ns.c_phi));
1270 cudaDeviceSynchronize();
1271 if (PROFILING == 1)
1272 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1273 &t_findPorositiesDev);
1274 checkForCudaErrorsIter("Post findPorositiesDev", iter);
1275
1276 #ifdef CFD_DEM_COUPLING
1277 /*if (params.nu <= 0.0) {
1278 std::cerr << "Error! The fluid needs a positive viscosity "
1279 "value in order to simulate particle-fluid interaction."
1280 << std::endl;
1281 exit(1);
1282 }*/
1283 if (iter == 0) {
1284 // set cell center ghost nodes
1285 KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
1286 (dev_ns_v, ns.bc_bot, ns.bc_top));
1287
1288 // find cell face velocities
1289 KERNEL_LAUNCH(interpolateCenterToFace, dimGridFluidFace, dimBlockFluidFace, 0,
1290 (dev_ns_v,
1291 dev_ns_v_x,
1292 dev_ns_v_y,
1293 dev_ns_v_z));
1294 cudaDeviceSynchronize();
1295 checkForCudaErrors("Post interpolateCenterToFace");
1296 }
1297
1298 KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
1299 (dev_ns_v_x,
1300 dev_ns_v_y,
1301 dev_ns_v_z,
1302 ns.bc_bot,
1303 ns.bc_top));
1304 cudaDeviceSynchronize();
1305 checkForCudaErrorsIter("Post setNSghostNodesFace", iter);
1306
1307 KERNEL_LAUNCH(findFaceDivTau, dimGridFluidFace, dimBlockFluidFace, 0,
1308 (dev_ns_v_x,
1309 dev_ns_v_y,
1310 dev_ns_v_z,
1311 ns.mu,
1312 dev_ns_div_tau_x,
1313 dev_ns_div_tau_y,
1314 dev_ns_div_tau_z));
1315 cudaDeviceSynchronize();
1316 checkForCudaErrorsIter("Post findFaceDivTau", iter);
1317
1318 KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluid, 0,
1319 (dev_ns_div_tau_x,
1320 dev_ns_div_tau_y,
1321 dev_ns_div_tau_z,
1322 ns.bc_bot,
1323 ns.bc_top));
1324 cudaDeviceSynchronize();
1325 checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_div_tau)",
1326 iter);
1327
1328 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1329 (dev_ns_p, ns.bc_bot, ns.bc_top));
1330 cudaDeviceSynchronize();
1331 checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)", iter);
1332
1333 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1334 (dev_ns_phi, ns.bc_bot, ns.bc_top));
1335 cudaDeviceSynchronize();
1336 checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)", iter);
1337
1338
1339 if (np > 0) {
1340
1341 // Per particle, find the fluid-particle interaction force
1342 // f_pf and apply it to the particle
1343 KERNEL_LAUNCH(findInteractionForce, dimGrid, dimBlock, 0,
1344 (dev_x,
1345 dev_vel,
1346 dev_ns_phi,
1347 dev_ns_p,
1348 dev_ns_v_x,
1349 dev_ns_v_y,
1350 dev_ns_v_z,
1351 dev_ns_div_tau_x,
1352 dev_ns_div_tau_y,
1353 dev_ns_div_tau_z,
1354 //ns.c_v,
1355 ns.mu,
1356 ns.rho_f,
1357 dev_ns_f_pf,
1358 dev_force,
1359 dev_ns_f_d,
1360 dev_ns_f_p,
1361 dev_ns_f_v,
1362 dev_ns_f_sum));
1363 cudaDeviceSynchronize();
1364 checkForCudaErrorsIter("Post findInteractionForce", iter);
1365
1366 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1367 (dev_ns_p, ns.bc_bot, ns.bc_top));
1368 cudaDeviceSynchronize();
1369 checkForCudaErrorsIter("Post setNSghostNodes(dev_ns_p)",
1370 iter);
1371
1372 // Apply fluid-particle interaction force to the fluid
1373 KERNEL_LAUNCH(applyInteractionForceToFluid, dimGridFluid, dimBlockFluid, 0,
1374 (dev_gridParticleIndex,
1375 dev_cellStart,
1376 dev_cellEnd,
1377 dev_ns_f_pf,
1378 dev_ns_F_pf));
1379 //dev_ns_F_pf_x,
1380 //dev_ns_F_pf_y,
1381 //dev_ns_F_pf_z);
1382 cudaDeviceSynchronize();
1383 checkForCudaErrorsIter("Post applyInteractionForceToFluid",
1384 iter);
1385
1386 KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
1387 (dev_ns_F_pf, ns.bc_bot, ns.bc_top));
1388 cudaDeviceSynchronize();
1389 checkForCudaErrorsIter("Post setNSghostNodes(F_pf)", iter);
1390 }
1391 #endif
1392
1393 if ((iter % ns.ndem) == 0) {
1394 // Initial guess for the top epsilon values. These may be
1395 // changed in setUpperPressureNS
1396 // TODO: Check if this should only be set when top bc=Dirichlet
1397 Float pressure = ns.p[idx(0,0,ns.nz-1)];
1398 Float pressure_new = pressure; // Dirichlet
1399 Float epsilon_value = pressure_new - ns.beta*pressure;
1400 KERNEL_LAUNCH(setNSepsilonTop, dimGridFluid, dimBlockFluid, 0,
1401 (dev_ns_epsilon,
1402 dev_ns_epsilon_new,
1403 epsilon_value));
1404 cudaDeviceSynchronize();
1405 checkForCudaErrorsIter("Post setNSepsilonTop", iter);
1406
1407 #if defined(REPORT_EPSILON) || defined(REPORT_V_P_COMPONENTS) || defined(REPORT_V_C_COMPONENTS)
1408 std::cout
1409 << "\n\n@@@@@@ TIME STEP " << iter << " @@@"
1410 << std::endl;
1411 #endif
1412
1413 // find cell containing top wall
1414 if (walls.nw > 0 &&
1415 (walls.wmode[0] == 1 || walls.wmode[0] == 3)) {
1416 wall0_iz = walls.nx->w/(grid.L[2]/grid.num[2]);
1417 KERNEL_LAUNCH(setNSepsilonAtTopWall, dimGridFluid, dimBlockFluid, 0,
1418 (dev_ns_epsilon,
1419 dev_ns_epsilon_new,
1420 wall0_iz,
1421 epsilon_value,
1422 dp_dz));
1423 cudaDeviceSynchronize();
1424 checkForCudaErrorsIter("Post setNSepsilonAtTopWall",
1425 iter);
1426
1427 #ifdef REPORT_EPSILON
1428 std::cout
1429 << "\n###### EPSILON setNSepsilonAtTopWall "
1430 << "######" << std::endl;
1431 transferNSepsilonFromGlobalDeviceMemory();
1432 printNSarray(stdout, ns.epsilon, "epsilon");
1433 #endif
1434 }
1435
1436 // Modulate the pressures at the upper boundary cells
1437 if ((ns.p_mod_A > 1.0e-5 || ns.p_mod_A < -1.0e-5) &&
1438 ns.p_mod_f > 1.0e-7) {
1439 // original pressure
1440 Float new_pressure = ns.p[idx(0,0,ns.nz-1)]
1441 + ns.p_mod_A*sin(2.0*M_PI*ns.p_mod_f*time.current
1442 + ns.p_mod_phi);
1443 KERNEL_LAUNCH(setUpperPressureNS, dimGridFluid, dimBlockFluid, 0,
1444 (dev_ns_p,
1445 dev_ns_epsilon,
1446 dev_ns_epsilon_new,
1447 ns.beta,
1448 new_pressure));
1449 cudaDeviceSynchronize();
1450 checkForCudaErrorsIter("Post setUpperPressureNS", iter);
1451
1452 #ifdef REPORT_MORE_EPSILON
1453 std::cout
1454 << "\n@@@@@@ TIME STEP " << iter << " @@@@@@"
1455 << "\n###### EPSILON AFTER setUpperPressureNS "
1456 << "######" << std::endl;
1457 transferNSepsilonFromGlobalDeviceMemory();
1458 printNSarray(stdout, ns.epsilon, "epsilon");
1459 #endif
1460 }
1461
1462 // Set the values of the ghost nodes in the grid
1463 if (PROFILING == 1)
1464 startTimer(&kernel_tic);
1465
1466 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1467 (dev_ns_p, ns.bc_bot, ns.bc_top));
1468
1469 //setNSghostNodes<Float3><<<dimGridFluid, dimBlockFluid>>>(
1470 //dev_ns_v, ns.bc_bot, ns.bc_top);
1471
1472 KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
1473 (dev_ns_v_p_x,
1474 dev_ns_v_p_y,
1475 dev_ns_v_p_z,
1476 ns.bc_bot, ns.bc_top));
1477
1478 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1479 (dev_ns_phi, ns.bc_bot, ns.bc_top));
1480
1481 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1482 (dev_ns_dphi, ns.bc_bot, ns.bc_top));
1483
1484 cudaDeviceSynchronize();
1485 if (PROFILING == 1)
1486 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1487 &t_setNSghostNodesDev);
1488 checkForCudaErrorsIter("Post setNSghostNodesDev", iter);
1489 /*std::cout
1490 << "\n###### EPSILON AFTER setNSghostNodesDev #####"
1491 << std::endl;
1492 transferNSepsilonFromGlobalDeviceMemory();
1493 printNSarray(stdout, ns.epsilon, "epsilon");*/
1494
1495 // interpolate velocities to cell centers which makes
1496 // velocity prediction easier
1497 KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
1498 (dev_ns_v_x,
1499 dev_ns_v_y,
1500 dev_ns_v_z,
1501 dev_ns_v));
1502 cudaDeviceSynchronize();
1503 checkForCudaErrorsIter(
1504 "Post interpolateFaceToCenter", iter);
1505
1506 // Set cell-center velocity ghost nodes
1507 KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
1508 (dev_ns_v, ns.bc_bot, ns.bc_top));
1509 cudaDeviceSynchronize();
1510 checkForCudaErrorsIter("Post setNSghostNodes(v)", iter);
1511
1512 // Find the divergence of phi*vi*v, needed for predicting
1513 // the fluid velocities
1514 if (PROFILING == 1)
1515 startTimer(&kernel_tic);
1516 KERNEL_LAUNCH(findNSdivphiviv, dimGridFluid, dimBlockFluid, 0,
1517 (dev_ns_phi,
1518 dev_ns_v,
1519 dev_ns_div_phi_vi_v));
1520 cudaDeviceSynchronize();
1521 if (PROFILING == 1)
1522 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1523 &t_findNSdivphiviv);
1524 checkForCudaErrorsIter("Post findNSdivphiviv", iter);
1525
1526 // Set cell-center ghost nodes
1527 KERNEL_LAUNCH(setNSghostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
1528 (dev_ns_div_phi_vi_v, ns.bc_bot, ns.bc_top));
1529 cudaDeviceSynchronize();
1530 checkForCudaErrorsIter("Post setNSghostNodes(div_phi_vi_v)",
1531 iter);
1532
1533 // Predict the fluid velocities on the base of the old
1534 // pressure field and ignoring the incompressibility
1535 // constraint
1536 if (PROFILING == 1)
1537 startTimer(&kernel_tic);
1538 KERNEL_LAUNCH(findPredNSvelocities, dimGridFluidFace, dimBlockFluidFace, 0,
1539 (dev_ns_p,
1540 dev_ns_v_x,
1541 dev_ns_v_y,
1542 dev_ns_v_z,
1543 dev_ns_phi,
1544 dev_ns_dphi,
1545 dev_ns_div_tau_x,
1546 dev_ns_div_tau_y,
1547 dev_ns_div_tau_z,
1548 dev_ns_div_phi_vi_v,
1549 ns.bc_bot,
1550 ns.bc_top,
1551 ns.beta,
1552 dev_ns_F_pf,
1553 ns.ndem,
1554 wall0_iz,
1555 ns.c_v,
1556 ns.mu,
1557 ns.rho_f,
1558 dev_ns_v_p_x,
1559 dev_ns_v_p_y,
1560 dev_ns_v_p_z));
1561 cudaDeviceSynchronize();
1562 if (PROFILING == 1)
1563 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1564 &t_findPredNSvelocities);
1565 checkForCudaErrorsIter("Post findPredNSvelocities", iter);
1566
1567 KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
1568 (dev_ns_v_p_x,
1569 dev_ns_v_p_y,
1570 dev_ns_v_p_z,
1571 ns.bc_bot, ns.bc_top));
1572 cudaDeviceSynchronize();
1573 checkForCudaErrorsIter(
1574 "Post setNSghostNodesFace(dev_ns_v_p)", iter);
1575
1576 KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
1577 (dev_ns_v_p_x,
1578 dev_ns_v_p_y,
1579 dev_ns_v_p_z,
1580 dev_ns_v_p));
1581 cudaDeviceSynchronize();
1582 checkForCudaErrorsIter(
1583 "Post interpolateFaceToCenter", iter);
1584
1585
1586 // In the first iteration of the sphere program, we'll need
1587 // to manually estimate the values of epsilon. In the
1588 // subsequent iterations, the previous values are used.
1589 if (iter == 0) {
1590
1591 // Define the first estimate of the values of epsilon.
1592 // The initial guess depends on the value of ns.beta.
1593 Float pressure = ns.p[idx(2,2,2)];
1594 Float pressure_new = pressure; // Guess p_curr = p_new
1595 Float epsilon_value = pressure_new - ns.beta*pressure;
1596 if (PROFILING == 1)
1597 startTimer(&kernel_tic);
1598 KERNEL_LAUNCH(setNSepsilonInterior, dimGridFluid, dimBlockFluid, 0,
1599 (dev_ns_epsilon, epsilon_value));
1600 cudaDeviceSynchronize();
1601
1602 KERNEL_LAUNCH(setNSnormZero, dimGridFluid, dimBlockFluid, 0,
1603 (dev_ns_norm));
1604 cudaDeviceSynchronize();
1605
1606 if (PROFILING == 1)
1607 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1608 &t_setNSepsilon);
1609 checkForCudaErrorsIter("Post setNSepsilonInterior",
1610 iter);
1611
1612 #ifdef REPORT_MORE_EPSILON
1613 std::cout
1614 << "\n###### EPSILON AFTER setNSepsilonInterior "
1615 << "######" << std::endl;
1616 transferNSepsilonFromGlobalDeviceMemory();
1617 printNSarray(stdout, ns.epsilon, "epsilon");
1618 #endif
1619
1620 // Set the epsilon values at the lower boundary
1621 pressure = ns.p[idx(0,0,0)];
1622 pressure_new = pressure; // Guess p_current = p_new
1623 epsilon_value = pressure_new - ns.beta*pressure;
1624 if (PROFILING == 1)
1625 startTimer(&kernel_tic);
1626 KERNEL_LAUNCH(setNSepsilonBottom, dimGridFluid, dimBlockFluid, 0,
1627 (dev_ns_epsilon,
1628 dev_ns_epsilon_new,
1629 epsilon_value));
1630 cudaDeviceSynchronize();
1631 if (PROFILING == 1)
1632 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1633 &t_setNSdirichlet);
1634 checkForCudaErrorsIter("Post setNSepsilonBottom", iter);
1635
1636 #ifdef REPORT_MORE_EPSILON
1637 std::cout
1638 << "\n###### EPSILON AFTER setNSepsilonBottom "
1639 << "######" << std::endl;
1640 transferNSepsilonFromGlobalDeviceMemory();
1641 printNSarray(stdout, ns.epsilon, "epsilon");
1642 #endif
1643
1644 /*setNSghostNodes<Float>
1645 <<<dimGridFluid, dimBlockFluid>>>(
1646 dev_ns_epsilon);
1647 cudaDeviceSynchronize();
1648 checkForCudaErrors(
1649 "Post setNSghostNodesFloat(dev_ns_epsilon)",
1650 iter);*/
1651 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1652 (dev_ns_epsilon,
1653 ns.bc_bot, ns.bc_top));
1654 cudaDeviceSynchronize();
1655 checkForCudaErrorsIter("Post setNSghostNodesEpsilon(1)",
1656 iter);
1657
1658 #ifdef REPORT_MORE_EPSILON
1659 std::cout <<
1660 "\n###### EPSILON AFTER setNSghostNodes(epsilon) "
1661 << "######" << std::endl;
1662 transferNSepsilonFromGlobalDeviceMemory();
1663 printNSarray(stdout, ns.epsilon, "epsilon");
1664 #endif
1665 }
1666
1667 // Solve the system of epsilon using a Jacobi iterative
1668 // solver. The average normalized residual is initialized
1669 // to a large value.
1670 //double avg_norm_res;
1671 double max_norm_res;
1672
1673 // Write a log file of the normalized residuals during the
1674 // Jacobi iterations
1675 std::ofstream reslog;
1676 if (write_res_log == 1)
1677 reslog.open("max_res_norm.dat");
1678
1679 // transfer normalized residuals from GPU to CPU
1680 #ifdef REPORT_MORE_EPSILON
1681 std::cout << "\n###### BEFORE FIRST JACOBI ITERATION ######"
1682 << "\n@@@@@@ TIME STEP " << iter << " @@@@@@"
1683 << std::endl;
1684 transferNSepsilonFromGlobalDeviceMemory();
1685 printNSarray(stdout, ns.epsilon, "epsilon");
1686 #endif
1687
1688 for (unsigned int nijac = 0; nijac<ns.maxiter; ++nijac) {
1689
1690 // Only grad(epsilon) changes during the Jacobi
1691 // iterations. The remaining terms of the forcing
1692 // function are only calculated during the first
1693 // iteration.
1694 if (PROFILING == 1)
1695 startTimer(&kernel_tic);
1696 KERNEL_LAUNCH(findNSforcing, dimGridFluid, dimBlockFluid, 0,
1697 (dev_ns_epsilon,
1698 dev_ns_phi,
1699 dev_ns_dphi,
1700 dev_ns_v_p,
1701 dev_ns_v_p_x,
1702 dev_ns_v_p_y,
1703 dev_ns_v_p_z,
1704 nijac,
1705 ns.ndem,
1706 ns.c_v,
1707 ns.rho_f,
1708 dev_ns_f1,
1709 dev_ns_f2,
1710 dev_ns_f));
1711 cudaDeviceSynchronize();
1712 if (PROFILING == 1)
1713 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1714 &t_findNSforcing);
1715 checkForCudaErrorsIter("Post findNSforcing", iter);
1716 /*setNSghostNodesForcing
1717 <<dimGridFluid, dimBlockFluid>>>(
1718 dev_ns_f1,
1719 dev_ns_f2,
1720 dev_ns_f,
1721 nijac);
1722 cudaDeviceSynchronize();
1723 checkForCudaErrors("Post setNSghostNodesForcing",
1724 iter);*/
1725
1726 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1727 (dev_ns_epsilon,
1728 ns.bc_bot, ns.bc_top));
1729 cudaDeviceSynchronize();
1730 checkForCudaErrorsIter("Post setNSghostNodesEpsilon(2)",
1731 iter);
1732
1733 #ifdef REPORT_EPSILON
1734 std::cout << "\n###### JACOBI ITERATION "
1735 << nijac
1736 << " after setNSghostNodes(epsilon,2) ######"
1737 << std::endl;
1738 transferNSepsilonFromGlobalDeviceMemory();
1739 printNSarray(stdout, ns.epsilon, "epsilon");
1740 #endif
1741
1742 // Perform a single Jacobi iteration
1743 if (PROFILING == 1)
1744 startTimer(&kernel_tic);
1745 KERNEL_LAUNCH(jacobiIterationNS, dimGridFluid, dimBlockFluid, 0,
1746 (dev_ns_epsilon,
1747 dev_ns_epsilon_new,
1748 dev_ns_norm,
1749 dev_ns_f,
1750 ns.bc_bot,
1751 ns.bc_top,
1752 ns.theta,
1753 wall0_iz,
1754 dp_dz));
1755 cudaDeviceSynchronize();
1756 if (PROFILING == 1)
1757 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1758 &t_jacobiIterationNS);
1759 checkForCudaErrorsIter("Post jacobiIterationNS", iter);
1760
1761 // set Dirichlet and Neumann BC at cells containing top
1762 // wall
1763 /*if (walls.nw > 0 && walls.wmode[0] == 1) {
1764 setNSepsilonAtTopWall
1765 <<<dimGridFluid, dimBlockFluid>>>(
1766 dev_ns_epsilon,
1767 dev_ns_epsilon_new,
1768 wall0_iz,
1769 epsilon_value,
1770 dp_dz);
1771 cudaDeviceSynchronize();
1772 checkForCudaErrorsIter("Post setNSepsilonAtTopWall",
1773 iter);
1774 }*/
1775
1776 // Copy new values to current values
1777 KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
1778 (dev_ns_epsilon_new,
1779 dev_ns_epsilon));
1780 cudaDeviceSynchronize();
1781 checkForCudaErrorsIter
1782 ("Post copyValues (epsilon_new->epsilon)", iter);
1783
1784 #ifdef REPORT_EPSILON
1785 std::cout << "\n###### JACOBI ITERATION "
1786 << nijac << " after jacobiIterationNS ######"
1787 << std::endl;
1788 transferNSepsilonFromGlobalDeviceMemory();
1789 printNSarray(stdout, ns.epsilon, "epsilon");
1790 #endif
1791
1792 if (nijac % nijacnorm == 0) {
1793
1794 // Read the normalized residuals from the device
1795 transferNSnormFromGlobalDeviceMemory();
1796
1797 // Write the normalized residuals to the terminal
1798 //printNSarray(stdout, ns.norm, "norm");
1799
1800 // Find the maximum value of the normalized
1801 // residuals
1802 max_norm_res = maxNormResNS();
1803
1804 // Write the Jacobi iteration number and maximum
1805 // value of the normalized residual to the log file
1806 if (write_res_log == 1)
1807 reslog << nijac << '\t' << max_norm_res
1808 << std::endl;
1809 }
1810
1811 if (max_norm_res < ns.tolerance) {
1812
1813 if (write_conv_log == 1
1814 && iter % conv_log_interval == 0)
1815 convlog << iter+1 << '\t' << nijac << std::endl;
1816
1817 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1818 (dev_ns_epsilon,
1819 ns.bc_bot, ns.bc_top));
1820 cudaDeviceSynchronize();
1821 checkForCudaErrorsIter
1822 ("Post setNSghostNodesEpsilon(4)", iter);
1823
1824 // Apply smoothing if requested
1825 if (ns.gamma > 0.0) {
1826
1827 KERNEL_LAUNCH(smoothing, dimGridFluid, dimBlockFluid, 0,
1828 (dev_ns_epsilon,
1829 ns.gamma,
1830 ns.bc_bot, ns.bc_top));
1831 cudaDeviceSynchronize();
1832 checkForCudaErrorsIter("Post smoothing", iter);
1833
1834 KERNEL_LAUNCH(setNSghostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1835 (dev_ns_epsilon,
1836 ns.bc_bot, ns.bc_top));
1837 cudaDeviceSynchronize();
1838 checkForCudaErrorsIter
1839 ("Post setNSghostNodesEpsilon(4)", iter);
1840 }
1841
1842 #ifdef REPORT_EPSILON
1843 std::cout << "\n###### JACOBI ITERATION "
1844 << nijac << " after smoothing ######"
1845 << std::endl;
1846 transferNSepsilonFromGlobalDeviceMemory();
1847 printNSarray(stdout, ns.epsilon, "epsilon");
1848 #endif
1849
1850 break; // solution has converged, exit Jacobi loop
1851 }
1852
1853 if (nijac >= ns.maxiter-1) {
1854
1855 if (write_conv_log == 1)
1856 convlog << iter+1 << '\t' << nijac << std::endl;
1857
1858 std::cerr << "\nIteration " << iter << ", time "
1859 << iter*time.dt << " s: "
1860 "Error, the epsilon solution in the fluid "
1861 "calculations did not converge. Try increasing "
1862 "the value of 'ns.maxiter' (" << ns.maxiter
1863 << ") or increase 'ns.tolerance' ("
1864 << ns.tolerance << ")." << std::endl;
1865 }
1866 //break; // end after Jacobi first iteration
1867 } // end Jacobi iteration loop
1868
1869 if (write_res_log == 1)
1870 reslog.close();
1871
1872 // Find the new pressures and velocities
1873 if (PROFILING == 1)
1874 startTimer(&kernel_tic);
1875
1876 KERNEL_LAUNCH(updateNSpressure, dimGridFluid, dimBlockFluid, 0,
1877 (dev_ns_epsilon,
1878 ns.beta,
1879 dev_ns_p));
1880 cudaDeviceSynchronize();
1881 checkForCudaErrorsIter("Post updateNSpressure", iter);
1882
1883 KERNEL_LAUNCH(updateNSvelocity, dimGridFluidFace, dimBlockFluidFace, 0,
1884 (dev_ns_v_p_x,
1885 dev_ns_v_p_y,
1886 dev_ns_v_p_z,
1887 dev_ns_phi,
1888 dev_ns_epsilon,
1889 ns.beta,
1890 ns.bc_bot,
1891 ns.bc_top,
1892 ns.ndem,
1893 ns.c_v,
1894 ns.rho_f,
1895 wall0_iz,
1896 iter,
1897 dev_ns_v_x,
1898 dev_ns_v_y,
1899 dev_ns_v_z));
1900 cudaDeviceSynchronize();
1901 if (PROFILING == 1)
1902 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1903 &t_updateNSvelocityPressure);
1904 checkForCudaErrorsIter("Post updateNSvelocity", iter);
1905
1906 KERNEL_LAUNCH(setNSghostNodesFace<Float>, dimGridFluidFace, dimBlockFluidFace, 0,
1907 (dev_ns_v_p_x,
1908 dev_ns_v_p_y,
1909 dev_ns_v_p_z,
1910 ns.bc_bot, ns.bc_top));
1911 cudaDeviceSynchronize();
1912 checkForCudaErrorsIter(
1913 "Post setNSghostNodesFace(dev_ns_v)", iter);
1914
1915 KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
1916 (dev_ns_v_x,
1917 dev_ns_v_y,
1918 dev_ns_v_z,
1919 dev_ns_v));
1920 cudaDeviceSynchronize();
1921 checkForCudaErrorsIter("Post interpolateFaceToCenter",
1922 iter);
1923 } // end iter % ns.dem == 0
1924 } // end cfd_solver == 0
1925
1926 // Darcy solution
1927 else if (cfd_solver == 1) {
1928
1929 #if defined(REPORT_EPSILON) || defined(REPORT_FORCING_TERMS)
1930 std::cout << "\n\n@@@@@@ TIME STEP " << iter << " @@@"
1931 << std::endl;
1932 #endif
1933
1934 if (walls.nw > 0 &&
1935 (walls.wmode[0] == 1 || walls.wmode[0] == 3)) {
1936 wall0_iz = walls.nx->w/(grid.L[2]/grid.num[2]);
1937 }
1938
1939 if (np > 0) {
1940
1941 if (PROFILING == 1)
1942 startTimer(&kernel_tic);
1943 KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
1944 (dev_darcy_p,
1945 darcy.bc_xn, darcy.bc_xp,
1946 darcy.bc_yn, darcy.bc_yp,
1947 darcy.bc_bot, darcy.bc_top));
1948 cudaDeviceSynchronize();
1949 if (PROFILING == 1)
1950 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1951 &t_setDarcyGhostNodes);
1952 checkForCudaErrorsIter("Post setDarcyGhostNodes("
1953 "dev_darcy_p) before findDarcyPressureForce", iter);
1954
1955 if (PROFILING == 1)
1956 startTimer(&kernel_tic);
1957 KERNEL_LAUNCH(findDarcyPressureGradient, dimGridFluid, dimBlockFluid, 0,
1958 (dev_darcy_p,
1959 dev_darcy_grad_p));
1960 cudaDeviceSynchronize();
1961 checkForCudaErrorsIter("After findDarcyPressureGradient",
1962 iter);
1963
1964 KERNEL_LAUNCH(setDarcyGhostNodes<Float3>, dimGridFluid, dimBlockFluid, 0,
1965 (dev_darcy_grad_p,
1966 darcy.bc_xn, darcy.bc_xp,
1967 darcy.bc_yn, darcy.bc_yp,
1968 darcy.bc_bot, darcy.bc_top));
1969 cudaDeviceSynchronize();
1970 if (PROFILING == 1)
1971 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1972 &t_setDarcyGhostNodes);
1973 checkForCudaErrorsIter("Post setDarcyGhostNodes("
1974 "dev_darcy_grad_p)", iter);
1975
1976 /*if (PROFILING == 1)
1977 startTimer(&kernel_tic);
1978 findDarcyPorositiesLinear<<<dimGridFluid, dimBlockFluid>>>(
1979 dev_cellStart,
1980 dev_cellEnd,
1981 dev_x_sorted,
1982 dev_vel_sorted,
1983 iter,
1984 darcy.ndem,
1985 np,
1986 darcy.c_phi,
1987 dev_darcy_phi,
1988 dev_darcy_dphi,
1989 dev_darcy_div_v_p);
1990 cudaDeviceSynchronize();
1991 if (PROFILING == 1)
1992 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
1993 &t_findDarcyPorosities);
1994 checkForCudaErrorsIter("Post findDarcyPorosities", iter);*/
1995
1996 /*findDarcyPressureForce<<<dimGrid, dimBlock>>>(
1997 dev_x,
1998 dev_darcy_p,
1999 wall0_iz,
2000 darcy.rho_f,
2001 dev_force,
2002 dev_darcy_f_p);*/
2003 KERNEL_LAUNCH(findDarcyPressureForceLinear, dimGrid, dimBlock, 0,
2004 (dev_x,
2005 dev_darcy_grad_p,
2006 dev_darcy_phi,
2007 wall0_iz,
2008 darcy.rho_f,
2009 darcy.bc_top,
2010 dev_force,
2011 dev_darcy_f_p));
2012 cudaDeviceSynchronize();
2013 if (PROFILING == 1)
2014 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2015 &t_findDarcyPressureForce);
2016 checkForCudaErrorsIter("Post findDarcyPressureForce",
2017 iter);
2018 }
2019
2020 if ((iter % darcy.ndem) == 0) {
2021
2022 if (PROFILING == 1)
2023 startTimer(&kernel_tic);
2024 /*findDarcyPorosities<<<dimGridFluid, dimBlockFluid>>>(
2025 dev_cellStart,
2026 dev_cellEnd,
2027 dev_x_sorted,
2028 dev_vel_sorted,
2029 iter,
2030 darcy.ndem,
2031 np,
2032 darcy.c_phi,
2033 dev_darcy_phi,
2034 dev_darcy_dphi);*/
2035 KERNEL_LAUNCH(findDarcyPorositiesLinear, dimGridFluid, dimBlockFluid, 0,
2036 (dev_cellStart,
2037 dev_cellEnd,
2038 dev_x_sorted,
2039 dev_vel_sorted,
2040 iter,
2041 darcy.ndem,
2042 np,
2043 darcy.c_phi,
2044 dev_darcy_phi,
2045 dev_darcy_dphi,
2046 dev_darcy_div_v_p,
2047 dev_darcy_vp_avg));
2048 cudaDeviceSynchronize();
2049 if (PROFILING == 1)
2050 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2051 &t_findDarcyPorosities);
2052 checkForCudaErrorsIter("Post findDarcyPorosities", iter);
2053
2054 // copy porosities to the frictionless Y boundaries
2055 if (grid.periodic == 2) {
2056 KERNEL_LAUNCH(copyDarcyPorositiesToEdges, dimGridFluid, dimBlockFluid, 0,
2057 (dev_darcy_phi,
2058 dev_darcy_dphi,
2059 dev_darcy_div_v_p,
2060 dev_darcy_vp_avg));
2061 cudaDeviceSynchronize();
2062 }
2063
2064 // copy porosities to the frictionless lower Z boundary
2065 if (grid.periodic == 2) {
2066 KERNEL_LAUNCH(copyDarcyPorositiesToBottom, dimGridFluid, dimBlockFluid, 0,
2067 (dev_darcy_phi,
2068 dev_darcy_dphi,
2069 dev_darcy_div_v_p,
2070 dev_darcy_vp_avg));
2071 cudaDeviceSynchronize();
2072 }
2073
2074 // Modulate the pressures at the upper boundary cells
2075 if ((darcy.p_mod_A > 1.0e-5 || darcy.p_mod_A < -1.0e-5) &&
2076 darcy.p_mod_f > 1.0e-7) {
2077 // original pressure
2078 Float new_pressure =
2079 darcy.p_top_orig + darcy.p_mod_A
2080 *sin(2.0*M_PI*darcy.p_mod_f*time.current
2081 + darcy.p_mod_phi);
2082 if (PROFILING == 1)
2083 startTimer(&kernel_tic);
2084 KERNEL_LAUNCH(setDarcyTopPressure, dimGridFluid, dimBlockFluid, 0,
2085 (new_pressure,
2086 dev_darcy_p,
2087 wall0_iz));
2088 cudaDeviceSynchronize();
2089 checkForCudaErrorsIter("Post setUpperPressureNS", iter);
2090
2091 // Modulate the pressures at the top wall
2092 KERNEL_LAUNCH(setDarcyTopWallPressure, dimGridFluid, dimBlockFluid, 0,
2093 (new_pressure,
2094 wall0_iz,
2095 dev_darcy_p));
2096 cudaDeviceSynchronize();
2097 checkForCudaErrorsIter("Post setDarcyTopWallPressure",
2098 iter);
2099
2100 if (PROFILING == 1)
2101 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2102 &t_setDarcyTopPressure);
2103 }
2104
2105 if (PROFILING == 1)
2106 startTimer(&kernel_tic);
2107 KERNEL_LAUNCH(findDarcyPermeabilities, dimGridFluid, dimBlockFluid, 0,
2108 (darcy.k_c, dev_darcy_phi, dev_darcy_k));
2109 cudaDeviceSynchronize();
2110 if (PROFILING == 1)
2111 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2112 &t_findDarcyPermeabilities);
2113 checkForCudaErrorsIter("Post findDarcyPermeabilities",
2114 iter);
2115
2116 if (PROFILING == 1)
2117 startTimer(&kernel_tic);
2118 KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
2119 (dev_darcy_phi,
2120 darcy.bc_xn, darcy.bc_xp,
2121 darcy.bc_yn, darcy.bc_yp,
2122 darcy.bc_bot, darcy.bc_top));
2123 cudaDeviceSynchronize();
2124 if (PROFILING == 1)
2125 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2126 &t_setDarcyGhostNodes);
2127 checkForCudaErrorsIter(
2128 "Post setDarcyGhostNodes(dev_darcy_phi)", iter);
2129
2130 if (PROFILING == 1)
2131 startTimer(&kernel_tic);
2132 KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
2133 (dev_darcy_k,
2134 darcy.bc_xn, darcy.bc_xp,
2135 darcy.bc_yn, darcy.bc_yp,
2136 darcy.bc_bot, darcy.bc_top));
2137 cudaDeviceSynchronize();
2138 if (PROFILING == 1)
2139 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2140 &t_setDarcyGhostNodes);
2141 checkForCudaErrorsIter(
2142 "Post setDarcyGhostNodes(dev_darcy_k)", iter);
2143
2144 if (PROFILING == 1)
2145 startTimer(&kernel_tic);
2146 KERNEL_LAUNCH(findDarcyPermeabilityGradients, dimGridFluid, dimBlockFluid, 0,
2147 (dev_darcy_k, dev_darcy_grad_k));
2148 cudaDeviceSynchronize();
2149 if (PROFILING == 1)
2150 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2151 &t_findDarcyPermeabilityGradients);
2152 checkForCudaErrorsIter(
2153 "Post findDarcyPermeabilityGradients", iter);
2154
2155 if (iter == 0) {
2156 KERNEL_LAUNCH(setDarcyNormZero, dimGridFluid, dimBlockFluid, 0,
2157 (dev_darcy_norm));
2158 cudaDeviceSynchronize();
2159 checkForCudaErrorsIter("Post setDarcyNormZero", iter);
2160
2161 if (PROFILING == 1)
2162 startTimer(&kernel_tic);
2163 KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
2164 (dev_darcy_p,
2165 dev_darcy_p_old));
2166 cudaDeviceSynchronize();
2167 if (PROFILING == 1)
2168 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2169 &t_copyValues);
2170 checkForCudaErrorsIter("Post copyValues(p -> p_old)",
2171 iter);
2172 }
2173
2174 if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
2175 if (PROFILING == 1)
2176 startTimer(&kernel_tic);
2177 KERNEL_LAUNCH(setDarcyGhostNodesFlux<Float>, dimGridFluid, dimBlockFluid, 0,
2178 (dev_darcy_p,
2179 darcy.bc_bot,
2180 darcy.bc_top,
2181 darcy.bc_bot_flux,
2182 darcy.bc_top_flux,
2183 dev_darcy_k,
2184 darcy.mu));
2185 cudaDeviceSynchronize();
2186 if (PROFILING == 1)
2187 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2188 &t_setDarcyGhostNodes);
2189 checkForCudaErrorsIter(
2190 "Post setDarcyGhostNodesFlux", iter);
2191 }
2192
2193 // Solve the system of epsilon using a Jacobi iterative
2194 // solver. The average normalized residual is initialized
2195 // to a large value.
2196 //double avg_norm_res;
2197 double max_norm_res;
2198
2199 // Write a log file of the normalized residuals during the
2200 // Jacobi iterations
2201 std::ofstream reslog;
2202 if (write_res_log == 1)
2203 reslog.open("max_res_norm.dat");
2204
2205 for (unsigned int nijac = 0; nijac<darcy.maxiter; ++nijac) {
2206
2207 #if defined(REPORT_EPSILON) || defined(REPORT_FORCING_TERMS)
2208 std::cout << "\n\n### Jacobi iteration " << nijac << std::endl;
2209 #endif
2210
2211 if (nijac == 0) {
2212 if (PROFILING == 1)
2213 startTimer(&kernel_tic);
2214 KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
2215 (dev_darcy_p,
2216 dev_darcy_p_old));
2217 cudaDeviceSynchronize();
2218 if (PROFILING == 1)
2219 stopTimer(&kernel_tic, &kernel_toc,
2220 &kernel_elapsed, &t_copyValues);
2221 checkForCudaErrorsIter(
2222 "Post copyValues(p -> p_old)", iter);
2223 }
2224
2225 if (PROFILING == 1)
2226 startTimer(&kernel_tic);
2227 KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
2228 (dev_darcy_p,
2229 darcy.bc_xn, darcy.bc_xp,
2230 darcy.bc_yn, darcy.bc_yp,
2231 darcy.bc_bot, darcy.bc_top));
2232 cudaDeviceSynchronize();
2233 if (PROFILING == 1)
2234 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2235 &t_setDarcyGhostNodes);
2236 checkForCudaErrorsIter("Post setDarcyGhostNodes("
2237 "dev_darcy_p) in Jacobi loop", iter);
2238
2239 if (nijac == 0) {
2240 if (PROFILING == 1)
2241 startTimer(&kernel_tic);
2242 KERNEL_LAUNCH(firstDarcySolution, dimGridFluid, dimBlockFluid, 0,
2243 (dev_darcy_p,
2244 dev_darcy_k,
2245 dev_darcy_phi,
2246 dev_darcy_dphi,
2247 dev_darcy_div_v_p,
2248 dev_darcy_vp_avg,
2249 dev_darcy_grad_k,
2250 darcy.beta_f,
2251 darcy.mu,
2252 darcy.bc_xn,
2253 darcy.bc_xp,
2254 darcy.bc_yn,
2255 darcy.bc_yp,
2256 darcy.bc_bot,
2257 darcy.bc_top,
2258 darcy.ndem,
2259 wall0_iz,
2260 dev_darcy_p_constant,
2261 dev_darcy_dp_expl));
2262 cudaDeviceSynchronize();
2263 if (PROFILING == 1)
2264 stopTimer(&kernel_tic, &kernel_toc,
2265 &kernel_elapsed,
2266 &t_updateDarcySolution);
2267 checkForCudaErrorsIter("Post updateDarcySolution",
2268 iter);
2269 }
2270
2271 if (PROFILING == 1)
2272 startTimer(&kernel_tic);
2273 KERNEL_LAUNCH(updateDarcySolution, dimGridFluid, dimBlockFluid, 0,
2274 (dev_darcy_p_old,
2275 //dev_darcy_dpdt,
2276 dev_darcy_dp_expl,
2277 dev_darcy_p,
2278 dev_darcy_k,
2279 dev_darcy_phi,
2280 dev_darcy_dphi,
2281 dev_darcy_div_v_p,
2282 dev_darcy_vp_avg,
2283 dev_darcy_grad_k,
2284 darcy.beta_f,
2285 darcy.mu,
2286 darcy.bc_xn,
2287 darcy.bc_xp,
2288 darcy.bc_yn,
2289 darcy.bc_yp,
2290 darcy.bc_bot,
2291 darcy.bc_top,
2292 darcy.ndem,
2293 wall0_iz,
2294 dev_darcy_p_constant,
2295 dev_darcy_p_new,
2296 dev_darcy_norm));
2297 cudaDeviceSynchronize();
2298 if (PROFILING == 1)
2299 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2300 &t_updateDarcySolution);
2301 checkForCudaErrorsIter("Post updateDarcySolution",
2302 iter);
2303
2304 if (darcy.bc_top == 1) {
2305 if (PROFILING == 1)
2306 startTimer(&kernel_tic);
2307 KERNEL_LAUNCH(setDarcyTopWallFixedFlow, dimGridFluid, dimBlockFluid, 0,
2308 (wall0_iz, dev_darcy_p));
2309 cudaDeviceSynchronize();
2310 if (PROFILING == 1)
2311 stopTimer(&kernel_tic, &kernel_toc,
2312 &kernel_elapsed,
2313 &t_updateDarcySolution);
2314 checkForCudaErrorsIter(
2315 "Post setDarcyTopWallFixedFlow", iter);
2316 }
2317
2318 if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
2319 if (PROFILING == 1)
2320 startTimer(&kernel_tic);
2321 KERNEL_LAUNCH(setDarcyGhostNodesFlux<Float>, dimGridFluid, dimBlockFluid, 0,
2322 (dev_darcy_p,
2323 darcy.bc_bot,
2324 darcy.bc_top,
2325 darcy.bc_bot_flux,
2326 darcy.bc_top_flux,
2327 dev_darcy_k,
2328 darcy.mu));
2329 cudaDeviceSynchronize();
2330 if (PROFILING == 1)
2331 stopTimer(&kernel_tic, &kernel_toc,
2332 &kernel_elapsed,
2333 &t_setDarcyGhostNodes);
2334 checkForCudaErrorsIter(
2335 "Post setDarcyGhostNodesFlux", iter);
2336 }
2337
2338 // Copy new values to current values
2339 if (PROFILING == 1)
2340 startTimer(&kernel_tic);
2341 KERNEL_LAUNCH(copyValues<Float>, dimGridFluid, dimBlockFluid, 0,
2342 (dev_darcy_p_new,
2343 dev_darcy_p));
2344 cudaDeviceSynchronize();
2345 if (PROFILING == 1)
2346 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2347 &t_copyValues);
2348 checkForCudaErrorsIter("Post copyValues(p_new -> p)",
2349 iter);
2350
2351 #ifdef REPORT_EPSILON
2352 std::cout << "\n###### JACOBI ITERATION "
2353 << nijac << " after copyValues ######" << std::endl;
2354 transferDarcyPressuresFromGlobalDeviceMemory();
2355 printDarcyArray(stdout, darcy.p, "p");
2356 #endif
2357
2358 if (nijac % nijacnorm == 0) {
2359 // Read the normalized residuals from the device
2360 transferDarcyNormFromGlobalDeviceMemory();
2361
2362 // Write the normalized residuals to the terminal
2363 //printDarcyArray(stdout, darcy.norm, "norm");
2364
2365 // Find the maximum value of the normalized
2366 // residuals
2367 max_norm_res = maxNormResDarcy();
2368
2369 // Write the Jacobi iteration number and maximum
2370 // value of the normalized residual to the log file
2371 if (write_res_log == 1)
2372 reslog << nijac << '\t' << max_norm_res
2373 << std::endl;
2374
2375 if (max_norm_res <= darcy.tolerance) {
2376 if (write_conv_log == 1
2377 && iter % conv_log_interval == 0)
2378 convlog << iter+1 << '\t' << nijac
2379 << std::endl;
2380
2381 break; // solution has converged
2382 }
2383 }
2384
2385 if (nijac == darcy.maxiter-1) {
2386
2387 if (write_conv_log == 1)
2388 convlog << iter+1 << '\t' << nijac << std::endl;
2389
2390 std::cerr << "\nIteration " << iter << ", time "
2391 << iter*time.dt << " s: "
2392 "Error, the pressure solution in the fluid "
2393 "calculations did not converge. Try increasing "
2394 "the value of 'darcy.maxiter' ("
2395 << darcy.maxiter
2396 << ") or increase 'darcy.tolerance' ("
2397 << darcy.tolerance << ")." << std::endl;
2398 }
2399
2400 if (write_res_log == 1)
2401 reslog.close();
2402
2403 //break; // end after first iteration
2404 }
2405
2406 // Zero all dphi values right after they are used in fluid
2407 // solution, unless a file is written in this step.
2408 if (filetimeclock + time.dt < time.file_dt) {
2409 KERNEL_LAUNCH(setDarcyZeros<Float>, dimGridFluid, dimBlockFluid, 0,
2410 (dev_darcy_dphi));
2411 cudaDeviceSynchronize();
2412 checkForCudaErrorsIter(
2413 "After setDarcyZeros(dev_darcy_dphi)", iter);
2414 }
2415
2416 if (PROFILING == 1)
2417 startTimer(&kernel_tic);
2418 KERNEL_LAUNCH(setDarcyGhostNodes<Float>, dimGridFluid, dimBlockFluid, 0,
2419 (dev_darcy_p,
2420 darcy.bc_xn, darcy.bc_xp,
2421 darcy.bc_yn, darcy.bc_yp,
2422 darcy.bc_bot, darcy.bc_top));
2423 cudaDeviceSynchronize();
2424 if (PROFILING == 1)
2425 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2426 &t_setDarcyGhostNodes);
2427 checkForCudaErrorsIter("Post setDarcyGhostNodes("
2428 "dev_darcy_p) after Jacobi loop", iter);
2429
2430 if (PROFILING == 1)
2431 startTimer(&kernel_tic);
2432 KERNEL_LAUNCH(findDarcyVelocities, dimGridFluid, dimBlockFluid, 0,
2433 (dev_darcy_p,
2434 dev_darcy_phi,
2435 dev_darcy_k,
2436 darcy.mu,
2437 dev_darcy_v));
2438 cudaDeviceSynchronize();
2439 if (PROFILING == 1)
2440 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2441 &t_findDarcyVelocities);
2442 checkForCudaErrorsIter("Post findDarcyVelocities", iter);
2443 }
2444 }
2445 }
2446 //break; // end after first iteration
2447
2448 if (np > 0) {
2449
2450 // Find shear stresses on upper fixed particles if a shear stress BC
2451 // is specified (wmode[0] == 3)
2452 if (walls.nw > 0 && walls.wmode[0] == 3) {
2453
2454 if (PROFILING == 1)
2455 startTimer(&kernel_tic);
2456 KERNEL_LAUNCH(findShearStressOnFixedMovingParticles, dimGrid, dimBlock, 0,
2457 (dev_x,
2458 dev_vel,
2459 dev_force,
2460 dev_walls_tau_eff_x_pp));
2461 cudaDeviceSynchronize();
2462 if (PROFILING == 1)
2463 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2464 &t_summation);
2465 checkForCudaErrorsIter(
2466 "Post findShearStressOnFixedMovingParticles", iter);
2467
2468 if (PROFILING == 1)
2469 startTimer(&kernel_tic);
2470 KERNEL_LAUNCH(summation, dimGrid, dimBlock, 0,
2471 (dev_walls_tau_eff_x_pp,
2472 dev_walls_tau_eff_x_partial));
2473 cudaDeviceSynchronize();
2474 if (PROFILING == 1)
2475 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2476 &t_summation);
2477 checkForCudaErrorsIter("Post shear stress summation", iter);
2478 }
2479
2480 // Determine whether it is time to step the velocity
2481 if (time.current >= v2_start && time.current < v2_end &&
2482 velocity_state == 1) {
2483 change_velocity_state = 1.0;
2484 velocity_state = 2;
2485 } else if (time.current >= v2_end && velocity_state == 2) {
2486 change_velocity_state = -1.0;
2487 velocity_state = 1;
2488 }
2489
2490 // Update particle kinematics
2491 if (PROFILING == 1)
2492 startTimer(&kernel_tic);
2493 KERNEL_LAUNCH(integrate, dimGrid, dimBlock, 0,
2494 (dev_x_sorted,
2495 dev_vel_sorted,
2496 dev_angvel_sorted,
2497 dev_x,
2498 dev_vel,
2499 dev_angvel,
2500 dev_force,
2501 dev_torque,
2502 dev_angpos,
2503 dev_acc,
2504 dev_angacc,
2505 dev_vel0,
2506 dev_angvel0,
2507 dev_xyzsum,
2508 dev_gridParticleIndex,
2509 iter,
2510 dev_walls_wmode,
2511 dev_walls_mvfd,
2512 dev_walls_tau_eff_x_partial,
2513 dev_walls_tau_x,
2514 walls.tau_x[0],
2515 change_velocity_state,
2516 velocity_factor,
2517 blocksPerGrid));
2518 cudaDeviceSynchronize();
2519 checkForCudaErrorsIter("Post integrate", iter);
2520 if (PROFILING == 1)
2521 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2522 &t_integrate);
2523
2524 if (change_velocity_state != 0)
2525 change_velocity_state = 0;
2526
2527 // Summation of forces on wall
2528 if (PROFILING == 1)
2529 startTimer(&kernel_tic);
2530 if (walls.nw > 0) {
2531 KERNEL_LAUNCH(summation, dimGrid, dimBlock, 0,
2532 (dev_walls_force_pp,
2533 dev_walls_force_partial));
2534 }
2535 cudaDeviceSynchronize();
2536 if (PROFILING == 1)
2537 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2538 &t_summation);
2539 checkForCudaErrorsIter("Post wall force summation", iter);
2540
2541 // Update wall kinematics
2542 if (PROFILING == 1)
2543 startTimer(&kernel_tic);
2544 if (walls.nw > 0) {
2545 KERNEL_LAUNCH(integrateWalls, 1, walls.nw, 0,
2546 (dev_walls_nx,
2547 dev_walls_mvfd,
2548 dev_walls_wmode,
2549 dev_walls_force_partial,
2550 dev_walls_acc,
2551 blocksPerGrid,
2552 time.current,
2553 iter));
2554 }
2555 cudaDeviceSynchronize();
2556 if (PROFILING == 1)
2557 stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
2558 &t_integrateWalls);
2559 checkForCudaErrorsIter("Post integrateWalls", iter);
2560 }
2561
2562 // Update timers and counters
2563 //time.current = iter*time.dt;
2564 time.current += time.dt;
2565 filetimeclock += time.dt;
2566 ++iter;
2567
2568 // Make sure all preceding tasks are complete
2569 if (cudaDeviceSynchronize() != cudaSuccess) {
2570 cerr << "Error during cudaDeviceSynchronize()" << endl;
2571 }
2572
2573 // Report time to console
2574 if (verbose == 1 && (iter % stdout_report == 0)) {
2575
2576 toc = clock();
2577 time_spent = (toc - tic)/(CLOCKS_PER_SEC); // real time spent
2578
2579 // Real time it takes to compute a second of model time
2580 t_ratio = time_spent/(time.current - t_start);
2581 time_t estimated_seconds_left(t_ratio*(time.total - time.current));
2582 tm *time_eta = gmtime(&estimated_seconds_left);
2583
2584 cout << "\r Current time: " << time.current << "/"
2585 << time.total << " s. ("
2586 << t_ratio << " s_real/s_sim, ETA: "
2587 << time_eta->tm_yday << "d "
2588 << std::setw(2) << std::setfill('0') << time_eta->tm_hour << ":"
2589 << std::setw(2) << std::setfill('0') << time_eta->tm_min << ":"
2590 << std::setw(2) << std::setfill('0') << time_eta->tm_sec
2591 << ") "; // << std::flush;
2592 }
2593
2594
2595 // Produce output binary if the time interval
2596 // between output files has been reached
2597 if (filetimeclock >= time.file_dt) {
2598
2599 // Pause the CPU thread until all CUDA calls previously issued are
2600 // completed
2601 cudaDeviceSynchronize();
2602 checkForCudaErrorsIter("Beginning of file output section", iter);
2603
2604 // v_x, v_y, v_z -> v
2605 if (fluid == 1 && cfd_solver == 0) {
2606 KERNEL_LAUNCH(interpolateFaceToCenter, dimGridFluid, dimBlockFluid, 0,
2607 (dev_ns_v_x,
2608 dev_ns_v_y,
2609 dev_ns_v_z,
2610 dev_ns_v));
2611 cudaDeviceSynchronize();
2612 checkForCudaErrorsIter("Post interpolateFaceToCenter", iter);
2613 }
2614
2615 //// Copy device data to host memory
2616 transferFromGlobalDeviceMemory();
2617 checkForCudaErrorsIter("After transferFromGlobalDeviceMemory()",
2618 iter);
2619
2620 // Empty the dphi values after device to host transfer
2621 if (fluid == 1) {
2622 if (cfd_solver == 1) {
2623 KERNEL_LAUNCH(setDarcyZeros<Float>, dimGridFluid, dimBlockFluid, 0,
2624 (dev_darcy_dphi));
2625 cudaDeviceSynchronize();
2626 checkForCudaErrorsIter(
2627 "After setDarcyZeros(dev_darcy_dphi) after transfer",
2628 iter);
2629 }
2630 }
2631
2632 // Pause the CPU thread until all CUDA calls previously issued are
2633 // completed
2634 cudaDeviceSynchronize();
2635
2636 // Check the numerical stability of the NS solver
2637 if (fluid == 1)
2638 if (cfd_solver == 0)
2639 checkNSstability();
2640
2641 // Write binary output file
2642 time.step_count += 1;
2643 snprintf(file, sizeof(file), "output/%s.output%05d.bin",
2644 sid.c_str(), time.step_count);
2645 writebin(file);
2646
2647 /*std::cout
2648 << "\n###### OUTPUT FILE " << time.step_count << " ######"
2649 << std::endl;
2650 transferNSepsilonFromGlobalDeviceMemory();
2651 printNSarray(stdout, ns.epsilon, "epsilon");*/
2652
2653 // Write fluid arrays
2654 /*if (fluid == 1 && cfd_solver == 0) {
2655 sprintf(file,"output/%s.ns_phi.output%05d.bin", sid.c_str(),
2656 time.step_count);
2657 writeNSarray(ns.phi, file);
2658 }*/
2659
2660 if (CONTACTINFO == 1) {
2661 // Write contact information to stdout
2662 cout << "\n\n---------------------------\n"
2663 << "t = " << time.current << " s.\n"
2664 << "---------------------------\n";
2665
2666 for (int n = 0; n < np; ++n) {
2667 cout << "\n## Particle " << n << " ##\n";
2668
2669 cout << "- contacts:\n";
2670 for (int nc = 0; nc < NC; ++nc)
2671 cout << "[" << nc << "]=" << k.contacts[nc+NC*n] <<
2672 '\n';
2673
2674 cout << "\n- delta_t:\n";
2675 for (int nc = 0; nc < NC; ++nc)
2676 cout << k.delta_t[nc+NC*n].x << '\t'
2677 << k.delta_t[nc+NC*n].y << '\t'
2678 << k.delta_t[nc+NC*n].z << '\t'
2679 << k.delta_t[nc+NC*n].w << '\n';
2680
2681 cout << "\n- distmod:\n";
2682 for (int nc = 0; nc < NC; ++nc)
2683 cout << k.distmod[nc+NC*n].x << '\t'
2684 << k.distmod[nc+NC*n].y << '\t'
2685 << k.distmod[nc+NC*n].z << '\t'
2686 << k.distmod[nc+NC*n].w << '\n';
2687 }
2688 cout << '\n';
2689 }
2690
2691 // Update status.dat at the interval of filetime
2692 outfile = "output/" + sid + ".status.dat";
2693 fp = fopen(outfile.c_str(), "w");
2694 fprintf(fp,"%2.4e %2.4e %d\n",
2695 time.current,
2696 100.0*time.current/time.total,
2697 time.step_count);
2698 fclose(fp);
2699
2700 filetimeclock = 0.0;
2701 }
2702
2703 // Uncomment break command to stop after the first iteration
2704 //break;
2705 }
2706
2707 if (write_conv_log == 1)
2708 convlog.close();
2709
2710
2711 // Stop clock and display calculation time spent
2712 toc = clock();
2713 cudaEventRecord(dev_toc, 0);
2714 cudaEventSynchronize(dev_toc);
2715
2716 time_spent = (toc - tic)/(CLOCKS_PER_SEC);
2717 cudaEventElapsedTime(&dev_time_spent, dev_tic, dev_toc);
2718
2719 if (verbose == 1) {
2720 cout << "\nSimulation ended. Statistics:\n"
2721 << " - Last output file number: "
2722 << time.step_count << "\n"
2723 << " - GPU time spent: "
2724 << dev_time_spent/1000.0f << " s\n"
2725 << " - CPU time spent: "
2726 << time_spent << " s\n"
2727 << " - Mean duration of iteration:\n"
2728 << " " << dev_time_spent/((double)iter*1000.0f) << " s"
2729 << std::endl;
2730 }
2731
2732 cudaEventDestroy(dev_tic);
2733 cudaEventDestroy(dev_toc);
2734
2735 cudaEventDestroy(kernel_tic);
2736 cudaEventDestroy(kernel_toc);
2737
2738 // Report time spent on each kernel
2739 if (PROFILING == 1 && verbose == 1) {
2740 double t_sum = t_calcParticleCellID + t_thrustsort + t_reorderArrays +
2741 t_topology + t_interact + t_bondsLinear + t_latticeBoltzmannD3Q19 +
2742 t_integrate + t_summation + t_integrateWalls + t_findPorositiesDev +
2743 t_findNSstressTensor +
2744 t_findNSdivphiviv + t_findNSdivtau + t_findPredNSvelocities +
2745 t_setNSepsilon + t_setNSdirichlet + t_setNSghostNodesDev +
2746 t_findNSforcing + t_jacobiIterationNS + t_updateNSvelocityPressure +
2747 t_findDarcyPorosities + t_setDarcyGhostNodes +
2748 t_findDarcyPressureForce + t_setDarcyTopPressure +
2749 t_findDarcyPermeabilities + t_findDarcyPermeabilityGradients +
2750 //t_findDarcyPressureChange +
2751 t_updateDarcySolution + t_copyValues + t_findDarcyVelocities;
2752
2753 cout << "\nKernel profiling statistics:\n"
2754 << " - calcParticleCellID:\t\t" << t_calcParticleCellID/1000.0
2755 << " s"
2756 << "\t(" << 100.0*t_calcParticleCellID/t_sum << " %)\n"
2757 << " - thrustsort:\t\t\t" << t_thrustsort/1000.0 << " s"
2758 << "\t(" << 100.0*t_thrustsort/t_sum << " %)\n"
2759 << " - reorderArrays:\t\t" << t_reorderArrays/1000.0 << " s"
2760 << "\t(" << 100.0*t_reorderArrays/t_sum << " %)\n";
2761 if (params.contactmodel == 2 || params.contactmodel == 3) {
2762 cout
2763 << " - topology:\t\t\t" << t_topology/1000.0 << " s"
2764 << "\t(" << 100.0*t_topology/t_sum << " %)\n";
2765 }
2766 cout << " - interact:\t\t\t" << t_interact/1000.0 << " s"
2767 << "\t(" << 100.0*t_interact/t_sum << " %)\n";
2768 if (params.nb0 > 0) {
2769 cout << " - bondsLinear:\t\t" << t_bondsLinear/1000.0 << " s"
2770 << "\t(" << 100.0*t_bondsLinear/t_sum << " %)\n";
2771 }
2772 cout << " - integrate:\t\t\t" << t_integrate/1000.0 << " s"
2773 << "\t(" << 100.0*t_integrate/t_sum << " %)\n"
2774 << " - summation:\t\t\t" << t_summation/1000.0 << " s"
2775 << "\t(" << 100.0*t_summation/t_sum << " %)\n"
2776 << " - integrateWalls:\t\t" << t_integrateWalls/1000.0 << " s"
2777 << "\t(" << 100.0*t_integrateWalls/t_sum << " %)\n";
2778 if (fluid == 1 && cfd_solver == 0) {
2779 cout << " - findPorositiesDev:\t\t" << t_findPorositiesDev/1000.0
2780 << " s" << "\t(" << 100.0*t_findPorositiesDev/t_sum << " %)\n"
2781 << " - findNSstressTensor:\t\t" << t_findNSstressTensor/1000.0
2782 << " s" << "\t(" << 100.0*t_findNSstressTensor/t_sum << " %)\n"
2783 << " - findNSdivphiviv:\t\t" << t_findNSdivphiviv/1000.0
2784 << " s" << "\t(" << 100.0*t_findNSdivphiviv/t_sum << " %)\n"
2785 << " - findNSdivtau:\t\t" << t_findNSdivtau/1000.0
2786 << " s" << "\t(" << 100.0*t_findNSdivtau/t_sum << " %)\n"
2787 << " - findPredNSvelocities:\t" <<
2788 t_findPredNSvelocities/1000.0 << " s" << "\t(" <<
2789 100.0*t_findPredNSvelocities/t_sum << " %)\n"
2790 << " - setNSepsilon:\t\t" << t_setNSepsilon/1000.0
2791 << " s" << "\t(" << 100.0*t_setNSepsilon/t_sum << " %)\n"
2792 << " - setNSdirichlet:\t\t" << t_setNSdirichlet/1000.0
2793 << " s" << "\t(" << 100.0*t_setNSdirichlet/t_sum << " %)\n"
2794 << " - setNSghostNodesDev:\t\t" << t_setNSghostNodesDev/1000.0
2795 << " s" << "\t(" << 100.0*t_setNSghostNodesDev/t_sum << " %)\n"
2796 << " - findNSforcing:\t\t" << t_findNSforcing/1000.0 << " s"
2797 << "\t(" << 100.0*t_findNSforcing/t_sum << " %)\n"
2798 << " - jacobiIterationNS:\t\t" << t_jacobiIterationNS/1000.0
2799 << " s"
2800 << "\t(" << 100.0*t_jacobiIterationNS/t_sum << " %)\n"
2801 << " - updateNSvelocityPressure:\t"
2802 << t_updateNSvelocityPressure/1000.0 << " s"
2803 << "\t(" << 100.0*t_updateNSvelocityPressure/t_sum << " %)\n";
2804 } else if (fluid == 1 && cfd_solver == 1) {
2805 cout << " - findDarcyPorosities:\t" <<
2806 t_findDarcyPorosities/1000.0 << " s" << "\t(" <<
2807 100.0*t_findDarcyPorosities/t_sum << " %)\n"
2808 << " - setDarcyGhostNodes:\t\t" <<
2809 t_setDarcyGhostNodes/1000.0 << " s" << "\t(" <<
2810 100.0*t_setDarcyGhostNodes/t_sum << " %)\n"
2811 << " - findDarcyPressureForce:\t" <<
2812 t_findDarcyPressureForce/1000.0 << " s" << "\t(" <<
2813 100.0*t_findDarcyPressureForce/t_sum << " %)\n"
2814 << " - setDarcyTopPressure:\t" <<
2815 t_setDarcyTopPressure/1000.0 << " s" << "\t(" <<
2816 100.0*t_setDarcyTopPressure/t_sum << " %)\n"
2817 << " - findDarcyPermeabilities:\t" <<
2818 t_findDarcyPermeabilities/1000.0 << " s" << "\t(" <<
2819 100.0*t_findDarcyPermeabilities/t_sum << " %)\n"
2820 << " - findDarcyPermeabilityGrads:\t" <<
2821 t_findDarcyPermeabilityGradients/1000.0 << " s" << "\t(" <<
2822 100.0*t_findDarcyPermeabilityGradients/t_sum << " %)\n"
2823 //<< " - findDarcyPressureChange:\t" <<
2824 //t_findDarcyPressureChange/1000.0 << " s" << "\t(" <<
2825 //100.0*t_findDarcyPressureChange/t_sum << " %)\n"
2826 << " - updateDarcySolution:\t" <<
2827 t_updateDarcySolution/1000.0 << " s" << "\t(" <<
2828 100.0*t_updateDarcySolution/t_sum << " %)\n"
2829 << " - copyValues:\t\t\t" <<
2830 t_copyValues/1000.0 << " s" << "\t(" <<
2831 100.0*t_copyValues/t_sum << " %)\n"
2832 << " - findDarcyVelocities:\t" <<
2833 t_findDarcyVelocities/1000.0 << " s" << "\t(" <<
2834 100.0*t_findDarcyVelocities/t_sum << " %)" << std::endl;
2835 }
2836 }
2837
2838 // Free GPU device memory
2839 freeGlobalDeviceMemory();
2840 checkForCudaErrorsIter("After freeGlobalDeviceMemory()", iter);
2841
2842 // Free contact info arrays
2843 delete[] k.contacts;
2844 delete[] k.distmod;
2845 delete[] k.delta_t;
2846
2847 if (fluid == 1) {
2848 if (cfd_solver == 0)
2849 endNS();
2850 else if (cfd_solver == 1)
2851 endDarcy();
2852 }
2853
2854 cudaDeviceReset();
2855 }
2856 // vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4