ttesting new porosity method, set fluid=true in initFluid - 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 7b3d83163c41681083609c0eda3027c040a27611
 (DIR) parent 782153575a3498d8d86e9e4c5deda9c5d573b848
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Fri, 25 Apr 2014 09:07:50 +0200
       
       ttesting new porosity method, set fluid=true in initFluid
       
       Diffstat:
         M python/sphere.py                    |       1 +
         M src/device.cu                       |       3 ++-
         M src/navierstokes.cuh                |      46 +++++++++++++++----------------
       
       3 files changed, 25 insertions(+), 25 deletions(-)
       ---
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -2328,6 +2328,7 @@ class sim:
                    has been specified
                :type hydrostatic: bool
                '''
       +        self.fluid = True
                self.mu = numpy.ones(1, dtype=numpy.float64) * mu
                self.rho_f = numpy.ones(1, dtype=numpy.float64) * rho
        
 (DIR) diff --git a/src/device.cu b/src/device.cu
       t@@ -844,7 +844,8 @@ __host__ void DEM::startTime()
                    // velocities
                    if (PROFILING == 1)
                        startTimer(&kernel_tic);
       -            findPorositiesVelocitiesDiametersSpherical
       +            //findPorositiesVelocitiesDiametersSpherical
       +            findPorositiesVelocitiesDiametersSphericalGradient
                        <<<dimGridFluid, dimBlockFluid>>>(
                                dev_cellStart,
                                dev_cellEnd,
 (DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
       t@@ -1123,8 +1123,12 @@ __global__ void findPorositiesVelocitiesDiametersSphericalGradient(
                    // Normal vector pointing from cell center towards particle center
                    Float3 n_p;
        
       +            // Normalized sphere-particle distance
       +            Float q;
       +
                    // Kernel function (2D disc)
       -            const Float dw_q = -1.0;
       +            //const Float dw_q = -1.0;
       +            Float dw_q;
        
                    // Iterate over 27 neighbor cells, R = 2*cell width
                    for (int z_dim=-2; z_dim<3; ++z_dim) { // z-axis
       t@@ -1169,33 +1173,27 @@ __global__ void findPorositiesVelocitiesDiametersSphericalGradient(
                                                    xr.y - X.y,
                                                    xr.z - X.z);
                                            d = length(x_p);
       -                                    n_p = x_p/length(x_p);
       +                                    n_p = x_p/d;
       +                                    q = d/R;
        
       -                                    // Lens shaped intersection
       -                                    if ((R - r) < d && d < (R + r)) {
       -                                        /*void_volume -=
       -                                            1.0/(12.0*d) * (
       -                                                    M_PI*(R + r - d)*(R + r - d)
       -                                                    *(d*d + 2.0*d*r - 3.0*r*r
       -                                                        + 2.0*d*R + 6.0*r*R
       -                                                        - 3.0*R*R) );*/
       -                                        v_avg += MAKE_FLOAT3(v.x, v.y, v.z);
       -                                        d_avg += 2.0*r;
       -                                        epsilon_ii +=
       -                                            dw_q*MAKE_FLOAT3(v.x, v.y, v.z)*n_p;
       -                                        n++;
       -                                    }
        
       -                                    // Particle fully contained in cell sphere
       -                                    if (d <= R - r) {
       -                                        //void_volume -= 4.0/3.0*M_PI*r*r*r;
       -                                        v_avg += MAKE_FLOAT3(v.x, v.y, v.z);
       -                                        d_avg += 2.0*r;
       -                                        epsilon_ii +=
       -                                            dw_q*MAKE_FLOAT3(v.x, v.y, v.z)*n_p;
       -                                        n++;
       +                                    dw_q = 0.0;
       +                                    if (0.0 < q && q < 1.0) {
       +                                        // 2d disc approximation
       +                                        //dw_q = -1.0;
       +
       +                                        // 3d sphere approximation
       +                                        dw_q = -1.5*pow(-q + 1.0, 0.5)
       +                                            *pow(q + 1.0, 0.5)
       +                                            + 0.5*pow(-q + 1.0, 1.5)
       +                                            *pow(q + 1.0, -0.5);
                                            }
        
       +                                    v_avg += MAKE_FLOAT3(v.x, v.y, v.z);
       +                                    d_avg += 2.0*r;
       +                                    epsilon_ii +=
       +                                        dw_q*MAKE_FLOAT3(v.x, v.y, v.z)*n_p;
       +                                    n++;
        
                                        }
                                    }