compute_stepsize Subroutine

private pure subroutine compute_stepsize(me, 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
real(kind=wp), intent(in) :: h

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

abs error tolerance (>0)
real(kind=wp), intent(in) :: 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 rk_module_variable_step::stepsize_class%compute_stepsize proc~integrate rk_module_variable_step::rk_variable_step_class%integrate proc~integrate->proc~compute_stepsize proc~integrate_to_event rk_module_variable_step::rk_variable_step_class%integrate_to_event proc~integrate_to_event->proc~compute_stepsize proc~step_size_test rk_module_variable_step::step_size_test proc~step_size_test->proc~compute_stepsize proc~rk_test_variable_step rk_module_variable_step::rk_test_variable_step proc~rk_test_variable_step->proc~integrate proc~rk_test_variable_step->proc~integrate_to_event

Source Code

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

    !! Compute the new step size using the specific method.

    implicit none

    class(stepsize_class),intent(in) :: me
    real(wp),intent(in)              :: h      !! current step size (<>0)
    real(wp),intent(in)              :: tol    !! abs error tolerance (>0)
    real(wp),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) :: hfactor  !! step size factor (>0)
    real(wp),parameter :: small = ten * epsilon(one) !! small error value

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

        hfactor = me%hfactor_accept
        accept = .true.

    else

        ! compute base factor based on the selected formula:
        !hfactor = abs(me%compute_h_factor(h,tol,err,p))
        if (me%relative_err) then
            hfactor = abs( me%safety_factor*abs(tol*h/err)**(one/real(p+me%p_exponent_offset,wp)) )
        else
            hfactor = abs( me%safety_factor*abs(tol/err)**(one/real(p+me%p_exponent_offset,wp)) )
        end if

        ! if the step is to be accepted:
        select case (me%accept_mode)
        case(1) !algorithm 17.12
            accept = (hfactor>=one)
        case(2) !algorithm 17.13
            accept = (err<=tol)
        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 subroutine compute_stepsize