trevised velocity prediction, velocity correction and forcing function - 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 4c2467acf0646b9a49314578eba3b28a7fdc9dff
(DIR) parent d1db5adde78dc6c20a26447c353001e1cf35e2e7
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Sun, 12 Oct 2014 16:11:58 +0200
revised velocity prediction, velocity correction and forcing function
Diffstat:
M src/navierstokes.cuh | 27 ++++++++++++++++++---------
M tests/fluid_particle_interaction.py | 50 +++++++++++++++++++++++++++++--
2 files changed, 65 insertions(+), 12 deletions(-)
---
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -2391,9 +2391,10 @@ __global__ void findPredNSvelocities(
const Float3 grad_p = MAKE_FLOAT3(
(p - p_xn)/dx,
(p - p_yn)/dy,
- (p - p_zn)/dz) * c_grad_p;
+ (p - p_zn)/dz);
#ifdef SET_1
- pressure_term = -beta*dt/(rho*phi)*grad_p;
+ //pressure_term = -beta*dt/(rho*phi)*grad_p;
+ pressure_term = -beta*c_grad_p*dt/rho*grad_p;
#endif
#ifdef SET_2
pressure_term = -beta*dt/rho*grad_p;
t@@ -2408,7 +2409,8 @@ __global__ void findPredNSvelocities(
// Determine the predicted velocity
#ifdef SET_1
const Float3 interaction_term = -dt/(rho*phi)*f_i;
- const Float3 diffusion_term = dt/(rho*phi)*div_tau;
+ //const Float3 diffusion_term = dt/(rho*phi)*div_tau;
+ const Float3 diffusion_term = dt/rho*div_tau;
const Float3 gravity_term = MAKE_FLOAT3(
devC_params.g[0], devC_params.g[1], devC_params.g[2])*dt;
const Float3 porosity_term = -1.0*v*dphi/phi;
t@@ -2548,10 +2550,16 @@ __global__ void findNSforcing(
//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 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);
+
+ const Float t1 = 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);
+ devC_params.rho_f*dot(v_p, grad_phi)/(phi*dt*c_grad_p);
+ const Float t4 = devC_params.rho_f*dphi/(dt*dt*c_grad_p*phi);
#endif
#ifdef SET_2
t@@ -2560,8 +2568,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; // 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@@ -2972,7 +2980,8 @@ __global__ void updateNSvelocity(
// Find new velocity
#ifdef SET_1
- Float3 v = v_p - ndem*devC_dt/(devC_params.rho_f*phi)*grad_epsilon;
+ //Float3 v = v_p - ndem*devC_dt/(devC_params.rho_f*phi)*grad_epsilon;
+ Float3 v = v_p - ndem*devC_dt/devC_params.rho_f*grad_epsilon;
#endif
#ifdef SET_2
Float3 v = v_p - ndem*devC_dt/devC_params.rho_f*grad_epsilon;
(DIR) diff --git a/tests/fluid_particle_interaction.py b/tests/fluid_particle_interaction.py
t@@ -19,7 +19,7 @@ sim.p_f[:,:,-1] = 1.0
sim.addParticle([0.5, 0.5, 0.5], 0.05)
sim.initTemporal(total=0.01, file_dt=0.001)
-sim.run(verbose=False)
+sim.run(device=1, verbose=False)
#sim.run(dry=True)
#sim.run(cudamemcheck=True)
#sim.writeVTKall()
t@@ -43,7 +43,7 @@ sim.addParticle([0.5, 0.5, 0.25], 0.05)
sim.initTemporal(total=0.0001, file_dt=0.00001)
-sim.run(verbose=False)
+sim.run(device=1, verbose=False)
#sim.writeVTKall()
sim.readlast()
t@@ -52,4 +52,48 @@ test(sim.vel[0,0] > 0.0, 'Particle 0 velocity:')
test(sim.vel[1,0] > 0.0, 'Particle 1 velocity:')
test(sim.vel[2,0] > 0.0, 'Particle 2 velocity:')
-sim.cleanup()
+
+'''
+print('# Test 3: Test pressure gradient force, c = 0.1')
+sim.p_f[:,:,0] = 10.0
+sim.p_f[:,:,-1] = 1.0
+sim.addParticle([0.5, 0.5, 0.5], 0.05)
+sim.initTemporal(total=0.01, file_dt=0.001)
+sim.c_grad_p[0] = 0.1
+
+sim.run(device=1, verbose=False)
+#sim.run(dry=True)
+#sim.run(cudamemcheck=True)
+#sim.writeVTKall()
+
+sim.readlast()
+test(sim.vel[0,2] > 0.0, 'Particle velocity:')
+
+
+
+# Sidewards gravity, homogenous pressure, Neumann boundaries.
+# Fluid should flow towards +x and drag particles in the same direction
+print('# Test 4: Test fluid drag force, c = 0.1')
+sim.initFluid()
+sim.zeroKinematics()
+sim.g[0] = 10.0
+sim.c_grad_p[0] = 0.1
+
+sim.deleteParticle(0)
+sim.addParticle([0.5, 0.5, 0.75], 0.05)
+sim.addParticle([0.5, 0.5, 0.50], 0.05)
+sim.addParticle([0.5, 0.5, 0.25], 0.05)
+
+sim.initTemporal(total=0.0001, file_dt=0.00001)
+
+sim.run(device=1, verbose=False)
+#sim.writeVTKall()
+
+sim.readlast()
+test((sim.v_f[:,:,:,0] > 0.0).all(), 'Fluid velocity:')
+test(sim.vel[0,0] > 0.0, 'Particle 0 velocity:')
+test(sim.vel[1,0] > 0.0, 'Particle 1 velocity:')
+test(sim.vel[2,0] > 0.0, 'Particle 2 velocity:')
+'''
+
+#sim.cleanup()