tstokes_law.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
       ---
       tstokes_law.py (2405B)
       ---
            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("### Stokes test - single sphere sedimentation ###")
           10 ## Stokes drag
           11 orig = sphere.sim(sid = "stokes_law_set2", fluid = True)
           12 cleanup(orig)
           13 orig.defaultParams()
           14 orig.addParticle([0.5,0.5,1.46], 0.05)
           15 #orig.defineWorldBoundaries([1.0,1.0,5.0], dx=0.1)
           16 #orig.defineWorldBoundaries([1.0,1.0,5.0])
           17 orig.defineWorldBoundaries([1.0,1.0,2.0])
           18 orig.initFluid(mu = 8.9e-4)
           19 #orig.initTemporal(total = 1.0, file_dt = 0.01)
           20 #orig.initTemporal(total = 1.0e-1, file_dt = 5.0e-3)
           21 #orig.initTemporal(total = 5.0, file_dt = 1.0e-2)
           22 orig.initTemporal(total = 1.0, file_dt = 1.0e-2)
           23 #orig.time_file_dt = orig.time_dt
           24 #orig.time_total = orig.time_dt*200
           25 #orig.time_file_dt = orig.time_dt*10
           26 #orig.time_total = orig.time_dt*1000
           27 #orig.g[2] = -10.0
           28 #orig.bc_bot[0] = 1      # No-flow BC at bottom
           29 #orig.setGamma(0.0)
           30 orig.vel[0,2] = -0.1
           31 #orig.vel[0,2] = -0.001
           32 #orig.setBeta(0.5)
           33 orig.setTolerance(1.0e-4)
           34 orig.setDEMstepsPerCFDstep(100)
           35 orig.run(dry=True)
           36 orig.run(verbose=True)
           37 py = sphere.sim(sid = orig.sid, fluid = True)
           38 
           39 ones = numpy.ones((orig.num))
           40 py.readlast(verbose = False)
           41 py.plotConvergence()
           42 py.writeVTKall()
           43 #compareNumpyArrays(ones, py.p_f, "Conservation of pressure:")
           44 
           45 it = numpy.loadtxt("../output/" + orig.sid + "-conv.log")
           46 #test((it[:,1] < 2000).all(), "Convergence rate:\t\t")
           47 
           48 t = numpy.empty(py.status())
           49 acc = numpy.empty(py.status())
           50 vel = numpy.empty(py.status())
           51 pos = numpy.empty(py.status())
           52 for i in range(py.status()):
           53     py.readstep(i+1, verbose=False)
           54     t[i] = py.time_current[0]
           55     acc[i] = py.force[0,2]/(V_sphere(py.radius[0])*py.rho[0]) + py.g[2]
           56     vel[i] = py.vel[0,2]
           57     pos[i] = py.x[0,2]
           58 
           59 fig = plt.figure()
           60 #plt.title('Convergence evolution in CFD solver in "' + self.sid + '"')
           61 plt.xlabel('Time [s]')
           62 plt.ylabel('$z$ value')
           63 plt.plot(t, acc, label='Acceleration')
           64 plt.plot(t, vel, label='Velocity')
           65 plt.plot(t, pos, label='Position')
           66 plt.grid()
           67 plt.legend()
           68 format = 'png'
           69 plt.savefig('./' + py.sid + '-stokes.' + format)
           70 plt.clf()
           71 plt.close(fig)
           72 #cleanup(orig)
           73 
           74 # Print 
           75 print('Final z-acceleration of particle: ' + str(acc[-1]) + ' m/s^2')
           76 print('Final z-velocity of particle:     ' + str(py.vel[0,2]) + ' m/s')
           77 print('Lowest fluid pressure:  ' + str(numpy.min(py.p_f)) + ' Pa')
           78 print('Highest fluid pressure: ' + str(numpy.max(py.p_f)) + ' Pa')