From b31fde779070a275364c43e5d18aeb571e553c2d Mon Sep 17 00:00:00 2001 From: Mike McKague Date: Fri, 22 Dec 2023 12:21:19 -0500 Subject: [PATCH 1/2] Added more integrator classes to current _scipy.py --- openpnm/integrators/_scipy.py | 43 +++++++++++++++++++++++++++++------ 1 file changed, 36 insertions(+), 7 deletions(-) diff --git a/openpnm/integrators/_scipy.py b/openpnm/integrators/_scipy.py index d48e25ee5a..cd5f538d4e 100644 --- a/openpnm/integrators/_scipy.py +++ b/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) From a6cc209085aaf0f02ab9074ba68a9eec689ea759 Mon Sep 17 00:00:00 2001 From: Mike McKague Date: Mon, 6 May 2024 11:14:15 -0400 Subject: [PATCH 2/2] Added unit tests for integrator classes --- tests/unit/algorithms/IntegratorsTest.py | 69 ++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 tests/unit/algorithms/IntegratorsTest.py diff --git a/tests/unit/algorithms/IntegratorsTest.py b/tests/unit/algorithms/IntegratorsTest.py new file mode 100644 index 0000000000..30bf45d760 --- /dev/null +++ b/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)()