Source code for porespy.filters._size_seq_satn

import numpy as np
from scipy.stats import rankdata

from porespy.tools import make_contiguous

__all__ = [
    'size_to_seq',
    'size_to_satn',
    'seq_to_satn',
    'pc_to_satn',
    'pc_to_seq',
    'satn_to_seq',
    'size_to_pc',
]


[docs] def size_to_seq(size, im=None, bins=None, mode='drainage'): r""" Converts a size map to a sequence map Parameters ---------- size : ndarray The image containing invasion size in each voxel. Values of 0 are assumed to be solid (if ``im`` is not given) and values of -1 are assumed to be uninvaded. im : ndarray, optional A binary image of the porous media, with ``True`` indicating the void space and ``False`` indicating the solid phase. If not given then it is assumed that the solid is identified as ``size == 0``. bins : array_like or int (optional) The bins to use when converting sizes to sequence. The default is to create 1 bin for each unique value in ``size`` (except for -1 and 0). If an **int** is supplied, it is interpreted as the number of bins between 1 and the maximum value in ``size``. If an array is supplied it is used as the bins directly. mode : str Controls how the sizes are converted to a sequence. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' The sizes are assumed to have been filled from largest to smallest, ignoring 0's and -1's 'imbibition' The sizes are assumed to have been filled from smallest to largest, ignoring 0's and -1's ============= ============================================================== Returns ------- seq : ndarray An ndarray the same shape as ``size`` with invasion size values replaced by the invasion sequence, according to the specified `mode`. Any uninvaded voxels, indicated by -1 in ``size`` will be indicated by -1 in ``seq``. Examples -------- `Click here <https://porespy.org/examples/filters/reference/size_to_seq.html>`_ to view online example. """ if im is None: solid = size == 0 else: solid = im == 0 uninvaded = size == -1 if bins is None: bins = np.unique(size) elif isinstance(bins, int): bins = np.linspace(0, size.max(), bins) vals = np.digitize(size, bins=bins, right=True) if mode.startswith('im'): vals[solid] = 0 vals[uninvaded] = -1 vals = make_contiguous(vals, mode='symmetric') if mode.startswith('dr'): vals = make_contiguous(vals, mode='symmetric') vals = vals.max() + 1 - vals vals[solid] = 0 vals[uninvaded] = -1 return vals
[docs] def size_to_satn(size, im=None, bins=None, mode='drainage'): r""" Converts size map to a saturation map Parameters ---------- size : ndarray The image containing invasion size values in each voxel. Solid should be indicated as 0's and uninvaded voxels as -1. im : ndarray, optional A binary image of the porous media, with ``True`` indicating the void space and ``False`` indicating the solid phase. If not given then it is assumed that the solid is identified as ``size == 0``. bins : array_like or int (optional) The bins to use when converting sizes to saturation. The default is to create 1 bin for each unique value in ``size``. If an **int** is supplied it is interpreted as the number of bins between 0 and the maximum value in ``size``. If an array is supplied it is used as the bins directly. mode : str Controls how the sizes are converted to saturations. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' The sizes are assumed to have been filled from largest to smallest, ignoring 0's and -1's 'imbibition' The sizes are assumed to have been filled from smallest to largest, ignoring 0's and -1's ============= ============================================================== Returns ------- satn : ndarray An ndarray the same shape as ``size`` but with size values replaced by the fraction of void space invaded at each size, according to the specified `mode`. Solid voxels and uninvaded voxels are represented by 0 and -1, respectively. Notes ----- If any ``-1`` values are present in `size` the maximum saturation will be less than 1.0 since this means that not all wetting phase was displaced. Examples -------- `Click here <https://porespy.org/examples/filters/reference/size_to_satn.html>`_ to view online example. """ if bins is None: bins = np.unique(size[size > 0]) elif isinstance(bins, int): bins = np.linspace(0, size.max(), bins) if im is None: im = ~(size == 0) void_vol = im.sum() satn = -np.ones_like(size, dtype=float) if mode.startswith('im'): for r in bins: hits = (size <= r) * (size > 0) temp = hits.sum()/void_vol satn[hits * (satn == -1)] = temp elif mode.startswith('dr'): for r in bins[-1::-1]: hits = (size >= r) * (size > 0) temp = hits.sum()/void_vol satn[hits * (satn == -1)] = temp satn *= (im > 0) return satn
[docs] def seq_to_satn(seq, im, mode='drainage'): r""" Converts a sequence map to a saturation map Parameters ---------- seq : ndarray The image containing invasion sequence values in each voxel. Residual phase should be indicated as 0's and uninvaded (i.e. trapped) voxels as -1. `im` is used to mask to remove solid voxels from computation im : ndarray A binary image of the porous media, with ``True`` indicating the void space and ``False`` indicating the solid phase. mode : str Controls how the sequences are converted to saturations. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' The saturation is assumed to increase with increasing sequence 'imbibition' The saturation is assumed to decrease with increasing sequence ============= ============================================================== Returns ------- satn : ndarray An ndarray the same shape as ``seq`` but with sequence values replaced by the fraction of void space invaded at the sequence number, accounting for the specified `mode`. Residual fluid and uninvaded voxels are represented by 0 and -1, respectively. Solie phase is also represented as 0, so `im` must be used as a mask. Notes ----- If any `-1` values are present in `seq` the maximum saturation will be less than 1.0 since this means that not all defending phase was displaced. If any `0` are present in `seq` then the minimum saturation will be greater than 0, since this means some invading phase was already present. Examples -------- `Click here <https://porespy.org/examples/filters/reference/seq_to_satn.html>`_ to view online example. """ seq = np.copy(seq).astype(int) trapped_mask = seq == -1 # Store uninvaded locations residual_mask = (seq == 0)*im seq[trapped_mask] = seq.max() + 1 # Assign trapped voxels with highest seq seq += 1 # Increment by 1 so residual voxels are 1 instead of 0 seq[~im] = 0 # Set solid voxels to 0 seq[residual_mask] = 1 # Set residual to 1 b = np.bincount(seq.flatten()) # Count number of voxels with each sequence if mode.startswith('imb'): c = 1 - np.cumsum(b[1:])/im.sum(dtype=np.int64) # Convert to saturation elif mode.startswith('drain'): c = np.cumsum(b[1:])/im.sum(dtype=np.int64) # Convert to saturation else: raise Exception('Unrecognized mode') c = np.hstack(([0], c)) # Add 0 for solid satn = c[seq] satn[~im] = 0.0 # Ensure solids are 0 satn[trapped_mask] = -1.0 # Update trapped voxels to have -1 still return satn
[docs] def pc_to_seq(pc, im, mode='drainage'): r""" Converts a capillary pressure map to a sequence map Parameters ---------- pc : ndarray A Numpy array with the value in each voxel indicating the capillary pressure at which it was invaded. In order to accommodate the possibility of both positive and negative capillary pressure values, trapped or residual wetting phase should be indicated by ``+inf`` and residual or trapped non-wetting phase by ``-inf``. im : ndarray A Numpy array with ``True`` values indicating the void space. This is used to mask the ``pc`` image so it does not matter what values are placed in the solid voxels of ``pc``. mode : str Controls how the pressures are converted to sequence. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' The pressures are assumed to have been filled from smallest to largest. Voxels with $-inf$ are treated as though they are invaded by non-wetting fluid at the start of the process, and voxels with $+inf$ are treated as though they are never invaded. 'imbibition' The pressures are assumed to have been filled from largest to smallest. Voxels with $+inf$ are treated as though they are already occupied by wetting fluid at the start of the process, and voxels with $*inf$ are treated as though they are never filled with wetting phase. ============= ============================================================== Returns ------- seq : ndarray A Numpy array the same shape as `pc`, with each voxel value indicating the sequence at which it was invaded, according to the specified `mode`. Uninvaded voxels (either trapped wp during drainage or trapped nwp during imbibion) are set to -1. Preinvaded voxels (either residual nwp during drainage or residual wp during imbition) are set to 0. Examples -------- `Click here <https://porespy.org/examples/filters/reference/pc_to_seq.html>`_ to view online example. """ pc = np.copy(pc) if mode == 'drainage': bins = np.unique(pc[im]) a = np.digitize(pc, bins=bins, right=False) + 1 # To ensure no voxels are 0 a[pc == -np.inf] = 0 # Set residual nwp to lowest sequence (i.e. 0) a[pc == np.inf] = -1 # Set trapped wp to -1 (only option) elif mode == 'imbibition': bins = np.flip(np.unique(pc[im])) a = np.digitize(pc, bins=bins, right=True) a[pc == np.inf] = 0 a[pc == -np.inf] = -1 a[~im] = 0 a = make_contiguous(a, mode='symmetric') return a
[docs] def pc_to_satn(pc, im, mode='drainage'): r""" Converts a capillary pressure map to a saturation map Parameters ---------- pc : ndarray A Numpy array with the value in each voxel indicating the capillary pressure at which it was invaded. In order to accommodate the possibility of both positive and negative capillary pressure values, uninvaded voxels should be indicated by ``+inf`` and residual phase by ``-inf``. Solid vs void phase is defined by ``im`` which is mandatory. im : ndarray A Numpy array with ``True`` values indicating the void space mode : str Controls how the pressures are converted to sequence. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' The pressures are assumed to have been filled from smallest to largest. 'imbibition' The pressures are assumed to have been filled from largest to smallest ============= ============================================================== Returns ------- satn : ndarray A Numpy array the same shape as `pc`, with each voxel value indicating the global saturation at which it was invaded, according to the specified `mode`. Voxels with `-inf` are treated as though they were invaded at the start of the simulation so are given a sequence number of 1 for both mode `drainage` and `imbibition`. Notes ----- If any ``+inf`` values are present the maximum saturation will be less than 1.0 since not all wetting phase was displaced. Examples -------- `Click here <https://porespy.org/examples/filters/reference/pc_to_satn.html>`_ to view online example. """ a = np.digitize(pc, bins=np.unique(pc)) a[~im] = 0 a[np.where(pc == np.inf)] = -1 satn = seq_to_satn(seq=a, im=im, mode=mode) return satn
[docs] def satn_to_seq(satn, im, residual=None, mode='drainage'): r""" Converts a saturaiton map to a sequence map Parameters ---------- satn : ndarray A Numpy array with the value in each voxel indicating the global saturation at the point it was invaded. -1 indicates a voxel that not invaded, and 0 indicates solid phase. im : ndarray A Numpy array with ``True`` values indicating the void space. residual : ndarray An `ndarray` with `True` values indicating the locations of any residual phase. If this is not provided, it will be assumed that no residual was present, which will result in an incorrect conversion to sequence if this assumption is not correct. mode : str Controls how the saturations are converted to sequence. The options are: ============= ============================================================== `mode` Description ============= ============================================================== 'drainage' Non-wetting fluid displaces wetting fluid 'imbibition' Wetting fluid dispalces the non-wetting fluid ============= ============================================================== Returns ------- seq : ndarray A Numpy array the same shape as `satn` with each voxel value indicating the sequence in which it was invaded, according to the specified `mode`. Solid voxels are indicated by 0 and uninvaded by -1. Examples -------- `Click here <https://porespy.org/examples/filters/reference/satn_to_seq.html>`_ to view online example. """ uninvaded = satn == -1 values = np.unique(satn) seq = np.digitize(satn, bins=values) # Set uninvaded to 0 seq[uninvaded] = -1 # Not sure this is needed? # Set solids back to 0 seq[~im] = 0 # Ensure values are contiguous while keeping -1 and 0 seq = make_contiguous(im=seq, mode='symmetric') if mode.startswith('im'): seq = (seq.max() + 1) - seq seq[~im] = 0 seq[uninvaded] = -1 if residual is not None: seq[residual] = 0 seq = make_contiguous(im=seq, mode='symmetric') return seq
[docs] def size_to_pc(im, size, f=None, **kwargs): r""" Converts size map into capillary pressure map Parameters ---------- im : ndarray A Numpy array with ``True`` values indicating the void space. size : ndarray The image containing invasion size values in each voxel. Solid should be indicated as 0's and uninvaded voxels as -1. f : function handle, optional A function to compute the capillary pressure which receives `size` as the first argument, followed by any additional `**kwargs`. If not provided then the Washburn equation is used, which requires `theta` and `sigma` to be specified as additional keyword arguments. **kwargs : Key word arguments All additional keyword arguments are passed on to `f`. Returns ------- pc : ndarray An image with each voxel containing the capillary pressure at which it was invaded. Any uninvaded voxels in `size` are set to `np.inf` which is meant to indicate that these voxels are never invaded. Notes ----- The function `f` should be of the form: .. code-block:: def func(r, a, b): pc = ... # Some equation for capillary pressure using r, a and b return pc """ if f is None: def f(r, sigma, theta, voxel_size): pc = -2*sigma*np.cos(np.deg2rad(theta))/(r*voxel_size) return pc pc = f(size, **kwargs) pc[~im] = 0.0 pc[im == -1] = np.inf return pc
# Seems to color in wrong order for imbibition # def satn_to_time(im, satn, flow_rate, voxel_size, mode='imbibition'): # f = 1 if mode == 'imbibition' else 0 # bins = np.digitize(satn.flatten()-f, bins=np.unique(satn-f), right=True) # counts = np.cumsum(np.bincount(bins)) # time = counts*(voxel_size**im.ndim)/flow_rate # time_map = np.reshape(time[bins], im.shape) # plt.imshow(time_map/im) if __name__ == "__main__": import porespy as ps import numpy as np import matplotlib.pyplot as plt im = np.ones([20, 20], dtype=bool) im[-1, :] = False im[0, :] = False pc = ps.filters.capillary_transform(im) inlets = ps.generators.faces(im.shape, inlet=1) inv = ps.simulations.invasion(im, pc, inlets=inlets) # Check drainage with trapping tmp = np.copy(inv.im_seq) tmp[tmp >= 12] = -1 # Set some to uninvaded/trapped satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='drainage') assert ((satn >= 0)*im).sum()/im.sum() == satn.max() # Now convert it back! seq = satn_to_seq(satn, im=im, mode='drainage') assert np.all(seq == tmp) # Check drainage wtih trapping and residual tmp = np.copy(inv.im_seq) tmp[tmp >= 12] = -1 # Set some to uninvaded/trapped tmp[tmp > 0] -= 1 residual = (tmp == 0)*im satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='drainage') assert ((satn >= 0)*im).sum()/im.sum() == satn.max() # Now convert it back! seq = satn_to_seq(satn, im=im, residual=residual, mode='drainage') assert np.all(seq == tmp) # Now check imbibition with trapping tmp = np.copy(inv.im_seq) tmp[tmp == 1] = -1 # Set some to uninvaded/trapped tmp[tmp > 0] = tmp.max() - tmp[tmp > 0] + 1 # Invert sequence numbers satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='imbibition') assert (1 - (satn > 0).sum()/im.sum()) == satn[satn > 0].min() # Now convert it back! seq = satn_to_seq(satn, im=im, mode='imbibition') assert np.all(seq == tmp) # Check drainage wtih trapping and residual tmp = np.copy(inv.im_seq) tmp[tmp == 1] = -1 # Set some to uninvaded/trapped tmp[tmp > 0] = tmp.max() - tmp[tmp > 0] + 1 # Invert sequence numbers tmp[(tmp > 0)*(tmp < 6)] = 0 tmp = ps.tools.make_contiguous(seq, mode='symmetric') residual = (tmp == 0)*im satn = ps.filters.seq_to_satn(seq=tmp, im=im, mode='drainage') assert ((satn >= 0)*im).sum()/im.sum() == satn.max() # Now convert it back! seq = satn_to_seq(satn, im=im, residual=residual, mode='drainage') assert np.all(seq == tmp)