tDarcy kernels disabled, IO verified - 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 7fff983d1864e9f0ca6fb35e4c0ac01dfce83e57
(DIR) parent 503e3db3d4c4ee215350fcb48da7956ab08817d5
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 9 Oct 2013 13:43:56 +0200
Darcy kernels disabled, IO verified
Diffstat:
M src/darcy.cpp | 83 ++++++++++++++++---------------
M src/darcy.cuh | 8 ++++----
M src/device.cu | 28 ++++++++++++++--------------
M src/file_io.cpp | 7 ++++---
M src/sphere.cpp | 5 +++++
M src/sphere.h | 25 ++++++-------------------
6 files changed, 76 insertions(+), 80 deletions(-)
---
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -14,10 +14,16 @@
#define PERIODIC_XY
// Initialize memory
-void DEM::initDarcyMem()
+void DEM::initDarcyMem(const Float cellsizemultiplier)
{
+ // Number of cells
+ d_nx = floor(grid.num[0]*cellsizemultiplier);
+ d_ny = floor(grid.num[1]*cellsizemultiplier);
+ d_nz = floor(grid.num[2]*cellsizemultiplier);
+
//unsigned int ncells = d_nx*d_ny*d_nz; // without ghost nodes
unsigned int ncells = (d_nx+2)*(d_ny+2)*(d_nz+2); // with ghost nodes
+
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
t@@ -45,9 +51,9 @@ void DEM::freeDarcyMem()
// 3D index to 1D index
unsigned int DEM::idx(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z)
+ const int x,
+ const int y,
+ const int z)
{
// without ghost nodes
//return x + d_nx*y + d_nx*d_ny*z;
t@@ -66,10 +72,12 @@ void DEM::initDarcyVals()
// Density of the fluid [kg/m^3]
const Float rho = 1000.0;
- unsigned int ix, iy, iz, cellidx;
- for (ix=0; ix<d_nx; ++ix) {
- for (iy=0; iy<d_ny; ++iy) {
- for (iz=0; iz<d_nz; ++iz) {
+ int ix, iy, iz, cellidx;
+
+ // Set values for all cells, including ghost nodes
+ for (ix=-1; ix<=d_nx; ++ix) {
+ for (iy=-1; iy<=d_ny; ++iy) {
+ for (iz=-1; iz<=d_nz; ++iz) {
cellidx = idx(ix,iy,iz);
t@@ -119,7 +127,7 @@ void DEM::copyDarcyVals(unsigned int read, unsigned int write)
// The edge (diagonal) cells are not written since they are not read
void DEM::setDarcyGhostNodes()
{
- unsigned int ix, iy, iz;
+ int ix, iy, iz;
// The x-normal plane
for (iy=0; iy<d_ny; ++iy) {
t@@ -188,7 +196,7 @@ void DEM::findDarcyTransmissivities()
// Grain size factor for Kozeny-Carman relationship
Float d_factor = r_bar2*r_bar2/180.0;
- unsigned int ix, iy, iz, cellidx;
+ int ix, iy, iz, cellidx;
Float K, k;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
t@@ -225,10 +233,10 @@ void DEM::findDarcyTransmissivities()
void DEM::setDarcyBCNeumannZero()
{
Float3 z3 = MAKE_FLOAT3(0.0, 0.0, 0.0);
- unsigned int ix, iy, iz;
- unsigned int nx = d_nx-1;
- unsigned int ny = d_ny-1;
- unsigned int nz = d_nz-1;
+ int ix, iy, iz;
+ int nx = d_nx-1;
+ int ny = d_ny-1;
+ int nz = d_nz-1;
// I don't care that the values at four edges are written twice
t@@ -248,7 +256,7 @@ void DEM::setDarcyBCNeumannZero()
}
}
- // y-z plane at x=0 and x=d_dx-1
+ // y-z plane at x=0 and x=d_nx-1
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
d_dH[idx( 0,iy,iz)] = z3;
t@@ -270,7 +278,7 @@ void DEM::findDarcyGradients()
const Float dz2 = 2.0*d_dz;
//Float H;
- unsigned int ix, iy, iz, cellidx;
+ int ix, iy, iz, cellidx;
// Without ghost-nodes
/*for (ix=1; ix<d_nx-1; ++ix) {
t@@ -367,7 +375,7 @@ void DEM::explDarcyStep()
// Explicit 3D finite difference scheme
// new = old + production*timestep + gradient*timestep
- unsigned int ix, iy, iz, cellidx;
+ int ix, iy, iz, cellidx;
Float K, H, deltaH;
Float Tx, Ty, Tz, S;
//Float Tx_n, Tx_p, Ty_n, Ty_p, Tz_n, Tz_p;
t@@ -493,7 +501,7 @@ void DEM::explDarcyStep()
// Print array values to file stream (stdout, stderr, other file)
void DEM::printDarcyArray(FILE* stream, Float* arr)
{
- unsigned int x, y, z;
+ int x, y, z;
for (z=0; z<d_nz; z++) {
for (y=0; y<d_ny; y++) {
for (x=0; x<d_nx; x++) {
t@@ -515,7 +523,7 @@ void DEM::printDarcyArray(FILE* stream, Float* arr, std::string desc)
// Print array values to file stream (stdout, stderr, other file)
void DEM::printDarcyArray(FILE* stream, Float3* arr)
{
- unsigned int x, y, z;
+ int x, y, z;
for (z=0; z<d_nz; z++) {
for (y=0; y<d_ny; y++) {
for (x=0; x<d_nx; x++) {
t@@ -552,7 +560,7 @@ void DEM::findDarcyVelocities()
// Find cell gradients
findDarcyGradients();
- unsigned int ix, iy, iz, cellidx;
+ int ix, iy, iz, cellidx;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
t@@ -582,9 +590,9 @@ void DEM::findDarcyVelocities()
// Return the lower corner coordinates of a cell
Float3 DEM::cellMinBoundaryDarcy(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z)
+ const int x,
+ const int y,
+ const int z)
{
const Float3 x_min = {x*d_dx, y*d_dy, z*d_dz};
return x_min;
t@@ -592,9 +600,9 @@ Float3 DEM::cellMinBoundaryDarcy(
// Return the upper corner coordinates of a cell
Float3 DEM::cellMaxBoundaryDarcy(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z)
+ const int x,
+ const int y,
+ const int z)
{
const Float3 x_max = {(x+1)*d_dx, (y+1)*d_dy, (z+1)*d_dz};
return x_max;
t@@ -609,9 +617,9 @@ Float DEM::cellVolumeDarcy()
// Find the porosity of a target cell
Float DEM::cellPorosity(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z)
+ const int x,
+ const int y,
+ const int z)
{
const Float3 x_min = cellMinBoundaryDarcy(x,y,z);
const Float3 x_max = cellMaxBoundaryDarcy(x,y,z);
t@@ -640,7 +648,7 @@ Float DEM::cellPorosity(
// Calculate the porosity for each cell
void DEM::findPorosities()
{
- unsigned int ix, iy, iz, cellidx;
+ int ix, iy, iz, cellidx;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
t@@ -700,7 +708,7 @@ void DEM::fluidDragDarcy()
Float DEM::getTmax()
{
Float max = -1.0e13; // initialize with a small number
- unsigned int ix,iy,iz;
+ int ix,iy,iz;
Float3 val;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
t@@ -721,7 +729,7 @@ Float DEM::getTmax()
Float DEM::getSsmin()
{
Float min = 1.0e13; // initialize with a small number
- unsigned int ix,iy,iz;
+ int ix,iy,iz;
Float val;
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
t@@ -748,7 +756,7 @@ void DEM::checkDarcyTimestep()
* (time.dt/(d_dx*d_dx) + time.dt/(d_dy*d_dy) + time.dt/(d_dz*d_dz));
if (value > 0.5) {
- std::cerr << "Error! The explicit darcy solution will be unstable.\n"
+ std::cerr << "Error! The explicit Darcy solution will be unstable.\n"
<< "This happens due to a combination of the following:\n"
<< " - The transmissivity T (i.e. hydraulic conductivity, K)"
<< " is too large (" << T_max << ")\n"
t@@ -764,7 +772,7 @@ void DEM::checkDarcyTimestep()
}
// Initialize darcy arrays, their values, and check the time step length
-void DEM::initDarcy(const Float cellsizemultiplier)
+void DEM::initDarcy()
{
if (params.nu <= 0.0) {
std::cerr << "Error in initDarcy. The dymamic viscosity (params.nu), "
t@@ -772,11 +780,6 @@ void DEM::initDarcy(const Float cellsizemultiplier)
exit(1);
}
- // Number of cells
- d_nx = floor(grid.num[0]*cellsizemultiplier);
- d_ny = floor(grid.num[1]*cellsizemultiplier);
- d_nz = floor(grid.num[2]*cellsizemultiplier);
-
// Cell size
d_dx = grid.L[0]/d_nx;
d_dy = grid.L[1]/d_ny;
t@@ -793,7 +796,7 @@ void DEM::initDarcy(const Float cellsizemultiplier)
<< d_dz << std::endl;
}
- initDarcyMem();
+ //initDarcyMem(); // done in readbin
initDarcyVals();
findDarcyTransmissivities();
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -32,9 +32,9 @@ void DEM::initDarcyMemDev(void)
cudaMalloc((void**)&dev_d_dH, memSizeF*3); // hydraulic pressure gradient
cudaMalloc((void**)&dev_d_K, memSizeF); // hydraulic conductivity
cudaMalloc((void**)&dev_d_T, memSizeF*3); // hydraulic transmissivity
- cudaMalloc((void**)&dev_d_Ss, memSizeF); // hydraulic storativi
+ cudaMalloc((void**)&dev_d_Ss, memSizeF); // hydraulic storativi
cudaMalloc((void**)&dev_d_W, memSizeF); // hydraulic recharge
- cudaMalloc((void**)&dev_d_phi, memSizeF); // cell porosity
+ cudaMalloc((void**)&dev_d_phi, memSizeF); // cell porosity
checkForCudaErrors("End of initDarcyMemDev");
}
t@@ -354,8 +354,8 @@ __global__ void findDarcyTransmissivitiesDev(
__syncthreads();
// Save hydraulic conductivity [m/s]
- const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
- //K = 0.5;
+ //const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
+ const Float K = 0.5;
dev_d_K[cellidx] = K;
// Hydraulic transmissivity [m2/s]
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -307,12 +307,14 @@ __host__ void DEM::allocateGlobalDeviceMemory(void)
// dev_walls_force_partial allocated later
// Fluid arrays
+#ifdef LBM_GPU
cudaMalloc((void**)&dev_f,
sizeof(Float)*grid.num[0]*grid.num[1]*grid.num[2]*19);
cudaMalloc((void**)&dev_f_new,
sizeof(Float)*grid.num[0]*grid.num[1]*grid.num[2]*19);
cudaMalloc((void**)&dev_v_rho,
sizeof(Float4)*grid.num[0]*grid.num[1]*grid.num[2]);
+#endif
checkForCudaErrors("End of allocateGlobalDeviceMemory");
if (verbose == 1)
t@@ -366,9 +368,11 @@ __host__ void DEM::freeGlobalDeviceMemory()
cudaFree(dev_walls_vel0);
// Fluid arrays
+#ifdef LBM_GPU
cudaFree(dev_f);
cudaFree(dev_f_new);
cudaFree(dev_v_rho);
+#endif
checkForCudaErrors("During cudaFree calls");
t@@ -460,10 +464,12 @@ __host__ void DEM::transferToGlobalDeviceMemory(int statusmsg)
sizeof(Float4)*grid.num[0]*grid.num[1]*grid.num[2],
cudaMemcpyHostToDevice);
#endif
- } //else {
- //transferDarcyToGlobalDeviceMemory(1);
- //}
- // Darcy arrays aren't ready yet
+ } else if (darcy == 1) {
+ transferDarcyToGlobalDeviceMemory(1);
+ } else {
+ std::cerr << "darcy value not understood ("
+ << darcy << ")" << std::endl;
+ }
}
checkForCudaErrors("End of transferToGlobalDeviceMemory");
t@@ -652,18 +658,10 @@ __host__ void DEM::startTime()
#endif
} else if (darcy == 1) {
#ifdef DARCY_GPU
- const Float cellsizemultiplier = 1.0;
- initDarcy(cellsizemultiplier);
- initDarcyMemDev();
- transferDarcyToGlobalDeviceMemory(1);
-
// Representative grain radius
const Float r_bar2 = meanRadius()*2.0;
// Grain size factor for Kozeny-Carman relationship
d_factor = r_bar2*r_bar2/180.0;
-#else
- const Float cellsizemultiplier = 1.0;
- initDarcy(cellsizemultiplier);
#endif
} else {
std::cerr << "Error, darcy value (" << darcy
t@@ -899,6 +897,7 @@ __host__ void DEM::startTime()
if (params.nu > 0.0 && darcy == 1) {
#ifdef DARCY_GPU
+ /*
checkForCudaErrors("Before findPorositiesDev", iter);
// Find cell porosities
if (PROFILING == 1)
t@@ -990,6 +989,7 @@ __host__ void DEM::startTime()
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
&t_findDarcyVelocitiesDev);
checkForCudaErrors("Post findDarcyVelocitiesDev", iter);
+ */
#else
// Copy device data to host memory
t@@ -1239,12 +1239,12 @@ __host__ void DEM::startTime()
delete[] k.delta_t;
#ifndef DARCY_GPU
- if (darcy == 1) {
+ if (darcy == 1 && params.nu > 0.0)
endDarcyDev();
endDarcy();
}
#else
- if (darcy == 1)
+ if (darcy == 1 && params.nu > 0.0)
endDarcy();
#endif
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -248,8 +248,8 @@ void DEM::readbin(const char *target)
if (verbose == 1)
cout << " - Reading LBM values:\t\t\t\t ";
- f = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
- f_new = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
+ //f = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
+ //f_new = new Float[grid.num[0]*grid.num[1]*grid.num[2]*19];
v_rho = new Float4[grid.num[0]*grid.num[1]*grid.num[2]];
for (z = 0; z<grid.num[2]; ++z) {
t@@ -269,7 +269,8 @@ void DEM::readbin(const char *target)
} else if (params.nu > 0.0 && darcy == 1) { // Darcy flow
- initDarcy();
+ const Float cellsizemultiplier = 1.0;
+ initDarcyMem(cellsizemultiplier);
if (verbose == 1)
cout << " - Reading Darcy values:\t\t\t ";
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -59,6 +59,11 @@ DEM::DEM(const std::string inputbin,
transferToConstantDeviceMemory();
}
+ if (params.nu > 0.0 && darcy == 1) {
+ initDarcy();
+ initDarcyMemDev();
+ }
+
// Allocate device memory for particle variables,
// tied to previously declared pointers in structures
allocateGlobalDeviceMemory();
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -174,7 +174,7 @@ class DEM {
//// Darcy functions
// Memory allocation
- void initDarcyMem();
+ void initDarcyMem(const Float cellsizemultiplier = 1.0);
void freeDarcyMem();
// Set some values for the Darcy parameters
t@@ -201,16 +201,10 @@ class DEM {
const Float3 min, const Float3 max);
// Return the lower corner coordinates of a cell
- Float3 cellMinBoundaryDarcy(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z);
+ Float3 cellMinBoundaryDarcy(const int x, const int y, const int z);
// Return the upper corner coordinates of a cell
- Float3 cellMaxBoundaryDarcy(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z);
+ Float3 cellMaxBoundaryDarcy(const int x, const int y, const int z);
// Returns the cell volume
Float cellVolumeDarcy();
t@@ -219,10 +213,7 @@ class DEM {
void fluidDragDarcy();
// Find porosity of cell
- Float cellPorosity(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z);
+ Float cellPorosity(const int x, const int y, const int z);
// Find and save all cell porosities
void findPorosities();
t@@ -234,14 +225,10 @@ class DEM {
Float meanRadius();
// Get linear (1D) index from 3D coordinate
- unsigned int idx(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z);
+ unsigned int idx(const int x, const int y, const int z);
// Initialize Darcy values and arrays
- void initDarcy(const Float cellsizemultiplier = 1.0);
- void initDarcyDev(const Float cellsizemultiplier = 1.0);
+ void initDarcy();
// Clean up Darcy arrays
void endDarcy();