import pnmlib as pnm
import numpy as np
import matplotlib.pyplot as plt
from pnmlib.inspect import tree/Users/jeffgostick/Library/CloudStorage/Dropbox/Flash Sync/Code/Git/pnmlib/.venv/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdmGenerate a network:
np.random.seed(0)
Nx, Ny, Nz = 6, 4, 1
Lx = 3.0
pn = pnm.generators.cubic(
shape=[Nx, Ny, Nz],
spacing=Lx,
node_prefix='pore',
edge_prefix='throat',
)
pn['pore.left'] = pn['pore.coords'][:, 0] < (Nx/2)*Lx
pn['pore.right'] = pn['pore.coords'][:, 0] > (Nx/2)*Lx
pn['throat.spacing'] = np.ones(pn['throat.conns'].shape[0])*LxCreate an empty dict and add the network:
sim = {}
pnm.core.add_network(sim, pn)
tree(sim)├── 'pore.all': 24/24, bool ├── 'pore.coords': (24, 3), float64 ├── 'pore.left': 12/24, bool ├── 'pore.right': 12/24, bool ├── 'throat.all': 38/38, bool ├── 'throat.conns': (38, 2), int64 └── 'throat.spacing': (38,), float64
Add some parameters and properties:
pnm.core.set_data(sim, 'param.T', 298.0)
pnm.core.set_data(sim, 'param.P', 101325)
pnm.core.set_data(sim, 'throat.diameter', 1.0)
pnm.core.set_data(sim, 'phase1/pore.temperature', 1.0)
tree(sim)├── 'param.P': 101325
├── 'param.T': 298.0
├── 'pore.all': 24/24, bool
├── 'pore.coords': (24, 3), float64
├── 'pore.left': 12/24, bool
├── 'pore.right': 12/24, bool
├── 'throat.all': 38/38, bool
├── 'throat.conns': (38, 2), int64
├── 'throat.diameter': (38,), float64
├── 'throat.spacing': (38,), float64
└── phase1
└── 'pore.temperature': (24,), float64
Now let’s define some model dicts:
onetime_models = {
'pore.seed': {
'model': pnm.models.geometry.random_seeds,
'num_range': [0.5, 1.0],
},
'network/pore.volume@left': {
'model': pnm.models.geometry.constant,
'value': 3.3,
},
'network/pore.volume@right': {
'model': pnm.models.geometry.constant,
'value': 1.1,
},
}
geo_models1 = {
'pore.diameter': {
'model': pnm.models.geometry.product,
'props': [3, 'pore.seed'],
},
}
geo_models2 = {
'pore.seed2': {
'model': pnm.models.geometry.random_seeds,
'num_range': [0.2, 0.8],
},
}Now let’s run apply the models. The apply_models function runs the models using the associated parameters, and sends the results to the sim dictionary:
pnm.models.apply_models(sim, onetime_models)
pnm.models.apply_models(sim, geo_models1)
pnm.models.apply_models(sim, geo_models2, group='network')pnm.inspect.tree(sim)├── 'param.P': 101325
├── 'param.T': 298.0
├── 'pore.all': 24/24, bool
├── 'pore.coords': (24, 3), float64
├── 'pore.diameter': (24,), float64
├── 'pore.left': 12/24, bool
├── 'pore.right': 12/24, bool
├── 'pore.seed': (24,), float64
├── 'throat.all': 38/38, bool
├── 'throat.conns': (38, 2), int64
├── 'throat.diameter': (38,), float64
├── 'throat.spacing': (38,), float64
├── network
│ ├── 'pore.seed2': (24,), float64
│ └── 'pore.volume': (24,), float64
└── phase1
└── 'pore.temperature': (24,), float64
phase_models = {
'phase1/pore.viscosity': {
'model': pnm.models.geometry.constant,
'value': 0.001,
},
'phase1/throat.viscosity': {
'model': pnm.models.geometry.constant,
'value': 0.001,
},
}
pnm.models.apply_models(sim, phase_models)
phys_models = {
'phase1/throat.hydraulic_conductance': {
'model': pnm.models.physics.hydraulic_conductance2,
'pore_viscosity': 'phase1/pore.viscosity',
'throat_viscosity': 'phase1/throat.viscosity',
'pore_diameter': 'pore.diameter',
'throat_diameter': 'throat.diameter',
},
}
pnm.models.apply_models(sim, phys_models)phase = pnm.core.get_data(sim, 'phase1/*')
pnm.inspect.info(phase)
══════════════════════════════════════════════════════════════════════════════
# Properties Valid Values
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
2 phase1/pore.temperature 24 / 24
3 phase1/pore.viscosity 24 / 24
4 phase1/throat.hydraulic_conductance 38 / 38
5 phase1/throat.viscosity 38 / 38
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
# Labels Assigned Locations
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――
――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――pnm.simulations.init_coefficient_matrix(
project=sim,
conductance_values='phase1/throat.hydraulic_conductance',
name='stokes_flow1',
)
pnm.simulations.init_rhs(
project=sim,
name='stokes_flow1',
)
pnm.simulations.set_value_bc(
project=sim,
values=1.5,
locs=[0],
name='stokes_flow1',
)
pnm.simulations.set_value_bc(
project=sim,
values=2.0,
locs=[20],
name='stokes_flow1',
)
pnm.simulations.set_rate_bc(
project=sim,
rates=-2.0,
locs=[12],
name='stokes_flow1',
)A = pnm.simulations.get_A(project=sim, name='stokes_flow1')
b = pnm.simulations.get_rhs(project=sim, name='stokes_flow1')
x = pnm.simulations.solve_sp(A, b)
fig, ax = plt.subplots(1, 2)
ax[0].imshow(x.reshape([Nx, Ny]).T, origin='lower', vmin=0, vmax=2)
ax[0].set_title("No reaction")
func = 'a*(x**b)'
s = pnm.simulations.get_source_term(func, a=-500, b=2)
s['locs'] = [12]
A, b, x = pnm.simulations.solve_reactive(
A, b, rxns=s, f=[0, 0], maxiter=100, rtol=1e-12, atol=1e-12, solver='sp')
ax[1].imshow(x.reshape([Nx, Ny]).T, origin='lower', vmin=0, vmax=2)
ax[1].set_title("With reaction");Solution converged after 6 iterations
Current residual is 3.5353936830706666e-13
pnm.core.set_data(sim, 'pore.concentration', x)
pnm.core.set_data(sim, 'throat.concentration', 1)
fig, ax = plt.subplots(figsize=(Nx, Ny))
pnm.visualization.draw_network(sim, cmap=plt.cm.viridis, color_by='concentration', ax=ax);