timplemented correction for spatial advection of porosity - sphere - GPU-based … | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit f3f8a7ca3858ddfea3d16e2c20abb1b983478231 | |
parent 6613f7d7b2f2e3a9e149fadc5296e484f617496a | |
Author: Anders Damsgaard <[email protected]> | |
Date: Tue, 31 Mar 2015 13:22:32 +0200 | |
implemented correction for spatial advection of porosity | |
Diffstat: | |
M python/halfshear-darcy-combined.py | 22 ++++++++++++---------- | |
M python/halfshear-darcy-rate-depend… | 32 +++++++++++++++++++++------… | |
M python/halfshear-darcy-sigma0=8000… | 6 +++--- | |
M src/darcy.cuh | 69 ++++++++++++++++++++++++-----… | |
M src/device.cu | 5 ++++- | |
M src/sphere.h | 7 ++++--- | |
6 files changed, 99 insertions(+), 42 deletions(-) | |
--- | |
diff --git a/python/halfshear-darcy-combined.py b/python/halfshear-darcy-combin… | |
t@@ -18,7 +18,8 @@ sid = 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-0… | |
outformat = 'pdf' | |
fluid = True | |
#threshold = 100.0 # [N] | |
-calculateforcechains = False | |
+calculateforcechains = True | |
+calculateforcechainhistory = False | |
legend_alpha=0.7 | |
linewidth=0.5 | |
t@@ -106,7 +107,7 @@ for i in numpy.arange(nsteps): | |
#''' | |
if calculateforcechains: | |
- if i > 0: | |
+ if i > 0 and calculateforcechainhistory: | |
loaded_contacts_prev = numpy.copy(loaded_contacts) | |
pairs_prev = numpy.copy(sim.pairs) | |
t@@ -118,17 +119,18 @@ for i in numpy.arange(nsteps): | |
sim.findCoordinationNumber() | |
coordinationnumber[i] = sim.findMeanCoordinationNumber() | |
- nfound = 0 | |
- if i > 0: | |
- for a in loaded_contacts[0]: | |
- for b in loaded_contacts_prev[0]: | |
- if (sim.pairs[:,a] == pairs_prev[:,b]).all(): | |
- nfound += 1; | |
+ if calculateforcechainhistory: | |
+ nfound = 0 | |
+ if i > 0: | |
+ for a in loaded_contacts[0]: | |
+ for b in loaded_contacts_prev[0]: | |
+ if (sim.pairs[:,a] == pairs_prev[:,b]).all(): | |
+ nfound += 1; | |
- nkept[i] = nfound | |
+ nkept[i] = nfound | |
+ print nfound | |
print coordinationnumber[i] | |
- print nfound | |
if calculateforcechains: | |
diff --git a/python/halfshear-darcy-rate-dependence.py b/python/halfshear-darcy… | |
t@@ -22,10 +22,14 @@ fluid = True | |
#threshold = 100.0 # [N] | |
-def creep_rheology(friction, n): | |
+def creep_rheology1(friction, n): | |
''' return strain rate from friction (tau/N) value ''' | |
return friction**n | |
+def creep_rheology2(friction, n, A): | |
+ ''' return strain rate from friction (tau/N) value ''' | |
+ return A*friction**n | |
+ | |
for sid in sids: | |
t@@ -79,12 +83,14 @@ for sid in sids: | |
#popt, pvoc = scipy.optimize.curve_fit( | |
#creep_rheology, tau/N, shearstrainrate) | |
popt, pvoc = scipy.optimize.curve_fit( | |
- creep_rheology, tau_nonzero[idxfit]/N_nonzero[idxfit], | |
+ creep_rheology1, tau_nonzero[idxfit]/N_nonzero[idxfit], | |
shearstrainrate_nonzero[idxfit]) | |
print '# Critical state' | |
print popt | |
print pvoc | |
n = popt[0] # stress exponent | |
+ #A = popt[1] # stress exponent | |
+ A = 1. | |
friction = tau_nonzero/N_nonzero | |
x_min = numpy.floor(numpy.min(friction)) | |
t@@ -93,27 +99,30 @@ for sid in sids: | |
numpy.linspace(numpy.min(tau_nonzero[idxfit]/N_nonzero[idxfit]), | |
numpy.max(tau_nonzero[idxfit]/N_nonzero[idxfit]), | |
100) | |
- strainrate_fit = friction_fit**n | |
+ strainrate_fit = A*friction_fit**n | |
### Consolidated state fit | |
- idxfit2 = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) & | |
- (shearstrainrate_nonzero < 0.1) & | |
- ((t_nonzero > 0.0) & (t_nonzero < 3.0))) | |
- #popt, pvoc = scipy.optimize.curve_fit( | |
- #creep_rheology, tau/N, shearstrainrate) | |
+ #idxfit2 = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) & | |
+ #(shearstrainrate_nonzero < 0.1) & | |
+ #((t_nonzero > 0.0) & (t_nonzero < 2.5))) | |
+ idxfit2 = numpy.nonzero((shearstrain_nonzero < 0.1) & | |
+ (tau_nonzero/N_nonzero < 0.38)) | |
+ ''' | |
popt2, pvoc2 = scipy.optimize.curve_fit( | |
- creep_rheology, tau_nonzero[idxfit2]/N_nonzero[idxfit2], | |
+ creep_rheology2, tau_nonzero[idxfit2]/N_nonzero[idxfit2], | |
shearstrainrate_nonzero[idxfit2]) | |
print '# Consolidated state' | |
print popt2 | |
print pvoc2 | |
n2 = popt2[0] # stress exponent | |
+ A2 = popt2[1] # prefactor | |
friction_fit2 =\ | |
numpy.linspace(numpy.min(tau_nonzero[idxfit2]/N_nonzero[idxfit2]), | |
numpy.max(tau_nonzero[idxfit2]/N_nonzero[idxfit2]), | |
100) | |
- strainrate_fit2 = friction_fit2**n2 | |
+ strainrate_fit2 = A2*friction_fit2**n2 | |
+ ''' | |
### Plotting | |
t@@ -182,6 +191,9 @@ for sid in sids: | |
fig = plt.figure(figsize=(3.5,2.5)) | |
ax1 = plt.subplot(111) | |
CS = ax1.scatter(friction, dilation[idx], | |
+ #tau_nonzero[idxfit2]/N_nonzero[idxfit2], | |
+ shearstrainrate_nonzero[idxfit2], | |
+ c=shearstrain_nonzero[idxfit2], linewidth=0.05, | |
c=shearstrain_nonzero, linewidth=0.05, | |
cmap=matplotlib.cm.get_cmap('afmhot')) | |
diff --git a/python/halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=1… | |
t@@ -1,3 +1,3 @@ | |
-3.210000000000000000e+03 3.285000000000000000e+03 3.140000000000000000e+03 2.8… | |
-0.000000000000000000e+00 3.101000000000000000e+03 3.097000000000000000e+03 2.8… | |
-4.956041666666666679e+00 5.008958333333333179e+00 4.964583333333333570e+00 4.9… | |
+3.210000000000000000e+03 3.285000000000000000e+03 3.140000000000000000e+03 2.8… | |
+6.903866600074103356e-310 1.014250921842801434e-315 1.158246314797981072e-315 … | |
+4.956041666666666679e+00 5.008958333333333179e+00 4.964583333333333570e+00 4.9… | |
diff --git a/src/darcy.cuh b/src/darcy.cuh | |
t@@ -28,7 +28,7 @@ void DEM::initDarcyMemDev(void) | |
cudaMalloc((void**)&dev_darcy_p, memSizeF); // hydraulic pressure | |
cudaMalloc((void**)&dev_darcy_p_new, memSizeF); // updated pressure | |
cudaMalloc((void**)&dev_darcy_v, memSizeF*3); // cell hydraulic velocity | |
- //cudaMalloc((void**)&dev_darcy_vp_avg, memSizeF*3); // avg. particle velo… | |
+ cudaMalloc((void**)&dev_darcy_vp_avg, memSizeF*3); // avg. particle veloci… | |
//cudaMalloc((void**)&dev_darcy_d_avg, memSizeF); // avg. particle diameter | |
cudaMalloc((void**)&dev_darcy_phi, memSizeF); // cell porosity | |
cudaMalloc((void**)&dev_darcy_dphi, memSizeF); // cell porosity change | |
t@@ -54,7 +54,7 @@ void DEM::freeDarcyMemDev() | |
cudaFree(dev_darcy_p); | |
cudaFree(dev_darcy_p_new); | |
cudaFree(dev_darcy_v); | |
- //cudaFree(dev_darcy_vp_avg); | |
+ cudaFree(dev_darcy_vp_avg); | |
//cudaFree(dev_darcy_d_avg); | |
cudaFree(dev_darcy_phi); | |
cudaFree(dev_darcy_dphi); | |
t@@ -323,7 +323,8 @@ __global__ void findDarcyPorositiesLinear( | |
const Float c_phi, // in | |
Float* __restrict__ dev_darcy_phi, // in + out | |
Float* __restrict__ dev_darcy_dphi, // in + out | |
- Float* __restrict__ dev_darcy_div_v_p) // out | |
+ Float* __restrict__ dev_darcy_div_v_p, // out | |
+ Float3* __restrict__ dev_darcy_vp_avg) // out | |
{ | |
// 3D thread index | |
const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x; | |
t@@ -359,6 +360,8 @@ __global__ void findDarcyPorositiesLinear( | |
Float s, vol_p; | |
Float phi = 1.00; | |
Float4 v; | |
+ Float3 vp_avg_num = MAKE_FLOAT3(0.0, 0.0, 0.0); | |
+ Float vp_avg_denum = 0.0; | |
//Float3 x3, v3; | |
//unsigned int n = 0; | |
t@@ -451,6 +454,11 @@ __global__ void findDarcyPorositiesLinear( | |
solid_volume += s*vol_p; | |
+ // Summation procedure for average velocity | |
+ vp_avg_num += | |
+ s*vol_p*MAKE_FLOAT3(v.x, v.y, v.z); | |
+ vp_avg_denum += s*vol_p; | |
+ | |
// Add particle contribution to cell face | |
// nodes of component-wise velocity | |
/*x3 += distmod; | |
t@@ -561,6 +569,7 @@ __global__ void findDarcyPorositiesLinear( | |
dev_darcy_phi[cellidx] = phi*c_phi; | |
dev_darcy_dphi[cellidx] = dphi*c_phi; | |
//dev_darcy_div_v_p[cellidx] = div_v_p; | |
+ dev_darcy_vp_avg[cellidx] = vp_avg_num/fmax(1.0e-16, vp_avg_denum); | |
//if (phi < 1.0 || div_v_p != 0.0) | |
//if (phi < 1.0) | |
t@@ -1459,6 +1468,7 @@ __global__ void firstDarcySolution( | |
const Float* __restrict__ dev_darcy_phi, // in | |
const Float* __restrict__ dev_darcy_dphi, // in | |
const Float* __restrict__ dev_darcy_div_v_p, // in | |
+ const Float3* __restrict__ dev_darcy_vp_avg, // in | |
const Float3* __restrict__ dev_darcy_grad_k, // in | |
const Float beta_f, // in | |
const Float mu, // in | |
t@@ -1490,11 +1500,18 @@ __global__ void firstDarcySolution( | |
// read values | |
__syncthreads(); | |
- 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 k = dev_darcy_k[cellidx]; | |
+ const Float3 grad_k = dev_darcy_grad_k[cellidx]; | |
+ const Float phi_xn = dev_darcy_phi[d_idx(x-1,y,z)]; | |
+ const Float phi = dev_darcy_phi[cellidx]; | |
+ const Float phi_xp = dev_darcy_phi[d_idx(x+1,y,z)]; | |
+ const Float phi_yn = dev_darcy_phi[d_idx(x,y-1,z)]; | |
+ const Float phi_yp = dev_darcy_phi[d_idx(x,y+1,z)]; | |
+ const Float phi_zn = dev_darcy_phi[d_idx(x,y,z-1)]; | |
+ const Float phi_zp = dev_darcy_phi[d_idx(x,y,z+1)]; | |
+ const Float dphi = dev_darcy_dphi[cellidx]; | |
//const Float div_v_p = dev_darcy_div_v_p[cellidx]; | |
+ const Float3 vp_avg = dev_darcy_vp_avg[cellidx]; | |
const Float p_xn = dev_darcy_p[d_idx(x-1,y,z)]; | |
const Float p = dev_darcy_p[cellidx]; | |
t@@ -1532,6 +1549,11 @@ __global__ void firstDarcySolution( | |
(p_yp - p_yn)/(dy + dy), | |
(p_zp - p_zn)/(dz + dz)); | |
+ const Float3 grad_phi = MAKE_FLOAT3( | |
+ (phi_xp - phi_xn)/(dx + dx), | |
+ (phi_yp - phi_yn)/(dy + dy), | |
+ (phi_zp - phi_zn)/(dz + dz)); | |
+ | |
// laplacian approximated by second-order central differences | |
const Float laplace_p = | |
(p_xp - (p + p) + p_xn)/(dx*dx) + | |
t@@ -1539,9 +1561,11 @@ __global__ void firstDarcySolution( | |
(p_zp - (p + p) + p_zn)/(dz*dz); | |
Float dp_expl = | |
- + (ndem*devC_dt)/(beta_f*phi*mu)*(k*laplace_p + dot(grad_k, grad_p… | |
+ + 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)); | |
+ //- dphi/(beta_f*phi*(1.0 - phi)); | |
+ - ndem*devC_dt/(beta_f*phi*(1.0 - phi)) | |
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi)); | |
// Dirichlet BC at dynamic top wall. wall0_iz will be larger than the | |
// grid if the wall isn't dynamic | |
t@@ -1607,6 +1631,7 @@ __global__ void updateDarcySolution( | |
const Float* __restrict__ dev_darcy_phi, // in | |
const Float* __restrict__ dev_darcy_dphi, // in | |
const Float* __restrict__ dev_darcy_div_v_p, // in | |
+ const Float3* __restrict__ dev_darcy_vp_avg, // in | |
const Float3* __restrict__ dev_darcy_grad_k, // in | |
const Float beta_f, // in | |
const Float mu, // in | |
t@@ -1639,11 +1664,18 @@ __global__ void updateDarcySolution( | |
// read values | |
__syncthreads(); | |
- 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 k = dev_darcy_k[cellidx]; | |
+ const Float3 grad_k = dev_darcy_grad_k[cellidx]; | |
+ const Float phi_xn = dev_darcy_phi[d_idx(x-1,y,z)]; | |
+ const Float phi = dev_darcy_phi[cellidx]; | |
+ const Float phi_xp = dev_darcy_phi[d_idx(x+1,y,z)]; | |
+ const Float phi_yn = dev_darcy_phi[d_idx(x,y-1,z)]; | |
+ const Float phi_yp = dev_darcy_phi[d_idx(x,y+1,z)]; | |
+ const Float phi_zn = dev_darcy_phi[d_idx(x,y,z-1)]; | |
+ const Float phi_zp = dev_darcy_phi[d_idx(x,y,z+1)]; | |
+ const Float dphi = dev_darcy_dphi[cellidx]; | |
//const Float div_v_p = dev_darcy_div_v_p[cellidx]; | |
+ const Float3 vp_avg = dev_darcy_vp_avg[cellidx]; | |
const Float p_old = dev_darcy_p_old[cellidx]; | |
const Float dp_expl = dev_darcy_dp_expl[cellidx]; | |
t@@ -1684,6 +1716,11 @@ __global__ void updateDarcySolution( | |
(p_yp - p_yn)/(dy + dy), | |
(p_zp - p_zn)/(dz + dz)); | |
+ const Float3 grad_phi = MAKE_FLOAT3( | |
+ (phi_xp - phi_xn)/(dx + dx), | |
+ (phi_yp - phi_yn)/(dy + dy), | |
+ (phi_zp - phi_zn)/(dz + dz)); | |
+ | |
// laplacian approximated by second-order central differences | |
const Float laplace_p = | |
(p_xp - (p + p) + p_xn)/(dx*dx) + | |
t@@ -1692,9 +1729,11 @@ __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… | |
+ + 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)); | |
+ //- dphi/(beta_f*phi*(1.0 - phi)); | |
+ - ndem*devC_dt/(beta_f*phi*(1.0 - phi)) | |
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi)); | |
// Dirichlet BC at dynamic top wall. wall0_iz will be larger than the | |
// grid if the wall isn't dynamic | |
diff --git a/src/device.cu b/src/device.cu | |
t@@ -1896,7 +1896,8 @@ __host__ void DEM::startTime() | |
darcy.c_phi, | |
dev_darcy_phi, | |
dev_darcy_dphi, | |
- dev_darcy_div_v_p); | |
+ dev_darcy_div_v_p, | |
+ dev_darcy_vp_avg); | |
cudaThreadSynchronize(); | |
if (PROFILING == 1) | |
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed, | |
t@@ -2041,6 +2042,7 @@ __host__ void DEM::startTime() | |
dev_darcy_phi, | |
dev_darcy_dphi, | |
dev_darcy_div_v_p, | |
+ dev_darcy_vp_avg, | |
dev_darcy_grad_k, | |
darcy.beta_f, | |
darcy.mu, | |
t@@ -2069,6 +2071,7 @@ __host__ void DEM::startTime() | |
dev_darcy_phi, | |
dev_darcy_dphi, | |
dev_darcy_div_v_p, | |
+ dev_darcy_vp_avg, | |
dev_darcy_grad_k, | |
darcy.beta_f, | |
darcy.mu, | |
diff --git a/src/sphere.h b/src/sphere.h | |
t@@ -305,14 +305,15 @@ class DEM { | |
Float* dev_darcy_phi; // Cell porosity | |
Float* dev_darcy_dphi; // Cell porosity change | |
Float* dev_darcy_div_v_p; // Cell particle velocity divergence | |
- Float* dev_darcy_v_p_x; // Cell particle velocity | |
- Float* dev_darcy_v_p_y; // Cell particle velocity | |
- Float* dev_darcy_v_p_z; // Cell particle velocity | |
+ //Float* dev_darcy_v_p_x; // Cell particle velocity | |
+ //Float* dev_darcy_v_p_y; // Cell particle velocity | |
+ //Float* dev_darcy_v_p_z; // Cell particle velocity | |
Float* dev_darcy_norm; // Normalized residual of epsilon values | |
Float4* dev_darcy_f_p; // Pressure gradient force on particles | |
Float* dev_darcy_k; // Cell hydraulic permeability | |
Float3* dev_darcy_grad_k; // Spatial gradient of permeability | |
Float3* dev_darcy_grad_p; // Spatial gradient of fluid pressure | |
+ Float3* dev_darcy_vp_avg; // Average particle velocity in cell | |
// Darcy functions | |
void initDarcyMem(); |