tdebugging neumann test - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit 07ad1b9a07bcb030f82bfe8813fd7f4e8ad800eb
DIR parent e2b07236328876b4eb81c92d4e1beb3ed3646206
HTML Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 21 Oct 2014 14:28:11 +0200
debugging neumann test
Diffstat:
M python/sphere.py | 46 +++++++++++++++++++++++++++++--
M src/debug.h | 2 +-
M src/device.cu | 18 ++++++++++++++++++
M src/navierstokes.cuh | 23 +++++++++++++++++------
M tests/cfd_tests_neumann.py | 12 ++++++++++--
5 files changed, 90 insertions(+), 11 deletions(-)
---
DIR diff --git a/python/sphere.py b/python/sphere.py
t@@ -2545,6 +2545,48 @@ class sim:
self.mu_ws[0] = 0.0
self.mu_wd[0] = 0.0
+ def largestFluidTimeStep(self, safety=0.01):
+ '''
+ Finds and returns the largest time step in the fluid phase by von
+ Neumann and Courant-Friedrichs-Lewy analysis given the current
+ velocities. This ensures stability in the diffusive and advective parts
+ of the momentum equation.
+
+ The value of the time step decreases with increasing fluid viscosity
+ (`self.mu`), and increases with fluid cell size (`self.L/self.num`)
+
+ and fluid velocities (`self.v_f`)
+
+ :param safety: Safety factor which is multiplied to the largest time
+ step.
+ :type safety: float
+ :returns: The largest timestep stable for the current fluid state.
+ :return type: float
+ '''
+
+ dx_min = numpy.min(self.L/self.num)
+ dt_min_von_neumann = 0.5*dx_min**2/self.mu[0]
+
+ # Normalized velocities
+ v_norm = numpy.empty(self.num[0]*self.num[1]*self.num[2])
+ idx = 0
+ for x in numpy.arange(self.num[0]):
+ for y in numpy.arange(self.num[1]):
+ for z in numpy.arange(self.num[2]):
+ v_norm[idx] = numpy.sqrt(self.v_f[x,y,z,:].dot(\
+ self.v_f[x,y,z,:]))
+ idx += 1
+
+ # Advection term. This term has to be reevaluated during the
+ # computations, as the fluid velocity changes.
+ v_max = numpy.amax(v_norm)
+ if v_max == 0:
+ v_max = 1.0e-7
+
+ dt_min_cfl = dx_min/v_max
+
+ return numpy.min([dt_min_von_neumann, dt_min_cfl])*safety
+
def initTemporal(self, total,
current = 0.0,
file_dt = 0.05,
t@@ -2604,7 +2646,7 @@ class sim:
self.L[2]/self.num[2]))
# Diffusion term
- if (self.mu[0]*self.time_dt[0]/(dx*dx) > 0.5):
+ if (self.mu[0]*self.time_dt[0]/(dx**2) > 0.5):
raise Exception("Error: The time step is too large to ensure "
+ "stability in the diffusive term of the fluid "
+ "momentum equation.")
t@@ -2617,7 +2659,7 @@ class sim:
for z in numpy.arange(self.num[2]):
v_norm[idx] = numpy.sqrt(self.v_f[x,y,z,:].dot(\
self.v_f[x,y,z,:]))
- idx = idx + 1
+ idx += 1
# Advection term. This term has to be reevaluated during the
# computations, as the fluid velocity changes.
DIR diff --git a/src/debug.h b/src/debug.h
t@@ -47,7 +47,7 @@ const int conv_log_interval = 1;
#define REPORT_V_C_COMPONENTS
// Enable reporting of forcing function terms to stdout
-#define REPORT_FORCING_TERMS
+//#define REPORT_FORCING_TERMS
// Choose solver model (see Zhou et al. 2010 "Discrete particle simulation of
// particle-fluid flow: model formulations and their applicability", table. 1.
DIR diff --git a/src/device.cu b/src/device.cu
t@@ -1650,6 +1650,24 @@ __host__ void DEM::startTime()
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
&t_updateNSvelocityPressure);
checkForCudaErrorsIter("Post updateNSvelocity", iter);
+
+ setNSghostNodesFace<Float>
+ <<<dimGridFluidFace, dimBlockFluidFace>>>(
+ dev_ns_v_p_x,
+ dev_ns_v_p_y,
+ dev_ns_v_p_z,
+ ns.bc_bot, ns.bc_top);
+ cudaThreadSynchronize();
+ checkForCudaErrorsIter(
+ "Post setNSghostNodesFace(dev_ns_v)", iter);
+
+ interpolateFaceToCenter<<<dimGridFluid, dimBlockFluid>>>(
+ dev_ns_v_x,
+ dev_ns_v_y,
+ dev_ns_v_z,
+ dev_ns_v);
+ cudaThreadSynchronize();
+ checkForCudaErrorsIter("Post interpolateFaceToCenter", iter);
} // end iter % ns.dem == 0
} // end navierstokes == 1
DIR diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2954,7 +2954,8 @@ __global__ void updateNSvelocity(
const unsigned int cellidx = vidx(x,y,z);
// Check that we are not outside the fluid grid
- if (x <= nx && y <= ny && z <= nz) {
+ //if (x <= nx && y <= ny && z <= nz) {
+ if (x < nx && y < ny && z < nz) {
// Read values
__syncthreads();
t@@ -2995,11 +2996,6 @@ __global__ void updateNSvelocity(
//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;
#endif
-
-#ifdef REPORT_V_C_COMPONENTS
- printf("[%d,%d,%d] v_c = %f\t%f\t%f\n",
- v.x-v_p.x, v.y-v_p.y, v.z-v_p.z);
-#endif
// Print values for debugging
/* if (z == 0) {
Float e_up = dev_ns_epsilon[idx(x,y,z+1)];
t@@ -3039,6 +3035,21 @@ __global__ void updateNSvelocity(
if (z >= wall0_iz)
v = MAKE_FLOAT3(0.0, 0.0, 0.0);
+#ifdef REPORT_V_C_COMPONENTS
+ printf("[%d,%d,%d]\n"
+ "\tv_p = % f\t% f\t% f\n"
+ "\tv_c = % f\t% f\t% f\n"
+ "\tv = % f\t% f\t% f\n"
+ "\tgrad(epsilon) = % f\t% f\t% f\n"
+ "\tphi = % f\t% f\t% f\n",
+ x, y, z,
+ v_p.x, v_p.y, v_p.z,
+ v.x-v_p.x, v.y-v_p.y, v.z-v_p.z,
+ v.x, v.y, v.z,
+ grad_epsilon.x, grad_epsilon.y, grad_epsilon.z,
+ phi.x, phi.y, phi.z);
+#endif
+
// Check the advection term using the Courant-Friedrichs-Lewy condition
__syncthreads();
if (v.x*ndem*devC_dt/dx
DIR diff --git a/tests/cfd_tests_neumann.py b/tests/cfd_tests_neumann.py
t@@ -45,16 +45,24 @@ print('''# Neumann bottom, Dirichlet top BC.
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.g[2] = -10.0
orig.initFluid(mu = 8.9e-4)
#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
-orig.initTemporal(total = 1.0e-3, file_dt = 1.0e-4, dt = 1.0e-4)
+#orig.initTemporal(total = 1.0e-2, file_dt = 1.0e-4, dt = 1.0e-4)
+orig.initTemporal(total = 2.0e-4, file_dt = 1.0e-4, dt = 1.0e-4)
+#print(orig.largestFluidTimeStep())
+#orig.initTemporal(total = orig.largestFluidTimeStep()*10.0,
+ #file_dt = orig.largestFluidTimeStep(),
+ #dt = orig.largestFluidTimeStep())
py = sphere.sim(sid = orig.sid, fluid = True)
orig.g[2] = -10.0
orig.bc_bot[0] = 1 # No-flow BC at bottom (Neumann)
+orig.setTolerance(1.0e-12)
#orig.run(dry=True)
orig.run(verbose=False)
orig.writeVTKall()
py.readlast(verbose = False)
+print(py.v_f)
#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])
t@@ -72,4 +80,4 @@ else:
print("Flow field:\t\t" + failed())
raise Exception("Failed")
-orig.cleanup()
+#orig.cleanup()