thansen-zoet-plots.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
       ---
       thansen-zoet-plots.py (1249B)
       ---
            1 #!/usr/bin/env python
            2 
            3 # Import sphere functionality
            4 import sphere
            5 import numpy
            6 import matplotlib.pyplot as plt
            7 from sklearn.linear_model import LinearRegression
            8 
            9 # Number of particles
           10 np = 1e4
           11 
           12 # Common simulation id
           13 sim_id = "hz"
           14 
           15 # Deviatoric stress [Pa]
           16 Nlist = numpy.array([51e3, 101e3, 202e3, 303e3, 404e3])
           17 
           18 ### INITIALIZATION ###
           19 
           20 # New class
           21 sim = sphere.sim(np = np, nd = 3, nw = 0, sid = sim_id + "-init")
           22 
           23 N_avg = 50
           24 tau_c = []
           25 
           26 for N in Nlist:
           27     sim.id("{}-shear-N{}".format(sim_id, N))
           28     lastf = sim.status()
           29     #sim.readlast()
           30     #sim.sheardisp()
           31 
           32     # average shear stress over last files
           33     tau_c_curr = 0.0
           34 
           35     for i in range(N_avg):
           36         sim.readstep(lastf - i)
           37         tau_c_curr += sim.shearStress('effective')
           38     tau_c.append(tau_c_curr/N_avg)
           39 
           40 fig = plt.figure()
           41 tau_c = numpy.array(tau_c)
           42 tau_c /= 1e3
           43 Nlist /= 1e3
           44 plt.plot(Nlist, tau_c, '+')
           45 mc = LinearRegression().fit(Nlist.reshape((-1, 1)), tau_c)
           46 Nspace = numpy.linspace(0.0, Nlist[-1], endpoint=True)
           47 plt.plot(Nspace, mc.predict(Nspace.reshape(-1,1)), '-')
           48 plt.title("friction = {:.2}, cohesion = {:.2} kPa".format(mc.coef_[0], mc.intercept_))
           49 plt.xlabel("Effective normal stress [kPa]")
           50 plt.ylabel("Shear stress [kPa]")
           51 plt.savefig("hz-mohr-coulomb.pdf")