How to get phase occupancy in Mixed Invasion Percolation algorithm?

[minzhangut]

I wrote a script to calculate the Kr curves during the drainage process using Ordinary Percolation algorithm. Now I want to do the same but using Mixed Invasion Percolation. But I am having a hard time to update the phase configuration to calculate phase occupancy and eventually to get the phase saturation and use StokesFlow to calculate effective permeability as well at each capillary pressure point.

I was trying to use the same code as the ones when I was using the Ordinary Percolation algorithm, but it is not working...

OP_1 = op.algorithms.MixedInvasionPercolation(network=pn)
OP_1.setup(phase=oil)
#OP_1.setup(late_pore_filling='pore.late_filling') 
OP_1.set_inlets(pn.pores('zmin'))
OP_1.run()
oil.update(OP_1.results(Pc=1000))
water['pore.occupancy'] = 1 - oil['pore.occupancy']
water['throat.occupancy'] = 1 - oil['throat.occupancy']
print(water['pore.occupancy'])

[jgostick]

Can you post a more full script, with perhaps a cubic network and simple geometry? That makes it easier for me to debug.

[mzdreambig]

import numpy as np
import pandas as pd
import scipy as sp
import openpnm as op
import porespy as ps
import os
import imageio

%config InlineBackend.figure_formats = ['svg']
import matplotlib.pyplot as plt
%matplotlib inline
ws = op.Workspace()
np.random.seed(10)
ws.settings["loglevel"] = 40
np.set_printoptions(precision=5)

N = 100
pn = op.network.Cubic(shape=[N, N, 1], spacing=2.5e-5)
geom = op.geometry.SpheresAndCylinders(network=pn, pores=pn.Ps, throats=pn.Ts)

water = op.phase.Water(network=pn, name='water')
water['pore.viscosity'] = 1.05e-3
water['throat.viscosity'] = 1.05e-3
water['pore.contact_angle'] = 0
water['throat.contact_angle'] = 0
water['pore.surface_tension'] = 30e-3
water['throat.surface_tension'] = 30e-3
print(water['throat.contact_angle'])

oil = op.phase.Air(network=pn, name='oil')
oil['pore.viscosity'] = 1.39e-3
oil['throat.viscosity'] = 1.39e-3
oil['pore.contact_angle'] = 180
oil['throat.contact_angle'] = 180
oil['pore.surface_tension'] = 30e-3
oil['throat.surface_tension'] = 30e-3
print(oil['throat.contact_angle'])

phys_water = op.physics.Standard(network=pn, phase=water, geometry=geom)
phys_oil = op.physics.Standard(network=pn, phase=oil, geometry=geom)

phys_water['pore.entry_pressure'] = 0 # since this is non-wetting phase displaces wetting phase, so it is bond percolation
phys_water.add_model(propname='throat.entry_pressure',
                   model=op.models.physics.capillary_pressure.washburn)
phys_water.add_model(propname='pore.pc_star',
                   model=op.models.misc.from_neighbor_throats,
                   prop='throat.entry_pressure',
                   mode='min')
phys_water.add_model(propname='pore.late_filling',
                   model=op.models.physics.multiphase.late_filling,
                   pressure='pore.pressure',
                   Pc_star='pore.pc_star',
                   #eta=0.85, Swp_star=0.21,
                   eta=1, Swp_star=0.91, # this value is not reasonable, I think?
                   regen_mode='deferred')

phys_oil['pore.entry_pressure'] = 0 # since this is non-wetting phase displaces wetting phase, so it is bond percolation
phys_oil.add_model(propname='throat.entry_pressure',
                   model=op.models.physics.capillary_pressure.washburn)
phys_oil.add_model(propname='pore.pc_star',
                   model=op.models.misc.from_neighbor_throats,
                   prop='throat.entry_pressure',
                   mode='min')
phys_oil.add_model(propname='pore.late_filling',
                   model=op.models.physics.multiphase.late_filling,
                   pressure='pore.pressure',
                   Pc_star='pore.pc_star',
                   #eta=0.85, Swp_star=0.21,
                   eta=1, Swp_star=0.91,
                   regen_mode='deferred')

OP_1 = op.algorithms.MixedInvasionPercolation(network=pn)
OP_1.setup(phase=oil)
OP_1.setup(late_pore_filling='pore.late_filling') 
OP_1.set_inlets(pn.pores('left'))
OP_1.run()

oil.update(OP_1.results(Pc=100000))
print(oil['pore.occupancy'])
#print(water['pore.occupancy'])
#print(pn['pore.volume'])

OP_1.results(Pc=100000)

I think the reason I got very small 'pore.occupancy' values instead of "0/1" is that I considered "pore.late_filling" mechanism? These small values are the actual volume of phase occupancy?
In addition, it seems that I need to convert it to be "True/False" oil['pore.occupancy'] information in order to run Stokes Flow and get the effective permeability, am I correct?

[mzdreambig]

data = OP_1.get_intrusion_data()
filter_pc = [data.Pcap[0]]
sat = [data.S_tot[0]]

print(np.array(data.Pcap).shape)
print(np.array(data.S_tot).shape)

for i, pc in enumerate(data.Pcap):
    if data.S_tot[i] - sat[-1] > 0.01:
        filter_pc.append(pc)
        sat.append(data.S_tot[i])
filter_pc.append(data.Pcap[-1])
sat.append(data.S_tot[-1])

import openpnm.models.physics as pm
phys_oil.add_model(model=pm.multiphase.conduit_conductance,
                   propname='throat.conduit_hydraulic_conductance',
                   throat_conductance='throat.hydraulic_conductance',
                  mode = 'medium')
phys_water.add_model(model=pm.multiphase.conduit_conductance,
                     propname='throat.conduit_hydraulic_conductance',
                     throat_conductance='throat.hydraulic_conductance',
                  mode = 'medium')

print('This is the total number of pores in the network:')
print(pn['pore.volume'].shape)
oil.update(OP_1.results(Pc=filter_pc[9]))
print('I "thought" this is the volume fraction that is occupied by the oil phase:')
print(oil[OP_1.settings['late_pore_filling']])
print(oil[OP_1.settings['late_pore_filling']].shape)
print('This is the actual volume that is occupied by the oil phase?')
print(oil['pore.occupancy'])
print(OP_1.results(Pc=filter_pc[9]))
print('The estimated volume fraction that is occupied by the oil phase?')
print(oil['pore.occupancy']/pn['pore.volume'])

print('The oil phase saturation from get_intrusion_data() function?')
print(sat[9])
print('The estimated oil phase saturation:')
print(np.sum(oil['pore.occupancy'])/(np.sum(pn["pore.volume"]) + np.sum(pn["throat.volume"])))
print(np.sum(oil['pore.occupancy'])/(np.sum(pn["pore.volume"]) ))

This is what I got:

This is the total number of pores in the network:
(12551,)
I "thought" this is the volume fraction that is occupied by the oil phase:
[0.73731 0.56285 0.65605 ... 0.57875 0.47461 0.29512]
(12551,)
This is the actual volume that is occupied by the oil phase?
[2.03558e-13 0.00000e+00 4.10739e-14 ... 0.00000e+00 0.00000e+00
 0.00000e+00]
{'pore.occupancy': array([2.03558e-13, 0.00000e+00, 4.10739e-14, ..., 0.00000e+00,
       0.00000e+00, 0.00000e+00]), 'throat.occupancy': array([3.31401e-24, 3.88539e-24, 4.11394e-24, ..., 7.01742e-16,
       9.14209e-25, 9.14209e-25])}
The estimated volume fraction that is occupied by the oil phase?
[0.73731 0.      0.65605 ... 0.      0.      0.     ]
The oil phase saturation from get_intrusion_data() function?
0.5958578745455175
The estimated oil phase saturation:
0.5452183017659206
0.5854535254627582

I got confused. I am wondering which calculation I should use to estimate the oil phase saturation? Did I make some mistakes when I defined the oil phase?
oil['pore.occupancy']/pn['pore.volume'] is the same as oil[OP_1.settings['late_pore_filling']] when its value is larger than 0.6?

[jgostick]

Hi @minzhangut, I am sorry, but I don't have the bandwidth to help you with this right now. If you're still having trouble in a week or so, please ping me again either by email or by posting here.

[EnochMayor]

Hi @minzhangut, were you able to figure it out and resolve this issue?
I have a similar situation.