Source code for openpnm.algorithms._advection_diffusion

import logging
import numpy as np
from openpnm.algorithms import ReactiveTransport
from openpnm.utils import Docorator


__all__ = ['AdvectionDiffusion']


docstr = Docorator()
logger = logging.getLogger(__name__)


@docstr.get_sections(base='AdvectionDiffusionSettings',
                     sections=['Parameters'])
@docstr.dedent
class AdvectionDiffusionSettings:
    r"""
    Parameters
    ----------
    %(ReactiveTransportSettings.parameters)s
    diffusive_conductance : str
        The name of the diffusive conductance values to be used by the
        specified ``conductance`` model to find the advective-diffusive
        conductance.
    hydraulic_conductance : str, optional
        The name of the hydraulic conductance values to be used by the
        specified ``conductance`` model to find the advective-diffusive
        conductance.
    pressure : str, optional
        The name of the pressure values calculated by the ``StokesFlow``
        algorithm.

    """
    quantity = 'pore.concentration'
    conductance = 'throat.ad_dif_conductance'
    diffusive_conductance = 'throat.diffusive_conductance'
    hydraulic_conductance = 'throat.hydraulic_conductance'
    pressure = 'pore.pressure'
    cache = False


[docs] class AdvectionDiffusion(ReactiveTransport): r""" A subclass of ReactiveTransport to simulate advection-diffusion """ def __init__(self, name='ad_?', **kwargs): super().__init__(name=name, **kwargs) self.settings._update(AdvectionDiffusionSettings()) self['pore.bc.outflow'] = np.nan
[docs] def set_outflow_BC(self, pores, mode='add'): r""" Adds outflow boundary condition to the selected pores Parameters ---------- pores : array_like The pore indices where the condition should be applied mode : str, optional Controls how the boundary conditions are applied. The default value is 'merge'. For definition of various modes, see the docstring for ``set_BC``. force : bool, optional If ``True`` then the ``'mode'`` is applied to all other bctypes as well. The default is ``False``. Notes ----- Outflow condition means that the gradient of the solved quantity does not change, i.e. is 0. """ pores = self._parse_indices(pores) # Calculating A[i,i] values to ensure the outflow condition network = self.project.network phase = self.project[self.settings['phase']] throats = network.find_neighbor_throats(pores=pores) C12 = network.conns[throats] P12 = phase[self.settings['pressure']][C12] gh = phase[self.settings['hydraulic_conductance']][throats] Q12 = -gh * np.diff(P12, axis=1).squeeze() Qp = np.zeros(self.Np) np.add.at(Qp, C12[:, 0], -Q12) np.add.at(Qp, C12[:, 1], Q12) self.set_BC(pores=pores, bcvalues=Qp[pores], bctype='outflow', mode=mode)
def _apply_BCs(self): r""" Applies Dirichlet, Neumann, and outflow BCs in order """ # Apply Dirichlet and rate BCs super()._apply_BCs() if 'pore.bc.outflow' not in self.keys(): return # Apply outflow BC diag = self.A.diagonal() ind = np.isfinite(self['pore.bc.outflow']) diag[ind] += self['pore.bc.outflow'][ind] self.A.setdiag(diag)
if __name__ == "__main__": import openpnm as op pn = op.network.Cubic(shape=[10, 10, 1]) pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders) pn.regenerate_models() air = op.phase.Air(network=pn) air.add_model_collection(op.models.collections.physics.standard) air.regenerate_models() del air.models['throat.diffusive_conductance'] del air.models['throat.hydraulic_conductance'] air['throat.diffusive_conductance'] = 1e-15 air['throat.hydraulic_conductance'] = 1e-15 flow = op.algorithms.StokesFlow(network=pn, phase=air) flow.set_value_BC(pores=pn.pores('left'), values=1) flow.set_value_BC(pores=pn.pores('right'), values=0) flow.run() ad = op.algorithms.AdvectionDiffusion(network=pn, phase=air) ad.settings['cache'] = False ad.set_value_BC(pores=pn.pores('front'), values=1) ad.set_outflow_BC(pores=pn.pores('back'), mode='overwrite', force=False) ad.run() import matplotlib.pyplot as plt fig, ax = plt.subplots(1, 1) ax.imshow(ad['pore.concentration'].reshape([10, 10]))