Size factor models return negative numbers for conduits with "engulfing" pores

[ma-sadeghi]

When one pore center lies within another pore, our size factor models break down, including cones_and_cylinders.

[LS-Maxwell]

hi, I also have a similar problem. And I write a small code to recover it.:

import numpy as np
import scipy as sp
import openpnm as op
import matplotlib.pyplot as plt
import math
np.random.seed(10)
ws = op.Workspace()
ws.settings["loglevel"] = 40

length_x=0.1
length_y=0.1
length_z=0.1
points_num=33
pore_scale=0.0001       #---pore scale---#
pn = op.network.Voronoi(shape=[length_x, length_y, length_z], points=points_num)  

fig, ax = plt.subplots(figsize=(5, 5))
op.topotools.plot_connections(network=pn, ax=ax)
op.topotools.plot_coordinates(network=pn, c='r', s=75, ax=ax)
op.topotools.plot_coordinates(pn, pn.pores('back'), c='blue', marker='*', 
                              markersize=50, ax=ax)  
op.topotools.plot_coordinates(pn, pn.pores('front'), c='cyan', marker='o', 
                              markersize=50, ax=ax) 

#------------------------------to cut unhealth pore----------------------------#
array_k = {}
pn2=pn
k2=0
for j in range(0,60,1):       #---60 is fine, it can enlarge---#
    if pn.Np==k2 :
        break
    pn2=pn
    k2=pn2.Np
    pn.regenerate_models()
    for i in range(0,pn.Np,1):    #----to check unhealth pore, to avoid pores too short----#
        too_close_pore = pn.find_nearby_pores(pores=[i], r=pore_scale)
        if len(too_close_pore[0]) > 0:
            too_close_pore=np.append(too_close_pore,[i])        #---put the points into---#
            array_k=too_close_pore
            #print(too_close_pore[0])
            op.topotools.merge_pores(network=pn, pores=array_k, labels=['merged'])
            pn.regenerate_models()
            k1=pn.Np
            break
pn.regenerate_models()
#------------------------------to cut unhealth pore----------------------------#

[LS-Maxwell]

I also have another problem. For example, if I want to set different water viscosity in model, where the top viscosity is 0.01 ,the bottom viscosity is 0.02, and the other location is 0.03. How can I bulit this heterogeneous property?

[LS-Maxwell]

I also have another problem. For example, if I want to set different water viscosity in model, where the top viscosity is 0.01 ,the bottom viscosity is 0.02, and the other location is 0.03. How can I bulit this heterogeneous property?

Well，I just write a studip for-if cyclic sentence to meet my needs. But if there are so many pores, this code will seriously slow me program down.

for i in range(0,pn.Np,1):
    z_coord=pn['pore.coords'][i][2]
    if z_coord==0.0:
        water['pore.pressure'][i]=44444.0
    else:
        water['pore.pressure'][i]=11111.0

[ma-sadeghi]

I also have another problem. For example, if I want to set different water viscosity in model, where the top viscosity is 0.01 ,the bottom viscosity is 0.02, and the other location is 0.03. How can I bulit this heterogeneous property?

You don't need the for loop. You can achieve that by creating masks.

Let's say that you want to set viscosity in the following way:

• z = 0.0 -> nu = 0.1
• z = 5.0 -> nu = 0.5
• otherwise -> nu = 0.25

X, Y, Z = net.coords.T

m1 = Z == 0.0
m2 = Z == 5.0
m3 = ~(m1 + m2)

water["pore.viscosity"][m1] = 0.1
water["pore.viscosity"][m2] = 0.5
water["pore.viscosity"][m3] = 0.25

Regarding the other issue (the center of one pore lies within the other pore), we'll try to address it soon (if we come up with a solution).

[LS-Maxwell]

I also have another problem. For example, if I want to set different water viscosity in model, where the top viscosity is 0.01 ,the bottom viscosity is 0.02, and the other location is 0.03. How can I bulit this heterogeneous property?

You don't need the for loop. You can achieve that by creating masks.

Let's say that you want to set viscosity in the following way:

• z = 0.0 -> nu = 0.1
• z = 5.0 -> nu = 0.5
• otherwise -> nu = 0.25

X, Y, Z = net.coords.T

m1 = Z == 0.0
m2 = Z == 5.0
m3 = ~(m1 + m2)

water["pore.viscosity"][m1] = 0.1
water["pore.viscosity"][m2] = 0.5
water["pore.viscosity"][m3] = 0.25

Regarding the other issue (the center of one pore lies within the other pore), we'll try to address it soon (if we come up with a solution).

Wow, thanks for the easier way . I will work overtime tonight writing code to solve my calculations.

[ma-sadeghi]

@jgostick Here's a suggestions: I think this situation mostly (if not always) occurs when working with extracted networks. Why not calculating conduit lengths directly in snow2 (as the distance between a pore center and throat centroid)?

[jgostick]

Perhaps we could add a new pore.safe_diameter which is the largest sphere that can safely be used in any conductance model?

[ma-sadeghi]

Note to self: here's a workaround until we formally patch it up:

def cones_and_cylinders(
    target,
    pore_diameter="pore.diameter",
    throat_diameter="throat.diameter"
):
    L_ctc = target['throat.spacing']
    D1, Dt, D2 = target.get_conduit_data(
        poreprop=pore_diameter,
        throatprop=throat_diameter
    ).T

    L1 = D1 / 2
    L2 = D2 / 2

    # Handle throats w/ overlapping pores
    _L1 = (4 * L_ctc**2 + D1**2 - D2**2) / (8 * L_ctc)
    mask = L_ctc - 0.5 * (D1 + D2) < 0
    if mask.any():
        L1[mask] = _L1[mask]
        L2[mask] = (L_ctc - L1)[mask]

    # Handle engulfing pores
    XYp = target.network.coords[target.network.conns][target.throats(to_global=True)]
    XYt = target["throat.global_peak"]
    _L1, _L2 = _np.linalg.norm(XYp - XYt[:, None], axis=2).T
    mask = L_ctc**2 <= (D1/2)**2 - (D2/2)**2
    if mask.any():
        L1[mask] = _L1[mask]
        L2[mask] = _L2[mask]
    mask = L_ctc**2 <= (D2/2)**2 - (D1/2)**2
    if mask.any():
        L1[mask] = _L1[mask]
        L2[mask] = _L2[mask]

    Lt = _np.maximum(L_ctc - (L1 + L2), 1e-15)
    return _np.vstack((L1, Lt, L2)).T