Repository: PMEAL/OpenPNM

Issue:
Implement Solution class for steady state simulations
Currently, `Solution` object is only implemented for transient simulations.

Discussion:

Patch:
diff --git a/openpnm/__init__.py b/openpnm/__init__.py
index 9f4758cfb..d5f833e2d 100644
--- a/openpnm/__init__.py
+++ b/openpnm/__init__.py
@@ -22,9 +22,9 @@ from . import geometry
 from . import phases
 from . import physics
 from . import models
+from . import algorithms
 from . import solvers
 from . import integrators
-from . import algorithms
 from . import materials
 from . import topotools
 from . import io
diff --git a/openpnm/algorithms/_generic_transport.py b/openpnm/algorithms/_generic_transport.py
index 84aef6658..4352e3ee9 100644
--- a/openpnm/algorithms/_generic_transport.py
+++ b/openpnm/algorithms/_generic_transport.py
@@ -1,16 +1,12 @@
-import numpy as np
-import scipy.sparse.linalg
 import warnings
+import numpy as np
 import scipy.sparse.csgraph as spgr
-from copy import deepcopy
-from scipy.spatial import ConvexHull
-from scipy.spatial import cKDTree
-from openpnm.topotools import iscoplanar, is_fully_connected, dimensionality
+from openpnm.topotools import is_fully_connected
 from openpnm.algorithms import GenericAlgorithm
 from openpnm.utils import logging, prettify_logger_message
 from openpnm.utils import Docorator, SettingsAttr
-from openpnm.utils import is_symmetric
 from openpnm.solvers import PardisoSpsolve
+from ._solution import SteadyStateSolution
 docstr = Docorator()
 logger = logging.getLogger(__name__)
 
@@ -433,6 +429,9 @@ class GenericTransport(GenericAlgorithm):
         in the ``settings`` attribute.
 
         """
+        # Initialize the solution object
+        self.soln = SteadyStateSolution(np.full(self.Np, fill_value=np.nan))
+        self.soln.is_converged = False
         logger.info('Running GenericTransport')
         solver = PardisoSpsolve() if solver is None else solver
         # Perform pre-solve validations
@@ -445,6 +444,7 @@ class GenericTransport(GenericAlgorithm):
         # Build A and b, then solve the system of equations
         self._update_A_and_b()
         self._run_special(solver=solver, x0=x0)
+        return self.soln
 
     def _run_special(self, solver, x0, w=1):
         # Make sure A,b are STILL well-defined
@@ -454,6 +454,9 @@ class GenericTransport(GenericAlgorithm):
         self.x = w * x_new + (1 - w) * self.x
         # Update A and b using the recent solution
         self._update_A_and_b()
+        # Update SteadyStateSolution object on algorithm
+        self.soln[:] = self.x
+        self.soln.is_converged = not bool(exit_code)
 
     def _update_A_and_b(self):
         r"""
diff --git a/openpnm/algorithms/_reactive_transport.py b/openpnm/algorithms/_reactive_transport.py
index 6a2f6cf38..13d54ff78 100644
--- a/openpnm/algorithms/_reactive_transport.py
+++ b/openpnm/algorithms/_reactive_transport.py
@@ -4,6 +4,7 @@ from scipy.optimize.nonlin import TerminationCondition
 from openpnm.algorithms import GenericTransport
 from openpnm.utils import logging
 from openpnm.utils import TypedList, Docorator, SettingsAttr
+from ._solution import SteadyStateSolution
 docstr = Docorator()
 logger = logging.getLogger(__name__)
 
@@ -238,15 +239,19 @@ class ReactiveTransport(GenericTransport):
         condition = TerminationCondition(f_rtol=f_rtol, x_rtol=x_rtol)
 
         for i in range(maxiter):
+            self.soln.num_iter = i + 1
             res = self._get_residual()
-            self.is_converged = condition.check(f=res, x=xold, dx=dx)
-            if self.is_converged:
+            is_converged = bool(condition.check(f=res, x=xold, dx=dx))
+            if is_converged:
+                self.soln.is_converged = is_converged
                 logger.info(f'Solution converged, residual norm: {norm(res):.4e}')
                 return
             super()._run_special(solver=solver, x0=xold, w=w)
             dx = self.x - xold
             xold = self.x
             logger.info(f'Iteration #{i:<4d} | Residual norm: {norm(res):.4e}')
+
+        self.soln.is_converged = False
         logger.critical(f"{self.name} didn't converge after {maxiter} iterations")
 
     def _update_A_and_b(self):
diff --git a/openpnm/algorithms/_solution.py b/openpnm/algorithms/_solution.py
new file mode 100644
index 000000000..2fa0bcf77
--- /dev/null
+++ b/openpnm/algorithms/_solution.py
@@ -0,0 +1,53 @@
+import numpy as np
+from scipy.interpolate import interp1d
+
+
+class Solution(np.ndarray):
+    """Brief description of 'Solution'"""
+    pass
+
+
+class SteadyStateSolution(np.ndarray):
+    r"""Brief description of 'SteadyStateSolution'"""
+
+    def __new__(cls, x):
+        obj = np.asarray(x).view(cls)
+        return obj
+
+
+class TransientSolution(Solution):
+    """Brief description of 'TransientSolution'"""
+
+    def __new__(cls, t, x):
+        obj = np.asarray(x).view(cls)
+        obj.t = np.asarray(t).view(cls)
+        return obj
+
+    def _create_interpolant(self):
+        self._interpolant = interp1d(self.t, self, bounds_error=True)
+
+    def interpolate(self, t):
+        """
+        Interpolates solution at time 't'.
+
+        Parameters
+        ----------
+        t : float
+            Time at which the solution is to be interpolated
+
+        Returns
+        -------
+        ndarray
+            Transient solution interpolated at the given time 't'
+
+        Notes
+        -----
+        't' must reside inside the 'tspan' used during integration.
+
+        """
+        # Cache interpolant to avoid overhead
+        if not hasattr(self, "_interpolant"):
+            self._create_interpolant()
+        return self._interpolant(t)
+
+    __call__ = interpolate
diff --git a/openpnm/integrators/_scipy.py b/openpnm/integrators/_scipy.py
index 407f5a27f..74ec21287 100644
--- a/openpnm/integrators/_scipy.py
+++ b/openpnm/integrators/_scipy.py
@@ -1,9 +1,8 @@
-import numpy as np
 from scipy.integrate import solve_ivp
-from scipy.interpolate import interp1d
 from openpnm.integrators import Integrator
+from openpnm.algorithms._solution import TransientSolution
 
-__all__ = ['ScipyRK45', 'Solution', 'TransientSolution']
+__all__ = ['ScipyRK45']
 
 
 class ScipyRK45(Integrator):
@@ -54,46 +53,3 @@ class ScipyRK45(Integrator):
         if sol.success:
             return TransientSolution(sol.t, sol.y)
         raise Exception(sol.message)
-
-
-class Solution(np.ndarray):
-    """Brief description of 'Solution'"""
-    ...
-
-
-class TransientSolution(Solution):
-    """Brief description of 'TransientSolution'"""
-
-    def __new__(cls, t, x):
-        obj = np.asarray(x).view(cls)
-        obj.t = np.asarray(t).view(cls)
-        return obj
-
-    def _create_interpolant(self):
-        self._interpolant = interp1d(self.t, self, bounds_error=True)
-
-    def interpolate(self, t):
-        """
-        Interpolates solution at time 't'.
-
-        Parameters
-        ----------
-        t : float
-            Time at which the solution is to be interpolated
-
-        Returns
-        -------
-        ndarray
-            Transient solution interpolated at the given time 't'
-
-        Notes
-        -----
-        't' must reside inside the 'tspan' used during integration.
-
-        """
-        # Cache interpolant to avoid overhead
-        if not hasattr(self, "_interpolant"):
-            self._create_interpolant()
-        return self._interpolant(t)
-
-    __call__ = interpolate


Test patch:
diff --git a/tests/unit/algorithms/GenericTransportTest.py b/tests/unit/algorithms/GenericTransportTest.py
index 66bf320ef..95c2c23d5 100644
--- a/tests/unit/algorithms/GenericTransportTest.py
+++ b/tests/unit/algorithms/GenericTransportTest.py
@@ -54,7 +54,12 @@ class GenericTransportTest:
         alg.settings['quantity'] = 'pore.mole_fraction'
         alg.set_value_BC(pores=self.net.pores('top'), values=1)
         alg.set_value_BC(pores=self.net.pores('bottom'), values=0)
-        alg.run()
+        soln = alg.run()
+        # Check returned data type
+        from openpnm.algorithms._solution import SteadyStateSolution
+        assert isinstance(soln, SteadyStateSolution)
+        # Ensure solution object is attached to the algorithm
+        assert isinstance(alg.soln, SteadyStateSolution)
 
     def test_two_value_conditions(self):
         alg = op.algorithms.GenericTransport(network=self.net,
diff --git a/tests/unit/algorithms/ReactiveTransportTest.py b/tests/unit/algorithms/ReactiveTransportTest.py
index 4646908c7..b1f6847d8 100644
--- a/tests/unit/algorithms/ReactiveTransportTest.py
+++ b/tests/unit/algorithms/ReactiveTransportTest.py
@@ -169,7 +169,7 @@ class ReactiveTransportTest:
         self.alg.settings._update({'relaxation_quantity': 20.0,
                                    'newton_maxiter': 25})
         self.alg.run()
-        assert not self.alg.is_converged
+        assert not self.alg.soln.is_converged
 
     # FIXME: we no longer want to throw exception when maxiter is reached
     def test_check_divergence_if_maxiter_reached(self):
diff --git a/tests/unit/algorithms/TransientReactiveTransportTest.py b/tests/unit/algorithms/TransientReactiveTransportTest.py
index ec3258e61..9c2dba520 100644
--- a/tests/unit/algorithms/TransientReactiveTransportTest.py
+++ b/tests/unit/algorithms/TransientReactiveTransportTest.py
@@ -46,9 +46,10 @@ class TransientReactiveTransportTest:
     def test_transient_solution(self):
         out = self.alg.run(x0=0, tspan=(0, 1), saveat=0.1)
         # Test datatype
-        assert isinstance(out, op.integrators.TransientSolution)
+        from openpnm.algorithms._solution import TransientSolution
+        assert isinstance(out, TransientSolution)
         # Ensure solution object is attached to the algorithm
-        assert isinstance(self.alg.soln, op.integrators.TransientSolution)
+        assert isinstance(self.alg.soln, TransientSolution)
         # Test shape
         nt.assert_array_equal(self.alg.soln.shape, (self.alg.Np, 11))
         # Test stored time points