tavoid nan/inf values when |v_rel| = 0 - 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 2a0eb86cc9df0ac9e97eef97b44d6c91f09334d2
 (DIR) parent ad6044c2f75bb6416249b65fc33496a4005c6fd8
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Fri, 23 May 2014 12:50:08 +0200
       
       avoid nan/inf values when |v_rel| = 0
       
       Diffstat:
         M src/navierstokes.cuh                |      11 ++++++++++-
       
       1 file changed, 10 insertions(+), 1 deletion(-)
       ---
 (DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
       t@@ -2827,9 +2827,12 @@ __global__ void findInteractionForce(
        
                const Float V_p = dx*dy*dz - phi*dx*dy*dz;
                const Float Re  = devC_params.rho_f*d*phi*v_rel_length/devC_params.mu;
       -        const Float Cd  = pow(0.63 + 4.8/pow(Re, 0.5), 2.0);
       +        Float Cd  = pow(0.63 + 4.8/pow(Re, 0.5), 2.0);
                const Float chi = 3.7 - 0.65*exp(-pow(1.5 - log10(Re), 2.0)/2.0);
        
       +        if (v_rel_length < 1.0e-6) // avoid Re=0 -> Cd=inf
       +            Cd = 0.0;
       +
                // Drag force
                const Float3 f_d = 0.125*Cd*devC_params.rho_f*M_PI*d*d*phi*phi
                    *v_rel_length*v_rel*pow(phi, -chi);
       t@@ -2852,6 +2855,12 @@ __global__ void findInteractionForce(
                const Float3 f_pf = f_d + f_p + f_v;
        
        #ifdef CHECK_NS_FINITE
       +        printf("%d [%d,%d,%d]\tV_p=%f Re=%f Cd=%f chi=%f\n",
       +                i, i_x, i_y, i_z, V_p, Re, Cd, chi);
       +        //checkFiniteFloat("V_p", i_x, i_y, i_z, V_p);
       +        //checkFiniteFloat("Re", i_x, i_y, i_z, Re);
       +        //checkFiniteFloat("Cd", i_x, i_y, i_z, Cd);
       +        //checkFiniteFloat("chi", i_x, i_y, i_z, chi);
                checkFiniteFloat3("f_d", i_x, i_y, i_z, f_d);
                checkFiniteFloat3("f_p", i_x, i_y, i_z, f_p);
                checkFiniteFloat3("f_v", i_x, i_y, i_z, f_v);