Added more Scipy integrators to handle stiff problems

openpnm/integrators/_scipy.py

@@ -3,13 +3,22 @@
 from openpnm.algorithms._solution import TransientSolution
 from openpnm.integrators import Integrator
 
-__all__ = ['ScipyRK45']
+__all__ = ['ScipyIntegrator',
+           'ScipyRK45',
+           'ScipyBDF',
+           'ScipyLSODA']
 
 
-class ScipyRK45(Integrator):
-    """Integrator class based on SciPy's implementation of RK45"""
+class ScipyIntegrator(Integrator):
+    """Integrator class based on SciPy's ODE solvers"""
 
-    def __init__(self, atol=1e-6, rtol=1e-6, verbose=False, linsolver=None):
+    def __init__(self,
+                 method="RK45",
+                 atol=1e-6,
+                 rtol=1e-6,
+                 verbose=False,
+                 linsolver=None):
+        self.method = method
         self.atol = atol
         self.rtol = rtol
         self.verbose = verbose
@@ -47,10 +56,30 @@ def solve(self, rhs, x0, tspan, saveat, **kwargs):
             "atol": self.atol,
             "rtol": self.rtol,
             "t_eval": saveat,
-            # FIXME: uncomment next line when/if scipy#11815 is merged
-            # "verbose": self.verbose,
+            # "verbose": self.verbose,  # Uncomment if supported in scipy
         }
-        sol = solve_ivp(rhs, tspan, x0, method="RK45", **options)
+        sol = solve_ivp(rhs, tspan, x0, method=self.method, **options)
         if sol.success:
             return TransientSolution(sol.t, sol.y)
         raise Exception(sol.message)
+
+
+class ScipyRK45(ScipyIntegrator):
+    """Integrator class using SciPy's RK45 method"""
+
+    def __init__(self, method="RK45", **kwargs):
+        super().__init__(method=method, **kwargs)
+
+
+class ScipyBDF(ScipyIntegrator):
+    """Integrator class using SciPy's BDF method"""
+
+    def __init__(self, method="BDF", **kwargs):
+        super().__init__(method=method, **kwargs)
+
+
+class ScipyLSODA(ScipyIntegrator):
+    """Integrator class using SciPy's LSODA method"""
+
+    def __init__(self, method="LSODA", **kwargs):
+        super().__init__(method=method, **kwargs)

tests/unit/algorithms/IntegratorsTest.py

@@ -0,0 +1,69 @@
+import openpnm as op
+import numpy as np
+import numpy.testing as nt
+
+
+class IntegratorsTest:
+
+    def setup_class(self):
+        np.random.seed(10)
+        self.net = op.network.Cubic([4, 4, 1], spacing=1e-4)
+        geo_mods = op.models.collections.geometry.spheres_and_cylinders.copy()
+        self.net.add_model_collection(geo_mods)
+        self.net.regenerate_models()
+        self.phase = op.phase.Water(network=self.net)
+        self.phase['pore.diffusivity'] = 1e-10
+        self.phase['pore.concentration'] = 100
+        phys_mods = op.models.collections.physics.basic.copy()
+        self.phase.add_model_collection(phys_mods)
+        self.phase.regenerate_models()
+        rxn = op.models.physics.source_terms.standard_kinetics
+        self.phase['pore.reaction_constant'] = 1e-14
+        self.phase['pore.reaction_order'] = 1
+        self.phase.add_model(propname='pore.source',
+                             model=rxn,
+                             X='pore.concentration',
+                             prefactor='pore.reaction_constant',
+                             exponent='pore.reaction_order')
+        self.alg = op.algorithms.TransientReactiveTransport(network=self.net,
+                                                            phase=self.phase)
+        self.alg.settings['quantity'] = 'pore.concentration'
+        self.alg.settings['conductance'] = 'throat.diffusive_conductance'
+        self.alg.set_source(pores=[0], propname='pore.source', mode='add')
+
+    def test_scipy_RK45(self):
+        x0 = np.ones(self.net.Np)*100
+        tspan = (0, 10)
+        saveat = 1.0
+        integrator = op.integrators.ScipyRK45()
+        self.alg.run(x0, tspan, saveat, integrator)
+        x = self.alg.x
+        nt.assert_allclose(x.mean(), 110.97095388, rtol=1e-5)
+
+    def test_scipy_BDF(self):
+        x0 = np.ones(self.net.Np)*100
+        tspan = (0, 10)
+        saveat = 1.0
+        integrator = op.integrators.ScipyBDF()
+        self.alg.run(x0, tspan, saveat, integrator)
+        x = self.alg.x
+        nt.assert_allclose(x.mean(), 110.97112100, rtol=1e-5)
+
+    def test_scipy_LSODA(self):
+        x0 = np.ones(self.net.Np)*100
+        tspan = (0, 10)
+        saveat = 1.0
+        integrator = op.integrators.ScipyLSODA()
+        self.alg.run(x0, tspan, saveat, integrator)
+        x = self.alg.x
+        nt.assert_allclose(x.mean(), 110.97097456, rtol=1e-5)
+
+
+if __name__ == '__main__':
+    t = IntegratorsTest()
+    t.setup_class()
+    self = t
+    for item in t.__dir__():
+        if item.startswith('test'):
+            print(f'Running test: {item}')
+            t.__getattribute__(item)()