tremove `c_a` parameter and reuse the slot for `dt_dem_fac` - 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 e6181dcce46e3632a3f4a18ab0aa4c5db4a008a1
(DIR) parent a9080ee20972d9541be77129557cfd414f69a3a1
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 28 Oct 2014 09:35:59 +0100
remove `c_a` parameter and reuse the slot for `dt_dem_fac`
Diffstat:
M python/halfshear-starter-rate.py | 5 +----
M python/shear-results-velfac.py | 8 ++++----
M python/shear-results.py | 23 ++++++++++++++---------
M python/sphere.py | 20 ++++++++++++--------
M src/constants.h | 2 +-
M src/datatypes.h | 2 +-
M src/device.cu | 1 -
M src/file_io.cpp | 4 ++--
M src/navierstokes.cuh | 9 ++++-----
M tests/CMakeLists.txt | 4 ++--
D tests/cfd_tests_neumann-c_a=0.0-c_… | 80 -------------------------------
A tests/cfd_tests_neumann-c_v=0.1.py | 78 +++++++++++++++++++++++++++++++
12 files changed, 119 insertions(+), 117 deletions(-)
---
(DIR) diff --git a/python/halfshear-starter-rate.py b/python/halfshear-starter-rate.py
t@@ -4,7 +4,7 @@ import numpy
import sys
# launch with:
-# $ python halfshear-starter-rate.py <DEVICE> <FLUID> <C_PHI> <C_V> <SIGMA_0> <VELOCITY FACTOR> <C_A>
+# $ python halfshear-starter-rate.py <DEVICE> <FLUID> <C_PHI> <C_V> <SIGMA_0> <VELOCITY FACTOR>
device = int(sys.argv[1])
wet = int(sys.argv[2])
t@@ -12,7 +12,6 @@ c_phi = float(sys.argv[3])
c_v = float(sys.argv[4])
sigma0 = float(sys.argv[5])
velfac = float(sys.argv[6])
-c_a = float(sys.argv[7])
#sim = sphere.sim('diffusivity-sigma0=' + str(sigma0) + '-c_phi=' + \
# str(c_phi) + '-c_grad_p=' + str(c_grad_p), fluid=True)
t@@ -30,7 +29,6 @@ sim.readlast()
sim.fluid = fluid
if fluid:
sim.id('halfshear-sigma0=' + str(sigma0) + '-c_v=' + str(c_v) + \
- '-c_a=' + str(c_a) + \
'-velfac=' + str(velfac) + '-shear')
else:
sim.id('halfshear-sigma0=' + str(sigma0) + \
t@@ -54,7 +52,6 @@ if fluid:
sim.setMaxIterations(2e5)
sim.c_phi[0] = c_phi
sim.c_v[0] = c_v
- sim.c_a[0] = c_a
sim.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
sim.w_devs[0] = sigma0
(DIR) diff --git a/python/shear-results-velfac.py b/python/shear-results-velfac.py
t@@ -206,17 +206,17 @@ for c in numpy.arange(len(velfacvals)):
if smoothed_results:
ax1.plot(shear_strain[c][1:], friction_smooth[c][1:], \
- label='$c_\\dot{\\gamma}$ = %.2f' % (velfacvals[c]), linewidth=1, alpha=0.5)
+ label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
else:
ax1.plot(shear_strain[c][1:], friction[c][1:], \
- label='$c_\\dot{\\gamma}$ = %.2f' % (velfacvals[c]), linewidth=1, alpha=0.5)
+ label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
ax2.plot(shear_strain[c][1:], dilation[c][1:], \
- label='$c_\\dot{\\gamma}$ = %.2f' % (velfacvals[c]), linewidth=1, alpha=0.5)
+ label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
if zflow and fluid:
ax3.plot(shear_strain[c][1:], v_f_z_mean[c][1:],
- label='$c_\\dot{\\gamma}$ = %.2f' % (velfacvals[c]), linewidth=1, alpha=0.5)
+ label='velfac = %.1f' % (velfacvals[c]), linewidth=1, alpha=0.5)
'''
(DIR) diff --git a/python/shear-results.py b/python/shear-results.py
t@@ -16,7 +16,7 @@ import matplotlib.pyplot as plt
smoothed_results = False
contact_forces = False
pressures = False
-zflow = True
+zflow = False
#sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
#sigma0 = 10.0e3
t@@ -195,8 +195,10 @@ for c in numpy.arange(1,len(cvals)+1):
c += 1
-#fig = plt.figure(figsize=(8,8)) # (w,h)
-fig = plt.figure(figsize=(8,10))
+if zflow:
+ fig = plt.figure(figsize=(8,10))
+else:
+ fig = plt.figure(figsize=(8,8)) # (w,h)
#fig = plt.figure(figsize=(8,12))
#fig = plt.figure(figsize=(8,16))
fig.subplots_adjust(hspace=0.0)
t@@ -204,11 +206,13 @@ fig.subplots_adjust(hspace=0.0)
#plt.subplot(3,1,1)
#plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
-#ax1 = plt.subplot(211)
-#ax2 = plt.subplot(212, sharex=ax1)
-ax1 = plt.subplot(311)
-ax2 = plt.subplot(312, sharex=ax1)
-ax3 = plt.subplot(313, sharex=ax1)
+if zflow:
+ ax1 = plt.subplot(311)
+ ax2 = plt.subplot(312, sharex=ax1)
+ ax3 = plt.subplot(313, sharex=ax1)
+else:
+ ax1 = plt.subplot(211)
+ ax2 = plt.subplot(212, sharex=ax1)
#ax3 = plt.subplot(413, sharex=ax1)
#ax4 = plt.subplot(414, sharex=ax1)
alpha = 0.5
t@@ -260,7 +264,8 @@ for c in numpy.arange(1,len(cvals)+1):
ax1.set_ylabel('Shear friction $\\tau/\\sigma\'$ [-]')
ax2.set_ylabel('Dilation $\\Delta h/(2r)$ [-]')
-ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
+if zflow:
+ ax3.set_ylabel('$\\boldsymbol{v}_\\text{f}^z h$ [ms$^{-1}$]')
#ax3.set_ylabel('Fluid pressure $p_\\text{f}$ [kPa]')
#ax4.set_ylabel('Particle contact force $||\\boldsymbol{f}_\\text{p}||$ [N]')
(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.06
+VERSION=1.07
class sim:
'''
t@@ -337,8 +337,8 @@ class sim:
# 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)
+ # DEM-CFD time scaling factor
+ self.dt_dem_fac = numpy.ones(1, dtype=numpy.float64)
## Interaction forces
self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
t@@ -608,7 +608,7 @@ class sim:
return(84)
elif (self.c_v != other.c_v):
print(85)
- elif (self.c_a != other.c_a):
+ elif (self.dt_dem_fac != other.dt_dem_fac):
print(85)
return(85)
elif (self.f_d != other.f_d).any():
t@@ -1031,12 +1031,16 @@ class sim:
self.ndem = 1
if self.version >= 1.04:
- self.c_phi = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.c_phi = \
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.c_v =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- if self.version >= 1.06:
+ if self.version == 1.06:
self.c_a =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ elif self.version >= 1.07:
+ self.dt_dem_fac =\
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
else:
self.c_a = numpy.ones(1, dtype=numpy.float64)
else:
t@@ -1230,7 +1234,7 @@ class sim:
fh.write(self.c_phi.astype(numpy.float64))
fh.write(self.c_v.astype(numpy.float64))
- fh.write(self.c_a.astype(numpy.float64))
+ fh.write(self.dt_dem_fac.astype(numpy.float64))
for i in numpy.arange(self.np):
fh.write(self.f_d[i,:].astype(numpy.float64))
t@@ -2776,7 +2780,7 @@ class sim:
self.c_phi = 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.dt_dem_fac = 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)
(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.06;
+const double VERSION = 1.07;
// Max. number of contacts per particle
//const int NC = 16;
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -141,7 +141,7 @@ struct NavierStokes {
unsigned int ndem; // Solver parameter: DEM time steps per CFD step
Float c_phi; // Porosity scaling coefficient
Float c_v; // Fluid velocity scaling coefficient
- Float c_a; // Fluid advection scaling coefficient
+ Float dt_dem_fac; // DEM-CFD time 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@@ -1321,7 +1321,6 @@ __host__ void DEM::startTime()
ns.ndem,
wall0_iz,
ns.c_v,
- ns.c_a,
dev_ns_v_p_x,
dev_ns_v_p_y,
dev_ns_v_p_z);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -296,7 +296,7 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.c_phi), sizeof(Float));
ifs.read(as_bytes(ns.c_v), sizeof(Float));
- ifs.read(as_bytes(ns.c_a), sizeof(Float));
+ ifs.read(as_bytes(ns.dt_dem_fac), sizeof(Float));
for (i = 0; i<np; ++i) {
ifs.read(as_bytes(ns.f_d[i].x), sizeof(Float));
t@@ -530,7 +530,7 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(ns.c_phi), sizeof(Float));
ofs.write(as_bytes(ns.c_v), sizeof(Float));
- ofs.write(as_bytes(ns.c_a), sizeof(Float));
+ ofs.write(as_bytes(ns.dt_dem_fac), 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@@ -2290,7 +2290,6 @@ __global__ void findPredNSvelocities(
const unsigned int ndem, // in
const unsigned int wall0_iz, // 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@@ -2429,7 +2428,7 @@ __global__ void findPredNSvelocities(
+ diffusion_term
+ gravity_term
+ porosity_term
- + c_a*advection_term);
+ + advection_term);
//// Neumann BCs
if ((z == 0 && bc_bot == 1) || (z > nz-1 && bc_top == 1))
t@@ -2471,9 +2470,9 @@ __global__ void findPredNSvelocities(
c_v*porosity_term.x,
c_v*porosity_term.y,
c_v*porosity_term.z,
- c_v*c_a*advection_term.x,
- c_v*c_a*advection_term.y,
- c_v*c_a*advection_term.z);
+ c_v*advection_term.x,
+ c_v*advection_term.y,
+ c_v*advection_term.z);
#endif
// Save the predicted velocity
(DIR) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
t@@ -27,8 +27,8 @@ add_test(cfd_tests ${PYTHON_EXECUTABLE}
add_test(cfd_tests_neumann ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/cfd_tests_neumann.py)
-add_test(cfd_tests_neumann-c_a=0.0-c_v=0.1 ${PYTHON_EXECUTABLE}
- ${CMAKE_CURRENT_BINARY_DIR}/cfd_tests_neumann-c_a=0.0-c_v=0.1.py)
+add_test(cfd_tests_neumann-c_v=0.1 ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/cfd_tests_neumann-c_v=0.1.py)
add_test(fluid_particle_interaction ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/fluid_particle_interaction.py)
(DIR) diff --git a/tests/cfd_tests_neumann-c_a=0.0-c_v=0.1.py b/tests/cfd_tests_neumann-c_a=0.0-c_v=0.1.py
t@@ -1,80 +0,0 @@
-#!/usr/bin/env python
-from pytestutils import *
-
-import sphere
-import sys
-import numpy
-import matplotlib.pyplot as plt
-
-print('### CFD tests - Dirichlet/Neumann BCs ###')
-
-print('''# Neumann bottom, Dirichlet top BC.
-# No gravity, no pressure gradients => no flow''')
-orig = sphere.sim("neumann", fluid = True)
-cleanup(orig)
-orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
-orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
-orig.initFluid(mu = 8.9e-4)
-#orig.initFluid(mu = 0.0)
-orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
-orig.c_v[0] = 0.1
-orig.c_a[0] = 0.0
-#orig.c_phi[0] = 0.1
-py = sphere.sim(sid = orig.sid, fluid = True)
-orig.bc_bot[0] = 1 # No-flow BC at bottom (Neumann)
-#orig.run(dry=True)
-orig.run(verbose=False)
-#orig.run(device=2)
-#orig.writeVTKall()
-py.readlast(verbose = False)
-ones = numpy.ones((orig.num))
-py.readlast(verbose = False)
-compareNumpyArraysClose(ones, py.p_f, "Conservation of pressure:",
- tolerance = 1.0e-5)
-
-# Fluid flow along z should be very small
-if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
- print("Flow field:\t\t" + passed())
-else:
- print("Flow field:\t\t" + failed())
- print(numpy.min(py.v_f))
- print(numpy.mean(py.v_f))
- print(numpy.max(py.v_f))
- raise Exception("Failed")
-
-print('''# Neumann bottom, Dirichlet top BC.
-# Gravity, pressure gradients => transient flow''')
-orig = sphere.sim("neumann", fluid = True)
-orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
-orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
-orig.initFluid(mu = 8.9e-4)
-#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
-#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 0.001)
-orig.initTemporal(total = 0.5, file_dt = 0.05, dt = 1.0e-4)
-py = sphere.sim(sid = orig.sid, fluid = True)
-orig.g[2] = -10.0
-orig.c_v[0] = 0.1
-orig.c_a[0] = 0.0
-orig.bc_bot[0] = 1 # No-flow BC at bottom (Neumann)
-#orig.run(dry=True)
-orig.run(verbose=False)
-#orig.run(device=2)
-orig.writeVTKall()
-py.readlast(verbose = False)
-#ideal_grad_p_z = numpy.linspace(
-# orig.p_f[0,0,0] + orig.L[2]*orig.rho_f*numpy.abs(orig.g[2]),
-# orig.p_f[0,0,-1], orig.num[2])
-ideal_grad_p_z = numpy.linspace(
- orig.p_f[0,0,0] + (orig.L[2]-orig.L[2]/orig.num[2])*orig.rho_f*numpy.abs(orig.g[2]),
- orig.p_f[0,0,-1], orig.num[2])
-compareNumpyArraysClose(ideal_grad_p_z, py.p_f[0,0,:],
- "Pressure gradient:\t", tolerance=1.0)
-
-# Fluid flow along z should be very small
-if ((numpy.abs(py.v_f[:,:,:,2]) < 1.0e-4).all()):
- print("Flow field:\t\t" + passed())
-else:
- print("Flow field:\t\t" + failed())
- raise Exception("Failed")
-
-#orig.cleanup()
(DIR) diff --git a/tests/cfd_tests_neumann-c_v=0.1.py b/tests/cfd_tests_neumann-c_v=0.1.py
t@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+from pytestutils import *
+
+import sphere
+import sys
+import numpy
+import matplotlib.pyplot as plt
+
+print('### CFD tests - Dirichlet/Neumann BCs ###')
+
+print('''# Neumann bottom, Dirichlet top BC.
+# No gravity, no pressure gradients => no flow''')
+orig = sphere.sim("neumann", fluid = True)
+cleanup(orig)
+orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
+orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
+orig.initFluid(mu = 8.9e-4)
+#orig.initFluid(mu = 0.0)
+orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
+orig.c_v[0] = 0.1
+#orig.c_phi[0] = 0.1
+py = sphere.sim(sid = orig.sid, fluid = True)
+orig.bc_bot[0] = 1 # No-flow BC at bottom (Neumann)
+#orig.run(dry=True)
+orig.run(verbose=False)
+#orig.run(device=2)
+#orig.writeVTKall()
+py.readlast(verbose = False)
+ones = numpy.ones((orig.num))
+py.readlast(verbose = False)
+compareNumpyArraysClose(ones, py.p_f, "Conservation of pressure:",
+ tolerance = 1.0e-5)
+
+# Fluid flow along z should be very small
+if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
+ print("Flow field:\t\t" + passed())
+else:
+ print("Flow field:\t\t" + failed())
+ print(numpy.min(py.v_f))
+ print(numpy.mean(py.v_f))
+ print(numpy.max(py.v_f))
+ raise Exception("Failed")
+
+print('''# Neumann bottom, Dirichlet top BC.
+# Gravity, pressure gradients => transient flow''')
+orig = sphere.sim("neumann", fluid = True)
+orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
+orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
+orig.initFluid(mu = 8.9e-4)
+#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
+#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 0.001)
+orig.initTemporal(total = 0.5, file_dt = 0.05, dt = 1.0e-4)
+py = sphere.sim(sid = orig.sid, fluid = True)
+orig.g[2] = -10.0
+orig.c_v[0] = 0.1
+orig.bc_bot[0] = 1 # No-flow BC at bottom (Neumann)
+#orig.run(dry=True)
+orig.run(verbose=False)
+#orig.run(device=2)
+orig.writeVTKall()
+py.readlast(verbose = False)
+#ideal_grad_p_z = numpy.linspace(
+# orig.p_f[0,0,0] + orig.L[2]*orig.rho_f*numpy.abs(orig.g[2]),
+# orig.p_f[0,0,-1], orig.num[2])
+ideal_grad_p_z = numpy.linspace(
+ orig.p_f[0,0,0] + (orig.L[2]-orig.L[2]/orig.num[2])*orig.rho_f*numpy.abs(orig.g[2]),
+ orig.p_f[0,0,-1], orig.num[2])
+compareNumpyArraysClose(ideal_grad_p_z, py.p_f[0,0,:],
+ "Pressure gradient:\t", tolerance=1.0)
+
+# Fluid flow along z should be very small
+if ((numpy.abs(py.v_f[:,:,:,2]) < 1.0e-4).all()):
+ print("Flow field:\t\t" + passed())
+else:
+ print("Flow field:\t\t" + failed())
+ raise Exception("Failed")
+
+#orig.cleanup()