tdiv(tau) will give wrong values on first iteration - 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 ad6044c2f75bb6416249b65fc33496a4005c6fd8
(DIR) parent a1210156b64db147b3827d400a4f0ce3dd2059cb
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 23 May 2014 12:42:35 +0200
div(tau) will give wrong values on first iteration
Diffstat:
M src/device.cu | 2 +-
M src/navierstokes.cuh | 29 +++++++++++++++++++----------
M tests/io_tests_fluid.py | 4 ++--
3 files changed, 22 insertions(+), 13 deletions(-)
---
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -833,7 +833,6 @@ __host__ void DEM::startTime()
// Solve Navier Stokes flow through the grid
if (navierstokes == 1) {
-
checkForCudaErrorsIter("Before findPorositiesDev", iter);
// Find cell porosities, average particle velocities, and average
// particle diameters. These are needed for predicting the fluid
t@@ -901,6 +900,7 @@ __host__ void DEM::startTime()
dev_ns_p,
dev_ns_v,
dev_ns_tau,
+ iter,
dev_ns_f_pf,
dev_force);
cudaThreadSynchronize();
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2782,14 +2782,15 @@ __device__ Float dragCoefficient(Float re)
// Find particle-fluid interaction force as outlined by Zhou et al. 2010, and
// originally by Gidaspow 1992.
__global__ void findInteractionForce(
- Float4* dev_x, // in
- Float4* dev_vel, // in
- Float* dev_ns_phi, // in
- Float* dev_ns_p, // in
- Float3* dev_ns_v, // in
- Float* dev_ns_tau, // in
- Float3* dev_ns_f_pf, // out
- Float4* dev_force) // out
+ Float4* dev_x, // in
+ Float4* dev_vel, // in
+ Float* dev_ns_phi, // in
+ Float* dev_ns_p, // in
+ Float3* dev_ns_v, // in
+ Float* dev_ns_tau, // in
+ const unsigned int iter, // in
+ Float3* dev_ns_f_pf, // out
+ Float4* dev_force) // out
{
unsigned int i = threadIdx.x + blockIdx.x*blockDim.x; // Particle index
t@@ -2838,14 +2839,22 @@ __global__ void findInteractionForce(
-1.0*gradient(dev_ns_p, i_x, i_y, i_z, dx, dy, dz)*V_p;
// Viscous force
- const Float3 f_v =
- -1.0*divergence_tensor(dev_ns_tau, i_x, i_y, i_z, dx, dy, dz)*V_p;
+ Float3 f_v;
+ if (iter == 0)
+ f_v = MAKE_FLOAT3(0.0, 0.0, 0.0);
+ else
+ f_v =
+ -1.0*divergence_tensor(dev_ns_tau, i_x, i_y, i_z, dx, dy, dz)
+ *V_p;
// Interaction force on particle (force) and fluid (f_pf)
__syncthreads();
const Float3 f_pf = f_d + f_p + f_v;
#ifdef CHECK_NS_FINITE
+ checkFiniteFloat3("f_d", i_x, i_y, i_z, f_d);
+ checkFiniteFloat3("f_p", i_x, i_y, i_z, f_p);
+ checkFiniteFloat3("f_v", i_x, i_y, i_z, f_v);
checkFiniteFloat3("f_pf", i_x, i_y, i_z, f_pf);
#endif
(DIR) diff --git a/tests/io_tests_fluid.py b/tests/io_tests_fluid.py
t@@ -22,8 +22,8 @@ py.readbin("../input/" + orig.sid + ".bin", verbose=False)
compare(orig, py, "Python IO:")
# Test C++ IO routines
-orig.run(verbose=True, hideinputfile=True)
-#orig.run(verbose=True, hideinputfile=False, cfd=True)
+#orig.run(verbose=True, hideinputfile=True)
+orig.run(verbose=True, hideinputfile=False)
cpp = sphere.sim(fluid=True)
cpp.readbin("../output/" + orig.sid + ".output00000.bin", verbose=False)
compare(orig, cpp, "C++ IO: ")