import openpnm as op
ws = op.Workspace()
proj = ws.new_project()
export = False

pn = op.network.Cubic(shape=[10, 10, 10], spacing=1e-4, project=proj)
geo = op.geometry.StickAndBall(network=pn, pores=pn.Ps, throats=pn.Ts)
air = op.phases.Air(network=pn, name='air')
water = op.phases.Water(network=pn, name='h2o')
hg = op.phases.Mercury(network=pn, name='hg')
phys_air = op.physics.Standard(network=pn, phase=air, geometry=geo)
phys_water = op.physics.Standard(network=pn, phase=water, geometry=geo)
phys_hg = op.physics.Standard(network=pn, phase=hg, geometry=geo)

mip = op.algorithms.Porosimetry(network=pn)
mip.setup(phase=hg)
mip.set_inlets(pores=pn.pores(['top', 'bottom']))
mip.run()
hg.update(mip.results(Pc=70000))
# mip.plot_intrusion_curve()

perm = op.algorithms.StokesFlow(network=pn)
perm.setup(phase=water)
perm.set_value_BC(pores=pn.pores('right'), values=0)
perm.set_value_BC(pores=pn.pores('left'), values=101325)
perm.run()
water.update(perm.results())
# print(perm.calc_effective_permeability())

# Add reaction to phys_air
mod = op.models.physics.source_terms.standard_kinetics
phys_air['pore.n'] = 2
phys_air['pore.A'] = -1e-5
phys_air.add_model(propname='pore.2nd_order_rxn', model=mod,
                   X='pore.concentration',
                   prefactor='pore.A', exponent='pore.n',
                   regen_mode='deferred')
rxn = op.algorithms.FickianDiffusion(network=pn)
rxn.setup(phase=air, solver='spsolve')
Ps = pn.find_nearby_pores(pores=50, r=5e-4, flatten=True)
rxn.set_source(propname='pore.2nd_order_rxn', pores=Ps)
rxn.set_value_BC(pores=pn.pores('top'), values=1)
rxn.run()
air.update(rxn.results())

fd = op.algorithms.FickianDiffusion(network=pn)
fd.setup(phase=air)
fd.set_value_BC(pores=pn.pores('left'), values=1)
fd.set_value_BC(pores=pn.pores('right'), values=0)
fd.run()
fd.calc_effective_diffusivity()

# Output network and the phases to a VTP file for visualization in Paraview
if export:
    proj.export_data(network=pn, phases=[hg, air, water], filename='output.vtp')

@@ -1,8 +1,27 @@
+import numpy as np%0A
 import o
@@ -37,67 +37,236 @@
  op%0A
-ws = op.Workspace()%0Aproj = ws.new_project()%0Aexport = False%0A
+from openpnm.models.physics import source_terms%0A%0A%0A# %25%25 Initialization: create Workspace and project objects%0Aws = op.Workspace()%0Aproj = ws.new_project()%0Anp.random.seed(9)%0A%0A# %25%25 Create network, geometry, phase, and physics objects
 %0Apn 
@@ -748,16 +748,52 @@
 y=geo)%0A%0A
+# %25%25 Perform porosimetry simulation%0A
 mip = op
@@ -941,18 +941,16 @@
 70000))%0A
-# 
 mip.plot
@@ -969,16 +969,52 @@
 urve()%0A%0A
+# %25%25 Perform Stokes flow simulation%0A
 perm = o
@@ -1227,16 +1227,15 @@
 ())%0A
-# print(
+Keff = 
 perm
@@ -1268,92 +1268,127 @@
 ty()
-)%0A%0A# Add reaction to phys_air%0Amod = op.models.physics.source_terms.standard_kinetics
+%5B0%5D%0Aprint(f%22Effective permeability: %7BKeff:.2e%7D%22)%0A%0A# %25%25 Perform reaction-diffusion simulation%0A# Add reaction to phys_air
 %0Aphy
@@ -1453,16 +1453,21 @@
 d_model(
+%0A    
 propname
@@ -1492,35 +1492,51 @@
 xn',
+%0A   
  model=
-mod,%0A               
+source_terms.standard_kinetics,%0A
     
@@ -1558,35 +1558,16 @@
 ration',
-%0A                  
  prefact
@@ -1598,31 +1598,16 @@
 ore.n',%0A
-               
     rege
@@ -1623,17 +1623,56 @@
 eferred'
-)
+%0A)%0A# Set up Fickian diffusion simulation
 %0Arxn = o
@@ -1953,16 +1953,84 @@
 lts())%0A%0A
+# %25%25 Perform pure diffusion simulation to get effective diffusivity%0A
 fd = op.
@@ -2199,16 +2199,23 @@
 d.run()%0A
+Deff = 
 fd.calc_
@@ -2237,19 +2237,69 @@
 sivity()
-%0A%0A#
+%5B0%5D%0Aprint(f%22Effective diffusivity: %7BDeff:.2e%7D%22)%0A%0A# %25%25
  Output 
@@ -2369,23 +2369,8 @@
 iew%0A
-if export:%0A    
 proj
@@ -2382,28 +2382,16 @@
 rt_data(
-network=pn, 
 phases=%5B

Refactor example.py