tdevs should still be set, add function to set pressure at top wall - sphere - … | |
git clone git://src.adamsgaard.dk/sphere | |
Log | |
Files | |
Refs | |
LICENSE | |
--- | |
commit bebe373425e17094afa919530ad3c562435e06e0 | |
parent 0f3a2baafbb052d268a67b86531d6ae1d25d447d | |
Author: Anders Damsgaard <[email protected]> | |
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(-) | |
--- | |
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 | |
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) | |
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) { | |
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 | |
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 valu… | |
-// 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… | |
+// 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, |