tfind porosities before pressure force - sphere - GPU-based 3D discrete element… | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit e040779c930d67d2c56fd5f843f02be9773cffe4 | |
parent 1e778bdb6456a0e499e4a03063dddefb3378d933 | |
Author: Anders Damsgaard <[email protected]> | |
Date: Mon, 23 Mar 2015 12:38:35 +0100 | |
find porosities before pressure force | |
Diffstat: | |
M src/darcy.cuh | 28 ++++++++++++++++------------ | |
M src/device.cu | 40 ++++++++++++++++-------------… | |
M tests/fluid_particle_interaction_d… | 2 +- | |
3 files changed, 37 insertions(+), 33 deletions(-) | |
--- | |
diff --git a/src/darcy.cuh b/src/darcy.cuh | |
t@@ -340,7 +340,7 @@ __global__ void findDarcyPorositiesLinear( | |
const Float dy = devC_grid.L[1]/ny; | |
const Float dz = devC_grid.L[2]/nz; | |
- Float void_volume = dx*dy*dz; // current void volume | |
+ Float solid_volume = 0.0; | |
Float4 xr; // particle pos. and radius | |
// check that we are not outside the fluid grid | |
t@@ -447,8 +447,7 @@ __global__ void findDarcyPorositiesLinear( | |
s = weightDist(dist, dx, dy, dz); | |
vol_p = sphereVolume(xr.w); | |
- // Subtract particle volume times weight | |
- void_volume -= s*vol_p; | |
+ solid_volume += s*vol_p; | |
// Add particle contribution to cell face | |
// nodes of component-wise velocity | |
t@@ -509,7 +508,7 @@ __global__ void findDarcyPorositiesLinear( | |
// Make sure that the porosity is in the interval [0.0;1.0] | |
//phi = fmin(0.9, fmax(0.1, void_volume/(dx*dy*dz))); | |
- phi = fmin(1.0, fmax(0.01, void_volume/(dx*dy*dz))); | |
+ phi = fmin(1.0, fmax(0.01, 1.0 - solid_volume/(dx*dy*dz))); | |
// Determine particle velocity divergence | |
/*const Float div_v_p = | |
t@@ -518,11 +517,11 @@ __global__ void findDarcyPorositiesLinear( | |
(v_p_zp - v_p_zn)/dz;*/ | |
const Float div_v_p = | |
(xp_num/fmax(1.e-12, xp_denum) | |
- - xn_num/fmax(1.e-12, xn_denum)) /dx + | |
+ - xn_num/fmax(1.e-12, xn_denum))/dx + | |
(yp_num/fmax(1.e-12, yp_denum) | |
- - yn_num/fmax(1.e-12, yn_denum)) /dy + | |
+ - yn_num/fmax(1.e-12, yn_denum))/dy + | |
(zp_num/fmax(1.e-12, zp_denum) | |
- - zn_num/fmax(1.e-12, zn_denum)) /dz; | |
+ - zn_num/fmax(1.e-12, zn_denum))/dz; | |
// Save porosity and porosity change | |
__syncthreads(); | |
t@@ -530,10 +529,12 @@ __global__ void findDarcyPorositiesLinear( | |
dev_darcy_phi[cellidx] = phi*c_phi; | |
dev_darcy_div_v_p[cellidx] = div_v_p; | |
- //if (phi < 1.0 || div_v_p != 0.0) | |
- if (div_v_p >= 1.0e-12 || div_v_p <= -1.0e-12) | |
+ if (phi < 1.0 || div_v_p != 0.0) | |
+ //if (div_v_p >= 1.0e-12 || div_v_p <= -1.0e-12) | |
printf("\n%d,%d,%d: findDarcyPorosities\n" | |
"\tphi = %f\n" | |
+ "\tsol_vol = %f\n" | |
+ "\tvol_p = %f\n" | |
"\tX = %.2e, %.2e, %.2e\n" | |
"\txr = %.2e, %.2e, %.2e\n" | |
"\tdiv_v_p = %.2e\n" | |
t@@ -543,6 +544,8 @@ __global__ void findDarcyPorositiesLinear( | |
//"\tv_p_z = %.2e, %.2e\n" | |
, x,y,z, | |
phi, | |
+ solid_volume, | |
+ vol_p, | |
X.x, X.y, X.z, | |
xr.x, xr.y, xr.z, | |
div_v_p, | |
t@@ -1046,7 +1049,7 @@ __global__ void findDarcyPressureForceLinear( | |
// find pressure gradient force plus buoyancy force. | |
// buoyancy force = weight of displaced fluid | |
// f_b = -rho_f*V*g | |
- Float3 f_p = -1.0*grad_p*v/(1.0-phi) | |
+ Float3 f_p = -1.0*grad_p*v/(1.0 - phi) | |
- rho_f*v*MAKE_FLOAT3( | |
devC_params.g[0], | |
devC_params.g[1], | |
t@@ -1057,14 +1060,15 @@ __global__ void findDarcyPressureForceLinear( | |
if (i_z >= wall0_iz) | |
f_p.z = 0.0; | |
- /*if (length(f_p) > 1.0e-12) | |
+ //if (length(f_p) > 1.0e-12) | |
printf("%d,%d,%d findPF:\n" | |
- //"\tphi = %f\n" | |
+ "\tphi = %f\n" | |
"\tx = %f, %f, %f\n" | |
"\tX = %f, %f, %f\n" | |
"\tgrad_p = %.2e, %.2e, %.2e\n" | |
"\tf_p = %.2e, %.2e, %.2e\n", | |
i_x, i_y, i_z, | |
+ phi, | |
x3.x, x3.y, x3.z, | |
X.x, X.y, X.z, | |
grad_p.x, grad_p.y, grad_p.z, | |
diff --git a/src/device.cu b/src/device.cu | |
t@@ -1827,6 +1827,26 @@ __host__ void DEM::startTime() | |
checkForCudaErrorsIter("Post setDarcyGhostNodes(" | |
"dev_darcy_grad_p)", iter); | |
+ if (PROFILING == 1) | |
+ startTimer(&kernel_tic); | |
+ findDarcyPorositiesLinear<<<dimGridFluid, dimBlockFluid>>>( | |
+ dev_cellStart, | |
+ dev_cellEnd, | |
+ dev_x_sorted, | |
+ dev_vel_sorted, | |
+ iter, | |
+ darcy.ndem, | |
+ np, | |
+ darcy.c_phi, | |
+ dev_darcy_phi, | |
+ dev_darcy_dphi, | |
+ dev_darcy_div_v_p); | |
+ cudaThreadSynchronize(); | |
+ if (PROFILING == 1) | |
+ stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, | |
+ &t_findDarcyPorosities); | |
+ checkForCudaErrorsIter("Post findDarcyPorosities", iter); | |
+ | |
/*findDarcyPressureForce<<<dimGrid, dimBlock>>>( | |
dev_x, | |
dev_darcy_p, | |
t@@ -1852,8 +1872,6 @@ __host__ void DEM::startTime() | |
if ((iter % darcy.ndem) == 0) { | |
- if (PROFILING == 1) | |
- startTimer(&kernel_tic); | |
/*findDarcyPorosities<<<dimGridFluid, dimBlockFluid>>>( | |
dev_cellStart, | |
dev_cellEnd, | |
t@@ -1865,24 +1883,6 @@ __host__ void DEM::startTime() | |
darcy.c_phi, | |
dev_darcy_phi, | |
dev_darcy_dphi);*/ | |
- findDarcyPorositiesLinear<<<dimGridFluid, dimBlockFluid>>>( | |
- dev_cellStart, | |
- dev_cellEnd, | |
- dev_x_sorted, | |
- dev_vel_sorted, | |
- iter, | |
- darcy.ndem, | |
- np, | |
- darcy.c_phi, | |
- dev_darcy_phi, | |
- dev_darcy_dphi, | |
- dev_darcy_div_v_p); | |
- cudaThreadSynchronize(); | |
- if (PROFILING == 1) | |
- stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, | |
- &t_findDarcyPorosities); | |
- checkForCudaErrorsIter("Post findDarcyPorosities", iter); | |
- | |
// Modulate the pressures at the upper boundary cells | |
if ((darcy.p_mod_A > 1.0e-5 || darcy.p_mod_A < -1.0e-5) && | |
darcy.p_mod_f > 1.0e-7) { | |
diff --git a/tests/fluid_particle_interaction_darcy.py b/tests/fluid_particle_i… | |
t@@ -15,7 +15,7 @@ sim.initFluid(cfd_solver = 1) | |
print('# Test 1: Test pressure gradient force') | |
#sim.p_f[:,:,-1] = 1.0 | |
#sim.addParticle([0.5, 0.5, 0.5], 0.01) | |
-sim.addParticle([0.55, 0.55, 0.55], 0.03) | |
+sim.addParticle([0.55, 0.55, 0.55], 0.05) | |
sim.vel[0,2] = 1.0e-2 | |
sim.initTemporal(total=0.001, file_dt=0.0001) | |
sim.time_file_dt[0] = sim.time_dt[0]*10 |