tconsolidation-curve.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
       ---
       tconsolidation-curve.py (3080B)
       ---
            1 #!/usr/bin/env python
            2 import matplotlib
            3 matplotlib.use('Agg')
            4 matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
            5 import os
            6 import shutil
            7 
            8 import sphere
            9 import numpy
           10 import matplotlib.pyplot as plt
           11 
           12 c_phi = 1.0
           13 #c_grad_p_list = [1.0, 0.1, 0.01, 0.001]
           14 #c_grad_p_list = [1.0, 0.1, 0.01]
           15 c_grad_p_list = [1.0, 0.1]
           16 #c_grad_p_list = [1.0]
           17 sigma0 = 10.0e3
           18 #sigma0 = 5.0e3
           19 
           20 t = [[], []]
           21 H = [[], []]
           22 #t = [[], [], []]
           23 #H = [[], [], []]
           24 #t = [[], [], [], []]
           25 #H = [[], [], [], []]
           26 
           27 c = 0
           28 for c_grad_p in c_grad_p_list:
           29 
           30     sid = 'cons-sigma0=' + str(sigma0) + '-c_phi=' + \
           31                      str(c_phi) + '-c_grad_p=' + str(c_grad_p)
           32     if c_grad_p != 1.0:
           33         sid += '-tall'
           34 
           35     if os.path.isfile('../output/' + sid + '.status.dat'):
           36         sim = sphere.sim(sid, fluid=True)
           37         t[c] = numpy.ones(sim.status()-1)
           38         H[c] = numpy.ones(sim.status()-1)
           39 
           40         #sim.visualize('walls')
           41         #sim.writeVTKall()
           42 
           43         #sim.plotLoadCurve()
           44         #sim.readfirst(verbose=True)
           45         #for i in numpy.arange(1, sim.status()+1):
           46         for i in numpy.arange(1, sim.status()):
           47             sim.readstep(i, verbose=False)
           48             t[c][i-1] = sim.time_current[0]
           49             H[c][i-1] = sim.w_x[0]
           50 
           51         '''
           52         # find consolidation parameters
           53         self.H0 = H[0]
           54         self.H100 = H[-1]
           55         self.H50 = (self.H0 + self.H100)/2.0
           56         T50 = 0.197 # case I
           57 
           58         # find the time where 50% of the consolidation (H50) has happened by
           59         # linear interpolation. The values in H are expected to be
           60         # monotonically decreasing. See Numerical Recipies p. 115
           61         i_lower = 0
           62         i_upper = self.status()-1
           63         while (i_upper - i_lower > 1):
           64             i_midpoint = int((i_upper + i_lower)/2)
           65             if (self.H50 < H[i_midpoint]):
           66                 i_lower = i_midpoint
           67             else:
           68                 i_upper = i_midpoint
           69         self.t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
           70                 (self.H50 - H[i_lower])/(H[i_upper] - H[i_lower])
           71 
           72         self.c_v = T50*self.H50**2.0/(self.t50)
           73         if self.fluid == True:
           74             e = numpy.mean(sb.phi[:,:,3:-8]) # ignore boundaries
           75         else:
           76             e = sb.voidRatio()
           77         '''
           78 
           79         H[c] -= H[c][0]
           80 
           81     c += 1
           82 
           83 # Normalize the thickness change
           84 #min_H = 0.0
           85 #for c in range(len(c_grad_p_list)):
           86     #min_H_c = numpy.min(H[c])
           87     #if min_H_c < min_H:
           88         #min_H = min_H_c
           89 
           90 plt.xlabel('Time [s]')
           91 #plt.ylabel('Normalized thickness change [-]')
           92 plt.ylabel('Thickness change [m]')
           93 #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
           94 #for c in range(len(c_grad_p_list)):
           95     #H[c] /= -min_H_c
           96 plt.semilogx(t[0], H[0], '-k', label='$c$ = %.2f' % (c_grad_p_list[0]))
           97 plt.semilogx(t[1], H[1], '--k', label='$c$ = %.2f' % (c_grad_p_list[1]))
           98 #plt.grid()
           99 
          100 plt.legend(loc='best', prop={'size':18}, fancybox=True, framealpha=0.5)
          101 plt.tight_layout()
          102 filename = 'cons-curves.pdf'
          103 plt.savefig(filename)
          104 #print(os.getcwd() + '/' + filename)
          105 shutil.copyfile(filename, '/home/adc/articles/own/2-org/' + filename)
          106 print(filename)