tc_grad_p replaced by velocity scaling term c_v and advection scaling term c_a - 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 356ed98dc6fbcc6150173c8698fe3bdc15c523b6
(DIR) parent 6985cb64fb0ffd0b5c9376ccd8ee4b2f05183bcf
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 23 Oct 2014 11:52:11 +0200
c_grad_p replaced by velocity scaling term c_v and advection scaling term c_a
Diffstat:
M python/sphere.py | 36 +++++++++++++++++++++-----------
M src/constants.h | 2 +-
M src/datatypes.h | 3 ++-
M src/device.cu | 9 +++++----
M src/file_io.cpp | 6 ++++--
M src/navierstokes.cuh | 78 +++++++++++++------------------
6 files changed, 68 insertions(+), 66 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -19,7 +19,7 @@ numpy.seterr(all='warn', over='raise')
# Sphere version number. This field should correspond to the value in
# `../src/constants.h`.
-VERSION=1.05
+VERSION=1.06
class sim:
'''
t@@ -334,8 +334,11 @@ class sim:
# Porosity scaling factor
self.c_phi = numpy.ones(1, dtype=numpy.float64)
- # Permeability scaling factor
- self.c_grad_p = numpy.ones(1, dtype=numpy.float64)
+ # Fluid velocity scaling factor
+ self.c_v = numpy.ones(1, dtype=numpy.float64)
+
+ # Advection scaling factor
+ self.c_a = numpy.ones(1, dtype=numpy.float64)
## Interaction forces
self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
t@@ -603,7 +606,9 @@ class sim:
elif (self.c_phi != other.c_phi):
print(84)
return(84)
- elif (self.c_grad_p != other.c_grad_p):
+ elif (self.c_v != other.c_v):
+ print(85)
+ elif (self.c_a != other.c_a):
print(85)
return(85)
elif (self.f_d != other.f_d).any():
t@@ -1020,18 +1025,23 @@ class sim:
self.tolerance =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.maxiter = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
- if (self.version >= 1.01):
+ if self.version >= 1.01:
self.ndem = numpy.fromfile(fh, dtype=numpy.uint32, count=1)
else:
self.ndem = 1
- if (self.version >= 1.04):
+ if self.version >= 1.04:
self.c_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
- self.c_grad_p =\
+ self.c_v =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ if self.version >= 1.06:
+ self.c_a =\
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ else:
+ self.c_a = numpy.ones(1, dtype=numpy.float64)
else:
self.c_phi = numpy.ones(1, dtype=numpy.float64)
- self.c_grad_p = numpy.ones(1, dtype=numpy.float64)
+ self.c_v = numpy.ones(1, dtype=numpy.float64)
if self.version >= 1.05:
self.f_d = numpy.empty_like(self.x)
t@@ -1219,7 +1229,8 @@ class sim:
fh.write(self.ndem.astype(numpy.uint32))
fh.write(self.c_phi.astype(numpy.float64))
- fh.write(self.c_grad_p.astype(numpy.float64))
+ fh.write(self.c_v.astype(numpy.float64))
+ fh.write(self.c_a.astype(numpy.float64))
for i in numpy.arange(self.np):
fh.write(self.f_d[i,:].astype(numpy.float64))
t@@ -2764,7 +2775,8 @@ class sim:
self.ndem = numpy.array(1)
self.c_phi = numpy.ones(1, dtype=numpy.float64)
- self.c_grad_p = numpy.ones(1, dtype=numpy.float64)
+ self.c_v = numpy.ones(1, dtype=numpy.float64)
+ self.c_a = numpy.ones(1, dtype=numpy.float64)
self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
t@@ -4622,7 +4634,7 @@ class sim:
self.t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
(self.H50 - H[i_lower])/(H[i_upper] - H[i_lower])
- self.c_v = T50*self.H50**2.0/(self.t50)
+ self.c_coeff = T50*self.H50**2.0/(self.t50)
if self.fluid:
e = numpy.mean(sb.phi[:,:,3:-8]) # ignore boundaries
else:
t@@ -4633,7 +4645,7 @@ class sim:
plt.xlabel('Time [s]')
plt.ylabel('Height [m]')
plt.title('$c_v$ = %.2e m$^2$ s$^{-1}$ at %.1f kPa and $e$ = %.2f' \
- % (self.c_v, sb.w_devs[0]/1000.0, e))
+ % (self.c_coeff, sb.w_devs[0]/1000.0, e))
plt.semilogx(t, H, '+-')
plt.axhline(y = self.H0, color='gray')
plt.axhline(y = self.H50, color='gray')
(DIR) diff --git a/src/constants.h b/src/constants.h
t@@ -16,7 +16,7 @@ const Float PI = 3.14159265358979;
const unsigned int ND = 3;
// Define source code version
-const double VERSION = 1.05;
+const double VERSION = 1.06;
// Max. number of contacts per particle
//const int NC = 16;
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -140,7 +140,8 @@ struct NavierStokes {
unsigned int maxiter; // Solver parameter: Max iterations to perform
unsigned int ndem; // Solver parameter: DEM time steps per CFD step
Float c_phi; // Porosity scaling coefficient
- Float c_grad_p; // Fluid pressure gradient scaling coefficient
+ Float c_v; // Fluid velocity scaling coefficient
+ Float c_a; // Fluid advection scaling coefficient
Float4* f_d; // Drag force on particles
Float4* f_p; // Pressure force on particles
Float4* f_v; // Viscous force on particles
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1126,7 +1126,7 @@ __host__ void DEM::startTime()
dev_ns_div_tau_x,
dev_ns_div_tau_y,
dev_ns_div_tau_z,
- //ns.c_grad_p,
+ //ns.c_v,
dev_ns_f_pf,
dev_force,
dev_ns_f_d,
t@@ -1320,7 +1320,8 @@ __host__ void DEM::startTime()
dev_ns_F_pf,
ns.ndem,
wall0_iz,
- ns.c_grad_p,
+ ns.c_v,
+ ns.c_a,
dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z);
t@@ -1466,7 +1467,7 @@ __host__ void DEM::startTime()
dev_ns_v_p_z,
nijac,
ns.ndem,
- ns.c_grad_p,
+ ns.c_v,
dev_ns_f1,
dev_ns_f2,
dev_ns_f);
t@@ -1648,7 +1649,7 @@ __host__ void DEM::startTime()
ns.bc_bot,
ns.bc_top,
ns.ndem,
- ns.c_grad_p,
+ ns.c_v,
wall0_iz,
iter,
dev_ns_v_x,
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -295,7 +295,8 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.ndem), sizeof(unsigned int));
ifs.read(as_bytes(ns.c_phi), sizeof(Float));
- ifs.read(as_bytes(ns.c_grad_p), sizeof(Float));
+ ifs.read(as_bytes(ns.c_v), sizeof(Float));
+ ifs.read(as_bytes(ns.c_a), sizeof(Float));
for (i = 0; i<np; ++i) {
ifs.read(as_bytes(ns.f_d[i].x), sizeof(Float));
t@@ -528,7 +529,8 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(ns.ndem), sizeof(unsigned int));
ofs.write(as_bytes(ns.c_phi), sizeof(Float));
- ofs.write(as_bytes(ns.c_grad_p), sizeof(Float));
+ ofs.write(as_bytes(ns.c_v), sizeof(Float));
+ ofs.write(as_bytes(ns.c_a), sizeof(Float));
for (i = 0; i<np; ++i) {
ofs.write(as_bytes(ns.f_d[i].x), sizeof(Float));
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2289,7 +2289,8 @@ __global__ void findPredNSvelocities(
const Float3* __restrict__ dev_ns_F_pf, // in
const unsigned int ndem, // in
const unsigned int wall0_iz, // in
- const Float c_grad_p, // in
+ const Float c_v, // in
+ const Float c_a, // in
Float* __restrict__ dev_ns_v_p_x, // out
Float* __restrict__ dev_ns_v_p_y, // out
Float* __restrict__ dev_ns_v_p_z) // out
t@@ -2397,12 +2398,12 @@ __global__ void findPredNSvelocities(
(p - p_zn)/dz);
#ifdef SET_1
//pressure_term = -beta*dt/(rho*phi)*grad_p;
- //pressure_term = -beta*c_grad_p*dt/rho*grad_p;
- pressure_term = -beta*c_grad_p*dt/(rho*phi)*grad_p;
+ //pressure_term = -beta*dt/rho*grad_p;
+ pressure_term = -beta*dt/(rho*phi)*grad_p;
#endif
#ifdef SET_2
//pressure_term = -beta*dt/rho*grad_p;
- pressure_term = -beta*c_grad_p*dt/rho*grad_p;
+ pressure_term = -beta*dt/rho*grad_p;
#endif
}
t@@ -2413,20 +2414,20 @@ __global__ void findPredNSvelocities(
// Determine the predicted velocity
#ifdef SET_1
- const Float3 interaction_term = -dt/(rho*phi)*f_i;
- const Float3 diffusion_term = dt/(rho*phi)*div_tau;
- const Float3 gravity_term = MAKE_FLOAT3(
+ const Float3 interaction_term = -c_v*dt/(rho*phi)*f_i;
+ const Float3 diffusion_term = c_v*dt/(rho*phi)*div_tau;
+ const Float3 gravity_term = c_v*MAKE_FLOAT3(
devC_params.g[0], devC_params.g[1], devC_params.g[2])*dt;
- const Float3 porosity_term = -1.0*v*dphi/phi;
- const Float3 advection_term = -1.0*div_phi_vi_v*dt/phi;
+ const Float3 porosity_term = -c_v*v*dphi/phi;
+ const Float3 advection_term = -c_v*c_a*div_phi_vi_v*dt/phi;
#endif
#ifdef SET_2
- const Float3 interaction_term = -dt/(rho*phi)*f_i;
- const Float3 diffusion_term = dt/rho*div_tau;
- const Float3 gravity_term = MAKE_FLOAT3(
+ const Float3 interaction_term = -c_v*dt/(rho*phi)*f_i;
+ const Float3 diffusion_term = c_v*dt/rho*div_tau;
+ const Float3 gravity_term = c_v*MAKE_FLOAT3(
devC_params.g[0], devC_params.g[1], devC_params.g[2])*dt;
- const Float3 porosity_term = -1.0*v*dphi/phi;
- const Float3 advection_term = -1.0*div_phi_vi_v*dt/phi;
+ const Float3 porosity_term = -c_v*v*dphi/phi;
+ const Float3 advection_term = -c_v*c_a*div_phi_vi_v*dt/phi;
#endif
Float3 v_p = v
t@@ -2511,7 +2512,7 @@ __global__ void findNSforcing(
const Float* __restrict__ dev_ns_v_p_z, // in
const unsigned int nijac, // in
const unsigned int ndem, // in
- const Float c_grad_p, // in
+ const Float c_v, // in
Float* __restrict__ dev_ns_f1, // out
Float3* __restrict__ dev_ns_f2, // out
Float* __restrict__ dev_ns_f) // out
t@@ -2563,30 +2564,20 @@ __global__ void findNSforcing(
// Find forcing function terms
#ifdef SET_1
- //const Float t1 = phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
- //const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(c_grad_p*dt);
- //const Float t4 = dphi*devC_params.rho_f/(c_grad_p*dt*dt);
-
- //const Float t1 = phi*phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
- //const Float t2 =
- //devC_params.rho_f*phi*dot(v_p, grad_phi)/(c_grad_p*dt);
- //const Float t4 = dphi*devC_params.rho_f*phi/(c_grad_p*dt*dt);
-
- //const Float t1 = devC_params.rho_f*div_v_p/(c_grad_p*dt);
- //const Float t2 =
- //devC_params.rho_f*dot(v_p, grad_phi)/(phi*dt*c_grad_p);
- //const Float t4 = devC_params.rho_f*dphi/(dt*dt*c_grad_p*phi);
-
- const Float t1 = devC_params.rho_f*phi/(c_grad_p*dt)*div_v_p;
- const Float t2 = devC_params.rho_f/(c_grad_p*dt)*dot(grad_phi, v_p);
- const Float t4 = devC_params.rho_f*dphi/(dt*dt*c_grad_p);
-
+ //const Float t1 = devC_params.rho_f*phi/dt*div_v_p;
+ //const Float t2 = devC_params.rho_f/dt*dot(grad_phi, v_p);
+ //const Float t4 = devC_params.rho_f*dphi/(dt*dt);
+ const Float t1 = devC_params.rho_f*phi/(c_v*dt)*div_v_p;
+ const Float t2 = devC_params.rho_f/(c_v*dt)*dot(grad_phi, v_p);
+ const Float t4 = devC_params.rho_f*dphi/(c_v*dt*dt);
#endif
#ifdef SET_2
- const Float t1 = devC_params.rho_f*div_v_p/(c_grad_p*dt);
- const Float t2 =
- devC_params.rho_f*dot(v_p, grad_phi)/(c_grad_p*phi*dt);
- const Float t4 = dphi*devC_params.rho_f/(c_grad_p*dt*dt*phi);
+ //const Float t1 = devC_params.rho_f*div_v_p/dt;
+ //const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(phi*dt);
+ //const Float t4 = dphi*devC_params.rho_f/(dt*dt*phi);
+ const Float t1 = devC_params.rho_f*div_v_p/(c_v*dt);
+ const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(c_v*phi*dt);
+ const Float t4 = dphi*devC_params.rho_f/(c_v*dt*dt*phi);
#endif
f1 = t1 + t2 + t4;
f2 = grad_phi/phi; // t3/grad(epsilon)
t@@ -2944,7 +2935,7 @@ __global__ void updateNSvelocity(
const int bc_bot, // in
const int bc_top, // in
const unsigned int ndem, // in
- const Float c_grad_p, // in
+ const Float c_v, // in
const unsigned int wall0_iz, // in
const unsigned int iter, // in
Float* __restrict__ dev_ns_v_x, // out
t@@ -3003,15 +2994,10 @@ __global__ void updateNSvelocity(
// Find new velocity
#ifdef SET_1
- //Float3 v = v_p - ndem*devC_dt/(devC_params.rho_f*phi)*grad_epsilon;
- //Float3 v = v_p - ndem*devC_dt*c_grad_p/devC_params.rho_f*grad_epsilon;
- Float3 v = v_p
- -ndem*devC_dt*c_grad_p/(phi*devC_params.rho_f)*grad_epsilon;
+ Float3 v = v_p - c_v*ndem*devC_dt/(phi*devC_params.rho_f)*grad_epsilon;
#endif
#ifdef SET_2
- //Float3 v = v_p - ndem*devC_dt/devC_params.rho_f*grad_epsilon;
- Float3 v = v_p
- -ndem*devC_dt*c_grad_p/devC_params.rho_f*grad_epsilon;
+ Float3 v = v_p - c_v*ndem*devC_dt/(devC_params.rho_f*grad_epsilon;
#endif
// Print values for debugging
t@@ -3189,7 +3175,7 @@ __global__ void findInteractionForce(
const Float* __restrict__ dev_ns_div_tau_x,// in
const Float* __restrict__ dev_ns_div_tau_y,// in
const Float* __restrict__ dev_ns_div_tau_z,// in
- //const Float c_grad_p, // in
+ //const Float c_v, // in
Float3* __restrict__ dev_ns_f_pf, // out
Float4* __restrict__ dev_force, // out
Float4* __restrict__ dev_ns_f_d, // out