tdebugging neumann test - sphere - GPU-based 3D discrete element method algorit… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit 07ad1b9a07bcb030f82bfe8813fd7f4e8ad800eb | |
parent e2b07236328876b4eb81c92d4e1beb3ed3646206 | |
Author: Anders Damsgaard <[email protected]> | |
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(-) | |
--- | |
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. | |
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. | |
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 | |
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 | |
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() |