tplot results from variable c, dpdz and K together - 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 f505b33d5273d15b321358d1acf0fd3852842800
 (DIR) parent 0e8957887ab4ae4b2a08669f4a7e5898164ccd49
 (HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
       Date:   Thu, 21 Aug 2014 11:31:28 +0200
       
       plot results from variable c, dpdz and K together
       
       Diffstat:
         M python/permeability-results-c=1.py  |       8 ++++----
         M python/permeability-results.py      |      99 +++++++++++++++++++++-----------
       
       2 files changed, 70 insertions(+), 37 deletions(-)
       ---
 (DIR) diff --git a/python/permeability-results-c=1.py b/python/permeability-results-c=1.py
       t@@ -9,14 +9,14 @@ import sphere
        from permeabilitycalculator import *
        import matplotlib.pyplot as plt
        
       -sigma0_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
       +dp_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
        
        sids = []
       -for sigma0 in sigma0_list:
       -    sids.append('permeability-dp=' + str(sigma0))
       +for dp in dp_list:
       +    sids.append('permeability-dp=' + str(dp))
        
        K = numpy.empty(len(sids))
       -dpdz = numpy.empty_like(K)
       +dp = numpy.empty_like(K)
        Q = numpy.empty_like(K)
        phi_bar = numpy.empty_like(K)
        i = 0
 (DIR) diff --git a/python/permeability-results.py b/python/permeability-results.py
       t@@ -2,48 +2,81 @@
        import matplotlib
        matplotlib.use('Agg')
        matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
       -import os
        
       +import os
        import numpy
        import sphere
        from permeabilitycalculator import *
        import matplotlib.pyplot as plt
        
       -sids = [
       -    'permeability-dp=1000.0-c_phi=1.0-c_grad_p=0.01',
       -    'permeability-dp=1000.0-c_phi=1.0-c_grad_p=0.5',
       -    'permeability-dp=1000.0',
       -    'permeability-dp=20000.0-c_phi=1.0-c_grad_p=0.01',
       -    'permeability-dp=20000.0-c_phi=1.0-c_grad_p=0.1',
       -    'permeability-dp=20000.0-c_phi=1.0-c_grad_p=0.5',
       -    'permeability-dp=4000.0-c_phi=1.0-c_grad_p=0.01',
       -    'permeability-dp=4000.0-c_phi=1.0-c_grad_p=0.1',
       -    'permeability-dp=4000.0-c_phi=1.0-c_grad_p=0.5',
       -    'permeability-dp=4000.0']
       -
       -K = numpy.empty(len(sids))
       -c_grad_p = numpy.empty_like(K)
       -i = 0
       -
       -for sid in sids:
       -    pc = PermeabilityCalc(sid, plot_evolution=False, print_results=False,
       -            verbose=False)
       -    K[i] = pc.conductivity()
       -    c_grad_p[i] = pc.c_grad_p()
       -    i += 1
       -        
       -
       -# produce VTK files
       -#for sid in sids:
       -    #sim = sphere.sim(sid, fluid=True)
       -    #sim.writeVTKall()
       +dp_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
       +cvals = [1.0, 0.1, 0.01]
       +c_phi = 1.0
       +
       +K = [[], [], []]
       +dpdz = [[], [], []]
       +Q = [[], [], []]
       +phi_bar = [[], [], []]
       +
       +
       +c = 0
       +for c_grad_p in cvals:
       +
       +    sids = []
       +    for dp in dp_list:
       +        sids.append('permeability-dp=' + str(dp) + '-c_phi=' + \
       +                str(c_phi) + '-c_grad_p=' + str(c_grad_p))
       +
       +    K[c] = numpy.empty(len(sids))
       +    dpdz[c] = numpy.empty_like(K)
       +    Q[c] = numpy.empty_like(K)
       +    phi_bar[c] = numpy.empty_like(K)
       +    i = 0
       +
       +    for sid in sids:
       +        pc = PermeabilityCalc(sid, plot_evolution=False, print_results=False,
       +                verbose=False)
       +        K[c][i] = pc.conductivity()
       +        pc.findPressureGradient()
       +        pc.findCrossSectionalFlux()
       +        dpdz[c][i] = pc.dPdL[2]
       +        Q[c][i] = pc.Q[2]
       +        pc.findMeanPorosity()
       +        phi_bar[c][i] = pc.phi_bar
       +
       +        i += 1
       +
       +    # produce VTK files
       +    #for sid in sids:
       +        #sim = sphere.sim(sid, fluid=True)
       +        #sim.writeVTKall()
       +    c += 1
        
        fig = plt.figure()
       -plt.xlabel('Pressure gradient coefficient $c$ [-]')
       -plt.ylabel('Hydraulic conductivity $K$ [m/s]')
       -plt.plot(c_grad_p, K, '+')
       +
       +#plt.subplot(3,1,1)
       +plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
       +plt.ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
       +plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
       +for c in range(len(c_vals)):
       +    dpdz /= 1000.0
       +    plt.plot(dpdz[c], K[c], 'o-k', label='$c$ = %.2f' % (c_vals[c]))
        plt.grid()
       +
       +#plt.subplot(3,1,2)
       +#plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
       +#plt.ylabel('Hydraulic flux $Q$ [m$^3$s$^{-1}$]')
       +#plt.plot(dpdz, Q, '+')
       +#plt.grid()
       +
       +#plt.subplot(3,1,3)
       +#plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
       +#plt.ylabel('Mean porosity $\\bar{\\phi}$ [-]')
       +#plt.plot(dpdz, phi_bar, '+')
       +#plt.grid()
       +
        plt.tight_layout()
       -filename = 'c_grad_p-vs-K.pdf'
       +filename = 'permeability-dpdz-vs-K-vs-c.pdf'
        plt.savefig(filename)
        print(os.getcwd() + '/' + filename)
       +plt.savefig(filename)