tBetter integration of darcy flow into sphere - 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 aa36e4b293f63604e6cd32a36ece1742a156fcc4
(DIR) parent 20e7d61d9aa9111c21e0106a432ac8ac7d4f9ebd
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Mon, 3 Jun 2013 15:32:17 +0200
Better integration of darcy flow into sphere
Diffstat:
M src/CMakeLists.txt | 6 +++---
M src/darcy.cpp | 184 +++++++++++++++++--------------
M src/device.cu | 16 +++++++++++-----
M src/file_io.cpp | 61 ++++++++++++++++++++++++++-----
M src/porousflow.cpp | 9 +++------
M src/sphere.cpp | 7 +++----
M src/sphere.h | 78 ++++++++++++++++++-------------
7 files changed, 217 insertions(+), 144 deletions(-)
---
(DIR) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
t@@ -17,11 +17,11 @@ SET(CUDA_NVCC_FLAGS
# Rule to build executable program
CUDA_ADD_EXECUTABLE(../sphere
- main.cpp file_io.cpp sphere.cpp device.cu utility.cu)
+ main.cpp file_io.cpp sphere.cpp device.cu utility.cu darcy.cpp)
CUDA_ADD_EXECUTABLE(../porosity
- porosity.cpp file_io.cpp sphere.cpp device.cu utility.cu)
+ porosity.cpp file_io.cpp sphere.cpp device.cu utility.cu darcy.cpp)
CUDA_ADD_EXECUTABLE(../forcechains
- forcechains.cpp file_io.cpp sphere.cpp device.cu utility.cu)
+ forcechains.cpp file_io.cpp sphere.cpp device.cu utility.cu darcy.cpp)
CUDA_ADD_EXECUTABLE(../porousflow
porousflow.cpp darcy.cpp file_io.cpp sphere.cpp device.cu utility.cu)
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -14,8 +14,9 @@
void DEM::initDarcyMem()
{
unsigned int ncells = d_nx*d_ny*d_nz;
- d_P = new Float[ncells]; // hydraulic pressure matrix
- d_dP = new Float3[ncells]; // Cell spatial gradient in hydraulic pressures
+ d_H = 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
d_S = new Float[ncells]; // hydraulic storativity matrix
d_W = new Float[ncells]; // hydraulic recharge
t@@ -24,8 +25,9 @@ void DEM::initDarcyMem()
// Free memory
void DEM::freeDarcyMem()
{
- free(d_P);
- free(d_dP);
+ free(d_H);
+ free(d_V);
+ free(d_dH);
free(d_K);
free(d_S);
free(d_W);
t@@ -47,9 +49,18 @@ void DEM::initDarcyVals()
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
- d_P[idx(ix,iy,iz)] = 1.0;
+
+ // Initial hydraulic head [m]
+ //d_H[idx(ix,iy,iz)] = 1.0;
+ // read from input binary
+
+ // Hydraulic conductivity [m/s]
d_K[idx(ix,iy,iz)] = 1.5;
+
+ // Hydraulic storativity [-]
d_S[idx(ix,iy,iz)] = 7.5e-3;
+
+ // Hydraulic recharge [Pa/s]
d_W[idx(ix,iy,iz)] = 0.0;
}
}
t@@ -75,43 +86,36 @@ Float DEM::minVal3dArr(Float* arr)
// Find the spatial gradient in pressures per cell
void DEM::findDarcyGradients()
{
-
- std::cout << "dx,dy,dz: "
- << d_dx << ","
- << d_dy << ","
- << d_dz << std::endl;
+ // Cell sizes squared
const Float dx2 = d_dx*d_dx;
const Float dy2 = d_dy*d_dy;
const Float dz2 = d_dz*d_dz;
- std::cout << "dx2,dy2,dz2: "
- << dx2 << ","
- << dy2 << ","
- << dz2 << std::endl;
- Float localP2;
- unsigned int ix, iy, iz;
+ 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) {
- localP2 = 2.0*d_P[idx(ix,iy,iz)];
+ cellidx = idx(ix,iy,iz);
+
+ H = d_H[cellidx]; // cell hydraulic pressure
- d_dP[idx(ix,iy,iz)].x
- = (d_P[idx(ix+1,iy,iz)] - localP2
- + d_P[idx(ix-1,iy,iz)])/dx2;
+ d_dH[cellidx].x
+ = (d_H[idx(ix+1,iy,iz)] - 2.0*H + d_H[idx(ix-1,iy,iz)])/dx2;
- d_dP[idx(ix,iy,iz)].y
- = (d_P[idx(ix,iy+1,iz)] - localP2
- + d_P[idx(ix,iy-1,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_dP[idx(ix,iy,iz)].z
- = (d_P[idx(ix,iy,iz+1)] - localP2
- + d_P[idx(ix,iy,iz-1)])/dz2;
+ 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@@ -126,77 +130,60 @@ void DEM::setDarcyBCNeumannZero()
// x-y plane at z=0 and z=d_dz-1
for (ix=0; ix<d_nx; ++ix) {
for (iy=0; iy<d_ny; ++iy) {
- d_dP[idx(ix,iy, 0)] = z3;
- d_dP[idx(ix,iy,nz)] = z3;
+ d_dH[idx(ix,iy, 0)] = z3;
+ d_dH[idx(ix,iy,nz)] = z3;
}
}
// x-z plane at y=0 and y=d_dy-1
for (ix=0; ix<d_nx; ++ix) {
for (iz=0; iz<d_nz; ++iz) {
- d_dP[idx(ix, 0,iz)] = z3;
- d_dP[idx(ix,ny,iz)] = z3;
+ d_dH[idx(ix, 0,iz)] = z3;
+ d_dH[idx(ix,ny,iz)] = z3;
}
}
// y-z plane at x=0 and x=d_dx-1
for (iy=0; iy<d_ny; ++iy) {
for (iz=0; iz<d_nz; ++iz) {
- d_dP[idx( 0,iy,iz)] = z3;
- d_dP[idx(nx,iy,iz)] = z3;
+ d_dH[idx( 0,iy,iz)] = z3;
+ d_dH[idx(nx,iy,iz)] = z3;
}
}
}
-void DEM::explDarcyStep(const Float dt)
+// Perform an explicit step.
+// Boundary conditions are fixed gradient values (Neumann)
+void DEM::explDarcyStep()
{
- // Find spatial gradients in all cells
- //findDarcyGradients();
- //printDarcyArray3(stdout, d_dP, "d_dP, after findDarcyGradients");
-
- // Set boundary conditions
- //setDarcyBCNeumannZero();
- //printDarcyArray3(stdout, d_dP, "d_dP, after setDarcyBCNeumannZero");
-
// Cell dims squared
const Float dx2 = d_dx*d_dx;
const Float dy2 = d_dy*d_dy;
const Float dz2 = d_dz*d_dz;
- std::cout << "dx2,dy2,dz2: "
- << dx2 << ","
- << dy2 << ","
- << dz2 << std::endl;
+
+ // Find cell gradients
+ findDarcyGradients();
+ setDarcyBCNeumannZero();
// Explicit 3D finite difference scheme
// new = old + gradient*timestep
unsigned int ix, iy, iz, cellidx;
- Float K, P;
- for (ix=1; ix<d_nx-1; ++ix) {
- for (iy=1; iy<d_ny-1; ++iy) {
- for (iz=1; iz<d_nz-1; ++iz) {
+ Float K, H;
+ Float dt = time.dt;
+ for (ix=0; ix<d_nx; ++ix) {
+ for (iy=0; iy<d_ny; ++iy) {
+ for (iz=0; iz<d_nz; ++iz) {
cellidx = idx(ix,iy,iz);
- /*d_P[cellidx]
- -= (d_W[cellidx]*dt
- + d_K[cellidx]*d_dP[cellidx].x/d_dx
- + d_K[cellidx]*d_dP[cellidx].y/d_dy
- + d_K[cellidx]*d_dP[cellidx].z/d_dz) / d_S[cellidx];*/
-
K = d_K[cellidx]; // cell hydraulic conductivity
- P = d_P[cellidx]; // cell hydraulic pressure
+ H = d_H[cellidx]; // cell hydraulic pressure
- d_P[cellidx]
+ d_H[cellidx]
+= d_W[cellidx]*dt // cell recharge
+ K*dt * // diffusivity term
- (
- (d_P[idx(ix+1,iy,iz)] - 2.0*P + d_P[idx(ix-1,iy,iz)])/dx2 +
- (d_P[idx(ix,iy+1,iz)] - 2.0*P + d_P[idx(ix,iy-1,iz)])/dy2 +
- (d_P[idx(ix,iy,iz+1)] - 2.0*P + d_P[idx(ix,iy,iz-1)])/dz2
- );
-
-
+ (d_dH[cellidx].x + d_dH[cellidx].y + d_dH[cellidx].z);
}
}
}
t@@ -249,23 +236,61 @@ void DEM::printDarcyArray3(FILE* stream, Float3* arr, std::string desc)
printDarcyArray3(stream, arr);
}
+// Find cell velocity
+void DEM::findDarcyVelocities()
+{
+ // Flux [m/s]: q = -k/nu * dH
+ // Pore velocity [m/s]: v = q/n
+ Float3 q, v, dH;
+
+ // Dynamic viscosity
+ Float nu = params.nu;
+
+ // Porosity [-]: n
+ Float n = 0.5; // CHANGE THIS
+
+ 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) {
+
+ cellidx = idx(ix,iy,iz);
+ dH = d_dH[cellidx];
+
+ // Calculate flux
+ q.x = -d_K[cellidx]/nu * dH.x;
+ q.y = -d_K[cellidx]/nu * dH.y;
+ q.z = -d_K[cellidx]/nu * dH.z;
+
+ // Calculate velocity
+ v.x = q.x/n;
+ v.y = q.y/n;
+ v.z = q.z/n;
-// Find hydraulic conductivities for each cell by finding the particle contents
-//
+ d_V[cellidx] = v;
+ }
+ }
+ }
+}
// Solve Darcy flow on a regular, cubic grid
-void DEM::startDarcy(
- const Float cellsizemultiplier)
+void DEM::initDarcy(const Float cellsizemultiplier)
{
+ if (params.nu <= 0.0) {
+ std::cerr << "Error in initDarcy. The dymamic viscosity (params.nu), "
+ << "should be larger than 0.0, but is " << params.nu << std::endl;
+ 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
- Float d_dx = grid.L[0]/d_nx;
- Float d_dy = grid.L[1]/d_ny;
- Float d_dz = grid.L[2]/d_nz;
+ d_dx = grid.L[0]/d_nx;
+ d_dy = grid.L[1]/d_ny;
+ d_dz = grid.L[2]/d_nz;
if (verbose == 1) {
std::cout << " - Fluid grid dimensions: "
t@@ -280,20 +305,11 @@ void DEM::startDarcy(
initDarcyMem();
initDarcyVals();
+}
- // Temporal loop
- //while(time.current <= time.total) {
-
- explDarcyStep(time.dt);
- time.current += time.dt;
- //}
-
-
- printDarcyArray(stdout, d_P, "d_P");
- //printDarcyArray3(stdout, d_dP, "d_dP");
- //printDarcyArray(stdout, d_K, "d_K");
- //printDarcyArray(stdout, d_S, "d_S");
- //printDarcyArray(stdout, d_W, "d_W");
+void DEM::endDarcy()
+{
+ printDarcyArray(stdout, d_H, "d_H");
freeDarcyMem();
}
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -448,7 +448,7 @@ __host__ void DEM::transferToGlobalDeviceMemory()
}
// Fluid arrays
- if (params.nu > 0.0) {
+ if (params.nu > 0.0 && darcy == 0) {
#ifdef LBM_GPU
cudaMemcpy( dev_f, f,
sizeof(Float)*grid.num[0]*grid.num[1]*grid.num[2]*19,
t@@ -529,7 +529,7 @@ __host__ void DEM::transferFromGlobalDeviceMemory()
// Fluid arrays
#ifdef LBM_GPU
- if (params.nu > 0.0) {
+ if (params.nu > 0.0 && darcy == 0) {
cudaMemcpy( f, dev_f,
sizeof(Float)*grid.num[0]*grid.num[1]*grid.num[2]*19,
cudaMemcpyDeviceToHost);
t@@ -605,7 +605,7 @@ __host__ void DEM::startTime()
<< dimBlock.x << "*" << dimBlock.y << "*" << dimBlock.z << "\n"
<< " - Shared memory required per block: " << smemSize << " bytes"
<< endl;
- if (params.nu > 0.0) {
+ if (params.nu > 0.0 && darcy == 0) {
cout << " - Blocks per fluid grid: "
<< dimGridFluid.x << "*" << dimGridFluid.y << "*" <<
dimGridFluid.z << "\n"
t@@ -630,7 +630,7 @@ __host__ void DEM::startTime()
fclose(fp);
// Initialize fluid distribution array
- if (params.nu > 0.0) {
+ if (params.nu > 0.0 && darcy == 0) {
#ifdef LBM_GPU
initFluid<<< dimGridFluid, dimBlockFluid >>>(dev_v_rho, dev_f);
cudaThreadSynchronize();
t@@ -825,7 +825,7 @@ __host__ void DEM::startTime()
}
// Process fluid and particle interaction in each cell
- if (params.nu > 0.0 && grid.periodic == 1) {
+ if (params.nu > 0.0 && darcy == 0 && grid.periodic == 1) {
#ifdef LBM_GPU
if (PROFILING == 1)
startTimer(&kernel_tic);
t@@ -856,6 +856,9 @@ __host__ void DEM::startTime()
}
+ // Solve darcy flow through grid
+ if (darcy == 1)
+ explDarcyStep();
// Update particle kinematics
t@@ -1064,5 +1067,8 @@ __host__ void DEM::startTime()
delete[] k.distmod;
delete[] k.delta_t;
+ if (darcy == 1)
+ endDarcy();
+
} /* EOF */
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -238,11 +238,20 @@ void DEM::readbin(const char *target)
// Read fluid parameters
ifs.read(as_bytes(params.nu), sizeof(params.nu));
- 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]];
- if (params.nu > 0.0) {
- unsigned int x, y, z;
+ unsigned int x, y, z;
+
+ if (verbose == 1)
+ cout << "Done\n";
+
+ if (params.nu > 0.0 && darcy == 0) { // Lattice-Boltzmann flow
+
+ if (verbose == 1)
+ cout << " - Reading LBM values:\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];
+ v_rho = new Float4[grid.num[0]*grid.num[1]*grid.num[2]];
+
for (z = 0; z<grid.num[2]; ++z) {
for (y = 0; y<grid.num[1]; ++y) {
for (x = 0; x<grid.num[0]; ++x) {
t@@ -254,15 +263,37 @@ void DEM::readbin(const char *target)
}
}
}
- }
+ if (verbose == 1)
+ cout << "Done" << std::endl;
+
+ } else if (params.nu > 0.0 && darcy == 1) { // Darcy flow
+
+ initDarcy();
+
+ if (verbose == 1)
+ cout << " - Reading Darcy values:\t\t\t ";
+
+ for (z = 0; z<grid.num[2]; ++z) {
+ for (y = 0; y<grid.num[1]; ++y) {
+ for (x = 0; x<grid.num[0]; ++x) {
+ i = idx(x,y,z);
+ ifs.read(as_bytes(d_V[i].x), sizeof(Float));
+ ifs.read(as_bytes(d_V[i].y), sizeof(Float));
+ ifs.read(as_bytes(d_V[i].z), sizeof(Float));
+ ifs.read(as_bytes(d_H[i]), sizeof(Float));
+ }
+ }
+ }
+
+ if (verbose == 1)
+ cout << "Done" << std::endl;
+ }
// Close file if it is still open
if (ifs.is_open())
ifs.close();
- if (verbose == 1)
- cout << "Done\n";
}
t@@ -417,7 +448,7 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(params.nu), sizeof(params.nu));
unsigned int x, y, z;
- if (params.nu > 0.0) {
+ if (params.nu > 0.0 && darcy == 0) { // Lattice Boltzmann flow
for (z = 0; z<grid.num[2]; ++z) {
for (y = 0; y<grid.num[1]; ++y) {
for (x = 0; x<grid.num[0]; ++x) {
t@@ -429,6 +460,18 @@ void DEM::writebin(const char *target)
}
}
}
+ } else if (params.nu > 0.0 && darcy == 1) { // Darcy flow
+ for (z = 0; z<d_dz; ++z) {
+ for (y = 0; y<d_dy; ++y) {
+ for (x = 0; x<d_dx; ++x) {
+ i = idx(x,y,z);
+ ofs.write(as_bytes(d_V[i].x), sizeof(Float));
+ ofs.write(as_bytes(d_V[i].y), sizeof(Float));
+ ofs.write(as_bytes(d_V[i].z), sizeof(Float));
+ ofs.write(as_bytes(d_H[i]), sizeof(Float));
+ }
+ }
+ }
}
// Close file if it is still open
(DIR) diff --git a/src/porousflow.cpp b/src/porousflow.cpp
t@@ -60,13 +60,10 @@ int main(const int argc, const char *argv[])
std::endl;
// Create DEM class, read data from input binary,
- // do not check values, do not init cuda, do not transfer const
- // mem
- DEM dem(argvi, verbose, 0, dry, 0, 0);
-
- // Start iterating through time
- dem.startDarcy();
+ // check values, init cuda, transfer const mem, solve Darcy flow
+ DEM dem(argvi, verbose, 1, dry, 1, 1, 1);
+ dem.startTime();
}
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -18,8 +18,9 @@ DEM::DEM(const std::string inputbin,
const int checkVals,
const int dry,
const int initCuda,
- const int transferConstMem)
-: verbose(verbosity)
+ const int transferConstMem,
+ const int darcyflow)
+: verbose(verbosity), darcy(darcyflow)
{
using std::cout;
using std::cerr;
t@@ -65,8 +66,6 @@ DEM::DEM(const std::string inputbin,
// Transfer data from host to gpu device memory
transferToGlobalDeviceMemory();
}
-
-
}
// Destructor: Liberates dynamically allocated host memory
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -144,16 +144,54 @@ class DEM {
Float4 *v_rho; // Fluid velocity v (xyz), and pressure rho (w)
Float4 *dev_v_rho; // Device equivalent
- // Darcy-flow values
+ //// Darcy-flow
+ int darcy; // 0: no, 1: yes
+
+ // Darcy values
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_P; // Cell hydraulic pressures
- Float3* d_dP; // Cell spatial gradient in pressures
+ Float* d_H; // Cell hydraulic heads
+ Float3* d_dH; // Cell spatial gradient in heads
Float* d_K; // Cell hydraulic conductivities (anisotropic)
Float* d_S; // Cell hydraulic storativity
Float* d_W; // Cell hydraulic recharge
- Float mu; // Fluid viscosity
+ Float3* d_V; // Cell fluid velocity
+
+ // Darcy functions
+
+ // Memory allocation
+ void initDarcyMem();
+ void freeDarcyMem();
+
+ // Set some values for the Darcy parameters
+ void initDarcyVals();
+
+ // Finds central difference gradients
+ void findDarcyGradients();
+
+ // Set gradient to zero at grid edges
+ void setDarcyBCNeumannZero();
+ // Find darcy flow velocities from specific flux (q)
+ void findDarcyVelocities();
+
+ // Get linear (1D) index from 3D coordinate
+ unsigned int idx(
+ const unsigned int x,
+ const unsigned int y,
+ const unsigned int z);
+
+ // Get minimum value in 1D array
+ Float minVal3dArr(Float* arr);
+
+ // Initialize Darcy values and arrays
+ void initDarcy(const Float cellsizemultiplier = 1.0);
+
+ // Clean up Darcy arrays
+ void endDarcy();
+
+ // Perform a single time step, explicit integration
+ void explDarcyStep();
public:
t@@ -165,7 +203,8 @@ class DEM {
const int checkVals = 1,
const int dry = 0,
const int initCuda = 1,
- const int transferConstMem = 1);
+ const int transferConstMem = 1,
+ const int darcyflow = 0);
// Destructor
~DEM(void);
t@@ -223,40 +262,13 @@ class DEM {
///// Darcy flow functions
- // Memory allocation and destruction
- void initDarcyMem();
- void freeDarcyMem();
-
- // Set some values for the Darcy parameters
- void initDarcyVals();
-
- // Get linear (1D) index from 3D coordinate
- unsigned int idx(
- const unsigned int x,
- const unsigned int y,
- const unsigned int z);
-
- // Get minimum value in 1D array
- Float minVal3dArr(Float* arr);
-
- // Finds central difference gradients
- void findDarcyGradients();
-
- // Set gradient to zero at grid edges
- void setDarcyBCNeumannZero();
-
- // Perform a single time step, explicit integration
- void explDarcyStep(const Float dt);
-
- // Calculate Darcy fluid flow through material
- void startDarcy(
- const Float cellsizemultiplier = 1.0);
// Print Darcy arrays to file stream
void printDarcyArray(FILE* stream, Float* arr);
void printDarcyArray(FILE* stream, Float* arr, std::string desc);
void printDarcyArray3(FILE* stream, Float3* arr);
void printDarcyArray3(FILE* stream, Float3* arr, std::string desc);
+
};
#endif