Event-finding integration routine for the rk_variable_step_class. Integrates until g(t,x)=0, or until t=tf (whichever happens first).
Note
There are some efficiency improvements that could be made here. This is a work in progress.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(rk_variable_step_class), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(in) | :: | t0 |
initial time |
||
| real(kind=wp), | intent(in), | dimension(me%n) | :: | x0 |
initial state |
|
| real(kind=wp), | intent(in) | :: | h |
abs(time step) |
||
| real(kind=wp), | intent(in) | :: | tmax |
max final time if event not located |
||
| real(kind=wp), | intent(in) | :: | tol |
function tolerance for root finding |
||
| real(kind=wp), | intent(out) | :: | tf |
actual final time reached |
||
| real(kind=wp), | intent(out), | dimension(me%n) | :: | xf |
final state (at tf) |
|
| real(kind=wp), | intent(out) | :: | gf |
g value at tf |
||
| integer, | intent(out), | optional | :: | ierr |
0 = no errors, <0 = error. if not present, an error will stop program. |
subroutine integrate_to_event(me,t0,x0,h,tmax,tol,tf,xf,gf,ierr) use brent_module implicit none class(rk_variable_step_class),intent(inout) :: me real(wp),intent(in) :: t0 !! initial time real(wp),dimension(me%n),intent(in) :: x0 !! initial state real(wp),intent(in) :: h !! abs(time step) real(wp),intent(in) :: tmax !! max final time if event not located real(wp),intent(in) :: tol !! function tolerance for root finding real(wp),intent(out) :: tf !! actual final time reached real(wp),dimension(me%n),intent(out) :: xf !! final state (at tf) real(wp),intent(out) :: gf !! g value at tf integer,intent(out),optional :: ierr !! 0 = no errors, !! <0 = error. !! if not present, an error will stop program. real(wp),dimension(me%n) :: etol,xp0 real(wp),dimension(me%n) :: x,g_xf real(wp),dimension(me%n) :: terr !! truncation error estimate integer :: i,p,iflag real(wp) :: t,dt,t2,ga,gb,dt_root,dum,err,dt_new,stol logical :: first,last,export,accept procedure(report_func),pointer :: report type(brent_class) :: solver if (present(ierr)) ierr = 0 if (.not. associated(me%f)) then if (present(ierr)) then ierr = -1 return else error stop 'Error in integrate_to_event: f is not associated.' end if end if if (.not. associated(me%g)) then if (present(ierr)) then ierr = -2 return else error stop 'Error in integrate_to_event: g is not associated.' end if end if me%num_rejected_steps = 0 export = associated(me%report) if (export) call me%report(t0,x0) !first point if (t0==tmax) then xf = x0 tf = t0 call me%g(t0,x0,gf) else first = .true. t = t0 x = x0 call me%g(t,x,ga) !evaluate event function if (h==zero) then ! compute an appropriate initial step size: ! WARNING: this may not be working in all cases ..... etol = me%rtol * me%stepsize_method%norm(x0) + me%atol call me%f(t0,x0,xp0) ! get initial dx/dt select case (me%hinit_method) case(1) call me%hstart(t0,tmax,x0,xp0,etol,dt) case(2) dt = me%hinit(t0,x0,sign(one,tmax-t0),xp0,me%stepsize_method%hmax,me%atol,me%rtol) case default if (present(ierr)) then ierr = -3 return else error stop 'invalid hinit_method selection' end if end select else ! user-specified initial step size: dt = sign(h,tmax-t0) ! (correct sign) end if p = me%order() !order of the method do t2 = t + dt last = ((dt>=zero .and. t2>=tmax) .or. & !adjust last time step (dt<zero .and. t2<=tmax)) ! if (last) then dt = tmax-t t2 = tmax end if do i=0,me%stepsize_method%max_attempts ! take a step: call me%step(t,x,dt,xf,terr) ! evaluate error and compute new step size: err = me%stepsize_method%norm(terr) stol = me%stepsize_method%norm( me%rtol * xf + me%atol ) call me%stepsize_method%compute_stepsize(dt,stol,err,p,dt_new,accept) dt = dt_new if (accept) then !accept this step exit else !step is rejected, repeat step with new dt me%num_rejected_steps = me%num_rejected_steps + 1 !note: if we have reached the min step size, and the error !is still too large, we can't proceed. if (i>=me%stepsize_method%max_attempts) then if (present(ierr)) then ierr = -4 return else error stop 'error: too many attempts to reduce step size.' end if end if if (abs(dt) <= abs(me%stepsize_method%hmin)) then if (present(ierr)) then ierr = -5 return else error stop 'warning: min step size.' end if end if !...... !... if we have two rejected steps and the step size hasn't changed.. ! then we need to abort, since no progress is being made... !...... last = ((dt>=zero .and. t2>=tmax) .or. & !adjust last time step (dt<zero .and. t2<=tmax)) ! if (last) then dt = tmax-t t2 = tmax else t2 = t + dt end if end if end do call me%g(t2,xf,gb) !evaluate event function if (first .and. abs(ga)<=tol) then !we ignore a root at t0 after the first step if (abs(gb)<=tol) then !check this one since it could have landed on a root gf = gb tf = t2 exit else if (last) then !exiting without having found a root tf = t2 gf = gb exit end if if (export) call me%report(t2,xf) !intermediate point x = xf t = t2 ga = gb end if elseif (ga*gb<=zero) then !there is a root somewhere on [t,t+dt] !find the root: call solver%set_function(solver_func) call solver%find_zero(zero,dt,tol,dt_root,dum,iflag,ga,gb) t2 = t + dt_root gf = solver_func(solver,dt_root) tf = t2 xf = g_xf !computed in the solver function exit else !no root yet, continue if (last) then !exiting without having found a root tf = t2 gf = gb exit end if if (export) call me%report(t2,xf) !intermediate point x = xf t = t2 ga = gb end if if (first) first = .false. if (last) exit x = xf t = t2 end do end if if (export) call me%report(tf,xf) !last point contains function solver_func(this,delt) result(g) !! root solver function. The input is the dt offset from time t. implicit none class(brent_class),intent(inout) :: this real(wp),intent(in) :: delt !! from [0 to dt] real(wp) :: g real(wp),dimension(me%n) :: terr !! truncation error estimate !take a step from t to t+delt and evaluate g function: ! [we don't check the error because we are within a ! step that was already accepted, so it should be ok] call me%step(t,x,delt,g_xf,terr) call me%g(t+delt,g_xf,g) end function solver_func end subroutine integrate_to_event