Repository: PMEAL/OpenPNM

Issue:
Outflow BC and other BC types must be prevented to coexist in a given set of pores

Discussion:

Patch:
diff --git a/openpnm/algorithms/AdvectionDiffusion.py b/openpnm/algorithms/AdvectionDiffusion.py
index 7554b0429..5373467f1 100644
--- a/openpnm/algorithms/AdvectionDiffusion.py
+++ b/openpnm/algorithms/AdvectionDiffusion.py
@@ -74,7 +74,7 @@ class AdvectionDiffusion(ReactiveTransport):
             **kwargs
     ):
         r"""
-        Setup method for setting/modifying algorithm settings.
+        Setup method for setting/modifying algorithm settings
         """
         if phase:
             self.settings['phase'] = phase.name
@@ -92,17 +92,31 @@ class AdvectionDiffusion(ReactiveTransport):
 
     def set_outflow_BC(self, pores, mode='merge'):
         r"""
-        Adds outflow boundary condition to the selected pores.
+        Adds outflow boundary condition to the selected pores
+
+        Parameters
+        ----------
+        pores : array_like
+            The pore indices where the condition should be applied
+        mode : string, optional
+            Controls how the boundary conditions are applied.  Options are:
+
+            'merge' - (Default) Adds supplied boundary conditions to already
+            existing conditions, and also overwrites any existing values.
+            If at rate or value BC exists at the given locations, these
+            are deleted, and outflow conditions are given priority.
+            'overwrite' - Deletes all boundary conditions of the given type
+            then adds the specified new ones.
 
         Notes
         -----
-        Outflow condition simply means that the gradient of the solved
-        quantity does not change, i.e. is 0.
+        Outflow condition means that the gradient of the solved quantity
+        does not change, i.e. is 0.
 
         """
         # Hijack the parse_mode function to verify mode/pores argument
-        allowed_modes = ['merge', 'overwrite', 'remove']
-        mode = self._parse_mode(mode, allowed=allowed_modes, single=True)
+        mode = self._parse_mode(mode, allowed=['merge', 'overwrite'],
+                                single=True)
         pores = self._parse_indices(pores)
 
         # Calculating A[i,i] values to ensure the outflow condition
@@ -117,11 +131,32 @@ class AdvectionDiffusion(ReactiveTransport):
         np.add.at(Qp, C12[:, 0], -Q12)
         np.add.at(Qp, C12[:, 1], Q12)
 
+        # Ensure other BCs are not already applied at given pores
+        hits = ~np.isnan(self['pore.bc_rate'][pores])
+        if np.any(hits):
+            self['pore.bc_rate'][pores] = np.nan
+            logger.warning('Rate boundary conditions found in some of the '
+                           + 'specified pores will be overwritten')
+        hits = ~np.isnan(self['pore.bc_value'][pores])
+        if np.any(hits):
+            self['pore.bc_value'][pores] = np.nan
+            logger.warning('Value boundary conditions found in some of the '
+                           + 'specified pores will be overwritten')
         # Store boundary values
         if ('pore.bc_outflow' not in self.keys()) or (mode == 'overwrite'):
             self['pore.bc_outflow'] = np.nan
         self['pore.bc_outflow'][pores] = Qp[pores]
 
+    def remove_BC(self, pores=None, bctype='all'):
+        # parse bctype argument
+        if isinstance(bctype, str):
+            bctype = [bctype]
+        if 'all' in bctype:
+            bctype = ['value', 'rate', 'outflow']
+        if ('pore.bc_outflow' in self.keys()) and ('outflow' in bctype):
+            self['pore.bc_outflow'][pores] = np.nan
+        super().remove_BC(pores=pores, bctype=bctype)
+
     def _apply_BCs(self):
         r"""
         Applies Dirichlet, Neumann, and outflow BCs in order
@@ -135,3 +170,15 @@ class AdvectionDiffusion(ReactiveTransport):
         ind = np.isfinite(self['pore.bc_outflow'])
         diag[ind] += self['pore.bc_outflow'][ind]
         self.A.setdiag(diag)
+
+    def _set_BC(self, pores, bctype, bcvalues=None, mode='merge'):
+        pores = self._parse_indices(pores)
+        # First check that given pores outflow BCs already applied
+        if 'pore.bc_outflow' in self.keys():
+            hits = ~np.isnan(self['pore.bc_outflow'][pores])
+            if np.any(hits):
+                raise Exception('Cannot apply BCs to the following pores '
+                                + 'which already have an outflow BC '
+                                + 'specified', pores[np.where(hits)])
+        # Then call parent class function if above check passes
+        super()._set_BC(pores=pores, bctype=bctype, bcvalues=bcvalues, mode=mode)
diff --git a/openpnm/algorithms/GenericTransport.py b/openpnm/algorithms/GenericTransport.py
index ef8440c11..dab16482e 100644
--- a/openpnm/algorithms/GenericTransport.py
+++ b/openpnm/algorithms/GenericTransport.py
@@ -266,13 +266,13 @@ class GenericTransport(GenericAlgorithm):
         mode : string, optional
             Controls how the boundary conditions are applied.  Options are:
 
-            +-------------+--------------------------------------------------+
-            | 'merge'     | (Default) Adds supplied boundary conditions to   |
-            |             | already existing conditions                      |
-            +-------------+--------------------------------------------------+
-            | 'overwrite' | Deletes all boundary condition on object then    |
-            |             | adds the given ones                              |
-            +-------------+--------------------------------------------------+
+            'merge' - (Default) Adds supplied boundary conditions to already
+            existing conditions, and also overwrites any existing values.
+            If BCs of the complementary type already exist in the given
+            locations, those values are kept.
+            'overwrite' - Deletes all boundary conditions of the given type
+            then adds the specified new ones (unless locations already have
+            BCs of the other type).
 
         Notes
         -----
@@ -280,7 +280,6 @@ class GenericTransport(GenericAlgorithm):
         ``settings``, e.g. ``alg.settings['quantity'] = 'pore.pressure'``.
 
         """
-        mode = self._parse_mode(mode, allowed=['merge', 'overwrite'], single=True)
         self._set_BC(pores=pores, bctype='value', bcvalues=values, mode=mode)
 
     def set_rate_BC(self, pores, rates=None, total_rate=None, mode='merge',
@@ -303,13 +302,13 @@ class GenericTransport(GenericAlgorithm):
         mode : str, optional
             Controls how the boundary conditions are applied.  Options are:
 
-            +-------------+--------------------------------------------------+
-            | 'merge'     | (Default) Adds supplied boundary conditions to   |
-            |             | already existing conditions                      |
-            +-------------+--------------------------------------------------+
-            | 'overwrite' | Deletes all boundary condition on object then    |
-            |             | adds the given ones                              |
-            +-------------+--------------------------------------------------+
+            'merge' - (Default) Adds supplied boundary conditions to already
+            existing conditions, and also overwrites any existing values.
+            If BCs of the complementary type already exist in the given
+            locations, these values are kept.
+            'overwrite' - Deletes all boundary conditions of the given type
+            then adds the specified new ones (unless locations already have
+            BCs of the other type).
 
         Notes
         -----
@@ -331,7 +330,6 @@ class GenericTransport(GenericAlgorithm):
                                 + 'total_rate')
             pores = self._parse_indices(pores)
             rates = total_rate/pores.size
-        mode = self._parse_mode(mode, allowed=['merge', 'overwrite'], single=True)
         self._set_BC(pores=pores, bctype='rate', bcvalues=rates, mode=mode)
 
     @docstr.get_sectionsf(base='GenericTransport._set_BC',
@@ -349,12 +347,8 @@ class GenericTransport(GenericAlgorithm):
             Specifies the type or the name of boundary condition to apply. The
             types can be one one of the following:
 
-            +-------------+--------------------------------------------------+
-            | 'value'     | Specify the value of the quantity in each        |
-            |             | location                                         |
-            +-------------+--------------------------------------------------+
-            | 'rate'      | Specify the flow rate into each location         |
-            +-------------+--------------------------------------------------+
+            'value' - Specify the value of the quantity in each location
+            'rate' - Specify the flow rate into each location
 
         bcvalues : int or array_like
             The boundary value to apply, such as concentration or rate.  If
@@ -366,13 +360,13 @@ class GenericTransport(GenericAlgorithm):
         mode : string, optional
             Controls how the boundary conditions are applied.  Options are:
 
-            +-------------+--------------------------------------------------+
-            | 'merge'     | (Default) Adds supplied boundary conditions to   |
-            |             | already existing conditions                      |
-            +-------------+--------------------------------------------------+
-            | 'overwrite' | Deletes all boundary condition on object then    |
-            |             | adds the given ones                              |
-            +-------------+--------------------------------------------------+
+            'merge' - (Default) Adds supplied boundary conditions to already
+            existing conditions, and also overwrites any existing values.
+            If BCs of the complementary type already exist in the given
+            locations, these values are kept.
+            'overwrite' - Deletes all boundary conditions of the given type
+            then adds the specified new ones (unless locations already have
+            BCs of the other type).
 
         Notes
         -----
@@ -385,6 +379,7 @@ class GenericTransport(GenericAlgorithm):
         # Hijack the parse_mode function to verify bctype argument
         bctype = self._parse_mode(bctype, allowed=['value', 'rate'],
                                   single=True)
+        othertype = list(set(['value', 'rate']).difference(set([bctype])))[0]
         mode = self._parse_mode(mode, allowed=['merge', 'overwrite'],
                                 single=True)
         pores = self._parse_indices(pores)
@@ -394,16 +389,25 @@ class GenericTransport(GenericAlgorithm):
             raise Exception('The number of boundary values must match the '
                             + 'number of locations')
 
-        # Warn the user that another boundary condition already exists
-        value_BC_mask = np.isfinite(self["pore.bc_value"])
-        rate_BC_mask = np.isfinite(self["pore.bc_rate"])
-        BC_locs = self.Ps[rate_BC_mask + value_BC_mask]
-        if np.intersect1d(pores, BC_locs).size:
-            logger.info('Another boundary condition detected in some locations!')
+        # Create boundary array if needed (though these are created on init)
+        if 'pore.bc_' + bctype not in self.keys():
+            self['pore.bc_' + bctype] = np.nan
 
-        # Clear old boundary values if needed
-        if ('pore.bc_' + bctype not in self.keys()) or (mode == 'overwrite'):
+        # Catch pores with existing BCs
+        if mode == 'merge':  # remove offenders, and warn user
+            existing_bcs = np.isfinite(self["pore.bc_" + othertype])
+            inds = pores[existing_bcs[pores]]
+        elif mode == 'overwrite':  # Remove existing BCs and write new ones
             self['pore.bc_' + bctype] = np.nan
+            existing_bcs = np.isfinite(self["pore.bc_" + othertype])
+            inds = pores[existing_bcs[pores]]
+        # Now drop any pore indices which have BCs that should be kept
+        if len(inds) > 0:
+            msg = r'Boundary conditions are already specified in ' + \
+                  r'the following given pores, so these will be skipped: '
+            msg = '\n'.join((msg, inds.__repr__()))
+            logger.warning(msg)
+            pores = np.array(list(set(pores).difference(set(inds))), dtype=int)
 
         # Store boundary values
         self['pore.bc_' + bctype][pores] = values
@@ -424,21 +428,18 @@ class GenericTransport(GenericAlgorithm):
             -*'all'*: (default) Removes all boundary conditions
             -*'value'*: Removes only value conditions
             -*'rate'*: Removes only rate conditions
-            -*'outflow'*: Removes only outflow conditions
 
         """
         if isinstance(bctype, str):
             bctype = [bctype]
         if 'all' in bctype:
-            bctype = ['value', 'rate', 'outflow']
+            bctype = ['value', 'rate']
         if pores is None:
             pores = self.Ps
         if ('pore.bc_value' in self.keys()) and ('value' in bctype):
             self['pore.bc_value'][pores] = np.nan
         if ('pore.bc_rate' in self.keys()) and ('rate' in bctype):
             self['pore.bc_rate'][pores] = np.nan
-        if ('pore.bc_outflow' in self.keys()) and ('outflow' in bctype):
-            self['pore.bc_outflow'][pores] = np.nan
 
     def _build_A(self):
         r"""


Test patch:
diff --git a/tests/unit/algorithms/AdvectionDiffusionTest.py b/tests/unit/algorithms/AdvectionDiffusionTest.py
index e4dc21e0c..46404808c 100644
--- a/tests/unit/algorithms/AdvectionDiffusionTest.py
+++ b/tests/unit/algorithms/AdvectionDiffusionTest.py
@@ -2,6 +2,7 @@ import numpy as np
 import openpnm as op
 from numpy.testing import assert_allclose
 from openpnm.algorithms import AdvectionDiffusion, StokesFlow
+import pytest
 
 
 class AdvectionDiffusionTest:
@@ -156,6 +157,30 @@ class AdvectionDiffusionTest:
             y = ad[ad.settings['quantity']].mean()
             assert_allclose(actual=y, desired=2, rtol=1e-5)
 
+    def test_add_outflow_overwrite_rate_and_value_BC(self):
+        ad = AdvectionDiffusion(network=self.net, phase=self.phase)
+        ad.set_rate_BC(pores=[0, 1], total_rate=1)
+        ad.set_value_BC(pores=[2, 3], values=1)
+        assert np.sum(np.isfinite(ad['pore.bc_rate'])) == 2
+        assert np.sum(np.isfinite(ad['pore.bc_value'])) == 2
+        ad.set_outflow_BC(pores=[0, 1, 2, 3])
+        assert np.sum(np.isfinite(ad['pore.bc_rate'])) == 0
+        assert np.sum(np.isfinite(ad['pore.bc_value'])) == 0
+
+    def test_value_BC_does_not_overwrite_outflow(self):
+        ad = AdvectionDiffusion(network=self.net, phase=self.phase)
+        ad.set_outflow_BC(pores=[0, 1])
+        with pytest.raises(Exception):
+            ad.set_value_BC(pores=[0, 1], values=1)
+
+    def test_add_rate_BC_fails_when_outflow_BC_present(self):
+        ad = AdvectionDiffusion(network=self.net, phase=self.phase)
+        ad.set_outflow_BC(pores=[0, 1])
+        with pytest.raises(Exception):
+            ad.set_rate_BC(pores=[0, 1], total_rate=1)
+        ad.set_rate_BC(pores=[2, 3], total_rate=1)
+        assert np.all(ad['pore.bc_rate'][[2, 3]] == 0.5)
+
     def test_outflow_BC_rigorous(self):
         ad = AdvectionDiffusion(network=self.net, phase=self.phase)
         ad.settings["cache_A"] = False
@@ -220,8 +245,8 @@ class AdvectionDiffusionTest:
 if __name__ == '__main__':
     t = AdvectionDiffusionTest()
     t.setup_class()
+    self = t
     for item in t.__dir__():
         if item.startswith('test'):
             print(f'Running test: {item}')
             t.__getattribute__(item)()
-    self = t
diff --git a/tests/unit/algorithms/GenericTransportTest.py b/tests/unit/algorithms/GenericTransportTest.py
index f8cf4858d..aa0b83ed4 100644
--- a/tests/unit/algorithms/GenericTransportTest.py
+++ b/tests/unit/algorithms/GenericTransportTest.py
@@ -98,6 +98,27 @@ class GenericTransportTest:
         y = np.unique(np.around(alg['pore.mole_fraction'], decimals=3))
         assert np.all(x == y)
 
+    def test_set_value_bc_where_rate_is_already_set_mode_merge(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_rate_BC(pores=[0, 1], values=1, mode='merge')
+        alg.set_value_BC(pores=[1, 2], values=0, mode='merge')
+        assert np.isfinite(alg['pore.bc_rate']).sum() == 2
+        assert np.isfinite(alg['pore.bc_value']).sum() == 1
+
+    def test_set_value_bc_where_rate_is_already_set_mode_overwrite(self):
+        alg = op.algorithms.GenericTransport(network=self.net,
+                                             phase=self.phase)
+        alg.settings['conductance'] = 'throat.diffusive_conductance'
+        alg.settings['quantity'] = 'pore.mole_fraction'
+        alg.set_rate_BC(pores=[0, 1], values=1, mode='merge')
+        alg.set_value_BC(pores=[2, 3, 4], values=0, mode='merge')
+        alg.set_value_BC(pores=[1, 2], values=0, mode='overwrite')
+        assert np.isfinite(alg['pore.bc_rate']).sum() == 2
+        assert np.isfinite(alg['pore.bc_value']).sum() == 1
+
     def test_cache_A(self):
         alg = op.algorithms.GenericTransport(network=self.net,
                                              phase=self.phase)
diff --git a/tests/unit/algorithms/SolversTest.py b/tests/unit/algorithms/SolversTest.py
index 4fd53af69..3214ec661 100644
--- a/tests/unit/algorithms/SolversTest.py
+++ b/tests/unit/algorithms/SolversTest.py
@@ -145,8 +145,8 @@ class SolversTest:
 if __name__ == '__main__':
     t = SolversTest()
     t.setup_class()
+    self = t
     for item in t.__dir__():
         if item.startswith('test'):
             print('running test: '+item)
             t.__getattribute__(item)()
-    self = t