tincorporate shear stress BC in integrate function - 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 cbd624d82308f27331fa8781ba282eff5661645f
   DIR parent d836c03f0c503c6ecde5e9d3cf3f2dcd65c89532
  HTML Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Wed, 14 Jan 2015 13:39:42 +0100
       
       incorporate shear stress BC in integrate function
       
       Diffstat:
         M src/integration.cuh                 |      33 +++++++++++++++++++++++++++++--
       
       1 file changed, 31 insertions(+), 2 deletions(-)
       ---
   DIR diff --git a/src/integration.cuh b/src/integration.cuh
       t@@ -27,7 +27,11 @@ __global__ void integrate(
            Float4* __restrict__ dev_angvel0,
            Float4* __restrict__ dev_xyzsum,
            const unsigned int* __restrict__ dev_gridParticleIndex,
       -    const unsigned int iter)
       +    const unsigned int iter,
       +    const int* __restrict__ dev_walls_wmode,
       +    const Float* __restrict__ dev_walls_tau_eff_x_partial,
       +    const Float* __restrict__ dev_walls_tau_x,
       +    const unsigned int blocksPerGrid)
        {
            unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x; // Thread id
        
       t@@ -45,6 +49,16 @@ __global__ void integrate(
                const Float4 angvel   = dev_angvel_sorted[idx];
                Float4 xyzsum = dev_xyzsum[orig_idx];
        
       +        // Read the mode of the top wall
       +        const int wall0mode = dev_walls_wmode[0];
       +        const Float wall0mass = dev_walls_mvfd[0].x;
       +
       +        // Find the final sum of shear stresses on the top particles
       +        Float tau_eff_x = 0.0;
       +        if (wall0mode == 3)
       +            for (int i=0; i<blocksPerGrid; ++i)
       +                tau_eff_x += dev_walls_tau_eff_x_partial[i];
       +
                // Get old accelerations for three-term Taylor expansion. These values
                // don't exist in the first time step
        #ifdef TY3
       t@@ -91,11 +105,26 @@ __global__ void integrate(
                angacc.y = torque.y * 1.0 / (2.0/5.0 * m * radius*radius);
                angacc.z = torque.z * 1.0 / (2.0/5.0 * m * radius*radius);
        
       +        // Fixed shear stress BC
       +        if (wall0mode == 3 && vel.w > 0.0001 && x.z > devC_grid.L[2]*0.5) {
       +
       +            // the force should be positive when abs(tau) > abs(tau_eff_x)
       +            const Float f_tau_x =
       +                (tau + tau_eff_x)*devC_grid.L[0]*devC_grid.L[1];
       +
       +            acc.x = f_tau_x/wall0mass;  // acceleration = force/mass
       +            acc.y = 0.0;
       +            acc.z -= devC_params.g[2];
       +
       +            // disable rotation
       +            angacc = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
       +        }
       +
                // Modify the acceleration if the particle is marked as having a fixed
                // velocity. In that case, zero the horizontal acceleration and disable
                // gravity to counteract segregation. Particles may move in the
                // z-dimension, to allow for dilation.
       -        if (vel.w > 0.0001 && vel.w < 10.0) {
       +        else if (vel.w > 0.0001 && vel.w < 10.0) {
        
                    acc.x = 0.0;
                    acc.y = 0.0;