tnumerous changes, added sphere based porosity est method - 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 8e970aa86be88c5c45b37b1bcb6e2b2a2edfde02
(DIR) parent 4cc89f81f1466e8ad89bd8e9d2dd2bca89eafa20
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 16 Oct 2013 10:50:47 +0200
numerous changes, added sphere based porosity est method
Diffstat:
M src/darcy.cpp | 75 +++++++++++++++++--------------
M src/darcy.cuh | 337 +++++++++++++++++++------------
M src/device.cu | 43 +++++++++++++++++++------------
M src/file_io.cpp | 9 +++++----
M src/sphere.cpp | 4 ++++
M src/sphere.h | 3 +++
6 files changed, 288 insertions(+), 183 deletions(-)
---
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -81,17 +81,9 @@ void DEM::initDarcyVals()
cellidx = idx(ix,iy,iz);
- // Initial hydraulic head [m]
- //d_H[cellidx] = 1.0;
- // read from input binary
-
- // Hydraulic permeability [m^2]
- //d_K[cellidx] = k*rho*-params.g[2]/params.nu;
- d_K[cellidx] = 0.5;
-
// Hydraulic storativity [-]
- d_Ss[cellidx] = 8.0e-3;
//d_Ss[cellidx] = 1.0;
+ d_Ss[cellidx] = 8.0e-3;
// Hydraulic recharge [s^-1]
d_W[cellidx] = 0.0;
t@@ -100,12 +92,12 @@ void DEM::initDarcyVals()
}
// Extract water from all cells in center
- ix = d_nx/2-1; iy = d_ny/2-1;
+ /*ix = d_nx/2; iy = d_ny/2;
Float cellvolume = d_dx*d_dy*d_dz;
for (iz=0; iz<d_nz; ++iz) {
- //d_W[idx(ix,iy,iz)] = -0.1*cellvolume;
- d_W[idx(ix,iy,iz)] = -1.0;
- }
+ //d_W[idx(ix,iy,iz)] = -1.0e-4/cellvolume;
+ d_W[idx(ix,iy,iz)] = -2.0e-3;
+ }*/
}
t@@ -352,7 +344,6 @@ Float hmean(Float a, Float b) {
}
// Perform an explicit step.
-// Boundary conditions are fixed values (Dirichlet)
void DEM::explDarcyStep()
{
t@@ -398,13 +389,8 @@ void DEM::explDarcyStep()
} else {*/
- // Cell hydraulic conductivity
- K = d_K[cellidx];
-
// Cell hydraulic transmissivities
- Tx = K*d_dx;
- Ty = K*d_dy;
- Tz = K*d_dz;
+ const Float3 T = d_T[cellidx];
// Cell hydraulic head
H = d_H[cellidx];
t@@ -443,14 +429,14 @@ void DEM::explDarcyStep()
* (d_H[idx(ix,iy+1,iz)] - H)/dy2;
*/
- gradx_n = hmean(Tx, d_T[idx(ix-1,iy,iz)].x)
+ gradx_n = hmean(T.x, d_T[idx(ix-1,iy,iz)].x)
* (d_H[idx(ix-1,iy,iz)] - H)/dxdx;
- gradx_p = hmean(Tx, d_T[idx(ix+1,iy,iz)].x)
+ gradx_p = hmean(T.x, d_T[idx(ix+1,iy,iz)].x)
* (d_H[idx(ix+1,iy,iz)] - H)/dxdx;
- grady_n = hmean(Ty, d_T[idx(ix,iy-1,iz)].y)
+ grady_n = hmean(T.y, d_T[idx(ix,iy-1,iz)].y)
* (d_H[idx(ix,iy-1,iz)] - H)/dydy;
- grady_p = hmean(Ty, d_T[idx(ix,iy+1,iz)].y)
+ grady_p = hmean(T.y, d_T[idx(ix,iy+1,iz)].y)
* (d_H[idx(ix,iy+1,iz)] - H)/dydy;
// Neumann (no-flow) boundary condition at +z and -z boundaries
t@@ -458,12 +444,12 @@ void DEM::explDarcyStep()
if (iz == 0)
gradz_n = 0.0;
else
- gradz_n = hmean(Tz, d_T[idx(ix,iy,iz-1)].z)
+ gradz_n = hmean(T.z, d_T[idx(ix,iy,iz-1)].z)
* (d_H[idx(ix,iy,iz-1)] - H)/dzdz;
if (iz == d_nz-1)
gradz_p = 0.0;
else
- gradz_p = hmean(Tz, d_T[idx(ix,iy,iz+1)].z)
+ gradz_p = hmean(T.z, d_T[idx(ix,iy,iz+1)].z)
* (d_H[idx(ix,iy,iz+1)] - H)/dzdz;
/*std::cerr << ix << ',' << iy << ',' << iz << '\t'
t@@ -803,19 +789,40 @@ void DEM::initDarcy()
checkDarcyTimestep();
}
+// Write values in scalar field to file
+void DEM::writeDarcyArray(Float* array, const char* filename)
+{
+ FILE* file;
+ if ((file = fopen(filename,"w"))) {
+ printDarcyArray(file, array);
+ fclose(file);
+ } else {
+ fprintf(stderr, "Error, could not open %s.\n", filename);
+ }
+}
+
+// Write values in vector field to file
+void DEM::writeDarcyArray(Float3* array, const char* filename)
+{
+ FILE* file;
+ if ((file = fopen(filename,"w"))) {
+ printDarcyArray(file, array);
+ fclose(file);
+ } else {
+ fprintf(stderr, "Error, could not open %s.\n", filename);
+ }
+}
+
+
// Print final heads and free memory
void DEM::endDarcy()
{
- /*FILE* Kfile;
- if ((Kfile = fopen("d_K.txt","w"))) {
- printDarcyArray(Kfile, d_K);
- fclose(Kfile);
- } else {
- fprintf(stderr, "Error, could not open d_K.txt\n");
- }*/
- //printDarcyArray(stdout, d_phi, "d_phi");
+ writeDarcyArray(d_phi, "d_phi.txt");
+ writeDarcyArray(d_K, "d_K.txt");
+
//printDarcyArray(stdout, d_K, "d_K");
//printDarcyArray(stdout, d_H, "d_H");
+ //printDarcyArray(stdout, d_H_new, "d_H_new");
//printDarcyArray(stdout, d_V, "d_V");
freeDarcyMem();
}
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -2,7 +2,7 @@
// CUDA implementation of Darcy flow
// Enable line below to perform Darcy flow computations on the GPU, disable for
-// CPU computation
+// Enable GPU computation
#define DARCY_GPU
#include <iostream>
t@@ -244,8 +244,10 @@ __global__ void setDarcyGhostNodesDev(
}
}
-// Find the porosity in each cell
-__global__ void findPorositiesDev(
+// Find the porosity in each cell on the base of a cubic grid, binning particles
+// into the cells containing their centers. This approximation causes
+// non-continuous porosities through time.
+__global__ void findPorositiesCubicDev(
unsigned int* dev_cellStart,
unsigned int* dev_cellEnd,
Float4* dev_x_sorted,
t@@ -280,26 +282,138 @@ __global__ void findPorositiesDev(
// Lowest particle index in cell
const unsigned int startIdx = dev_cellStart[cellID];
- // Highest particle index in cell
- const unsigned int endIdx = dev_cellEnd[cellID];
+ Float phi = 0.99;
- // Iterate over cell particles
- for (unsigned int i = startIdx; i<endIdx; ++i) {
+ // Make sure cell is not empty
+ if (startIdx != 0xffffffff) {
- // Read particle position and radius
- __syncthreads();
- xr = dev_x_sorted[i];
+ // Highest particle index in cell
+ const unsigned int endIdx = dev_cellEnd[cellID];
+
+ // Iterate over cell particles
+ for (unsigned int i = startIdx; i<endIdx; ++i) {
+
+ // Read particle position and radius
+ __syncthreads();
+ xr = dev_x_sorted[i];
- // Subtract particle volume from void volume
- void_volume -= 4.0/3.0*M_PI*xr.w*xr.w*xr.w;
+ // Subtract particle volume from void volume
+ void_volume -= 4.0/3.0*M_PI*xr.w*xr.w*xr.w;
+ }
+
+ // Make sure that the porosity is in the interval ]0.0;1.0[
+ phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
+ }
+
+ // Save porosity
+ __syncthreads();
+ dev_d_phi[idx(x,y,z)] = phi;
+ }
+}
+
+
+// Find the porosity in each cell on the base of a sphere, centered at the cell
+// center. This approximation is continuous through time and generally
+// preferable to findPorositiesCubicDev, although it's slower.
+__global__ void findPorositiesSphericalDev(
+ unsigned int* dev_cellStart,
+ unsigned int* dev_cellEnd,
+ Float4* dev_x_sorted,
+ Float* dev_d_phi)
+{
+ // 3D thread index
+ const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
+ const unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
+ const unsigned int z = blockDim.z * blockIdx.z + threadIdx.z;
+
+ // Grid dimensions
+ const unsigned int nx = devC_grid.num[0];
+ const unsigned int ny = devC_grid.num[1];
+ const unsigned int nz = devC_grid.num[2];
+
+ // Cell dimensions
+ const Float dx = devC_grid.L[0]/nx;
+ const Float dy = devC_grid.L[1]/ny;
+ const Float dz = devC_grid.L[2]/nz;
+
+ // Cell sphere radius
+ const Float R = fmin(dx, fmin(dy,dz)) * 0.5;
+ const Float cell_volume = 4.0/3.0*M_PI*R*R*R;
+
+ Float void_volume = cell_volume;
+ Float4 xr; // particle pos. and radius
+
+ // check that we are not outside the fluid grid
+ if (x < nx && y < ny && z < nz) {
+
+ // Cell sphere center position
+ const Float3 X = MAKE_FLOAT3(
+ nx*dx + 0.5*dx,
+ ny*dy + 0.5*dy,
+ nz*dz + 0.5*dz);
+
+ Float d, r;
+ Float phi = 0.99;
+
+ // Iterate over 27 neighbor cells
+ for (int z_dim=-1; z_dim<2; ++z_dim) { // z-axis
+ for (int y_dim=-1; y_dim<2; ++y_dim) { // y-axis
+ for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
+
+
+ // Calculate linear cell ID
+ const unsigned int cellID = x + y*devC_grid.num[0]
+ + (devC_grid.num[0] * devC_grid.num[1])*z;
+
+ // Lowest particle index in cell
+ const unsigned int startIdx = dev_cellStart[cellID];
+
+
+ // Make sure cell is not empty
+ if (startIdx != 0xffffffff) {
+
+ // Highest particle index in cell
+ const unsigned int endIdx = dev_cellEnd[cellID];
+
+ // Iterate over cell particles
+ for (unsigned int i = startIdx; i<endIdx; ++i) {
+
+ // Read particle position and radius
+ __syncthreads();
+ xr = dev_x_sorted[i];
+ r = xr.w;
+
+ // Find center distance
+ d = length(MAKE_FLOAT3(
+ X.x - xr.x,
+ X.y - xr.y,
+ X.z - xr.z));
+
+ // if ((R + r) <= d) -> no common volume
+
+ // Lens shaped intersection
+ if (((R - r) < d) && (d < (R + r))) {
+ void_volume -=
+ 1.0/(12.0*d) * (
+ M_PI*(R + r - d)*(R + r - d) *
+ (d*d + 2.0*d*r - 3.0*r*r + 2.0*d*R
+ + 6.0*r*R - 3.0*R*R) );
+
+ // Particle fully contained in cell sphere
+ } else if (d <= (R - r)) {
+ void_volume -= 4.0/3.0*M_PI*r*r*r;
+ }
+ }
+ }
+ }
+ }
}
// Make sure that the porosity is in the interval ]0.0;1.0[
- const Float phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
+ phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
// Save porosity
__syncthreads();
-
dev_d_phi[idx(x,y,z)] = phi;
}
}
t@@ -325,9 +439,9 @@ __global__ void findDarcyTransmissivitiesDev(
const unsigned int nz = devC_grid.num[2];
// Grid sizes
- const Float d_dx = devC_grid.L[0]/nx;
- const Float d_dy = devC_grid.L[1]/ny;
- const Float d_dz = devC_grid.L[2]/nz;
+ const Float dx = devC_grid.L[0]/nx;
+ const Float dy = devC_grid.L[1]/ny;
+ const Float dz = devC_grid.L[2]/nz;
// Density of the fluid [kg/m^3]
const Float rho = 1000.0;
t@@ -340,7 +454,7 @@ __global__ void findDarcyTransmissivitiesDev(
__syncthreads();
- // Cell porosity [-]
+ // Read the cell porosity [-]
const Float phi = dev_d_phi[cellidx];
// Calculate permeability from the Kozeny-Carman relationship
t@@ -350,16 +464,17 @@ __global__ void findDarcyTransmissivitiesDev(
// Schwartz and Zhang 2003
Float k = phi*phi*phi/((1.0-phi)*(1.0-phi)) * d_factor;
-
- __syncthreads();
-
// Save hydraulic conductivity [m/s]
- //const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
- const Float K = 0.5;
- dev_d_K[cellidx] = K;
+ const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
+ //const Float K = 0.5;
+ //const Float K = 1.0e-2;
// Hydraulic transmissivity [m2/s]
- Float3 T = {K*d_dx, K*d_dy, K*d_dz};
+ Float3 T = {K*dx, K*dy, K*dz};
+
+ // Save values. Note! The K values are unused
+ __syncthreads();
+ dev_d_K[cellidx] = K;
dev_d_T[cellidx] = T;
}
t@@ -391,25 +506,29 @@ __global__ void findDarcyGradientsDev(
// Check that we are not outside the fluid grid
Float3 gradient;
+ Float xp, xn, yp, yn, zp, zn;
if (x < nx && y < ny && z < nz) {
+ // Read 6 neighbor cells
__syncthreads();
-
+ xp = dev_scalarfield[idx(x+1,y,z)];
+ xn = dev_scalarfield[idx(x-1,y,z)];
+ yp = dev_scalarfield[idx(x,y+1,z)];
+ yn = dev_scalarfield[idx(x,y-1,z)];
+ zp = dev_scalarfield[idx(x,y,z+1)];
+ zn = dev_scalarfield[idx(x,y,z-1)];
+
+ // Calculate central-difference gradients
// x
- gradient.x =
- (dev_scalarfield[idx(x+1,y,z)] - dev_scalarfield[idx(x-1,y,z)])
- /(2.0*dx);
+ gradient.x = (xp - xn)/(2.0*dx);
// y
- gradient.y =
- (dev_scalarfield[idx(x,y+1,z)] - dev_scalarfield[idx(x,y-1,z)])
- /(2.0*dy);
+ gradient.y = (yp - yn)/(2.0*dy);
// z
- gradient.z =
- (dev_scalarfield[idx(x,y,z+1)] - dev_scalarfield[idx(x,y,z-1)])
- /(2.0*dz);
+ gradient.z = (zp - zn)/(2.0*dz);
+ // Write gradient
__syncthreads();
dev_vectorfield[cellidx] = gradient;
}
t@@ -443,7 +562,7 @@ __global__ void explDarcyStepDev(
const unsigned int ny = devC_grid.num[1];
const unsigned int nz = devC_grid.num[2];
- // Grid sizes
+ // Cell sizes
const Float dx = devC_grid.L[0]/nx;
const Float dy = devC_grid.L[1]/ny;
const Float dz = devC_grid.L[2]/nz;
t@@ -458,52 +577,64 @@ __global__ void explDarcyStepDev(
// new = old + production*timestep + gradient*timestep
// Enforce Dirichlet BC
- if (x == 0 || y == 0 || z == 0 ||
+ /*if (x == 0 || y == 0 || z == 0 ||
x == nx-1 || y == ny-1 || z == nz-1) {
__syncthreads();
dev_d_H_new[cellidx] = dev_d_H[cellidx];
- } else {
-
- // Cell hydraulic conductivity
- __syncthreads();
- //const Float K = dev_d_K[cellidx];
+ } else {*/
- // Cell hydraulic transmissivities
- const Float3 T = dev_d_T[cellidx];
- //const Float Tx = K*dx;
- //const Float Ty = K*dy;
- //const Float Tz = K*dz;
-
- // Harmonic mean of transmissivity
- // (in neg. and pos. direction along axis from cell)
- __syncthreads();
- const Float Tx_n = hmeanDev(T.x, dev_d_T[idx(x-1,y,z)].x);
- const Float Tx_p = hmeanDev(T.x, dev_d_T[idx(x+1,y,z)].x);
- const Float Ty_n = hmeanDev(T.y, dev_d_T[idx(x,y-1,z)].y);
- const Float Ty_p = hmeanDev(T.y, dev_d_T[idx(x,y+1,z)].y);
- const Float Tz_n = hmeanDev(T.z, dev_d_T[idx(x,y,z-1)].z);
- const Float Tz_p = hmeanDev(T.z, dev_d_T[idx(x,y,z+1)].z);
-
- // Cell hydraulic storativity
- const Float S = dev_d_Ss[cellidx]*dx*dy*dz;
-
- // Cell hydraulic head
- const Float H = dev_d_H[cellidx];
-
- // Laplacian operator
- const Float deltaH = devC_dt/S *
- ( Tx_n * (dev_d_H[idx(x-1,y,z)] - H)/(dx*dx)
- + Tx_p * (dev_d_H[idx(x+1,y,z)] - H)/(dx*dx)
- + Ty_n * (dev_d_H[idx(x,y-1,z)] - H)/(dy*dy)
- + Ty_p * (dev_d_H[idx(x,y+1,z)] - H)/(dy*dy)
- + Tz_n * (dev_d_H[idx(x,y,z-1)] - H)/(dy*dz)
- + Tz_p * (dev_d_H[idx(x,y,z+1)] - H)/(dy*dz)
- + dev_d_W[cellidx] );
-
- // Calculate new hydraulic pressure in cell
- __syncthreads();
- dev_d_H_new[cellidx] = H + deltaH;
- }
+ // Read cell and the six neighbor cell hydraulic transmissivities
+ __syncthreads();
+ const Float3 T = dev_d_T[cellidx];
+ const Float3 T_xn = dev_d_T[idx(x-1,y,z)];
+ const Float3 T_xp = dev_d_T[idx(x+1,y,z)];
+ const Float3 T_yn = dev_d_T[idx(x,y-1,z)];
+ const Float3 T_yp = dev_d_T[idx(x,y+1,z)];
+ const Float3 T_zn = dev_d_T[idx(x,y,z-1)];
+ const Float3 T_zp = dev_d_T[idx(x,y,z+1)];
+
+ // Read the cell hydraulic specific storativity
+ const Float Ss = dev_d_Ss[cellidx];
+
+ // Read the cell hydraulic volumetric flux
+ const Float W = dev_d_W[cellidx];
+
+ // Read the cell and the six neighbor cell hydraulic pressures
+ const Float H = dev_d_H[cellidx];
+ const Float H_xn = dev_d_H[idx(x-1,y,z)];
+ const Float H_xp = dev_d_H[idx(x+1,y,z)];
+ const Float H_yn = dev_d_H[idx(x,y-1,z)];
+ const Float H_yp = dev_d_H[idx(x,y+1,z)];
+ const Float H_zn = dev_d_H[idx(x,y,z-1)];
+ const Float H_zp = dev_d_H[idx(x,y,z+1)];
+
+ // Calculate the gradients in the pressure between
+ // the cell and it's six neighbors
+ const Float dH_xn = hmeanDev(T.x, T_xn.x) * (H_xn - H)/(dx*dx);
+ const Float dH_xp = hmeanDev(T.x, T_xp.x) * (H_xp - H)/(dx*dx);
+ const Float dH_yn = hmeanDev(T.y, T_yn.y) * (H_yn - H)/(dy*dy);
+ const Float dH_yp = hmeanDev(T.y, T_yp.y) * (H_yp - H)/(dy*dy);
+ Float dH_zn = hmeanDev(T.z, T_zn.z) * (H_zn - H)/(dz*dz);
+ Float dH_zp = hmeanDev(T.z, T_zp.z) * (H_zp - H)/(dz*dz);
+
+ // Neumann (no-flow) boundary condition at +z and -z boundaries
+ // enforced by a gradient value of 0.0
+ if (z == 0)
+ dH_zn = 0.0;
+ if (z == nz-1)
+ dH_zp = 0.0;
+
+ // Determine the Laplacian operator
+ const Float deltaH = devC_dt/(Ss*dx*dy*dz) *
+ ( dH_xn + dH_xp
+ + dH_yn + dH_yp
+ + dH_zn + dH_zp
+ + W );
+
+ // Write the new hydraulic pressure in cell
+ __syncthreads();
+ dev_d_H_new[cellidx] = H + deltaH;
+ //}
}
}
t@@ -557,59 +688,9 @@ __global__ void findDarcyVelocitiesDev(
}
}
-// Solve Darcy flow on a regular, cubic grid
-/*void DEM::initDarcyDev(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
- 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: "
- << d_nx << "*"
- << d_ny << "*"
- << d_nz << std::endl;
- std::cout << " - Fluid grid cell size: "
- << d_dx << "*"
- << d_dy << "*"
- << d_dz << std::endl;
- }
-
- initDarcyMemDev();
- initDarcyVals();
- findDarcyTransmissivities();
-
- checkDarcyTimestep();
-
- transferDarcyToGlobalDeviceMemory(1);
-}*/
-
// Print final heads and free memory
void DEM::endDarcyDev()
{
- /*FILE* Kfile;
- if ((Kfile = fopen("d_K.txt","w"))) {
- printDarcyArray(Kfile, d_K);
- fclose(Kfile);
- } else {
- fprintf(stderr, "Error, could not open d_K.txt\n");
- }*/
- //printDarcyArray(stdout, d_phi, "d_phi");
- //printDarcyArray(stdout, d_K, "d_K");
- //printDarcyArray(stdout, d_H, "d_H");
- //printDarcyArray(stdout, d_V, "d_V");
freeDarcyMemDev();
}
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -324,7 +324,7 @@ __host__ void DEM::allocateGlobalDeviceMemory(void)
__host__ void DEM::freeGlobalDeviceMemory()
{
if (verbose == 1)
- printf("\nLiberating device memory: ");
+ printf("\nFreeing device memory: ");
// Particle arrays
cudaFree(dev_x);
cudaFree(dev_xysum);
t@@ -467,7 +467,7 @@ __host__ void DEM::transferToGlobalDeviceMemory(int statusmsg)
} else if (darcy == 1) {
transferDarcyToGlobalDeviceMemory(1);
} else {
- std::cerr << "darcy value not understood ("
+ std::cerr << "Error: Darcy value not understood ("
<< darcy << ")" << std::endl;
}
}
t@@ -552,7 +552,9 @@ __host__ void DEM::transferFromGlobalDeviceMemory()
cudaMemcpyDeviceToHost);
#endif
} else {
+#ifdef DARCY_GPU
transferDarcyFromGlobalDeviceMemory(0);
+#endif
}
}
t@@ -698,7 +700,7 @@ __host__ void DEM::startTime()
double t_summation = 0.0;
double t_integrateWalls = 0.0;
- double t_findPorositiesDev = 0.0;
+ double t_findPorositiesCubicDev = 0.0;
double t_findDarcyTransmissivitiesDev = 0.0;
double t_setDarcyGhostNodesDev = 0.0;
double t_explDarcyStepDev = 0.0;
t@@ -898,20 +900,25 @@ __host__ void DEM::startTime()
#ifdef DARCY_GPU
//*
- checkForCudaErrors("Before findPorositiesDev", iter);
+ checkForCudaErrors("Before findPorositiesCubicDev", iter);
// Find cell porosities
if (PROFILING == 1)
startTimer(&kernel_tic);
- findPorositiesDev<<<dimGridFluid, dimBlockFluid>>>(
+ findPorositiesCubicDev<<<dimGridFluid, dimBlockFluid>>>(
dev_cellStart,
dev_cellEnd,
dev_x_sorted,
dev_d_phi);
+ /*findPorositiesSphericalDev<<<dimGridFluid, dimBlockFluid>>>(
+ dev_cellStart,
+ dev_cellEnd,
+ dev_x_sorted,
+ dev_d_phi);*/
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
- &t_findPorositiesDev);
- checkForCudaErrors("Post findPorositiesDev", iter);
+ &t_findPorositiesCubicDev);
+ checkForCudaErrors("Post findPorositiesCubicDev", iter);
// Find resulting cell transmissivities
if (PROFILING == 1)
t@@ -978,7 +985,7 @@ __host__ void DEM::startTime()
&t_findDarcyGradientsDev);
checkForCudaErrors("Post findDarcyGradientsDev", iter);
- // Find the pressure gradients
+ // Find the velocities caused by the pressure gradients
if (PROFILING == 1)
startTimer(&kernel_tic);
findDarcyVelocitiesDev<<<dimGridFluid, dimBlockFluid>>>(
t@@ -1111,6 +1118,11 @@ __host__ void DEM::startTime()
sprintf(file,"output/%s.output%05d.bin", sid.c_str(), time.step_count);
writebin(file);
+ // Write Darcy arrays
+ sprintf(file,"output/%s.d_phi.output%05d.bin", sid.c_str(), time.step_count);
+ writeDarcyArray(d_phi, file);
+ sprintf(file,"output/%s.d_K.output%05d.bin", sid.c_str(), time.step_count);
+ writeDarcyArray(d_K, file);
if (CONTACTINFO == 1) {
// Write contact information to stdout
t@@ -1187,7 +1199,7 @@ __host__ void DEM::startTime()
if (PROFILING == 1 && verbose == 1) {
double t_sum = t_calcParticleCellID + t_thrustsort + t_reorderArrays +
t_topology + t_interact + t_bondsLinear + t_latticeBoltzmannD3Q19 +
- t_integrate + t_summation + t_integrateWalls + t_findPorositiesDev +
+ t_integrate + t_summation + t_integrateWalls + t_findPorositiesCubicDev +
t_findDarcyTransmissivitiesDev + t_setDarcyGhostNodesDev +
t_explDarcyStepDev + t_findDarcyGradientsDev +
t_findDarcyVelocitiesDev;
t@@ -1227,8 +1239,8 @@ __host__ void DEM::startTime()
<< "\t(" << 100.0*t_integrateWalls/t_sum << " %)\n";
if (params.nu > 0.0 && darcy == 1) {
cout
- << " - findPorositiesDev:\t\t" << t_findPorositiesDev/1000.0
- << " s" << "\t(" << 100.0*t_findPorositiesDev/t_sum << " %)\n"
+ << " - findPorositiesCubicDev:\t\t" << t_findPorositiesCubicDev/1000.0
+ << " s" << "\t(" << 100.0*t_findPorositiesCubicDev/t_sum << " %)\n"
<< " - findDarcyTransmis.Dev:\t" <<
t_findDarcyTransmissivitiesDev/1000.0 << " s"
<< "\t(" << 100.0*t_findDarcyTransmissivitiesDev/t_sum << " %)\n"
t@@ -1254,15 +1266,12 @@ __host__ void DEM::startTime()
delete[] k.distmod;
delete[] k.delta_t;
-#ifndef DARCY_GPU
- if (darcy == 1 && params.nu > 0.0)
+ if (darcy == 1 && params.nu > 0.0) {
+#ifdef DARCY_GPU
endDarcyDev();
+#endif
endDarcy();
}
-#else
- if (darcy == 1 && params.nu > 0.0)
- endDarcy();
-#endif
} /* EOF */
// vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -448,7 +448,7 @@ void DEM::writebin(const char *target)
}
ofs.write(as_bytes(params.nu), sizeof(params.nu));
- unsigned int x, y, z;
+ int x, y, z;
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) {
t@@ -462,14 +462,15 @@ void DEM::writebin(const char *target)
}
}
} else if (params.nu > 0.0 && darcy == 1) { // Darcy flow
- for (z = 0; z<d_nz; ++z) {
- for (y = 0; y<d_ny; ++y) {
- for (x = 0; x<d_nx; ++x) {
+ for (z=0; z<d_nz; z++) {
+ for (y=0; y<d_ny; y++) {
+ for (x=0; x<d_nx; 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));
+ //printf("%d,%d,%d: d_H[%d] = %f\n", x,y,z, i, d_H[i]);
}
}
}
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -76,6 +76,8 @@ DEM::DEM(const std::string inputbin,
// Destructor: Liberates dynamically allocated host memory
DEM::~DEM(void)
{
+ if (verbose == 1)
+ std::cout << "Freeing host memory: ";
delete[] k.x;
delete[] k.xysum;
delete[] k.vel;
t@@ -90,6 +92,8 @@ DEM::~DEM(void)
delete[] e.p;
delete[] walls.nx;
delete[] walls.mvfd;
+ if (verbose == 1)
+ std::cout << "Done" << std::endl;
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -325,6 +325,9 @@ class DEM {
void printDarcyArray(FILE* stream, Float3* arr);
void printDarcyArray(FILE* stream, Float3* arr, std::string desc);
+ // Write Darcy arrays to file
+ void writeDarcyArray(Float* array, const char* filename);
+ void writeDarcyArray(Float3* array, const char* filename);
};
#endif