How to handle extra boundary pores added by SNOW extraction method from porespy?

test.py

def generate_network(shape, radius, ncyl, phimax, thetamax):

    ws = op.Workspace()
    ws.clear()

    proj = ws.new_project(name='proj1')

    im = psp.generators.cylinders(shape=shape, radius=radius, ncylinders=ncyl, phi_max=phimax, theta_max=thetamax)
    im = psp.filters.fill_blind_pores(im)
    im = psp.filters.trim_floating_solid(im)

    # extract network
    net = psp.networks.snow(im=im, voxel_size=1e-6)

    pn = op.network.GenericNetwork(name='net1', project=proj)
    pn.update(net)

    op.utils.Workspace.save_workspace(ws, filename='pnm_networks/testws.pnm')


def add_boundary_pores(pn):

    op.topotools.add_boundary_pores(pn, pores=pn.pores('top'), offset=[0, 0, max(pn['pore.diameter'])],
                                    apply_label='top_inlet')
    op.topotools.add_boundary_pores(pn, pores=pn.pores('bottom'), offset=[0, 0, -max(pn['pore.diameter'])],
                                    apply_label='bottom_inlet')
    op.topotools.add_boundary_pores(pn, pores=pn.pores('front'), offset=[0, -max(pn['pore.diameter']), 0],
                                    apply_label='front_inlet')
    op.topotools.add_boundary_pores(pn, pores=pn.pores('back'), offset=[0, max(pn['pore.diameter']), 0],
                                    apply_label='back_inlet')
    op.topotools.add_boundary_pores(pn, pores=pn.pores('top'), offset=[-max(pn['pore.diameter']), 0, 0],
                                    apply_label='left_inlet')
    op.topotools.add_boundary_pores(pn, pores=pn.pores('right'), offset=[max(pn['pore.diameter']), 0, 0],
                                    apply_label='right_inlet')
    
    # make geometry obect for inlet pores
    PsB = pn.pores('*inlet')
    TsB = pn.find_neighbor_throats(pores=PsB, mode='xor', flatten=True)
    geom_b = op.geometry.GenericGeometry(network=pn, pores=PsB, throats=TsB, name='boun1')

    # add properties to boundary pores and throats
    geom_b['pore.volume'] = 0.0
    geom_b['pore.inscribed_diameter'] = max(pn['pore.inscribed_diameter'])
    geom_b['pore.diameter'] = max(pn['pore.diameter'])
    geom_b['pore.equivalent_diameter'] = max(pn['pore.equivalent_diameter'])
    geom_b['pore.extended_diameter'] = max(pn['pore.extended_diameter'])
    geom_b['pore.surface_area'] = max(pn['pore.surface_area'])
    geom_b['pore.area'] = max(pn['pore.area'])

    geom_b['throat.volume'] = 0.0
    geom_b['throat.diameter'] = np.nanmax(pn['throat.diameter'])
    geom_b['throat.inscribed_diameter'] = max(pn['throat.inscribed_diameter'])
    geom_b['throat.area'] = max(pn['throat.area'])
    geom_b['throat.perimeter'] = max(pn['throat.perimeter'])
    geom_b['throat.equivalent_diameter'] = max(pn['throat.equivalent_diameter'])
    geom_b['throat.total_length'] = max(pn['throat.total_length'])
    geom_b['throat.length'] = max(pn['throat.length'])
    geom_b['throat.direct_length'] = max(pn['throat.direct_length'])
    geom_b['throat.conduit_lengths.pore1'] = max(pn['throat.conduit_lengths.pore1'])
    geom_b['throat.conduit_lengths.pore2'] = max(pn['throat.conduit_lengths.pore2'])
    geom_b['throat.conduit_lengths.throat'] = max(pn['throat.conduit_lengths.throat'])

    return geom_b


if __name__ == '__main__':

    if not os.path.isfile('pnm_networks/testws.pnm'):  # check if a file is present 
        generate_network([1000, 1000, 250], 5, 200, 90, 90)  # if no then generate a network 

    ws = op.Workspace()
    ws.clear()
    ws.load_project('pnm_networks/testws.pnm')

    pn = ws['proj1']['net1']  # load project 

    geom = op.geometry.Imported(network=pn, name='geo1')
    geomb = add_boundary_pores(pn)  # <------------------------------ adds boundary pores to network 

    dis_pore = pn.check_network_health()

    if dis_pore['trim_pores'] is not None:
        op.topotools.trim(pn, pores=dis_pore['trim_pores'])
        dis_pore = pn.check_network_health()
    
    # phases
    air = op.phases.Air(network=pn, name='air')
    water = op.phases.Water(network=pn, name='water')
    water['pore.contact_angle'] = 110
    
    # physicss 
    phys_water = op.physics.Standard(network=pn, phase=water, geometry=geom)
    phys_air = op.physics.Standard(network=pn, phase=air, geometry=geom)
    phys_water_b = op.physics.Standard(network=pn, phase=water, geometry=geomb)
    phys_air_b = op.physics.Standard(network=pn, phase=air, geometry=geomb)
    
    # run ordinary percolation 
    OP1 = op.algorithms.OrdinaryPercolation(network=pn)
    OP1.set_inlets(pores=pn.pores('top_inlet'))
    OP1.set_outlets(pores=pn.pores('bottom'))
    OP1.settings['trapping'] = True
    OP1.setup(phase=water, pore_volume='pore.volume', throat_volume='throat.volume', access_limited=True)
    OP1.run(points=100)

    # plt.figure()
    OP1.plot_intrusion_curve()

    data = OP1.get_intrusion_data()
    sat = data.Snwp

    def update_phys_phase(results):
        water['pore.occupancy'] = results['pore.occupancy']
        air['pore.occupancy'] = 1 - results['pore.occupancy']
        water['throat.occupancy'] = results['throat.occupancy']
        air['throat.occupancy'] = 1 - results['throat.occupancy']

        # internal physics
        phys_air.add_model(model=pm.multiphase.conduit_conductance,
                           propname='throat.conduit_hydraulic_conductance',
                           throat_conductance='throat.hydraulic_conductance',
                           mode='strict')
        phys_water.add_model(model=pm.multiphase.conduit_conductance,
                             propname='throat.conduit_hydraulic_conductance',
                             throat_conductance='throat.hydraulic_conductance',
                             mode='strict')
        phys_air.add_model(model=pm.multiphase.conduit_conductance,
                           propname='throat.conduit_diffusive_conductance',
                           throat_conductance='throat.diffusive_conductance',
                           mode='strict')
        phys_water.add_model(model=pm.multiphase.conduit_conductance,
                             propname='throat.conduit_diffusive_conductance',
                             throat_conductance='throat.diffusive_conductance',
                             mode='strict')

        # boundary physics
        phys_air_b.add_model(model=pm.multiphase.conduit_conductance,
                             propname='throat.conduit_hydraulic_conductance',
                             throat_conductance='throat.diffusive_conductance',
                             mode='strict')
        phys_water_b.add_model(model=pm.multiphase.conduit_conductance,
                               propname='throat.conduit_hydraulic_conductance',
                               throat_conductance='throat.diffusive_conductance',
                               mode='strict')
        phys_air_b.add_model(model=pm.multiphase.conduit_conductance,
                             propname='throat.conduit_diffusive_conductance',
                             throat_conductance='throat.diffusive_conductance',
                             mode='strict')
        phys_water_b.add_model(model=pm.multiphase.conduit_conductance,
                               propname='throat.conduit_diffusive_conductance',
                               throat_conductance='throat.diffusive_conductance',
                               mode='strict')


    bounds = [['front', 'back'], ['left', 'right'], ['top', 'bottom']]
    diff_air = {'0': [], '1': [], '2': []}
    diff_water = {'0': [], '1': [], '2': []}
    perm_air = {'0': [], '1': [], '2': []}
    perm_water = {'0': [], '1': [], '2': []}

    max_Pc = max(OP1['throat.invasion_pressure'])

    pore_volumes = pn['pore.volume']
    throat_volumes = pn['throat.volume']
    tot_vol = np.sum(pore_volumes) + np.sum(throat_volumes)

    SPP_air = [None, None, None]
    SPP_water = [None, None, None]
    SPD_air = [None, None, None]
    SPD_water = [None, None, None]

    atm1 = 1
    atm2 = 0

    for bound_increment in range(len(bounds)):
        BC1_pores = pn.pores(labels=bounds[bound_increment][0])
        BC2_pores = pn.pores(labels=bounds[bound_increment][1])

        # setup and run stokes and fickian diffusion
        SPST_air = op.algorithms.StokesFlow(network=pn, phase=air)
        SPST_air.setup(conductance='throat.hydraulic_conductance')
        SPST_air.set_value_BC(values=atm1, pores=BC1_pores)
        SPST_air.set_value_BC(values=atm2, pores=BC2_pores)
        SPST_air.run()

        SPST_water = op.algorithms.StokesFlow(network=pn, phase=water)
        SPST_water.setup(conductance='throat.hydraulic_conductance')
        SPST_water.set_value_BC(values=atm1, pores=BC1_pores)
        SPST_water.set_value_BC(values=atm2, pores=BC2_pores)
        SPST_water.run()

        SPFD_air = op.algorithms.FickianDiffusion(network=pn, phase=air)
        SPFD_air.setup(conductance='throat.diffusive_conductance')
        SPFD_air.set_value_BC(values=atm1, pores=BC1_pores)
        SPFD_air.set_value_BC(values=atm2, pores=BC2_pores)
        SPFD_air.run()

        SPFD_water = op.algorithms.FickianDiffusion(network=pn, phase=water)
        SPFD_water.setup(conductance='throat.diffusive_conductance')
        SPFD_water.set_value_BC(values=atm1, pores=BC1_pores)
        SPFD_water.set_value_BC(values=atm2, pores=BC2_pores)
        SPFD_water.run()

        # multiplying the one length by 1e-6 has the same impact as multiplying all lengths by the scale factor
        SPP_air[bound_increment] = SPST_air.calc_effective_permeability(
            domain_area=1000*1000*1e-6, domain_length=250
        )
        SPP_water[bound_increment] = SPST_water.calc_effective_permeability(
            domain_area=1000*1000*1e-6, domain_length=250
        )
        SPD_air[bound_increment] = SPFD_air.calc_effective_diffusivity(
            domain_area=1000*1000*1e-6, domain_length=250
        )
        SPD_water[bound_increment] = SPFD_water.calc_effective_diffusivity(
            domain_area=1000*1000*1e-6, domain_length=250
        )

        pn.project.purge_object(obj=SPST_air)
        pn.project.purge_object(obj=SPST_water)
        pn.project.purge_object(obj=SPFD_air)
        pn.project.purge_object(obj=SPFD_water)

    for Pc in data.Pcap:
        update_phys_phase(OP1.results(Pc=Pc))

        for bound_increment in range(len(bounds)):
            BC1_pores = pn.pores(labels=bounds[bound_increment][0])
            BC2_pores = pn.pores(labels=bounds[bound_increment][1])

            # setup and run multiphase stokes flow and fickian diffusion
            MPST_air = op.algorithms.StokesFlow(network=pn, phase=air)
            MPST_air.setup(conductance='throat.conduit_hydraulic_conductance')
            MPST_air.set_value_BC(values=atm1, pores=BC1_pores)
            MPST_air.set_value_BC(values=atm2, pores=BC2_pores)
            MPST_air.run()

            MPST_water = op.algorithms.StokesFlow(network=pn, phase=water)
            MPST_water.setup(conductance='throat.conduit_hydraulic_conductance')
            MPST_water.set_value_BC(values=atm1, pores=BC1_pores)
            MPST_water.set_value_BC(values=atm2, pores=BC2_pores)
            MPST_water.run()

            MPFD_air = op.algorithms.FickianDiffusion(network=pn, phase=air)
            MPFD_air.setup(conductance='throat.conduit_diffusive_conductance')
            MPFD_air.set_value_BC(values=atm1, pores=BC1_pores)
            MPFD_air.set_value_BC(values=atm2, pores=BC2_pores)
            MPFD_air.run()

            MPFD_water = op.algorithms.FickianDiffusion(network=pn, phase=water)
            MPFD_water.setup(conductance='throat.conduit_diffusive_conductance')
            MPFD_water.set_value_BC(values=atm1, pores=BC1_pores)
            MPFD_water.set_value_BC(values=atm2, pores=BC2_pores)
            MPFD_water.run()

            MPair_eff_perm = MPST_air.calc_effective_permeability(
                domain_area=1000*1000*1e-6, domain_length=250
            )
            MPwater_eff_perm = MPST_water.calc_effective_permeability(
                domain_area=1000*1000*1e-6, domain_length=250
            )
            MPair_eff_dif = MPFD_air.calc_effective_diffusivity(
                domain_area=1000*1000*1e-6, domain_length=250
            )
            MPwater_eff_dif = MPFD_water.calc_effective_diffusivity(
                domain_area=1000*1000*1e-6, domain_length=250
            )

            rel_eff_perm_air = MPair_eff_perm / SPP_air[bound_increment]
            rel_eff_perm_water = MPwater_eff_perm / SPP_water[bound_increment]
            rel_eff_dif_air = MPair_eff_dif / SPD_air[bound_increment]
            rel_eff_dif_water = MPwater_eff_dif / SPD_water[bound_increment]

            perm_air[str(bound_increment)].append(rel_eff_perm_air)
            perm_water[str(bound_increment)].append(rel_eff_perm_water)
            diff_air[str(bound_increment)].append(rel_eff_dif_air)
            diff_water[str(bound_increment)].append(rel_eff_dif_water)

            pn.project.purge_object(obj=MPFD_air)
            pn.project.purge_object(obj=MPFD_water)
            pn.project.purge_object(obj=MPST_air)
            pn.project.purge_object(obj=MPST_water)