tuse central differences for porosity change - 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 ef0c5c667eeab770245e87f4d4664aacee5c9081
(DIR) parent 5cf55424f3bd21f92b22998c12cccb14a42c807b
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Wed, 26 Nov 2014 11:31:44 +0100
use central differences for porosity change
Diffstat:
M src/darcy.cuh | 49 ++++++++++++++++++-------------
M src/device.cu | 8 ++++----
2 files changed, 32 insertions(+), 25 deletions(-)
---
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -277,10 +277,10 @@ __global__ void findDarcyPorosities(
const unsigned int np, // in
const Float c_phi, // in
Float* __restrict__ dev_darcy_phi, // in + out
- //Float* __restrict__ dev_darcy_dphi) // in + out
- Float* __restrict__ dev_darcy_dphi, // in + out
+ Float* __restrict__ dev_darcy_dphi) // in + out
+ //Float* __restrict__ dev_darcy_dphi, // in + out
//Float* __restrict__ dev_darcy_dphi) // out
- Float* __restrict__ dev_darcy_div_v_p) // out
+ //Float* __restrict__ dev_darcy_div_v_p) // out
{
// 3D thread index
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
t@@ -434,8 +434,14 @@ __global__ void findDarcyPorosities(
//phi = fmin(0.99, fmax(0.01, void_volume/cell_volume));
//phi = void_volume/cell_volume;
- /*Float dphi = phi - phi_0;*/
- Float dphi = phi_new - phi;
+ // Backwards Euler
+ //Float dphi = phi - phi_0;
+
+ // Forwards Euler
+ //Float dphi = phi_new - phi;
+
+ // Central difference
+ Float dphi = 0.5*(phi_new - phi_0);
/*if (iteration == 0)
dphi = 0.0;*/
t@@ -454,7 +460,7 @@ __global__ void findDarcyPorosities(
dev_darcy_dphi[cellidx] += dphi*c_phi;
//dev_darcy_vp_avg[cellidx] = v_avg;
//dev_darcy_d_avg[cellidx] = d_avg;
- dev_darcy_div_v_p[cellidx] = dot_epsilon_kk*c_phi;
+ //dev_darcy_div_v_p[cellidx] = dot_epsilon_kk*c_phi;
/*printf("\n%d,%d,%d: findDarcyPorosities\n"
"\tphi = %f\n"
t@@ -475,7 +481,7 @@ __global__ void findDarcyPorosities(
#ifdef CHECK_FLUID_FINITE
(void)checkFiniteFloat("phi", x, y, z, phi);
(void)checkFiniteFloat("dphi", x, y, z, dphi);
- (void)checkFiniteFloat("dot_epsilon_kk", x, y, z, dot_epsilon_kk);
+ //(void)checkFiniteFloat("dot_epsilon_kk", x, y, z, dot_epsilon_kk);
//(void)checkFiniteFloat3("v_avg", x, y, z, v_avg);
//(void)checkFiniteFloat("d_avg", x, y, z, d_avg);
#endif
t@@ -857,7 +863,7 @@ __global__ void firstDarcySolution(
const Float* __restrict__ dev_darcy_k, // in
const Float* __restrict__ dev_darcy_phi, // in
const Float* __restrict__ dev_darcy_dphi, // in
- const Float* __restrict__ dev_darcy_div_v_p, // in
+ //const Float* __restrict__ dev_darcy_div_v_p, // in
const Float3* __restrict__ dev_darcy_grad_k, // in
const Float beta_f, // in
const Float mu, // in
t@@ -903,8 +909,8 @@ __global__ void firstDarcySolution(
const Float k = dev_darcy_k[cellidx];
const Float3 grad_k = dev_darcy_grad_k[cellidx];
const Float phi = dev_darcy_phi[cellidx];
- //const Float dphi = dev_darcy_dphi[cellidx];
- const Float div_v_p = dev_darcy_div_v_p[cellidx];
+ const Float dphi = dev_darcy_dphi[cellidx];
+ //const Float div_v_p = dev_darcy_div_v_p[cellidx];
const Float p_xn = dev_darcy_p[d_idx(x-1,y,z)];
const Float p = dev_darcy_p[cellidx];
t@@ -949,8 +955,8 @@ __global__ void firstDarcySolution(
Float dp_expl =
+ (ndem*devC_dt)/(beta_f*phi*mu)*(k*laplace_p + dot(grad_k, grad_p))
- - div_v_p/(beta_f*phi);
- //- dphi/(beta_f*phi*(1.0 - phi));
+ //- div_v_p/(beta_f*phi);
+ - dphi/(beta_f*phi*(1.0 - phi));
// Dirichlet BC at dynamic top wall. wall0_iz will be larger than the
// grid if the wall isn't dynamic
t@@ -974,8 +980,8 @@ __global__ void firstDarcySolution(
"grad_k = %e, %e, %e\n"
"dp_diff = %e\n"
"dp_forc = %e\n"
- "div_v_p = %e\n"
- //"dphi = %e\n"
+ //"div_v_p = %e\n"
+ "dphi = %e\n"
//"dphi/dt = %e\n"
,
x,y,z,
t@@ -988,8 +994,9 @@ __global__ void firstDarcySolution(
grad_p.x, grad_p.y, grad_p.z,
grad_k.x, grad_k.y, grad_k.z,
dp_diff, dp_forc,
- div_v_p//,
- //dphi, dphi/(ndem*devC_dt)
+ //div_v_p//,
+ dphi//,
+ //dphi/(ndem*devC_dt)
);
#endif
t@@ -1013,7 +1020,7 @@ __global__ void updateDarcySolution(
const Float* __restrict__ dev_darcy_k, // in
const Float* __restrict__ dev_darcy_phi, // in
const Float* __restrict__ dev_darcy_dphi, // in
- const Float* __restrict__ dev_darcy_div_v_p, // in
+ //const Float* __restrict__ dev_darcy_div_v_p, // in
const Float3* __restrict__ dev_darcy_grad_k, // in
const Float beta_f, // in
const Float mu, // in
t@@ -1060,8 +1067,8 @@ __global__ void updateDarcySolution(
const Float k = dev_darcy_k[cellidx];
const Float3 grad_k = dev_darcy_grad_k[cellidx];
const Float phi = dev_darcy_phi[cellidx];
- //const Float dphi = dev_darcy_dphi[cellidx];
- const Float div_v_p = dev_darcy_div_v_p[cellidx];
+ const Float dphi = dev_darcy_dphi[cellidx];
+ //const Float div_v_p = dev_darcy_div_v_p[cellidx];
const Float p_old = dev_darcy_p_old[cellidx];
const Float dp_expl = dev_darcy_dp_expl[cellidx];
t@@ -1110,8 +1117,8 @@ __global__ void updateDarcySolution(
//Float p_new = p_old
Float dp_impl =
+ (ndem*devC_dt)/(beta_f*phi*mu)*(k*laplace_p + dot(grad_k, grad_p))
- - div_v_p/(beta_f*phi);
- //- dphi/(beta_f*phi*(1.0 - phi));
+ //- div_v_p/(beta_f*phi);
+ - dphi/(beta_f*phi*(1.0 - phi));
// Dirichlet BC at dynamic top wall. wall0_iz will be larger than the
// grid if the wall isn't dynamic
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -1796,8 +1796,8 @@ __host__ void DEM::startTime()
np,
darcy.c_phi,
dev_darcy_phi,
- dev_darcy_dphi,
- dev_darcy_div_v_p);
+ dev_darcy_dphi);//,
+ //dev_darcy_div_v_p);
cudaThreadSynchronize();
if (PROFILING == 1)
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
t@@ -1976,7 +1976,7 @@ __host__ void DEM::startTime()
dev_darcy_k,
dev_darcy_phi,
dev_darcy_dphi,
- dev_darcy_div_v_p,
+ //dev_darcy_div_v_p,
dev_darcy_grad_k,
darcy.beta_f,
darcy.mu,
t@@ -2003,7 +2003,7 @@ __host__ void DEM::startTime()
dev_darcy_k,
dev_darcy_phi,
dev_darcy_dphi,
- dev_darcy_div_v_p,
+ //dev_darcy_div_v_p,
dev_darcy_grad_k,
darcy.beta_f,
darcy.mu,