compute_stepsize Subroutine

private subroutine compute_stepsize(me, n, h, tol, err, p, hnew, accept)

Compute the new step size using the specific method.

Type Bound

stepsize_class

Arguments

Type IntentOptional Attributes Name
class(stepsize_class), intent(in) :: me
integer, intent(in) :: n

number of variables
real(kind=wp), intent(in) :: h

current step size (<>0)
real(kind=wp), intent(in), dimension(n) :: tol

abs error tolerance (>0)
real(kind=wp), intent(in), dimension(n) :: err

truncation error estimate (>0)
integer, intent(in) :: p

order of the method
real(kind=wp), intent(out) :: hnew

new step size (<>0)
logical, intent(out) :: accept

if the step is accepted

Called by

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

Source Code

    subroutine compute_stepsize(me,n,h,tol,err,p,hnew,accept)

    implicit none

    class(stepsize_class),intent(in) :: me
    integer,intent(in)               :: n      !! number of variables
    real(wp),intent(in)              :: h      !! current step size (<>0)
    real(wp),dimension(n),intent(in) :: tol    !! abs error tolerance (>0)
    real(wp),dimension(n),intent(in) :: err    !! truncation error estimate (>0)
    integer,intent(in)               :: p      !! order of the method
    real(wp),intent(out)             :: hnew   !! new step size (<>0)
    logical,intent(out)              :: accept !! if the step is accepted

    real(wp) :: e        !! exponent
    real(wp) :: hfactor  !! step size factor (>0)
    real(wp) :: max_err  !! max error for all the elements

    real(wp),parameter :: small = 10.0_wp * epsilon(1.0_wp) !! small error value

    if (me%fixed_step_mode) then
        ! do not adjust the step size
        accept = .true.
        hnew = h
    else

        if (all(err<=small)) then ! the error is extremely small

            hfactor = me%hfactor_accept
            accept = .true.

        else

            ! compute base factor based on the selected formula:
            if (me%relative_err) then
                max_err = me%norm(err/tol*abs(h))
            else
                max_err = me%norm(err/tol)
            end if

            e = 1.0_wp / real(p+me%p_exponent_offset,wp)
            hfactor = abs( me%safety_factor * (1.0_wp/max_err)**e )

            ! if the step is to be accepted:
            select case (me%accept_mode)
            case(1) !algorithm 17.12
                accept = hfactor >= 1.0_wp
            case(2) !algorithm 17.13
                accept = me%norm(err/tol) <= 1.0_wp
            end select

            !...notes:
            ! see: L. Shampine "Some Practical Runge-Kutta Formulas",
            !      Mathematics of Computation, 46(173), Jan 1986.
            ! different conditions for satisfying error conditions:
            !  ||err|| <= tol   -- Error per step (EPS)
            !  ||err|| <= h*tol -- Error per unit step (EPUS)

            !compute the actual hfactor based on the limits:
            if (accept) then
                hfactor = min(me%hfactor_accept, hfactor)
            else
                hfactor = max(me%hfactor_reject, hfactor)
            end if

        end if

        ! compute the new step size (enforce min/max bounds & add sign):
        hnew = sign(max(me%hmin,min(me%hmax,abs(h)*hfactor)),h)

    end if

    end subroutine compute_stepsize