<s> """<STR_LIT>""" <EOL> import scipy as sp <EOL> import numpy as np <EOL> import matplotlib . pyplot as plt <EOL> from OpenPNM . Algorithms import GenericAlgorithm <EOL> from OpenPNM . Base import logging <EOL> logger = logging . getLogger ( __name__ ) <EOL> class Drainage ( GenericAlgorithm ) : <EOL> r"""<STR_LIT>""" <EOL> def __init__ ( self , network , name = None ) : <EOL> super ( ) . __init__ ( network = network , name = name ) <EOL> def setup ( self , <EOL> invading_phase , <EOL> defending_phase , <EOL> entry_pressure = '<STR_LIT>' , <EOL> trapping = False , <EOL> pore_filling = None , <EOL> throat_filling = None , <EOL> pore_volume = '<STR_LIT>' , <EOL> throat_volume = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> self [ '<STR_LIT>' ] = invading_phase [ entry_pressure ] <EOL> self [ '<STR_LIT>' ] = sp . inf <EOL> self [ '<STR_LIT>' ] = sp . inf <EOL> self [ '<STR_LIT>' ] = sp . inf <EOL> self [ '<STR_LIT>' ] = sp . inf <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> self . _inv_phase = invading_phase <EOL> self . _def_phase = defending_phase <EOL> self . _trapping = trapping <EOL> self . _pore_filling = pore_filling <EOL> self . _throat_filling = throat_filling <EOL> self . _throat_volume = '<STR_LIT>' <EOL> self . _pore_volume = '<STR_LIT>' <EOL> def set_inlets ( self , pores = None , mode = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> Ps = self . _parse_locations ( pores ) <EOL> if mode in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> self [ '<STR_LIT>' ] = False <EOL> if sum ( self [ '<STR_LIT>' ] [ Ps ] ) > <NUM_LIT:0> : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> bool_val = True <EOL> if mode is '<STR_LIT>' : <EOL> bool_val = False <EOL> self [ '<STR_LIT>' ] [ Ps ] = bool_val <EOL> def set_outlets ( self , pores = None , mode = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> if self . _trapping is False : <EOL> raise Exception ( '<STR_LIT>' + <EOL> '<STR_LIT>' ) <EOL> Ps = self . _parse_locations ( pores ) <EOL> if mode in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> self [ '<STR_LIT>' ] = False <EOL> if sum ( self [ '<STR_LIT>' ] [ Ps ] ) > <NUM_LIT:0> : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> bool_val = True <EOL> if mode is '<STR_LIT>' : <EOL> bool_val = False <EOL> self [ '<STR_LIT>' ] [ Ps ] = bool_val <EOL> def set_residual ( self , pores = None , throats = None , mode = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> bool_val = True <EOL> if mode is '<STR_LIT>' : <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> return <EOL> if pores is not None : <EOL> Ps = self . _parse_locations ( pores ) <EOL> if mode in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> self [ '<STR_LIT>' ] = False <EOL> if mode is '<STR_LIT>' : <EOL> bool_val = False <EOL> self [ '<STR_LIT>' ] [ Ps ] = bool_val <EOL> if throats is not None : <EOL> Ts = self . _parse_locations ( throats ) <EOL> if mode in [ '<STR_LIT>' , '<STR_LIT>' ] : <EOL> self [ '<STR_LIT>' ] = False <EOL> if mode is '<STR_LIT>' : <EOL> bool_val = False <EOL> self [ '<STR_LIT>' ] [ Ts ] = bool_val <EOL> def set_boundary_conditions ( self , bc_type = None , pores = None , throats = None , <EOL> mode = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> if bc_type is None : <EOL> if mode == '<STR_LIT>' : <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> self [ '<STR_LIT>' ] = False <EOL> return <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> if bc_type == '<STR_LIT>' : <EOL> self . set_residual ( pores = pores , throats = throats , mode = mode ) <EOL> elif bc_type == '<STR_LIT>' : <EOL> self . set_inlets ( pores = pores , mode = mode ) <EOL> elif bc_type == '<STR_LIT>' : <EOL> self . set_outlets ( pores = pores , mode = mode ) <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> def run ( self , npts = <NUM_LIT> , inv_pressures = None ) : <EOL> r"""<STR_LIT>""" <EOL> if inv_pressures is None : <EOL> logger . info ( '<STR_LIT>' ) <EOL> min_p = sp . amin ( self [ '<STR_LIT>' ] ) * <NUM_LIT> <EOL> max_p = sp . amax ( self [ '<STR_LIT>' ] ) * <NUM_LIT> <EOL> inv_points = sp . logspace ( sp . log10 ( min_p ) , <EOL> sp . log10 ( max_p ) , <EOL> npts ) <EOL> else : <EOL> inv_points = sp . unique ( inv_pressures ) <EOL> self . _inv_points = inv_points <EOL> if sp . sum ( self [ '<STR_LIT>' ] ) == <NUM_LIT:0> : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> if self . _trapping : <EOL> if sp . sum ( self [ '<STR_LIT>' ] ) == <NUM_LIT:0> : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> for inv_val in self . _inv_points : <EOL> logger . info ( '<STR_LIT>' + str ( inv_val ) ) <EOL> self . _apply_percolation ( inv_val ) <EOL> if self . _trapping : <EOL> logger . info ( '<STR_LIT>' ) <EOL> self . _check_trapping ( inv_val ) <EOL> Pinv = self [ '<STR_LIT>' ] <EOL> self [ '<STR_LIT>' ] = sp . searchsorted ( sp . unique ( Pinv ) , Pinv ) <EOL> Tinv = self [ '<STR_LIT>' ] <EOL> self [ '<STR_LIT>' ] = sp . searchsorted ( sp . unique ( Tinv ) , Tinv ) <EOL> def _check_trapping ( self , inv_val ) : <EOL> r"""<STR_LIT>""" <EOL> Tinvaded = self [ '<STR_LIT>' ] < sp . inf <EOL> Tinvaded = Tinvaded + self [ '<STR_LIT>' ] <EOL> Tdefended = ~ Tinvaded <EOL> [ pclusters , tclusters ] = self . _net . find_clusters2 ( mask = Tdefended , <EOL> t_labels = True ) <EOL> outlets = self [ '<STR_LIT>' ] * ( self [ '<STR_LIT>' ] == sp . inf ) <EOL> def_clusters = sp . unique ( pclusters [ outlets ] ) <EOL> temp = sp . in1d ( sp . unique ( pclusters ) , def_clusters , invert = True ) <EOL> trapped_clusters = sp . unique ( pclusters ) [ temp ] <EOL> trapped_clusters = trapped_clusters [ trapped_clusters >= <NUM_LIT:0> ] <EOL> pmask = np . in1d ( pclusters , trapped_clusters ) <EOL> pinds = ( self [ '<STR_LIT>' ] == sp . inf ) * ( pmask ) <EOL> self [ '<STR_LIT>' ] [ pinds ] = inv_val <EOL> tinds = self . _net . find_neighbor_throats ( pores = pinds , <EOL> mode = '<STR_LIT>' ) <EOL> self [ '<STR_LIT>' ] [ tinds ] = inv_val <EOL> self [ '<STR_LIT>' ] [ tinds ] = <NUM_LIT> <EOL> def _apply_percolation ( self , inv_val ) : <EOL> r"""<STR_LIT>""" <EOL> Tinvaded = self [ '<STR_LIT>' ] <= inv_val <EOL> Tinvaded = Tinvaded + self [ '<STR_LIT>' ] <EOL> [ pclusters , tclusters ] = self . _net . find_clusters2 ( mask = Tinvaded , <EOL> t_labels = True ) <EOL> inv_clusters = sp . unique ( pclusters [ self [ '<STR_LIT>' ] ] ) <EOL> inv_clusters = inv_clusters [ inv_clusters >= <NUM_LIT:0> ] <EOL> pmask = np . in1d ( pclusters , inv_clusters ) <EOL> pinds = ( self [ '<STR_LIT>' ] == sp . inf ) * ( pmask ) <EOL> self [ '<STR_LIT>' ] [ pinds ] = inv_val <EOL> tmask = np . in1d ( tclusters , inv_clusters ) <EOL> tinds = ( self [ '<STR_LIT>' ] == sp . inf ) * ( tmask ) <EOL> self [ '<STR_LIT>' ] [ tinds ] = inv_val <EOL> if sp . any ( self [ '<STR_LIT>' ] ) : <EOL> self [ '<STR_LIT>' ] [ self [ '<STR_LIT>' ] ] = <NUM_LIT:0> <EOL> if sp . any ( self [ '<STR_LIT>' ] ) : <EOL> self [ '<STR_LIT>' ] [ self [ '<STR_LIT>' ] ] = <NUM_LIT:0> <EOL> def get_drainage_data ( self ) : <EOL> r"""<STR_LIT>""" <EOL> PcPoints = self . _inv_points <EOL> if PcPoints [ - <NUM_LIT:1> ] == sp . inf : <EOL> PcPoints = PcPoints [ : - <NUM_LIT:1> ] <EOL> Pvol = self . _net [ self . _pore_volume ] <EOL> Tvol = self . _net [ self . _throat_volume ] <EOL> Total_vol = sp . sum ( Pvol ) + sp . sum ( Tvol ) <EOL> Vnwp_t = [ ] <EOL> Vnwp_p = [ ] <EOL> Vnwp_all = [ ] <EOL> for Pc in PcPoints : <EOL> p_inv = self [ '<STR_LIT>' ] <= Pc <EOL> if self . _pore_filling is None : <EOL> Vp = sp . sum ( Pvol [ p_inv ] ) <EOL> else : <EOL> Vp = self . _calc_fractional_filling ( pressure = Pc , <EOL> element = '<STR_LIT>' ) <EOL> Vp = sp . sum ( Vp [ p_inv ] ) <EOL> t_inv = self [ '<STR_LIT>' ] <= Pc <EOL> if self . _throat_filling is None : <EOL> Vt = sp . sum ( Tvol [ t_inv ] ) <EOL> else : <EOL> Vt = self . _calc_fractional_filling ( pressure = Pc , <EOL> element = '<STR_LIT>' ) <EOL> Vt = sp . sum ( Vt [ t_inv ] ) <EOL> Vnwp_p . append ( Vp ) <EOL> Vnwp_t . append ( Vt ) <EOL> Vnwp_all . append ( Vp + Vt ) <EOL> Snwp_all = [ V / Total_vol for V in Vnwp_all ] <EOL> data = { } <EOL> data [ '<STR_LIT>' ] = PcPoints <EOL> data [ '<STR_LIT>' ] = Snwp_all <EOL> data [ '<STR_LIT>' ] = [ <NUM_LIT:1> - s for s in Snwp_all ] <EOL> return data <EOL> def _calc_fractional_filling ( self , element , pressure ) : <EOL> r"""<STR_LIT>""" <EOL> if element == '<STR_LIT>' : <EOL> key = self . _pore_filling <EOL> vol = self . _pore_volume <EOL> elif element == '<STR_LIT>' : <EOL> key = self . _throat_filling <EOL> vol = self . _throat_volume <EOL> else : <EOL> raise Exception ( '<STR_LIT>' ) <EOL> Snwp = sp . zeros ( ( self . _count ( element ) , ) ) <EOL> for phys in self . _inv_phase . _physics : <EOL> temp_Pc = phys . models [ key ] [ '<STR_LIT>' ] <EOL> phys . models [ key ] [ '<STR_LIT>' ] = pressure <EOL> Snwp [ phys . Pnet ] = phys . models [ key ] . run ( ) <EOL> if sp . any ( self [ element + '<STR_LIT>' ] ) : <EOL> Snwp [ self [ element + '<STR_LIT>' ] ] = <NUM_LIT:1.0> <EOL> phys . models [ key ] [ '<STR_LIT>' ] = temp_Pc <EOL> V = self . _net [ vol ] * Snwp <EOL> return V <EOL> def plot_drainage_curve ( self , <EOL> data = None , <EOL> x_values = '<STR_LIT>' , <EOL> y_values = '<STR_LIT>' ) : <EOL> r"""<STR_LIT>""" <EOL> if data is None : <EOL> data = self . get_drainage_data ( ) <EOL> xdata = data [ x_values ] <EOL> ydata = data [ y_values ] <EOL> fig = plt . figure ( ) <EOL> plt . plot ( xdata , ydata , '<STR_LIT>' ) <EOL> plt . ylabel ( y_values ) <EOL> plt . xlabel ( x_values ) <EOL> plt . grid ( True ) <EOL> if sp . amax ( xdata ) <= <NUM_LIT:1> : <EOL> plt . xlim ( xmin = <NUM_LIT:0> , xmax = <NUM_LIT:1> ) <EOL> if sp . amax ( ydata ) <= <NUM_LIT:1> : <EOL> plt . ylim ( ymin = <NUM_LIT:0> , ymax = <NUM_LIT:1> ) <EOL> return fig <EOL> def return_results ( self , Pc ) : <EOL> r"""<STR_LIT>""" <EOL> Psatn = self [ '<STR_LIT>' ] <= Pc <EOL> Tsatn = self [ '<STR_LIT>' ] <= Pc <EOL> self . _inv_phase [ '<STR_LIT>' ] = sp . array ( Psatn , dtype = float ) <EOL> self . _inv_phase [ '<STR_LIT>' ] = sp . array ( Tsatn , dtype = float ) <EOL> self . _def_phase [ '<STR_LIT>' ] = sp . array ( ~ Psatn , dtype = float ) <EOL> self . _def_phase [ '<STR_LIT>' ] = sp . array ( ~ Tsatn , dtype = float ) <EOL> self . _inv_phase [ '<STR_LIT>' ] = self [ '<STR_LIT>' ] <EOL> self . _inv_phase [ '<STR_LIT>' ] = self [ '<STR_LIT>' ] <EOL> self . _inv_phase [ '<STR_LIT>' ] = self [ '<STR_LIT>' ] <EOL> self . _inv_phase [ '<STR_LIT>' ] = self [ '<STR_LIT>' ] <EOL> if self . _pore_filling : <EOL> Vp = self . _calc_fractional_filling ( element = '<STR_LIT>' , pressure = Pc ) <EOL> Sp = Vp / self . _net [ self . _pore_volume ] <EOL> zero_ps = self . _net [ self . _pore_volume ] == <NUM_LIT:0.0> <EOL> Sp [ zero_ps ] = self . _inv_phase [ '<STR_LIT>' ] [ zero_ps ] <EOL> self . _inv_phase [ '<STR_LIT>' ] = Sp <EOL> self . _def_phase [ '<STR_LIT>' ] = <NUM_LIT:1> - Sp <EOL> if self . _throat_filling : <EOL> Vt = self . _calc_fractional_filling ( element = '<STR_LIT>' , pressure = Pc ) <EOL> St = Vt / self . _net [ self . _throat_volume ] <EOL> zero_ts = self . _net [ self . _throat_volume ] == <NUM_LIT:0.0> <EOL> St [ zero_ts ] = self . _inv_phase [ '<STR_LIT>' ] [ zero_ts ] <EOL> self . _inv_phase [ '<STR_LIT>' ] = St <EOL> self . _def_phase [ '<STR_LIT>' ] = <NUM_LIT:1> - St </s>
