tdarcy flow working - 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 ef0e8ba6bae7051abf968c42029f18fa8352ba10
(DIR) parent 1b27cb6a26bebce2f0f2d81abfa7daf395471e19
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 29 Oct 2013 14:18:47 +0100
darcy flow working
Diffstat:
M src/contactsearch.cuh | 306 ++++++++++++++++----------------
M src/darcy.cpp | 39 ++++++++++++++++++++-----------
M src/darcy.cuh | 236 +++++++++++++++++++++----------
M src/datatypes.h | 1 +
M src/device.cu | 89 ++++++++++++++++---------------
M src/sphere.cpp | 5 ++++-
M src/sphere.h | 15 ++++++++-------
7 files changed, 398 insertions(+), 293 deletions(-)
---
(DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
t@@ -99,77 +99,77 @@ __device__ void findAndProcessContactsInCell(int3 targetCell,
// Get distance modifier for interparticle
// vector, if it crosses a periodic boundary
Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
- if (findDistMod(&targetCell, &distmod) == -1)
- return; // Target cell lies outside the grid
-
-
- //// Check and process particle-particle collisions
-
- // Calculate linear cell ID
- unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
- + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
-
- // Lowest particle index in cell
- unsigned int startIdx = dev_cellStart[cellID];
-
- // Make sure cell is not empty
- if (startIdx != 0xffffffff) {
-
- // Highest particle index in cell + 1
- unsigned int endIdx = dev_cellEnd[cellID];
-
- // Iterate over cell particles
- for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
- if (idx_b != idx_a) { // Do not collide particle with itself
-
- // Fetch position and velocity of particle B.
- Float4 x_b = dev_x_sorted[idx_b];
- Float radius_b = x_b.w;
- Float kappa = devC_params.kappa;
-
- // Distance between particle centers (Float4 -> Float3)
- Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
- x_a.y - x_b.y,
- x_a.z - x_b.z);
-
- // Adjust interparticle vector if periodic boundary/boundaries
- // are crossed
- x_ab += distmod;
- Float x_ab_length = length(x_ab);
-
- // Distance between particle perimeters
- Float delta_ab = x_ab_length - (radius_a + radius_b);
-
- // Check for particle overlap
- if (delta_ab < 0.0f) {
- contactLinearViscous(F, T, es_dot, ev_dot, p,
- idx_a, idx_b,
- dev_vel_sorted,
- dev_angvel_sorted,
- radius_a, radius_b,
- x_ab, x_ab_length,
- delta_ab, kappa);
- } else if (delta_ab < devC_params.db) {
- // Check wether particle distance satisfies the capillary bond distance
- capillaryCohesion_exp(F, radius_a, radius_b, delta_ab,
- x_ab, x_ab_length, kappa);
- }
+ if (findDistMod(&targetCell, &distmod) != -1) {
+
+
+ //// Check and process particle-particle collisions
+
+ // Calculate linear cell ID
+ unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
+ + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
+
+ // Lowest particle index in cell
+ unsigned int startIdx = dev_cellStart[cellID];
+
+ // Make sure cell is not empty
+ if (startIdx != 0xffffffff) {
+
+ // Highest particle index in cell + 1
+ unsigned int endIdx = dev_cellEnd[cellID];
+
+ // Iterate over cell particles
+ for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
+ if (idx_b != idx_a) { // Do not collide particle with itself
+
+ // Fetch position and velocity of particle B.
+ Float4 x_b = dev_x_sorted[idx_b];
+ Float radius_b = x_b.w;
+ Float kappa = devC_params.kappa;
+
+ // Distance between particle centers (Float4 -> Float3)
+ Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
+ x_a.y - x_b.y,
+ x_a.z - x_b.z);
+
+ // Adjust interparticle vector if periodic boundary/boundaries
+ // are crossed
+ x_ab += distmod;
+ Float x_ab_length = length(x_ab);
+
+ // Distance between particle perimeters
+ Float delta_ab = x_ab_length - (radius_a + radius_b);
+
+ // Check for particle overlap
+ if (delta_ab < 0.0f) {
+ contactLinearViscous(F, T, es_dot, ev_dot, p,
+ idx_a, idx_b,
+ dev_vel_sorted,
+ dev_angvel_sorted,
+ radius_a, radius_b,
+ x_ab, x_ab_length,
+ delta_ab, kappa);
+ } else if (delta_ab < devC_params.db) {
+ // Check wether particle distance satisfies the capillary bond distance
+ capillaryCohesion_exp(F, radius_a, radius_b, delta_ab,
+ x_ab, x_ab_length, kappa);
+ }
- // Check wether particles are bonded together
- /*if (bonds.x == idx_b || bonds.y == idx_b ||
- bonds.z == idx_b || bonds.w == idx_b) {
- bondLinear(F, T, es_dot, p, % ev_dot missing
- idx_a, idx_b,
- dev_x_sorted, dev_vel_sorted,
- dev_angvel_sorted,
- radius_a, radius_b,
- x_ab, x_ab_length,
- delta_ab);
- }*/
-
- } // Do not collide particle with itself end
- } // Iterate over cell particles end
- } // Check wether cell is empty end
+ // Check wether particles are bonded together
+ /*if (bonds.x == idx_b || bonds.y == idx_b ||
+ bonds.z == idx_b || bonds.w == idx_b) {
+ bondLinear(F, T, es_dot, p, % ev_dot missing
+ idx_a, idx_b,
+ dev_x_sorted, dev_vel_sorted,
+ dev_angvel_sorted,
+ radius_a, radius_b,
+ x_ab, x_ab_length,
+ delta_ab);
+ }*/
+
+ } // Do not collide particle with itself end
+ } // Iterate over cell particles end
+ } // Check wether cell is empty end
+ } // Periodic boundary and distance adjustment end
} // End of overlapsInCell(...)
t@@ -192,118 +192,118 @@ __device__ void findContactsInCell(int3 targetCell,
// Get distance modifier for interparticle
// vector, if it crosses a periodic boundary
Float3 distmod = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
- if (findDistMod(&targetCell, &distmod) == -1)
- return; // Target cell lies outside the grid
+ if (findDistMod(&targetCell, &distmod) != -1) {
- __syncthreads();
+ __syncthreads();
- //// Check and process particle-particle collisions
+ //// Check and process particle-particle collisions
- // Calculate linear cell ID
- unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
- + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
+ // Calculate linear cell ID
+ unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
+ + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
- // Lowest particle index in cell
- unsigned int startIdx = dev_cellStart[cellID];
+ // Lowest particle index in cell
+ unsigned int startIdx = dev_cellStart[cellID];
- // Make sure cell is not empty
- if (startIdx != 0xffffffff) {
+ // Make sure cell is not empty
+ if (startIdx != 0xffffffff) {
- __syncthreads();
+ __syncthreads();
- // Highest particle index in cell + 1
- unsigned int endIdx = dev_cellEnd[cellID];
+ // Highest particle index in cell + 1
+ unsigned int endIdx = dev_cellEnd[cellID];
- // Read the original index of particle A
- unsigned int idx_a_orig = dev_gridParticleIndex[idx_a];
+ // Read the original index of particle A
+ unsigned int idx_a_orig = dev_gridParticleIndex[idx_a];
- // Iterate over cell particles
- for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
- if (idx_b != idx_a) { // Do not collide particle with itself
+ // Iterate over cell particles
+ for (unsigned int idx_b = startIdx; idx_b<endIdx; ++idx_b) {
+ if (idx_b != idx_a) { // Do not collide particle with itself
- // Fetch position and radius of particle B.
- Float4 x_b = dev_x_sorted[idx_b];
- Float radius_b = x_b.w;
+ // Fetch position and radius of particle B.
+ Float4 x_b = dev_x_sorted[idx_b];
+ Float radius_b = x_b.w;
- // Read the original index of particle B
- unsigned int idx_b_orig = dev_gridParticleIndex[idx_b];
+ // Read the original index of particle B
+ unsigned int idx_b_orig = dev_gridParticleIndex[idx_b];
- //__syncthreads();
+ //__syncthreads();
- // Distance between particle centers (Float4 -> Float3)
- Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
- x_a.y - x_b.y,
- x_a.z - x_b.z);
+ // Distance between particle centers (Float4 -> Float3)
+ Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
+ x_a.y - x_b.y,
+ x_a.z - x_b.z);
- // Adjust interparticle vector if periodic boundary/boundaries
- // are crossed
- x_ab += distmod;
+ // Adjust interparticle vector if periodic boundary/boundaries
+ // are crossed
+ x_ab += distmod;
- Float x_ab_length = length(x_ab);
+ Float x_ab_length = length(x_ab);
- // Distance between particle perimeters
- Float delta_ab = x_ab_length - (radius_a + radius_b);
+ // Distance between particle perimeters
+ Float delta_ab = x_ab_length - (radius_a + radius_b);
- // Check for particle overlap
- if (delta_ab < 0.0f) {
+ // Check for particle overlap
+ if (delta_ab < 0.0f) {
- // If the particle is not yet registered in the contact list,
- // use the next position in the array
- int cpos = *nc;
- unsigned int cidx;
+ // If the particle is not yet registered in the contact list,
+ // use the next position in the array
+ int cpos = *nc;
+ unsigned int cidx;
- // Find out, if particle is already registered in contact list
- for (int i=0; i < devC_nc; ++i) {
- __syncthreads();
- cidx = dev_contacts[(unsigned int)(idx_a_orig*devC_nc+i)];
- if (cidx == devC_np) // Write to position of now-deleted contact
- cpos = i;
- else if (cidx == idx_b_orig) { // Write to position of same contact
- cpos = i;
- break;
+ // Find out, if particle is already registered in contact list
+ for (int i=0; i < devC_nc; ++i) {
+ __syncthreads();
+ cidx = dev_contacts[(unsigned int)(idx_a_orig*devC_nc+i)];
+ if (cidx == devC_np) // Write to position of now-deleted contact
+ cpos = i;
+ else if (cidx == idx_b_orig) { // Write to position of same contact
+ cpos = i;
+ break;
+ }
}
- }
- __syncthreads();
+ __syncthreads();
- // Write the particle index to the relevant position,
- // no matter if it already is there or not (concurrency of write)
- dev_contacts[(unsigned int)(idx_a_orig*devC_nc+cpos)] = idx_b_orig;
+ // Write the particle index to the relevant position,
+ // no matter if it already is there or not (concurrency of write)
+ dev_contacts[(unsigned int)(idx_a_orig*devC_nc+cpos)] = idx_b_orig;
- // Write the interparticle vector and radius of particle B
- //dev_x_ab_r_b[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(x_ab, radius_b);
- dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = MAKE_FLOAT4(distmod.x, distmod.y, distmod.z, radius_b);
+ // Write the interparticle vector and radius of particle B
+ //dev_x_ab_r_b[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(x_ab, radius_b);
+ dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = MAKE_FLOAT4(distmod.x, distmod.y, distmod.z, radius_b);
- // Increment contact counter
- ++*nc;
+ // Increment contact counter
+ ++*nc;
- // Check if the number of contacts of particle A
- // exceeds the max. number of contacts per particle
- if (*nc > devC_nc)
- return; // I would like to throw an error, but it doesn't seem possible...
+ // Check if the number of contacts of particle A
+ // exceeds the max. number of contacts per particle
+ if (*nc > devC_nc)
+ return; // I would like to throw an error, but it doesn't seem possible...
- }
+ }
- // Write the inter-particle position vector correction and radius of particle B
- //dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(distmod, radius_b);
-
- // Check wether particles are bonded together
- /*if (bonds.x == idx_b || bonds.y == idx_b ||
- bonds.z == idx_b || bonds.w == idx_b) {
- bondLinear(F, T, es_dot, p, % ev_dot missing
- idx_a, idx_b,
- dev_x_sorted, dev_vel_sorted,
- dev_angvel_sorted,
- radius_a, radius_b,
- x_ab, x_ab_length,
- delta_ab);
- }*/
-
- } // Do not collide particle with itself end
- } // Iterate over cell particles end
- } // Check wether cell is empty end
+ // Write the inter-particle position vector correction and radius of particle B
+ //dev_distmod[(unsigned int)(idx_a_orig*devC_nc+cpos)] = make_Float4(distmod, radius_b);
+
+ // Check wether particles are bonded together
+ /*if (bonds.x == idx_b || bonds.y == idx_b ||
+ bonds.z == idx_b || bonds.w == idx_b) {
+ bondLinear(F, T, es_dot, p, % ev_dot missing
+ idx_a, idx_b,
+ dev_x_sorted, dev_vel_sorted,
+ dev_angvel_sorted,
+ radius_a, radius_b,
+ x_ab, x_ab_length,
+ delta_ab);
+ }*/
+
+ } // Do not collide particle with itself end
+ } // Iterate over cell particles end
+ } // Check wether cell is empty end
+ } // Periodic boundary and distmod end
} // End of findContactsInCell(...)
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -33,6 +33,7 @@ void DEM::initDarcyMem(const Float cellsizemultiplier)
d.Ss = new Float[ncells]; // hydraulic storativity matrix
d.W = new Float[ncells]; // hydraulic recharge
d.phi = new Float[ncells]; // cell porosity
+ d.dphi = new Float[ncells]; // cell porosity change
}
// Free memory
t@@ -47,6 +48,7 @@ void DEM::freeDarcyMem()
delete[] d.Ss;
delete[] d.W;
delete[] d.phi;
+ delete[] d.dphi;
}
// 3D index to 1D index
t@@ -82,8 +84,8 @@ void DEM::initDarcyVals()
cellidx = idx(ix,iy,iz);
// Hydraulic storativity [-]
- //d.Ss[cellidx] = 1.0;
- 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@@ -113,6 +115,7 @@ void DEM::copyDarcyVals(unsigned int read, unsigned int write)
d.Ss[write] = d.Ss[read];
d.W[write] = d.W[read];
d.phi[write] = d.phi[read];
+ d.dphi[write] = d.dphi[read];
}
// Update ghost nodes from their parent cell values
t@@ -179,12 +182,8 @@ void DEM::findDarcyTransmissivities()
// Density of the fluid [kg/m^3]
const Float rho = 1000.0;
- // Kozeny-Carman parameter
- //Float a = 1.0e-8;
- //Float a = 1.0;
-
// Representative grain radius
- Float r_bar2 = meanRadius()*2.0;
+ Float r_bar2 = meanRadius()*2.0 * 1.0e-6;
// Grain size factor for Kozeny-Carman relationship
Float d_factor = r_bar2*r_bar2/180.0;
t@@ -208,8 +207,7 @@ void DEM::findDarcyTransmissivities()
// Save hydraulic conductivity [m/s]
//K = d.K[cellidx];
- //K = k*rho*-params.g[2]/params.nu;
- K = 0.5;
+ K = k*rho*-params.g[2]/params.nu;
d.K[cellidx] = K;
// Hydraulic transmissivity [m2/s]
t@@ -463,7 +461,7 @@ void DEM::explDarcyStep()
( gradx_n + gradx_p
+ grady_n + grady_p
+ gradz_n + gradz_p
- + d.W[cellidx] );
+ + d.W[cellidx] - d.dphi[cellidx]/time.dt);
// Calculate new hydraulic pressure in cell
d.H_new[cellidx] = H + deltaH;
t@@ -624,7 +622,6 @@ Float DEM::cellPorosity(
// Return the porosity, which should always be ]0.0;1.0[
Float phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
- phi = 0.5;
return phi;
}
t@@ -632,10 +629,14 @@ Float DEM::cellPorosity(
void DEM::findPorosities()
{
int ix, iy, iz, cellidx;
+ Float phi, phi_0;
for (ix=0; ix<d.nx; ++ix) {
for (iy=0; iy<d.ny; ++iy) {
for (iz=0; iz<d.nz; ++iz) {
- d.phi[idx(ix,iy,iz)] = cellPorosity(ix,iy,iz);
+ phi_0 = d.phi[idx(ix,iy,iz)];
+ phi = cellPorosity(ix,iy,iz);
+ d.phi[idx(ix,iy,iz)] = phi;
+ d.dphi[idx(ix,iy,iz)] = phi - phi_0;
}
}
}
t@@ -741,8 +742,9 @@ void DEM::checkDarcyTimestep()
if (value > 0.5) {
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"
+ << " - The transmissivity T (i.e. hydraulic conductivity, K ("
+ << T_max/(d.dx)
+ << ") is too large (" << T_max << ")\n"
<< " - The storativity S is too small"
<< " (" << S_min << ")\n"
<< " - The time step is too large"
t@@ -783,6 +785,15 @@ void DEM::initDarcy()
initDarcyVals();
findDarcyTransmissivities();
+ // set dphi values to zero
+ for (int ix=0; ix<d.nx; ++ix) {
+ for (int iy=0; iy<d.ny; ++iy) {
+ for (int iz=0; iz<d.nz; ++iz) {
+ d.dphi[idx(ix,iy,iz)] = 0.0;
+ }
+ }
+ }
+
checkDarcyTimestep();
}
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -35,6 +35,7 @@ void DEM::initDarcyMemDev(void)
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_dphi, memSizeF); // cell porosity change
checkForCudaErrors("End of initDarcyMemDev");
}
t@@ -51,6 +52,7 @@ void DEM::freeDarcyMemDev()
cudaFree(dev_d_Ss);
cudaFree(dev_d_W);
cudaFree(dev_d_phi);
+ cudaFree(dev_d_dphi);
}
// Transfer to device
t@@ -78,6 +80,7 @@ void DEM::transferDarcyToGlobalDeviceMemory(int statusmsg)
cudaMemcpy(dev_d_Ss, d.Ss, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_d_W, d.W, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_d_phi, d.phi, memSizeF, cudaMemcpyHostToDevice);
+ cudaMemcpy(dev_d_dphi, d.dphi, memSizeF, cudaMemcpyHostToDevice);
checkForCudaErrors("End of transferDarcyToGlobalDeviceMemory");
//if (verbose == 1 && statusmsg == 1)
t@@ -105,6 +108,7 @@ void DEM::transferDarcyFromGlobalDeviceMemory(int statusmsg)
cudaMemcpy(d.Ss, dev_d_Ss, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(d.W, dev_d_W, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(d.phi, dev_d_phi, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(d.dphi, dev_d_dphi, memSizeF, cudaMemcpyDeviceToHost);
checkForCudaErrors("End of transferDarcyFromGlobalDeviceMemory");
if (verbose == 1 && statusmsg == 1)
t@@ -130,7 +134,8 @@ __device__ void copyDarcyValsDev(
Float3* dev_d_V, Float3* dev_d_dH,
Float* dev_d_K, Float3* dev_d_T,
Float* dev_d_Ss, Float* dev_d_W,
- Float* dev_d_phi)
+ Float* dev_d_phi,
+ Float* dev_d_dphi)
{
// Coalesced read
const Float H = dev_d_H[read];
t@@ -142,6 +147,7 @@ __device__ void copyDarcyValsDev(
const Float Ss = dev_d_Ss[read];
const Float W = dev_d_W[read];
const Float phi = dev_d_phi[read];
+ const Float dphi = dev_d_dphi[read];
// Coalesced write
__syncthreads();
t@@ -154,6 +160,7 @@ __device__ void copyDarcyValsDev(
dev_d_Ss[write] = Ss;
dev_d_W[write] = W;
dev_d_phi[write] = phi;
+ dev_d_dphi[write] = dphi;
}
// Update ghost nodes from their parent cell values
t@@ -164,7 +171,8 @@ __global__ void setDarcyGhostNodesDev(
Float3* dev_d_V, Float3* dev_d_dH,
Float* dev_d_K, Float3* dev_d_T,
Float* dev_d_Ss, Float* dev_d_W,
- Float* dev_d_phi)
+ Float* dev_d_phi,
+ Float* dev_d_dphi)
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -192,7 +200,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
if (x == nx-1) {
writeidx = idx(-1,y,z);
t@@ -201,7 +210,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
if (y == 0) {
t@@ -211,7 +221,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
if (y == ny-1) {
writeidx = idx(x,-1,z);
t@@ -220,7 +231,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
if (z == 0) {
t@@ -230,7 +242,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
if (z == nz-1) {
writeidx = idx(x,y,-1);
t@@ -239,7 +252,8 @@ __global__ void setDarcyGhostNodesDev(
dev_d_V, dev_d_dH,
dev_d_K, dev_d_T,
dev_d_Ss, dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
}
}
}
t@@ -251,7 +265,8 @@ __global__ void findPorositiesCubicDev(
unsigned int* dev_cellStart,
unsigned int* dev_cellEnd,
Float4* dev_x_sorted,
- Float* dev_d_phi)
+ Float* dev_d_phi,
+ Float* dev_d_dphi)
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -282,6 +297,10 @@ __global__ void findPorositiesCubicDev(
// Lowest particle index in cell
const unsigned int startIdx = dev_cellStart[cellID];
+ // Read old porosity
+ __syncthreads();
+ Float phi_0 = dev_d_phi[idx(x,y,z)];
+
Float phi = 0.99;
// Make sure cell is not empty
t@@ -305,9 +324,10 @@ __global__ void findPorositiesCubicDev(
phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
}
- // Save porosity
+ // Save porosity and porosity change
__syncthreads();
- dev_d_phi[idx(x,y,z)] = phi;
+ dev_d_phi[idx(x,y,z)] = phi;
+ dev_d_dphi[idx(x,y,z)] = phi - phi_0;
}
}
t@@ -319,7 +339,9 @@ __global__ void findPorositiesSphericalDev(
unsigned int* dev_cellStart,
unsigned int* dev_cellEnd,
Float4* dev_x_sorted,
- Float* dev_d_phi)
+ Float* dev_d_phi,
+ Float* dev_d_dphi,
+ unsigned int iteration)
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -348,60 +370,85 @@ __global__ void findPorositiesSphericalDev(
// Cell sphere center position
const Float3 X = MAKE_FLOAT3(
- nx*dx + 0.5*dx,
- ny*dy + 0.5*dy,
- nz*dz + 0.5*dz);
+ x*dx + 0.5*dx,
+ y*dy + 0.5*dy,
+ z*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];
+ // Read old porosity
+ __syncthreads();
+ Float phi_0 = dev_d_phi[idx(x,y,z)];
- // Iterate over cell particles
- for (unsigned int i = startIdx; i<endIdx; ++i) {
+ // The cell 3d index
+ const int3 gridPos = make_int3((int)x,(int)y,(int)z);
- // Read particle position and radius
- __syncthreads();
- xr = dev_x_sorted[i];
- r = xr.w;
+ // The neighbor cell 3d index
+ int3 targetCell;
- // Find center distance
- d = length(MAKE_FLOAT3(
- X.x - xr.x,
- X.y - xr.y,
- X.z - xr.z));
+ // The distance modifier for particles across periodic boundaries
+ Float3 dist, distmod;
- // if ((R + r) <= d) -> no common volume
+ // 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
- // 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) );
+ // Index of neighbor cell this iteration is looking at
+ targetCell = gridPos + make_int3(x_dim, y_dim, z_dim);
+
+ // Get distance modifier for interparticle
+ // vector, if it crosses a periodic boundary
+ // DISTMOD IS NEVER NOT 0,0,0 !!!!!!
+ distmod = MAKE_FLOAT3(0.0, 0.0, 0.0);
+ if (findDistMod(&targetCell, &distmod) != -1) {
+
+ // Calculate linear cell ID
+ const unsigned int cellID =
+ targetCell.x + targetCell.y * devC_grid.num[0]
+ + (devC_grid.num[0] * devC_grid.num[1])
+ * targetCell.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
+ dist = MAKE_FLOAT3(
+ X.x - xr.x,
+ X.y - xr.y,
+ X.z - xr.z);
+ dist += distmod;
+ d = length(dist);
+
+ // 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;
+ if (d <= R - r) {
+ void_volume -= 4.0/3.0*M_PI*r*r*r;
+ }
}
}
}
t@@ -410,11 +457,19 @@ __global__ void findPorositiesSphericalDev(
}
// Make sure that the porosity is in the interval ]0.0;1.0[
- phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
+ phi = fmin(0.9, fmax(0.1, void_volume/cell_volume));
+ //phi = void_volume/cell_volume;
+
+ Float dphi = phi - phi_0;
+ if (iteration == 0) {
+ // Do not use the initial CPU porosity estimates
+ dphi = 0.0;
+ }
- // Save porosity
+ // Save porosity and porosity change
__syncthreads();
- dev_d_phi[idx(x,y,z)] = phi;
+ dev_d_phi[idx(x,y,z)] = phi;
+ dev_d_dphi[idx(x,y,z)] = dphi;
}
}
t@@ -465,6 +520,7 @@ __global__ void findDarcyTransmissivitiesDev(
Float k = phi*phi*phi/((1.0-phi)*(1.0-phi)) * d_factor;
// Save hydraulic conductivity [m/s]
+ //const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
const Float K = k*rho*-devC_params.g[2]/devC_params.nu;
//const Float K = 0.5;
//const Float K = 1.0e-2;
t@@ -550,7 +606,8 @@ __global__ void explDarcyStepDev(
Float* dev_d_H_new,
Float3* dev_d_T,
Float* dev_d_Ss,
- Float* dev_d_W)
+ Float* dev_d_W,
+ Float* dev_d_dphi)
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -596,8 +653,11 @@ __global__ void explDarcyStepDev(
// 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 hydraulic recharge
+ const Float Q = dev_d_W[cellidx];
+
+ // Read the cell porosity change
+ const Float dphi = dev_d_dphi[cellidx];
// Read the cell and the six neighbor cell hydraulic pressures
const Float H = dev_d_H[cellidx];
t@@ -610,30 +670,56 @@ __global__ void explDarcyStepDev(
// 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);
+ const Float TdH_xn = hmeanDev(T.x, T_xn.x) * (H_xn - H)/(dx*dx);
+ const Float TdH_xp = hmeanDev(T.x, T_xp.x) * (H_xp - H)/(dx*dx);
+ const Float TdH_yn = hmeanDev(T.y, T_yn.y) * (H_yn - H)/(dy*dy);
+ const Float TdH_yp = hmeanDev(T.y, T_yp.y) * (H_yp - H)/(dy*dy);
+ Float TdH_zn = hmeanDev(T.z, T_zn.z) * (H_zn - H)/(dz*dz);
+ Float TdH_zp = hmeanDev(T.z, T_zp.z) * (H_zp - H)/(dz*dz);
+
+ // Mean cell dimension
+ const Float dx_bar = (dx+dy+dz)/3.0;
// 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;
+ TdH_zn = 0.0;
if (z == nz-1)
- dH_zp = 0.0;
+ TdH_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 );
+ const Float laplacianH =
+ TdH_xn + TdH_xp
+ + TdH_yn + TdH_yp
+ + TdH_zn + TdH_zp;
+
+ const Float porevol = -dx_bar*dphi/devC_dt;
+
+ // The change in hydraulic pressure
+ const Float deltaH = devC_dt/(Ss*dx*dy*dz)*(laplacianH + Q + porevol);
+
+ /*printf("%d,%d,%d: H = %f, deltaH = %f,\tlaplacianH = %f,"
+ "dphi = %f,\tpore = %f, "
+ "TdH_xn = %f, TdH_xp = %f, "
+ "TdH_yn = %f, TdH_yp = %f, "
+ "TdH_zn = %f, TdH_zp = %f\n"
+ "\tT_xn = %f,\tT_xp = %f,\t"
+ "T_yn = %f,\tT_yp = %f,\t"
+ "T_zn = %f,\tT_zp = %f\n",
+ x,y,z, H, deltaH, laplacianH, dphi, porevol,
+ TdH_xn, TdH_xp,
+ TdH_yn, TdH_yp,
+ TdH_zn, TdH_zp,
+ T_xn.x, T_xp.x,
+ T_yn.y, T_yp.y,
+ T_zn.z, T_zp.z);*/
+
+ // The pressure should never be negative
+ const Float H_new = fmax(0.0, H + deltaH);
// Write the new hydraulic pressure in cell
__syncthreads();
- dev_d_H_new[cellidx] = H + deltaH;
+ dev_d_H_new[cellidx] = H_new;
//}
}
}
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -120,6 +120,7 @@ struct Darcy {
Float* Ss; // Cell hydraulic storativity
Float* W; // Cell hydraulic recharge
Float* phi; // Cell porosity
+ Float* dphi; // Cell porosity change
};
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -18,8 +18,6 @@
#include "constants.cuh"
#include "debug.h"
-//#include "cuPrintf.cu"
-
#include "sorting.cuh"
#include "contactmodels.cuh"
#include "cohesion.cuh"
t@@ -60,7 +58,8 @@ __host__ void DEM::initializeGPU(void)
cout << " System contains 1 CUDA compatible device.\n";
} else {
if (verbose == 1)
- cout << " System contains " << devicecount << " CUDA compatible devices.\n";
+ cout << " System contains " << devicecount
+ << " CUDA compatible devices.\n";
}
cudaGetDeviceProperties(&prop, cudadevice);
t@@ -77,7 +76,8 @@ __host__ void DEM::initializeGPU(void)
<< cudaRuntimeVersion%100 << std::endl;
}
- // Comment following line when using a system only containing exclusive mode GPUs
+ // Comment following line when using a system only containing exclusive mode
+ // GPUs
cudaChooseDevice(&cudadevice, &prop);
checkForCudaErrors("While initializing CUDA device");
t@@ -380,19 +380,10 @@ __host__ void DEM::freeGlobalDeviceMemory()
#ifdef DARCY_GPU
if (params.nu > 0.0 && darcy == 1) {
- cudaFree(dev_d_H);
- cudaFree(dev_d_H_new);
- cudaFree(dev_d_V);
- cudaFree(dev_d_dH);
- cudaFree(dev_d_K);
- cudaFree(dev_d_T);
- cudaFree(dev_d_Ss);
- cudaFree(dev_d_W);
- cudaFree(dev_d_phi);
+ freeDarcyMemDev();
}
#endif
-
checkForCudaErrors("During cudaFree calls");
if (verbose == 1)
t@@ -584,9 +575,10 @@ __host__ void DEM::transferFromGlobalDeviceMemory()
// Iterate through time by explicit time integration
__host__ void DEM::startTime()
{
- using std::cout; // Namespace directive
- using std::cerr; // Namespace directive
- using std::endl; // Namespace directive
+ using std::cout;
+ using std::cerr;
+ using std::endl;
+
std::string outfile;
char file[200];
FILE *fp;
t@@ -657,7 +649,7 @@ __host__ void DEM::startTime()
// Initialize counter variable values
filetimeclock = 0.0;
long iter = 0;
- const int stdout_report = 10; // the no of time steps between reporting to stdout
+ const int stdout_report = 10; // no of steps between reporting to stdout
// Create first status.dat
//sprintf(file,"output/%s.status.dat", sid);
t@@ -682,7 +674,7 @@ __host__ void DEM::startTime()
// 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;
+ d_factor = r_bar2*r_bar2/180.0 * 1.0e-6;
#endif
} else {
std::cerr << "Error, darcy value (" << darcy
t@@ -719,7 +711,7 @@ __host__ void DEM::startTime()
double t_summation = 0.0;
double t_integrateWalls = 0.0;
- double t_findPorositiesCubicDev = 0.0;
+ double t_findPorositiesDev = 0.0;
double t_findDarcyTransmissivitiesDev = 0.0;
double t_setDarcyGhostNodesDev = 0.0;
double t_explDarcyStepDev = 0.0;
t@@ -759,32 +751,36 @@ __host__ void DEM::startTime()
// Synchronization point
cudaThreadSynchronize();
if (PROFILING == 1)
- stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, &t_calcParticleCellID);
+ stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
+ &t_calcParticleCellID);
checkForCudaErrors("Post calcParticleCellID");
- // Sort particle (key, particle ID) pairs by hash key with Thrust radix sort
+ // Sort particle (key, particle ID) pairs by hash key with Thrust radix
+ // sort
if (PROFILING == 1)
startTimer(&kernel_tic);
thrust::sort_by_key(thrust::device_ptr<uint>(dev_gridParticleCellID),
thrust::device_ptr<uint>(dev_gridParticleCellID + np),
thrust::device_ptr<uint>(dev_gridParticleIndex));
- cudaThreadSynchronize(); // Needed? Does thrust synchronize threads implicitly?
+ cudaThreadSynchronize(); // Needed? Does thrust synchronize implicitly?
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, &t_thrustsort);
checkForCudaErrors("Post thrust::sort_by_key");
- // Zero cell array values by setting cellStart to its highest possible value,
- // specified with pointer value 0xffffffff, which for a 32 bit unsigned int
+ // Zero cell array values by setting cellStart to its highest possible
+ // value, specified with pointer value 0xffffffff, which for a 32 bit
+ // unsigned int
// is 4294967295.
cudaMemset(dev_cellStart, 0xffffffff,
grid.num[0]*grid.num[1]*grid.num[2]*sizeof(unsigned int));
cudaThreadSynchronize();
checkForCudaErrors("Post cudaMemset");
- // Use sorted order to reorder particle arrays (position, velocities, radii) to ensure
- // coherent memory access. Save ordered configurations in new arrays (*_sorted).
+ // Use sorted order to reorder particle arrays (position, velocities,
+ // radii) to ensure coherent memory access. Save ordered configurations
+ // in new arrays (*_sorted).
if (PROFILING == 1)
startTimer(&kernel_tic);
reorderArrays<<<dimGrid, dimBlock, smemSize>>>(dev_cellStart,
t@@ -919,25 +915,28 @@ __host__ void DEM::startTime()
#ifdef DARCY_GPU
- checkForCudaErrors("Before findPorositiesCubicDev", iter);
+ checkForCudaErrors("Before findPorositiesDev", iter);
// Find cell porosities
if (PROFILING == 1)
startTimer(&kernel_tic);
- findPorositiesCubicDev<<<dimGridFluid, dimBlockFluid>>>(
+ /*findPorositiesCubicDev<<<dimGridFluid, dimBlockFluid>>>(
+ dev_cellStart,
+ dev_cellEnd,
+ dev_x_sorted,
+ dev_d_phi,
+ dev_d_dphi);*/
+ findPorositiesSphericalDev<<<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);
+ dev_d_phi,
+ dev_d_dphi,
+ iter);
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
- &t_findPorositiesCubicDev);
- checkForCudaErrors("Post findPorositiesCubicDev", iter);
+ &t_findPorositiesDev);
+ checkForCudaErrors("Post findPorositiesDev", iter);
// Find resulting cell transmissivities
if (PROFILING == 1)
t@@ -965,7 +964,8 @@ __host__ void DEM::startTime()
dev_d_T,
dev_d_Ss,
dev_d_W,
- dev_d_phi);
+ dev_d_phi,
+ dev_d_dphi);
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
t@@ -980,7 +980,8 @@ __host__ void DEM::startTime()
dev_d_H_new,
dev_d_T,
dev_d_Ss,
- dev_d_W);
+ dev_d_W,
+ dev_d_dphi);
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
t@@ -1185,7 +1186,9 @@ __host__ void DEM::startTime()
filetimeclock = 0.0;
}
- //break; // Stop after the first iteration
+
+ // Uncomment break command to stop after the first iteration
+ //break;
}
// Stop clock and display calculation time spent
t@@ -1219,7 +1222,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_findPorositiesCubicDev +
+ t_integrate + t_summation + t_integrateWalls + t_findPorositiesDev +
t_findDarcyTransmissivitiesDev + t_setDarcyGhostNodesDev +
t_explDarcyStepDev + t_findDarcyGradientsDev +
t_findDarcyVelocitiesDev;
t@@ -1259,8 +1262,8 @@ __host__ void DEM::startTime()
<< "\t(" << 100.0*t_integrateWalls/t_sum << " %)\n";
if (params.nu > 0.0 && darcy == 1) {
cout
- << " - findPorositiesCubicDev:\t\t" << t_findPorositiesCubicDev/1000.0
- << " s" << "\t(" << 100.0*t_findPorositiesCubicDev/t_sum << " %)\n"
+ << " - findPorositiesDev:\t\t" << t_findPorositiesDev/1000.0
+ << " s" << "\t(" << 100.0*t_findPorositiesDev/t_sum << " %)\n"
<< " - findDarcyTransmis.Dev:\t" <<
t_findDarcyTransmissivitiesDev/1000.0 << " s"
<< "\t(" << 100.0*t_findDarcyTransmissivitiesDev/t_sum << " %)\n"
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -49,6 +49,10 @@ DEM::DEM(const std::string inputbin,
if (dry == 1)
exit(0);
+ if (params.nu > 0.0 && darcy == 1) {
+ initDarcy();
+ }
+
if (initCuda == 1) {
// Initialize CUDA
t@@ -60,7 +64,6 @@ DEM::DEM(const std::string inputbin,
}
if (params.nu > 0.0 && darcy == 1) {
- initDarcy();
initDarcyMemDev();
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -151,15 +151,16 @@ class DEM {
Darcy d;
// Darcy values, device
- Float* dev_d_H; // Cell hydraulic heads
+ Float* dev_d_H; // Cell hydraulic heads
Float* dev_d_H_new; // Cell hydraulic heads
- Float3* dev_d_V; // Cell fluid velocity
- Float3* dev_d_dH; // Cell spatial gradient in heads
- Float* dev_d_K; // Cell hydraulic conductivities
- Float3* dev_d_T; // Cell hydraulic transmissivity
- Float* dev_d_Ss; // Cell hydraulic storativity
- Float* dev_d_W; // Cell hydraulic recharge
+ Float3* dev_d_V; // Cell fluid velocity
+ Float3* dev_d_dH; // Cell spatial gradient in heads
+ Float* dev_d_K; // Cell hydraulic conductivities
+ Float3* dev_d_T; // Cell hydraulic transmissivity
+ Float* dev_d_Ss; // Cell hydraulic storativity
+ Float* dev_d_W; // Cell hydraulic recharge
Float* dev_d_phi; // Cell porosity
+ Float* dev_d_dphi; // Cell porosity change
//// Darcy functions