Throat occupancy does not change in MixedInvasionPercolation.result()

[saeidabd]

@jgostick @TomTranter Dear OpenPNM Team, I have problem in finding relative permeability. I used MixedInvasionPercolation.result() to update occupancy of pores and throat then set a little value for hydraulic conductance in desired place and run stock flow but throat occupancy remains zero in all capilarry pressures. what can I do?

`import porespy as ps import openpnm as op import numpy as np import matplotlib.pyplot as plt

#importing network and geometry from Berea Sandston----------------------------
fromraw=np.fromfile( path ,dtype=np.uint8)
im=fromraw.shape=(300,300,300)
im = ~np.array(fromraw, dtype=bool)
print(ps.metrics.porosity(im))
resolution = 4.9e-6
net = ps.networks.snow(im=im, voxel_size=resolution)
ws = op.Workspace()
proj = op.Project()
pn = op.network.GenericNetwork(name='Berea', project=proj)
pn.update(net)
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
h = pn.check_network_health()

water = op.phases.Water(network=pn,name='water')
oil = op.phases.GenericPhase(network=pn,name='oil')

oil['pore.viscosity'] = 0.001
oil['pore.contact_angle'] = 120
water['pore.contact_angle']= 60
oil['pore.surface_tension'] = 0.04
water["pore.surface_tension"]=0.04
water["pore.temperature"] = 350
water.regenerate_models()
oil.regenerate_models()

phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=pn)
phys_oil = op.physics.GenericPhysics(network=pn, phase=oil, geometry=pn)
mod = op.models.physics.capillary_pressure.washburn
water.add_model(propname='throat.entry_pressure', model=mod)
oil.add_model(propname='throat.entry_pressure', model=mod)

mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys_water.add_model(propname='throat.hydraulic_conductance', model=mod)
phys_oil.add_model(propname='throat.hydraulic_conductance', model=mod)

IP_withdrawal = op.algorithms.MixedInvasionPercolation(network=pn)
IP_withdrawal.setup(phase=oil)
IP_withdrawal.set_inlets(pores=pn.pores('right'))
IP_withdrawal.run()
IP_withdrawal.set_outlets(pores=pn.pores(['left']))
IP_withdrawal.apply_trapping()
injection_data = IP_withdrawal.get_intrusion_data()

Pcdata = {}
Pcdata['capillary_pressure'] = injection_data.Pcap  
Pcdata['oil_phase_saturation'] = injection_data.S_tot
Pcdata['water_phase_saturation'] = []
for i in injection_data.S_tot:
    Pcdata['water_phase_saturation'].append(1-i)

FlowAlg_water = op.algorithms.StokesFlow(network=pn, phase=water)
FlowAlg_water.set_value_BC(pores=pn.pores('left'), values=202650)
FlowAlg_water.set_value_BC(pores=pn.pores('right'), values=101325)
FlowAlg_water.run()

FlowAlg_oil = op.algorithms.StokesFlow(network=pn, phase=oil)
FlowAlg_oil.set_value_BC(pores=pn.pores('left'), values=202650)
FlowAlg_oil.set_value_BC(pores=pn.pores('right'), values=101325)
FlowAlg_oil.run()

Qtot_water = FlowAlg_water.rate(pores=pn.pores('left'))
Qtot_oil = FlowAlg_oil.rate(pores=pn.pores('left'))

pnVolume= np.sum(pn['pore.volume']) + np.sum(pn['throat.volume'])

Krdata= {}
Krdata['Krw'] = []
Krdata['Krnw'] = []
Krdata['Sw'] = []
for Pc in Pcdata['capillary_pressure']:
    oil.update(IP_withdrawal.results(Pc))
    Po=[]
    To=[]
    for i in oil['throat.occupancy']:
        if i==0:
            To.append(False)
        else:
            To.append(True)
    for i in oil['pore.occupancy']:
        if i==0:
            Po.append(False)
        else:
            Po.append(True)
    Po=np.array(Po)
    To=np.array(To)
    phys_oil['throat.hydraulic_conductance'][~To] = 1e-20
    FlowAlg_oil.run()
    Q_partial = FlowAlg_oil.rate(pores=pn.pores('left'))
    Q_partial[Q_partial<0] = 0
    Krdata['Krnw'].append(Q_partial[0]/Qtot_oil[0])
    Vt = np.sum(pn['throat.volume'][To])
    Vp = np.sum(pn['pore.volume'][Po])
    Krdata['Sw'].append(1-(Vp+Vt)/pnVolume)
    oil.regenerate_models()
    phys_oil.regenerate_models()

    phys_water['throat.hydraulic_conductance'][To] = 1e-20
    FlowAlg_water.run()
    Q_partial = FlowAlg_water.rate(pores=pn.pores('left'))
    Q_partial[Q_partial<0] = 0
    Krdata['Krw'].append(Q_partial[0]/Qtot_water[0])
    water.regenerate_models()
    phys_water.regenerate_models()

plt.figure()
plt.plot(Krdata['Sw'],Krdata['Krnw']  ,'y*-',
               label='Non-wet phase (oil)')
plt.plot(Krdata['Sw'],Krdata['Krw']  ,'b*-',
               label='Wet phase (water)')
plt.ylabel('Relative permeability')
plt.xlabel('Saturation \n Porous Volume Fraction occupied by wet phase')
plt.grid(True)
plt.legend()
plt.show()`

[TomTranter]

Hi @saeidabd, your problem is that the results are given as a fraction of the pore space that is filled at this pressure allowing for fractional pore filling and corner flow models. The following code should work. I also reduced number of inlets and number of points to evaluate.

import openpnm as op
import numpy as np
import matplotlib.pyplot as plt
import os
from pathlib import Path

INST_DIR = os.path.dirname(os.path.abspath(op.__file__))
ROOT_DIR = Path(INST_DIR).parent
BEREA_DIR = os.path.join(ROOT_DIR, 'tests\\fixtures\\ICL-Sandstone(Berea)')
project = op.io.Statoil.load(path=BEREA_DIR, prefix='Berea')

pn = project.network
pn['pore.diameter'] = pn['pore.radius'] * 2
pn['throat.diameter'] = pn['throat.radius'] * 2
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
h = pn.check_network_health()

water = op.phases.Water(network=pn, name='water')
oil = op.phases.GenericPhase(network=pn, name='oil')

oil['pore.viscosity'] = 0.001
oil['pore.contact_angle'] = 120
water['pore.contact_angle'] = 60
oil['pore.surface_tension'] = 0.04
water["pore.surface_tension"] = 0.04
water["pore.temperature"] = 350
water.regenerate_models()
oil.regenerate_models()

phys_water = op.physics.GenericPhysics(network=pn, phase=water, geometry=pn)
phys_oil = op.physics.GenericPhysics(network=pn, phase=oil, geometry=pn)
mod = op.models.physics.capillary_pressure.washburn
water.add_model(propname='throat.entry_pressure', model=mod)
oil.add_model(propname='throat.entry_pressure', model=mod)

mod = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys_water.add_model(propname='throat.hydraulic_conductance', model=mod)
phys_oil.add_model(propname='throat.hydraulic_conductance', model=mod)

IP_withdrawal = op.algorithms.MixedInvasionPercolation(network=pn)
IP_withdrawal.setup(phase=oil)
# Select every 4th pore as inlet to prevent closing it off straight away
IP_withdrawal.set_inlets(pores=pn.pores('inlets')[::4])
IP_withdrawal.run()
IP_withdrawal.set_outlets(pores=pn.pores(['outlets']))
IP_withdrawal.apply_trapping()
injection_data = IP_withdrawal.get_intrusion_data()

Pcdata = {}
Pcdata['capillary_pressure'] = injection_data.Pcap
Pcdata['oil_phase_saturation'] = injection_data.S_tot
Pcdata['water_phase_saturation'] = []
for i in injection_data.S_tot:
    Pcdata['water_phase_saturation'].append(1-i)

FlowAlg_water = op.algorithms.StokesFlow(network=pn, phase=water)
FlowAlg_water.set_value_BC(pores=pn.pores('inlets'), values=202650)
FlowAlg_water.set_value_BC(pores=pn.pores('outlets'), values=101325)
FlowAlg_water.run()

FlowAlg_oil = op.algorithms.StokesFlow(network=pn, phase=oil)
FlowAlg_oil.set_value_BC(pores=pn.pores('inlets'), values=202650)
FlowAlg_oil.set_value_BC(pores=pn.pores('outlets'), values=101325)
FlowAlg_oil.run()

Qtot_water = FlowAlg_water.rate(pores=pn.pores('inlets'))
Qtot_oil = FlowAlg_oil.rate(pores=pn.pores('inlets'))

pnVolume = np.sum(pn['pore.volume']) + np.sum(pn['throat.volume'])

Krdata = {}
Krdata['Krw'] = []
Krdata['Krnw'] = []
Krdata['Sw'] = []
for Pc in np.linspace(0, 20000, 21):
    res = IP_withdrawal.results(Pc)
    Po = res['pore.occupancy'] > 0
    To = res['throat.occupancy'] > 0
    oil.update({'pore.occupancy': Po,
                'throat.occupancy': To})
    phys_oil['throat.hydraulic_conductance'][~To] = 1e-20
    FlowAlg_oil.run()
    Q_partial = FlowAlg_oil.rate(pores=pn.pores('inlets'))
    Q_partial[Q_partial < 0] = 0
    Krdata['Krnw'].append(Q_partial[0]/Qtot_oil[0])
    Vt = np.sum(pn['throat.volume'][To])
    Vp = np.sum(pn['pore.volume'][Po])
    Krdata['Sw'].append(1-(Vp+Vt)/pnVolume)
    oil.regenerate_models()
    phys_oil.regenerate_models()
    phys_water['throat.hydraulic_conductance'][To] = 1e-20
    FlowAlg_water.run()
    Q_partial = FlowAlg_water.rate(pores=pn.pores('inlets'))
    Q_partial[Q_partial < 0] = 0
    Krdata['Krw'].append(Q_partial[0]/Qtot_water[0])
    water.regenerate_models()
    phys_water.regenerate_models()

plt.figure()
plt.plot(Krdata['Sw'], Krdata['Krnw'], 'y*-',
         label='Non-wet phase (oil)')
plt.plot(Krdata['Sw'], Krdata['Krw'], 'b*-',
         label='Wet phase (water)')
plt.ylabel('Relative permeability')
plt.xlabel('Saturation \n Porous Volume Fraction occupied by wet phase')
plt.grid(True)
plt.legend()
plt.show()

[saeidabd]

Dear @TomTranter I saw this part of source code of MixedInvasionPercolation that is for result function:

    def results(self, Pc):
        r"""
        Places the results of the IP simulation into the Phase object.

        Parameters
        ----------
        Pc : float
            Capillary Pressure at which phase configuration was reached

        """
        if Pc is None:
            results = {'pore.invasion_sequence':
                       self['pore.invasion_sequence'],
                       'throat.invasion_sequence':
                       self['throat.invasion_sequence']}
        else:
            phase = self.project.find_phase(self)
            net = self.project.network
            inv_p = self['pore.invasion_pressure'].copy()
            inv_t = self['throat.invasion_pressure'].copy()
            # Handle trapped pores and throats by moving their pressure up
            # to be ignored
            if np.sum(self['pore.invasion_sequence'] == -1) > 0:
                inv_p[self['pore.invasion_sequence'] == -1] = Pc + 1
            if np.sum(self['throat.invasion_sequence'] == -1) > 0:
                inv_t[self['throat.invasion_sequence'] == -1] = Pc + 1
            p_inv = inv_p <= Pc
            t_inv = inv_t <= Pc

            if self.settings['late_pore_filling']:
                # Set pressure on phase to current capillary pressure
                phase['pore.pressure'] = Pc
                # Regenerate corresponding physics model
                for phys in self.project.find_physics(phase=phase):
                    phys.regenerate_models(self.settings['late_pore_filling'])
                # Fetch partial filling fraction from phase object (0->1)
                frac = phase[self.settings['late_pore_filling']]
                p_vol = net['pore.volume']*frac
            else:
                p_vol = net['pore.volume']
            if self.settings['late_throat_filling']:
                # Set pressure on phase to current capillary pressure
                phase['throat.pressure'] = Pc
                # Regenerate corresponding physics model
                for phys in self.project.find_physics(phase=phase):
                    phys.regenerate_models(self.settings['late_throat_filling'])
                # Fetch partial filling fraction from phase object (0->1)
                frac = phase[self.settings['late_throat_filling']]
                t_vol = net['throat.volume']*frac
            else:
                t_vol = net['throat.volume']
            results = {'pore.occupancy': p_inv*p_vol,
                       'throat.occupancy': t_inv*t_vol}
        return results

It sounds throat occupancy should be update. I try statoil and snow algorithm both (without updating throat like you did for me), snow algorithm has problem and throat occupancy does not change but statoil one works correctly. I use snow algorithm because the statoil one has big error in calculating permeability. I really don't know what should I do, any Idea?

[TomTranter]

I used the statoil file as I don't have your raw file. I don't really understand what your problem is now but I think you should spend some time debugging your code. Maybe throat volume is zero in the snow algorithm so throat occupancy is zero?

[saeidabd]

@TomTranter I did check it and you are right, in the snow throat volumes are zero. I attach the file for snow algorithm here. Thank you for your responding.

[jgostick]

I will also add one comment here: The snow algorithm puts some questionable values into boundary pores. We are working to remedy this, but I would suggest either not including boundary pores in the snow algorithm, or deleting them once you're in openpnm, and seeing if that helps.

[saeidabd]

@TomTranter @jgostick Thank you dear @jgostick , it didn't really help. As Tom pointed to me finding out zero volume of throat, I add throat.volume dictionary manually again by their length and diameter and it got me some curves. I will close this issue.