Unify RHS and SOLOUT user interrupts under a single idid exit path (#52)

* Add RHS callback interrupt path

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/4b6afe5b-daf3-4c54-a604-ea10ab860a5b

Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>

* Verify RHS interrupt regression coverage

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/4b6afe5b-daf3-4c54-a604-ea10ab860a5b

Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>

* Finalize unified callback interrupt handling

Agent-Logs-Url: https://github.com/ivan-pi/stiff3/sessions/4b6afe5b-daf3-4c54-a604-ea10ab860a5b

Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>

---------

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Co-authored-by: ivan-pi <21085643+ivan-pi@users.noreply.github.com>

Author: Copilot

README.md

@@ -133,10 +133,11 @@ program vanpol
 
 contains
 
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     f(1) = y(2)
     f(2) = mu*(1.0_wp - y(1)**2)*y(2) - y(1)
   end subroutine
@@ -191,7 +192,7 @@ The integer exit flag `idid` reports solver status:
 
 - `0` — successful completion at `x1`
 - `-1` — LU factorization failed because the Jacobian system matrix was singular
-- `-2` — integration interrupted by the user `solout` callback (`irtrn < 0`)
+- `-2` — integration interrupted by a user callback (`ires < 0` in `fun` or `irtrn < 0` in `solout`)
 - `-3` — step-size underflow occurred during bisection (`abs(h) <= spacing(x_current)`)
 
 When using this explicit-workspace overload together with `solout`, accepted-step dense output is available through the generic interface `stiff3_interp`:

example/e5_chemical.f90

@@ -219,10 +219,11 @@ function ycoord(value, ymin, ymax, height, top, bottom) result(yplot)
       (value - ymin)*(real(height,wp) - top - bottom)/(ymax - ymin)
   end function
 
-  subroutine fun(n, y, f)
+  subroutine fun(n, y, f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     real(wp) :: prod1, prod2, prod3, prod4
 
     prod1 = 7.89e-10_wp*y(1)

example/hires.f90

@@ -107,11 +107,11 @@ subroutine write_gnuplot_data(filename, eps, cpu, err, stats)
     close(unit)
   end subroutine
 
-  subroutine fun(n, y, f)
+  subroutine fun(n, y, f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
-
+    integer, intent(inout) :: ires
     f(1) = -1.71_wp*y(1) + 0.43_wp*y(2) + 8.32_wp*y(3) + 0.0007_wp
     f(2) = 1.71_wp*y(1) - 8.75_wp*y(2)
     f(3) = -10.03_wp*y(3) + 0.43_wp*y(4) + 0.035_wp*y(5)

example/lorenz.f90

@@ -25,10 +25,11 @@ program lorenz
 
 contains
 
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     f(1) = sigma*(y(2) - y(1))
     f(2) = y(1)*(rho - y(3)) - y(2)
     f(3) = y(1)*y(2) - beta*y(3)

example/oregonator.f90

@@ -253,11 +253,11 @@ function ycoord(value, ymin, ymax, height, top, bottom) result(yplot)
       (value - ymin)*(real(height,wp) - top - bottom)/(ymax - ymin)
   end function
 
-  subroutine fun(n, y, f)
+  subroutine fun(n, y, f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
-
+    integer, intent(inout) :: ires
     f(1) = 77.27_wp*(y(2) + y(1)*(1.0_wp - 8.375e-6_wp*y(1) - y(2)))
     f(2) = (y(3) - (1.0_wp + y(1))*y(2))/77.27_wp
     f(3) = 0.161_wp*(y(1) - y(3))

example/pendant_drop.f90

@@ -36,10 +36,11 @@ program pendant_drop
 
 contains
 
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     f(1) = pL - drho*y(3) - sin(y(1))/y(2)
     f(2) = cos(y(1))
     f(3) = sin(y(1))

example/predator_prey.f90

@@ -26,10 +26,11 @@ program predator_prey
 
 contains
 
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     f(1) = alpha*y(1) - beta*y(1)*y(2)
     f(2) = delta*y(1)*y(2) - gamma*y(2)
   end subroutine

example/ring_modulator.f90

@@ -115,11 +115,11 @@ subroutine write_gnuplot_data(filename, eps, cpu, err, stats)
     close(unit)
   end subroutine
 
-  subroutine fun(n, y, f)
+  subroutine fun(n, y, f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
-
+    integer, intent(inout) :: ires
     real(wp), parameter :: c = 1.6e-8_wp
     real(wp), parameter :: cs = 2.0e-12_wp
     real(wp), parameter :: cp = 1.0e-8_wp

example/robertson.f90

@@ -9,10 +9,11 @@ module robertson_rate
   !
   ! Robertson rate equations
   !
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
+    integer, intent(inout) :: ires
     f(1) = -4.e-2_wp*y(1)+1.e4_wp*y(2)*y(3)
     f(2) = -f(1)-3.e7_wp*y(2)**2
     f(3) = 3.e7_wp*y(2)**2

example/three_equation_system.f90

@@ -29,11 +29,11 @@ program three_equation_system
 
 contains
 
-  subroutine fun(n,y,f)
+  subroutine fun(n,y,f, ires)
     integer, intent(in) :: n
     real(wp), intent(in) :: y(n)
     real(wp), intent(inout) :: f(n)
-
+    integer, intent(inout) :: ires
     f(1) = -k1*y(1) - k2*y(1)*y(3)
     f(2) = -k3*y(2)*y(3)
     f(3) = -k1*y(1) - k2*y(1)*y(3) - k3*y(2)*y(3)