tcfd_tests_neumann.py - 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
       ---
       tcfd_tests_neumann.py (2801B)
       ---
            1 #!/usr/bin/env python
            2 from pytestutils import *
            3 
            4 import sphere
            5 import sys
            6 import numpy
            7 import matplotlib.pyplot as plt
            8 
            9 print('### CFD tests - Dirichlet/Neumann BCs ###')
           10 
           11 print('''# Neumann bottom, Dirichlet top BC.
           12 # No gravity, no pressure gradients => no flow''')
           13 '''
           14 orig = sphere.sim("neumann", fluid = True)
           15 cleanup(orig)
           16 orig.defaultParams(mu_s = 0.4, mu_d = 0.4)
           17 orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
           18 #orig.initFluid(mu = 8.9e-4)
           19 orig.initFluid(mu = 0.0)
           20 orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
           21 py = sphere.sim(sid = orig.sid, fluid = True)
           22 orig.bc_bot[0] = 1      # No-flow BC at bottom (Neumann)
           23 #orig.run(dry=True)
           24 orig.run(verbose=False)
           25 #orig.writeVTKall()
           26 py.readlast(verbose = False)
           27 ones = numpy.ones((orig.num))
           28 py.readlast(verbose = False)
           29 compareNumpyArraysClose(ones, py.p_f, "Conservation of pressure:",
           30         tolerance = 1.0e-5)
           31 
           32 # Fluid flow along z should be very small
           33 if ((numpy.abs(py.v_f[:,:,:,:]) < 1.0e-6).all()):
           34     print("Flow field:\t\t" + passed())
           35 else:
           36     print("Flow field:\t\t" + failed())
           37     print(numpy.min(py.v_f))
           38     print(numpy.mean(py.v_f))
           39     print(numpy.max(py.v_f))
           40     raise Exception("Failed")
           41 '''
           42 
           43 print('''# Neumann bottom, Dirichlet top BC.
           44 # Gravity, pressure gradients => transient flow''')
           45 orig = sphere.sim("neumann", fluid = True)
           46 orig.cleanup()
           47 #orig.defineWorldBoundaries([0.4, 0.4, 1], dx = 0.1)
           48 orig.defineWorldBoundaries([0.3, 0.3, 0.3], dx = 0.1)
           49 #orig.g[2] = -10.0
           50 orig.initFluid(mu = 8.9e-4)
           51 orig.initTemporal(total = 0.05, file_dt = 0.005, dt = 1.0e-4)
           52 #orig.initTemporal(total = 1.0e-2, file_dt = 1.0e-4, dt = 1.0e-4)
           53 #orig.initTemporal(total = 1.0e-3, file_dt = 1.0e-4, dt = 1.0e-4)
           54 #print(orig.largestFluidTimeStep())
           55 #orig.initTemporal(total = orig.largestFluidTimeStep()*10.0,
           56         #file_dt = orig.largestFluidTimeStep(),
           57         #dt = orig.largestFluidTimeStep())
           58 py = sphere.sim(sid = orig.sid, fluid = True)
           59 orig.g[2] = -10.0
           60 orig.bc_bot[0] = 1      # No-flow BC at bottom (Neumann)
           61 #orig.run(dry=True)
           62 orig.run(verbose=False)
           63 #orig.writeVTKall()
           64 py.readlast(verbose = False)
           65 print(py.v_f)
           66 #ideal_grad_p_z = numpy.linspace(
           67 #        orig.p_f[0,0,0] + orig.L[2]*orig.rho_f*numpy.abs(orig.g[2]),
           68 #        orig.p_f[0,0,-1], orig.num[2])
           69 ideal_grad_p_z = numpy.linspace(
           70         orig.p_f[0,0,0] + (orig.L[2]-orig.L[2]/orig.num[2])*orig.rho_f*numpy.abs(orig.g[2]),
           71         orig.p_f[0,0,-1], orig.num[2])
           72 compareNumpyArraysClose(ideal_grad_p_z, py.p_f[0,0,:],
           73         "Pressure gradient:\t", tolerance=1.0e2)
           74 
           75 # Fluid flow along z should be very small
           76 #if ((numpy.abs(py.v_f[:,:,:,2]) < 5.0e-2).all()):
           77 if ((numpy.abs(py.v_f[:,:,:,2]) < 1.0e-4).all()):
           78     print("Flow field:\t\t" + passed())
           79 else:
           80     print("Flow field:\t\t" + failed())
           81     raise Exception("Failed")
           82 
           83 #orig.cleanup()