fix: Improved Radau5ODE handling of final time

CHANGELOG

@@ -1,4 +1,7 @@
 --- CHANGELOG ---
+--- Assimulo-3.7.3---
+    * Improved Radau5ODE handling of steps ending close to final time. 
+
 --- Assimulo-3.7.2---
     * Fixed issue with 'report_continuously = True' for backwards simulation.
 

tests/solvers/utils.py → tests/conftest.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python 
 # -*- coding: utf-8 -*-
 
-# Copyright (C) 2025 Modelon AB
+# Copyright (C) 2026 Modelon AB
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
@@ -15,7 +15,9 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program. If not, see <http://www.gnu.org/licenses/>.
 
-# Testing utilities shared in a variety of solver tests
+import pytest
+import dataclasses
+from typing import Union
 
 import numpy as np
 
@@ -111,11 +113,12 @@ def init_mode(self, solver):
         solver.y[1] = (-1.0 if solver.sw[1] else 3.0)
         solver.y[2] = (0.0 if solver.sw[2] else 2.0)
 
+
 class Eval_Failure(Explicit_Problem):
     """Problem for testing evaluation failures starting from a given time point and 
     aborting on BaseExceptions."""
     y0 = np.array([1.])
-    def __init__(self, t_failure = 0.5, max_evals = 1000):
+    def __init__(self, t_failure = 0.5, max_evals = -1):
         self.t_failure = t_failure
         self.max_evals = max_evals
         self.n_evals = 0
@@ -187,6 +190,7 @@ def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, even
         super().__init__(self.aux.f, y0, sw0 = sw0)
         self.jac = self.aux.jac
 
+
 class ImplicitProbBaseException(Implicit_Problem):
     def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, event_n = 5):
         self.aux = BaseExceptionAux(dim = dim, fcn = fcn, jac = jac, event = event, fcn_n = fcn_n, event_n = event_n)
@@ -199,6 +203,7 @@ def __init__(self, dim, fcn = False, jac = False, event = False, fcn_n = 5, even
             sw0 = np.array([1.])
         super().__init__(self.aux.f_impl, y0, yd0, sw0 = sw0)
 
+
 class ExplicitTimeEventCloseToFinalTime(Explicit_Problem):
     t0 = 0.0
     def __init__(self, tfinal):
@@ -213,7 +218,8 @@ def handle_event(self, solver, event_info):
 
     def time_events(self, t, y, sw = None):
         return self._time_event if self._time_event > t else None
-
+
+
 class ImplicitTimeEventCloseToFinalTime(Implicit_Problem):
     t0 = 0.0
     def __init__(self, tfinal):
@@ -228,4 +234,49 @@ def handle_event(self, solver, event_info):
 
     def time_events(self, t, y, yd, sw = None):
         return self._time_event if self._time_event > t else None
-
+
+@pytest.fixture
+def extended_problem():
+    return Extended_Problem()
+
+@pytest.fixture
+def eval_failure():
+    return Eval_Failure()
+
+@pytest.fixture
+def explicit_prob_base_exception_func_eval():
+    return ExplicitProbBaseException(dim = 2, fcn = True)
+
+@pytest.fixture
+def explicit_prob_base_exception_jac_eval():
+    return ExplicitProbBaseException(dim = 2, jac = True)
+
+@pytest.fixture
+def explicit_prob_base_exception_event():
+    return ExplicitProbBaseException(dim = 1, event = True, event_n = 3)
+
+@pytest.fixture
+def implicit_prob_base_exception_func_eval():
+    return ImplicitProbBaseException(dim = 2, fcn = True)
+
+
+@dataclasses.dataclass
+class TimeEventCloseToFinalTimeCase:
+    problem: Union[Explicit_Problem, Implicit_Problem]
+    tfinal: float
+
+@pytest.fixture(params = [1e-6, 1, 1e6])
+def explicit_prob_time_event_close_to_final_time(request):
+    tfinal = request.param
+    return TimeEventCloseToFinalTimeCase(
+        problem = ExplicitTimeEventCloseToFinalTime(tfinal = tfinal),
+        tfinal = tfinal
+        )
+
+@pytest.fixture(params = [1e-6, 1, 1e6])
+def implicit_prob_time_event_close_to_final_time(request):
+    tfinal = request.param
+    return TimeEventCloseToFinalTimeCase(
+        problem = ImplicitTimeEventCloseToFinalTime(tfinal = tfinal),
+        tfinal = tfinal
+        )

tests/solvers/test_euler.py

@@ -22,8 +22,6 @@
 import numpy as np
 import scipy.sparse as sps
 
-from .utils import Extended_Problem
-
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
 
@@ -41,10 +39,8 @@ def setup_class(cls):
         cls.problem = Explicit_Problem(f, y0)
         cls.simulator = ExplicitEuler(cls.problem)
 
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = ExplicitEuler(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = ExplicitEuler(extended_problem) #Create the solver
 
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True

tests/solvers/test_odepack.py

@@ -25,8 +25,6 @@
 import numpy as np
 import scipy.sparse as sps
 
-from .utils import Extended_Problem
-
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
 
@@ -92,10 +90,8 @@ def jac_sparse(t,y):
         cls.sim.rtol = 1e-6 #Default 1e-6
         cls.sim.usejac = False
 
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = LSODAR(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = LSODAR(extended_problem) #Create the solver
 
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True

tests/solvers/test_radau5.py

@@ -26,20 +26,9 @@
 import scipy.sparse as sps
 import numpy as np
 
-from .utils import (
-    Extended_Problem,
-    Eval_Failure,
-    ExplicitProbBaseException,
-    ImplicitProbBaseException,
-    ExplicitTimeEventCloseToFinalTime,
-    ImplicitTimeEventCloseToFinalTime,
-)
-
-
 import re
 float_regex = r"[\s]*[\d]*.[\d]*((e|E)(\+|\-)\d\d|)"
 
-
 class Test_Explicit_Radau5_Py:
     """
     Tests the explicit Radau solver (Python implementation).
@@ -90,10 +79,8 @@ def jac(t,y):
         cls.sim.usejac = False
 
     @pytest.mark.skip("Does not support state events")
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = _Radau5ODE(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = _Radau5ODE(extended_problem) #Create the solver
 
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -351,10 +338,8 @@ def jac_sparse(t,y):
         cls.sim.inith = 1.e-4 #Initial step-size
         cls.sim.usejac = False
 
-    def test_event_localizer(self):
-        exp_mod = Extended_Problem() #Create the problem
-
-        exp_sim = Radau5ODE(exp_mod) #Create the solver
+    def test_event_localizer(self, extended_problem):
+        exp_sim = Radau5ODE(extended_problem) #Create the solver
 
         exp_sim.verbosity = 0
         exp_sim.report_continuously = True
@@ -872,37 +857,33 @@ def test_linear_solver(self):
         with pytest.raises(Radau_Exception, match = err_msg.format('0', "<class 'int'>")):
             self.sim.linear_solver = 0
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, explicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation. Radau5 + C + explicit problem."""
-        prob = ExplicitProbBaseException(dim = 2, fcn = True)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_func_eval)
         with pytest.raises(BaseException, match = "f"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_jac(self):
+    def test_base_exception_interrupt_jac(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation. Radau5 + C + explicit problem."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_jac_eval)
         sim.usejac = True
 
         with pytest.raises(BaseException, match = "jac"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_jac_sparse(self):
+    def test_base_exception_interrupt_jac_sparse(self, explicit_prob_base_exception_jac_eval):
         """Test that BaseExceptions in jacobian terminate the simulation. Radau5 + C + explicit problem + sparse jac."""
-        prob = ExplicitProbBaseException(dim = 2, jac = True)
-        prob.jac_nnz = 2
-        sim = Radau5ODE(prob)
+        explicit_prob_base_exception_jac_eval.jac_nnz = 2
+        sim = Radau5ODE(explicit_prob_base_exception_jac_eval)
         sim.linear_solver = 'SPARSE'
         sim.usejac = True
 
         with pytest.raises(BaseException, match = "jac"):
             sim.simulate(1.)
 
-    def test_base_exception_interrupt_event_indicator(self):
+    def test_base_exception_interrupt_event_indicator(self, explicit_prob_base_exception_event):
         """Test that BaseExceptions in event indicator function resp. solout callback correctly terminate solution."""
-        prob = ExplicitProbBaseException(dim = 1, event = True, event_n = 3)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(explicit_prob_base_exception_event)
 
         with pytest.raises(BaseException, match = "event"):
             sim.simulate(1.)
@@ -946,20 +927,23 @@ def f(t, y):
             sim.simulate(1.0)
         assert any(sim.statistics[k] > 0 for k in sim.statistics.keys()), "No statistics was found to be stored"
 
-    def test_no_progress(self):
+    def test_no_progress(self, eval_failure):
         """Test example where solver cannot make progress past a given time."""
-        prob = Eval_Failure(t_failure = 0.5, max_evals = -1)
-        sim = Radau5ODE(prob)
+        sim = Radau5ODE(eval_failure)
 
         err_msg = "passed failure time"
         with pytest.raises(ValueError, match = re.escape(err_msg)):
             sim.simulate(1.0)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self, 
+                                       report_continuously, 
+                                       explicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = Radau5ODE(ExplicitTimeEventCloseToFinalTime(tfinal))
+        problem = explicit_prob_time_event_close_to_final_time.problem
+        tfinal = explicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = Radau5ODE(problem)
         exp_sim.report_continuously = report_continuously
         tt, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2)
         assert tt[-1] < tfinal # check final interval is skipped
@@ -972,6 +956,24 @@ def test_backwards_report_continuously(self):
         sim.report_continuously = True
         sim.simulate(0, ncp = 10)
 
+    @pytest.mark.parametrize("scale_factor", [1e-6, 1, 1e6])
+    def test_final_step_skip_due_to_t_final_rounding(self, scale_factor):
+        """Test the case of skipping a final step in case a step finishes an epsilon 
+        before the final time."""
+        mod = Explicit_Problem(lambda t, y: [0], y0 = [1], t0 = 0)
+
+        maxh = 1*scale_factor
+        inith = 1*scale_factor
+        eps = 1e-15*scale_factor
+        t_final = 2*scale_factor + eps
+
+        sim = Radau5ODE(mod)
+        sim.maxh = maxh
+        sim.inith = inith
+        t, _ = sim.simulate(t_final)
+
+        assert t[-1] == t_final
+        assert sim.get_statistics()["nsteps"] == 2
 
 class Test_Implicit_Radau5:
     """
@@ -1232,11 +1234,15 @@ def f(t, y, yd):
         with pytest.raises(Radau5Error):
             sim.simulate(1.)
 
-    @pytest.mark.parametrize("tfinal", [1e-6, 1, 1e6])
     @pytest.mark.parametrize("report_continuously", [True, False])
-    def test_event_close_to_final_time(self, tfinal, report_continuously):
+    def test_event_close_to_final_time(self,
+                                       report_continuously,
+                                       implicit_prob_time_event_close_to_final_time):
         """Test event close to final time."""
-        exp_sim = Radau5DAE(ImplicitTimeEventCloseToFinalTime(tfinal))
+        prob = implicit_prob_time_event_close_to_final_time.problem
+        tfinal = implicit_prob_time_event_close_to_final_time.tfinal
+
+        exp_sim = Radau5DAE(prob)
         exp_sim.report_continuously = report_continuously
         tt, _, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2)
         assert tt[-1] < tfinal # check final interval is skipped
@@ -1387,10 +1393,9 @@ def test_maxh(self):
         self.sim.simulate(0.5)
         assert max(np.diff(self.sim.t_sol))-np.finfo('double').eps <= 0.01
 
-    def test_base_exception_interrupt_fcn(self):
+    def test_base_exception_interrupt_fcn(self, implicit_prob_base_exception_func_eval):
         """Test that BaseExceptions in right-hand side terminate the simulation. Radau5 + C + implicit problem."""
-        prob = ImplicitProbBaseException(dim = 2, fcn = True)
-        sim = Radau5DAE(prob)
+        sim = Radau5DAE(implicit_prob_base_exception_func_eval)
 
         err_msg = "Unrecoverable exception encountered during callback to problem (right-hand side/jacobian)."
         with pytest.raises(Radau5Error, match = re.escape(err_msg)):