tsign error in fixvel evaluation, still debugging - 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 132d83ec86f4d11cf364241e0825c68106259e27
 (DIR) parent 478b4e1cdaff74aa7a3269eb7dfc95341d10f368
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Wed,  2 Apr 2014 15:58:45 +0200
       
       sign error in fixvel evaluation, still debugging
       
       Diffstat:
         M python/fluidshear.py                |       5 +----
         M python/sphere.py                    |       1 -
         M src/integration.cuh                 |       2 +-
         M src/navierstokes.cuh                |      10 +++++-----
       
       4 files changed, 7 insertions(+), 11 deletions(-)
       ---
 (DIR) diff --git a/python/fluidshear.py b/python/fluidshear.py
       t@@ -7,7 +7,7 @@ sid = 'fluidshear'
        ## Initialization from loose packing to a gravitationally collapsed state
        ## without fluids
        sim = sphere.sim(sid + '-init', np = 2400, fluid = False)
       -sim.cleanup()
       +#sim.cleanup()
        sim.radius[:] = 0.05
        sim.initRandomGridPos(gridnum = [12, 12, 9000])
        sim.initTemporal(total = 5.0, file_dt = 0.05)
       t@@ -28,9 +28,6 @@ sim.run()
        sim.writeVTKall()
        sim.visualize('walls')
        
       -import sys
       -sys.exit()
       -
        ## Shear with fluids
        sim.readlast()
        sim.sid = sid + '-shear'
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -2086,7 +2086,6 @@ class sim:
                #self.w_m[idx] = numpy.array([self.rho[0]*self.np*math.pi \
                #        *(cellsize/2.0)**3])
                self.w_m[idx] = numpy.array([self.rho*self.L[0]*self.L[1]*d_max])
       -        print(self.w_m[idx])
        
            def consolidate(self, normal_stress = 10e3):
                '''
 (DIR) diff --git a/src/integration.cuh b/src/integration.cuh
       t@@ -97,7 +97,7 @@ __global__ void integrate(Float4* dev_x_sorted, Float4* dev_vel_sorted, // Input
                    angacc = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
                }
        
       -        if (vel.w < 0.0001) {
       +        if (vel.w < -0.0001) {
                    acc.x = 0.0;
                    acc.y = 0.0;
                    acc.z = 0.0;
 (DIR) diff --git a/src/navierstokes.cuh b/src/navierstokes.cuh
       t@@ -998,8 +998,8 @@ __global__ void findPorositiesVelocitiesDiametersSpherical(
        
                    // Save porosity, porosity change, average velocity and average diameter
                    __syncthreads();
       -            phi = 0.5; dphi = 0.0; // disable porosity effects
                    const unsigned int cellidx = idx(x,y,z);
       +            //phi = 0.5; dphi = 0.0; // disable porosity effects const unsigned int cellidx = idx(x,y,z);
                    dev_ns_phi[cellidx]  = phi;
                    dev_ns_dphi[cellidx] = dphi;
                    dev_ns_vp_avg[cellidx] = v_avg;
       t@@ -1015,8 +1015,6 @@ __global__ void findPorositiesVelocitiesDiametersSpherical(
        
                    __syncthreads();
                    const unsigned int cellidx = idx(x,y,z);
       -            //dev_ns_phi[cellidx]  = 1.0;
       -            //dev_ns_dphi[cellidx] = 0.0;
        
                    dev_ns_phi[cellidx]  = 1.0;
                    dev_ns_dphi[cellidx] = 0.0;
       t@@ -1682,7 +1680,8 @@ __global__ void findPredNSvelocities(
                //const Float3 f_g = devC_params.rho_f*dx*dy*dz*phi*g;
                const Float3 f_g
                    = MAKE_FLOAT3(devC_params.g[0], devC_params.g[1], devC_params.g[2])
       -            * devC_params.rho_f * dx*dy*dz * phi;
       +            * devC_params.rho_f * dx*dy*dz * 0.5;
       +            //* devC_params.rho_f * dx*dy*dz * phi;
                //const Float3 f_g = MAKE_FLOAT3(0.0, 0.0, 0.0);
        
                // Find pressure gradient
       t@@ -2329,7 +2328,8 @@ __global__ void findInteractionForce(
                if (v_rel_length > 0.0) {
                    if (phi <= 0.8)       // Ergun equation
                        fi = (150.0*devC_params.mu*not_phi*not_phi/(phi*d_avg*d_avg)
       -                        + 1.75*not_phi*devC_params.rho_f*v_rel_length/d_avg)*v_rel;
       +                        + 1.75*not_phi*devC_params.rho_f*v_rel_length/d_avg)
       +                    *v_rel;
                    else if (phi < 0.999) // Wen and Yu equation
                        fi = (3.0/4.0*cd*not_phi*pow(phi,
                                    -2.65)*devC_params.mu*devC_params.rho_f