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 9ede2ec8e2b8f465001c16a96246af21101e61d7
(DIR) parent cdb27fe393e02eeefeca2d85bd1f584d9c0a8d2d
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 1 Apr 2016 00:14:02 +0200
Merge branch 'master' of github.com:anders-dc/sphere
Conflicts:
python/channel.py
Diffstat:
M README.rst | 4 +---
A demParticles.cpd | 184 +++++++++++++++++++++++++++++++
M output/cube-init.status.dat | 2 +-
M python/capillary-cohesion.py | 35 +++++++++++++++++++------------
M python/channel.py | 91 +++++++++++++++++++------------
M python/cube-init.py | 10 +++++-----
A python/rate-state.py | 103 +++++++++++++++++++++++++++++++
M python/sphere.py | 10 ++++++++--
M src/debug.h | 3 +++
M src/device.cu | 63 +++++++++++++++++++++++++++++++
M src/integration.cuh | 10 ++++++++++
M src/sphere.cpp | 3 +++
12 files changed, 459 insertions(+), 59 deletions(-)
---
(DIR) diff --git a/README.rst b/README.rst
t@@ -13,9 +13,6 @@ A powerful Nvidia GPU with proper support for double precision floating is
highly recommended. ``sphere`` has been tested with the Nvidia Tesla K20 and
Nvidia Tesla M2050 GPUs.
- **Note:** CUDA 6.5 is the recommended version. CUDA 7.0 may work but is not
- thoroughly tested yet.
-
License
-------
``sphere`` is licensed under the GNU General Public License, v.3.
t@@ -128,3 +125,4 @@ Author
Anders Damsgaard, `adamsgaard@ucsd.edu <mailto:adamsgaard@ucsd.edu>`_,
`webpage <https://cs.au.dk/~adc>`_.
+With contributions from Cooper Elsworth.
(DIR) diff --git a/demParticles.cpd b/demParticles.cpd
t@@ -0,0 +1,184 @@
+<CustomFilterDefinitions>
+ <CustomProxyDefinition name="demParticles" group="filters">
+ <CompoundSourceProxy id="4921" servers="1">
+ <Proxy group="filters" type="Glyph" id="4758" servers="1" compound_name="Glyph1">
+ <Property name="GlyphMode" id="4758.GlyphMode" number_of_elements="1">
+ <Element index="0" value="0"/>
+ <Domain name="enum" id="4758.GlyphMode.enum">
+ <Entry value="0" text="All Points"/>
+ <Entry value="1" text="Every Nth Point"/>
+ <Entry value="2" text="Uniform Spatial Distribution"/>
+ </Domain>
+ </Property>
+ <Property name="GlyphTransform" id="4758.GlyphTransform" number_of_elements="1">
+ <Proxy value="4680"/>
+ <Domain name="proxy_list" id="4758.GlyphTransform.proxy_list">
+ <Proxy value="4680"/>
+ </Domain>
+ </Property>
+ <Property name="Input" id="4758.Input" number_of_elements="1">
+ <Domain name="groups" id="4758.Input.groups"/>
+ <Domain name="input_array1" id="4758.Input.input_array1"/>
+ <Domain name="input_array2" id="4758.Input.input_array2"/>
+ <Domain name="input_type" id="4758.Input.input_type"/>
+ </Property>
+ <Property name="MaximumNumberOfSamplePoints" id="4758.MaximumNumberOfSamplePoints" number_of_elements="1">
+ <Element index="0" value="5000"/>
+ <Domain name="range" id="4758.MaximumNumberOfSamplePoints.range"/>
+ </Property>
+ <Property name="Orient" id="4758.Orient" number_of_elements="1">
+ <Element index="0" value="1"/>
+ <Domain name="bool" id="4758.Orient.bool"/>
+ </Property>
+ <Property name="Scalars" id="4758.Scalars" number_of_elements="5">
+ <Element index="0" value=""/>
+ <Element index="1" value=""/>
+ <Element index="2" value=""/>
+ <Element index="3" value="0"/>
+ <Element index="4" value="Diameter"/>
+ <Domain name="array_list" id="4758.Scalars.array_list">
+ <String text="Diameter"/>
+ <String text="FixedVel"/>
+ <String text="Pressure [Pa]"/>
+ <String text="Shear Energy Rate [J/s]"/>
+ <String text="Shear Energy [J]"/>
+ <String text="Type color"/>
+ <String text="Viscous Energy Rate [J/s]"/>
+ <String text="Viscous Energy [J]"/>
+ </Domain>
+ </Property>
+ <Property name="ScaleFactor" id="4758.ScaleFactor" number_of_elements="1">
+ <Element index="0" value="1"/>
+ <Domain name="bounds" id="4758.ScaleFactor.bounds"/>
+ <Domain name="scalar_range" id="4758.ScaleFactor.scalar_range"/>
+ <Domain name="vector_range" id="4758.ScaleFactor.vector_range"/>
+ </Property>
+ <Property name="ScaleMode" id="4758.ScaleMode" number_of_elements="1">
+ <Element index="0" value="0"/>
+ <Domain name="enum" id="4758.ScaleMode.enum">
+ <Entry value="0" text="scalar"/>
+ <Entry value="1" text="vector"/>
+ <Entry value="2" text="vector_components"/>
+ <Entry value="3" text="off"/>
+ </Domain>
+ </Property>
+ <Property name="Seed" id="4758.Seed" number_of_elements="1">
+ <Element index="0" value="10339"/>
+ <Domain name="range" id="4758.Seed.range"/>
+ </Property>
+ <Property name="Source" id="4758.Source" number_of_elements="1">
+ <Proxy value="4736" output_port="0"/>
+ <Domain name="groups" id="4758.Source.groups"/>
+ <Domain name="input_type" id="4758.Source.input_type"/>
+ <Domain name="proxy_list" id="4758.Source.proxy_list">
+ <Proxy value="4681"/>
+ <Proxy value="4692"/>
+ <Proxy value="4703"/>
+ <Proxy value="4714"/>
+ <Proxy value="4725"/>
+ <Proxy value="4736"/>
+ <Proxy value="4747"/>
+ </Domain>
+ </Property>
+ <Property name="Stride" id="4758.Stride" number_of_elements="1">
+ <Element index="0" value="1"/>
+ <Domain name="range" id="4758.Stride.range"/>
+ </Property>
+ <Property name="Vectors" id="4758.Vectors" number_of_elements="5">
+ <Element index="0" value="1"/>
+ <Element index="1" value=""/>
+ <Element index="2" value=""/>
+ <Element index="3" value="0"/>
+ <Element index="4" value="Angular position[rad]"/>
+ <Domain name="array_list" id="4758.Vectors.array_list">
+ <String text="Angular position[rad]"/>
+ <String text="Angular velocity [rad/s]"/>
+ <String text="Displacement [m]"/>
+ <String text="Fluid pressure force [N]"/>
+ <String text="Force [N]"/>
+ <String text="Torque [Nm]"/>
+ <String text="Velocity [m/s]"/>
+ </Domain>
+ </Property>
+ </Proxy>
+ <Proxy group="extended_sources" type="Transform2" id="4680" servers="1" compound_name="auto_4680">
+ <Property name="Position" id="4680.Position" number_of_elements="3">
+ <Element index="0" value="0"/>
+ <Element index="1" value="0"/>
+ <Element index="2" value="0"/>
+ <Domain name="range" id="4680.Position.range"/>
+ </Property>
+ <Property name="PositionInfo" id="4680.PositionInfo" number_of_elements="3">
+ <Element index="0" value="0"/>
+ <Element index="1" value="0"/>
+ <Element index="2" value="0"/>
+ </Property>
+ <Property name="Rotation" id="4680.Rotation" number_of_elements="3">
+ <Element index="0" value="0"/>
+ <Element index="1" value="0"/>
+ <Element index="2" value="0"/>
+ <Domain name="range" id="4680.Rotation.range"/>
+ </Property>
+ <Property name="RotationInfo" id="4680.RotationInfo" number_of_elements="3">
+ <Element index="0" value="0"/>
+ <Element index="1" value="0"/>
+ <Element index="2" value="0"/>
+ </Property>
+ <Property name="Scale" id="4680.Scale" number_of_elements="3">
+ <Element index="0" value="1"/>
+ <Element index="1" value="1"/>
+ <Element index="2" value="1"/>
+ <Domain name="range" id="4680.Scale.range"/>
+ </Property>
+ <Property name="ScaleInfo" id="4680.ScaleInfo" number_of_elements="3">
+ <Element index="0" value="1"/>
+ <Element index="1" value="1"/>
+ <Element index="2" value="1"/>
+ </Property>
+ </Proxy>
+ <Proxy group="sources" type="SphereSource" id="4736" servers="1" compound_name="auto_4736">
+ <Property name="Center" id="4736.Center" number_of_elements="3">
+ <Element index="0" value="0"/>
+ <Element index="1" value="0"/>
+ <Element index="2" value="0"/>
+ <Domain name="range" id="4736.Center.range"/>
+ </Property>
+ <Property name="EndPhi" id="4736.EndPhi" number_of_elements="1">
+ <Element index="0" value="180"/>
+ <Domain name="range" id="4736.EndPhi.range"/>
+ </Property>
+ <Property name="EndTheta" id="4736.EndTheta" number_of_elements="1">
+ <Element index="0" value="360"/>
+ <Domain name="range" id="4736.EndTheta.range"/>
+ </Property>
+ <Property name="PhiResolution" id="4736.PhiResolution" number_of_elements="1">
+ <Element index="0" value="8"/>
+ <Domain name="range" id="4736.PhiResolution.range"/>
+ </Property>
+ <Property name="Radius" id="4736.Radius" number_of_elements="1">
+ <Element index="0" value="0.5"/>
+ <Domain name="range" id="4736.Radius.range"/>
+ </Property>
+ <Property name="StartPhi" id="4736.StartPhi" number_of_elements="1">
+ <Element index="0" value="0"/>
+ <Domain name="range" id="4736.StartPhi.range"/>
+ </Property>
+ <Property name="StartTheta" id="4736.StartTheta" number_of_elements="1">
+ <Element index="0" value="0"/>
+ <Domain name="range" id="4736.StartTheta.range"/>
+ </Property>
+ <Property name="ThetaResolution" id="4736.ThetaResolution" number_of_elements="1">
+ <Element index="0" value="8"/>
+ <Domain name="range" id="4736.ThetaResolution.range"/>
+ </Property>
+ </Proxy>
+ <ExposedProperties>
+ <Property name="Input" proxy_name="Glyph1" exposed_name="Input"/>
+ </ExposedProperties>
+ <OutputPort name="Output" proxy="Glyph1" port_index="0"/>
+ <Hints>
+ <ShowInMenu/>
+ </Hints>
+ </CompoundSourceProxy>
+ </CustomProxyDefinition>
+</CustomFilterDefinitions>
(DIR) diff --git a/output/cube-init.status.dat b/output/cube-init.status.dat
t@@ -1 +1 @@
-5.8000e+00 9.6667e+01 29
+3.8000e+00 6.3333e+01 19
(DIR) diff --git a/python/capillary-cohesion.py b/python/capillary-cohesion.py
t@@ -13,28 +13,37 @@
# top of a flat wall.
import sphere
-import numpy
+#import numpy
import sys
-device = sys.argv[1]
-cohesion = sys.argv[2]
-gravity = sys.argv[3]
+device = int(sys.argv[1])
+cohesion = int(sys.argv[2])
+gravity = int(sys.argv[3])
# Create packing
sim = sphere.sim('cap-cohesion=' + str(cohesion) + '-init-grav=' \
+ str(gravity), np=2000)
-sim.mu_s[0] = 0.0
-sim.mu_d[0] = 0.0
-sim.k_n[0] = 1.0e7
-sim.k_t[0] = 1.0e7
+#sim.mu_s[0] = 0.0
+#sim.mu_d[0] = 0.0
+#sim.k_n[0] = 1.0e7
+#sim.k_t[0] = 1.0e7
sim.generateRadii(psd='uni', mean=1.0e-3, variance=1.0e-4)
sim.contactModel(1)
-sim.initRandomGridPos([12, 12, 10000])
+sim.initRandomGridPos(gridnum=[24, 24, 10000], padding=1.4)
+sim.defaultParams(gamma_t = 1.0e3, capillaryCohesion=1)
sim.initTemporal(5.0, file_dt=0.01, epsilon=0.07)
+#I = numpy.nonzero(sim.x[:,2] < sim.L[2]*0.5)
+#sim.vel[I[0], 2] = 0.01 # add a instability seeding perturbation
+#I = numpy.nonzero(sim.x[:,2] > sim.L[2]*0.5)
+#sim.vel[I[0], 2] = -0.01 # add a instability seeding perturbation
if gravity == 1:
sim.g[2] = -10.0
-sim.run()
+sim.run(dry=True)
+sim.run(device=device)
+sim.writeVTKall()
+# if gravity is enabled, read last output file, place the sediment in a large
+# box, and restart time
if gravity == 1:
sim.readlast()
sim.sid = 'cap-cohesion=' + str(cohesion)
t@@ -50,6 +59,6 @@ if gravity == 1:
sim.x[:,1] += 0.5*sim.L[1] - 0.5*init_ly
sim.initTemporal(2.0, file_dt=0.01, epsilon=0.07)
- sim.run()
-
-sim.writeVTKall()
+ sim.run(dry=True)
+ sim.run(device=device)
+ sim.writeVTKall()
(DIR) diff --git a/python/channel.py b/python/channel.py
t@@ -2,8 +2,9 @@
import sphere
import numpy
-relaxation = False
-consolidation = True
+relaxation = True
+consolidation = False
+water = False
id_prefix = 'channel3'
N = 10e3
t@@ -88,7 +89,8 @@ sim.x[:, 2] = sim.x[:, 2] - min_z
sim.defineWorldBoundaries(L=[numpy.max(sim.x[:, 0] + sim.radius[:]),
numpy.max(sim.x[:, 1] + sim.radius[:]),
- numpy.max(sim.x[:, 2] + sim.radius[:])*1.2])
+ numpy.max(sim.x[:, 2] + sim.radius[:])*10.2])
+ #numpy.max(sim.x[:, 2] + sim.radius[:])*1.2])
sim.k_t[0] = 2.0/3.0*sim.k_n[0]
# sim.cleanup()
t@@ -107,10 +109,17 @@ sim.normalBoundariesXY()
# sim.consolidate(normal_stress=0.0)
# assign automatic colors, overwriting values from grid array
-sim.checkerboardColors(nx=grid.shape[1]/2, ny=2, nz=grid.shape[0])
+sim.checkerboardColors(nx=grid.shape[1], ny=2, nz=grid.shape[0]/4)
+
+sim.contactmodel[0] = 2
+sim.mu_s[0] = 0.5
+sim.mu_d[0] = 0.5
# Set duration of simulation, automatically determine timestep, etc.
sim.initTemporal(total=3.0, file_dt=0.01, epsilon=0.07)
+sim.time_dt[0] = 1.0e-20
+sim.time_file_dt = sim.time_dt
+sim.time_total = sim.time_file_dt*5.
sim.zeroKinematics()
if relaxation:
t@@ -118,38 +127,50 @@ if relaxation:
sim.run()
sim.writeVTKall()
+exit()
# Consolidation under constant normal stress
-sim.readlast()
-sim.id(id_prefix + '-' + str(int(N/1000.)) + 'kPa')
-sim.cleanup()
-sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
-
-# fix lowest plane of particles
-I = numpy.nonzero(sim.x[:, 2] < 1.5*numpy.mean(sim.radius))
-sim.fixvel[I] = -1
-sim.color[I] = 0
-
-sim.zeroKinematics()
-
-# Wall parameters
-sim.mu_ws[0] = 0.5
-sim.mu_wd[0] = 0.5
-sim.gamma_wn[0] = 1.0e2
-sim.gamma_wt[0] = 1.0e2
-# sim.gamma_wn[0] = 0.0
-# sim.gamma_wt[0] = 0.0
-
-# Particle parameters
-sim.mu_s[0] = 0.5
-sim.mu_d[0] = 0.5
-sim.gamma_n[0] = 0.0
-sim.gamma_t[0] = 0.0
-
-# apply effective normal stress from upper wall
-sim.consolidate(normal_stress=N)
-
if consolidation:
+ sim.readlast()
+ sim.id(id_prefix + '-' + str(int(N/1000.)) + 'kPa')
+ sim.cleanup()
+ sim.initTemporal(current=0.0, total=10.0, file_dt=0.01, epsilon=0.07)
+
+ # fix lowest plane of particles
+ I = numpy.nonzero(sim.x[:, 2] < 1.5*numpy.mean(sim.radius))
+ sim.fixvel[I] = -1
+ sim.color[I] = 0
+
+ sim.zeroKinematics()
+
+ # Wall parameters
+ sim.mu_ws[0] = 0.5
+ sim.mu_wd[0] = 0.5
+ sim.gamma_wn[0] = 1.0e2
+ sim.gamma_wt[0] = 1.0e2
+ # sim.gamma_wn[0] = 0.0
+ # sim.gamma_wt[0] = 0.0
+
+ # Particle parameters
+ sim.mu_s[0] = 0.5
+ sim.mu_d[0] = 0.5
+ sim.gamma_n[0] = 0.0
+ sim.gamma_t[0] = 0.0
+
+ # apply effective normal stress from upper wall
+ sim.consolidate(normal_stress=N)
+
sim.run(dry=True)
- sim.run()
- sim.writeVTKall()
+ #sim.run()
+ #sim.writeVTKall()
+
+## Add water
+#if water:
+ #sim.readlast()
+ #sim.id(id_prefix + '-wet')
+ #sim.wet()
+ #sim.initTemporal(total=3.0, file_dt=0.01, epsilon=0.07)
+#
+ #sim.run(dry=True)
+ #sim.run()
+ #sim.writeVTKall()
(DIR) diff --git a/python/cube-init.py b/python/cube-init.py
t@@ -1,7 +1,7 @@
#!/usr/bin/env python
import sphere
-init = sphere.sim('cube-init', np=1e2)
+init = sphere.sim('cube-init', np=1e3)
init.generateRadii(psd='uni', mean=0.01, variance=0.002)
t@@ -12,14 +12,14 @@ init.initRandomGridPos(gridnum=(6, 6, 1e12))
# Disable friction to dissipate energy fast
init.k_n[0] = 1.0e8
-init.mu_s[0] = 0.0
-init.mu_d[0] = 0.0
+init.mu_s[0] = 0.5
+init.mu_d[0] = 0.5
# Choose the tangential contact model
# 1) Visco-frictional (somewhat incorrect, fast computations)
# 2) Elastic-viscous-frictional (more correct, slow computations in dense
# packings)
-init.contactmodel[0] = 1
+init.contactmodel[0] = 2
# Add gravitational acceleration
init.g[2] = -10.0
t@@ -28,6 +28,6 @@ init.g[2] = -10.0
init.initTemporal(total=6.0, file_dt=0.2)
print(init.num)
-init.run(dry = True)
+init.run(dry=True)
init.run()
init.writeVTKall()
(DIR) diff --git a/python/rate-state.py b/python/rate-state.py
t@@ -0,0 +1,103 @@
+#!/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])
+device = 0
+wet = 0
+c_phi = 1.0
+k_c = 3.5e-13
+sigma0 = 80000.0
+mu = 2.080e-7
+velfac = 1.0
+
+if wet == 1:
+ fluid = True
+else:
+ fluid = False
+
+# load consolidated granular assemblage
+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()
+
+sim.shearVelocityGradient(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)
+
+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.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
+
+I = numpy.nonzero(sim.fixvel > 0)
+sim.fixvel[I] = 8.0 # step-wise velocity change when fixvel in ]5.0; 10.0[
+
+sim.run(dry=True)
+
+sim.run(device=device)
+sim.writeVTKall()
+sim.visualize('shear')
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -2645,7 +2645,9 @@ class sim:
self.num[1] += 1
- def initRandomGridPos(self, gridnum = numpy.array([12, 12, 32])):
+ def initRandomGridPos(self,
+ gridnum = numpy.array([12, 12, 32]),
+ padding=2.1):
'''
Initialize particle positions in loose, cubic configuration with some
variance. ``gridnum`` is the number of cells in the x, y and z
t@@ -2656,6 +2658,10 @@ class sim:
:param gridnum: The number of particles in x, y and z directions
:type gridnum: numpy.array
+ :param padding: Increase distance between particles in x, y and z
+ directions with this multiplier. Large values create more random
+ packings.
+ :type padding: float
'''
# Calculate cells in grid
t@@ -2665,7 +2671,7 @@ class sim:
r_max = numpy.amax(self.radius)
# Cells in grid 2*size to make space for random offset
- cellsize = 2.1 * r_max * 2
+ cellsize = padding * r_max * 2
# Check whether there are enough grid cells
if ((coarsegrid[0]-1)*(coarsegrid[1]-1)*(coarsegrid[2]-1)) < self.np:
(DIR) diff --git a/src/debug.h b/src/debug.h
t@@ -37,6 +37,9 @@ const int conv_log_interval = 10;
// Enable drag force and particle fluid coupling
#define CFD_DEM_COUPLING
+// Check if initial particle positions are finite values
+#define CHECK_PARTICLES_FINITE
+
// Check for nan/inf values in fluid solver kernels
#define CHECK_FLUID_FINITE
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -233,6 +233,44 @@ __global__ void checkConstantValues(int* dev_equal,
*dev_equal = 29; // Not ok
}
+__global__ void checkParticlePositions(
+ const Float4* __restrict__ dev_x)
+{
+ unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x; // Thread id
+
+ if (idx < devC_np) { // Condition prevents block size error
+ Float4 x = dev_x[idx];
+
+ // make sure grain doesn't have NaN or Inf position
+ if (!isfinite(x.x) || !isfinite(x.y) || !isfinite(x.z)) {
+ __syncthreads();
+ printf("\nParticle %d has non-finite position: x = %f %f %f",
+ idx, x.x, x.y, x.z);
+ }
+
+ /*__syncthreads();
+ printf("\nParticle %d: x = %f %f %f",
+ idx, x.x, x.y, x.z);*/
+
+ // check that the particle is inside of the simulation domain
+ if (x.x < devC_grid.origo[0] ||
+ x.y < devC_grid.origo[1] ||
+ x.z < devC_grid.origo[2] ||
+ x.x > devC_grid.L[0] ||
+ x.y > devC_grid.L[1] ||
+ x.z > devC_grid.L[2]) {
+ __syncthreads();
+ printf("\nParticle %d is outside the computational domain "
+ "(%f %f %f to %f %f %f): x = %f %f %f",
+ idx,
+ devC_grid.origo[0], devC_grid.origo[1], devC_grid.origo[2],
+ devC_grid.L[0], devC_grid.L[1], devC_grid.L[2],
+ x.x, x.y, x.z);
+ }
+ }
+}
+
+
// Copy the constant data components to device memory,
// and check whether the values correspond to the
// values in constant memory.
t@@ -904,6 +942,12 @@ __host__ void DEM::startTime()
double t_start = time.current;
double t_ratio; // ration between time flow in model vs. reality
+ // Hard-coded parameters for stepwise velocity change
+ int change_velocity_state = 0; // 1: increase velocity, 2: decrease vel.
+ const Float velocity_factor = 1.1; // v2 = v1*velocity_factor
+ const Float v2_start = 5.0; // seconds
+ const Float v2_end = 10.0; // seconds
+
// Index of dynamic top wall (if it exists)
unsigned int wall0_iz = 10000000;
// weight of fluid between two cells in z direction
t@@ -933,9 +977,17 @@ __host__ void DEM::startTime()
cout << " Current simulation time: " << time.current << " s.";
+
// MAIN CALCULATION TIME LOOP
while (time.current <= time.total) {
+ // check if particle positions have finite values
+#ifdef CHECK_PARTICLES_FINITE
+ checkParticlePositions<<<dimGrid, dimBlock>>>(dev_x);
+ cudaThreadSynchronize();
+ checkForCudaErrorsIter("Post checkParticlePositions", iter);
+#endif
+
// Print current step number to terminal
//printf("\n\n@@@ DEM time step: %ld\n", iter);
t@@ -2293,6 +2345,12 @@ __host__ void DEM::startTime()
checkForCudaErrorsIter("Post shear stress summation", iter);
}
+ // Determine whether it is time to step the velocity
+ if (time.current >= v2_start && time.current < v2_end)
+ change_velocity_state = 1.0;
+ else if (time.current >= 10.0)
+ change_velocity_state = -1.0;
+
// Update particle kinematics
if (PROFILING == 1)
startTimer(&kernel_tic);
t@@ -2317,6 +2375,8 @@ __host__ void DEM::startTime()
dev_walls_tau_eff_x_partial,
dev_walls_tau_x,
walls.tau_x[0],
+ change_velocity_state,
+ velocity_factor,
blocksPerGrid);
cudaThreadSynchronize();
checkForCudaErrorsIter("Post integrate", iter);
t@@ -2324,6 +2384,9 @@ __host__ void DEM::startTime()
stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
&t_integrate);
+ if (change_velocity_state != 0)
+ change_velocity_state = 0;
+
// Summation of forces on wall
if (PROFILING == 1)
startTimer(&kernel_tic);
(DIR) diff --git a/src/integration.cuh b/src/integration.cuh
t@@ -33,6 +33,8 @@ __global__ void integrate(
const Float* __restrict__ dev_walls_tau_eff_x_partial,
const Float* __restrict__ dev_walls_tau_x,
const Float tau_x,
+ const int change_velocity_state, // 1: v *= vel_fac, -1: v /= vel_fac
+ const Float velocity_factor,
const unsigned int blocksPerGrid)
{
unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x; // Thread id
t@@ -280,6 +282,14 @@ __global__ void integrate(
x_new.x -= L.x;
}
+ // step-wise velocity change for rate-and-state experiments
+ if (vel.w > 5.0 && vel.w < 10.0) {
+ if (change_velocity_state == 1)
+ vel_new.x *= velocity_factor;
+ else if (change_velocity_state == -1)
+ vel_new.x /= velocity_factor;
+ }
+
// Hold threads for coalesced write
__syncthreads();
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -387,6 +387,9 @@ void DEM::reportValues()
exit(1);
}
+ if (params.kappa > 0.0 && params.V_b > 0.0)
+ cout << " - Capillary cohesion enabled\n";
+
cout << " - Number of dynamic walls: " << walls.nw << '\n';
if (grid.periodic == 1)