tadd scaling coefficients c_phi and c_grad_p - 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 9ace0175d14ee2a2ea10c135ac5598d07c059c61
(DIR) parent c4dcd878677647a252f73067a71dbb3aad50570a
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 28 Jul 2014 15:36:44 +0200
add scaling coefficients c_phi and c_grad_p
Diffstat:
M python/sphere.py | 23 +++++++++++++++++++++++
M src/datatypes.h | 18 ++++++++++--------
M src/device.cu | 5 ++++-
M src/file_io.cpp | 6 ++++++
M src/navierstokes.cuh | 13 ++++++++-----
5 files changed, 51 insertions(+), 14 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -331,6 +331,12 @@ class sim:
# The number of DEM time steps to perform between CFD updates
self.ndem = numpy.array(1)
+ # 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)
+
# Particle color marker
self.color = numpy.zeros(self.np, dtype=numpy.int32)
t@@ -588,6 +594,12 @@ class sim:
elif (self.ndem != other.ndem):
print(83)
return 83
+ elif (self.c_phi != other.c_phi):
+ print(84)
+ return(84)
+ elif (self.c_grad_p != other.c_grad_p):
+ print(85)
+ return(85)
if ((self.color != other.color)).any():
print(90)
t@@ -956,6 +968,11 @@ class sim:
else:
self.ndem = 1
+ if (self.version >= 1.04):
+ self.c_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.c_grad_p =\
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
+
if (self.version >= 1.02):
self.color =\
numpy.fromfile(fh, dtype=numpy.int32, count=self.np)
t@@ -1109,6 +1126,9 @@ class sim:
fh.write(self.maxiter.astype(numpy.uint32))
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.color.astype(numpy.int32))
finally:
t@@ -2484,6 +2504,9 @@ class sim:
self.maxiter = numpy.array(1e4)
self.ndem = numpy.array(1)
+ self.c_phi = numpy.array(1.0)
+ self.c_grad_p = numpy.array(1.0)
+
def setFluidBottomNoFlow(self):
'''
Set the lower boundary of the fluid domain to follow the no-flow
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -121,14 +121,14 @@ struct NavierStokes {
//Float* v_p_x; // Predicted fluid velocity in staggered grid
//Float* v_p_y; // Predicted fluid velocity in staggered grid
//Float* v_p_z; // Predicted fluid velocity in staggered grid
- Float* phi; // Cell porosity
- Float* dphi; // Cell porosity change
- Float* norm; // Normalized residual of epsilon updates
- Float* epsilon; // Iterative solution parameter
- Float* epsilon_new; // Updated value of iterative solution parameter
- Float p_mod_A; // Pressure modulation amplitude at top
- Float p_mod_f; // Pressure modulation frequency at top
- Float p_mod_phi; // Pressure modulation phase at top
+ Float* phi; // Cell porosity
+ Float* dphi; // Cell porosity change
+ Float* norm; // Normalized residual of epsilon updates
+ Float* epsilon; // Iterative solution parameter
+ Float* epsilon_new; // Updated value of iterative solution parameter
+ Float p_mod_A; // Pressure modulation amplitude at top
+ Float p_mod_f; // Pressure modulation frequency at top
+ Float p_mod_phi; // Pressure modulation phase at top
int bc_bot; // 0: Dirichlet, 1: Neumann
int bc_top; // 0: Dirichlet, 1: Neumann
int free_slip_bot; // 0: no, 1: yes
t@@ -139,6 +139,8 @@ struct NavierStokes {
Float tolerance; // Solver parameter: Max residual tolerance
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
};
// Image structure
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1031,7 +1031,8 @@ __host__ void DEM::startTime()
dev_ns_vp_avg,
dev_ns_d_avg,
iter,
- np);
+ np,
+ ns.c_phi);
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
t@@ -1278,6 +1279,7 @@ __host__ void DEM::startTime()
ns.beta,
dev_ns_F_pf,
ns.ndem,
+ ns.c_grad_p,
dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z);
t@@ -1641,6 +1643,7 @@ __host__ void DEM::startTime()
ns.bc_bot,
ns.bc_top,
ns.ndem,
+ ns.c_grad_p,
dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -294,6 +294,9 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.maxiter), sizeof(unsigned int));
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));
+
if (verbose == 1)
cout << "Done" << std::endl;
}
t@@ -502,6 +505,9 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(ns.tolerance), sizeof(Float));
ofs.write(as_bytes(ns.maxiter), sizeof(unsigned int));
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));
}
for (i = 0; i<np; ++i)
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -1003,7 +1003,8 @@ __global__ void findPorositiesVelocitiesDiametersSpherical(
Float3* dev_ns_vp_avg,
Float* dev_ns_d_avg,
const unsigned int iteration,
- const unsigned int np)
+ const unsigned int np,
+ const Float c_phi)
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -1161,8 +1162,8 @@ __global__ void findPorositiesVelocitiesDiametersSpherical(
__syncthreads();
//phi = 0.5; dphi = 0.0; // disable porosity effects
const unsigned int cellidx = idx(x,y,z);
- dev_ns_phi[cellidx] = phi;
- dev_ns_dphi[cellidx] = dphi;
+ dev_ns_phi[cellidx] = phi*c_phi;
+ dev_ns_dphi[cellidx] = dphi*c_phi;
dev_ns_vp_avg[cellidx] = v_avg;
dev_ns_d_avg[cellidx] = d_avg;
t@@ -2145,6 +2146,7 @@ __global__ void findPredNSvelocities(
Float beta, // in
Float3* dev_ns_F_pf, // in
unsigned int ndem, // in
+ Float c_grad_p, // in
Float* dev_ns_v_p_x, // out
Float* dev_ns_v_p_y, // out
Float* dev_ns_v_p_z) // out
t@@ -2248,7 +2250,7 @@ __global__ void findPredNSvelocities(
const Float3 grad_p = MAKE_FLOAT3(
(p - p_xn)/dx,
(p - p_yn)/dy,
- (p - p_zn)/dz);
+ (p - p_zn)/dz) * c_grad_p;
#ifdef SET_1
pressure_term = -beta*dt/(rho*phi)*grad_p;
#endif
t@@ -2754,6 +2756,7 @@ __global__ void updateNSvelocity(
int bc_bot, // in
int bc_top, // in
unsigned int ndem, // in
+ Float c_grad_p, // in
Float* dev_ns_v_x, // out
Float* dev_ns_v_y, // out
Float* dev_ns_v_z) // out
t@@ -2806,7 +2809,7 @@ __global__ void updateNSvelocity(
const Float3 grad_epsilon = MAKE_FLOAT3(
(epsilon_c - epsilon_xn)/dx,
(epsilon_c - epsilon_yn)/dy,
- (epsilon_c - epsilon_zn)/dz);
+ (epsilon_c - epsilon_zn)/dz) * c_grad_p;
// Find new velocity
#ifdef SET_1