tuse reformulated forcing function, remove c_grad_p from pressure gradient force - 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 d1db5adde78dc6c20a26447c353001e1cf35e2e7
(DIR) parent 5fb2209cc878db151655af8f0681a723063a5957
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 10 Oct 2014 10:15:40 +0200
use reformulated forcing function, remove c_grad_p from pressure gradient force
Diffstat:
M src/device.cu | 2 +-
M src/navierstokes.cuh | 17 +++++++++++------
2 files changed, 12 insertions(+), 7 deletions(-)
---
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1125,7 +1125,7 @@ __host__ void DEM::startTime()
dev_ns_div_tau_x,
dev_ns_div_tau_y,
dev_ns_div_tau_z,
- ns.c_grad_p,
+ //ns.c_grad_p,
dev_ns_f_pf,
dev_force,
dev_ns_f_d,
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2545,9 +2545,13 @@ __global__ void findNSforcing(
// Find forcing function terms
#ifdef SET_1
- const Float t1 = phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
- const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(c_grad_p*dt);
- const Float t4 = dphi*devC_params.rho_f/(c_grad_p*dt*dt);
+ //const Float t1 = phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
+ //const Float t2 = devC_params.rho_f*dot(v_p, grad_phi)/(c_grad_p*dt);
+ //const Float t4 = dphi*devC_params.rho_f/(c_grad_p*dt*dt);
+ const Float t1 = phi*phi*devC_params.rho_f*div_v_p/(c_grad_p*dt);
+ const Float t2 =
+ devC_params.rho_f*phi*dot(v_p, grad_phi)/(c_grad_p*dt);
+ const Float t4 = dphi*devC_params.rho_f*phi/(c_grad_p*dt*dt);
#endif
#ifdef SET_2
t@@ -2556,7 +2560,8 @@ __global__ void findNSforcing(
const Float t4 = dphi*devC_params.rho_f/(dt*dt*phi);
#endif
f1 = t1 + t2 + t4;
- f2 = grad_phi/phi; // t3/grad(epsilon)
+ //f2 = grad_phi/phi; // t3/grad(epsilon)
+ f2 = grad_phi; // t3/grad(epsilon)
#ifdef REPORT_FORCING_TERMS
// Report values terms in the forcing function for debugging
t@@ -3133,7 +3138,7 @@ __global__ void findInteractionForce(
const Float* __restrict__ dev_ns_div_tau_x,// in
const Float* __restrict__ dev_ns_div_tau_y,// in
const Float* __restrict__ dev_ns_div_tau_z,// in
- const Float c_grad_p, // in
+ //const Float c_grad_p, // in
Float3* __restrict__ dev_ns_f_pf, // out
Float4* __restrict__ dev_force, // out
Float4* __restrict__ dev_ns_f_d, // out
t@@ -3213,7 +3218,7 @@ __global__ void findInteractionForce(
// Pressure gradient force
const Float3 f_p =
- -c_grad_p*gradient(dev_ns_p, i_x, i_y, i_z, dx, dy, dz)*V_p;
+ -1.0*gradient(dev_ns_p, i_x, i_y, i_z, dx, dy, dz)*V_p;
// Viscous force
const Float3 f_v = -1.0*div_tau*V_p;