tset velocity BCs on predicted values as well - 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 0aaf70618d39783d3f0a86d1fa21bad0a7ebe8a9
(DIR) parent e4516a4cde0fcf1ea1f1eee8ad48c172e681f2f3
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 22 Oct 2014 14:19:52 +0200
set velocity BCs on predicted values as well
Diffstat:
M src/debug.h | 2 +-
M src/device.cu | 10 +++++++++-
M src/navierstokes.cuh | 27 +++++++++++++++++++++------
M tests/cfd_tests_neumann.py | 2 --
4 files changed, 31 insertions(+), 10 deletions(-)
---
(DIR) diff --git a/src/debug.h b/src/debug.h
t@@ -22,7 +22,7 @@ const int write_res_log = 0;
// Report epsilon values during Jacobi iterations to stdout
#define REPORT_EPSILON
-#define REPORT_MORE_EPSILON
+//#define REPORT_MORE_EPSILON
// Report the number of iterations it took before convergence to logfile
// 'output/<sid>-conv.dat'
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1176,6 +1176,12 @@ __host__ void DEM::startTime()
cudaThreadSynchronize();
checkForCudaErrorsIter("Post setNSepsilonTop", iter);
+#if defined(REPORT_EPSILON) || defined(REPORT_V_P_COMPONENTS) || defined(REPORT_V_C_COMPONENTS)
+ std::cout
+ << "\n\n@@@@@@ TIME STEP " << iter << " @@@"
+ << std::endl;
+#endif
+
// find cell containing top wall
if (walls.nw > 0 && walls.wmode[0] == 1) {
wall0_iz = walls.nx->w/(grid.L[2]/grid.num[2]);
t@@ -1187,9 +1193,9 @@ __host__ void DEM::startTime()
dp_dz);
cudaThreadSynchronize();
checkForCudaErrorsIter("Post setNSepsilonAtTopWall", iter);
+
#ifdef REPORT_EPSILON
std::cout
- << "\n@@@@@@ TIME STEP " << iter << " @@@@@@"
<< "\n###### EPSILON setNSepsilonAtTopWall "
<< "######" << std::endl;
transferNSepsilonFromGlobalDeviceMemory();
t@@ -1313,6 +1319,7 @@ __host__ void DEM::startTime()
ns.beta,
dev_ns_F_pf,
ns.ndem,
+ wall0_iz,
ns.c_grad_p,
dev_ns_v_p_x,
dev_ns_v_p_y,
t@@ -1643,6 +1650,7 @@ __host__ void DEM::startTime()
ns.ndem,
ns.c_grad_p,
wall0_iz,
+ iter,
dev_ns_v_x,
dev_ns_v_y,
dev_ns_v_z);
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2288,7 +2288,8 @@ __global__ void findPredNSvelocities(
const Float beta, // in
const Float3* __restrict__ dev_ns_F_pf, // in
const unsigned int ndem, // in
- const Float __restrict__ c_grad_p, // in
+ const unsigned int wall0_iz, // in
+ const Float c_grad_p, // 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@@ -2449,6 +2450,17 @@ __global__ void findPredNSvelocities(
v_p.z = 0.0;
}*/
+ //if ((z == 0 && bc_bot == 1) || (z == nz-1 && bc_top == 1))
+ if ((z == 0 && bc_bot == 1) || (z > nz-1 && bc_top == 1))
+ v_p.z = 0.0;
+
+ //if ((z == 0 && bc_bot == 2) || (z == nz-1 && bc_top == 2))
+ if ((z == 0 && bc_bot == 2) || (z > nz-1 && bc_top == 2))
+ v_p = MAKE_FLOAT3(0.0, 0.0, 0.0);
+
+ // Set velocities to zero above the dynamic wall
+ if (z >= wall0_iz)
+ v_p = MAKE_FLOAT3(0.0, 0.0, 0.0);
#ifdef REPORT_V_P_COMPONENTS
// Report velocity components to stdout for debugging
t@@ -2934,6 +2946,7 @@ __global__ void updateNSvelocity(
const unsigned int ndem, // in
const Float c_grad_p, // in
const unsigned int wall0_iz, // in
+ const unsigned int iter, // in
Float* __restrict__ dev_ns_v_x, // out
Float* __restrict__ dev_ns_v_y, // out
Float* __restrict__ dev_ns_v_z) // out
t@@ -2992,13 +3005,15 @@ __global__ void updateNSvelocity(
#ifdef SET_1
//Float3 v = v_p - ndem*devC_dt/(devC_params.rho_f*phi)*grad_epsilon;
//Float3 v = v_p - ndem*devC_dt*c_grad_p/devC_params.rho_f*grad_epsilon;
- Float3 v =
- v_p - ndem*devC_dt*c_grad_p/(phi*devC_params.rho_f)*grad_epsilon;
+ Float3 v = v_p
+ -ndem*devC_dt*c_grad_p/(phi*devC_params.rho_f)*grad_epsilon;
#endif
#ifdef SET_2
//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;
+ Float3 v = v_p
+ -ndem*devC_dt*c_grad_p/devC_params.rho_f*grad_epsilon;
#endif
+
// Print values for debugging
/* if (z == 0) {
Float e_up = dev_ns_epsilon[idx(x,y,z+1)];
t@@ -3396,7 +3411,7 @@ __global__ void interpolateCenterToFace(
const Float z_val = amean(center.z, zn.z);
__syncthreads();
- printf("c2f [%d,%d,%d] = %f,%f,%f\n", x,y,z, x_val, y_val, z_val);
+ //printf("c2f [%d,%d,%d] = %f,%f,%f\n", x,y,z, x_val, y_val, z_val);
dev_out_x[faceidx] = x_val;
dev_out_y[faceidx] = y_val;
dev_out_z[faceidx] = z_val;
t@@ -3438,7 +3453,7 @@ __global__ void interpolateFaceToCenter(
amean(z_n, z_p));
__syncthreads();
- printf("f2c [%d,%d,%d] = %f, %f, %f\n", x,y,z, val.x, val.y, val.z);
+ //printf("f2c [%d,%d,%d] = %f, %f, %f\n", x,y,z, val.x, val.y, val.z);
dev_out[idx(x,y,z)] = val;
}
}
(DIR) diff --git a/tests/cfd_tests_neumann.py b/tests/cfd_tests_neumann.py
t@@ -51,7 +51,6 @@ orig.initFluid(mu = 8.9e-4)
#orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
#orig.initTemporal(total = 1.0e-2, file_dt = 1.0e-4, dt = 1.0e-4)
orig.initTemporal(total = 1.0e-3, file_dt = 1.0e-4, dt = 1.0e-4)
-orig.v_f[:,:,:,2] = 0.0
#print(orig.largestFluidTimeStep())
#orig.initTemporal(total = orig.largestFluidTimeStep()*10.0,
#file_dt = orig.largestFluidTimeStep(),
t@@ -59,7 +58,6 @@ orig.v_f[:,:,:,2] = 0.0
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()