+
Skip to content

Specify root finder for event location #27

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

Merged
merged 1 commit into from
Jul 19, 2023
Merged

Specify root finder for event location #27

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
43 changes: 26 additions & 17 deletions src/rklib_module.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,7 @@ module rklib_module
procedure(deriv_func),pointer :: f => null() !! user-specified derivative function
procedure(report_func),pointer :: report => null() !! user-specified report function
procedure(event_func),pointer :: g => null() !! event function (stop when this is zero)
type(root_method) :: solver = root_method_brent !! the root solver method to use for even finding

real(wp),dimension(:,:),allocatable :: funcs !! matrix to store the function
!! evalutaions in the step function.
Expand Down Expand Up @@ -507,7 +508,8 @@ end subroutine set_fsal_cache
! Initialize the [[rk_class]].

subroutine initialize_rk_class(me,n,f,report,g,stop_on_errors,&
max_number_of_steps,report_rate)
max_number_of_steps,report_rate,&
solver)

implicit none

Expand All @@ -524,6 +526,8 @@ subroutine initialize_rk_class(me,n,f,report,g,stop_on_errors,&
!! `1` : report every point,
!! `2` : report every other point, etc.
!! The first and last point are always reported.
class(root_method),intent(in),optional :: solver !! the root-finding method to use for even finding.
!! if not present, then `brent_solver` is used.

type(rklib_properties) :: props !! to get the method properties

Expand All @@ -536,6 +540,7 @@ subroutine initialize_rk_class(me,n,f,report,g,stop_on_errors,&
if (present(stop_on_errors)) me%stop_on_errors = stop_on_errors
if (present(max_number_of_steps)) me%max_number_of_steps = abs(max_number_of_steps)
if (present(report_rate)) me%report_rate = abs(report_rate)
if (present(solver)) me%solver = solver

! allocate the registers:
props = me%properties()
Expand Down Expand Up @@ -577,7 +582,8 @@ end subroutine begin_integration_rk_fixed_step_class
! Initialize the [[rk_fixed_step_class]].

subroutine initialize_fixed_step(me,n,f,report,g,stop_on_errors,&
max_number_of_steps,report_rate)
max_number_of_steps,report_rate,&
solver)

implicit none

Expand All @@ -594,9 +600,11 @@ subroutine initialize_fixed_step(me,n,f,report,g,stop_on_errors,&
!! `1` : report every point,
!! `2` : report every other point, etc.
!! The first and last point are always reported.
class(root_method),intent(in),optional :: solver !! the root-finding method to use for even finding.
!! if not present, then `brent_solver` is used.

! base init all we need here:
call me%init(n,f,report,g,stop_on_errors,max_number_of_steps,report_rate)
call me%init(n,f,report,g,stop_on_errors,max_number_of_steps,report_rate,solver)

end subroutine initialize_fixed_step
!*****************************************************************************************
Expand Down Expand Up @@ -695,7 +703,6 @@ subroutine integrate_to_event_fixed_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
real(wp),dimension(me%n) :: g_xf !! state vector from the root finder
logical :: first !! it is the first step
logical :: last !! it is the last step
type(brent_solver) :: solver !! for the root finding problem
integer :: iflag !! return flag from `solver`

if (.not. associated(me%f)) then
Expand Down Expand Up @@ -767,11 +774,12 @@ subroutine integrate_to_event_fixed_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt]

!find the root:
call solver%initialize(solver_func)
call solver%solve(zero,dt,dt_root,dum,iflag,fax=ga,fbx=gb)
call root_scalar(me%solver,solver_func,zero,dt,dt_root,dum,iflag,&
fax=ga, fbx=gb, rtol=tol, atol=tol)
! ftol,maxiter,bisect_on_failure) ! other options
if (me%stopped) return
t2 = t + dt_root
gf = solver_func(solver,dt_root)
gf = solver_func(dt_root)
if (me%stopped) return
tf = t2
xf = g_xf !computed in the solver function
Expand Down Expand Up @@ -801,13 +809,12 @@ subroutine integrate_to_event_fixed_step(me,t0,x0,h,tmax,tol,tf,xf,gf)

contains

function solver_func(this,delt) result(g)
function solver_func(delt) result(g)

!! root solver function. The input is the `dt` offset from time `t`.

implicit none

class(root_solver),intent(inout) :: this
real(wp),intent(in) :: delt !! from [0 to `dt`]
real(wp) :: g

Expand Down Expand Up @@ -1022,7 +1029,8 @@ end subroutine begin_integration_rk_variable_step_class

subroutine initialize_variable_step(me,n,f,rtol,atol,stepsize_method,&
hinit_method,report,g,stop_on_errors,&
max_number_of_steps,report_rate)
max_number_of_steps,report_rate,&
solver)

implicit none

Expand Down Expand Up @@ -1052,11 +1060,13 @@ subroutine initialize_variable_step(me,n,f,rtol,atol,stepsize_method,&
!! `1` : report every point,
!! `2` : report every other point, etc.
!! The first and last point are always reported.
class(root_method),intent(in),optional :: solver !! the root-finding method to use for even finding.
!! if not present, then `brent_solver` is used.

real(wp),parameter :: default_tol = 100*epsilon(1.0_wp) !! if tols not specified

! base init:
call me%init(n,f,report,g,stop_on_errors,max_number_of_steps,report_rate)
call me%init(n,f,report,g,stop_on_errors,max_number_of_steps,report_rate,solver)

! variable-step specific inputs:
if (allocated(me%rtol)) deallocate(me%rtol)
Expand Down Expand Up @@ -1299,7 +1309,6 @@ subroutine integrate_to_event_variable_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
integer :: i,p,iflag
real(wp) :: t,dt,t2,ga,gb,dt_root,dum,dt_new
logical :: first,last,accept
type(brent_solver) :: solver

if (.not. associated(me%f)) then
call me%raise_exception(RKLIB_ERROR_F_NOT_ASSOCIATED)
Expand Down Expand Up @@ -1416,11 +1425,12 @@ subroutine integrate_to_event_variable_step(me,t0,x0,h,tmax,tol,tf,xf,gf)
elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt]

!find the root:
call solver%initialize(solver_func, rtol=tol, atol=tol)
call solver%solve(zero,dt,dt_root,dum,iflag,fax=ga,fbx=gb)
call root_scalar(me%solver,solver_func,zero,dt,dt_root,dum,iflag,&
fax=ga, fbx=gb, rtol=tol, atol=tol)
! ftol,maxiter,bisect_on_failure) ! other options
if (me%stopped) return
t2 = t + dt_root
gf = solver_func(solver,dt_root)
gf = solver_func(dt_root)
if (me%stopped) return
tf = t2
xf = g_xf !computed in the solver function
Expand Down Expand Up @@ -1452,13 +1462,12 @@ subroutine integrate_to_event_variable_step(me,t0,x0,h,tmax,tol,tf,xf,gf)

contains

function solver_func(this,delt) result(g)
function solver_func(delt) result(g)

!! root solver function. The input is the `dt` offset from time `t`.

implicit none

class(root_solver),intent(inout) :: this
real(wp),intent(in) :: delt !! from [0 to `dt`]
real(wp) :: g

Expand Down
4 changes: 3 additions & 1 deletion test/rk_test_variable_step.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ program rk_test_variable_step
use test_support
use pyplot_module
use iso_fortran_env
use root_module, only: root_method_brent

implicit none

Expand Down Expand Up @@ -428,7 +429,8 @@ subroutine run_test(name)

call s2%initialize(n=n,f=twobody,g=twobody_event,&
rtol=[1.0e-12_wp],atol=[1.0e-12_wp],&
stepsize_method=sz)
stepsize_method=sz,&
solver = root_method_brent) ! specify root solver

fevals = 0
first = .true.
Expand Down
点击 这是indexloc提供的php浏览器服务，不要输入任何密码和下载