Question on Porosimetry

""" last edited Thur Oct 3 :  achieving mercury curve of third model (a 100*80 2D) 
##----------------------------------------------------------------------------
import openpnm as op
import scipy as sp
from math import pi
import matplotlib.pyplot as plt
import pore_size_normal
##Network----------------------------------------------------------------------
pn=op.network.Cubic(shape=[100,80,1],spacing=0.00001) #SI unit: meter
##Geometry---------------------------------------------------------------------
geom=op.geometry.GenericGeometry(network=pn,pores=pn.Ps,throats=pn.Ts)
geom['pore.diameter']=pore_size_normal.normal(8000,0,0.00001,0.000005,0.0000025)
fig1 = plt.hist(geom['pore.diameter'],bins=20, facecolor='green', alpha=0.75)
plt.xlabel('Pore Diameter(um)')
plt.ylabel('Frequency')
plt.title('Pore Diameter Histogram')
plt.grid(True)
plt.show()
P12=pn['throat.conns'] # array with numbers of pair pore
D12=geom['pore.diameter'][P12]
Dt=sp.amin(D12,axis=1 ) #axis= 0 , 1
geom['throat.diameter']=Dt
Rp=geom['pore.diameter']/2
geom['pore.volume']=(4/3)*pi*(Rp)**3
totalporevolume=sum(geom['pore.volume'])
C2C=0.0001
Rp12=Rp[P12]
geom['throat.length']=C2C-sp.sum(Rp12,axis=1)
Rt=geom['throat.diameter']/2
Lt=geom['throat.length']
geom['throat.volume']=pi*Lt*Rt**2
##Phases ----------------------------------------------------------------------
##Phases & physics ------------------------------------------------------------
hg = op.phases.Mercury(network=pn)

phys_hg = op.physics.GenericPhysics(network=pn, phase=hg, geometry=geom)

model = op.models.physics.capillary_pressure.washburn
phys_hg.add_model(propname='throat.entry_pressure',
                   model=model, contact_angle='pore.contact_angle',
                    surface_tension='pore.surface_tension')

##Algorithm--------------------------------------------------------------------
drcur=op.algorithms.Porosimetry(network=pn)
drcur.setup(phase=hg)
drcur.set_inlets(pn.pores(['left', 'right', 'top', 'bottom', 'front','back']))
drcur.run(points=100)
fig = drcur.plot_intrusion_curve()
##Finall part importing for paraview-------------------------------------------

prj=pn.project
prj.export_data(filename='untitles11',filetype='vtp')