tdo not calculate all velocity components at the outer grid edges - 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 88c6e28ce4e3c26bbde24ea7275f2e6430bd3572
 (DIR) parent 837f320e197769680e0c4ac3a95ede30f4bb5072
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Tue, 26 Aug 2014 15:32:17 +0200
       
       do not calculate all velocity components at the outer grid edges
       
       Diffstat:
         M src/debug.h                         |       4 ++--
         M src/navierstokes.cuh                |      29 ++++++++++++++++++++++++-----
       
       2 files changed, 26 insertions(+), 7 deletions(-)
       ---
 (DIR) diff --git a/src/debug.h b/src/debug.h
       t@@ -21,8 +21,8 @@ const unsigned int nijacnorm = 10;
        const int write_res_log = 0;
        
        // Report epsilon values during Jacobi iterations to stdout
       -//#define REPORT_EPSILON
       -//#define REPORT_MORE_EPSILON
       +#define REPORT_EPSILON
       +#define REPORT_MORE_EPSILON
        
        // Report the number of iterations it took before convergence to logfile
        // 'output/<sid>-conv.dat'
 (DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
       t@@ -2894,7 +2894,6 @@ __global__ void updateNSvelocity(
                const Float epsilon_yn = dev_ns_epsilon[idx(x,y-1,z)];
                const Float epsilon_zn = dev_ns_epsilon[idx(x,y,z-1)];
        
       -
                const Float phi_xn = dev_ns_phi[idx(x-1,y,z)];
                const Float phi_c  = dev_ns_phi[idx(x,y,z)];
                const Float phi_yn = dev_ns_phi[idx(x,y-1,z)];
       t@@ -2932,20 +2931,40 @@ __global__ void updateNSvelocity(
                   e_down);
                   }*/
        
       -        if ((z == 0 && bc_bot == 1) || (z == nz-1 && bc_top == 1))
       +        //if ((z == 0 && bc_bot == 1) || (z == nz-1 && bc_top == 1))
       +        if ((z == 0 && bc_bot == 1) || (z > nz-1 && bc_top == 1))
                    v.z = 0.0;
        
       -        if ((z == 0 && bc_bot == 2) || (z == nz-1 && bc_top == 2))
       +        //if ((z == 0 && bc_bot == 2) || (z == nz-1 && bc_top == 2))
       +        if ((z == 0 && bc_bot == 2) || (z > nz-1 && bc_top == 2))
                    v = MAKE_FLOAT3(0.0, 0.0, 0.0);
        
       +        // Do not calculate all components at the outer grid edges (nx, ny, nz)
       +        if (x == nx) {
       +            v.y = 0.0;
       +            v.z = 0.0;
       +        }
       +        if (y == ny) {
       +            v.x = 0.0;
       +            v.z = 0.0;
       +        }
       +        if (z == nz) {
       +            v.x = 0.0;
       +            v.y = 0.0;
       +        }
       +
                // Check the advection term using the Courant-Friedrichs-Lewy condition
                __syncthreads();
                if (v.x*ndem*devC_dt/dx
                    + v.y*ndem*devC_dt/dy
                    + v.z*ndem*devC_dt/dz > 1.0) {
                    printf("[%d,%d,%d] Warning: Advection term in fluid may be "
       -                   "unstable (CFL condition), v = %f,%f,%f\n",
       -                   x,y,z, v.x, v.y, v.z);
       +                   "unstable (CFL condition)\n"
       +                   "\tv = %.2e,%.2e,%.2e\n"
       +                   "\te_c,e_xn,e_yn,e_zn = %.2e,%.2e,%.2e,%.2e\n",
       +                   x,y,z, v.x, v.y, v.z,
       +                   epsilon_c, epsilon_xn, epsilon_yn, epsilon_zn
       +                   );
                }
        
                // Write new values