tpermeability-results-c=1.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
       ---
       tpermeability-results-c=1.py (1690B)
       ---
            1 #!/usr/bin/env python
            2 import matplotlib
            3 matplotlib.use('Agg')
            4 matplotlib.rcParams.update({'font.size': 18, 'font.family': 'serif'})
            5 
            6 import os
            7 import numpy
            8 import sphere
            9 from permeabilitycalculator import *
           10 import matplotlib.pyplot as plt
           11 
           12 dp_list = numpy.array([1.0e3, 2.0e3, 4.0e3, 10.0e3, 20.0e3, 40.0e3])
           13 
           14 sids = []
           15 for dp in dp_list:
           16     sids.append('permeability-dp=' + str(dp))
           17 
           18 K = numpy.empty(len(sids))
           19 dp = numpy.empty_like(K)
           20 Q = numpy.empty_like(K)
           21 phi_bar = numpy.empty_like(K)
           22 i = 0
           23 
           24 for sid in sids:
           25     pc = PermeabilityCalc(sid, plot_evolution=False, print_results=False,
           26             verbose=False)
           27     K[i] = pc.conductivity()
           28     pc.findPressureGradient()
           29     pc.findCrossSectionalFlux()
           30     dpdz[i] = pc.dPdL[2]
           31     Q[i] = pc.Q[2]
           32     pc.findMeanPorosity()
           33     phi_bar[i] = pc.phi_bar
           34 
           35     i += 1
           36 
           37 # produce VTK files
           38 #for sid in sids:
           39     #sim = sphere.sim(sid, fluid=True)
           40     #sim.writeVTKall()
           41 
           42 fig = plt.figure()
           43 
           44 #plt.subplot(3,1,1)
           45 plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [kPa m$^{-1}$]')
           46 plt.ylabel('Hydraulic conductivity $K$ [ms$^{-1}$]')
           47 plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
           48 dpdz /= 1000.0
           49 plt.plot(dpdz, K, 'o-k')
           50 plt.grid()
           51 
           52 #plt.subplot(3,1,2)
           53 #plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
           54 #plt.ylabel('Hydraulic flux $Q$ [m$^3$s$^{-1}$]')
           55 #plt.plot(dpdz, Q, '+')
           56 #plt.grid()
           57 
           58 #plt.subplot(3,1,3)
           59 #plt.xlabel('Pressure gradient $\\Delta p/\\Delta z$ [Pa m$^{-1}$]')
           60 #plt.ylabel('Mean porosity $\\bar{\\phi}$ [-]')
           61 #plt.plot(dpdz, phi_bar, '+')
           62 #plt.grid()
           63 
           64 plt.tight_layout()
           65 filename = 'permeability-dpdz-vs-K.pdf'
           66 plt.savefig(filename)
           67 print(os.getcwd() + '/' + filename)
           68 plt.savefig(filename)