tdiffusivity-test.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-test.py (2193B)
       ---
            1 #!/usr/bin/env python
            2 import sphere
            3 import numpy
            4 import sys
            5 
            6 # launch with:
            7 # $ python diffusivity-test <DEVICE> <C_PHI> <C_GRAD_P>
            8 
            9 # Unique simulation parameters
           10 device = int(sys.argv[1])
           11 c_phi = float(sys.argv[2])
           12 c_grad_p = float(sys.argv[3])
           13 
           14 # Unconsolidated input
           15 sim = sphere.sim('diffusivity-relax')
           16 sim.readlast()
           17 
           18 # Load sequence as per Bowles 1992, p. 135, units: Pa. dp/p = 1
           19 sigma0_list = numpy.array([5.0e3, 10.0e3, 20.0e3, 40.0e3, 80.0e3, 160.0e3])
           20 
           21 i = 0
           22 for sigma0 in sigma0_list:
           23 
           24     if (i == 0):
           25         i += 1
           26         continue
           27 
           28     # Read previous output if not first load test
           29     if (i > 0):
           30         sim.sid = 'cons-sigma0=' + str(sigma0_list[i-1]) + '-c_phi=' + \
           31                 str(c_phi) + '-c_grad_p=' + str(c_grad_p)
           32         sim.readlast()
           33 
           34     sim.sid = 'cons-sigma0=' + str(sigma0) + '-c_phi=' + str(c_phi) + \
           35             '-c_grad_p=' + str(c_grad_p)
           36     print('\n###### ' + sim.sid + ' ######')
           37 
           38     # Checkerboard colors
           39     x_min = numpy.min(sim.x[:,0])
           40     x_max = numpy.max(sim.x[:,0])
           41     y_min = numpy.min(sim.x[:,1])
           42     y_max = numpy.max(sim.x[:,1])
           43     z_min = numpy.min(sim.x[:,2])
           44     z_max = numpy.max(sim.x[:,2])
           45     color_nx = 6
           46     color_ny = 6
           47     color_nz = 6
           48     for n in range(sim.np):
           49         ix = numpy.floor((sim.x[n,0] - x_min)/(x_max/color_nx))
           50         iy = numpy.floor((sim.x[n,1] - y_min)/(y_max/color_ny))
           51         iz = numpy.floor((sim.x[n,2] - z_min)/(z_max/color_nz))
           52         sim.color[n] = (-1)**ix + (-1)**iy + (-1)**iz
           53 
           54     sim.cleanup()
           55     sim.adjustUpperWall()
           56     sim.zeroKinematics()
           57 
           58     sim.consolidate(normal_stress = sigma0)
           59 
           60     sim.initFluid(mu = 17.87e-4, p = 1.0e5, hydrostatic = True)
           61     sim.setFluidBottomNoFlow()
           62     sim.setFluidTopFixedPressure()
           63     sim.setDEMstepsPerCFDstep(10)
           64     sim.setMaxIterations(2e5)
           65     sim.initTemporal(total = 20.0, file_dt = 0.01, epsilon=0.07)
           66     sim.c_grad_p[0] = c_grad_p
           67     sim.c_phi[0] = c_phi
           68 
           69     # Fix lowermost particles
           70     dz = sim.L[2]/sim.num[2]
           71     I = numpy.nonzero(sim.x[:,2] < 1.5*dz)
           72     sim.fixvel[I] = 1
           73     
           74     sim.run(dry=True)
           75     sim.run(device=device)
           76     #sim.writeVTKall()
           77     sim.visualize('walls')
           78     sim.visualize('fluid-pressure')
           79 
           80     i += 1