How to use the openpnm.algorithms.StokesFlow function in openpnm

# phase
phase = op.phases.Water(network=net)
phase['pore.consistency'] = 0.00089319
phase['pore.flow_index'] = 2

# physics
phys = op.physics.GenericPhysics(network=net,
                                 phase=phase,
                                 geometry=geo)

mod1 = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys.add_model(propname='throat.hydraulic_conductance',
               model=mod1, regen_mode='normal')

# algorithms: Newtonian Stokes flow
sf = op.algorithms.StokesFlow(network=net, phase=phase)
sf.set_value_BC(pores=net.pores('front'), values=1)
sf.set_value_BC(pores=net.pores('back'), values=2)
sf.run()
phase.update(sf.results())
phase['pore.pressure_sf'] = phase['pore.pressure']

mod2 = op.models.physics.hydraulic_conductance.hagen_poiseuille_power_law
phys.add_model(propname='throat.nonNewtonian_hydraulic_conductance',
               model=mod2, regen_mode='normal')

# algorithms: Non Newtonian Stokes flow
nnsf = op.algorithms.NonNewtonianStokesFlow(network=net, phase=phase)
nnsf.set_value_BC(pores=net.pores('front'), values=1)
nnsf.set_value_BC(pores=net.pores('back'), values=2)
nnsf.run()
phase.update(nnsf.results())

model=current, regen_mode='normal')

eA_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
phys.add_model(propname='throat.diffusive_conductance.' + Na.name,
               pore_diffusivity='pore.diffusivity.' + Na.name,
               throat_diffusivity='throat.diffusivity.' + Na.name,
               model=eA_dif, regen_mode='normal')

eB_dif = op.models.physics.diffusive_conductance.ordinary_diffusion
phys.add_model(propname='throat.diffusive_conductance.' + Cl.name,
               pore_diffusivity='pore.diffusivity.' + Cl.name,
               throat_diffusivity='throat.diffusivity.' + Cl.name,
               model=eB_dif, regen_mode='normal')

# algorithms
sf = op.algorithms.StokesFlow(network=net, phase=sw)
sf.set_value_BC(pores=net.pores('back'), values=0.01)
sf.set_value_BC(pores=net.pores('front'), values=0.00)
sf.settings['rxn_tolerance'] = 1e-12
sf.run()
sw.update(sf.results())

p = op.algorithms.ChargeConservation(network=net, phase=sw)
p.set_value_BC(pores=net.pores('left'), values=0.01)
p.set_value_BC(pores=net.pores('right'), values=0.00)
p.settings['rxn_tolerance'] = 1e-12
p.settings['charge_conservation'] = 'electroneutrality'

eA = op.algorithms.NernstPlanck(network=net, phase=sw, ion=Na.name)
eA.set_value_BC(pores=net.pores('back'), values=100)
eA.set_value_BC(pores=net.pores('front'), values=90)
eA.settings['rxn_tolerance'] = 1e-12

geom['pore.area'] = sp.pi*(geom['pore.diameter']**2)/4.0
mod2 = op.models.physics.diffusive_conductance.bulk_diffusion
phys_water.add_model(propname='throat.diffusive_conductance',
                     model=mod2, diffusivity='pore.diffusivity')

phys_water.regenerate_models()

inlet = pn.pores('back')  # pore inlet
outlet = pn.pores('front')  # pore outlet

inlet2 = pn.pores('left')  # pore inlet2
outlet2 = pn.pores('right')  # pore outlet2

# ALGORITHMS
alg1 = op.algorithms.StokesFlow(network=pn, phase=water)
alg1.set_dirichlet_BC(pores=inlet, values=10)
alg1.set_dirichlet_BC(pores=outlet, values=0)
alg1.run()
water['pore.pressure'] = alg1['pore.pressure']

alg2 = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water)
alg2.set_IC(0)
alg2.set_dirichlet_BC(pores=inlet, values=2)
alg2.set_dirichlet_BC(pores=outlet, values=0)
alg2.run()

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 term to phys_air
mod = op.models.physics.generic_source_term.standard_kinetics
phys_air['pore.n'] = 2
phys_air['pore.A'] = -1e-5
phys_air.add_model(propname='pore.2nd_order_rxn', model=mod,
                   quantity='pore.concentration',
                   prefactor='pore.A', exponent='pore.n',
                   regen_mode='deferred')
rxn = op.algorithms.FickianDiffusion(network=pn)

# %%Simulate capillary pressure curve
mod = op.models.physics.capillary_pressure.washburn_slit
phys_hg.add_model(propname='throat.entry_pressure', model=mod)

mip = op.algorithms.Porosimetry(network=net)
mip.setup(phase=hg)
mip.set_inlets(net.pores(['top', 'bottom']))
mip.run(points=25)
mip.plot_intrusion_curve()

# %%Calculating permeability
mod = op.models.physics.hydraulic_conductance.hagen_poiseuille_slit
phys_water.add_model(propname='throat.hydraulic_conductance', model=mod)

alg = op.algorithms.StokesFlow(network=net, phase=water)
BC1_pores = net.pores('pore.front')
alg.set_value_BC(values=202650, pores=BC1_pores)
BC2_pores = net.pores('pore.back')
alg.set_value_BC(values=101325, pores=BC2_pores)
alg.run()
Q = alg.rate(pores=net.pores('front'))

A = (Ly*Lz)*(Ny*Nz)
L = Lx*Nx
mu = sp.mean(water['throat.viscosity'])
Kxx = Q*mu*L/(A*101325)
print("The permeability coefficient is:", Kxx/1e-15, 'mD')

# %%Calculating porosity
Vp = geo['pore.volume'][net.Ps]
Vt = geo['throat.volume'][net.Ts]

Output: array [Kewp, Kenwp]
        The value of effective permeability of defending (if there is any) and
        invading phase in the direction that is defined by flow_pores.

        Note: To account for multiphase flow, multiphase physics model is added
        and updated in each saturation (saturation is related to
        the presence of another phase). Here, the conduit_hydraulic conductance
        is used as the conductance required by stokes flow algorithm.

        """
        network = self.project.network
        self._regenerate_models()
        if self.settings['wp'] is not None:
            wp = self.project[self.settings['wp']]
            St_mp_wp = StokesFlow(network=network, phase=wp)
            St_mp_wp.setup(conductance='throat.conduit_hydraulic_conductance')
            St_mp_wp.set_value_BC(pores=flow_pores[0], values=1)
            St_mp_wp.set_value_BC(pores=flow_pores[1], values=0)
            St_mp_wp.run()
            Kewp = np.sum(abs(St_mp_wp.rate(pores=flow_pores[1])))
            self.project.purge_object(obj=St_mp_wp)
        else:
            Kewp = None
            pass
        nwp = self.project[self.settings['nwp']]
        St_mp_nwp = StokesFlow(network=network, phase=nwp)
        St_mp_nwp.set_value_BC(pores=flow_pores[0], values=1)
        St_mp_nwp.set_value_BC(pores=flow_pores[1], values=0)
        St_mp_nwp.setup(conductance='throat.conduit_hydraulic_conductance')
        St_mp_nwp.run()
        Kenwp = np.sum(abs(St_mp_nwp.rate(pores=flow_pores[1])))

from openpnm.algorithms import TransientReactiveTransport, StokesFlow
from openpnm.utils import logging
logger = logging.getLogger(__name__)


class TransientStokesFlow(TransientReactiveTransport, StokesFlow):
    r"""
    A subclass of GenericTransport to simulate Stokes flow.

    """

    def __init__(self, settings={}, **kwargs):
        super().__init__(**kwargs)
        # Apply any received settings to overwrite defaults
        self.settings.update(settings)

geom['pore.area'] = sp.pi*(geom['pore.diameter']**2)/4.0
mod2 = op.models.physics.diffusive_conductance.bulk_diffusion
phys_water.add_model(propname='throat.diffusive_conductance',
                     model=mod2, diffusivity='pore.diffusivity')

phys_water.regenerate_models()

inlet = pn.pores('back')  # pore inlet
outlet = pn.pores('front')  # pore outlet

inlet2 = pn.pores('left')  # pore inlet2
outlet2 = pn.pores('right')  # pore outlet2

# ALGORITHMS
alg1 = op.algorithms.StokesFlow(network=pn, phase=water)
alg1.set_dirichlet_BC(pores=inlet, values=5)
alg1.set_dirichlet_BC(pores=outlet, values=0)
alg1.run()

alg1b = op.algorithms.TransientStokesFlow(network=pn, phase=water)
alg1b.set_IC(0)
alg1b.set_dirichlet_BC(pores=inlet, values=5)
alg1b.set_dirichlet_BC(pores=outlet, values=0)
alg1b.run()

flow_pores: numpy array
        Boundary pores that will have constant value boundary condition to in
        StokesFlow algorithm. First element is the inlet face (pores) for flow
        of invading phase through porous media. Second element is the outlet
        face (pores).

        Output: array [Kwp, Knwp]
        The value of absolute permeability of defending (if there is any) and
        invadin phase in the direction that is defined by flow_pores.

        Note: Absolute permeability is not dependent to the phase, but here
        we just need to calculate the rate instead of all variables that are
        contributing to the darcy's law.
        """
        network = self.project.network
        St_p = StokesFlow(network=network, phase=phase)
        St_p.set_value_BC(pores=flow_pores[0], values=1)
        St_p.set_value_BC(pores=flow_pores[1], values=0)
        St_p.run()
        val = np.sum(abs(St_p.rate(pores=flow_pores[1])))
        K_abs = val
        self.project.purge_object(obj=St_p)
        return K_abs