timplemented Neumann BCs (free slip, no slip) from Griebel et al 1998, UNTESTED - 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 903c2687b5df54b83f93bcceb6859b862633637e
(DIR) parent 02d6027fd3121cc0637cd39899d79af1e9320075
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Tue, 24 Jun 2014 14:35:23 +0200
implemented Neumann BCs (free slip, no slip) from Griebel et al 1998, UNTESTED
Diffstat:
M python/shortening.py | 15 ++++++++-------
M python/sphere.py | 26 +++++++++++++++++++++++---
M src/navierstokes.cuh | 55 +++++++++++++++++++-------------
3 files changed, 64 insertions(+), 32 deletions(-)
---
(DIR) diff --git a/python/shortening.py b/python/shortening.py
t@@ -137,16 +137,17 @@ for i in range(sim.np):
I = numpy.nonzero(sim.x[:,1] < 1.5*numpy.mean(sim.radius))
sim.fixvel[I] = -1
sim.color[I] = 0
+sim.x[I,1] = 0.0 # move into center into lower wall to avoid stuck particles
# fix left-most plane of particles
-I = numpy.nonzero(sim.x[:,2] < 1.5*numpy.mean(sim.radius))
-sim.fixvel[I] = -1
-sim.color[I] = 0
+#I = numpy.nonzero(sim.x[:,2] < 1.5*numpy.mean(sim.radius))
+#sim.fixvel[I] = -1
+#sim.color[I] = 0
# fix right-most plane of particles
-I = numpy.nonzero(sim.x[:,2] > z_max - 1.5*numpy.mean(sim.radius))
-sim.fixvel[I] = -1
-sim.color[I] = 0
+#I = numpy.nonzero(sim.x[:,2] > z_max - 1.5*numpy.mean(sim.radius))
+#sim.fixvel[I] = -1
+#sim.color[I] = 0
#sim.normalBoundariesXY()
#sim.periodicBoundariesX()
t@@ -170,7 +171,7 @@ sim.gamma_t[0] = 0.0
compressional_strain = 0.5
wall_velocity = -compressional_strain*(z_max - z_min)/sim.time_total[0]
sim.uniaxialStrainRate(wvel = wall_velocity)
-sim.vel[I,2] = wall_velocity
+#sim.vel[I,2] = wall_velocity
sim.run(dry=True)
sim.run()
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -300,7 +300,7 @@ 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
+ # 0: Dirichlet, 1: Neumann free slip, 2: Neumann no slip, 3: Periodic
self.bc_bot = numpy.zeros(1, dtype=numpy.int32)
self.bc_top = numpy.zeros(1, dtype=numpy.int32)
# Free slip boundaries? 1: yes
t@@ -2490,13 +2490,23 @@ class sim:
def setFluidBottomNoFlow(self):
'''
Set the lower boundary of the fluid domain to follow the no-flow
- (Neumann) boundary condition.
+ (Neumann) boundary condition with free slip parallel to the boundary.
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidBottomFixedPressure()`.
'''
self.bc_bot[0] = 1
+ def setFluidBottomNoFlowNoSlip(self):
+ '''
+ Set the lower boundary of the fluid domain to follow the no-flow
+ (Neumann) boundary condition with no slip parallel to the boundary.
+
+ The default behavior for the boundary is fixed value (Dirichlet), see
+ :func:`setFluidBottomFixedPressure()`.
+ '''
+ self.bc_bot[0] = 2
+
def setFluidBottomFixedPressure(self):
'''
Set the lower boundary of the fluid domain to follow the fixed pressure
t@@ -2510,13 +2520,23 @@ class sim:
def setFluidTopNoFlow(self):
'''
Set the upper boundary of the fluid domain to follow the no-flow
- (Neumann) boundary condition.
+ (Neumann) boundary condition with free slip parallel to the boundary.
The default behavior for the boundary is fixed value (Dirichlet), see
:func:`setFluidTopFixedPressure()`.
'''
self.bc_top[0] = 1
+ def setFluidTopNoFlowNoSlip(self):
+ '''
+ Set the upper boundary of the fluid domain to follow the no-flow
+ (Neumann) boundary condition with no slip parallel to the boundary.
+
+ The default behavior for the boundary is fixed value (Dirichlet), see
+ :func:`setFluidTopFixedPressure()`.
+ '''
+ self.bc_top[0] = 2
+
def setFluidTopFixedPressure(self):
'''
Set the upper boundary of the fluid domain to follow the fixed pressure
(DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -709,16 +709,27 @@ __global__ void setNSghostNodesFace(
dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Dirichlet -z
}
if (z == 0 && bc_bot == 1) {
- dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann -z
- dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Neumann -z
- dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Neumann -z
+ //dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann free slip -z
+ //dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Neumann free slip -z
+ //dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Neumann free slip -z
+ dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann free slip -z
+ dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Neumann free slip -z
+ dev_scalarfield_z[vidx(x,y,-1)] = 0.0; // Neumann free slip -z
}
if (z == 0 && bc_bot == 2) {
+ //dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Neumann no slip -z
+ //dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Neumann no slip -z
+ //dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Neumann no slip -z
+ dev_scalarfield_x[vidx(x,y,-1)] = -val_x; // Neumann no slip -z
+ dev_scalarfield_y[vidx(x,y,-1)] = -val_y; // Neumann no slip -z
+ dev_scalarfield_z[vidx(x,y,-1)] = 0.0; // Neumann no slip -z
+ }
+ if (z == 0 && bc_bot == 3) {
dev_scalarfield_x[vidx(x,y,nz)] = val_x; // Periodic -z
dev_scalarfield_y[vidx(x,y,nz)] = val_y; // Periodic -z
dev_scalarfield_z[vidx(x,y,nz)] = val_z; // Periodic -z
}
- if (z == 1 && bc_bot == 2) {
+ if (z == 1 && bc_bot == 3) {
dev_scalarfield_z[vidx(x,y,nz+1)] = val_z; // Periodic -z
}
t@@ -726,12 +737,26 @@ __global__ void setNSghostNodesFace(
dev_scalarfield_z[vidx(x,y,nz)] = val_z; // Dirichlet +z
}
if (z == nz-1 && bc_top == 1) {
- dev_scalarfield_x[vidx(x,y,nz)] = val_x; // Neumann +z
- dev_scalarfield_y[vidx(x,y,nz)] = val_y; // Neumann +z
- dev_scalarfield_z[vidx(x,y,nz)] = val_z; // Neumann +z
- dev_scalarfield_z[vidx(x,y,nz+1)] = val_z; // Neumann +z
+ //dev_scalarfield_x[vidx(x,y,nz)] = val_x; // Neumann free slip +z
+ //dev_scalarfield_y[vidx(x,y,nz)] = val_y; // Neumann free slip +z
+ //dev_scalarfield_z[vidx(x,y,nz)] = val_z; // Neumann free slip +z
+ //dev_scalarfield_z[vidx(x,y,nz+1)] = val_z; // Neumann free slip +z
+ dev_scalarfield_x[vidx(x,y,nz)] = val_x; // Neumann free slip +z
+ dev_scalarfield_y[vidx(x,y,nz)] = val_y; // Neumann free slip +z
+ dev_scalarfield_z[vidx(x,y,nz)] = 0.0; // Neumann free slip +z
+ dev_scalarfield_z[vidx(x,y,nz+1)] = 0.0; // Neumann free slip +z
}
if (z == nz-1 && bc_top == 2) {
+ //dev_scalarfield_x[vidx(x,y,nz)] = val_x; // Neumann no slip +z
+ //dev_scalarfield_y[vidx(x,y,nz)] = val_y; // Neumann no slip +z
+ //dev_scalarfield_z[vidx(x,y,nz)] = val_z; // Neumann no slip +z
+ //dev_scalarfield_z[vidx(x,y,nz+1)] = val_z; // Neumann no slip +z
+ dev_scalarfield_x[vidx(x,y,nz)] = -val_x; // Neumann no slip +z
+ dev_scalarfield_y[vidx(x,y,nz)] = -val_y; // Neumann no slip +z
+ dev_scalarfield_z[vidx(x,y,nz)] = 0.0; // Neumann no slip +z
+ dev_scalarfield_z[vidx(x,y,nz+1)] = 0.0; // Neumann no slip +z
+ }
+ if (z == nz-1 && bc_top == 3) {
dev_scalarfield_x[vidx(x,y,-1)] = val_x; // Periodic +z
dev_scalarfield_y[vidx(x,y,-1)] = val_y; // Periodic +z
dev_scalarfield_z[vidx(x,y,-1)] = val_z; // Periodic +z
t@@ -2285,20 +2310,6 @@ __global__ void findPredNSvelocities(
v_p.z = v.z;
//v_p.z = 0.0;
- // Do not set components placed in ghost node cells
- if (x == nx) {
- v_p.y = 0.0;
- v_p.z = 0.0;
- }
- if (y == ny) {
- v_p.x = 0.0;
- v_p.z = 0.0;
- }
- if (z == nz) {
- v_p.x = 0.0;
- v_p.y = 0.0;
- }
-
// Save the predicted velocity
__syncthreads();
dev_ns_v_p_x[fidx] = v_p.x;