foodie_integrator_runge_kutta_emd Module

  1. foodie_integrator_runge_kutta_embedded.f90
  2. foodie_integrator_runge_kutta_emd
module~~foodie_integrator_runge_kutta_emd~~UsesGraph module~foodie_integrator_runge_kutta_emd foodie_integrator_runge_kutta_emd module~foodie_integrator_multistage_object foodie_integrator_multistage_object module~foodie_integrator_multistage_object->module~foodie_integrator_runge_kutta_emd penf penf penf->module~foodie_integrator_runge_kutta_emd penf->module~foodie_integrator_multistage_object module~foodie_integrator_object foodie_integrator_object penf->module~foodie_integrator_object module~foodie_error_codes foodie_error_codes penf->module~foodie_error_codes module~foodie_integrand_object foodie_integrand_object penf->module~foodie_integrand_object module~foodie_integrator_object->module~foodie_integrator_runge_kutta_emd module~foodie_integrator_object->module~foodie_integrator_multistage_object module~foodie_error_codes->module~foodie_integrator_runge_kutta_emd module~foodie_integrand_object->module~foodie_integrator_runge_kutta_emd module~foodie_integrand_object->module~foodie_integrator_multistage_object iso_fortran_env iso_fortran_env iso_fortran_env->module~foodie_integrator_multistage_object iso_fortran_env->module~foodie_integrator_object
Help

FOODIE integrator: provide an explicit class of embedded Runge-Kutta schemes, from 2nd to 10th order accurate.

The integrators provided have the embedded pairs property allowing for automatic step size control. The schemes are explicit and defined through the extended Butcher's table syntax, see[1] .

Considering the following ODE system:

$$ U_t = R(t,U) $$

where \(U_t = \frac{dU}{dt}\), U is the vector of state variables being a function of the time-like independent variable t, R is the (vectorial) residual function, the class of schemes implemented are written in the form:

$$ U_p^{n+1} = U^n +\Delta t \sum_{s=1}^{Ns}\beta_p^s K^s $$ $$ U_{p+1}^{n+1} = U^n +\Delta t \sum_{s=1}^{Ns}\beta_{p+1}^s K^s $$

p is the lower accuracy order scheme and p+1 is the higher one; Ns is the number of stages used and \(K^s\) is the \(s^{th}\) stage computed as:

$$ K^s = R\left( t^n+\gamma^s \Delta t, U^n +\Delta t \sum_{i=1}^{s-1}\alpha^{s,i} K^i \right) $$

The schemes are explicit thus the above summation is up to \(s-1\). The coefficients \(\beta\), \(\alpha\) and \(\gamma\) are given in the extended Butcher table form:

  gamma^1    | alpha^{1,1}       alpha^{1,2}       ...        alpha^{1,Ns}
  gamma^2    | alpha^{2,1}       alpha^{2,2}       ...        alpha^{2,Ns}
  .          | .                 .                 .          .
  .          | .                 .                  .         .
  .          | .                 .                   .        .
  gamma^{Ns} | alpha^{Ns,1}      alpha^{Ns,2}      ...        alpha^{Ns,Ns}
 ------------|-------------------------------------------------------------
             | beta_{p+1}^1      beta_{p+1}^2      ...        beta_{p+1}^{Ns}
             | beta_p^1          beta_p^2          ...        beta_p^{Ns}

Because only explicit schemes are considered the Butcher table reduces to diagonal matrix:

  gamma^1    | 0                 0                 ...        0
  gamma^2    | alpha^{2,1}       0                 ...        0
  .          | .                 .                 .          .
  .          | .                 .                  .         .
  .          | .                 .                   .        .
  gamma^{Ns} | alpha^{Ns,1}      alpha^{Ns,2}      ...        0
 ------------|-------------------------------------------------------------
             | beta_{p+1}^1      beta_{p+1}^2      ...        beta_{p+1}^{Ns}
             | beta_p^1          beta_p^2          ...        beta_p^{Ns}

Moreover the following relation always holds: \( \gamma^s = \sum_{i=1}^{Ns}\alpha^{s,i} \)

The different schemes are selected accordingly to the number of stages used. Currently the following schemes are available:

2 stages, 2th order

This scheme is due to Heun-Euler.

  0  | 0
  1  | 1     0
 ----------------
     | 1/2   1/2
     | 1      0
6 stages, 5th order

This scheme is due to Cash and Karp, see [3].

  0    | 0
  1/5    | 1/5
  3/10 | 3/40            9/40
  3/5    | 3/10          -9/10        6/5
  1    | -11/54        5/2          -70/27      35/27
  7/8    | 1631/55296    175/512      575/13824   44275/110592     253/4096     0
 ----------------------------------------------------------------------------------------
       | 37/378        0           250/621      125/594          0            512/1771
       | 2825/27648    0           18575/48384  13525/55296      277/14336    1/4
7 stages, 4th order

This scheme is due to Dormand and Prince, see [1].

  0    | 0
  1/5  | 1/5
  3/10 | 3/40          9/40
  4/5  | 44/45        -56/15        32/9
  8/9  | 19372/6561   -25360/2187   64448/6561   -212/729
  1    | 9017/3168    -355/33       46732/5247    49/176      -5103/18656
  1    | 35/384        0            500/1113      125/192     -2187/6784      11/84      0
 --------------------------------------------------------------------------------------------
       | 5179/57600    0            7571/16695    393/640     -92097/339200   187/2100   1/40
       | 35/384        0            500/1113      125/192     -2187/6784      11/84      0
9 stages, 6th order

This scheme is due to Calvo et al., see [2].

  0                 | 0
  2/15              | 2/15
  1/5               | 1/20                  3/20
  3/10              | 3/40                  0                      9/40
  14/25             | 86727015/196851553    -60129073/52624712     957436434/1378352377    83886832/147842441
  19/25             | -86860849/45628967    111022885/25716487     108046682/101167669     -141756746/36005461
  35226607/35688279 | 77759591/16096467     -49252809/6452555      -381680111/51572984     879269579/66788831
  1                 | 237564263/39280295    -100523239/10677940    -265574846/27330247     317978411/18988713
  1                 | 17572349/289262523    0                      57513011/201864250      15587306/354501571
 --------------------------------------------------------------------------------------------------------------
                    | 17572349/289262523    0                      57513011/201864250      15587306/354501571
                    | 15231665/510830334    0                      59452991/116050448      -28398517/122437738
 ...continued...
  0                 |
  2/15              |
  1/5               |
  3/10              |
  14/25             |
  19/25             | 73139862/60170633
  35226607/35688279 | -90453121/33722162     111179552/157155827
  1                 | -124494385/35453627    86822444/100138635     -12873523/724232625
  1                 | 71783021/234982865     29672000/180480167     65567621/127060952     -79074570/210557597    0
 -----------------------------------------------------------------------------------------------------------------------
                    | 71783021/234982865     29672000/180480167     65567621/127060952     -79074570/210557597    0
                    | 56673824/137010559     68003849/426673583     7097631/37564021       -71226429/583093742    1/20
17 stages, 10th order

This scheme is due to Feagin, see [4].

  0                        |  0
  0.1                      |  0.1
  0.539357840802981787532  | -0.915176561375291440520  1.454534402178273228052
  0.809036761204472681298  |  0.202259190301118170324  0                        0.606777570903354510974
  0.309036761204472681298  |  0.184024714708643575149  0                        0.197966831227192369068
  0.981074190219795268254  |  0.087900734020668133731  0                        0
  0.833333333333333333333  |  0.085970050490246030218  0                        0
  0.354017365856802376329  |  0.120930449125333720660  0                        0
  0.882527661964732346425  |  0.110854379580391483508  0                        0
  0.642615758240322548157  |  0.112054414752879004829  0                        0
  0.357384241759677451842  |  0.113976783964185986138  0                        0
  0.117472338035267653574  |  0.079831452828019604635  0                        0
  0.833333333333333333333  |  0.985115610164857280120  0                        0
  0.309036761204472681298  |  0.895080295771632891049  0                        0.197966831227192369068
  0.539357840802981787532  | -0.915176561375291440520  1.454534402178273228052  0
  0.1                      |  0.1                      0                       -0.157178665799771163367
  1                        |  0.181781300700095283888  0.675                    0.342758159847189839942
 ------------------------------------------------------------------------------------------------------
                           |  0.033333333333333333333  0.025                    0.033333333333333333333
 ...continued...
  0                        |
  0.1                      |
  0.539357840802981787532  |
  0.809036761204472681298  |
  0.309036761204472681298  | -0.072954784731363262918
  0.981074190219795268254  |  0.410459702520260645318  0.482713753678866489204
  0.833333333333333333333  |  0.330885963040722183948  0.489662957309450192844 -0.073185637507085073678
  0.354017365856802376329  |  0                        0.260124675758295622809  0.032540262154909133015
  0.882527661964732346425  |  0                        0                       -0.060576148825500558762
  0.642615758240322548157  |  0                        0                       -0.144942775902865915672
  0.357384241759677451842  |  0                        0                       -0.076881336420335693858
  0.117472338035267653574  |  0                        0                       -0.052032968680060307651
  0.833333333333333333333  |  0.330885963040722183948  0.489662957309450192844 -1.378964865748435675821
  0.309036761204472681298  | -0.072954784731363262918  0                       -0.851236239662007619739
  0.539357840802981787532  |  0                       -0.777333643644968233538  0
  0.1                      |  0                        0                        0
  1                        |  0                        0.259111214548322744512 -0.358278966717952089048
 ------------------------------------------------------------------------------------------------------
                           |  0                        0.05                     0
 ...continued...
  0                        |
  0.1                      |
  0.539357840802981787532  |
  0.809036761204472681298  |
  0.309036761204472681298  |
  0.981074190219795268254  |
  0.833333333333333333333  |
  0.354017365856802376329  | -0.059578021181736100156
  0.882527661964732346425  |  0.321763705601778390100  0.510485725608063031577
  0.642615758240322548157  | -0.333269719096256706589  0.499269229556880061353  0.509504608929686104236
  0.357384241759677451842  |  0.239527360324390649107  0.397774662368094639047  0.010755895687360745555
  0.117472338035267653574  | -0.057695414616854888173  0.194781915712104164976  0.145384923188325069727
  0.833333333333333333333  | -0.861164195027635666673  5.784288136375372200229  3.288077619851035668904
  0.309036761204472681298  |  0.398320112318533301719  3.639372631810356060294  1.548228770398303223653
  0.539357840802981787532  | -0.091089566215517606959  0                        0
  0.1                      |  0                        0                        0
  1                        | -1.045948959408833060950  0.930327845415626983292  1.779509594317081024461
 ------------------------------------------------------------------------------------------------------
                           |  0.04                     0                        0.189237478148923490158
 ...continued...
  0                        |
  0.1                      |
  0.539357840802981787532  |
  0.809036761204472681298  |
  0.309036761204472681298  |
  0.981074190219795268254  |
  0.833333333333333333333  |
  0.354017365856802376329  |
  0.882527661964732346425  |
  0.642615758240322548157  |
  0.357384241759677451842  | -0.327769124164018874147
  0.117472338035267653574  | -0.078294271035167077755 -0.114503299361098912184
  0.833333333333333333333  | -2.386339050931363840134 -3.254793424836439186545 -2.163435416864229823539
  0.309036761204472681298  | -2.122217147040537160260 -1.583503985453261727133 -1.715616082859362649220
  0.539357840802981787532  |  0                        0                        0
  0.1                      |  0                        0                        0
  1                        |  0.1                     -0.282547569539044081612 -0.159327350119972549169
 ------------------------------------------------------------------------------------------------------
                           |  0.277429188517743176508  0.277429188517743176508  0.189237478148923490158
 ...continued...
  0                        |
  0.1                      |
  0.539357840802981787532  |
  0.809036761204472681298  |
  0.309036761204472681298  |
  0.981074190219795268254  |
  0.833333333333333333333  |
  0.354017365856802376329  |
  0.882527661964732346425  |
  0.642615758240322548157  |
  0.357384241759677451842  |
  0.117472338035267653574  |
  0.833333333333333333333  |
  0.309036761204472681298  | -0.024403640575012745213
  0.539357840802981787532  |  0.091089566215517606959  0.777333643644968233538
  0.1                      |  0                        0                        0.157178665799771163367
  1                        | -0.145515894647001510860 -0.259111214548322744512 -0.342758159847189839942 -0.675
 -------------------------------------------------------------------------------------------------------------
                           | -0.04                    -0.05                    -0.033333333333333333333 -0.025
 ...continued...
  0                        |
  0.1                      |
  0.539357840802981787532  |
  0.809036761204472681298  |
  0.309036761204472681298  |
  0.981074190219795268254  |
  0.833333333333333333333  |
  0.354017365856802376329  |
  0.882527661964732346425  |
  0.642615758240322548157  |
  0.357384241759677451842  |
  0.117472338035267653574  |
  0.833333333333333333333  |
  0.309036761204472681298  |
  0.539357840802981787532  |
  0.1                      |
  1                        |
 ---------------------------------------------------
                           | 0.033333333333333333333

Bibliography

[1] A family of embedded Runge-Kutta formulae, Dormand, J. R., Prince, P. J. (1980), Journal of Computational and Applied Mathematics 6 (1): 19–26, doi:10.1016/0771-050X(80)90013-3.

[2] A New Embedded Pair of Runge-Kutta Formulas of orders 5 and 6, M. Calvo, J.I. Montijano, L. Randez, Computers & Mathematics with Applications, Volume 20, Issue 1, 1990, Pages 15–24, ISSN 0898-1221, http://dx.doi.org/10.1016/0898-1221(90)90064-Q.

[3] A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides, J. R. Cash, A. H. Karp, ACM Transactions on Mathematical Software, vol. 16, pp. 201–222, 1990, doi:10.1145/79505.79507.

[4] A tenth-order Runge-Kutta method with error estimate, Feagin, T., Proceedings of the IAENG Conf. on Scientific Computing. 2007.

Used By

module~~foodie_integrator_runge_kutta_emd~~UsedByGraph module~foodie_integrator_runge_kutta_emd foodie_integrator_runge_kutta_emd module~foodie foodie module~foodie_integrator_runge_kutta_emd->module~foodie program~integrate_burgers integrate_burgers module~foodie->program~integrate_burgers program~integrate_lorenz integrate_lorenz module~foodie->program~integrate_lorenz module~type_euler_1d type_euler_1D module~foodie->module~type_euler_1d module~foodie_test_integrand_tester_object foodie_test_integrand_tester_object module~foodie->module~foodie_test_integrand_tester_object module~foodie_test_lcce_test foodie_test_lcce_test module~foodie->module~foodie_test_lcce_test module~type_euler_1d_openmp type_euler_1D_openmp module~foodie->module~type_euler_1d_openmp program~integrate_euler_1d_caf~2 integrate_euler_1D_caf module~foodie->program~integrate_euler_1d_caf~2 module~foodie_test_oscillation_test foodie_test_oscillation_test module~foodie->module~foodie_test_oscillation_test program~integrate_euler_1d integrate_euler_1D module~foodie->program~integrate_euler_1d module~foodie_test_object foodie_test_object module~foodie->module~foodie_test_object module~foodie_test_integrand_oscillation foodie_test_integrand_oscillation module~foodie->module~foodie_test_integrand_oscillation module~foodie_test_integrand_lcce foodie_test_integrand_lcce module~foodie->module~foodie_test_integrand_lcce module~type_euler_1d_caf type_euler_1D_caf module~foodie->module~type_euler_1d_caf module~type_burgers type_burgers module~foodie->module~type_burgers module~foodie_test_integrand_ladvection foodie_test_integrand_ladvection module~foodie->module~foodie_test_integrand_ladvection program~integrate_euler_1d_openmp integrate_euler_1D_openmp module~foodie->program~integrate_euler_1d_openmp module~type_lorenz type_lorenz module~foodie->module~type_lorenz module~type_euler_1d->program~integrate_euler_1d module~foodie_test_integrand_tester_object->module~foodie_test_object module~foodie_test_integrand_tester_object->module~foodie_test_integrand_oscillation module~foodie_test_integrand_tester_object->module~foodie_test_integrand_lcce module~foodie_test_integrand_tester_object->module~foodie_test_integrand_ladvection program~foodie_test_lcce foodie_test_lcce module~foodie_test_lcce_test->program~foodie_test_lcce module~type_euler_1d_openmp->program~integrate_euler_1d_openmp program~foodie_test_oscillation foodie_test_oscillation module~foodie_test_oscillation_test->program~foodie_test_oscillation program~foodie_tester foodie_tester module~foodie_test_object->program~foodie_tester module~foodie_test_integrand_oscillation->module~foodie_test_oscillation_test module~foodie_test_integrand_oscillation->module~foodie_test_object module~foodie_test_integrand_lcce->module~foodie_test_lcce_test module~foodie_test_integrand_lcce->module~foodie_test_object module~type_euler_1d_caf->program~integrate_euler_1d_caf~2 module~type_burgers->program~integrate_burgers module~foodie_test_integrand_ladvection->module~foodie_test_object module~type_lorenz->program~integrate_lorenz
Help

Variables

Derived Types

Functions

Subroutines


Variables

TypeVisibility AttributesNameInitial
character(len=99), private, parameter:: class_name_ ='runge_kutta_emd'

Name of the class of schemes.
logical, private, parameter:: has_fast_mode_ =.true.

Flag to check if integrator provides fast mode integrate.
character(len=99), private, parameter:: supported_schemes_(1:4) =[trim(class_name_)//'_stages_6_order_5  ', trim(class_name_)//'_stages_7_order_4  ', trim(class_name_)//'_stages_9_order_6  ', trim(class_name_)//'_stages_17_order_10']

List of supported schemes.

Derived Types

type, public, extends(integrator_multistage_object) :: integrator_runge_kutta_emd

FOODIE integrator: provide an explicit class of embedded Runge-Kutta schemes, from 2nd to 10th order accurate.

Components

TypeVisibility AttributesNameInitial
class(integrand_object), private, allocatable:: U1

First U evaluation.
class(integrand_object), private, allocatable:: U2

Second U evaluation.
real(kind=R_P), private, allocatable:: alph(:,:)

\(\alpha\) Butcher's coefficients.
real(kind=R_P), private, allocatable:: beta(:,:)

\(\beta\) Butcher's coefficients.
class(integrand_object), public, allocatable:: buffer

Buffer used for fast integration.
character(len=:), public, allocatable:: description_

Informative description of the integrator.
integer(kind=I_P), public :: error =0

Error status code.
character(len=:), public, allocatable:: error_message

Error message, hopefully meaningful.
real(kind=R_P), private, allocatable:: gamm(:)

\(\gamma\) Butcher's coefficients.
real(kind=R_P), private :: pp1_inv

\(1/(p+1)\) where p is the accuracy order of the lower accurate scheme.
integer(kind=I_P), public :: registers

Number of registers used for stages.
class(integrand_object), public, allocatable:: stage(:)

Stages.
integer(kind=I_P), public :: stages

Number of stages.
real(kind=R_P), private :: tolerance

Tolerance on the local truncation error.

Type-Bound Procedures

procedure, public, pass(self) :: allocate_integrand_members

Allocate integrand members.

procedure, public, pass(lhs) :: assign_abstract

Assign ony members of abstract integrator_object type.

procedure, public, pass(lhs) :: assign_multistage

Assign members of integrator_multistage_object and parents.

generic, public :: assignment(=) => integr_assign_integr

Overload =.

procedure, public, pass(self) :: check_error

Check for error occurrencies.

procedure, public, pass(self) :: class_name

Return the class name of schemes.

procedure, public, pass(self) :: description

Return informative integrator description.

procedure, public, pass(self) :: destroy

Destroy the integrator.

procedure, public, pass(self) :: destroy_abstract

Destroy only members of abstract integrator_object type.

procedure, public, pass(self) :: destroy_multistage

Destroy the integrator.

procedure, public, pass(self) :: has_fast_mode

Return .true. if the integrator class has fast mode integrate.

procedure, public, pass(self) :: initialize

Initialize (create) the integrator.

procedure, public, pass(lhs) :: integr_assign_integr

Operator =.

procedure, public, pass(self) :: integrate

Integrate integrand field.

procedure, public, pass(self) :: integrate_fast

Integrate integrand field, fast mode.

procedure, public, pass(self) :: is_multistage

Return .true. for multistage integrator.

procedure, public, pass(self) :: is_multistep

Return .true. for multistep integrator.

procedure, public, pass(self) :: is_supported

Return .true. if the integrator class support the given scheme.

procedure, private, pass(self) :: new_Dt

Compute new estimation of the time step Dt.

procedure, public, pass(self) :: stages_number

Return number of stages used.

procedure, public, pass(self) :: steps_number

Return number of steps used.

procedure, public, pass(self) :: supported_schemes

Return the list of supported schemes.

procedure, public, pass(self) :: trigger_error

Trigger an error.

Functions

private pure function class_name(self)

Return the class name of schemes.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(in) :: self

Integrator.

Return Value character(len=99)

Class name.

private elemental function has_fast_mode(self)

Return .true. if the integrator class has fast mode integrate.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(in) :: self

Integrator.

Return Value logical

Inquire result.

private elemental function is_supported(self, scheme)

Return .true. if the integrator class support the given scheme.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(in) :: self

Integrator.
character(len=*), intent(in) :: scheme

Selected scheme.

Return Value logical

Inquire result.

private pure function supported_schemes(self) result(schemes)

Return the list of supported schemes.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(in) :: self

Integrator.

Return Value character(len=99), allocatable, (:)

Queried scheme.


Subroutines

private subroutine allocate_integrand_members(self, U)

Allocate members of interpolator being of integrand_object class.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: self

Integrator.
class(integrand_object), intent(in) :: U

Integrand.

private elemental subroutine destroy(self)

Destroy the integrator.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: self

Integrator.

private subroutine initialize(self, scheme, U, tolerance, stop_on_fail)

Create the actual RK integrator: initialize the Butcher' table coefficients.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: self

Integrator.
character(len=*), intent(in) :: scheme

Selected scheme.
class(integrand_object), intent(in), optional :: U

Integrand molding prototype.
real(kind=R_P), intent(in), optional :: tolerance

Tolerance on the local truncation error (default 0.01).
logical, intent(in), optional :: stop_on_fail

Stop execution if initialization fail.

private subroutine integr_assign_integr(lhs, rhs)

Operator =.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: lhs

Left hand side.
class(integrator_object), intent(in) :: rhs

Right hand side.

private subroutine integrate(self, U, Dt, t, new_Dt)

Integrate field with explicit embedded Runge-Kutta scheme.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: self

Integrator.
class(integrand_object), intent(inout) :: U

Field to be integrated.
real(kind=R_P), intent(in) :: Dt

Time step.
real(kind=R_P), intent(in) :: t

Time.
real(kind=R_P), intent(out), optional :: new_Dt

New adapted time step.

private subroutine integrate_fast(self, U, Dt, t, new_Dt)

Integrate field with explicit embedded Runge-Kutta scheme, fast mode.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(inout) :: self

Integrator.
class(integrand_object), intent(inout) :: U

Field to be integrated.
real(kind=R_P), intent(in) :: Dt

Time step.
real(kind=R_P), intent(in) :: t

Time.
real(kind=R_P), intent(out), optional :: new_Dt

New adapted time step.

private elemental subroutine new_Dt(self, error, Dt)

Compute new estimation of the time step Dt.

Arguments

Type IntentOptional AttributesName
class(integrator_runge_kutta_emd), intent(in) :: self

Integrator.
real(kind=R_P), intent(in) :: error

Local truncation error estimation.
real(kind=R_P), intent(inout) :: Dt

Time step.

All Modules