tMerge branch 'master' of github.com:anders-dc/sphere - 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 9cf1ddce58685e261bac0c142e788d7eae415c85
(DIR) parent 7c7adcf33014a7a965f18f1b6285c7c85c6207f9
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 23 Feb 2017 15:05:00 -0800
Merge branch 'master' of github.com:anders-dc/sphere
Diffstat:
M README.rst | 4 ++++
M python/creep-master.py | 24 ++++++++++++++----------
M python/halfshear-darcy-combined.py | 25 +++++++++++++++----------
M python/halfshear-darcy-rate-depend… | 203 ++++++++++++++++++++++++++++---
M python/shortening.py | 12 ++++++------
M python/sphere.py | 55 ++++++++++++++++++++++++++++---
M src/darcy.cpp | 2 ++
M src/darcy.cuh | 411 ++++++++++++++++++++++---------
M src/datatypes.h | 2 ++
M src/device.cu | 29 ++++++++++++++---------------
M src/file_io.cpp | 26 +++++++++++++++++++++++++-
M src/navierstokes.cpp | 2 ++
M src/sphere.h | 1 +
M src/version.h | 2 +-
14 files changed, 617 insertions(+), 181 deletions(-)
---
(DIR) diff --git a/README.rst b/README.rst
t@@ -103,6 +103,10 @@ Publications
``sphere`` has been used to produce results in the following scientific
publications and presentations:
+- Damsgaard, A., D.L. Egholm, L.H. Beem, S. Tulaczyk, N.K. Larsen, J.A.
+ Piotrowski, and M.R. Siegfried (2016), Ice flow dynamics forced by water
+ pressure variations in subglacial granular beds, Geophysical Research Letters,
+ 43, `doi:10.1002/2016GL071579 <http://dx.doi.org/10.1002/2016GL071579>`_.
- Damsgaard, A., D.L. Egholm, J.A. Piotrowski, S. Tulaczyk, N.K. Larsen, and
C.F. Brædstrup (2015), A new methodology to simulate subglacial deformation of
water-saturated granular material, The Cryosphere, 9, 2183-2200,
(DIR) diff --git a/python/creep-master.py b/python/creep-master.py
t@@ -7,13 +7,13 @@ import numpy
### EXPERIMENT SETUP ###
initialization = False
consolidation = False
-shearing = True
+shearing = False
creeping = True
rendering = False
plots = True
# Common simulation id
-sim_id = "creep1"
+sim_id = "creep2"
# Fluid-pressure gradient [Pa/m]
dpdx = -100.0
t@@ -37,7 +37,8 @@ np = 1e4
### INITIALIZATION ###
# New class
-init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
+#init = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
+init = sphere.sim(np = np, nd = 3, nw = 0, sid = 'creep1' + "-init")
# Uniform radii from 0.8 cm to 1.2 cm
init.generateRadii(psd = 'uni', mean = 0.005, variance = 0.001)
t@@ -47,7 +48,7 @@ init.defaultParams(gamma_n = 100.0, mu_s = 0.6, mu_d = 0.6)
init.setYoungsModulus(1e8)
# Add gravity
-init.g[2] = -9.81
+init.g[2] = -g
# Periodic x and y boundaries
init.periodicBoundariesXY()
t@@ -84,7 +85,9 @@ cons = sphere.sim(np = init.np, nw = 1, sid = sim_id +
# Read last output file of initialization step
lastf = init.status()
-cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
+#cons.readbin("../output/" + sim_id + "-init.output{:0=5}.bin".format(lastf), verbose=False)
+cons.readbin("../output/" + 'creep1' + "-init.output{:0=5}.bin".format(lastf),
+ verbose=False)
# Periodic x and y boundaries
cons.periodicBoundariesXY()
t@@ -97,7 +100,7 @@ cons.gamma_wn[0] = 0.0
cons.gamma_wt[0] = 0.0
cons.rho[0] = rho_g
-cons.g[2] = g
+cons.g[2] = -g
# Set duration of simulation
cons.initTemporal(total = 1.5)
t@@ -136,7 +139,7 @@ shear.readbin("../output/" + sim_id +
shear.periodicBoundariesXY()
shear.rho[0] = rho_g
-shear.g[2] = g
+shear.g[2] = -g
# Disable particle viscosities
shear.gamma_n[0] = 0.0
t@@ -191,12 +194,13 @@ creep.setFluidCompressibility(4.6e-10) # water at 25 C
creep.setFluidTopNoFlow()
creep.setFluidBottomNoFlow()
creep.setFluidXFixedPressure()
+creep.adaptiveGrid()
# set fluid pressures at the boundaries and internally
-dx = sim.L[0]/sim.num[0]
+dx = creep.L[0]/creep.num[0]
for ix in range(creep.num[0]):
x = ix + 0.5*dx
- creep.p_f[ix,:,:] = numpy.abs(sim.L[0]*dpdx) - x*dpdx
+ creep.p_f[ix,:,:] = numpy.abs(creep.L[0]*dpdx) + x*dpdx
creep.zeroKinematics()
t@@ -211,7 +215,7 @@ creep.color[numpy.nonzero(creep.fixvel == 1)] == -1
creep.adaptiveGrid()
# Set duration of simulation
-creep.initTemporal(total = 60.0)
+creep.initTemporal(total = 20.0)
if (creeping == True):
(DIR) diff --git a/python/halfshear-darcy-combined.py b/python/halfshear-darcy-combined.py
t@@ -15,17 +15,18 @@ matplotlib.rcParams.update({'font.size': 7, 'font.family': 'sans-serif'})
matplotlib.rc('text', usetex=True)
matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
-if len(sys.argv) > 1:
- sid = sys.argv[1]
-else:
- raise Exception('Missing sid')
+#if len(sys.argv) > 1:
+ #sid = sys.argv[1]
+#else:
+ #raise Exception('Missing sid')
+sid = 'halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2'
outformat = 'pdf'
fluid = True
#threshold = 100.0 # [N]
if len(sys.argv) > 3:
calculateforcechains = bool(int(sys.argv[3]))
else:
- calculateforcechains = True
+ calculateforcechains = False
calculateforcechainhistory = False
legend_alpha=0.7
linewidth=0.5
t@@ -35,10 +36,10 @@ rasterized=True # rasterize colored areas (pcolormesh and colorbar)
izfactor = 1 # factor for vertical discretization in poros
-if len(sys.argv) > 2:
- t_DEM_to_t_real = float(sys.argv[2])
-else:
- raise Exception('Missing t_DEM_to_t_real')
+#if len(sys.argv) > 2:
+# t_DEM_to_t_real = float(sys.argv[2])
+#else:
+# raise Exception('Missing t_DEM_to_t_real')
###################
t@@ -149,6 +150,7 @@ else:
n, nkept, coordinationnumber = numpy.loadtxt(sid + '-fc.txt')
# Transform time from model time to real time [s]
+t_DEM_to_t_real = sim.mu/1.787e-3
t = t/t_DEM_to_t_real
## integrate velocities to displacement along x (xdispint)
t@@ -225,11 +227,14 @@ ax1.text(bbox_x, bbox_y, 'a',
## ax3: v, ax4: unused
ax3 = plt.subplot(5, 1, 2, sharex=ax1)
+#ax3.semilogy(t, v*t_DEM_to_t_real, 'k', linewidth=linewidth)
ax3.semilogy(t, v, 'k', linewidth=linewidth)
ax3.set_ylabel('Shear velocity [ms$^{-1}$]')
+#ax3.set_ylabel('Shear velocity [ma$^{-1}$]')
+#ax3.set_ylabel('Shear velocity [ma$^{-1}$]')
# shade stick periods
collection = matplotlib.collections.BrokenBarHCollection.span_where(
- t, ymin=1.0e-11, ymax=1.0,
+ t, ymin=1.0e-11, ymax=1.0e4,
where=numpy.isclose(v, 0.0),
facecolor='black', alpha=0.2,
linewidth=0)
(DIR) diff --git a/python/halfshear-darcy-rate-dependence.py b/python/halfshear-darcy-rate-dependence.py
t@@ -33,6 +33,8 @@ sids =\
'halfshear-darcy-sigma0=10000.0-k_c=2e-16-mu=2.08e-07-ss=3000.0-A=4500.0-f=0.2',
]
#['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.797e-06-ss=10000.0-A=70000.0-f=0.2']
+sids = \
+ ['halfshear-darcy-sigma0=80000.0-k_c=3.5e-13-mu=1.04e-07-ss=10000.0-A=70000.0-f=0.2']
outformat = 'pdf'
fluid = True
#threshold = 100.0 # [N]
t@@ -66,12 +68,13 @@ for sid in sids:
dilation = numpy.empty(sim.status())
for i in numpy.arange(sim.status()):
sim.readstep(i+1, verbose=False)
+ t_DEM_to_t_real = sim.mu[0]/1.787e-3
#tau = sim.shearStress()
tau[i] = sim.w_tau_x # defined shear stress
#tau[i] = sim.shearStress()[0] # measured shear stress along x
N[i] = sim.currentNormalStress() # defined normal stress
#v[i] = sim.shearVel()
- shearstrainrate[i] = sim.shearStrainRate()
+ shearstrainrate[i] = sim.shearStrainRate()*t_DEM_to_t_real
shearstrain[i] = sim.shearStrain()
t[i] = sim.currentTime()
t@@ -88,27 +91,168 @@ for sid in sids:
tau_nonzero = tau[idx]
N_nonzero = N[idx]
shearstrain_nonzero = shearstrain[idx]
+
t_nonzero = t[idx]
### "Critical" state fit
# The algorithm uses the Levenberg-Marquardt algorithm through leastsq
idxfit = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
- (shearstrainrate_nonzero < 0.1) &
+ (shearstrainrate_nonzero < 0.1*t_DEM_to_t_real) &
(((t_nonzero > 6.0) & (t_nonzero < 8.0)) |
((t_nonzero > 11.0) & (t_nonzero < 13.0)) |
- ((t_nonzero > 16.0) & (t_nonzero < 18.0))))
+ ((t_nonzero > 16.0) & (t_nonzero < 18.0)) |
+ ((t_nonzero > 21.0) & (t_nonzero < 23.0)) |
+ ((t_nonzero > 26.0) & (t_nonzero < 28.0)) |
+ ((t_nonzero > 31.0) & (t_nonzero < 33.0)) |
+ ((t_nonzero > 36.0) & (t_nonzero < 38.0)) |
+ ((t_nonzero > 41.0) & (t_nonzero < 43.0)) |
+ ((t_nonzero > 46.0) & (t_nonzero < 48.0)) |
+ ((t_nonzero > 51.0) & (t_nonzero < 53.0)) |
+ ((t_nonzero > 56.0) & (t_nonzero < 58.0))))
+
#popt, pvoc = scipy.optimize.curve_fit(
#creep_rheology, tau/N, shearstrainrate)
popt, pvoc = scipy.optimize.curve_fit(
creep_rheology1, tau_nonzero[idxfit]/N_nonzero[idxfit],
shearstrainrate_nonzero[idxfit])
- print '# Critical state'
+ popt2, pvoc2 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit]/N_nonzero[idxfit],
+ shearstrainrate_nonzero[idxfit])
+ print '# Creep rheology 1'
print popt
print pvoc
n = popt[0] # stress exponent
#A = popt[1] # stress exponent
A = 1.
+ print '# Creep rheology 2'
+ print popt2
+ print pvoc2
+ n = popt2[0] # stress exponent
+ A = popt2[1] # stress exponent
+ #A = 1.
+
+
+ '''
+ # Investigate n and A evolution from cycle to cycle
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 1.0) & (t_nonzero < 3.0))))
+ popt_0, pvoc_0 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 6.0) & (t_nonzero < 8.0))))
+ popt_1, pvoc_1 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 11.0) & (t_nonzero < 13.0))))
+ popt_2, pvoc_2 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 16.0) & (t_nonzero < 18.0))))
+ popt_3, pvoc_3 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 21.0) & (t_nonzero < 23.0))))
+ popt_4, pvoc_4 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 26.0) & (t_nonzero < 28.0))))
+ popt_5, pvoc_5 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 31.0) & (t_nonzero < 33.0))))
+ popt_6, pvoc_6 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 36.0) & (t_nonzero < 38.0))))
+ popt_7, pvoc_7 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 41.0) & (t_nonzero < 43.0))))
+ popt_8, pvoc_8 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 46.0) & (t_nonzero < 48.0))))
+ popt_9, pvoc_9 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 51.0) & (t_nonzero < 53.0))))
+ popt_10, pvoc_10 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ idxfit_ = numpy.nonzero((tau_nonzero/N_nonzero < 0.38) &
+ (shearstrainrate_nonzero < 0.1) &
+ (((t_nonzero > 56.0) & (t_nonzero < 58.0))))
+ popt_11, pvoc_11 = scipy.optimize.curve_fit(
+ creep_rheology2, tau_nonzero[idxfit_]/N_nonzero[idxfit_],
+ shearstrainrate_nonzero[idxfit_])
+
+ n_list = [
+ popt_0[0],
+ popt_1[0],
+ popt_2[0],
+ popt_3[0],
+ popt_4[0],
+ popt_5[0],
+ popt_6[0],
+ popt_7[0],
+ popt_8[0],
+ popt_9[0],
+ popt_10[0],
+ popt_11[0]]
+
+ A_list = [
+ popt_0[1],
+ popt_1[1],
+ popt_2[1],
+ popt_3[1],
+ popt_4[1],
+ popt_5[1],
+ popt_6[1],
+ popt_7[1],
+ popt_8[1],
+ popt_9[1],
+ popt_10[1],
+ popt_11[1]]
+ cycle = numpy.arange(1, len(n_list)+1)
+
+ print(n_list)
+ print(A_list)
+ '''
+
+
friction = tau_nonzero/N_nonzero
x_min = numpy.floor(numpy.min(friction))
x_max = numpy.ceil(numpy.max(friction)) + 0.05
t@@ -117,14 +261,14 @@ for sid in sids:
numpy.max(tau_nonzero[idxfit]/N_nonzero[idxfit]),
100)
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 < 2.5)))
idxfit2 = numpy.nonzero((shearstrain_nonzero < 0.1) &
(tau_nonzero/N_nonzero < 0.38))
- '''
+
popt2, pvoc2 = scipy.optimize.curve_fit(
creep_rheology2, tau_nonzero[idxfit2]/N_nonzero[idxfit2],
shearstrainrate_nonzero[idxfit2])
t@@ -159,11 +303,11 @@ for sid in sids:
#cmap=matplotlib.cm.get_cmap('afmhot'))
## plastic limit
- x = [0.28, 0.28]
+ x = [0.26, 0.26]
y = ax1.get_ylim()
limitcolor = '#333333'
ax1.plot(x, y, '--', linewidth=2, color=limitcolor)
- ax1.text(x[0]+0.03, 2.0e-4,
+ ax1.text(x[0]+0.03, 2.0e-4*t_DEM_to_t_real,
'Yield strength', rotation=90, color=limitcolor,
bbox={'fc':'#ffffff', 'pad':3, 'alpha':0.7})
t@@ -171,9 +315,10 @@ for sid in sids:
## Fit
ax1.plot(friction_fit, strainrate_fit)
#ax1.plot(friction_fit2, strainrate_fit2)
- ax1.annotate('$\\dot{{\\gamma}} = (\\tau/N)^{{{:.1f}}}$'.format(n),
+ ax1.annotate('$\\dot{{\\gamma}}$ = ' + '{:.2e} s$^{{-1}}$'.format(A) +
+ '$(\\tau/N)^{{{:.2f}}}$'.format(n),
xy = (friction_fit[40], strainrate_fit[40]),
- xytext = (0.32+0.05, 2.0e-9),
+ xytext = (0.31+0.05, 2.0e-9*t_DEM_to_t_real),
arrowprops=dict(facecolor='blue', edgecolor='blue', shrink=0.1,
width=1, headwidth=4, frac=0.2),)
#xytext = (friction_fit[50]+0.15, strainrate_fit[50]-1.0e-5))#,
t@@ -200,9 +345,9 @@ for sid in sids:
plt.tight_layout()
filename = sid + '-rate-dependence.' + outformat
plt.savefig(filename)
- plt.close()
- shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
print(filename)
+ plt.close()
+ #shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
## dilation vs. rate
'''
t@@ -223,11 +368,9 @@ for sid in sids:
ax1.text(x[0]+0.03, 2.0e-4,
'Yield strength', rotation=90, color=limitcolor,
bbox={'fc':'#ffffff', 'pad':3, 'alpha':0.7})
-
+ '''
## Fit
- '''
- '''
ax1.plot(friction_fit, strainrate_fit)
#ax1.plot(friction_fit2, strainrate_fit2)
ax1.annotate('$\\dot{\\gamma} = (\\tau/N)^{6.4}$',
t@@ -260,4 +403,32 @@ for sid in sids:
plt.close()
shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
print(filename)
- '''
+
+
+
+ fig = plt.figure(figsize=(3.5,2.5))
+ ax1 = plt.subplot(111)
+ #ax1.semilogy(N/1000., v)
+ #ax1.semilogy(tau_nonzero/N_nonzero, v_nonzero, '+k')
+ #ax1.plot(tau/N, v, '.')
+ #CS = ax1.scatter(friction, v_nonzero, c=shearstrain_nonzero,
+ #linewidth=0)
+ CS = ax1.plot(cycle, A_list, '-k')
+
+ ax1.set_xlabel('Forcing cycle [-]')
+ #ax1.set_ylabel('Shear velocity [m/s]')
+ ax1.set_ylabel('$A$ [s$^{-1}$]')
+
+ ax2 = ax1.twinx()
+ ax2color = 'blue'
+ lns = ax2.plot(cycle, n_list, ':', color=ax2color)
+ ax2.set_ylabel('$n$ [-]')
+ for tl in ax2.get_yticklabels():
+ tl.set_color(ax2color)
+
+ plt.tight_layout()
+ filename = sid + '-rate-dependence-A-n.' + outformat
+ plt.savefig(filename)
+ print(filename)
+ plt.close()
+ #shutil.copyfile(filename, '/home/adc/articles/own/3/graphics/' + filename)
(DIR) diff --git a/python/shortening.py b/python/shortening.py
t@@ -60,15 +60,15 @@ for z in range(nz):
# move to x=0
min_x = numpy.min(sim.x[:,0] - sim.radius[:])
-sim.x[:,0] = sim.x[:,0] - min_x
+sim.x[:,0] = sim.x[:,0] - min_x
# move to y=0
min_y = numpy.min(sim.x[:,1] - sim.radius[:])
-sim.x[:,1] = sim.x[:,1] - min_y
+sim.x[:,1] = sim.x[:,1] - min_y
# move to z=0
min_z = numpy.min(sim.x[:,2] - sim.radius[:])
-sim.x[:,2] = sim.x[:,2] - min_z
+sim.x[:,2] = sim.x[:,2] - min_z
#sim.defineWorldBoundaries(L=[Lx, Lz*3, Ly])
sim.defineWorldBoundaries(L=[numpy.max(sim.x[:,0] + sim.radius[:]), Lz*3, Ly])
t@@ -98,9 +98,9 @@ sim.uniaxialStrainRate(wvel = 0.0)
sim.initTemporal(total=3.0, file_dt = 0.01)
sim.zeroKinematics()
-#sim.run(dry=True)
-#sim.run()
-#sim.writeVTKall()
+sim.run(dry=True)
+sim.run()
+sim.writeVTKall()
## Shortening
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -24,7 +24,7 @@ numpy.seterr(all='warn', over='raise')
# Sphere version number. This field should correspond to the value in
# `../src/version.h`.
-VERSION = 2.14
+VERSION = 2.15
# Transparency on plot legends
legend_alpha = 0.5
t@@ -337,6 +337,8 @@ class sim:
self.p_mod_f = numpy.zeros(1, dtype=numpy.float64) # Frequency [Hz]
self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
+ ## Fluid solver parameters
+
if self.cfd_solver[0] == 1: # Darcy solver
# Boundary conditions at the sides of the fluid grid
# 0: Dirichlet
t@@ -363,7 +365,10 @@ class sim:
self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
- ## Solver parameters
+ # Hold pressures constant in fluid cell (0: True, 1: False)
+ self.p_f_constant = numpy.zeros(
+ (self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.int32)
# Navier-Stokes
if self.cfd_solver[0] == 0:
t@@ -440,7 +445,6 @@ class sim:
'''
Called when to sim objects are compared. Returns 0 if the values
are identical.
- TODO: Replace print(#) with print("field name")
'''
if self.version != other.version:
print('version')
t@@ -690,6 +694,9 @@ class sim:
elif self.bc_top_flux != other.bc_top_flux:
print('bc_top_flux')
return 91
+ elif (self.p_f_constant != other.p_f_constant).any():
+ print('p_f_constant')
+ return 96
if self.cfd_solver == 0:
if self.gamma != other.gamma:
t@@ -1190,7 +1197,7 @@ class sim:
self.p_mod_phi =\
numpy.fromfile(fh, dtype=numpy.float64, count=1)
- if self.version >= 2.12 and self.cfd_solver == 1:
+ if self.version >= 2.12 and self.cfd_solver[0] == 1:
self.bc_xn =\
numpy.fromfile(fh, dtype=numpy.int32, count=1)
self.bc_xp =\
t@@ -1217,6 +1224,23 @@ class sim:
self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
+ if self.version >= 2.15:
+ self.p_f_constant = \
+ numpy.empty((self.num[0],self.num[1],self.num[2]),
+ dtype=numpy.int32)
+
+ for z in numpy.arange(self.num[2]):
+ for y in numpy.arange(self.num[1]):
+ for x in numpy.arange(self.num[0]):
+ self.p_f_constant[x,y,z] = \
+ numpy.fromfile(fh,
+ dtype=numpy.int32,
+ count=1)
+ else:
+ self.p_f_constant = numpy.zeros(
+ (self.num[0], self.num[1], self.num[2]),
+ dtype=numpy.int32)
+
if self.version >= 2.0 and self.cfd_solver == 0:
self.gamma = \
numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1465,6 +1489,12 @@ class sim:
fh.write(self.bc_bot_flux.astype(numpy.float64))
fh.write(self.bc_top_flux.astype(numpy.float64))
+ for z in numpy.arange(self.num[2]):
+ for y in numpy.arange(self.num[1]):
+ for x in numpy.arange(self.num[0]):
+ fh.write(self.p_f_constant[x,y,z].astype(
+ numpy.int32))
+
# Navier Stokes
if self.cfd_solver[0] == 0:
fh.write(self.gamma.astype(numpy.float64))
t@@ -2087,6 +2117,14 @@ class sim:
else:
K.SetNumberOfTuples(grid.GetNumberOfPoints())
+ p_f_constant = vtk.vtkDoubleArray()
+ p_f_constant.SetName("Constant pressure [-]")
+ p_f_constant.SetNumberOfComponents(1)
+ if cell_centered:
+ p_f_constant.SetNumberOfTuples(grid.GetNumberOfCells())
+ else:
+ p_f_constant.SetNumberOfTuples(grid.GetNumberOfPoints())
+
# insert values
for z in range(self.num[2]):
for y in range(self.num[1]):
t@@ -2100,6 +2138,7 @@ class sim:
if self.cfd_solver[0] == 1:
k.SetValue(idx, self.k[x,y,z])
K.SetValue(idx, self.K[x,y,z])
+ p_f_constant.SetValue(idx, self.p_f_constant[x,y,z])
# add pres array to grid
if cell_centered:
t@@ -2111,6 +2150,7 @@ class sim:
if self.cfd_solver[0] == 1:
grid.GetCellData().AddArray(k)
grid.GetCellData().AddArray(K)
+ grid.GetCellData().AddArray(p_f_constant)
else:
grid.GetPointData().AddArray(pres)
grid.GetPointData().AddArray(vel)
t@@ -2120,6 +2160,7 @@ class sim:
if self.cfd_solver[0] == 1:
grid.GetPointData().AddArray(k)
grid.GetPointData().AddArray(K)
+ grid.GetPointData().AddArray(p_f_constant)
# write VTK XML image data file
writer = vtk.vtkXMLImageDataWriter()
t@@ -3403,6 +3444,10 @@ class sim:
self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
+ self.p_f_constant = numpy.zeros((self.num[0], self.num[1],
+ self.num[2]),
+ dtype=numpy.int32)
+
# Fluid solver type
# 0: Navier Stokes (fluid with inertia)
# 1: Stokes-Darcy (fluid without inertia)
t@@ -4459,7 +4504,7 @@ class sim:
if dry:
dryarg = "--dry "
if valgrind:
- valgrindbin = "valgrind -q "
+ valgrindbin = "valgrind -q --track-origins=yes "
if cudamemcheck:
cudamemchk = "cuda-memcheck --leak-check full "
if self.fluid:
(DIR) diff --git a/src/darcy.cpp b/src/darcy.cpp
t@@ -30,6 +30,7 @@ void DEM::initDarcyMem()
darcy.norm = new Float[ncells]; // normalized residual of epsilon
darcy.f_p = new Float4[np]; // pressure force on particles
darcy.k = new Float[ncells]; // hydraulic pressure
+ darcy.p_constant = new int[ncells]; // constant pressure (0: no, 1: yes)
}
unsigned int DEM::darcyCells()
t@@ -59,6 +60,7 @@ void DEM::freeDarcyMem()
delete[] darcy.norm;
delete[] darcy.f_p;
delete[] darcy.k;
+ delete[] darcy.p_constant;
}
// 3D index to 1D index
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -41,6 +41,8 @@ void DEM::initDarcyMemDev(void)
//cudaMalloc((void**)&dev_darcy_v_p_x, memSizeFace); // v_p.x
//cudaMalloc((void**)&dev_darcy_v_p_y, memSizeFace); // v_p.y
//cudaMalloc((void**)&dev_darcy_v_p_z, memSizeFace); // v_p.z
+ cudaMalloc((void**)&dev_darcy_p_constant,
+ sizeof(int)*darcyCells()); // grad(pressure)
checkForCudaErrors("End of initDarcyMemDev");
}
t@@ -67,6 +69,7 @@ void DEM::freeDarcyMemDev()
//cudaFree(dev_darcy_v_p_y);
//cudaFree(dev_darcy_v_p_z);
cudaFree(dev_darcy_grad_p);
+ cudaFree(dev_darcy_p_constant);
}
// Transfer to device
t@@ -91,6 +94,8 @@ void DEM::transferDarcyToGlobalDeviceMemory(int statusmsg)
cudaMemcpy(dev_darcy_dphi, darcy.dphi, memSizeF, cudaMemcpyHostToDevice);
cudaMemcpy(dev_darcy_f_p, darcy.f_p, sizeof(Float4)*np,
cudaMemcpyHostToDevice);
+ cudaMemcpy(dev_darcy_p_constant, darcy.p_constant, sizeof(int)*darcyCells(),
+ cudaMemcpyHostToDevice);
checkForCudaErrors("End of transferDarcyToGlobalDeviceMemory");
//if (verbose == 1 && statusmsg == 1)
t@@ -114,6 +119,8 @@ void DEM::transferDarcyFromGlobalDeviceMemory(int statusmsg)
cudaMemcpy(darcy.f_p, dev_darcy_f_p, sizeof(Float4)*np,
cudaMemcpyDeviceToHost);
cudaMemcpy(darcy.k, dev_darcy_k, memSizeF, cudaMemcpyDeviceToHost);
+ cudaMemcpy(darcy.p_constant, dev_darcy_p_constant, sizeof(int)*darcyCells(),
+ cudaMemcpyDeviceToHost);
checkForCudaErrors("End of transferDarcyFromGlobalDeviceMemory", 0);
if (verbose == 1 && statusmsg == 1)
t@@ -612,7 +619,7 @@ __global__ void findDarcyPorositiesLinear(
Float cell_volume = dx*dy*dz;
if (z == nz - 1)
- cell_volume *= 0.75;
+ cell_volume *= 0.875;
// 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)));
t@@ -710,7 +717,7 @@ __global__ void findDarcyPorositiesLinear(
//}
//dev_darcy_phi[cellidx] = phi;
//dev_darcy_dphi[cellidx] = dphi;
- dev_darcy_phi[cellidx] = 0.999;
+ dev_darcy_phi[cellidx] = 0.9;
dev_darcy_dphi[cellidx] = 0.0;
//dev_darcy_vp_avg[cellidx] = MAKE_FLOAT3(0.0, 0.0, 0.0);
t@@ -798,40 +805,6 @@ __global__ void copyDarcyPorositiesToBottom(
}
-// Copy the porosity, porosity change, div_v_p and vp_avg values to the grid top
-// from the grid interior if the grid is adaptive (grid.adaptive == 1).
-__global__ void copyDarcyPorositiesToTop(
- Float* __restrict__ dev_darcy_phi, // in + out
- Float* __restrict__ dev_darcy_dphi, // in + out
- Float* __restrict__ dev_darcy_div_v_p, // in + out
- Float3* __restrict__ dev_darcy_vp_avg) // in + out
-{
- // 3D thread index
- const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
- const unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
- const unsigned int z = blockDim.z * blockIdx.z + threadIdx.z;
-
- // Grid dimensions
- const unsigned int nx = devC_grid.num[0];
- const unsigned int ny = devC_grid.num[1];
- const unsigned int nz = devC_grid.num[2];
-
- // check that we are not outside the fluid grid
- if (devC_grid.adaptive == 1 &&
- x < nx && y < ny && z == nz - 1) {
-
- const unsigned int readidx = d_idx(x, y, nz - 2);
- const unsigned int writeidx = d_idx(x, y, z);
-
- __syncthreads();
- dev_darcy_phi[writeidx] = dev_darcy_phi[readidx];
- dev_darcy_dphi[writeidx] = dev_darcy_dphi[readidx];
- dev_darcy_div_v_p[writeidx] = dev_darcy_div_v_p[readidx];
- dev_darcy_vp_avg[writeidx] = dev_darcy_vp_avg[readidx];
- }
-}
-
-
// Find the porosity in each cell on the base of a sphere, centered at the cell
// center.
__global__ void findDarcyPorosities(
t@@ -1686,6 +1659,7 @@ __global__ void firstDarcySolution(
const int bc_top, // in
const unsigned int ndem, // in
const unsigned int wall0_iz, // in
+ const int* __restrict__ dev_darcy_p_constant, // in
Float* __restrict__ dev_darcy_dp_expl) // out
{
// 3D thread index
t@@ -1710,8 +1684,8 @@ __global__ void firstDarcySolution(
// read values
__syncthreads();
- const Float k = dev_darcy_k[cellidx];
- const Float3 grad_k = dev_darcy_grad_k[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)];
t@@ -1722,14 +1696,23 @@ __global__ void firstDarcySolution(
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 int p_constant = dev_darcy_p_constant[cellidx];
- Float p_xn = dev_darcy_p[d_idx(x-1,y,z)];
- const Float p = dev_darcy_p[cellidx];
- Float p_xp = dev_darcy_p[d_idx(x+1,y,z)];
- Float p_yn = dev_darcy_p[d_idx(x,y-1,z)];
- Float p_yp = dev_darcy_p[d_idx(x,y+1,z)];
- Float p_zn = dev_darcy_p[d_idx(x,y,z-1)];
- Float p_zp = dev_darcy_p[d_idx(x,y,z+1)];
+ Float p_xn = dev_darcy_p[d_idx(x-1,y,z)];
+ const Float p = dev_darcy_p[cellidx];
+ Float p_xp = dev_darcy_p[d_idx(x+1,y,z)];
+ Float p_yn = dev_darcy_p[d_idx(x,y-1,z)];
+ Float p_yp = dev_darcy_p[d_idx(x,y+1,z)];
+ Float p_zn = dev_darcy_p[d_idx(x,y,z-1)];
+ Float p_zp = dev_darcy_p[d_idx(x,y,z+1)];
+
+ const Float k_xn = dev_darcy_k[d_idx(x-1,y,z)];
+ const Float k = dev_darcy_k[cellidx];
+ const Float k_xp = dev_darcy_k[d_idx(x+1,y,z)];
+ const Float k_yn = dev_darcy_k[d_idx(x,y-1,z)];
+ const Float k_yp = dev_darcy_k[d_idx(x,y+1,z)];
+ const Float k_zn = dev_darcy_k[d_idx(x,y,z-1)];
+ const Float k_zp = dev_darcy_k[d_idx(x,y,z+1)];
// Neumann BCs
if (x == 0 && bc_xn == 1)
t@@ -1745,44 +1728,121 @@ __global__ void firstDarcySolution(
if (z == nz-1 && bc_top == 1)
p_zp = p;
- // upwind coefficients for grad(p) determined from values of k
- // k = 1.0: backwards difference
- // k = -1.0: forwards difference
- /*const Float3 e_k = MAKE_FLOAT3(
- copysign(1.0, grad_k.x),
- copysign(1.0, grad_k.y),
- copysign(1.0, grad_k.z));
-
- // gradient approximated by first-order forward differences
- const Float3 grad_p = MAKE_FLOAT3(
- ((1.0 + e_k.x)*(p - p_xn) + (1.0 - e_k.x)*(p_xp - p))/(dx + dx),
- ((1.0 + e_k.y)*(p - p_yn) + (1.0 - e_k.y)*(p_yp - p))/(dy + dy),
- ((1.0 + e_k.z)*(p - p_zn) + (1.0 - e_k.z)*(p_zp - p))/(dz + dz)
- );*/
-
- // gradient approximated by first-order central differences
- const Float3 grad_p = MAKE_FLOAT3(
+ /*
+ // gradients approximated by first-order central differences, order of
+ // approximation is O(dx*dx)
+ const Float3 grad_p_central = MAKE_FLOAT3(
(p_xp - p_xn)/(dx + dx),
(p_yp - p_yn)/(dy + dy),
(p_zp - p_zn)/(dz + dz));
- const Float3 grad_phi = MAKE_FLOAT3(
+ const Float3 grad_k_central = MAKE_FLOAT3(
+ (k_xp - k_xn)/(dx + dx),
+ (k_yp - k_yn)/(dy + dy),
+ (k_zp - k_zn)/(dz + dz));
+ */
+
+ const Float3 grad_phi_central = MAKE_FLOAT3(
(phi_xp - phi_xn)/(dx + dx),
(phi_yp - phi_yn)/(dy + dy),
(phi_zp - phi_zn)/(dz + dz));
+ /*
+ // upwind coefficients for grad(p) determined from values of k
+ // e_p = 1.0: backwards difference
+ // e_p = -1.0: forwards difference
+ const Float3 e_p = MAKE_FLOAT3(
+ copysign(1.0, -(p_xp - p_xn)),
+ copysign(1.0, -(p_yp - p_yn)),
+ copysign(1.0, -(p_zp - p_zn)));
+
+ // gradient approximated by first-order forward differences, order of
+ // approximation is O(dx)
+ const Float3 grad_p_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(p - p_xn) + (1.0 - e_p.x)*(p_xp - p))/dx,
+ ((1.0 + e_p.y)*(p - p_yn) + (1.0 - e_p.y)*(p_yp - p))/dy,
+ ((1.0 + e_p.z)*(p - p_zn) + (1.0 - e_p.z)*(p_zp - p))/dz
+ );
+
+ const Float3 grad_k_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(k - k_xn) + (1.0 - e_p.x)*(k_xp - k))/dx,
+ ((1.0 + e_p.y)*(k - k_yn) + (1.0 - e_p.y)*(k_yp - k))/dy,
+ ((1.0 + e_p.z)*(k - k_zn) + (1.0 - e_p.z)*(k_zp - k))/dz
+ );
+
+ const Float3 grad_phi_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(phi - phi_xn)
+ + (1.0 - e_p.x)*(phi_xp - phi))/dx,
+ ((1.0 + e_p.y)*(phi - phi_yn)
+ + (1.0 - e_p.y)*(phi_yp - phi))/dy,
+ ((1.0 + e_p.z)*(phi - phi_zn)
+ + (1.0 - e_p.z)*(phi_zp - phi))/dz
+ );
+
+ // Average central and upwind discretizations to get intermediate order
+ // of approximation
+ Float gamma = 0.5; // in [0;1], where 0: fully central, 1: fully upwind
+
+ const Float3 grad_p = MAKE_FLOAT3(
+ gamma*grad_p_upwind.x,
+ gamma*grad_p_upwind.y,
+ gamma*grad_p_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_p_central.x,
+ (1.0 - gamma) * grad_p_central.y,
+ (1.0 - gamma) * grad_p_central.z);
+
+ const Float3 grad_k = MAKE_FLOAT3(
+ gamma*grad_k_upwind.x,
+ gamma*grad_k_upwind.y,
+ gamma*grad_k_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_k_central.x,
+ (1.0 - gamma) * grad_k_central.y,
+ (1.0 - gamma) * grad_k_central.z);
+
+ const Float3 grad_phi = MAKE_FLOAT3(
+ gamma*grad_phi_upwind.x,
+ gamma*grad_phi_upwind.y,
+ gamma*grad_phi_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_phi_central.x,
+ (1.0 - gamma) * grad_phi_central.y,
+ (1.0 - gamma) * grad_phi_central.z);
+
// laplacian approximated by second-order central differences
const Float laplace_p =
(p_xp - (p + p) + p_xn)/(dx*dx) +
(p_yp - (p + p) + p_yn)/(dy*dy) +
(p_zp - (p + p) + p_zn)/(dz*dz);
+ */
+
+ // Solve div(k*grad(p)) as a single term, using harmonic mean for
+ // permeability. k_HM*grad(p) is found between the pressure nodes.
+ const Float div_k_grad_p =
+ (2.*k_xp*k/(k_xp + k) *
+ (p_xp - p)/dx
+ -
+ 2.*k_xn*k/(k_xn + k) *
+ (p - p_xn)/dx)/dx
+ +
+ (2.*k_yp*k/(k_yp + k) *
+ (p_yp - p)/dy
+ -
+ 2.*k_yn*k/(k_yn + k) *
+ (p - p_yn)/dy)/dy
+ +
+ (2.*k_zp*k/(k_zp + k) *
+ (p_zp - p)/dz
+ -
+ 2.*k_zn*k/(k_zn + k) *
+ (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))
+ ndem*devC_dt/(beta_f*phi*mu)*div_k_grad_p
//- div_v_p/(beta_f*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));
+ -(ndem*devC_dt/(beta_f*phi*(1.0 - phi)))
+ //*(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi));
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi_central));
// Dirichlet BC at fixed-pressure boundaries and at the dynamic top
// wall. wall0_iz will be larger than the grid if the wall isn't
t@@ -1790,29 +1850,42 @@ __global__ void firstDarcySolution(
if ((bc_bot == 0 && z == 0) || (bc_top == 0 && z == nz-1)
|| (z >= wall0_iz && bc_top == 0)
|| (bc_xn == 0 && x == 0) || (bc_xp == 0 && x == nx-1)
- || (bc_yn == 0 && y == 0) || (bc_yp == 0 && y == nx-1))
+ || (bc_yn == 0 && y == 0) || (bc_yp == 0 && y == nx-1)
+ || p_constant == 1)
dp_expl = 0.0;
#ifdef REPORT_FORCING_TERMS
- const Float dp_diff = (ndem*devC_dt)/(beta_f*phi*mu)
- *(k*laplace_p + dot(grad_k, grad_p));
- const Float dp_forc = -dphi/(beta_f*phi*(1.0 - phi));
+ const Float dp_diff =
+ ndem*devC_dt/(beta_f*phi*mu)
+ //*(k*laplace_p + dot(grad_k, grad_p));
+ *div_k_grad_p;
+ //const Float dp_forc = -dphi/(beta_f*phi*(1.0 - phi));
//const Float dp_forc = -div_v_p/(beta_f*phi);
+ const Float dp_forc =
+ -(ndem*devC_dt/(beta_f*phi*(1.0 - phi)))
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi));
+
printf("\n%d,%d,%d firstDarcySolution\n"
"p = %e\n"
"p_x = %e, %e\n"
"p_y = %e, %e\n"
"p_z = %e, %e\n"
"dp_expl = %e\n"
- "laplace_p = %e\n"
- "grad_p = %e, %e, %e\n"
- "grad_k = %e, %e, %e\n"
+ //"laplace_p = %e\n"
+ //"grad_p = %e, %e, %e\n"
+ //"grad_k = %e, %e, %e\n"
+ "div_k_grad_p= %e\n"
"dp_diff = %e\n"
"dp_forc = %e\n"
//"div_v_p = %e\n"
+ "phi = %e\n"
"dphi = %e\n"
"dphi/dt = %e\n"
"vp_avg = %e, %e, %e\n"
+ "grad_phi = %e, %e, %e\n"
+ //"ndem*dt/(beta*phi*(1-phi)) = %e\n"
+ //"dphi/(ndem*dt) = %e\n"
+ //"dot(v_p,grad_phi) = %e\n"
,
x,y,z,
p,
t@@ -1820,14 +1893,21 @@ __global__ void firstDarcySolution(
p_yn, p_yp,
p_zn, p_zp,
dp_expl,
- laplace_p,
- grad_p.x, grad_p.y, grad_p.z,
- grad_k.x, grad_k.y, grad_k.z,
- dp_diff, dp_forc,
+ //laplace_p,
+ //grad_p.x, grad_p.y, grad_p.z,
+ //grad_k.x, grad_k.y, grad_k.z,
+ div_k_grad_p,
+ dp_diff,
+ dp_forc,
//div_v_p//,
+ phi,
dphi,
dphi/(ndem*devC_dt),
- vp_avg.x, vp_avg.y, vp_avg.z
+ vp_avg.x, vp_avg.y, vp_avg.z,
+ grad_phi_central.x, grad_phi_central.y, grad_phi_central.z//,
+ //ndem*devC_dt/(beta_f*phi*(1.0 - phi)),
+ //dphi/(ndem*devC_dt),
+ //dot(vp_avg, grad_phi)
);
#endif
t@@ -1864,6 +1944,7 @@ __global__ void updateDarcySolution(
const int bc_top, // in
const unsigned int ndem, // in
const unsigned int wall0_iz, // in
+ const int* __restrict__ dev_darcy_p_constant, // in
Float* __restrict__ dev_darcy_p_new, // out
Float* __restrict__ dev_darcy_norm) // out
{
t@@ -1889,8 +1970,8 @@ __global__ void updateDarcySolution(
// read values
__syncthreads();
- const Float k = dev_darcy_k[cellidx];
- const Float3 grad_k = dev_darcy_grad_k[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)];
t@@ -1901,17 +1982,26 @@ __global__ void updateDarcySolution(
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 int p_constant = dev_darcy_p_constant[cellidx];
const Float p_old = dev_darcy_p_old[cellidx];
const Float dp_expl = dev_darcy_dp_expl[cellidx];
- Float p_xn = dev_darcy_p[d_idx(x-1,y,z)];
- const Float p = dev_darcy_p[cellidx];
- Float p_xp = dev_darcy_p[d_idx(x+1,y,z)];
- Float p_yn = dev_darcy_p[d_idx(x,y-1,z)];
- Float p_yp = dev_darcy_p[d_idx(x,y+1,z)];
- Float p_zn = dev_darcy_p[d_idx(x,y,z-1)];
- Float p_zp = dev_darcy_p[d_idx(x,y,z+1)];
+ Float p_xn = dev_darcy_p[d_idx(x-1,y,z)];
+ const Float p = dev_darcy_p[cellidx];
+ Float p_xp = dev_darcy_p[d_idx(x+1,y,z)];
+ Float p_yn = dev_darcy_p[d_idx(x,y-1,z)];
+ Float p_yp = dev_darcy_p[d_idx(x,y+1,z)];
+ Float p_zn = dev_darcy_p[d_idx(x,y,z-1)];
+ Float p_zp = dev_darcy_p[d_idx(x,y,z+1)];
+
+ const Float k_xn = dev_darcy_k[d_idx(x-1,y,z)];
+ const Float k = dev_darcy_k[cellidx];
+ const Float k_xp = dev_darcy_k[d_idx(x+1,y,z)];
+ const Float k_yn = dev_darcy_k[d_idx(x,y-1,z)];
+ const Float k_yp = dev_darcy_k[d_idx(x,y+1,z)];
+ const Float k_zn = dev_darcy_k[d_idx(x,y,z-1)];
+ const Float k_zp = dev_darcy_k[d_idx(x,y,z+1)];
// Neumann BCs
if (x == 0 && bc_xn == 1)
t@@ -1927,45 +2017,123 @@ __global__ void updateDarcySolution(
if (z == nz-1 && bc_top == 1)
p_zp = p;
- // upwind coefficients for grad(p) determined from values of k
- // k = 1.0: backwards difference
- // k = -1.0: forwards difference
- /*const Float3 e_k = MAKE_FLOAT3(
- copysign(1.0, grad_k.x),
- copysign(1.0, grad_k.y),
- copysign(1.0, grad_k.z));
-
- // gradient approximated by first-order forward differences
- const Float3 grad_p = MAKE_FLOAT3(
- ((1.0 + e_k.x)*(p - p_xn) + (1.0 - e_k.x)*(p_xp - p))/(dx + dx),
- ((1.0 + e_k.y)*(p - p_yn) + (1.0 - e_k.y)*(p_yp - p))/(dy + dy),
- ((1.0 + e_k.z)*(p - p_zn) + (1.0 - e_k.z)*(p_zp - p))/(dz + dz)
- );*/
-
- // gradient approximated by first-order central differences
- const Float3 grad_p = MAKE_FLOAT3(
+ /*
+ // gradients approximated by first-order central differences, order of
+ // approximation is O(dx*dx)
+ const Float3 grad_p_central = MAKE_FLOAT3(
(p_xp - p_xn)/(dx + dx),
(p_yp - p_yn)/(dy + dy),
(p_zp - p_zn)/(dz + dz));
- const Float3 grad_phi = MAKE_FLOAT3(
+ const Float3 grad_k_central = MAKE_FLOAT3(
+ (k_xp - k_xn)/(dx + dx),
+ (k_yp - k_yn)/(dy + dy),
+ (k_zp - k_zn)/(dz + dz));
+ */
+
+ const Float3 grad_phi_central = MAKE_FLOAT3(
(phi_xp - phi_xn)/(dx + dx),
(phi_yp - phi_yn)/(dy + dy),
(phi_zp - phi_zn)/(dz + dz));
+ /*
+ // upwind coefficients for grad(p) determined from values of k
+ // e_p = 1.0: backwards difference
+ // e_p = -1.0: forwards difference
+ const Float3 e_p = MAKE_FLOAT3(
+ copysign(1.0, -(p_xp - p_xn)),
+ copysign(1.0, -(p_yp - p_yn)),
+ copysign(1.0, -(p_zp - p_zn)));
+
+ // gradient approximated by first-order forward differences, order of
+ // approximation is O(dx)
+ const Float3 grad_p_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(p - p_xn) + (1.0 - e_p.x)*(p_xp - p))/dx,
+ ((1.0 + e_p.y)*(p - p_yn) + (1.0 - e_p.y)*(p_yp - p))/dy,
+ ((1.0 + e_p.z)*(p - p_zn) + (1.0 - e_p.z)*(p_zp - p))/dz
+ );
+
+ const Float3 grad_k_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(k - k_xn) + (1.0 - e_p.x)*(k_xp - k))/dx,
+ ((1.0 + e_p.y)*(k - k_yn) + (1.0 - e_p.y)*(k_yp - k))/dy,
+ ((1.0 + e_p.z)*(k - k_zn) + (1.0 - e_p.z)*(k_zp - k))/dz
+ );
+
+ const Float3 grad_phi_upwind = MAKE_FLOAT3(
+ ((1.0 + e_p.x)*(phi - phi_xn)
+ + (1.0 - e_p.x)*(phi_xp - phi))/dx,
+ ((1.0 + e_p.y)*(phi - phi_yn)
+ + (1.0 - e_p.y)*(phi_yp - phi))/dy,
+ ((1.0 + e_p.z)*(phi - phi_zn)
+ + (1.0 - e_p.z)*(phi_zp - phi))/dz
+ );
+
+ // Average central and upwind discretizations to get intermediate order
+ // of approximation
+ Float gamma = 1.0; // in [0;1], where 0: fully central, 1: fully upwind
+
+ const Float3 grad_p = MAKE_FLOAT3(
+ gamma*grad_p_upwind.x,
+ gamma*grad_p_upwind.y,
+ gamma*grad_p_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_p_central.x,
+ (1.0 - gamma) * grad_p_central.y,
+ (1.0 - gamma) * grad_p_central.z);
+
+ const Float3 grad_k = MAKE_FLOAT3(
+ gamma*grad_k_upwind.x,
+ gamma*grad_k_upwind.y,
+ gamma*grad_k_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_k_central.x,
+ (1.0 - gamma) * grad_k_central.y,
+ (1.0 - gamma) * grad_k_central.z);
+
+ const Float3 grad_phi = MAKE_FLOAT3(
+ gamma*grad_phi_upwind.x,
+ gamma*grad_phi_upwind.y,
+ gamma*grad_phi_upwind.z) + MAKE_FLOAT3(
+ (1.0 - gamma) * grad_phi_central.x,
+ (1.0 - gamma) * grad_phi_central.y,
+ (1.0 - gamma) * grad_phi_central.z);
+
// laplacian approximated by second-order central differences
const Float laplace_p =
(p_xp - (p + p) + p_xn)/(dx*dx) +
(p_yp - (p + p) + p_yn)/(dy*dy) +
(p_zp - (p + p) + p_zn)/(dz*dz);
+ */
+
+ // Solve div(k*grad(p)) as a single term, using harmonic mean for
+ // permeability. k_HM*grad(p) is found between the pressure nodes.
+ const Float div_k_grad_p =
+ (2.*k_xp*k/(k_xp + k) *
+ (p_xp - p)/dx
+ -
+ 2.*k_xn*k/(k_xn + k) *
+ (p - p_xn)/dx)/dx
+ +
+ (2.*k_yp*k/(k_yp + k) *
+ (p_yp - p)/dy
+ -
+ 2.*k_yn*k/(k_yn + k) *
+ (p - p_yn)/dy)/dy
+ +
+ (2.*k_zp*k/(k_zp + k) *
+ (p_zp - p)/dz
+ -
+ 2.*k_zn*k/(k_zn + k) *
+ (p - p_zn)/dz)/dz;
+
//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))
+ ndem*devC_dt/(beta_f*phi*mu)*div_k_grad_p
//- div_v_p/(beta_f*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));
+ -(ndem*devC_dt/(beta_f*phi*(1.0 - phi)))
+ //*(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi));
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi_central));
// Dirichlet BC at fixed-pressure boundaries and at the dynamic top
// wall. wall0_iz will be larger than the grid if the wall isn't
t@@ -1973,7 +2141,8 @@ __global__ void updateDarcySolution(
if ((bc_bot == 0 && z == 0) || (bc_top == 0 && z == nz-1)
|| (z >= wall0_iz && bc_top == 0)
|| (bc_xn == 0 && x == 0) || (bc_xp == 0 && x == nx-1)
- || (bc_yn == 0 && y == 0) || (bc_yp == 0 && y == nx-1))
+ || (bc_yn == 0 && y == 0) || (bc_yp == 0 && y == nx-1)
+ || p_constant == 1)
dp_impl = 0.0;
//p_new = p;
t@@ -1995,7 +2164,12 @@ __global__ void updateDarcySolution(
#ifdef REPORT_FORCING_TERMS_JACOBIAN
const Float dp_diff = (ndem*devC_dt)/(beta_f*phi*mu)
*(k*laplace_p + dot(grad_k, grad_p));
- const Float dp_forc = -dphi/(beta_f*phi*(1.0 - phi));
+ //const Float dp_forc = -dphi/(beta_f*phi*(1.0 - phi));
+ //const Float dp_forc = -div_v_p/(beta_f*phi);
+ const Float dp_forc =
+ -(ndem*devC_dt/(beta_f*phi*(1.0 - phi)))
+ *(dphi/(ndem*devC_dt) + dot(vp_avg, grad_phi));
+ //const Float dp_forc = -dphi/(beta_f*phi*(1.0 - phi));
//const Float dp_forc = -div_v_p/(beta_f*phi);
printf("\n%d,%d,%d updateDarcySolution\n"
"p_new = %e\n"
t@@ -2005,9 +2179,10 @@ __global__ void updateDarcySolution(
"p_z = %e, %e\n"
"dp_expl = %e\n"
"p_old = %e\n"
- "laplace_p = %e\n"
- "grad_p = %e, %e, %e\n"
- "grad_k = %e, %e, %e\n"
+ //"laplace_p = %e\n"
+ //"grad_p = %e, %e, %e\n"
+ //"grad_k = %e, %e, %e\n"
+ "div_k_grad_p= %e\n"
"dp_diff = %e\n"
"dp_forc = %e\n"
"div_v_p = %e\n"
t@@ -2022,10 +2197,12 @@ __global__ void updateDarcySolution(
p_zn, p_zp,
dp_expl,
p_old,
- laplace_p,
- grad_p.x, grad_p.y, grad_p.z,
- grad_k.x, grad_k.y, grad_k.z,
- dp_diff, dp_forc,
+ //laplace_p,
+ //grad_p.x, grad_p.y, grad_p.z,
+ //grad_k.x, grad_k.y, grad_k.z,
+ div_k_grad_p,
+ dp_diff,
+ dp_forc,
//div_v_p,
dphi, dphi/(ndem*devC_dt),
res_norm); //
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -136,6 +136,7 @@ struct NavierStokes {
int free_slip_top; // 0: no, 1: yes
Float bc_bot_flux; // Flux normal to boundary
Float bc_top_flux; // Flux normal to boundary
+ int* p_constant; // Keep pressure in cell constant (0: False, 1:True)
Float gamma; // Solver parameter: Smoothing
Float theta; // Solver parameter: Under-relaxation
Float beta; // Solver parameter: Solution method
t@@ -176,6 +177,7 @@ struct Darcy {
int free_slip_top; // 0: no, 1: yes
Float bc_bot_flux; // Flux normal to boundary
Float bc_top_flux; // Flux normal to boundary
+ int* p_constant; // Keep pressure in cell constant (0: False, 1:True)
Float tolerance; // Solver parameter: Max residual tolerance
unsigned int maxiter; // Solver parameter: Max iterations to perform
unsigned int ndem; // Solver parameter: DEM time steps per CFD step
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -5,6 +5,8 @@
#include <cstdio>
#include <cuda.h>
#include <helper_math.h>
+#include <iomanip>
+#include <time.h>
#include "vector_arithmetic.h" // for arbitrary prec. vectors
//#include <vector_functions.h> // for single prec. vectors
t@@ -1044,7 +1046,7 @@ __host__ void DEM::startTime()
// If the grid is adaptive, readjust the grid height to equal the
// positions of the dynamic walls
- if (grid.adaptive == 1) {
+ if (grid.adaptive == 1 && walls.nw > 0) {
updateGridSize();
}
t@@ -2037,17 +2039,6 @@ __host__ void DEM::startTime()
cudaThreadSynchronize();
}
- // copy porosities to the upper Z boundary
- /*if (grid.adaptive == 1) {
- copyDarcyPorositiesToTop<<<dimGridFluid,
- dimBlockFluid>>>(
- dev_darcy_phi,
- dev_darcy_dphi,
- dev_darcy_div_v_p,
- dev_darcy_vp_avg);
- cudaThreadSynchronize();
- }*/
-
// 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) {
t@@ -2238,6 +2229,7 @@ __host__ void DEM::startTime()
darcy.bc_top,
darcy.ndem,
wall0_iz,
+ dev_darcy_p_constant,
dev_darcy_dp_expl);
cudaThreadSynchronize();
if (PROFILING == 1)
t@@ -2271,6 +2263,7 @@ __host__ void DEM::startTime()
darcy.bc_top,
darcy.ndem,
wall0_iz,
+ dev_darcy_p_constant,
dev_darcy_p_new,
dev_darcy_norm);
cudaThreadSynchronize();
t@@ -2556,11 +2549,17 @@ __host__ void DEM::startTime()
// Real time it takes to compute a second of model time
t_ratio = time_spent/(time.current - t_start);
+ time_t estimated_seconds_left(t_ratio*(time.total - time.current));
+ tm *time_eta = gmtime(&estimated_seconds_left);
- cout << "\r Current simulation time: "
- << time.current << "/"
+ cout << "\r Current time: " << time.current << "/"
<< time.total << " s. ("
- << t_ratio << " s_real/s_sim) "; // << std::flush;
+ << t_ratio << " s_real/s_sim, ETA: "
+ << time_eta->tm_yday << "d "
+ << std::setw(2) << std::setfill('0') << time_eta->tm_hour << ":"
+ << std::setw(2) << std::setfill('0') << time_eta->tm_min << ":"
+ << std::setw(2) << std::setfill('0') << time_eta->tm_sec
+ << ") "; // << std::flush;
}
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -204,7 +204,7 @@ void DEM::readbin(const char *target)
// Wall mass (x), velocity (y), force (z), and deviatoric stress (w)
walls.nx = new Float4[walls.nw];
walls.mvfd = new Float4[walls.nw];
- walls.tau_x = new Float[walls.nw];
+ walls.tau_x = new Float[1];
for (i = 0; i<walls.nw; ++i)
ifs.read(as_bytes(walls.wmode[i]), sizeof(walls.wmode[i]));
t@@ -303,6 +303,12 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.bc_bot_flux), sizeof(Float));
ifs.read(as_bytes(ns.bc_top_flux), sizeof(Float));
+ for (z = 0; z<grid.num[2]; ++z)
+ for (y = 0; y<grid.num[1]; ++y)
+ for (x = 0; x<grid.num[0]; ++x)
+ ifs.read(as_bytes(ns.p_constant[idx(x,y,z)]),
+ sizeof(int));
+
ifs.read(as_bytes(ns.gamma), sizeof(Float));
ifs.read(as_bytes(ns.theta), sizeof(Float));
ifs.read(as_bytes(ns.beta), sizeof(Float));
t@@ -377,6 +383,12 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(darcy.bc_bot_flux), sizeof(Float));
ifs.read(as_bytes(darcy.bc_top_flux), sizeof(Float));
+ for (z = 0; z<darcy.nz; ++z)
+ for (y = 0; y<darcy.ny; ++y)
+ for (x = 0; x<darcy.nx; ++x)
+ ifs.read(as_bytes(darcy.p_constant[d_idx(x,y,z)]),
+ sizeof(int));
+
ifs.read(as_bytes(darcy.tolerance), sizeof(Float));
ifs.read(as_bytes(darcy.maxiter), sizeof(unsigned int));
ifs.read(as_bytes(darcy.ndem), sizeof(unsigned int));
t@@ -603,6 +615,12 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(ns.bc_bot_flux), sizeof(Float));
ofs.write(as_bytes(ns.bc_top_flux), sizeof(Float));
+ for (z = 0; z<ns.nz; ++z)
+ for (y = 0; y<ns.ny; ++y)
+ for (x = 0; x<ns.nx; ++x)
+ ofs.write(as_bytes(ns.p_constant[idx(x,y,z)]),
+ sizeof(int));
+
ofs.write(as_bytes(ns.gamma), sizeof(Float));
ofs.write(as_bytes(ns.theta), sizeof(Float));
ofs.write(as_bytes(ns.beta), sizeof(Float));
t@@ -678,6 +696,12 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(darcy.bc_bot_flux), sizeof(Float));
ofs.write(as_bytes(darcy.bc_top_flux), sizeof(Float));
+ for (z = 0; z<darcy.nz; ++z)
+ for (y = 0; y<darcy.ny; ++y)
+ for (x = 0; x<darcy.nx; ++x)
+ ofs.write(as_bytes(darcy.p_constant[d_idx(x,y,z)]),
+ sizeof(int));
+
ofs.write(as_bytes(darcy.tolerance), sizeof(Float));
ofs.write(as_bytes(darcy.maxiter), sizeof(unsigned int));
ofs.write(as_bytes(darcy.ndem), sizeof(unsigned int));
(DIR) diff --git a/src/navierstokes.cpp b/src/navierstokes.cpp
t@@ -41,6 +41,7 @@ void DEM::initNSmem()
ns.f_p = new Float4[np]; // pressure force on particles
ns.f_v = new Float4[np]; // viscous force on particles
ns.f_sum = new Float4[np]; // sum of fluid forces on particles
+ ns.p_constant = new int[ncells]; // unused
}
unsigned int DEM::NScells()
t@@ -81,6 +82,7 @@ void DEM::freeNSmem()
delete[] ns.f_p;
delete[] ns.f_v;
delete[] ns.f_sum;
+ delete[] ns.p_constant;
}
// 3D index to 1D index
(DIR) diff --git a/src/sphere.h b/src/sphere.h
t@@ -317,6 +317,7 @@ class DEM {
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
+ int* dev_darcy_p_constant; // Constant pressure (0: False, 1: True)
// Darcy functions
void initDarcyMem();
(DIR) diff --git a/src/version.h b/src/version.h
t@@ -2,6 +2,6 @@
#define VERSION_H_
// Define source code version
-const double VERSION = 2.14;
+const double VERSION = 2.15;
#endif