Repository: pmeal/openpnm
Issue: Solvers don't check for convergence when maximum iterations is reached

This is a corner-case that @mkaguer has recently found. Basically, when `max_iter` is reached, the solver gets tricked that the solution has converged. This is because unless `update_iterative_props` is called, the solution `x`, not surprisingly, satisfies `Ax-b~0`. It's only when we update the `A` matrix with the new solution (due to a nonlinear source term for instance), that `Ax-b~0` no longer necessarily holds. When `max_iter` is reached, the `for` loop is ended, and `update_iterative_props` which happens to be at the top of the loop never gets called again. 

Changed files:
tests/unit/algorithms/TransientMultiphysicsNernstPlanckSolverTest.py, openpnm/algorithms/ReactiveTransport.py, tests/unit/algorithms/ReactiveTransportTest.py

diff --git a/openpnm/algorithms/ReactiveTransport.py b/openpnm/algorithms/ReactiveTransport.py
index 0f8489626..4d4c5b4a3 100644
--- a/openpnm/algorithms/ReactiveTransport.py
+++ b/openpnm/algorithms/ReactiveTransport.py
@@ -135,7 +135,7 @@ class ReactiveTransport(GenericTransport):
         if conductance:
             self.settings['conductance'] = conductance
         if nlin_max_iter:
-            self.settings['max_iter'] = nlin_max_iter
+            self.settings['nlin_max_iter'] = nlin_max_iter
         if relaxation_source:
             self.settings['relaxation_source'] = relaxation_source
         if relaxation_quantity:
@@ -223,7 +223,7 @@ class ReactiveTransport(GenericTransport):
             the algorithm.
 
         """
-        if type(propnames) is str:  # Convert string to list if necessary
+        if isinstance(propnames, str):  # Convert string to list if necessary
             propnames = [propnames]
         d = self.settings["variable_props"]
         self.settings["variable_props"] = list(set(d) | set(propnames))
@@ -333,28 +333,38 @@ class ReactiveTransport(GenericTransport):
         max_it = self.settings['nlin_max_iter']
         # Write initial guess to algorithm obj (for _update_iterative_props to work)
         self[quantity] = x = x0
+        # Update A and b based on self[quantity]
+        self._update_A_and_b()
+        # Just in case you got a lucky guess, i.e. x0!
+        if self._is_converged():
+            logger.info(f'Solution converged: {self._get_residual():.4e}')
+            return x
 
         for itr in range(max_it):
-            # Update iterative properties on phase, geometries, and physics
-            self._update_iterative_props()
-            # Build A and b, apply BCs/source terms
-            self._build_A()
-            self._build_b()
-            self._apply_BCs()
-            self._apply_sources()
-            # Check solution convergence
-            res = self._get_residual()
-            if itr >= 1 and self._is_converged():
-                logger.info(f'Solution converged: {res:.4e}')
-                return x
-            logger.info(f'Tolerance not met: {res:.4e}')
             # Solve, use relaxation, and update solution on algorithm obj
             self[quantity] = x = self._solve(x0=x) * w + x * (1 - w)
+            self._update_A_and_b()
+            # Check solution convergence
+            if self._is_converged():
+                logger.info(f'Solution converged: {self._get_residual():.4e}')
+                return x
+            logger.info(f'Tolerance not met: {self._get_residual():.4e}')
 
-        # Check solution convergence after max_it iterations
         if not self._is_converged():
             raise Exception(f"Not converged after {max_it} iterations.")
 
+    def _update_A_and_b(self):
+        r"""
+        Updates A and b based on the most recent solution stored on algorithm object.
+        """
+        # Update iterative properties on phase, geometries, and physics
+        self._update_iterative_props()
+        # Build A and b, apply BCs/source terms
+        self._build_A()
+        self._build_b()
+        self._apply_BCs()
+        self._apply_sources()
+
     def _get_iterative_props(self):
         r"""
         Find and return properties that need to be iterated while running the
diff --git a/tests/unit/algorithms/ReactiveTransportTest.py b/tests/unit/algorithms/ReactiveTransportTest.py
index ce5572d29..ee8382ec7 100644
--- a/tests/unit/algorithms/ReactiveTransportTest.py
+++ b/tests/unit/algorithms/ReactiveTransportTest.py
@@ -26,6 +26,22 @@ class ReactiveTransportTest:
         self.alg = op.algorithms.ReactiveTransport(network=self.net,
                                                    phase=self.phase)
 
+    def test_setup(self):
+        temp = self.alg.settings.copy()
+        self.alg.setup(
+            conductance="throat.cond",
+            quantity="pore.test",
+            nlin_max_iter=123,
+            relaxation_source=1.23,
+            relaxation_quantity=3.21
+        )
+        assert self.alg.settings["conductance"] == "throat.cond"
+        assert self.alg.settings["quantity"] == "pore.test"
+        assert self.alg.settings["nlin_max_iter"] == 123
+        assert self.alg.settings["relaxation_source"] == 1.23
+        assert self.alg.settings["relaxation_quantity"] == 3.21
+        self.alg.settings = temp
+
     def test_set_variable_props(self):
         assert len(self.alg.settings["variable_props"]) == 0
         self.alg._set_variable_props(propnames="pore.pressure")
@@ -160,6 +176,19 @@ class ReactiveTransportTest:
         with pytest.raises(Exception):
             self.alg.run()
 
+    def test_check_divergence_if_maxiter_reached(self):
+        self.alg.reset(bcs=True, source_terms=True)
+        self.alg.setup(nlin_max_iter=2)
+        self.alg.settings.update({'conductance': 'throat.diffusive_conductance',
+                                  'quantity': 'pore.concentration'})
+        self.alg.set_source(pores=self.net.pores('bottom'), propname='pore.reaction')
+        self.alg.set_value_BC(pores=self.net.pores('top'), values=1.0)
+        self.alg.settings['relaxation_quantity'] = 1.0
+        self.alg.settings['relaxation_source'] = 1.0
+        with pytest.raises(Exception):
+            self.alg.run()
+        self.alg.setup(nlin_max_iter=5000)
+
     def test_variable_conductance(self):
         self.alg.reset(bcs=True, source_terms=True, variable_props=True)
 
diff --git a/tests/unit/algorithms/TransientMultiphysicsNernstPlanckSolverTest.py b/tests/unit/algorithms/TransientMultiphysicsNernstPlanckSolverTest.py
index 0025d0aac..c36d5009c 100644
--- a/tests/unit/algorithms/TransientMultiphysicsNernstPlanckSolverTest.py
+++ b/tests/unit/algorithms/TransientMultiphysicsNernstPlanckSolverTest.py
@@ -168,9 +168,9 @@ class TransientMultiphysicsNernstPlanckSolverTest:
             "pore.concentration.Na_mix_01@3000",
             "pore.concentration.Cl_mix_01@3000",
             "pore.potential@3000",
-            "pore.concentration.Na_mix_01@4000",
-            "pore.concentration.Cl_mix_01@4000",
-            "pore.potential@4000",
+            "pore.concentration.Na_mix_01@3500",
+            "pore.concentration.Cl_mix_01@3500",
+            "pore.potential@3500",
         ]
         assert set(times).issubset(set(self.sw.keys()))