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 9bc372451003f70474ea166830f10cf2a4a1763c
(DIR) parent d8d27df5b9343bd74b9cebfaf9d707252fdd7188
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Sun, 6 Sep 2015 08:47:22 +0200
Merge branch 'master' of github.com:anders-dc/sphere
Diffstat:
A output/shear-sigma0=10000.0.output… | 0
A output/shear-sigma0=10000.0.status… | 1 +
A python/halfshear-darcy-rate-streng… | 254 +++++++++++++++++++++++++++++++
M python/halfshear-darcy-strength-di… | 9 +++++----
M python/halfshear-darcy-strength-di… | 2 +-
A python/sigma-sim1-starter.py | 182 +++++++++++++++++++++++++++++++
M python/sphere.py | 71 +++++++++++++++++++++++++++----
M src/darcy.cuh | 48 +++++++++++++++++++++++++++++++
M src/datatypes.h | 8 ++++++--
M src/device.cu | 41 +++++++++++++++++++++++++++++++
M src/file_io.cpp | 8 ++++++++
M src/version.h | 2 +-
12 files changed, 609 insertions(+), 17 deletions(-)
---
(DIR) diff --git a/output/shear-sigma0=10000.0.output01999.bin b/output/shear-sigma0=10000.0.output01999.bin
Binary files differ.
(DIR) diff --git a/output/shear-sigma0=10000.0.status.dat b/output/shear-sigma0=10000.0.status.dat
t@@ -0,0 +1 @@
+1.9990e+01 9.9950e+01 1999
(DIR) diff --git a/python/halfshear-darcy-rate-strength.py b/python/halfshear-darcy-rate-strength.py
t@@ -0,0 +1,254 @@
+#!/usr/bin/env python
+import matplotlib
+matplotlib.use('Agg')
+matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
+matplotlib.rc('text', usetex=True)
+matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
+import shutil
+
+import os
+import sys
+import numpy
+import sphere
+from permeabilitycalculator import *
+import matplotlib.pyplot as plt
+import scipy.optimize
+
+import seaborn as sns
+#sns.set(style='ticks', palette='Set2')
+sns.set(style='ticks', palette='colorblind')
+#sns.set(style='ticks', palette='muted')
+#sns.set(style='ticks', palette='pastel')
+sns.despine() # remove right and top spines
+
+pressures = True
+zflow = False
+contact_forces = False
+smooth_friction = True
+
+
+sids = [
+ 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-06-velfac=1.0-shear',
+ 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-14-mu=1.797e-06-velfac=1.0-shear',
+ 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-06-velfac=1.0-shear',
+
+ #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=3.594e-07-velfac=1.0-shear',
+ #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=6.11e-07-velfac=1.0-shear',
+
+ 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-07-velfac=1.0-shear',
+ 'halfshear-darcy-sigma0=20000.0-k_c=3.5e-15-mu=1.797e-08-velfac=1.0-shear',
+
+
+ #'halfshear-sigma0=20000.0-shear'
+
+ # from halfshear-darcy-perm.sh
+ #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-07-velfac=1.0-shear',
+ #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-08-velfac=1.0-shear',
+ #'halfshear-darcy-sigma0=20000.0-k_c=3.5e-13-mu=1.797e-09-velfac=1.0-shear',
+ ]
+fluids = [
+ True,
+ True,
+ True,
+ True,
+ True,
+ #False,
+ True,
+ True,
+ True,
+ ]
+
+def smooth(x, window_len=10, window='hanning'):
+ """smooth the data using a window with requested size.
+
+ This method is based on the convolution of a scaled window with the signal.
+ The signal is prepared by introducing reflected copies of the signal
+ (with the window size) in both ends so that transient parts are minimized
+ in the begining and end part of the output signal.
+
+ input:
+ x: the input signal
+ window_len: the dimension of the smoothing window
+ window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
+ flat window will produce a moving average smoothing.
+
+ output:
+ the smoothed signal
+
+ example:
+
+ import numpy as np
+ t = np.linspace(-2,2,0.1)
+ x = np.sin(t)+np.random.randn(len(t))*0.1
+ y = smooth(x)
+
+ see also:
+
+ numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
+ scipy.signal.lfilter
+
+ TODO: the window parameter could be the window itself if an array instead of a string
+ """
+
+ if x.ndim != 1:
+ raise ValueError, "smooth only accepts 1 dimension arrays."
+
+ if x.size < window_len:
+ raise ValueError, "Input vector needs to be bigger than window size."
+
+ if window_len < 3:
+ return x
+
+ if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
+ raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
+
+ s=numpy.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
+ #print(len(s))
+
+ if window == 'flat': #moving average
+ w = numpy.ones(window_len,'d')
+ else:
+ w = getattr(numpy, window)(window_len)
+ y = numpy.convolve(w/w.sum(), s, mode='same')
+ return y[window_len-1:-window_len+1]
+
+
+max_step = 400
+friction = numpy.zeros(max_step)
+
+velratios = []
+peakfrictions = []
+
+sim = sphere.sim(sids[0], fluid=True)
+sim.readfirst(verbose=False)
+fig = plt.figure(figsize=(3.5, 3))
+
+if False:
+ it=0
+ for sid in sids:
+
+ print '\n' + sid
+ sim.id(sid)
+ sim.fluid=fluids[it]
+ it += 1
+
+ sim.readfirst(verbose=False)
+
+ velratio = 1.0
+ if sim.fluid:
+ if sim.k_c[0] > 3.5e-15 and sim.mu[0] < 1.797e-06:
+ print 'Large permeability and low viscosity'
+ velratio = 3.5e-15/sim.k_c[0] \
+ * sim.mu[0]/1.797e-06
+ elif sim.k_c[0] > 3.5e-15:
+ print 'Large permeability'
+ velratio = 3.5e-15/sim.k_c[0]
+ elif sim.mu < 1.797e-06:
+ print 'Low viscosity'
+ velratio = sim.mu[0]/1.797e-06
+ elif numpy.isclose(sim.k_c[0], 3.5e-15) and \
+ numpy.isclose(sim.mu[0], 1.797e-06):
+ print 'Normal permeability and viscosity'
+ velratio = 1.
+ else:
+ raise Exception('Could not determine velratio')
+ else:
+ velratio = 0.001
+
+ print 'velratio = ' + str(velratio)
+
+ for i in numpy.arange(max_step):
+ sim.readstep(i+1, verbose=False)
+
+ friction[i] = sim.shearStress('effective')/\
+ sim.currentNormalStress('defined')
+
+ smoothfriction = smooth(friction, 20, window='hanning')
+ peakfriction = numpy.amax(smoothfriction[14:])
+
+ plt.plot(numpy.arange(max_step), friction, alpha=0.3, color='b')
+ plt.plot(numpy.arange(max_step), smoothfriction, color='b')
+ plt.title(str(peakfriction))
+ plt.savefig(sid + '-friction.pdf')
+ plt.clf()
+
+ velratios.append(velratio)
+ peakfrictions.append(peakfriction)
+else:
+ velratios =\
+ [0.01,
+ 0.099999999999999992,
+ 1.0,
+ 0.10000000000000001,
+ #0.010000000000000002,
+ #0.001,
+ #0.0001,
+ #1.0000000000000001e-05
+ ]
+ peakfrictions = \
+ [0.619354807290315,
+ 0.6536161052814875,
+ 0.70810354077280957,
+ 0.64649301787774571,
+ #0.65265739261434697,
+ #0.85878138368962764,
+ #1.6263903846405066,
+ #0.94451353171977692
+ ]
+
+
+#def frictionfunction(velratio, a, b=numpy.min(peakfrictions)):
+def frictionfunction(x, a, b):
+ #return numpy.exp(a*velratio) + numpy.min(peakfrictions)
+ #return -numpy.log(a*velratio) + numpy.min(peakfrictions)
+ #return a*numpy.log(velratio) + b
+ return a*10**(x) + b
+
+popt, pvoc = scipy.optimize.curve_fit(
+ frictionfunction,
+ peakfrictions,
+ velratios)
+print popt
+print pvoc
+a=popt[0]
+b=popt[1]
+
+#a = 0.025
+#b = numpy.min(peakfrictions)
+#b = numpy.max(peakfrictions)
+
+shearvel = sim.shearVel()/1000.*numpy.asarray(velratios) * (3600.*24.*365.25)
+
+'''
+plt.semilogx(
+ #[1, numpy.min(shearvel)],
+ [1, 6],
+ #[1, 1000],
+ #[numpy.min(peakfrictions), numpy.min(peakfrictions)],
+ [0.615, 0.615],
+ '--', color='gray', linewidth=2)
+'''
+
+#plt.semilogx(velratios, peakfrictions, 'o')
+#plt.semilogx(shearvel, peakfrictions, 'o')
+plt.semilogx(shearvel, peakfrictions, 'o')
+#plt.plot(velratios, peakfrictions, 'o')
+
+xfit = numpy.linspace(numpy.min(velratios), numpy.max(velratios))
+#yfit = frictionfunction(xfit, popt[0], popt[1])
+yfit = frictionfunction(xfit, a, b)
+#plt.semilogx(xfit, yfit)
+#plt.plot(xfit, yfit)
+
+plt.xlabel('Shear velocity [m a$^{-1}$]')
+plt.ylabel('Peak shear friction, $\\max(\\tau/\\sigma_0)$ [-]')
+#plt.xlim([0.8, 110])
+#plt.xlim([0.008, 1.1])
+plt.xlim([1., 1000,])
+plt.ylim([0.6, 0.72])
+plt.tight_layout()
+sns.despine() # remove right and top spines
+filename = 'halfshear-darcy-rate-strength.pdf'
+plt.savefig(filename)
+plt.close()
+print(filename)
(DIR) diff --git a/python/halfshear-darcy-strength-dilation-rate.py b/python/halfshear-darcy-strength-dilation-rate.py
t@@ -219,11 +219,11 @@ for sigma0 in sigma0_list:
mu_f = mu_f_vals[c]
if numpy.isclose(mu_f, 1.797e-6):
- label = 'ref. shear velocity'
+ label = 'ref.\\ shear velocity'
#elif numpy.isclose(mu_f, 1.204e-6):
- #label = 'ref. shear velocity$\\times$0.67'
+ #label = 'ref.\\ shear velocity$\\times$0.67'
#elif numpy.isclose(mu_f, 1.797e-8):
- #label = 'ref. shear velocity$\\times$0.01'
+ #label = 'ref.\\ shear velocity$\\times$0.01'
else:
#label = '$\\mu_\\text{{f}}$ = {:.3e} Pa s'.format(mu_f)
label = 'shear velocity$\\times${:.2}'.format(mu_f/mu_f_vals[0])
t@@ -428,7 +428,8 @@ for sigma0 in sigma0_list:
#ax2.set_xlim([0.0, 0.2])
#ax1.set_ylim([-7, 45])
- ax1.set_xlim([0.0, 1.0])
+ #ax1.set_xlim([0.0, 1.0])
+ ax1.set_xlim([0.0, 2.0])
ax1.set_ylim([0.0, 1.0])
ax2.set_ylim([0.0, 1.0])
ax3.set_ylim([-150., 150.])
(DIR) diff --git a/python/halfshear-darcy-strength-dilation.py b/python/halfshear-darcy-strength-dilation.py
t@@ -220,7 +220,7 @@ for sigma0 in sigma0_list:
elif numpy.isclose(k_c, 3.5e-13, atol=1.0e-16):
label = 'high permeability'
elif numpy.isclose(k_c, 3.5e-14, atol=1.0e-16):
- label = 'interm. permeability'
+ label = 'interm.\\ permeability'
elif numpy.isclose(k_c, 3.5e-15, atol=1.0e-16):
label = 'low permeability'
else:
(DIR) diff --git a/python/sigma-sim1-starter.py b/python/sigma-sim1-starter.py
t@@ -0,0 +1,182 @@
+#!/usr/bin/env python
+import sphere
+import numpy
+import sys
+
+# launch with:
+# $ ipython sigma-sim1-starter.py <device> <fluid> <c_phi> <k_c> <sigma_0> <mu> <velfac>
+
+# start with
+# ipython sigma-sim1-starter.py 0 1 1.0 2.0e-16 10000.0 2.080e-7 1.0
+
+device = int(sys.argv[1])
+wet = int(sys.argv[2])
+c_phi = float(sys.argv[3])
+k_c = float(sys.argv[4])
+sigma0 = float(sys.argv[5])
+mu = float(sys.argv[6])
+velfac = float(sys.argv[7])
+
+if wet == 1:
+ fluid = True
+else:
+ fluid = False
+
+#sim = sphere.sim('halfshear-sigma0=' + str(sigma0), fluid=False)
+sim = sphere.sim('shear-sigma0=' + str(sigma0), fluid=False)
+sim.readlast()
+#sim.readbin('../input/shear-sigma0=10000.0-new.bin')
+#sim.scaleSize(0.01) # from 1 cm to 0.01 cm = 100 micro m (fine sand)
+
+sim.fluid = fluid
+if fluid:
+ #sim.id('halfshear-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
+ #'-mu=' + str(mu) + '-velfac=' + str(velfac) + '-shear')
+ sim.id('s1-darcy-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
+ '-mu=' + str(mu) + '-velfac=' + str(velfac) + '-noflux-shear')
+else:
+ sim.id('s1-darcy-sigma0=' + str(sigma0) + '-velfac=' + str(velfac) + \
+ '-noflux-shear')
+
+#sim.checkerboardColors(nx=6,ny=3,nz=6)
+sim.checkerboardColors(nx=6,ny=6,nz=6)
+sim.cleanup()
+sim.adjustUpperWall()
+sim.zeroKinematics()
+
+# customized shear function for linear velocity gradient
+def shearVelocityGradient(sim, shear_strain_rate = 1.0, shear_stress = False):
+ '''
+ Setup shear experiment either by a constant shear rate or a constant
+ shear stress. The shear strain rate is the shear velocity divided by
+ the initial height per second. The shear movement is along the positive
+ x axis. The function zeroes the tangential wall viscosity (gamma_wt) and
+ the wall friction coefficients (mu_ws, mu_wn).
+
+ :param shear_strain_rate: The shear strain rate [-] to use if
+ shear_stress isn't False.
+ :type shear_strain_rate: float
+ :param shear_stress: The shear stress value to use [Pa].
+ :type shear_stress: float or bool
+ '''
+
+ sim.nw[0] = 1
+
+ # Find lowest and heighest point
+ z_min = numpy.min(sim.x[:,2] - sim.radius)
+ z_max = numpy.max(sim.x[:,2] + sim.radius)
+
+ # the grid cell size is equal to the max. particle diameter
+ cellsize = sim.L[0] / sim.num[0]
+
+ # make grid one cell heigher to allow dilation
+ sim.num[2] += 1
+ sim.L[2] = sim.num[2] * cellsize
+
+ # zero kinematics
+ sim.zeroKinematics()
+
+ # Adjust grid and placement of upper wall
+ sim.wmode = numpy.array([1])
+
+ # Fix horizontal velocity to 0.0 of lowermost particles
+ d_max_below = numpy.max(sim.radius[numpy.nonzero(sim.x[:,2] <
+ (z_max-z_min)*0.3)])*2.0
+ I = numpy.nonzero(sim.x[:,2] < (z_min + d_max_below))
+ sim.fixvel[I] = 1
+ sim.angvel[I,0] = 0.0
+ sim.angvel[I,1] = 0.0
+ sim.angvel[I,2] = 0.0
+ sim.vel[I,0] = 0.0 # x-dim
+ sim.vel[I,1] = 0.0 # y-dim
+ sim.color[I] = -1
+
+ # Fix horizontal velocity to specific value of uppermost particles
+ d_max_top = numpy.max(sim.radius[numpy.nonzero(sim.x[:,2] >
+ (z_max-z_min)*0.7)])*2.0
+ I = numpy.nonzero(sim.x[:,2] > (z_max - d_max_top))
+ sim.fixvel[I] = 1
+ sim.angvel[I,0] = 0.0
+ sim.angvel[I,1] = 0.0
+ sim.angvel[I,2] = 0.0
+ if shear_stress == False:
+ prefactor = sim.x[I,1]/sim.L[1]
+ sim.vel[I,0] = prefactor*(z_max-z_min)*shear_strain_rate
+ else:
+ sim.vel[I,0] = 0.0
+ sim.wmode[0] = 3
+ sim.w_tau_x[0] = float(shear_stress)
+ sim.vel[I,1] = 0.0 # y-dim
+ sim.color[I] = -1
+
+ # Set wall tangential viscosity to zero
+ sim.gamma_wt[0] = 0.0
+
+ # Set wall friction coefficients to zero
+ sim.mu_ws[0] = 0.0
+ sim.mu_wd[0] = 0.0
+ return sim
+
+sim = shearVelocityGradient(sim, 1.0/20.0 * velfac)
+K_q_real = 36.4e9
+K_w_real = 2.2e9
+
+
+#K_q_sim = 1.16e9
+K_q_sim = 1.16e6
+sim.setStiffnessNormal(K_q_sim)
+sim.setStiffnessTangential(K_q_sim)
+K_w_sim = K_w_real/K_q_real * K_q_sim
+
+
+if fluid:
+ #sim.num[2] *= 2
+ sim.num[:] /= 2
+ #sim.L[2] *= 2.0
+ #sim.initFluid(mu = 1.787e-6, p = 600.0e3, cfd_solver = 1)
+ sim.initFluid(mu = mu, p = 0.0, cfd_solver = 1)
+ sim.setFluidTopFixedPressure()
+ #sim.setFluidTopFixedFlow()
+ sim.setFluidBottomNoFlow()
+ #sim.setFluidBottomFixedPressure()
+ #sim.setDEMstepsPerCFDstep(10)
+ sim.setMaxIterations(2e5)
+ sim.setPermeabilityPrefactor(k_c)
+ sim.setFluidCompressibility(1.0/K_w_sim)
+
+
+# frictionless side boundaries
+sim.periodicBoundariesX()
+
+# rearrange particle assemblage to accomodate frictionless side boundaries
+sim.x[:,1] += numpy.abs(numpy.min(sim.x[:,1] - sim.radius[:]))
+sim.L[1] = numpy.max(sim.x[:,1] + sim.radius[:])
+
+
+sim.w_sigma0[0] = sigma0
+sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
+
+#sim.setStiffnessNormal(36.4e9 * 0.1 / 2.0)
+#sim.setStiffnessTangential(36.4e9/3.0 * 0.1 / 2.0)
+sim.setStiffnessNormal(K_q_sim)
+sim.setStiffnessTangential(K_q_sim)
+sim.mu_s[0] = 0.5
+sim.mu_d[0] = 0.5
+sim.setDampingNormal(0.0)
+sim.setDampingTangential(0.0)
+#sim.deleteAllParticles()
+#sim.fixvel[:] = -1.0
+
+sim.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
+#sim.time_dt[0] *= 1.0e-2
+#sim.initTemporal(total = 1.0e-4, file_dt = 1.0e-5, epsilon=0.07)
+
+# Fix lowermost particles
+#dz = sim.L[2]/sim.num[2]
+#I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
+#sim.fixvel[I] = 1
+
+sim.run(dry=True)
+
+sim.run(device=device)
+sim.writeVTKall()
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
import matplotlib.collections
matplotlib.rcParams.update({'font.size': 7, 'font.family': 'serif'})
matplotlib.rc('text', usetex=True)
-matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{amsmath}"]
+matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
from matplotlib.font_manager import FontProperties
import subprocess
import pickle as pl
t@@ -24,11 +24,12 @@ numpy.seterr(all='warn', over='raise')
# Sphere version number. This field should correspond to the value in
# `../src/constants.h`.
-VERSION = 2.1
+VERSION = 2.11
# Transparency on plot legends
legend_alpha = 0.5
+
class sim:
'''
Class containing all ``sphere`` data.
t@@ -54,12 +55,12 @@ class sim:
'''
def __init__(self,
- sid = 'unnamed',
- np = 0,
- nd = 3,
- nw = 0,
- fluid = False,
- cfd_solver = 0):
+ sid='unnamed',
+ np=0,
+ nd=3,
+ nw=0,
+ fluid=False,
+ cfd_solver=0):
# Sphere version number
self.version = numpy.ones(1, dtype=numpy.float64)*VERSION
t@@ -329,13 +330,21 @@ class sim:
self.p_mod_phi = numpy.zeros(1, dtype=numpy.float64) # Shift [rad]
# Boundary conditions at the top and bottom of the fluid grid
- # 0: Dirichlet, 1: Neumann free slip, 2: Neumann no slip, 3: Periodic
+ # 0: Dirichlet
+ # 1: Neumann free slip
+ # 2: Neumann no slip
+ # 3: Periodic
+ # 4: Constant flux (Darcy solver only)
self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
self.bc_top = numpy.zeros(1, dtype=numpy.int32)
# Free slip boundaries? 1: yes
self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
+ # Boundary-normal flux (in case of bc_* = 4)
+ self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
+ self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
+
## Solver parameters
t@@ -652,6 +661,12 @@ class sim:
elif self.free_slip_top != other.free_slip_top:
print(77)
return 77
+ elif self.bc_bot_flux != other.bc_bot_flux:
+ print(91)
+ return 91
+ elif self.bc_top_flux != other.bc_top_flux:
+ print(91)
+ return 91
if self.cfd_solver == 0:
if self.gamma != other.gamma:
t@@ -1137,6 +1152,14 @@ class sim:
numpy.fromfile(fh, dtype=numpy.int32, count=1)
self.free_slip_top =\
numpy.fromfile(fh, dtype=numpy.int32, count=1)
+ if self.version >= 2.11:
+ self.bc_bot_flux =\
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ self.bc_top_flux =\
+ numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ else:
+ self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
+ self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
if self.version >= 2.0 and self.cfd_solver == 0:
self.gamma = \
t@@ -1375,6 +1398,8 @@ class sim:
fh.write(self.bc_top.astype(numpy.int32))
fh.write(self.free_slip_bot.astype(numpy.int32))
fh.write(self.free_slip_top.astype(numpy.int32))
+ fh.write(self.bc_bot_flux.astype(numpy.float64))
+ fh.write(self.bc_top_flux.astype(numpy.float64))
# Navier Stokes
if self.cfd_solver[0] == 0:
t@@ -3278,6 +3303,8 @@ class sim:
self.bc_top = numpy.zeros(1, dtype=numpy.int32)
self.free_slip_bot = numpy.ones(1, dtype=numpy.int32)
self.free_slip_top = numpy.ones(1, dtype=numpy.int32)
+ self.bc_bot_flux = numpy.zeros(1, dtype=numpy.float64)
+ self.bc_top_flux = numpy.zeros(1, dtype=numpy.float64)
# Fluid solver type
# 0: Navier Stokes (fluid with inertia)
t@@ -3362,6 +3389,19 @@ class sim:
'''
self.bc_bot[0] = 0
+ def setFluidBottomFixedFlux(self, specific_flux):
+ '''
+ Define a constant fluid flux normal to the boundary.
+
+ The default behavior for the boundary is fixed value (Dirichlet), see
+ :func:`setFluidBottomFixedPressure()`.
+
+ :param specific_flux: Specific flux values across boundary (positive
+ values upwards), [m/s]
+ '''
+ self.bc_bot[0] = 4
+ self.bc_bot_flux[0] = specific_flux
+
def setFluidTopNoFlow(self):
'''
Set the upper boundary of the fluid domain to follow the no-flow
t@@ -3392,6 +3432,19 @@ class sim:
'''
self.bc_top[0] = 0
+ def setFluidTopFixedFlux(self, specific_flux):
+ '''
+ Define a constant fluid flux normal to the boundary.
+
+ The default behavior for the boundary is fixed value (Dirichlet), see
+ :func:`setFluidBottomFixedPressure()`.
+
+ :param specific_flux: Specific flux values across boundary (positive
+ values upwards), [m/s]
+ '''
+ self.bc_top[0] = 4
+ self.bc_top_flux[0] = specific_flux
+
def setPermeabilityPrefactor(self, k_c, verbose=True):
'''
Set the permeability prefactor from Goren et al 2011, eq. 24. The
(DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
t@@ -251,6 +251,54 @@ __global__ void setDarcyGhostNodes(
}
}
+// Update a field in the ghost nodes from their parent cell values. The edge
+// (diagonal) cells are not written since they are not read. Launch this kernel
+// for all cells in the grid using
+// setDarcyGhostNodes<datatype><<<.. , ..>>>( .. );
+ template<typename T>
+__global__ void setDarcyGhostNodesFlux(
+ T* __restrict__ dev_scalarfield, // out
+ const int bc_bot, // in
+ const int bc_top, // in
+ const Float bc_bot_flux, // in
+ const Float bc_top_flux, // in
+ const Float* __restrict__ dev_darcy_k, // in
+ const Float mu) // in
+{
+ // 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 (x < nx && y < ny && z < nz && (bc_bot == 4 || bc_top == 4)) {
+
+ const T p = dev_scalarfield[d_idx(x,y,z)];
+ const Float k = dev_darcy_k[d_idx(x,y,z)];
+ const Float dz = devC_grid.L[2]/nz;
+
+ Float q_z = 0.;
+ if (z == 0)
+ q_z = bc_bot_flux;
+ else if (z == nz-1)
+ q_z = bc_top_flux;
+
+ const Float p_ghost = -mu/k*q_z * dz + p;
+
+ // z
+ if (z == 0 && bc_bot == 4)
+ dev_scalarfield[idx(x,y,-1)] = p_ghost;
+
+ if (z == nz-1 && bc_top == 4)
+ dev_scalarfield[idx(x,y,nz)] = p_ghost;
+ }
+}
+
/*
__global__ void findDarcyParticleVelocities(
const unsigned int* __restrict__ dev_cellStart, // in
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -128,10 +128,12 @@ struct NavierStokes {
Float p_mod_A; // Pressure modulation amplitude at top
Float p_mod_f; // Pressure modulation frequency at top
Float p_mod_phi; // Pressure modulation phase at top
- int bc_bot; // 0: Dirichlet, 1: Neumann
- int bc_top; // 0: Dirichlet, 1: Neumann
+ int bc_bot; // 0: Dirichlet, 1: Neumann, 3: Periodic, 4: Flux
+ int bc_top; // 0: Dirichlet, 1: Neumann, 3: Periodic, 4: Flux
int free_slip_bot; // 0: no, 1: yes
int free_slip_top; // 0: no, 1: yes
+ Float bc_bot_flux; // Flux normal to boundary
+ Float bc_top_flux; // Flux normal to boundary
Float gamma; // Solver parameter: Smoothing
Float theta; // Solver parameter: Under-relaxation
Float beta; // Solver parameter: Solution method
t@@ -166,6 +168,8 @@ struct Darcy {
int bc_top; // 0: Dirichlet, 1: Neumann
int free_slip_bot; // 0: no, 1: yes
int free_slip_top; // 0: no, 1: yes
+ Float bc_bot_flux; // Flux normal to boundary
+ Float bc_top_flux; // Flux normal to boundary
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@@ -2000,6 +2000,26 @@ __host__ void DEM::startTime()
iter);
}
+ if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
+ if (PROFILING == 1)
+ startTimer(&kernel_tic);
+ setDarcyGhostNodesFlux<Float>
+ <<<dimGridFluid, dimBlockFluid>>>(
+ dev_darcy_p,
+ darcy.bc_bot,
+ darcy.bc_top,
+ darcy.bc_bot_flux,
+ darcy.bc_top_flux,
+ dev_darcy_k,
+ darcy.mu);
+ cudaThreadSynchronize();
+ if (PROFILING == 1)
+ stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
+ &t_setDarcyGhostNodes);
+ checkForCudaErrorsIter(
+ "Post setDarcyGhostNodesFlux", iter);
+ }
+
// Solve the system of epsilon using a Jacobi iterative
// solver. The average normalized residual is initialized
// to a large value.
t@@ -2114,6 +2134,27 @@ __host__ void DEM::startTime()
"Post setDarcyTopWallFixedFlow", iter);
}
+ if (darcy.bc_bot == 4 || darcy.bc_top == 4) {
+ if (PROFILING == 1)
+ startTimer(&kernel_tic);
+ setDarcyGhostNodesFlux<Float>
+ <<<dimGridFluid, dimBlockFluid>>>(
+ dev_darcy_p,
+ darcy.bc_bot,
+ darcy.bc_top,
+ darcy.bc_bot_flux,
+ darcy.bc_top_flux,
+ dev_darcy_k,
+ darcy.mu);
+ cudaThreadSynchronize();
+ if (PROFILING == 1)
+ stopTimer(&kernel_tic, &kernel_toc,
+ &kernel_elapsed,
+ &t_setDarcyGhostNodes);
+ checkForCudaErrorsIter(
+ "Post setDarcyGhostNodesFlux", iter);
+ }
+
// Copy new values to current values
if (PROFILING == 1)
startTimer(&kernel_tic);
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -298,6 +298,8 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(ns.bc_top), sizeof(int));
ifs.read(as_bytes(ns.free_slip_bot), sizeof(int));
ifs.read(as_bytes(ns.free_slip_top), sizeof(int));
+ ifs.read(as_bytes(ns.bc_bot_flux), sizeof(Float));
+ ifs.read(as_bytes(ns.bc_top_flux), sizeof(Float));
ifs.read(as_bytes(ns.gamma), sizeof(Float));
ifs.read(as_bytes(ns.theta), sizeof(Float));
t@@ -366,6 +368,8 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(darcy.bc_top), sizeof(int));
ifs.read(as_bytes(darcy.free_slip_bot), sizeof(int));
ifs.read(as_bytes(darcy.free_slip_top), sizeof(int));
+ ifs.read(as_bytes(darcy.bc_bot_flux), sizeof(Float));
+ ifs.read(as_bytes(darcy.bc_top_flux), sizeof(Float));
ifs.read(as_bytes(darcy.tolerance), sizeof(Float));
ifs.read(as_bytes(darcy.maxiter), sizeof(unsigned int));
t@@ -588,6 +592,8 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(ns.bc_top), sizeof(int));
ofs.write(as_bytes(ns.free_slip_bot), sizeof(int));
ofs.write(as_bytes(ns.free_slip_top), sizeof(int));
+ ofs.write(as_bytes(ns.bc_bot_flux), sizeof(Float));
+ ofs.write(as_bytes(ns.bc_top_flux), sizeof(Float));
ofs.write(as_bytes(ns.gamma), sizeof(Float));
ofs.write(as_bytes(ns.theta), sizeof(Float));
t@@ -657,6 +663,8 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(darcy.bc_top), sizeof(int));
ofs.write(as_bytes(darcy.free_slip_bot), sizeof(int));
ofs.write(as_bytes(darcy.free_slip_top), sizeof(int));
+ ofs.write(as_bytes(darcy.bc_bot_flux), sizeof(Float));
+ ofs.write(as_bytes(darcy.bc_top_flux), sizeof(Float));
ofs.write(as_bytes(darcy.tolerance), sizeof(Float));
ofs.write(as_bytes(darcy.maxiter), sizeof(unsigned int));
(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.10;
+const double VERSION = 2.11;
#endif