stepsize_class

type, public :: stepsize_class

Algorithms for adjusting the step size for variable-step Runge-Kutta integrators.

The new step size factor is computed from the scaled error e = norm(err/tol) (or norm(err/(tol*|h|)) if relative_err), where tol = rtol*s + atol, and s is set by rtol_scaling:

$h_{new} = h \cdot SF \cdot e^{-k} \cdot e_{old}^{\beta}$ , where $k = 1/(p + p_{offset}) - \beta_{mult} \beta$ ,

limited to [hfactor_reject, hfactor_accept]. The $e_{old}^{\beta}$ term is only used for accepted steps ( $\beta$ = pi_beta, which is 0 by default).

The defaults are the original rklib settings. To use the step size control of Hairer's DOPRI5 and DOP853 codes, use:

norm = rms_func, rtol_scaling = 2
p_exponent_offset = 0 (for these methods, $1/p = 1/(q+1)$ , where q is the order of the error estimate)
hfactor_accept = 6, hfactor_reject = 1/3 (approximately: 0.333), safety_factor = 0.9
no_increase_after_reject = .true.
pi_beta = 0.04, pi_beta_mult = 0.75 for DOPRI5 (rkdp54), or pi_beta = 0 for DOP853 (rkdp85, also with error_combination = 2)

Components

Type	Visibility	Attributes		Name		Initial
logical,	private		::	fixed_step_mode	=	.false.	if true, then the method runs in fixed step mode with not error estimation
real(kind=wp),	private		::	hmax	=	1.0e+6_wp	maximum allowed step size
real(kind=wp),	private		::	hmin	=	1.0e-6_wp	minimum allowed step size
real(kind=wp),	private		::	hfactor_reject	=	0.5_wp	minimum allowed factor for decreasing step size after rejected step
real(kind=wp),	private		::	hfactor_accept	=	2.0_wp	maximum allowed factor for increasing step size after accepted step
integer,	private		::	accept_mode	=	2	method to determine if step is accepted [1,2]
integer,	private		::	max_attempts	=	10000	maximum number of attempts to decrease step size before giving up
logical,	private		::	relative_err	=	.false.	if true, use the error per unit step: $\|\|err\|\| \le \|h\| tol$ . otherwise, use the error per step: $\|\|err\|\| \le tol$ . (The scaled error is then proportional to one less power of `h`, so `p_exponent_offset` should be one less than for the error per step.)
real(kind=wp),	private		::	safety_factor	=	0.9_wp	for `hfactor` equation (>0)
integer,	private		::	p_exponent_offset	=	1	the exponent of the `hfactor` equation is `1/(p + p_exponent_offset)` (`-1`, `0`, or `1`; `p + p_exponent_offset` is limited to be at least 1). It should be one over the power of `h` that the scaled error is proportional to. For a `p(p-1)` pair, the error estimate is proportional to $h^{q+1} = h^p$ (where `q=p-1` is the order of the error estimate), so `0` for the error per step, and `-1` for the error per unit step (`relative_err`).
procedure(norm_func),	private,	nopass, pointer	::	norm	=>	maxval_func	routine for computing the norm of the scaled error (e.g., maxval_func, norm2_func, rms_func, or a user function)
integer,	private		::	rtol_scaling	=	1	the state used to scale the relative tolerance: `1`: `tol = rtol\|x(t+h)\| + atol`, `2`: `tol = rtolmax(\|x(t)\|,\|x(t+h)\|) + atol` (as in `DOPRI5` and `DOP853`)
integer,	private		::	error_combination	=	1	for methods with two error estimates that are combined (currently `rkdp85`, which uses the 5th and 3rd order estimates of `DOP853`, see rkdp85): `1`: combined for each component, `2`: combined using the scaled norms of the two estimates (as in `DOP853`)
real(kind=wp),	private		::	pi_beta	=	0.0_wp	the $\beta$ parameter of the PI (Lund) step size control (stabilized step size control, see Hairer & Wanner, Solving ODEs I, section IV.2). `0` (the default) is the standard controller. Values up to `0.04` are recommended for `DOPRI5`.
real(kind=wp),	private		::	pi_beta_mult	=	0.75_wp	the error exponent is reduced by `pi_beta_mult*pi_beta` (`0.75` is used in `DOPRI5`, `0.2` in `DOP853`)
logical,	private		::	no_increase_after_reject	=	.false.	if true, the step size is not increased for the step after a rejected step (as in `DOPRI5` and `DOP853`)
real(kind=wp),	private		::	err_old	=	1.0e-4_wp	scaled error of the last accepted step (for the PI control)
logical,	private		::	last_rejected	=	.false.	if the last step was rejected

Type-Bound Procedures

procedure, public :: initialize => stepsize_class_constructor

private subroutine stepsize_class_constructor(me, hmin, hmax, hfactor_reject, hfactor_accept, norm, accept_mode, relative_err, safety_factor, p_exponent_offset, max_attempts, fixed_step_mode, rtol_scaling, error_combination, pi_beta, pi_beta_mult, no_increase_after_reject)

Constructor for a stepsize_class.

Arguments

Type	Intent	Optional		Name
class(stepsize_class),	intent(out)		::	me	(all options not specified are set to their defaults)
real(kind=wp),	intent(in),	optional	::	hmin	minimum allowed step size (>0)
real(kind=wp),	intent(in),	optional	::	hmax	maximum allowed step size (>0)
real(kind=wp),	intent(in),	optional	::	hfactor_reject	minimum allowed factor for decreasing step size after rejected step (>0)
real(kind=wp),	intent(in),	optional	::	hfactor_accept	maximum allowed factor for decreasing step size after accepted step (>0)
procedure(norm_func),		optional	::	norm	the user-specified $\|\|x\|\|$ function
integer,	intent(in),	optional	::	accept_mode	method to determine if step is accepted [1,2]
logical,	intent(in),	optional	::	relative_err	if true, use the error per unit step: $\|\|err\|\| \le \|h\| tol$ . Otherwise (the default), use the error per step: $\|\|err\|\| \le tol$ . Note that, for the error per unit step, the scaled error is proportional to one less power of `h`, so `p_exponent_offset` should be one less than for the error per step (e.g., `-1` rather than `0` for a `p(p-1)` pair).
real(kind=wp),	intent(in),	optional	::	safety_factor	for `hfactor` equation (>0)
integer,	intent(in),	optional	::	p_exponent_offset	the exponent of the step size factor is `1/(p + p_exponent_offset)` (`-1`, `0`, or `1`). See stepsize_class for details.
integer,	intent(in),	optional	::	max_attempts	max step size change attempts after rejected step
logical,	intent(in),	optional	::	fixed_step_mode	if true, then the method runs in fixed step mode with not error estimation. All the other inputs are ignored. Note that this requires a `dt /= 0` input for the integrator.
integer,	intent(in),	optional	::	rtol_scaling	the state used to scale the relative tolerance: `1`: `\|x(t+h)\|` (the default), `2`: `max(\|x(t)\|,\|x(t+h)\|)` (as in `DOPRI5` and `DOP853`).
integer,	intent(in),	optional	::	error_combination	for methods with two error estimates (`rkdp85`): `1`: combined for each component (the default), `2`: combined using their scaled norms (as in `DOP853`).
real(kind=wp),	intent(in),	optional	::	pi_beta	the $\beta$ parameter for the PI (Lund) stabilized step size control (`0` to disable, which is the default). Values up to `0.04` are recommended for `DOPRI5`.
real(kind=wp),	intent(in),	optional	::	pi_beta_mult	the error exponent is reduced by `pi_beta_mult*pi_beta` (`0.75` is the default, and is used in `DOPRI5`, `0.2` is used in `DOP853`).
logical,	intent(in),	optional	::	no_increase_after_reject	if true, the step size is not increased for the step after a rejected step (default is false).

procedure, public :: compute_stepsize

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

Compute the new step size using the specific method.

Arguments

Type	Intent	Attributes		Name
class(stepsize_class),	intent(inout)		::	me	(`inout` since it has some state, see `pi_beta` and `no_increase_after_reject`)
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

procedure, public :: destroy => destroy_stepsize_class

private subroutine destroy_stepsize_class(me)

Destructor for stepsize_class.

Arguments

Type	Intent	Optional	Attributes		Name
class(stepsize_class),	intent(out)			::	me

procedure, public :: reset => stepsize_class_reset

private pure subroutine stepsize_class_reset(me)

Reset the internal state of the step size control. This is called at the start of each integration.

Arguments

Type	Intent	Optional	Attributes		Name
class(stepsize_class),	intent(inout)			::	me

procedure, public :: tolerance => stepsize_class_tolerance

private pure function stepsize_class_tolerance(me, rtol, atol, x, xf) result(tol)

The error tolerance vector for a step from x to xf: rtol*s + atol, where s depends on the rtol_scaling option.

Arguments

Type	Intent	Attributes		Name
class(stepsize_class),	intent(in)		::	me
real(kind=wp),	intent(in),	dimension(:)	::	rtol	relative tolerance
real(kind=wp),	intent(in),	dimension(:)	::	atol	absolute tolerance
real(kind=wp),	intent(in),	dimension(:)	::	x	state at the start of the step
real(kind=wp),	intent(in),	dimension(:)	::	xf	state at the end of the step

Return Value real(kind=wp), dimension(size(x))

the error tolerance for each component

Source Code

    type,public :: stepsize_class

        !! Algorithms for adjusting the step size for variable-step
        !! Runge-Kutta integrators.
        !!
        !! The new step size factor is computed from the scaled error
        !! `e = norm(err/tol)` (or `norm(err/(tol*|h|))` if `relative_err`), where
        !! `tol = rtol*s + atol`, and `s` is set by `rtol_scaling`:
        !!
        !! \( h_{new} = h \cdot SF \cdot e^{-k} \cdot e_{old}^{\beta} \), where
        !! \( k = 1/(p + p_{offset}) - \beta_{mult} \beta \),
        !!
        !! limited to `[hfactor_reject, hfactor_accept]`. The \( e_{old}^{\beta} \) term is
        !! only used for accepted steps (\( \beta \) = `pi_beta`, which is 0 by default).
        !!
        !! The defaults are the original rklib settings. To use the step size control
        !! of Hairer's `DOPRI5` and `DOP853` codes, use:
        !!
        !!  * `norm = rms_func`, `rtol_scaling = 2`
        !!  * `p_exponent_offset = 0` (for these methods, \( 1/p = 1/(q+1) \),
        !!    where `q` is the order of the error estimate)
        !!  * `hfactor_accept = 6`, `hfactor_reject = 1/3` (approximately: `0.333`), `safety_factor = 0.9`
        !!  * `no_increase_after_reject = .true.`
        !!  * `pi_beta = 0.04`, `pi_beta_mult = 0.75` for `DOPRI5` (`rkdp54`), or
        !!    `pi_beta = 0` for `DOP853` (`rkdp85`, also with `error_combination = 2`)

        private

        logical :: fixed_step_mode = .false. !! if true, then the method runs in
                                             !! fixed step mode with not error estimation

        real(wp) :: hmax           = 1.0e+6_wp  !! maximum allowed step size
        real(wp) :: hmin           = 1.0e-6_wp  !! minimum allowed step size
        real(wp) :: hfactor_reject = 0.5_wp     !! minimum allowed factor for decreasing step size after rejected step
        real(wp) :: hfactor_accept = 2.0_wp     !! maximum allowed factor for increasing step size after accepted step
        integer  :: accept_mode    = 2          !! method to determine if step is accepted [1,2]
        integer  :: max_attempts   = 10000      !! maximum number of attempts to decrease step size before giving up

        ! for the `hfactor` equation:
        logical  :: relative_err      = .false. !! if true, use the error per unit step: \( ||err|| \le |h| tol \).
                                                !! otherwise, use the error per step: \( ||err|| \le tol \).
                                                !! (The scaled error is then proportional to one less power of `h`, so
                                                !! `p_exponent_offset` should be one less than for the error per step.)
        real(wp) :: safety_factor     = 0.9_wp  !! for `hfactor` equation (>0)
        integer  :: p_exponent_offset = 1       !! the exponent of the `hfactor` equation is `1/(p + p_exponent_offset)`
                                                !! (`-1`, `0`, or `1`; `p + p_exponent_offset` is limited to be at least 1).
                                                !! It should be one over the power of `h` that the scaled error is
                                                !! proportional to. For a `p(p-1)` pair, the error estimate is
                                                !! proportional to \( h^{q+1} = h^p \) (where `q=p-1` is the order of the
                                                !! error estimate), so `0` for the error per step, and `-1` for the
                                                !! error per unit step (`relative_err`).

        procedure(norm_func),nopass,pointer :: norm => maxval_func
            !! routine for computing the norm of the scaled error
            !! (e.g., [[maxval_func]], [[norm2_func]], [[rms_func]], or a user function)

        integer  :: rtol_scaling = 1 !! the state used to scale the relative tolerance:
                                     !! `1`: `tol = rtol*|x(t+h)| + atol`,
                                     !! `2`: `tol = rtol*max(|x(t)|,|x(t+h)|) + atol` (as in `DOPRI5` and `DOP853`)
        integer  :: error_combination = 1 !! for methods with two error estimates that are combined
                                          !! (currently `rkdp85`, which uses the 5th and 3rd order
                                          !! estimates of `DOP853`, see [[rkdp85]]):
                                          !! `1`: combined for each component,
                                          !! `2`: combined using the scaled norms of the two estimates (as in `DOP853`)
        real(wp) :: pi_beta = 0.0_wp      !! the \( \beta \) parameter of the PI (Lund) step size control
                                          !! (stabilized step size control, see Hairer & Wanner,
                                          !! Solving ODEs I, section IV.2). `0` (the default) is the standard
                                          !! controller. Values up to `0.04` are recommended for `DOPRI5`.
        real(wp) :: pi_beta_mult = 0.75_wp !! the error exponent is reduced by `pi_beta_mult*pi_beta`
                                           !! (`0.75` is used in `DOPRI5`, `0.2` in `DOP853`)
        logical  :: no_increase_after_reject = .false. !! if true, the step size is not increased
                                                       !! for the step after a rejected step (as in `DOPRI5` and `DOP853`)

        ! internal state (reset at the start of each integration, see [[stepsize_class_reset]]):
        real(wp) :: err_old = 1.0e-4_wp !! scaled error of the last accepted step (for the PI control)
        logical  :: last_rejected = .false. !! if the last step was rejected

        contains

        private

        procedure,public :: initialize => stepsize_class_constructor
        procedure,public :: compute_stepsize
        procedure,public :: destroy => destroy_stepsize_class
        procedure,public :: reset => stepsize_class_reset
        procedure,public :: tolerance => stepsize_class_tolerance

    end type stepsize_class

stepsize_class Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public :: stepsize_class

Inherited by

Components

Type-Bound Procedures

procedure, public :: initialize => stepsize_class_constructor

private subroutine stepsize_class_constructor(me, hmin, hmax, hfactor_reject, hfactor_accept, norm, accept_mode, relative_err, safety_factor, p_exponent_offset, max_attempts, fixed_step_mode, rtol_scaling, error_combination, pi_beta, pi_beta_mult, no_increase_after_reject)

Arguments

procedure, public :: compute_stepsize

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

Arguments

procedure, public :: destroy => destroy_stepsize_class

private subroutine destroy_stepsize_class(me)

Arguments

procedure, public :: reset => stepsize_class_reset

private pure subroutine stepsize_class_reset(me)

Arguments

procedure, public :: tolerance => stepsize_class_tolerance

private pure function stepsize_class_tolerance(me, rtol, atol, x, xf) result(tol)

Arguments

Return Value real(kind=wp), dimension(size(x))

Source Code