tMoved CPU utility functions to separate utility.h - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit 819e1feb3c46bdb460501c4f73743ec7a9b5682c
(DIR) parent 72813c78d9d2c75a835a58382f4219361bb1c41b
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Wed, 12 Jun 2013 12:29:35 +0200
Moved CPU utility functions to separate utility.h
Diffstat:
M src/darcy.cpp | 124 ++++++++++++++++++++-----------
M src/device.cu | 9 +++++----
M src/latticeboltzmann.cuh | 10 ++--------
M src/sphere.h | 3 ++-
A src/utility.h | 21 +++++++++++++++++++++
5 files changed, 112 insertions(+), 55 deletions(-)
---
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -8,12 +8,14 @@
#include "datatypes.h"
#include "constants.h"
#include "sphere.h"
+#include "utility.h"
// Initialize memory
void DEM::initDarcyMem()
{
unsigned int ncells = d_nx*d_ny*d_nz;
d_H = new Float[ncells]; // hydraulic pressure matrix
+ d_H_new = new Float[ncells]; // hydraulic pressure matrix
d_V = new Float3[ncells]; // Cell hydraulic velocity
d_dH = new Float3[ncells]; // Cell spatial gradient in hydraulic pressures
d_K = new Float[ncells]; // hydraulic conductivity matrix
t@@ -26,6 +28,7 @@ void DEM::initDarcyMem()
void DEM::freeDarcyMem()
{
free(d_H);
+ free(d_H_new);
free(d_V);
free(d_dH);
free(d_K);
t@@ -34,6 +37,14 @@ void DEM::freeDarcyMem()
free(d_n);
}
+// Swap two arrays pointers
+void swapFloatArrays(Float* arr1, Float* arr2)
+{
+ Float* tmp = arr1;
+ arr1 = arr2;
+ arr2 = tmp;
+}
+
// 3D index to 1D index
unsigned int DEM::idx(
const unsigned int x,
t@@ -84,39 +95,6 @@ Float DEM::minVal3dArr(Float* arr)
}
}
-// Find the spatial gradient in pressures per cell
-void DEM::findDarcyGradients()
-{
- // Cell sizes squared
- const Float dx2 = d_dx*d_dx;
- const Float dy2 = d_dy*d_dy;
- const Float dz2 = d_dz*d_dz;
-
- Float H;
- unsigned int ix, iy, iz, cellidx;
-
- for (ix=1; ix<d_nx-1; ++ix) {
- for (iy=1; iy<d_ny-1; ++iy) {
- for (iz=1; iz<d_nz-1; ++iz) {
-
- cellidx = idx(ix,iy,iz);
-
- H = d_H[cellidx]; // cell hydraulic pressure
-
- // Second order central difference
- d_dH[cellidx].x
- = (d_H[idx(ix+1,iy,iz)] - 2.0*H + d_H[idx(ix-1,iy,iz)])/dx2;
-
- d_dH[cellidx].y
- = (d_H[idx(ix,iy+1,iz)] - 2.0*H + d_H[idx(ix,iy-1,iz)])/dy2;
-
- d_dH[cellidx].z
- = (d_H[idx(ix,iy,iz+1)] - 2.0*H + d_H[idx(ix,iy,iz-1)])/dz2;
- }
- }
- }
-}
-
// Set the gradient to 0.0 in all dimensions at the boundaries
void DEM::setDarcyBCNeumannZero()
t@@ -155,6 +133,40 @@ void DEM::setDarcyBCNeumannZero()
}
+// Find the spatial gradient in pressures per cell
+void DEM::findDarcyGradients()
+{
+ // Cell sizes squared
+ const Float dx2 = d_dx*d_dx;
+ const Float dy2 = d_dy*d_dy;
+ const Float dz2 = d_dz*d_dz;
+
+ Float H;
+ unsigned int ix, iy, iz, cellidx;
+
+ for (ix=1; ix<d_nx-1; ++ix) {
+ for (iy=1; iy<d_ny-1; ++iy) {
+ for (iz=1; iz<d_nz-1; ++iz) {
+
+ cellidx = idx(ix,iy,iz);
+
+ H = d_H[cellidx]; // cell hydraulic pressure
+
+ // Second order central difference
+ d_dH[cellidx].x
+ = (d_H[idx(ix+1,iy,iz)] - 2.0*H + d_H[idx(ix-1,iy,iz)])/dx2;
+
+ d_dH[cellidx].y
+ = (d_H[idx(ix,iy+1,iz)] - 2.0*H + d_H[idx(ix,iy-1,iz)])/dy2;
+
+ d_dH[cellidx].z
+ = (d_H[idx(ix,iy,iz+1)] - 2.0*H + d_H[idx(ix,iy,iz-1)])/dz2;
+ }
+ }
+ }
+}
+
+
// Perform an explicit step.
// Boundary conditions are fixed gradient values (Neumann)
void DEM::explDarcyStep()
t@@ -171,25 +183,53 @@ void DEM::explDarcyStep()
// Explicit 3D finite difference scheme
// new = old + gradient*timestep
unsigned int ix, iy, iz, cellidx;
- Float K, H;
+ Float K, H, d;
+ Float T, Tx, Ty, Tz, S;
Float dt = time.dt;
- for (ix=0; ix<d_nx; ++ix) {
+ /*for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
- for (iz=0; iz<d_nz; ++iz) {
+ for (iz=0; iz<d_nz; ++iz) {*/
+ for (ix=1; ix<d_nx-1; ++ix) {
+ for (iy=1; iy<d_ny-1; ++iy) {
+ for (iz=1; iz<d_nz-1; ++iz) {
cellidx = idx(ix,iy,iz);
K = d_K[cellidx]; // cell hydraulic conductivity
- H = d_H[cellidx]; // cell hydraulic pressure
+ //H = d_H[cellidx]; // cell hydraulic pressure
+
+ // Cell size
+ d = d_dx;
- d_H[cellidx]
- += d_W[cellidx]*dt // cell recharge
- + K*dt * // diffusivity term
- (d_dH[cellidx].x + d_dH[cellidx].y + d_dH[cellidx].z);
+ // Cell hydraulic transmissivity (as long as nx=ny=nz)
+ T = K*d;
+
+ // Transmissivity (arithmetic mean, harmonic mean may be better)
+ Tx = (d_K[idx(ix-1,iy,iz)]*d + T + d_K[idx(ix+1,iy,iz)])/3.0;
+ Ty = (d_K[idx(ix,iy-1,iz)]*d + T + d_K[idx(ix,iy+1,iz)])/3.0;
+ Tz = (d_K[idx(ix,iy,iz-1)]*d + T + d_K[idx(ix,iy,iz+1)])/3.0;
+
+ // Cell hydraulic storativity (as long as nx=ny=nz)
+ S = d_S[cellidx]*d; // not d^3 ?
+
+ d_H_new[cellidx] = d_H[cellidx] +
+ dt/S *
+ (Tx*d_dH[cellidx].x
+ + Ty*d_dH[cellidx].y
+ + Tz*d_dH[cellidx].z
+ + d_W[cellidx]); // Cell recharge
+
+ //d_H[cellidx]
+ //+= d_W[cellidx]*dt // cell recharge
+ //+ K*dt * // diffusivity term
+ //(d_dH[cellidx].x + d_dH[cellidx].y + d_dH[cellidx].z);
}
}
}
+ // Swap d_H and d_H_new
+ swapFloatArrays(d_H, d_H_new);
+
// Find macroscopic cell fluid velocities
findDarcyVelocities();
}
t@@ -411,7 +451,7 @@ void DEM::initDarcy(const Float cellsizemultiplier)
// Print final heads and free memory
void DEM::endDarcy()
{
- //printDarcyArray(stdout, d_H, "d_H");
+ printDarcyArray(stdout, d_H, "d_H");
//printDarcyArray(stdout, d_V, "d_V");
//printDarcyArray(stdout, d_n, "d_n");
freeDarcyMem();
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -14,6 +14,7 @@
#include "sphere.h"
#include "datatypes.h"
#include "utility.cuh"
+#include "utility.h"
#include "constants.cuh"
#include "debug.h"
t@@ -375,9 +376,9 @@ __host__ void DEM::freeGlobalDeviceMemory()
}
-__host__ void DEM::transferToGlobalDeviceMemory()
+__host__ void DEM::transferToGlobalDeviceMemory(int statusmsg)
{
- if (verbose == 1)
+ if (verbose == 1 && statusmsg == 1)
std::cout << " Transfering data to the device: ";
// Commonly-used memory sizes
t@@ -460,7 +461,7 @@ __host__ void DEM::transferToGlobalDeviceMemory()
}
checkForCudaErrors("End of transferToGlobalDeviceMemory");
- if (verbose == 1)
+ if (verbose == 1 && statusmsg == 1)
std::cout << "Done" << std::endl;
}
t@@ -869,7 +870,7 @@ __host__ void DEM::startTime()
explDarcyStep();
// Transfer data from host to device memory
- transferToGlobalDeviceMemory();
+ transferToGlobalDeviceMemory(0);
// Pause the CPU thread until all CUDA calls previously issued are completed
cudaThreadSynchronize();
(DIR) diff --git a/src/latticeboltzmann.cuh b/src/latticeboltzmann.cuh
t@@ -1,6 +1,8 @@
#ifndef LATTICEBOLTZMANN_CUH_
#define LATTICEBOLTZMANN_CUH_
+#include "utility.h"
+
// Enable line below to perform lattice Boltzmann computations on the
// GPU, disable for CPU computation
//#define LBM_GPU
t@@ -174,14 +176,6 @@ void initFluid(
#endif
}
-// Swap two arrays pointers
-void swapFloatArrays(Float* arr1, Float* arr2)
-{
- Float* tmp = arr1;
- arr1 = arr2;
- arr2 = tmp;
-}
-
// Combined streaming and collision step with particle coupling and optional
// periodic boundaries. Derived from A. Monitzer 2013
#ifdef LBM_GPU
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -120,7 +120,7 @@ class DEM {
void rt_freeGlobalDeviceMemory(void);
// Copy non-constant data to global GPU memory
- void transferToGlobalDeviceMemory(void);
+ void transferToGlobalDeviceMemory(int status = 1);
// Copy non-constant data from global GPU memory to host RAM
void transferFromGlobalDeviceMemory(void);
t@@ -151,6 +151,7 @@ class DEM {
int d_nx, d_ny, d_nz; // Number of cells in each dim
Float d_dx, d_dy, d_dz; // Cell length in each dim
Float* d_H; // Cell hydraulic heads
+ Float* d_H_new; // Cell hydraulic heads
Float3* d_V; // Cell fluid velocity
Float3* d_dH; // Cell spatial gradient in heads
Float* d_K; // Cell hydraulic conductivities (anisotropic)
(DIR) diff --git a/src/utility.h b/src/utility.h
t@@ -0,0 +1,21 @@
+#ifndef UTILITY_H_
+#define UTILITY_H_
+
+// MISC. UTILITY FUNCTIONS
+
+
+//Round a / b to nearest higher integer value
+unsigned int iDivUp(unsigned int a, unsigned int b) {
+ return (a % b != 0) ? (a / b + 1) : (a / b);
+}
+
+// Swap two arrays pointers
+void swapFloatArrays(Float* arr1, Float* arr2)
+{
+ Float* tmp = arr1;
+ arr1 = arr2;
+ arr2 = tmp;
+}
+
+
+#endif