tincrease velocity step factor to 10x - 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 8f10a0860a007df266703f23ba6a19a7846bc4ed
 (DIR) parent ad8cff4022033105dbd14ac7cc5a210420714f28
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Fri,  8 Apr 2016 12:38:22 -0700
       
       increase velocity step factor to 10x
       
       Diffstat:
         M output/cube-init.status.dat         |       2 +-
         M python/cube-init.py                 |       2 +-
         M python/rate-state.py                |     121 ++++++++++++++++++-------------
         M src/device.cu                       |       4 ++--
       
       4 files changed, 76 insertions(+), 53 deletions(-)
       ---
 (DIR) diff --git a/output/cube-init.status.dat b/output/cube-init.status.dat
       t@@ -1 +1 @@
       -5.0000e+00 8.3333e+01 25
       +1.8000e+00 3.0000e+01 9
 (DIR) diff --git a/python/cube-init.py b/python/cube-init.py
       t@@ -8,7 +8,7 @@ init.generateRadii(psd='uni', mean=0.01, variance=0.002)
        init.periodicBoundariesXY()
        
        # Initialize positions in random grid (also sets world size)
       -init.initRandomGridPos(gridnum=(6, 6, 1e12))
       +init.initRandomGridPos(gridnum=(7, 7, 1e12))
        
        # Disable friction to dissipate energy fast
        init.k_n[0] = 1.0e8
 (DIR) diff --git a/python/rate-state.py b/python/rate-state.py
       t@@ -9,7 +9,7 @@ import numpy
        # start with
        # ipython sigma-sim1-starter.py 0 1 1.0 2.0e-16 10000.0 2.080e-7 1.0
        
       -sid_prefix = 'ratestate2'
       +sid_prefix = 'ratestate3'
        
        # device = int(sys.argv[1])
        # wet = int(sys.argv[2])
       t@@ -22,10 +22,13 @@ device = 0
        wet = 0
        c_phi = 1.0
        k_c = 3.5e-13
       -sigma0 = 80000.0
       +#sigma0 = 80000.0
       +sigma0 = 40000.0
        mu = 2.080e-7
        velfac = 1.0
        
       +start_from_beginning = False
       +
        if wet == 1:
            fluid = True
        else:
       t@@ -33,11 +36,29 @@ else:
        
        # load consolidated granular assemblage
        #sim = sphere.sim('halfshear-sigma0=' + str(sigma0), fluid=False)
       -sim = sphere.sim('shear-sigma0=' + str(sigma0), fluid=False)
       -sim.readlast()
       +sim = sphere.sim(fluid=False)
       +if start_from_beginning:
       +    sim = sphere.sim('shear-sigma0=' + str(sigma0), fluid=False)
       +    sim.readlast()
       +else:
       +    if fluid:
       +        sim.id('ratestate2-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
       +                '-mu=' + str(mu) + '-shear')
       +    else:
       +        sim.id('ratestate2-sigma0=' + str(sigma0) + '-shear')
       +
       +    sim.readTime(4.9)
       +
       +    if fluid:
       +        sim.id(sid_prefix + '-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
       +                '-mu=' + str(mu) + '-shear')
       +    else:
       +        sim.id(sid_prefix + '-sigma0=' + str(sigma0) + '-shear')
       +
        #sim.readbin('../input/shear-sigma0=10000.0-new.bin')
        #sim.scaleSize(0.01) # from 1 cm to 0.01 cm = 100 micro m (fine sand)
        
       +
        sim.fluid = fluid
        if fluid:
            sim.id(sid_prefix + '-sigma0=' + str(sigma0) + '-k_c=' + str(k_c) + \
       t@@ -45,52 +66,54 @@ if fluid:
        else:
            sim.id(sid_prefix + '-sigma0=' + str(sigma0) + '-shear')
        
       -sim.checkerboardColors(nx=6,ny=3,nz=6)
       -#sim.checkerboardColors(nx=6,ny=6,nz=6)
       -sim.cleanup()
       -sim.adjustUpperWall()
       -sim.zeroKinematics()
       +if start_from_beginning:
       +    sim.checkerboardColors(nx=6,ny=3,nz=6)
       +    #sim.checkerboardColors(nx=6,ny=6,nz=6)
       +    sim.cleanup()
       +    sim.adjustUpperWall()
       +    sim.zeroKinematics()
       +
       +    sim.shear(1.0/20.0 * velfac)
       +    K_q_real = 36.4e9
       +    K_w_real =  2.2e9
       +
       +
       +    K_q_sim  = 1.16e9
       +    #K_q_sim = 1.16e6
       +    sim.setStiffnessNormal(K_q_sim)
       +    sim.setStiffnessTangential(K_q_sim)
       +    K_w_sim  = K_w_real/K_q_real * K_q_sim
       +
       +
       +    if fluid:
       +        #sim.num[2] *= 2
       +        sim.num[:] /= 2
       +        #sim.L[2] *= 2.0
       +        #sim.initFluid(mu = 1.787e-6, p = 600.0e3, cfd_solver = 1)
       +        sim.initFluid(mu = mu, p = 0.0, cfd_solver = 1)
       +        sim.setFluidTopFixedPressure()
       +        #sim.setFluidTopFixedFlow()
       +        sim.setFluidBottomNoFlow()
       +        #sim.setFluidBottomFixedPressure()
       +        #sim.setDEMstepsPerCFDstep(10)
       +        sim.setMaxIterations(2e5)
       +        sim.setPermeabilityPrefactor(k_c)
       +        sim.setFluidCompressibility(1.0/K_w_sim)
       +
       +    sim.w_sigma0[0] = sigma0
       +    sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
       +
       +    #sim.setStiffnessNormal(36.4e9 * 0.1 / 2.0)
       +    #sim.setStiffnessTangential(36.4e9/3.0 * 0.1 / 2.0)
       +    sim.setStiffnessNormal(K_q_sim)
       +    sim.setStiffnessTangential(K_q_sim)
       +    sim.mu_s[0] = 0.5
       +    sim.mu_d[0] = 0.5
       +    sim.setDampingNormal(0.0)
       +    sim.setDampingTangential(0.0)
       +
       +    sim.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
        
       -sim.shear(1.0/20.0 * velfac)
       -K_q_real = 36.4e9
       -K_w_real =  2.2e9
       -
       -
       -K_q_sim  = 1.16e9
       -#K_q_sim = 1.16e6
       -sim.setStiffnessNormal(K_q_sim)
       -sim.setStiffnessTangential(K_q_sim)
       -K_w_sim  = K_w_real/K_q_real * K_q_sim
       -
       -
       -if fluid:
       -    #sim.num[2] *= 2
       -    sim.num[:] /= 2
       -    #sim.L[2] *= 2.0
       -    #sim.initFluid(mu = 1.787e-6, p = 600.0e3, cfd_solver = 1)
       -    sim.initFluid(mu = mu, p = 0.0, cfd_solver = 1)
       -    sim.setFluidTopFixedPressure()
       -    #sim.setFluidTopFixedFlow()
       -    sim.setFluidBottomNoFlow()
       -    #sim.setFluidBottomFixedPressure()
       -    #sim.setDEMstepsPerCFDstep(10)
       -    sim.setMaxIterations(2e5)
       -    sim.setPermeabilityPrefactor(k_c)
       -    sim.setFluidCompressibility(1.0/K_w_sim)
       -
       -sim.w_sigma0[0] = sigma0
       -sim.w_m[0] = numpy.abs(sigma0*sim.L[0]*sim.L[1]/sim.g[2])
       -
       -#sim.setStiffnessNormal(36.4e9 * 0.1 / 2.0)
       -#sim.setStiffnessTangential(36.4e9/3.0 * 0.1 / 2.0)
       -sim.setStiffnessNormal(K_q_sim)
       -sim.setStiffnessTangential(K_q_sim)
       -sim.mu_s[0] = 0.5
       -sim.mu_d[0] = 0.5
       -sim.setDampingNormal(0.0)
       -sim.setDampingTangential(0.0)
       -
       -sim.initTemporal(total = 20.0/velfac, file_dt = 0.01/velfac, epsilon=0.07)
        
        I = numpy.nonzero(sim.fixvel > 0)
        sim.fixvel[I] = 8.0 # step-wise velocity change when fixvel in ]5.0; 10.0[
 (DIR) diff --git a/src/device.cu b/src/device.cu
       t@@ -942,10 +942,10 @@ __host__ void DEM::startTime()
            double t_start = time.current;
            double t_ratio;     // ration between time flow in model vs. reality
        
       -    // Hard-coded parameters for stepwise velocity change
       +    // Hard-coded parameters for stepwise velocity change (rate-state exp)
            int velocity_state = 1;  // 1: v1, 2: v2
            int change_velocity_state = 0;  // 1: increase velocity, 2: decrease vel.
       -    const Float velocity_factor = 1.1;  // v2 = v1*velocity_factor
       +    const Float velocity_factor = 10.0;  // v2 = v1*velocity_factor
            const Float v2_start = 5.0; // seconds
            const Float v2_end = 10.0;  // seconds