integrate_to_event_variable_step Subroutine

private subroutine integrate_to_event_variable_step(me, t0, x0, h, tmax, tol, tf, xf, gf)

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 Bound

rk_variable_step_class

Arguments

Type IntentOptional 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

Calls

proc~~integrate_to_event_variable_step~~CallsGraph proc~integrate_to_event_variable_step rklib_module::rk_variable_step_class%integrate_to_event_variable_step g g proc~integrate_to_event_variable_step->g proc~begin_integration_rk_variable_step_class rklib_module::rk_variable_step_class%begin_integration_rk_variable_step_class proc~integrate_to_event_variable_step->proc~begin_integration_rk_variable_step_class proc~compute_initial_step rklib_module::rk_variable_step_class%compute_initial_step proc~integrate_to_event_variable_step->proc~compute_initial_step proc~compute_stepsize rklib_module::stepsize_class%compute_stepsize proc~integrate_to_event_variable_step->proc~compute_stepsize proc~export_point rklib_module::rk_class%export_point proc~integrate_to_event_variable_step->proc~export_point proc~order rklib_module::rk_variable_step_class%order proc~integrate_to_event_variable_step->proc~order proc~raise_exception rklib_module::rk_class%raise_exception proc~integrate_to_event_variable_step->proc~raise_exception root_scalar root_scalar proc~integrate_to_event_variable_step->root_scalar step step proc~integrate_to_event_variable_step->step begin begin proc~begin_integration_rk_variable_step_class->begin proc~destroy_fsal_cache rklib_module::rk_variable_step_fsal_class%destroy_fsal_cache proc~begin_integration_rk_variable_step_class->proc~destroy_fsal_cache f f proc~compute_initial_step->f proc~hinit rklib_module::rk_variable_step_class%hinit proc~compute_initial_step->proc~hinit proc~hstart rklib_module::rk_variable_step_class%hstart proc~compute_initial_step->proc~hstart properties properties proc~order->properties proc~hinit->proc~order proc~hinit->f proc~hstart->proc~order proc~hstart->f

Source Code

    subroutine integrate_to_event_variable_step(me,t0,x0,h,tmax,tol,tf,xf,gf)

    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

    real(wp),dimension(me%n) :: x,g_xf
    real(wp),dimension(me%n) :: xerr !! truncation error estimate
    real(wp),dimension(me%n) :: stol
    integer :: i,p,iflag
    real(wp) :: t,dt,t2,ga,gb,dt_root,dum,dt_new
    logical :: first,last,accept

    if (.not. associated(me%f)) then
        call me%raise_exception(RKLIB_ERROR_F_NOT_ASSOCIATED)
        return
    end if
    if (.not. associated(me%g)) then
        call me%raise_exception(RKLIB_ERROR_G_NOT_ASSOCIATED)
        return
    end if

    call me%begin_integration()

    call me%export_point(t0,x0,.true.)  !first point

    if (abs(t0-tmax)<=zero) 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
        dt = me%compute_initial_step(t0,tmax,x0,h)
        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,xerr)
                if (me%stopped) return

                if (me%stepsize_method%fixed_step_mode) then
                    ! 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
                end if

                if (accept) then
                    !accept this step
                    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)
                        return
                    end if
                    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
                        call me%raise_exception(RKLIB_ERROR_TOO_MANY_REDUCTIONS)
                        return
                    end if
                    if (abs(dt) < abs(me%stepsize_method%hmin)) then
                        call me%raise_exception(RKLIB_ERROR_MIN_STEP_SIZE)
                        return
                    end if

                    last = ((dt>=zero .and. (t+dt)>=tmax) .or. &  !adjust last time step
                            (dt<zero .and. (t+dt)<=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
                    call me%export_point(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 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(dt_root)
                if (me%stopped) return
                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
                call me%export_point(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

    call me%export_point(tf,xf,.true.)   !last point

    contains

        function solver_func(delt) result(g)

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

        implicit none

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

        real(wp),dimension(me%n) :: xerr !! 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,xerr)
        if (me%stopped) return
        call me%g(t+delt,g_xf,g)

        end function solver_func

    end subroutine integrate_to_event_variable_step