From a66fad3fe128386edf0b40fd745966bed7dde95c Mon Sep 17 00:00:00 2001 From: petermeisrimelmodelon Date: Mon, 26 Jan 2026 13:30:47 +0000 Subject: [PATCH 1/3] fix: Improved Radau5ODE handling of final time --- CHANGELOG | 3 +++ tests/solvers/test_radau5.py | 18 ++++++++++++++++++ tests/solvers/test_sundials.py | 18 ++++++++++++++++++ thirdparty/radau5/radau5.c | 7 +++++++ 4 files changed, 46 insertions(+) diff --git a/CHANGELOG b/CHANGELOG index 23d0a16b..3775e495 100644 --- a/CHANGELOG +++ b/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. diff --git a/tests/solvers/test_radau5.py b/tests/solvers/test_radau5.py index 7493fe9b..14819e72 100644 --- a/tests/solvers/test_radau5.py +++ b/tests/solvers/test_radau5.py @@ -972,6 +972,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: """ diff --git a/tests/solvers/test_sundials.py b/tests/solvers/test_sundials.py index dc56ed96..8f65ae10 100644 --- a/tests/solvers/test_sundials.py +++ b/tests/solvers/test_sundials.py @@ -901,6 +901,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 = CVode(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_IDA: @classmethod diff --git a/thirdparty/radau5/radau5.c b/thirdparty/radau5/radau5.c index 122c4737..74e9340e 100644 --- a/thirdparty/radau5/radau5.c +++ b/thirdparty/radau5/radau5.c @@ -626,6 +626,13 @@ static int _radcor(radau_mem_t *rmem, int n, FP_CB_f fcn, void *fcn_EXT, for (i = 1; i <= n; ++i) { scal[i] = atol[i] + rtol[i] * radau5_abs(y[i]); } + /* possibly round to h_output if sufficiently close; Note: this bypasses hmax; + * This avoid a too small stepsize in the following step that would case an error. + * Needs to be done before reporting solution via callback. */ + if (radau5_abs(*xend - xph) < 100*rmem->input->uround*radau5_abs(xph)) { + *x = *xend; + last = TRUE_; + } if (iout) { rmem->xsol = *x; for (i = 1; i <= n; ++i) { From d85146b2e511a868c7f2122adf95c9da2f1decb7 Mon Sep 17 00:00:00 2001 From: petermeisrimelmodelon Date: Mon, 26 Jan 2026 16:16:07 +0000 Subject: [PATCH 2/3] chore: refactoring some testing to start using more fixture --- tests/{solvers/utils.py => conftest.py} | 61 ++++++++++++++++-- tests/solvers/test_euler.py | 8 +-- tests/solvers/test_odepack.py | 8 +-- tests/solvers/test_radau5.py | 75 +++++++++------------- tests/solvers/test_sundials.py | 84 +++++++++++-------------- 5 files changed, 127 insertions(+), 109 deletions(-) rename tests/{solvers/utils.py => conftest.py} (84%) diff --git a/tests/solvers/utils.py b/tests/conftest.py similarity index 84% rename from tests/solvers/utils.py rename to tests/conftest.py index 3625673c..95e8aed7 100644 --- a/tests/solvers/utils.py +++ b/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 . -# 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 + ) diff --git a/tests/solvers/test_euler.py b/tests/solvers/test_euler.py index 05b6128e..e799c798 100644 --- a/tests/solvers/test_euler.py +++ b/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 diff --git a/tests/solvers/test_odepack.py b/tests/solvers/test_odepack.py index 1d91af7c..df89ff0e 100644 --- a/tests/solvers/test_odepack.py +++ b/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 diff --git a/tests/solvers/test_radau5.py b/tests/solvers/test_radau5.py index 14819e72..79e396b5 100644 --- a/tests/solvers/test_radau5.py +++ b/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', "")): 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 @@ -1250,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 @@ -1405,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)): diff --git a/tests/solvers/test_sundials.py b/tests/solvers/test_sundials.py index 8f65ae10..436205e7 100644 --- a/tests/solvers/test_sundials.py +++ b/tests/solvers/test_sundials.py @@ -23,15 +23,6 @@ from assimulo.exception import AssimuloException, TimeLimitExceeded, TerminateSimulation import numpy as np import scipy.sparse as sps -from .utils import ( - Extended_Problem, - Eval_Failure, - ExplicitProbBaseException, - ImplicitProbBaseException, - ExplicitTimeEventCloseToFinalTime, - ImplicitTimeEventCloseToFinalTime, -) - class Test_CVode: @@ -67,11 +58,9 @@ def f(t, y): assert np.all(t == np.arange(0,11)[::-1]) - def test_event_localizer(self): + def test_event_localizer(self, extended_problem): """Test that CVode internal event localization works correctly.""" - exp_mod = Extended_Problem() # Create the problem - - exp_sim = CVode(exp_mod) # Create the solver + exp_sim = CVode(extended_problem) # Create the solver exp_sim.verbosity = 0 exp_sim.report_continuously = True @@ -85,11 +74,9 @@ def test_event_localizer(self): assert y[-1][1] == pytest.approx(3.0) assert y[-1][2] == pytest.approx(2.0) - def test_event_localizer_external(self): + def test_event_localizer_external(self, extended_problem): """Test that CVode with Assimulo event localization works correctly.""" - exp_mod = Extended_Problem() # Create the problem - - exp_sim = CVode(exp_mod) # Create the solver + exp_sim = CVode(extended_problem) # Create the solver exp_sim.verbosity = 0 exp_sim.report_continuously = True @@ -784,11 +771,10 @@ def test_rtol_vector_sense(self): sim.rtol = [1e-6] sim.rtol = np.array([1e-6]) - def test_no_progress_force_min_h(self): + def test_no_progress_force_min_h(self, eval_failure): """Test example where CVode fails to make progress and minimal stepsize will be forced.""" - prob = Eval_Failure(t_failure = 0.5, max_evals = -1) - sim = CVode(prob) + sim = CVode(eval_failure) sim.maxstepshnil = 10 # default assert sim.minh == 0 err_msg = "The right-hand side function had repeated recoverable errors." @@ -796,11 +782,11 @@ def test_no_progress_force_min_h(self): sim.simulate(1.) assert sim.minh > 0 - def test_no_progress_maxstepshnil_not_active(self): + def test_no_progress_maxstepshnil_not_active(self, eval_failure): """Test example where CVode fails to make progress, but forcing minh is deactivated.""" - prob = Eval_Failure(t_failure = 0.5, max_evals = 1000) - sim = CVode(prob) + eval_failure.max_evals = 1000 + sim = CVode(eval_failure) sim.maxstepshnil = 0 assert sim.minh == 0 err_msg = "failed in an unrecoverable manner" @@ -814,11 +800,10 @@ def test_no_progress_maxstepshnil_not_active(self): (0, [0., 0.7, 1.]) ] ) - def test_maxstepshnil_not_enforced(self, ncp, ncp_list): + def test_maxstepshnil_not_enforced(self, ncp, ncp_list, eval_failure): """Test example where CVode fails to make progress, but minh will not be forced. CVode fails in ordinary ways.""" - prob = Eval_Failure(t_failure = 0.5, max_evals = -1) - sim = CVode(prob) + sim = CVode(eval_failure) assert sim.minh == 0 sim.report_continuously = False sim.maxsteps = 200 @@ -844,18 +829,16 @@ def test_maxstepshnil_set_invalid(self, val): with pytest.raises(TypeError, match = re.escape(msg)): sim.maxstepshnil = val - 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.""" - prob = ExplicitProbBaseException(dim = 2, fcn = True) - sim = CVode(prob) + sim = CVode(explicit_prob_base_exception_func_eval) msg = "The user-provided rhs function failed in an unrecoverable manner" with pytest.raises(CVodeError, match = re.escape(msg)): 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.""" - prob = ExplicitProbBaseException(dim = 2, jac = True) - sim = CVode(prob) + sim = CVode(explicit_prob_base_exception_jac_eval) sim.usejac = True msg = "The linear solvers setup function failed in an unrecoverable manner" @@ -863,11 +846,10 @@ def test_base_exception_interrupt_jac(self): sim.simulate(1.) @pytest.mark.skip(reason = "This works, but causes a segfault in the deconstructor") - 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.""" - prob = ExplicitProbBaseException(dim = 2, jac = True) - prob.jac_nnz = 2 - sim = CVode(prob) + explicit_prob_base_exception_jac_eval.jac_nnz = 2 + sim = CVode(explicit_prob_base_exception_jac_eval) sim.linear_solver = 'SPARSE' sim.usejac = True @@ -875,20 +857,23 @@ def test_base_exception_interrupt_jac_sparse(self): with pytest.raises(CVodeError, match = re.escape(msg)): 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 = CVode(prob) + sim = CVode(explicit_prob_base_exception_event) msg = "The rootfinding function failed in an unrecoverable manner" with pytest.raises(CVodeError, match = re.escape(msg)): 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, + explicit_prob_time_event_close_to_final_time): """Test event close to final time.""" - exp_sim = CVode(ExplicitTimeEventCloseToFinalTime(tfinal)) + problem = explicit_prob_time_event_close_to_final_time.problem + tfinal = explicit_prob_time_event_close_to_final_time.tfinal + + exp_sim = CVode(problem) exp_sim.report_continuously = report_continuously tt, _ = exp_sim.simulate(tfinal = tfinal, ncp = 2) assert tt[-1] < tfinal # check final interval is skipped @@ -1353,20 +1338,23 @@ def completed_step(solver): assert len(sim.t_sol) == sim.statistics["nsteps"] + 1 assert nsteps == sim.statistics["nsteps"] - 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 = IDA(prob) + sim = IDA(implicit_prob_base_exception_func_eval) err_msg = "The user-provided residual function failed in an unrecoverable manner." with pytest.raises(IDAError, match = re.escape(err_msg)): 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 = IDA(ImplicitTimeEventCloseToFinalTime(tfinal)) + prob = implicit_prob_time_event_close_to_final_time.problem + tfinal = implicit_prob_time_event_close_to_final_time.tfinal + + exp_sim = IDA(prob) exp_sim.report_continuously = report_continuously tt, _, _ = exp_sim.simulate(tfinal = tfinal, ncp = 10) assert tt[-1] < tfinal # check final interval is skipped From 1f8abcd8559d0e5a0336fa6b08eaedd6801730de Mon Sep 17 00:00:00 2001 From: petermeisrimelmodelon Date: Tue, 27 Jan 2026 09:44:23 +0000 Subject: [PATCH 3/3] improving comments --- thirdparty/radau5/radau5.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/thirdparty/radau5/radau5.c b/thirdparty/radau5/radau5.c index 74e9340e..5583fd65 100644 --- a/thirdparty/radau5/radau5.c +++ b/thirdparty/radau5/radau5.c @@ -626,9 +626,9 @@ static int _radcor(radau_mem_t *rmem, int n, FP_CB_f fcn, void *fcn_EXT, for (i = 1; i <= n; ++i) { scal[i] = atol[i] + rtol[i] * radau5_abs(y[i]); } - /* possibly round to h_output if sufficiently close; Note: this bypasses hmax; - * This avoid a too small stepsize in the following step that would case an error. - * Needs to be done before reporting solution via callback. */ + /* possibly round to xend if sufficiently close; Note: this bypasses hmax; + * This avoids an error due to too small stepsize in the following (final) step. + * Done before solout callback such that solution is reported for x = xend. */ if (radau5_abs(*xend - xph) < 100*rmem->input->uround*radau5_abs(xph)) { *x = *xend; last = TRUE_;