diff --git a/src/rklib_module.F90 b/src/rklib_module.F90 index 3b40cd54..b178c695 100644 --- a/src/rklib_module.F90 +++ b/src/rklib_module.F90 @@ -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. @@ -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 @@ -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 @@ -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() @@ -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 @@ -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 !***************************************************************************************** @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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) @@ -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 @@ -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 diff --git a/test/rk_test_variable_step.f90 b/test/rk_test_variable_step.f90 index 6e073cca..945807b4 100644 --- a/test/rk_test_variable_step.f90 +++ b/test/rk_test_variable_step.f90 @@ -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 @@ -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.