tcorrected grid boundary check, debugging cfd - 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 90c2fb7d69cb5dfcfdf6f426af69b9ecf80c38de
 (DIR) parent 86fdfad131a636d5ce80131892b66d157250f0c5
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Wed, 18 Jun 2014 08:58:08 +0200
       
       corrected grid boundary check, debugging cfd
       
       Diffstat:
         M python/sphere.py                    |       3 ++-
         M tests/cfd_tests.py                  |      10 ++++++++++
         M tests/cfd_tests_neumann.py          |      12 ++++++++----
       
       3 files changed, 20 insertions(+), 5 deletions(-)
       ---
 (DIR) diff --git a/python/sphere.py b/python/sphere.py
       t@@ -1783,7 +1783,8 @@ class sim:
                self.num[2] = numpy.ceil((self.L[2]-self.origo[2])/cellsize_min)
        
                #if (self.num.any() < 4):
       -        if (self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4):
       +        #if (self.num[0] < 4 or self.num[1] < 4 or self.num[2] < 4):
       +        if (self.num[0] < 3 or self.num[1] < 3 or self.num[2] < 3):
                    raise Exception("Error: The grid must be at least 3 cells in each "
                    + "direction\nGrid: x={}, y={}, z={}\n".format(\
                            self.num[0], self.num[1], self.num[2])
 (DIR) diff --git a/tests/cfd_tests.py b/tests/cfd_tests.py
       t@@ -34,6 +34,16 @@ compareNumpyArrays(ones, py.p_f, "Conservation of pressure:")
        it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
        compare(it[:,1].sum(), 0.0, "Convergence rate (1/2):\t")
        
       +# Fluid flow should be very small
       +if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
       +    print("Flow field:\t\t" + passed())
       +else:
       +    print("Flow field:\t\t" + failed())
       +    print(numpy.min(py.v_f))
       +    print(numpy.mean(py.v_f))
       +    print(numpy.max(py.v_f))
       +    raise Exception("Failed")
       +
        # Add pressure gradient
        # This test passes with BETA=0.0 and tolerance=1.0e-9
        orig.p_f[:,:,-1] = 1.1
 (DIR) diff --git a/tests/cfd_tests_neumann.py b/tests/cfd_tests_neumann.py
       t@@ -14,7 +14,8 @@ orig = sphere.sim("neumann", fluid = True)
        cleanup(orig)
        orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
        orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
       -orig.initFluid(mu = 8.9e-4)
       +#orig.initFluid(mu = 8.9e-4)
       +orig.initFluid(mu = 0.0)
        orig.initTemporal(total = 0.5, file_dt = 0.05, dt = 1.0e-4)
        py = sphere.sim(sid = orig.sid, fluid = True)
        orig.bc_bot[0] = 1      # No-flow BC at bottom (Neumann)
       t@@ -25,13 +26,16 @@ py.readlast(verbose = False)
        ones = numpy.ones((orig.num))
        py.readlast(verbose = False)
        compareNumpyArraysClose(ones, py.p_f, "Conservation of pressure:",
       -        tolerance = 1.0e-1)
       +        tolerance = 1.0e-3)
        
        # Fluid flow along z should be very small
       -if ((numpy.abs(py.v_f[:,:,:,2]) < 1.0e-4).all()):
       +if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
            print("Flow field:\t\t" + passed())
        else:
            print("Flow field:\t\t" + failed())
       +    print(numpy.min(py.v_f))
       +    print(numpy.mean(py.v_f))
       +    print(numpy.max(py.v_f))
            raise Exception("Failed")
        
        print('''# Neumann bottom, Dirichlet top BC.
       t@@ -52,7 +56,7 @@ ideal_grad_p_z = numpy.linspace(
                orig.p_f[0,0,0] + orig.L[2]*orig.rho_f*numpy.abs(orig.g[2]),
                orig.p_f[0,0,-1], orig.num[2])
        compareNumpyArraysClose(ideal_grad_p_z, py.p_f[0,0,:],
       -        "Pressure gradient:\t", tolerance=1.0e3)
       +        "Pressure gradient:\t", tolerance=1.0e2)
        
        # Fluid flow along z should be very small
        if ((numpy.abs(py.v_f[:,:,:,2]) < 5.0e-2).all()):