From 43b36fdbc4e09190319fe29506b5ddb5fdbc8a9d Mon Sep 17 00:00:00 2001
From: Gabriel Gerlero <gerlero@users.noreply.github.com>
Date: Tue, 10 Oct 2023 22:13:17 -0300
Subject: [PATCH] Work around solve_ivp issue

---
 CHANGELOG.md            |  7 +++++++
 fronts/__init__.py      |  2 +-
 fronts/_semiinfinite.py | 42 ++++++++++++++++++++++++++---------------
 3 files changed, 35 insertions(+), 16 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 184b257..e451596 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/fronts/__init__.py b/fronts/__init__.py
index 5ff113f..530d192 100644
--- a/fronts/__init__.py
+++ b/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, 
diff --git a/fronts/_semiinfinite.py b/fronts/_semiinfinite.py
index 05e1597..5938e35 100644
--- a/fronts/_semiinfinite.py
+++ b/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: