Work around solve_ivp issue

CHANGELOG.md

@@ -5,6 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## [1.1.2] - Unreleased
+
+### Fixed
+
+- Work around a SciPy bug that can cause `solve()` to fail with an unexpected `ValueError` (https://github.com/scipy/scipy/issues/17066).
+
 ## [1.1.1] - 2023-10-09
 
 ### Added
@@ -198,6 +204,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 First public pre-release version.
 
+[1.1.2]: https://github.com/gerlero/fronts/compare/v1.1.1...v1.1.2
 [1.1.1]: https://github.com/gerlero/fronts/compare/v1.1.0...v1.1.1
 [1.1.0]: https://github.com/gerlero/fronts/compare/v1.0.3...v1.1.0
 [1.0.3]: https://github.com/gerlero/fronts/compare/v1.0.2...v1.0.3

fronts/__init__.py

@@ -3,7 +3,7 @@
 transformation.
 """
 
-__version__ = '1.1.1'
+__version__ = '1.1.2'
 
 from ._boltzmann import ode, BaseSolution, o, do_dr, do_dt, r, t, as_o
 from ._semiinfinite import (solve, solve_flowrate, solve_from_guess, Solution, 

fronts/_semiinfinite.py

@@ -300,21 +300,33 @@ def integrate(self, b, d_dob):
 
         with np.errstate(divide='ignore', invalid='ignore'):
 
-            if self._method == 'explicit':
-                ivp_result = solve_ivp(self._fun,
-                                       t_span=(self._ob, np.inf),
-                                       y0=(b, d_dob),
-                                       method='DOP853',
-                                       events=self._events,
-                                       dense_output=True)
-            else:
-                ivp_result = solve_ivp(self._fun,
-                                       t_span=(self._ob, np.inf),
-                                       y0=(b, d_dob),
-                                       method='Radau',
-                                       jac=self._jac,
-                                       events=self._events,
-                                       dense_output=True)
+            try:
+                if self._method == 'explicit':
+                    ivp_result = solve_ivp(self._fun,
+                                            t_span=(self._ob, np.inf),
+                                            y0=(b, d_dob),
+                                            method='DOP853',
+                                            events=self._events,
+                                            dense_output=True)
+                else:
+                    ivp_result = solve_ivp(self._fun,
+                                            t_span=(self._ob, np.inf),
+                                            y0=(b, d_dob),
+                                            method='Radau',
+                                            jac=self._jac,
+                                            events=self._events,
+                                            dense_output=True)
+            except ValueError as e:
+                # Workaround for https://github.com/scipy/scipy/issues/17066
+                if str(e) == "`ts` must be strictly increasing or decreasing.":
+                    return self.Result(b=b,
+                                       d_dob=d_dob,
+                                       i_residual=self._theta_direction*np.inf,
+                                       D_calls=None,
+                                       o=None,
+                                       sol=None)
+                else:
+                    raise
 
         if ivp_result.success and ivp_result.t_events[0].size == 1: