Skip to content

Allow users to flag a time step (x values) as invalid within the derivative function #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mind-bending-forks wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Allow users to flag a time step (x values) as invalid within the derivative function #44

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
0168ffc
Add an error message corresponding to the need to stop due to encount…
mind-bending-forks Apr 21, 2026
faac74c
Add and initialise a flag to track step validity, similar to the stop…
mind-bending-forks Apr 21, 2026
29d81ce
Add a user-callable method to allow flagging of a timestep as invalid
mind-bending-forks Apr 21, 2026
8f7936f
Raise exception if step is marked as invalid in scenrio where timeste…
mind-bending-forks Apr 21, 2026
98bfe21
Bisect the timestep and try again in scenarios where the timestep can…
mind-bending-forks Apr 21, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
80 changes: 67 additions & 13 deletions src/rklib_module.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ module rklib_module
real(wp),parameter :: zero = 0.0_wp

integer,parameter :: max_error_len = 100 !! max size of error message strings
integer,parameter,public :: RKLIB_ERROR_INVALID_F = -11
integer,parameter,public :: RKLIB_ERROR_TOO_MANY_STEPS = -10
integer,parameter,public :: RKLIB_ERROR_INVALID_RTOL_SIZE = -9
integer,parameter,public :: RKLIB_ERROR_INVALID_ATOL_SIZE = -8
Expand All @@ -50,8 +51,9 @@ module rklib_module
integer,parameter,public :: RKLIB_ERROR_G_NOT_ASSOCIATED = -2
integer,parameter,public :: RKLIB_ERROR_F_NOT_ASSOCIATED = -1
integer,parameter,public :: RKLIB_ERROR_NONE = 0
character(len=max_error_len),dimension(RKLIB_ERROR_TOO_MANY_STEPS:RKLIB_ERROR_NONE),parameter :: &
character(len=max_error_len),dimension(RKLIB_ERROR_INVALID_F:RKLIB_ERROR_NONE),parameter :: &
rklib_error_messages = [&
'Cannot avoid invalid timestep ', & ! -11
'Too many steps ', & ! -10
'Invalid size for rtol array ', & ! -9
'Invalid size for atol array ', & ! -8
Expand Down Expand Up @@ -122,6 +124,7 @@ module rklib_module

integer :: istatus = 0 !! status code
logical :: stopped = .false. !! if user has stopped the integration in `f` or `report`.
logical :: invalid = .false. !! if user has flagged x as invalid at the specified time in `f`
integer :: num_steps = 0 !! number of accepted steps taken
integer :: max_number_of_steps = huge(1) !! maximum number of steps to take
integer :: report_rate = 1 !! how often to call report function:
Expand All @@ -147,6 +150,7 @@ module rklib_module

procedure,public :: destroy !! destructor
procedure,public :: stop => rk_class_stop !! user-callable method to stop the integration
procedure,public :: flag_invalid => rk_class_flag_invalid !! user-callable method to flag a timestep as invalid
procedure,public :: status => rk_class_status !! get status code and message
procedure,public :: failed

Expand Down Expand Up @@ -388,6 +392,16 @@ subroutine rk_class_stop(me)
end subroutine rk_class_stop
!*****************************************************************************************

!*****************************************************************************************
!>
! User-callable method to flag a timestep as invalid

subroutine rk_class_flag_invalid(me)
class(rk_class),intent(inout) :: me
me%invalid = .true.
end subroutine rk_class_flag_invalid
!*****************************************************************************************

!*****************************************************************************************
!>
! Returns true if there was an error.
Expand Down Expand Up @@ -551,6 +565,7 @@ subroutine initialize_rk_class(me,n,f,report,g,stop_on_errors,&
! reset internal variables:
me%num_steps = 0
me%stopped = .false.
me%invalid = .false.

end subroutine initialize_rk_class
!*****************************************************************************************
Expand All @@ -564,6 +579,7 @@ subroutine begin_integration_rk_class(me)
call me%clear_exception()
me%num_steps = 0
me%stopped = .false.
me%invalid = .false.
end subroutine begin_integration_rk_class
!*****************************************************************************************

Expand Down Expand Up @@ -655,6 +671,10 @@ subroutine integrate_fixed_step(me,t0,x0,h,tf,xf)
if (last) dt = tf-t !
call me%step(t,x,dt,xf)
if (me%stopped) return
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
me%num_steps = me%num_steps + 1
if (me%num_steps > me%max_number_of_steps) then
call me%raise_exception(RKLIB_ERROR_TOO_MANY_STEPS)
Expand Down Expand Up @@ -745,6 +765,10 @@ subroutine integrate_to_event_fixed_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
end if
call me%step(t,x,dt,xf)
if (me%stopped) return
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
me%num_steps = me%num_steps + 1
if (me%num_steps > me%max_number_of_steps) then
call me%raise_exception(RKLIB_ERROR_TOO_MANY_STEPS)
Expand Down Expand Up @@ -821,6 +845,10 @@ function solver_func(delt) result(g)
!take a step from t to t+delt and evaluate g function:
call me%step(t,x,delt,g_xf)
if (me%stopped) return
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
call me%g(t+delt,g_xf,g)

end function solver_func
Expand Down Expand Up @@ -1224,16 +1252,27 @@ subroutine integrate_variable_step(me,t0,x0,h,tf,xf)
if (me%stopped) return

if (me%stepsize_method%fixed_step_mode) then
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
! don't adjust the step size
accept = .true.
me%last_accepted_step_size = dt ! save it [really only needs to be done once]
else
! evaluate error and compute new step size:
xerr = abs(xerr)
tol = me%rtol * abs(xf) + me%atol
call me%stepsize_method%compute_stepsize(me%n,dt,tol,xerr,p,dt_new,accept)
if (accept) me%last_accepted_step_size = dt ! save it
dt = dt_new
if (me%invalid) then
! Halve the timestep and try again...
me%invalid = .false.
accept = .false.
dt = dt * 0.5_wp
else
! evaluate error and compute new step size:
xerr = abs(xerr)
tol = me%rtol * abs(xf) + me%atol
call me%stepsize_method%compute_stepsize(me%n,dt,tol,xerr,p,dt_new,accept)
if (accept) me%last_accepted_step_size = dt ! save it
dt = dt_new
end if
end if

if (accept) then
Expand Down Expand Up @@ -1353,16 +1392,27 @@ subroutine integrate_to_event_variable_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
if (me%stopped) return

if (me%stepsize_method%fixed_step_mode) then
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
! don't adjust the step size
accept = .true.
me%last_accepted_step_size = dt ! save it [really only needs to be done once]
else
! evaluate error and compute new step size:
xerr = abs(xerr)
stol = me%rtol * abs(xf) + me%atol
call me%stepsize_method%compute_stepsize(me%n,dt,stol,xerr,p,dt_new,accept)
if (accept) me%last_accepted_step_size = dt ! save it
dt = dt_new
if (me%invalid) then
! Halve the timestep and try again...
me%invalid = .false.
accept = .false.
dt = dt * 0.5_wp
else
! evaluate error and compute new step size:
xerr = abs(xerr)
stol = me%rtol * abs(xf) + me%atol
call me%stepsize_method%compute_stepsize(me%n,dt,stol,xerr,p,dt_new,accept)
if (accept) me%last_accepted_step_size = dt ! save it
dt = dt_new
end if
end if

if (accept) then
Expand Down Expand Up @@ -1478,6 +1528,10 @@ function solver_func(delt) result(g)
! step that was already accepted, so it should be ok]
call me%step(t,x,delt,g_xf,xerr)
if (me%stopped) return
if (me%invalid) then
call me%raise_exception(RKLIB_ERROR_INVALID_F)
return
end if
call me%g(t+delt,g_xf,g)

end function solver_func
Expand Down