Source code for openpnm.algorithms._transient_reactive_transport
import logging
import numpy as np
from openpnm.algorithms import ReactiveTransport
from openpnm.utils import Docorator
from openpnm.integrators import ScipyRK45
from openpnm.algorithms._solution import SolutionContainer
__all__ = ['TransientReactiveTransport']
docstr = Docorator()
logger = logging.getLogger(__name__)
@docstr.get_sections(base='TransientReactiveTransportSettings',
sections=['Parameters', 'Other Parameters'])
@docstr.dedent
class TransientReactiveTransportSettings:
r"""
Parameters
----------
%(ReactiveTransportSettings.parameters)s
"""
pore_volume = 'pore.volume'
[docs]
class TransientReactiveTransport(ReactiveTransport):
r"""
A subclass of ReactiveTransport for transient simulations.
Parameters
----------
network : Network
The Network with which this algorithm is associated.
Notes
-----
Either a Network or a Project must be specified.
"""
def __init__(self, phase, name='trans_react_?', **kwargs):
super().__init__(phase=phase, name=name, **kwargs)
self.settings._update(TransientReactiveTransportSettings())
self.settings['phase'] = phase.name
self["pore.ic"] = np.nan
[docs]
def run(self, x0, tspan, saveat=None, integrator=None):
"""
Runs the transient algorithm and returns the solution.
Parameters
----------
x0 : ndarray or float
Array (or scalar) containing initial condition values.
tspan : array_like
Tuple (or array) containing the integration time span.
saveat : array_like or float, optional
If an array is passed, it signifies the time points at which
the solution is to be stored, and if a scalar is passed, it
refers to the interval at which the solution is to be stored.
integrator : Integrator, optional
Integrator object which will be used to to the time stepping.
Can be instantiated using openpnm.integrators module.
Returns
-------
TransientSolution
The solution object, which is basically a numpy array with
the added functionality that it can be called to return the
solution at intermediate times (i.e., those not stored in the
solution object).
"""
logger.info('Running TransientTransport')
if np.isscalar(saveat):
saveat = np.arange(*tspan, saveat)
# FIXME: why do we forcibly add tspan[1] to saveat even if the user
# didn't want to?
if (saveat is not None) and (tspan[1] not in saveat):
saveat = np.hstack((saveat, [tspan[1]]))
integrator = ScipyRK45() if integrator is None else integrator
# Perform pre-solve validations
self._validate_settings()
self._validate_topology_health()
self._validate_linear_system()
# Write x0 to algorithm the obj (needed by _update_iterative_props)
self['pore.ic'] = x0 = np.ones(self.Np, dtype=float) * x0
self._merge_inital_and_boundary_values()
# Build RHS (dx/dt = RHS), then integrate the system of ODEs
rhs = self._build_rhs()
# Integrate RHS using the given solver
soln = integrator.solve(rhs, x0, tspan, saveat)
# Return solution as dictionary
self.soln = SolutionContainer()
self.soln[self.settings['quantity']] = soln
def _run_special(self, x0):
pass
def _build_rhs(self):
"""
Returns a function handle, which calculates dy/dt = rhs(y, t).
Notes
-----
``y`` is the variable that the algorithms solves for, e.g., for
``TransientFickianDiffusion``, it would be concentration.
"""
def ode_func(t, y):
# TODO: add a cache mechanism
self.x = y
self._update_A_and_b()
A = self.A.tocsc()
b = self.b
V = self.network[self.settings["pore_volume"]]
return (-A.dot(y) + b) / V # much faster than A*y
return ode_func
def _merge_inital_and_boundary_values(self):
x0 = self['pore.ic']
bc_pores = ~np.isnan(self['pore.bc.value'])
x0[bc_pores] = self['pore.bc.value'][bc_pores]
quantity = self.settings['quantity']
self[quantity] = x0
