tdiffusivity.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
       ---
       tdiffusivity.py (2490B)
       ---
            1 #!/usr/bin/env python
            2 import sphere
            3 import numpy
            4 
            5 sid = 'diffusivity'
            6 
            7 ## Initialization from loose packing to a gravitationally collapsed state
            8 ## without fluids
            9 #sim = sphere.sim(sid + '-init', np = 2400, fluid = False)
           10 sim = sphere.sim(sid + '-init', np = 10000, fluid = False)
           11 #sim.cleanup()
           12 sim.radius[:] = 0.01
           13 #sim.initRandomGridPos(gridnum = [12, 12, 9000])
           14 sim.initRandomGridPos(gridnum = [24, 24, 9000])
           15 sim.initTemporal(total = 5.0, file_dt = 0.05, epsilon=0.07)
           16 sim.g[2] = -9.81
           17 sim.periodicBoundariesXY()
           18 #sim.run(dry=True)
           19 #sim.run()
           20 #sim.writeVTKall()
           21 
           22 # Stack two init assemblages on top of each other
           23 sim.readlast()
           24 max_z = numpy.max(sim.x[:,2] + sim.radius[:])
           25 for i in numpy.arange(sim.np):
           26     sim.addParticle(
           27             [sim.x[i,0], sim.x[i,1], sim.x[i,2] + max_z],
           28             radius=sim.radius[i])
           29 
           30 cellsize_min = 2.1*numpy.max(sim.radius)
           31 sim.L[2] = numpy.max(sim.x[:,2] + sim.radius[:])
           32 #sim.num[0] = numpy.ceil((sim.L[0]-sim.origo[0])/cellsize_min)
           33 #sim.num[1] = numpy.ceil((sim.L[1]-sim.origo[1])/cellsize_min)
           34 #sim.num[2] = numpy.ceil((sim.L[2]-sim.origo[2])/cellsize_min)
           35 sim.initGrid()
           36 
           37 ## Relaxation step
           38 sim.sid = sid + '-relax'
           39 sim.initTemporal(total = 2.0, file_dt = 0.05, epsilon=0.07)
           40 sim.mu_s[0] = 0.3
           41 sim.mu_d[0] = 0.3
           42 sim.zeroKinematics()
           43 sim.periodicBoundariesXY()
           44 #sim.run(dry=True)
           45 #sim.run()
           46 #sim.writeVTKall()
           47 
           48 ## Consolidation from a top wall with fluids
           49 sim.readlast()
           50 sim.sid = sid + '-cons-wet'
           51 
           52 # Checkerboard colors
           53 x_min = numpy.min(sim.x[:,0])
           54 x_max = numpy.max(sim.x[:,0])
           55 y_min = numpy.min(sim.x[:,1])
           56 y_max = numpy.max(sim.x[:,1])
           57 z_min = numpy.min(sim.x[:,2])
           58 z_max = numpy.max(sim.x[:,2])
           59 color_nx = 6
           60 color_ny = 6
           61 color_nz = 6
           62 for i in range(sim.np):
           63     ix = numpy.floor((sim.x[i,0] - x_min)/(x_max/color_nx))
           64     iy = numpy.floor((sim.x[i,1] - y_min)/(y_max/color_ny))
           65     iz = numpy.floor((sim.x[i,2] - z_min)/(z_max/color_nz))
           66     sim.color[i] = (-1)**ix + (-1)**iy + (-1)**iz
           67 
           68 sim.cleanup()
           69 sim.adjustUpperWall()
           70 sim.zeroKinematics()
           71 sim.consolidate(normal_stress = 10.0e3)
           72 sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True)  # mu = water at 0 deg C
           73 #sim.initFluid(mu = 8.9e-4, p = 1.0e5, hydrostatic = True)  # mu = water at 20 deg C
           74 sim.setFluidBottomNoFlow()
           75 sim.setFluidTopFixedPressure()
           76 #sim.setDEMstepsPerCFDstep(100)
           77 sim.setDEMstepsPerCFDstep(10)
           78 sim.initTemporal(total = 5.0, file_dt = 0.01, epsilon=0.07)
           79 sim.run(dry=True)
           80 sim.run()
           81 sim.writeVTKall()
           82 sim.visualize('walls')
           83 sim.visualize('fluid-pressure')