timprove diffusivity plot, time in continue sim additional length - 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 040315a879ff217c2fcde6a4a64e60a87ca3eea4
 (DIR) parent cb202e45887b3246e155760cef1ecfef296ab0c2
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Thu,  4 Sep 2014 13:42:38 +0200
       
       improve diffusivity plot, time in continue sim additional length
       
       Diffstat:
         M python/continue_sim.py              |       8 ++++----
         M python/diffusivity-results.py       |      25 +++++++++++++++----------
       
       2 files changed, 19 insertions(+), 14 deletions(-)
       ---
 (DIR) diff --git a/python/continue_sim.py b/python/continue_sim.py
       t@@ -4,11 +4,11 @@ import sys
        
        def print_usage():
            print('Usage: ' + sys.argv[0]
       -            + ' <simulation id> <fluid> <device> [end_time]')
       +            + ' <simulation id> <fluid> <device> [duration]')
            print('where "simulation id" is a string and "fluid" is either 0 or 1.')
            print('"device" is the number of the GPU device.')
       -    print('The total simulation can optionally be changed with the end_time '
       -            'parameter')
       +    print('The total simulation can optionally be defined to continue from the '
       +            'current time and "duration" seconds more.')
        
        if len(sys.argv) < 2:
            print_usage()
       t@@ -18,5 +18,5 @@ else:
            sim = sphere.sim(sys.argv[1], fluid = int(sys.argv[2]))
            sim.readlast()
            if len(sys.argv) == 5:
       -        sim.time_total = float(sys.argv[4])
       +        sim.time_total[0] = sim.time_current[0] + float(sys.argv[4])
            sim.run(device=sys.argv[3])
 (DIR) diff --git a/python/diffusivity-results.py b/python/diffusivity-results.py
       t@@ -13,8 +13,11 @@ c_phi = 1.0
        c_grad_p = 1.0
        #sigma0_list = numpy.array([5.0e3, 10.0e3, 20.0e3, 40.0e3, 80.0e3, 160.0e3])
        sigma0_list = numpy.array([10.0e3, 20.0e3, 40.0e3, 80.0e3, 160.0e3])
       -alpha = numpy.empty_like(sigma0_list)
       -phi_bar = numpy.empty_like(sigma0_list)
       +#alpha = numpy.empty_like(sigma0_list)
       +#phi_bar = numpy.empty_like(sigma0_list)
       +load = numpy.array([])
       +alpha = numpy.array([])
       +phi_bar = numpy.array([])
        
        #dc = diffusivitycalc.DiffusivityCalc()
        
       t@@ -26,10 +29,11 @@ for sigma0 in sigma0_list:
            if os.path.isfile('../output/' + sid + '.status.dat'):
                sim = sphere.sim(sid, fluid=True)
        
       -        sim.visualize('walls')
       +        #sim.visualize('walls')
                sim.plotLoadCurve()
       -        alpha[i] = sim.c_v
       -        phi_bar[i] = sim.phi_bar
       +        load = numpy.append(load, sigma0)
       +        alpha = numpy.append(alpha, sim.c_v)
       +        phi_bar = numpy.append(phi_bar, sim.phi_bar)
                #sim.writeVTKall()
        
            else:
       t@@ -38,19 +42,20 @@ for sigma0 in sigma0_list:
            i += 1
        
        fig, ax1 = plt.subplots()
       -sigma0_list /= 1000.0
       -
       +load /= 1000.0
        
        #plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
       -ax1.plot(sigma0_list, alpha, '.-k')
       +ax1.plot(load, alpha, 'o-k')
        ax1.set_xlabel('Normal stress $\\sigma_0$ [kPa]')
        ax1.set_ylabel('Hydraulic diffusivity $\\alpha$ [m$^2$s$^{-1}$]')
       +#ax1.ticklabel_format(style='plain', axis='y')
       +ax1.get_xaxis().get_major_formatter().set_useOffset(False)
        #ax1.grid()
        
        ax2 = ax1.twinx()
        color = 'b'
       -ax2.plot(sigma0_list, phi_bar, '.--' + color)
       -ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]')
       +ax2.plot(load, phi_bar, 'o--' + color)
       +ax2.set_ylabel('Mean porosity $\\bar{\\phi}$ [-]', color=color)
        for tl in ax2.get_yticklabels():
            tl.set_color(color)