Repository: PMEAL/OpenPNM
Pull request #1188 (2019-01-04 16:19:18)

Exclusive rate method for AdvectionDiffusion algorithm
Addressing issue #1150. This PR adds an exclusive rate method for AdvectionDiffusion algorithm. This method overrides the rate method on GenericTransport since it only accounted for diffusive flux.


Patch:
diff --git a/openpnm/algorithms/AdvectionDiffusion.py b/openpnm/algorithms/AdvectionDiffusion.py
index e3a76bf9b..bcfe534cd 100644
--- a/openpnm/algorithms/AdvectionDiffusion.py
+++ b/openpnm/algorithms/AdvectionDiffusion.py
@@ -136,3 +136,92 @@ class AdvectionDiffusion(ReactiveTransport):
         ind = np.isfinite(self['pore.bc_outflow'])
         diag[ind] += self['pore.bc_outflow'][ind]
         self.A.setdiag(diag)
+
+    def rate(self, pores=[], throats=[], mode='group'):
+        r"""
+        Calculates the net rate of material moving into a given set of pores or
+        throats
+
+        Parameters
+        ----------
+        pores : array_like
+            The pores for which the rate should be calculated
+
+        throats : array_like
+            The throats through which the rate should be calculated
+
+        mode : string, optional
+            Controls how to return the rate.  Options are:
+
+            *'group'*: (default) Returns the cumulative rate of material
+            moving into the given set of pores
+
+            *'single'* : Calculates the rate for each pore individually
+
+        Returns
+        -------
+        If ``pores`` are specified, then the returned values indicate the
+        net rate of material exiting the pore or pores.  Thus a positive
+        rate indicates material is leaving the pores, and negative values
+        mean material is entering.
+
+        If ``throats`` are specified the rate is calculated in the direction of
+        the gradient, thus is always positive.
+
+        If ``mode`` is 'single' then the cumulative rate through the given
+        pores (or throats) are returned as a vector, if ``mode`` is 'group'
+        then the individual rates are summed and returned as a scalar.
+
+        """
+        pores = self._parse_indices(pores)
+        throats = self._parse_indices(throats)
+
+        network = self.project.network
+        C12 = network["throat.conns"]
+        phase = self.project.phases()[self.settings['phase']]
+
+        P = phase[self.settings['pressure']]
+        X = self[self.settings['quantity']]
+        gh = phase[self.settings['hydraulic_conductance']]
+        gd = phase[self.settings['diffusive_conductance']]
+
+        Q12 = -gh * np.diff(P[C12], axis=1).squeeze()
+        Pe12 = Q12 / gd
+        Pe12 = np.abs(Pe12).clip(min=1e-10) * np.sign(Pe12)
+        Q12 = Pe12 * gd
+
+        X12 = X[C12]
+        dX = -np.diff(X12, axis=1).squeeze()
+
+        s_dis = self.settings['s_scheme']
+
+        if s_dis == 'upwind':
+            w = gd + np.maximum(0, -Q12)
+            Qt = -Q12 * X12[:, 0] - dX * w
+        elif s_dis == 'hybrid':
+            w = np.maximum(0, np.maximum(-Q12, gd-Q12/2))
+            Qt = -Q12 * X12[:, 0] - dX * w
+        elif s_dis == 'powerlaw':
+            w = gd * np.maximum(0, (1 - 0.1*np.abs(Pe12))**5) + \
+                np.maximum(0, -Q12)
+            Qt = -Q12 * X12[:, 0] - dX * w
+        elif s_dis == 'exponential':
+            w = -Q12 / (1 - np.exp(Pe12))
+            Qt = -Q12 * X12[:, 0] - dX * w
+
+        if len(throats) and len(pores):
+            raise Exception('Must specify either pores or throats, not both')
+
+        if len(throats):
+            R = np.absolute(Qt[throats])
+
+        if len(pores):
+            Qp = np.zeros((self.Np,))
+            np.add.at(Qp, C12[:, 0], -Qt)
+            np.add.at(Qp, C12[:, 1], Qt)
+            R = Qp[pores]
+
+        if mode == 'group':
+            R = np.sum(R)
+
+        return np.array(R, ndmin=1)


Test patch:
diff --git a/tests/unit/algorithms/AdvectionDiffusionTest.py b/tests/unit/algorithms/AdvectionDiffusionTest.py
index 880f70295..abdacdaa7 100644
--- a/tests/unit/algorithms/AdvectionDiffusionTest.py
+++ b/tests/unit/algorithms/AdvectionDiffusionTest.py
@@ -72,13 +72,15 @@ class AdvectionDiffusionTest:
         assert_allclose(actual=y, desired=x)
 
     def test_outflow_BC(self):
-        for scheme in ['upwind', 'hybrid', 'powerlaw']:
+        for scheme in ['upwind', 'hybrid', 'powerlaw', 'exponential']:
             ad = op.algorithms.AdvectionDiffusion(network=self.net,
                                                   phase=self.phase)
             ad.setup(quantity='pore.concentration',
                      diffusive_conductance='throat.diffusive_conductance',
                      hydraulic_conductance='throat.hydraulic_conductance',
-                     pressure='pore.pressure')
+                     pressure='pore.pressure',
+                     s_scheme=scheme)
+
             ad.set_value_BC(pores=self.net.pores('back'), values=2)
             ad.set_outflow_BC(pores=self.net.pores('front'))
             ad.run()
@@ -86,6 +88,31 @@ class AdvectionDiffusionTest:
             y = ad[ad.settings['quantity']].mean()
             assert_allclose(actual=y, desired=2.0)
 
+    def test_rate(self):
+        for scheme in ['upwind', 'hybrid', 'powerlaw', 'exponential']:
+            ad = op.algorithms.AdvectionDiffusion(network=self.net,
+                                                  phase=self.phase)
+            ad.setup(quantity='pore.concentration',
+                     diffusive_conductance='throat.diffusive_conductance',
+                     hydraulic_conductance='throat.hydraulic_conductance',
+                     pressure='pore.pressure',
+                     s_scheme=scheme)
+
+            ad.set_value_BC(pores=self.net.pores('back'), values=2)
+            ad.set_value_BC(pores=self.net.pores('front'), values=0)
+            ad.run()
+
+            mdot_inlet = ad.rate(pores=self.net.pores("back"))[0]
+            mdot_outlet = ad.rate(pores=self.net.pores("front"))[0]
+            temp = sp.random.choice(self.net.pores(["back", "front"],
+                                                    mode="not"),
+                                    size=3, replace=False)
+            mdot_internal = ad.rate(pores=temp)[0]
+            # Ensure no mass is generated within the network
+            assert_allclose(mdot_inlet - mdot_internal, mdot_inlet)
+            # Ensure mass is not conserved
+            assert_allclose(mdot_inlet, -mdot_outlet)
+
     def teardown_class(self):
         ws = op.Workspace()
         ws.clear()