tincrease wall mass instead of applying normal stress - 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 6b29476057a9cab955d7fdaeed45eb80177fc0bf
(DIR) parent 0b6294bb9552807ee7c7671c2ca4a073b9f13076
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 15 Sep 2014 15:02:12 +0200
increase wall mass instead of applying normal stress
Diffstat:
M python/permeability-results.py | 28 ++++++++++++++++++++++++++--
M python/shear-dry.sh | 2 +-
M python/shear-starter.py | 8 +++++---
M src/integration.cuh | 8 ++++++++
4 files changed, 40 insertions(+), 6 deletions(-)
---
(DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
t@@ -21,6 +21,7 @@ dpdz = [[], [], []]
Q = [[], [], []]
phi_bar = [[], [], []]
Re = [[], [], []]
+fp_fsum = [[], [], []]
c = 0
t@@ -39,12 +40,13 @@ for c_grad_p in cvals:
Q[c] = numpy.zeros_like(K[c])
phi_bar[c] = numpy.zeros_like(K[c])
Re[c] = numpy.zeros_like(K[c])
+ fp_fsum[c] = numpy.zeros_like(K[c])
i = 0
for sid in sids:
if os.path.isfile('../output/' + sid + '.status.dat'):
- pc = PermeabilityCalc(sid, plot_evolution=False, print_results=False,
- verbose=False)
+ pc = PermeabilityCalc(sid, plot_evolution=False,
+ print_results=False, verbose=False)
K[c][i] = pc.conductivity()
pc.findPressureGradient()
pc.findCrossSectionalFlux()
t@@ -57,6 +59,20 @@ for c_grad_p in cvals:
sim = sphere.sim(sid, fluid=True)
sim.readlast(verbose=False)
Re[c][i] = numpy.mean(sim.ReynoldsNumber())
+
+ # find magnitude of fluid pressure force and total interaction force
+ '''
+ fp_magn = numpy.empty(sim.np)
+ fsum_magn = numpy.empty(sim.np)
+ for i in numpy.arange(sim.np):
+ fp_magn[i] = sim.f_p[i,:].dot(sim.f_p[i,:])
+ fsum_magn[i] = sim.f_sum[i,:].dot(sim.f_sum[i,:])
+
+ fp_fsum[c][i] = numpy.mean(fp_magn/fsum_magn)
+ # interaction forces not written in these old output files!
+ '''
+
+
else:
print(sid + ' not found')
t@@ -73,6 +89,7 @@ fig = plt.figure(figsize=(8,12))
ax1 = plt.subplot(3,1,1)
ax2 = plt.subplot(3,1,2, sharex=ax1)
ax3 = plt.subplot(3,1,3, sharex=ax1)
+#ax4 = plt.subplot(4,1,4, sharex=ax1)
lines = ['-', '--', '-.', ':']
markers = ['o', 'x', '^', '+']
for c in range(len(cvals)):
t@@ -86,6 +103,8 @@ for c in range(len(cvals)):
linestyle=lines[c], marker=markers[c], color='black')
ax3.loglog(dpdz[c], Re[c], label='$c$ = %.2f' % (cvals[c]),
linestyle=lines[c], marker=markers[c], color='black')
+ #ax4.loglog(dpdz[c], fp_fsum[c], label='$c$ = %.2f' % (cvals[c]),
+ #linestyle=lines[c], marker=markers[c], color='black')
ax1.set_ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
#ax1.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
t@@ -94,12 +113,17 @@ ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
ax3.set_ylabel('Mean Reynolds number $\\bar{Re}$ [-]')
+#ax4.set_ylabel('$\\bar{\\boldsymbol{f}_{\\Delta p}/\\bar{\\boldsymbol{f}_\\text{pf}}$ [-]')
+
ax3.set_xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
+
plt.setp(ax1.get_xticklabels(), visible=False)
+plt.setp(ax2.get_xticklabels(), visible=False)
ax1.grid()
ax2.grid()
ax3.grid()
+#ax4.grid()
#plt.subplot(3,1,2)
#plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
(DIR) diff --git a/python/shear-dry.sh b/python/shear-dry.sh
t@@ -1,5 +1,5 @@
#!/bin/sh
-#PBS -N shear-dry-hi_mu
+#PBS -N shear-dry-hi_mu-hw
#PBS -l nodes=1:ppn=3
#PBS -l walltime=19200:00:00
#PBS -q qfermi
(DIR) diff --git a/python/shear-starter.py b/python/shear-starter.py
t@@ -27,9 +27,9 @@ sim.readlast()
if fluid:
sim.sid = 'shear-sigma0=' + str(sigma0) + '-c_phi=' + str(c_phi) + \
- '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc'
+ '-c_grad_p=' + str(c_grad_p) + '-hi_mu-lo_visc-hw'
else:
- sim.sid = 'shear-sigma0=' + str(sigma0)
+ sim.sid = 'shear-sigma0=' + str(sigma0) + '-hw'
print(sim.sid)
sim.fluid = fluid
t@@ -53,7 +53,9 @@ sim.setMaxIterations(2e5)
sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
sim.c_phi[0] = c_phi
sim.c_grad_p[0] = c_grad_p
-sim.w_devs[0] = sigma0
+#sim.w_devs[0] = sigma0
+sim.w_devs[0] = 0.0
+sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
sim.mu_s[0] = 0.5
sim.mu_d[0] = 0.5
(DIR) diff --git a/src/integration.cuh b/src/integration.cuh
t@@ -344,6 +344,14 @@ __global__ void integrateWalls(
// (Wall mass is stored in x component of position Float4)
Float acc = (w_mvfd.z + N)/w_mvfd.x;
+ // Apply gravity
+ if (idx == 0)
+ acc += devC_params.g[2];
+ else if (idx == 1)
+ acc += devC_params.g[0];
+ else if (idx == 2)
+ acc += devC_params.g[1];
+
// If Wall BC is controlled by velocity, it should not change
if (wmode == 2) {
acc = 0.0;