tplot results for all c values - 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 ae993c5f50daa01d742ab7681b7215da75581145
 (DIR) parent 75f00d3a5b69b8cb76c535aa57bcf66fbcbba7f9
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Wed, 27 Aug 2014 10:16:39 +0200
       
       plot results for all c values
       
       Diffstat:
         M python/consolidation-curve.py       |      82 ++++++++++++++++---------------
       
       1 file changed, 42 insertions(+), 40 deletions(-)
       ---
 (DIR) diff --git a/python/consolidation-curve.py b/python/consolidation-curve.py
       t@@ -9,8 +9,8 @@ import numpy
        import matplotlib.pyplot as plt
        
        c_phi = 1.0
       -#c_grad_p_list = [1.0, 0.1, 0.01]
       -c_grad_p_list = [1.0]
       +c_grad_p_list = [1.0, 0.1, 0.01]
       +#c_grad_p_list = [1.0]
        sigma0 = 10.0e3
        #sigma0 = 5.0e3
        
       t@@ -20,50 +20,52 @@ H = [[], [], []]
        c = 0
        for c_grad_p in c_grad_p_list:
        
       -    sim = sphere.sim('cons-sigma0=' + str(sigma0) + '-c_phi=' + \
       -                     str(c_phi) + '-c_grad_p=' + str(c_grad_p), fluid=True)
       -    t[c] = numpy.empty(sim.status())
       -    H[c] = numpy.empty(sim.status())
       +    sid = 'cons-sigma0=' + str(sigma0) + '-c_phi=' + \
       +                     str(c_phi) + '-c_grad_p=' + str(c_grad_p)
       +    if os.path.isfile('../output/' + sid + '.status.dat'):
       +        sim = sphere.sim(sid, fluid=True)
       +        t[c] = numpy.empty(sim.status())
       +        H[c] = numpy.empty(sim.status())
        
       -    #sim.visualize('walls')
       -    #sim.writeVTKall()
       +        #sim.visualize('walls')
       +        #sim.writeVTKall()
        
       -    #sim.plotLoadCurve()
       -    #sim.readfirst(verbose=True)
       -    for i in numpy.arange(1, sim.status()+1):
       -        sim.readstep(i, verbose=False)
       -        t[c][i-1] = sim.time_current[0]
       -        H[c][i-1] = sim.w_x[0]
       +        #sim.plotLoadCurve()
       +        #sim.readfirst(verbose=True)
       +        for i in numpy.arange(1, sim.status()+1):
       +            sim.readstep(i, verbose=False)
       +            t[c][i-1] = sim.time_current[0]
       +            H[c][i-1] = sim.w_x[0]
        
       -    '''
       -    # find consolidation parameters
       -    self.H0 = H[0]
       -    self.H100 = H[-1]
       -    self.H50 = (self.H0 + self.H100)/2.0
       -    T50 = 0.197 # case I
       +        '''
       +        # find consolidation parameters
       +        self.H0 = H[0]
       +        self.H100 = H[-1]
       +        self.H50 = (self.H0 + self.H100)/2.0
       +        T50 = 0.197 # case I
        
       -    # find the time where 50% of the consolidation (H50) has happened by
       -    # linear interpolation. The values in H are expected to be
       -    # monotonically decreasing. See Numerical Recipies p. 115
       -    i_lower = 0
       -    i_upper = self.status()-1
       -    while (i_upper - i_lower > 1):
       -        i_midpoint = int((i_upper + i_lower)/2)
       -        if (self.H50 < H[i_midpoint]):
       -            i_lower = i_midpoint
       -        else:
       -            i_upper = i_midpoint
       -    self.t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
       -            (self.H50 - H[i_lower])/(H[i_upper] - H[i_lower])
       +        # find the time where 50% of the consolidation (H50) has happened by
       +        # linear interpolation. The values in H are expected to be
       +        # monotonically decreasing. See Numerical Recipies p. 115
       +        i_lower = 0
       +        i_upper = self.status()-1
       +        while (i_upper - i_lower > 1):
       +            i_midpoint = int((i_upper + i_lower)/2)
       +            if (self.H50 < H[i_midpoint]):
       +                i_lower = i_midpoint
       +            else:
       +                i_upper = i_midpoint
       +        self.t50 = t[i_lower] + (t[i_upper] - t[i_lower]) * \
       +                (self.H50 - H[i_lower])/(H[i_upper] - H[i_lower])
        
       -    self.c_v = T50*self.H50**2.0/(self.t50)
       -    if self.fluid == True:
       -        e = numpy.mean(sb.phi[:,:,3:-8]) # ignore boundaries
       -    else:
       -        e = sb.voidRatio()
       -    '''
       +        self.c_v = T50*self.H50**2.0/(self.t50)
       +        if self.fluid == True:
       +            e = numpy.mean(sb.phi[:,:,3:-8]) # ignore boundaries
       +        else:
       +            e = sb.voidRatio()
       +        '''
        
       -    H[c] -= H[c][0]
       +        H[c] -= H[c][0]
        
            c += 1