tdevs should still be set, add function to set pressure at top wall - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
HTML git clone git://src.adamsgaard.dk/sphere
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
DIR commit bebe373425e17094afa919530ad3c562435e06e0
DIR parent 0f3a2baafbb052d268a67b86531d6ae1d25d447d
HTML Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 19 Sep 2014 11:16:39 +0200
devs should still be set, add function to set pressure at top wall
Diffstat:
M python/shear-starter.py | 6 +++---
M python/shear2.py | 2 +-
M src/device.cu | 19 +++++++++++++++++++
M src/integration.cuh | 6 ++++--
M src/navierstokes.cuh | 47 ++++++++++++++++++++++++++++---
5 files changed, 70 insertions(+), 10 deletions(-)
---
DIR diff --git a/python/shear-starter.py b/python/shear-starter.py
t@@ -44,7 +44,7 @@ sim.shear(1.0/20.0)
if fluid:
sim.num[2] *= 2
sim.L[2] *= 2.0
- sim.initFluid(mu = 1.787e-6, p = 1.0e5, hydrostatic = True)
+ sim.initFluid(mu = 1.787e-6, p = 600.0e3, hydrostatic = True)
#sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True)
sim.setFluidBottomNoFlow()
sim.setFluidTopFixedPressure()
t@@ -53,8 +53,8 @@ sim.setMaxIterations(2e5)
sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
sim.c_phi[0] = c_phi
sim.c_grad_p[0] = c_grad_p
-#sim.w_devs[0] = sigma0
-sim.w_devs[0] = 0.0
+sim.w_devs[0] = sigma0
+#sim.w_devs[0] = 0.0
sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
sim.mu_s[0] = 0.5
sim.mu_d[0] = 0.5
DIR diff --git a/python/shear2.py b/python/shear2.py
t@@ -11,7 +11,7 @@ sim.shear(1.0/20.0)
if fluid:
sim.num[2] *= 2
sim.L[2] *= 2.0
- sim.initFluid(mu=1.797e-6, p=600.0e3, hydrostatic=True)
+ sim.initFluid(mu=1.787e-6, p=600.0e3, hydrostatic=True)
sim.setFluidBottomNoFlow()
sim.setFluidTopFixedPressure()
sim.setDEMstepsPerCFDstep(100)
DIR diff --git a/src/device.cu b/src/device.cu
t@@ -840,6 +840,11 @@ __host__ void DEM::startTime()
double t_start = time.current;
double t_ratio; // ration between time flow in model vs. reality
+ // Index of dynamic top wall (if it exists)
+ unsigned int wall0_iz;
+ // weight of fluid between two cells in z direction
+ const Float dp_dz = fabs(params.rho_f*params.g[2]*grid.L[2]/grid.num[2]);
+
// Write a log file of the number of iterations it took before
// convergence in the fluid solver
std::ofstream convlog;
t@@ -1157,6 +1162,7 @@ __host__ void DEM::startTime()
if ((iter % ns.ndem) == 0) {
// Initial guess for the top epsilon values. These may be
// changed in setUpperPressureNS
+ // TODO: Check if this should only be set when top bc=Dirichlet
Float pressure = ns.p[idx(0,0,ns.nz-1)];
Float pressure_new = pressure; // Dirichlet
Float epsilon_value = pressure_new - ns.beta*pressure;
t@@ -1167,6 +1173,19 @@ __host__ void DEM::startTime()
cudaThreadSynchronize();
checkForCudaErrorsIter("Post setNSepsilonTop", iter);
+ // find cell containing top wall
+ if (walls.nw > 0 && walls.wmode[0] == 1) {
+ wall0_iz = walls.nx->w/(grid.L[2]/grid.num[2]);
+ setNSepsilonAtTopWall<<<dimGridFluid, dimBlockFluid>>>(
+ dev_ns_epsilon,
+ dev_ns_epsilon_new,
+ wall0_iz,
+ epsilon_value,
+ dp_dz);
+ cudaThreadSynchronize();
+ checkForCudaErrorsIter("Post setNSepsilonAtTopWall", iter);
+ }
+
// Modulate the pressures at the upper boundary cells
if ((ns.p_mod_A > 1.0e-5 || ns.p_mod_A < -1.0e-5) &&
ns.p_mod_f > 1.0e-7) {
DIR diff --git a/src/integration.cuh b/src/integration.cuh
t@@ -334,10 +334,12 @@ __global__ void integrateWalls(
const Float dt = devC_dt;
- // Normal load = Deviatoric stress times wall surface area,
- // directed downwards.
+ // Apply sinusoidal variation if specified
const Float sigma_0 = w_mvfd.w
+ devC_params.devs_A*sin(2.0*PI*devC_params.devs_f * t_current);
+
+ // Normal load = Normal stress times wall surface area,
+ // directed downwards.
const Float N = -sigma_0*devC_grid.L[0]*devC_grid.L[1];
// Calculate resulting acceleration of wall
DIR diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
t@@ -417,11 +417,9 @@ __global__ void setNSepsilonBottom(
}
}
-// Set the constant values of epsilon at the lower boundary. Since the
+// Set the constant values of epsilon at the upper boundary. Since the
// Dirichlet boundary values aren't transfered during array swapping, the values
-// also need to be written to the new array of epsilons. A value of 0 equals
-// the Dirichlet boundary condition: the new value should be identical to the
-// old value, i.e. the temporal gradient is 0
+// also need to be written to the new array of epsilons.
__global__ void setNSepsilonTop(
Float* __restrict__ dev_ns_epsilon,
Float* __restrict__ dev_ns_epsilon_new,
t@@ -442,6 +440,47 @@ __global__ void setNSepsilonTop(
dev_ns_epsilon_new[cellidx] = value;
}
}
+
+// Set the constant values of epsilon and grad_z(epsilon) at the upper wall, if
+// it is dynamic (Dirichlet+Neumann). Since the Dirichlet boundary values aren't
+// transfered during array swapping, the values also need to be written to the
+// new array of epsilons.
+__global__ void setNSepsilonAtTopWall(
+ Float* __restrict__ dev_ns_epsilon,
+ Float* __restrict__ dev_ns_epsilon_new,
+ const unsigned int iz,
+ const Float value,
+ const Float dp_dz)
+{
+ // 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;
+
+ const unsigned int cellidx = idx(x,y,z);
+
+ // cells containing the wall
+ if (x < devC_grid.num[0] && y < devC_grid.num[1] && z == iz) {
+ __syncthreads();
+ dev_ns_epsilon[cellidx] = value;
+ dev_ns_epsilon_new[cellidx] = value;
+ }
+
+ // cells above the wall
+ if (x < devC_grid.num[0] && y < devC_grid.num[1] && z == iz+1) {
+
+ // Pressure value in order to obtain hydrostatic pressure distribution
+ // for Neumann BC. The pressure should equal the value at the top wall
+ // minus the contribution due to the fluid weight.
+ // p_iz+1 = p_iz - rho_f*g*dz
+ const Float p = value - dp_dz;
+
+ __syncthreads();
+ dev_ns_epsilon[cellidx] = p;
+ dev_ns_epsilon_new[cellidx] = p;
+ }
+}
+
__device__ void copyNSvalsDev(
const unsigned int read,
const unsigned int write,