tcorrected BCs at top wall - 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 5e234fe2c18c005b32bd0747fc74c2485ec09bb8
 (DIR) parent ef4d24ef5091809277658fe911f92035efd89e4e
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Fri, 19 Sep 2014 15:32:26 +0200
       
       corrected BCs at top wall
       
       Diffstat:
         M src/device.cu                       |       7 ++++---
         M src/navierstokes.cuh                |      20 +++++++++++++++-----
       
       2 files changed, 19 insertions(+), 8 deletions(-)
       ---
 (DIR) diff --git a/src/device.cu b/src/device.cu
       t@@ -1500,7 +1500,8 @@ __host__ void DEM::startTime()
                                    ns.bc_bot,
                                    ns.bc_top,
                                    ns.theta,
       -                            wall0_iz);
       +                            wall0_iz,
       +                            dp_dz);
                            cudaThreadSynchronize();
                            if (PROFILING == 1)
                                stopTimer(&kernel_tic, &kernel_toc, &kernel_elapsed,
       t@@ -1508,7 +1509,7 @@ __host__ void DEM::startTime()
                            checkForCudaErrorsIter("Post jacobiIterationNS", iter);
        
                            // set Dirichlet and Neumann BC at cells containing top wall
       -                    if (walls.nw > 0 && walls.wmode[0] == 1) {
       +                    /*if (walls.nw > 0 && walls.wmode[0] == 1) {
                                setNSepsilonAtTopWall<<<dimGridFluid, dimBlockFluid>>>(
                                        dev_ns_epsilon,
                                        dev_ns_epsilon_new,
       t@@ -1518,7 +1519,7 @@ __host__ void DEM::startTime()
                                cudaThreadSynchronize();
                                checkForCudaErrorsIter("Post setNSepsilonAtTopWall",
                                        iter);
       -                    }
       +                    }*/
        
                            // Copy new values to current values
                            copyValues<Float><<<dimGridFluid, dimBlockFluid>>>(
 (DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
       t@@ -459,15 +459,17 @@ __global__ void setNSepsilonAtTopWall(
        
            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) {
       +    // cells containing the wall (Dirichlet BC)
       +    if (x < devC_grid.num[0] && y < devC_grid.num[1] && z < devC_grid.num[2] &&
       +            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) {
       +    // cells above the wall (Neumann BC)
       +    if (x < devC_grid.num[0] && y < devC_grid.num[1] && z < devC_grid.num[2] &&
       +            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
       t@@ -2670,7 +2672,8 @@ __global__ void jacobiIterationNS(
            const int bc_bot,
            const int bc_top,
            const Float theta,
       -    const unsigned int wall0_iz)
       +    const unsigned int wall0_iz,
       +    const Float dp_dz)
        {
            // 3D thread index
            const unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
       t@@ -2739,9 +2742,16 @@ __global__ void jacobiIterationNS(
                       + dxdx*dydy*(e_zn + e_zp))
                    /(2.0*(dxdx*dydy + dxdx*dzdz + dydy*dzdz));
        
       +
       +        // Dirichlet BC at dynamic top wall. wall0_iz will be larger than the
       +        // grid if the wall isn't dynamic
                if (z == wall0_iz)
                    e_new = e;
        
       +        // Neumann BC at dynamic top wall
       +        if (z == wall0_iz + 1)
       +            e_new = e_zn + dp_dz;
       +
                // New value of epsilon in 1D update
                //const Float e_new = (e_zp + e_zn - dz*dz*f)/2.0;