integrate_variable_step Subroutine

private subroutine integrate_variable_step(me, t0, x0, h, tf, xf)

Main integration routine for the rk_variable_step_class.

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(:) :: x0

initial state
real(kind=wp), intent(in) :: h

initial abs(time step)
real(kind=wp), intent(in) :: tf

final time
real(kind=wp), intent(out), dimension(:) :: xf

final state

Calls

proc~~integrate_variable_step~~CallsGraph proc~integrate_variable_step rklib_module::rk_variable_step_class%integrate_variable_step proc~begin_integration_rk_variable_step_class rklib_module::rk_variable_step_class%begin_integration_rk_variable_step_class proc~integrate_variable_step->proc~begin_integration_rk_variable_step_class proc~compute_initial_step rklib_module::rk_variable_step_class%compute_initial_step proc~integrate_variable_step->proc~compute_initial_step proc~compute_stepsize rklib_module::stepsize_class%compute_stepsize proc~integrate_variable_step->proc~compute_stepsize proc~export_point rklib_module::rk_class%export_point proc~integrate_variable_step->proc~export_point proc~order rklib_module::rk_variable_step_class%order proc~integrate_variable_step->proc~order proc~raise_exception rklib_module::rk_class%raise_exception proc~integrate_variable_step->proc~raise_exception step step proc~integrate_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

Called by

proc~~integrate_variable_step~~CalledByGraph proc~integrate_variable_step rklib_module::rk_variable_step_class%integrate_variable_step program~rklib_example rklib_example program~rklib_example->proc~integrate_variable_step

Source Code

    subroutine integrate_variable_step(me,t0,x0,h,tf,xf)

    implicit none

    class(rk_variable_step_class),intent(inout)     :: me
    real(wp),intent(in)               :: t0    !! initial time
    real(wp),dimension(:),intent(in)  :: x0    !! initial state
    real(wp),intent(in)               :: h     !! initial abs(time step)
    real(wp),intent(in)               :: tf    !! final time
    real(wp),dimension(:),intent(out) :: xf    !! final state

    real(wp) :: t,dt,t2,dt_new
    real(wp),dimension(me%n) :: x,xerr,tol
    logical :: last !! it is the last step
    logical :: accept !! the step is accepted
    integer :: i !! max step size reduction attempts counter
    integer :: p !! order of the method

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

    call me%begin_integration()

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

    if (abs(t0-tf)<=zero) then
        xf = x0
    else

        t = t0
        x = x0
        dt = me%compute_initial_step(t0,tf,x0,h)
        p = me%order()     !order of the method

        do
            t2 = t + dt
            last = ((dt>=zero .and. t2>=tf) .or. &  !adjust last time step
                    (dt<zero .and. t2<=tf))         !
            if (last) dt = tf-t                     !

            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)
                    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

                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)>=tf) .or. &  !adjust last time step
                            (dt<zero .and. (t+dt)<=tf))         !
                    if (last) dt = tf-t                         !
                    t2 = t + dt

                end if

            end do

            if (last) exit
            call me%export_point(t2,xf)   !intermediate point
            x = xf
            t = t2
        end do

    end if

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

    end subroutine integrate_variable_step