tstoring fluid interaction force components in output - 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 6ce8fd65468c3a399a4b131a3bce345a4cc3fdac
(DIR) parent 096a799126d8262917e89ae78079f57e350fbaeb
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 15 Aug 2014 12:01:20 +0200
storing fluid interaction force components in output
Diffstat:
M python/sphere.py | 139 +++++++++++++++++++++++++++++--
M src/constants.h | 2 +-
M src/datatypes.h | 4 ++++
M src/device.cu | 6 +++++-
M src/file_io.cpp | 44 ++++++++++++++++++++++++++++++-
M src/navierstokes.cpp | 8 ++++++++
M src/navierstokes.cuh | 29 ++++++++++++++++++++++++++++-
M src/sphere.h | 4 ++++
M tests/io_tests_fluid.py | 2 +-
9 files changed, 224 insertions(+), 14 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.04
+VERSION=1.05
class sim:
'''
t@@ -203,15 +203,15 @@ class sim:
# Wall normals
self.w_n = numpy.zeros((self.nw, self.nd), dtype=numpy.float64)
- if (self.nw >= 1):
+ if self.nw >= 1:
self.w_n[0,2] = -1.0
- if (self.nw >= 2):
+ if self.nw >= 2:
self.w_n[1,0] = -1.0
- if (self.nw >= 3):
+ if self.nw >= 3:
self.w_n[2,0] = 1.0
- if (self.nw >= 4):
+ if self.nw >= 4:
self.w_n[3,1] = -1.0
- if (self.nw >= 5):
+ if self.nw >= 5:
self.w_n[4,1] = 1.0
# Wall positions on the axes that are parallel to the wall normal [m]
t@@ -270,7 +270,7 @@ class sim:
# Simulate fluid? True: Yes, False: no
self.fluid = fluid
- if (self.fluid == True):
+ if self.fluid == True:
# Fluid dynamic viscosity [N/(m/s)]
self.mu = numpy.zeros(1, dtype=numpy.float64)
t@@ -337,6 +337,12 @@ class sim:
# Permeability scaling factor
self.c_grad_p = numpy.ones(1, dtype=numpy.float64)
+ ## Interaction forces
+ self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+
# Particle color marker
self.color = numpy.zeros(self.np, dtype=numpy.int32)
t@@ -600,6 +606,18 @@ class sim:
elif (self.c_grad_p != other.c_grad_p):
print(85)
return(85)
+ elif (self.f_d != other.f_d).any():
+ print(86)
+ return(86)
+ elif (self.f_p != other.f_p).any():
+ print(87)
+ return(87)
+ elif (self.f_v != other.f_v).any():
+ print(88)
+ return(88)
+ elif (self.f_sum != other.f_sum).any():
+ print(89)
+ return(89)
if ((self.color != other.color)).any():
print(90)
t@@ -673,6 +691,10 @@ class sim:
self.ev = numpy.append(self.ev, ev)
self.p = numpy.append(self.p, p)
self.color = numpy.append(self.color, color)
+ self.f_d = numpy.append(self.f_d, [numpy.zeros(3)], axis=0)
+ self.f_p = numpy.append(self.f_p, [numpy.zeros(3)], axis=0)
+ self.f_v = numpy.append(self.f_v, [numpy.zeros(3)], axis=0)
+ self.f_sum = numpy.append(self.f_sum, [numpy.zeros(3)], axis=0)
def deleteParticle(self, i):
'''
t@@ -699,6 +721,10 @@ class sim:
self.ev = numpy.delete(self.ev, i)
self.p = numpy.delete(self.p, i)
self.color = numpy.delete(self.color, i)
+ self.f_d = numpy.delete(self.f_d, i, axis=0)
+ self.f_p = numpy.delete(self.f_p, i, axis=0)
+ self.f_v = numpy.delete(self.f_v, i, axis=0)
+ self.f_sum = numpy.delete(self.f_sum, i, axis=0)
def deleteAllParticles(self):
'''
t@@ -720,6 +746,10 @@ class sim:
self.ev = numpy.zeros(self.np, dtype=numpy.float64)
self.p = numpy.zeros(self.np, dtype=numpy.float64)
self.color = numpy.zeros(self.np, dtype=numpy.int32)
+ self.f_d = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
def readbin(self, targetbin, verbose = True, bonds = True, devsmod = True,
esysparticle = False):
t@@ -747,7 +777,7 @@ class sim:
'''
fh = None
- try :
+ try:
if (verbose == True):
print("Input file: {0}".format(targetbin))
fh = open(targetbin, "rb")
t@@ -819,7 +849,7 @@ class sim:
self.torque = numpy.fromfile(fh, dtype=numpy.float64,\
count=self.np*self.nd).reshape(self.np, self.nd)
- if (esysparticle == True):
+ if esysparticle == True:
return
# Per-particle single-value parameters
t@@ -976,6 +1006,38 @@ class sim:
self.c_phi = numpy.ones(1, dtype=numpy.float64)
self.c_grad_p = numpy.ones(1, dtype=numpy.float64)
+ if self.version >= 1.05:
+ self.f_d = numpy.empty_like(self.x)
+ self.f_p = numpy.empty_like(self.x)
+ self.f_v = numpy.empty_like(self.x)
+ self.f_sum = numpy.empty_like(self.x)
+
+ for i in range(self.np[0]):
+ self.f_d[i,:] = \
+ numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nd)
+ for i in range(self.np[0]):
+ self.f_p[i,:] = \
+ numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nd)
+ for i in range(self.np[0]):
+ self.f_v[i,:] = \
+ numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nd)
+ for i in range(self.np[0]):
+ self.f_sum[i,:] = \
+ numpy.fromfile(fh, dtype=numpy.float64,
+ count=self.nd)
+ else:
+ self.f_d = numpy.zeros((self.np, self.nd),
+ dtype=numpy.float64)
+ self.f_p = numpy.zeros((self.np, self.nd),
+ dtype=numpy.float64)
+ self.f_v = numpy.zeros((self.np, self.nd),
+ dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd),
+ dtype=numpy.float64)
+
if (self.version >= 1.02):
self.color =\
numpy.fromfile(fh, dtype=numpy.int32, count=self.np)
t@@ -1132,6 +1194,15 @@ class sim:
fh.write(self.c_phi.astype(numpy.float64))
fh.write(self.c_grad_p.astype(numpy.float64))
+ for i in numpy.arange(self.np):
+ fh.write(self.f_d[i,:].astype(numpy.float64))
+ for i in numpy.arange(self.np):
+ fh.write(self.f_p[i,:].astype(numpy.float64))
+ for i in numpy.arange(self.np):
+ fh.write(self.f_v[i,:].astype(numpy.float64))
+ for i in numpy.arange(self.np):
+ fh.write(self.f_sum[i,:].astype(numpy.float64))
+
fh.write(self.color.astype(numpy.int32))
finally:
t@@ -1309,6 +1380,51 @@ class sim:
fh.write('\n')
fh.write(' </DataArray>\n')
+ if self.fluid == True:
+ # Fluid interaction force
+ fh.write(' <DataArray type="Float32" '
+ + 'Name="Fluid force total" '
+ + 'NumberOfComponents="3" format="ascii">\n')
+ fh.write(' ')
+ for i in range(self.np):
+ fh.write('%f %f %f ' % \
+ (self.f_sum[i,0], self.f_sum[i,1], self.f_sum[i,2]))
+ fh.write('\n')
+ fh.write(' </DataArray>\n')
+
+ # Fluid drag force
+ fh.write(' <DataArray type="Float32" '
+ + 'Name="Fluid drag force" '
+ + 'NumberOfComponents="3" format="ascii">\n')
+ fh.write(' ')
+ for i in range(self.np):
+ fh.write('%f %f %f ' % \
+ (self.f_d[i,0], self.f_d[i,1], self.f_d[i,2]))
+ fh.write('\n')
+ fh.write(' </DataArray>\n')
+
+ # Fluid pressure force
+ fh.write(' <DataArray type="Float32" '
+ + 'Name="Fluid pressure force" '
+ + 'NumberOfComponents="3" format="ascii">\n')
+ fh.write(' ')
+ for i in range(self.np):
+ fh.write('%f %f %f ' % \
+ (self.f_p[i,0], self.f_p[i,1], self.f_p[i,2]))
+ fh.write('\n')
+ fh.write(' </DataArray>\n')
+
+ # Fluid viscous force
+ fh.write(' <DataArray type="Float32" '
+ + 'Name="Fluid viscous force" '
+ + 'NumberOfComponents="3" format="ascii">\n')
+ fh.write(' ')
+ for i in range(self.np):
+ fh.write('%f %f %f ' % \
+ (self.f_v[i,0], self.f_v[i,1], self.f_v[i,2]))
+ fh.write('\n')
+ fh.write(' </DataArray>\n')
+
# fixvel
fh.write(' <DataArray type="Float32" Name="FixedVel" '
+ 'format="ascii">\n')
t@@ -2508,6 +2624,11 @@ class sim:
self.c_phi = numpy.ones(1, dtype=numpy.float64)
self.c_grad_p = 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)
+ self.f_v = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+ self.f_sum = numpy.zeros((self.np, self.nd), dtype=numpy.float64)
+
def setFluidBottomNoFlow(self):
'''
Set the lower boundary of the fluid domain to follow the no-flow
(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.04;
+const double VERSION = 1.05;
// Max. number of contacts per particle
//const int NC = 16;
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -141,6 +141,10 @@ struct NavierStokes {
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
+ Float4* f_d; // Drag force on particles
+ Float4* f_p; // Pressure force on particles
+ Float4* f_v; // Viscous force on particles
+ Float4* f_sum; // Viscous force on particles
};
// Image structure
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1120,7 +1120,11 @@ __host__ void DEM::startTime()
dev_ns_div_tau_y,
dev_ns_div_tau_z,
dev_ns_f_pf,
- dev_force);
+ dev_force,
+ dev_ns_f_d,
+ dev_ns_f_p,
+ dev_ns_f_v,
+ dev_ns_f_sum);
cudaThreadSynchronize();
checkForCudaErrorsIter("Post findInteractionForce", iter);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -12,7 +12,7 @@
// Get the address of the first byte of an object's representation
// See Stroustrup (2008) p. 388
template<class T>
-char* as_bytes(T& i) // treat a T as a sequence of bytes
+char* as_bytes(T& i) // treat a T as a sequence of bytes
{
// get the address of the first byte of memory used
// to store the object
t@@ -297,6 +297,27 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.c_phi), sizeof(Float));
ifs.read(as_bytes(ns.c_grad_p), sizeof(Float));
+ for (i = 0; i<np; ++i) {
+ ifs.read(as_bytes(ns.f_d[i].x), sizeof(Float));
+ ifs.read(as_bytes(ns.f_d[i].y), sizeof(Float));
+ ifs.read(as_bytes(ns.f_d[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ifs.read(as_bytes(ns.f_p[i].x), sizeof(Float));
+ ifs.read(as_bytes(ns.f_p[i].y), sizeof(Float));
+ ifs.read(as_bytes(ns.f_p[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ifs.read(as_bytes(ns.f_v[i].x), sizeof(Float));
+ ifs.read(as_bytes(ns.f_v[i].y), sizeof(Float));
+ ifs.read(as_bytes(ns.f_v[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ifs.read(as_bytes(ns.f_sum[i].x), sizeof(Float));
+ ifs.read(as_bytes(ns.f_sum[i].y), sizeof(Float));
+ ifs.read(as_bytes(ns.f_sum[i].z), sizeof(Float));
+ }
+
if (verbose == 1)
cout << "Done" << std::endl;
}
t@@ -508,6 +529,27 @@ void DEM::writebin(const char *target)
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) {
+ ofs.write(as_bytes(ns.f_d[i].x), sizeof(Float));
+ ofs.write(as_bytes(ns.f_d[i].y), sizeof(Float));
+ ofs.write(as_bytes(ns.f_d[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ofs.write(as_bytes(ns.f_p[i].x), sizeof(Float));
+ ofs.write(as_bytes(ns.f_p[i].y), sizeof(Float));
+ ofs.write(as_bytes(ns.f_p[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ofs.write(as_bytes(ns.f_v[i].x), sizeof(Float));
+ ofs.write(as_bytes(ns.f_v[i].y), sizeof(Float));
+ ofs.write(as_bytes(ns.f_v[i].z), sizeof(Float));
+ }
+ for (i = 0; i<np; ++i) {
+ ofs.write(as_bytes(ns.f_sum[i].x), sizeof(Float));
+ ofs.write(as_bytes(ns.f_sum[i].y), sizeof(Float));
+ ofs.write(as_bytes(ns.f_sum[i].z), sizeof(Float));
+ }
}
for (i = 0; i<np; ++i)
(DIR) diff --git a/src/navierstokes.cpp b/src/navierstokes.cpp
t@@ -37,6 +37,10 @@ void DEM::initNSmem()
ns.norm = new Float[ncells]; // normalized residual of epsilon
ns.epsilon = new Float[ncells]; // normalized residual of epsilon
ns.epsilon_new = new Float[ncells]; // normalized residual of epsilon
+ ns.f_d = new Float4[np]; // drag force on particles
+ ns.f_p = new Float4[np]; // pressure force on particles
+ ns.f_v = new Float4[np]; // viscous force on particles
+ ns.f_sum = new Float4[np]; // sum of fluid forces on particles
}
unsigned int DEM::NScells()
t@@ -73,6 +77,10 @@ void DEM::freeNSmem()
delete[] ns.norm;
delete[] ns.epsilon;
delete[] ns.epsilon_new;
+ delete[] ns.f_d;
+ delete[] ns.f_p;
+ delete[] ns.f_v;
+ delete[] ns.f_sum;
}
// 3D index to 1D index
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -97,6 +97,10 @@ void DEM::initNSmemDev(void)
cudaMalloc((void**)&dev_ns_div_phi_vi_v, memSizeF*3); // div(phi*vi*v)
//cudaMalloc((void**)&dev_ns_div_phi_tau, memSizeF*3); // div(phi*tau)
cudaMalloc((void**)&dev_ns_f_pf, sizeof(Float3)*np); // particle fluid force
+ cudaMalloc((void**)&dev_ns_f_d, sizeof(Float4)*np); // drag force
+ cudaMalloc((void**)&dev_ns_f_p, sizeof(Float4)*np); // pressure force
+ cudaMalloc((void**)&dev_ns_f_v, sizeof(Float4)*np); // viscous force
+ cudaMalloc((void**)&dev_ns_f_sum, sizeof(Float4)*np); // sum of int. forces
checkForCudaErrors("End of initNSmemDev");
}
t@@ -137,6 +141,10 @@ void DEM::freeNSmemDev()
cudaFree(dev_ns_div_tau_y);
cudaFree(dev_ns_div_tau_z);
cudaFree(dev_ns_f_pf);
+ cudaFree(dev_ns_f_d);
+ cudaFree(dev_ns_f_p);
+ cudaFree(dev_ns_f_v);
+ cudaFree(dev_ns_f_sum);
}
// Transfer to device
t@@ -184,6 +192,11 @@ void DEM::transferNSfromGlobalDeviceMemory(int statusmsg)
cudaMemcpy(ns.phi, dev_ns_phi, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(ns.dphi, dev_ns_dphi, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(ns.norm, dev_ns_norm, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.f_d, dev_ns_f_d, sizeof(Float4)*np, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.f_p, dev_ns_f_p, sizeof(Float4)*np, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.f_v, dev_ns_f_v, sizeof(Float4)*np, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.f_sum, dev_ns_f_sum, sizeof(Float4)*np,
+ cudaMemcpyDeviceToHost);
checkForCudaErrors("End of transferNSfromGlobalDeviceMemory", 0);
if (verbose == 1 && statusmsg == 1)
t@@ -3041,7 +3054,11 @@ __global__ void findInteractionForce(
const Float* __restrict__ dev_ns_div_tau_y,// in
const Float* __restrict__ dev_ns_div_tau_z,// in
Float3* __restrict__ dev_ns_f_pf, // out
- Float4* __restrict__ dev_force) // out
+ Float4* __restrict__ dev_force, // out
+ Float4* __restrict__ dev_ns_f_d, // out
+ Float4* __restrict__ dev_ns_f_p, // out
+ Float4* __restrict__ dev_ns_f_v, // out
+ Float4* __restrict__ dev_ns_f_sum) // out
{
unsigned int i = threadIdx.x + blockIdx.x*blockDim.x; // Particle index
t@@ -3141,6 +3158,7 @@ __global__ void findInteractionForce(
checkFiniteFloat3("f_pf", i_x, i_y, i_z, f_pf);
#endif
+ __syncthreads();
#ifdef SET_1
dev_ns_f_pf[i] = f_pf;
#endif
t@@ -3149,7 +3167,16 @@ __global__ void findInteractionForce(
dev_ns_f_pf[i] = f_d;
#endif
+ __syncthreads();
dev_force[i] += MAKE_FLOAT4(f_pf.x, f_pf.y, f_pf.z, 0.0);
+ dev_ns_f_d[i] = MAKE_FLOAT4(f_d.x, f_d.y, f_d.z, 0.0);
+ dev_ns_f_p[i] = MAKE_FLOAT4(f_p.x, f_p.y, f_p.z, 0.0);
+ dev_ns_f_v[i] = MAKE_FLOAT4(f_v.x, f_v.y, f_v.z, 0.0);
+ dev_ns_f_sum[i] = MAKE_FLOAT4(
+ f_d.x + f_p.x + f_v.x,
+ f_d.y + f_p.y + f_v.y,
+ f_d.z + f_p.z + f_v.z,
+ 0.0);
}
}
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -213,6 +213,10 @@ class DEM {
Float* dev_ns_div_tau_y; // div(tau) on y-face
Float* dev_ns_div_tau_z; // div(tau) on z-face
Float3* dev_ns_f_pf; // Interaction force on particles
+ Float4* dev_ns_f_d; // Drag force on particles
+ Float4* dev_ns_f_p; // Pressure force on particles
+ Float4* dev_ns_f_v; // Viscous force on particles
+ Float4* dev_ns_f_sum; // Total force on particles
// Helper device arrays, input
unsigned int** hdev_gridParticleIndex;
(DIR) diff --git a/tests/io_tests_fluid.py b/tests/io_tests_fluid.py
t@@ -22,7 +22,7 @@ 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()
#orig.run(dry=True)
#orig.run(verbose=True, hideinputfile=False, cudamemcheck=True)
cpp = sphere.sim(fluid=True)