twrite cell face velocities to files - 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 e07f904c139196389023592dc85028e8cd657a3b
(DIR) parent 0c32bfa91146574058a5c91800efebde8e957ec0
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Thu, 26 Jun 2014 11:17:19 +0200
write cell face velocities to files
Diffstat:
M python/shortening.py | 6 +++---
M src/datatypes.h | 6 +++---
M src/file_io.cpp | 9 ++++++---
M src/navierstokes.cpp | 12 ++++++------
M src/navierstokes.cuh | 27 ++++++++++++++++++++-------
M tests/cfd_inclined.py | 3 ++-
M tests/cfd_tests_neumann.py | 2 ++
7 files changed, 42 insertions(+), 23 deletions(-)
---
(DIR) diff --git a/python/shortening.py b/python/shortening.py
t@@ -108,7 +108,7 @@ sim = sphere.sim('shortening-relaxation', nw=1)
sim.readlast()
sim.sid = 'shortening'
sim.cleanup()
-sim.initTemporal(current=0.0, total=20.0, file_dt = 0.01, epsilon=0.01)
+sim.initTemporal(current=0.0, total=20.0, file_dt = 0.01, epsilon=0.07)
# set colors again
y_min = numpy.min(sim.x[:,1])
t@@ -165,8 +165,8 @@ sim.gamma_wt[0] = 0.0
# Particle parameters
sim.mu_s[0] = 0.5
sim.mu_d[0] = 0.5
-sim.gamma_n[0] = 0.0
-sim.gamma_t[0] = 0.0
+sim.gamma_n[0] = 100.0
+sim.gamma_t[0] = 100.0
# push down upper wall
compressional_strain = 0.5
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -114,9 +114,9 @@ struct NavierStokes {
Float dx, dy, dz; // Cell length in each dim
Float* p; // Cell hydraulic pressures
Float3* v; // Cell fluid velocity
- //Float* v_x; // Fluid velocity in staggered grid
- //Float* v_y; // Fluid velocity in staggered grid
- //Float* v_z; // Fluid velocity in staggered grid
+ Float* v_x; // Fluid velocity in staggered grid
+ Float* v_y; // Fluid velocity in staggered grid
+ Float* v_z; // Fluid velocity in staggered grid
//Float3* v_p; // Predicted fluid velocity
//Float* v_p_x; // Predicted fluid velocity in staggered grid
//Float* v_p_y; // Predicted fluid velocity in staggered grid
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -468,9 +468,12 @@ void DEM::writebin(const char *target)
for (y=0; y<ns.ny; y++) {
for (x=0; x<ns.nx; x++) {
i = idx(x,y,z);
- ofs.write(as_bytes(ns.v[i].x), sizeof(Float));
- ofs.write(as_bytes(ns.v[i].y), sizeof(Float));
- ofs.write(as_bytes(ns.v[i].z), sizeof(Float));
+ //ofs.write(as_bytes(ns.v[i].x), sizeof(Float));
+ //ofs.write(as_bytes(ns.v[i].y), sizeof(Float));
+ //ofs.write(as_bytes(ns.v[i].z), sizeof(Float));
+ ofs.write(as_bytes(ns.v_x[i]), sizeof(Float));
+ ofs.write(as_bytes(ns.v_y[i]), sizeof(Float));
+ ofs.write(as_bytes(ns.v_z[i]), sizeof(Float));
ofs.write(as_bytes(ns.p[i]), sizeof(Float));
ofs.write(as_bytes(ns.phi[i]), sizeof(Float));
ofs.write(as_bytes(ns.dphi[i]), sizeof(Float));
(DIR) diff --git a/src/navierstokes.cpp b/src/navierstokes.cpp
t@@ -25,9 +25,9 @@ void DEM::initNSmem()
ns.p = new Float[ncells]; // hydraulic pressure
ns.v = new Float3[ncells]; // hydraulic velocity
- //ns.v_x = new Float[ncells_st]; // hydraulic velocity in staggered grid
- //ns.v_y = new Float[ncells_st]; // hydraulic velocity in staggered grid
- //ns.v_z = new Float[ncells_st]; // hydraulic velocity in staggered grid
+ ns.v_x = new Float[ncells_st]; // hydraulic velocity in staggered grid
+ ns.v_y = new Float[ncells_st]; // hydraulic velocity in staggered grid
+ ns.v_z = new Float[ncells_st]; // hydraulic velocity in staggered grid
//ns.v_p = new Float3[ncells]; // predicted hydraulic velocity
//ns.v_p_x = new Float[ncells_st]; // pred. hydraulic velocity in st. grid
//ns.v_p_y = new Float[ncells_st]; // pred. hydraulic velocity in st. grid
t@@ -61,9 +61,9 @@ void DEM::freeNSmem()
{
delete[] ns.p;
delete[] ns.v;
- //delete[] ns.v_x;
- //delete[] ns.v_y;
- //delete[] ns.v_z;
+ delete[] ns.v_x;
+ delete[] ns.v_y;
+ delete[] ns.v_z;
//delete[] ns.v_p;
//delete[] ns.v_p_x;
//delete[] ns.v_p_y;
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -177,6 +177,9 @@ void DEM::transferNSfromGlobalDeviceMemory(int statusmsg)
cudaMemcpy(ns.p, dev_ns_p, memSizeF, cudaMemcpyDeviceToHost);
checkForCudaErrors("In transferNSfromGlobalDeviceMemory, dev_ns_p", 0);
cudaMemcpy(ns.v, dev_ns_v, memSizeF*3, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.v_x, dev_ns_v_x, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.v_y, dev_ns_v_y, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(ns.v_z, dev_ns_v_z, memSizeF, cudaMemcpyDeviceToHost);
//cudaMemcpy(ns.v_p, dev_ns_v_p, memSizeF*3, cudaMemcpyDeviceToHost);
cudaMemcpy(ns.phi, dev_ns_phi, memSizeF, cudaMemcpyDeviceToHost);
cudaMemcpy(ns.dphi, dev_ns_dphi, memSizeF, cudaMemcpyDeviceToHost);
t@@ -714,7 +717,7 @@ __global__ void setNSghostNodesFace(
//dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Neumann free slip -z
dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann free slip -z
dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Neumann free slip -z
- dev_scalarfield_z[vidx(x,y,-1)] = 0.0; // Neumann free slip -z
+ dev_scalarfield_z[vidx(x,y,-1)] = 0.0; // Neumann free slip -z
}
if (z == 0 && bc_bot == 2) {
//dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann no slip -z
t@@ -2285,8 +2288,23 @@ __global__ void findPredNSvelocities(
+ porosity_term
+ advection_term;
+ //// Neumann BCs
+
+ // Free slip
+ if ((z == 0 && bc_bot == 1) || (z == nz-1 && bc_top == 1))
+ v_p.z = v.z;
+
+ // No slip
+ if ((z == 0 && bc_bot == 2) || (z == nz-1 && bc_top == 2)) {
+ v_p.x = 0.0;
+ v_p.y = 0.0;
+ v_p.z = 0.0;
+ }
+
+
#ifdef REPORT_V_P_COMPONENTS
// Report velocity components to stdout for debugging
+ if (z==0)
printf("\n[%d,%d,%d]"
"\tv_p = %+e %+e %+e\n"
"\tpres = %+e %+e %+e\n"
t@@ -2294,7 +2312,7 @@ __global__ void findPredNSvelocities(
"\tdiff = %+e %+e %+e\n"
"\tgrav = %+e %+e %+e\n"
"\tporos = %+e %+e %+e\n"
- "\tadv = %+e %+e %+e\n"
+ "\tadv = %+e %+e %+e\n",
x, y, z,
v_p.x, v_p.y, v_p.z,
pressure_term.x, pressure_term.y, pressure_term.z,
t@@ -2305,11 +2323,6 @@ __global__ void findPredNSvelocities(
advection_term.x, advection_term.y, advection_term.z);
#endif
- // Enforce Neumann BC if specified
- if ((z == 0 && bc_bot == 1) || (z == nz-1 && bc_top == 1))
- v_p.z = v.z;
- //v_p.z = 0.0;
-
// Save the predicted velocity
__syncthreads();
dev_ns_v_p_x[fidx] = v_p.x;
(DIR) diff --git a/tests/cfd_inclined.py b/tests/cfd_inclined.py
t@@ -9,7 +9,8 @@ orig.defineWorldBoundaries([0.3, 0.3, 0.3], dx = 0.06)
orig.initFluid(mu=8.9e-4) # inviscid "fluids" (mu=0) won't work!
#orig.initTemporal(total = 0.5, file_dt = 0.05, dt = 1.0e-4)
orig.initTemporal(total = 1.0e-0, file_dt = 1.0e-1, dt = 1.0e-3)
-orig.bc_bot[0] = 2 # No-flow, no slip BC at bottom (Neumann)
+orig.bc_bot[0] = 1 # No-flow, free slip BC at bottom (Neumann)
+#orig.bc_bot[0] = 2 # No-flow, no slip BC at bottom (Neumann)
#orig.bc_top[0] = 1 # No-flow, free slip BC at top (Neumann)
angle = 10.0 # slab inclination in degrees
(DIR) diff --git a/tests/cfd_tests_neumann.py b/tests/cfd_tests_neumann.py
t@@ -10,6 +10,7 @@ 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)
t@@ -37,6 +38,7 @@ else:
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''')