FOODIE integrator: provide an explicit class of embedded Runge-Kutta schemes, from 2nd to 10th order accurate.
!< FOODIE integrator: provide an explicit class of embedded Runge-Kutta schemes, from 2nd to 10th order accurate. module foodie_integrator_runge_kutta_emd !< 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) $$ !< !< @note The value of \(\Delta t\) must be provided, it not being computed by the integrator. !< !< 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: !< !< @bug Presently, the 2 stages Heun-Euler seems to not work, **do not use it**. !< !<##### 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. use foodie_error_codes, only : ERROR_UNSUPPORTED_SCHEME use foodie_integrand_object, only : integrand_object use foodie_integrator_multistage_object, only : integrator_multistage_object use foodie_integrator_object, only : integrator_object use penf, only : I_P, R_P implicit none private public :: integrator_runge_kutta_emd character(len=99), parameter :: class_name_='runge_kutta_emd' !< Name of the class of schemes. character(len=99), 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. logical, parameter :: has_fast_mode_=.true. !< Flag to check if integrator provides *fast mode* integrate. type, 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. !< !< @note The integrator must be created or initialized (initialize the RK coefficients) before used. private real(R_P) :: tolerance !< Tolerance on the local truncation error. real(R_P) :: pp1_inv !< \(1/(p+1)\) where p is the accuracy order of the lower accurate scheme. real(R_P), allocatable :: alph(:,:) !< \(\alpha\) Butcher's coefficients. real(R_P), allocatable :: beta(:,:) !< \(\beta\) Butcher's coefficients. real(R_P), allocatable :: gamm(:) !< \(\gamma\) Butcher's coefficients. class(integrand_object), allocatable :: U1 !< First U evaluation. class(integrand_object), allocatable :: U2 !< Second U evaluation. contains ! deferred methods procedure, pass(self) :: class_name !< Return the class name of schemes. procedure, pass(self) :: has_fast_mode !< Return .true. if the integrator class has *fast mode* integrate. procedure, pass(lhs) :: integr_assign_integr !< Operator `=`. procedure, pass(self) :: integrate !< Integrate integrand field. procedure, pass(self) :: integrate_fast !< Integrate integrand field, fast mode. procedure, pass(self) :: is_supported !< Return .true. if the integrator class support the given scheme. procedure, pass(self) :: supported_schemes !< Return the list of supported schemes. ! public methods procedure, pass(self) :: destroy !< Destroy the integrator. procedure, pass(self) :: initialize !< Initialize (create) the integrator. ! overridden public methods procedure, pass(self) :: allocate_integrand_members !< Allocate integrand members. ! private methods procedure, pass(self), private :: new_Dt !< Compute new estimation of the time step Dt. endtype integrator_runge_kutta_emd contains ! deferred methods pure function class_name(self) !< Return the class name of schemes. class(integrator_runge_kutta_emd), intent(in) :: self !< Integrator. character(len=99) :: class_name !< Class name. class_name = trim(adjustl(class_name_)) endfunction class_name elemental function has_fast_mode(self) !< Return .true. if the integrator class has *fast mode* integrate. class(integrator_runge_kutta_emd), intent(in) :: self !< Integrator. logical :: has_fast_mode !< Inquire result. has_fast_mode = has_fast_mode_ endfunction has_fast_mode subroutine integr_assign_integr(lhs, rhs) !< Operator `=`. class(integrator_runge_kutta_emd), intent(inout) :: lhs !< Left hand side. class(integrator_object), intent(in) :: rhs !< Right hand side. call lhs%assign_abstract(rhs=rhs) select type(rhs) class is (integrator_runge_kutta_emd) lhs%tolerance = rhs%tolerance lhs%pp1_inv = rhs%pp1_inv lhs%stages = rhs%stages if (allocated(rhs%alph)) lhs%alph = rhs%alph if (allocated(rhs%beta)) lhs%beta = rhs%beta if (allocated(rhs%gamm)) lhs%gamm = rhs%gamm endselect endsubroutine integr_assign_integr subroutine integrate(self, U, Dt, t, new_Dt) !< Integrate field with explicit embedded Runge-Kutta scheme. !< !< The time steps is adaptively resized using the local truncation error estimation by means of the embedded pairs of RK schemes. class(integrator_runge_kutta_emd), intent(inout) :: self !< Integrator. class(integrand_object), intent(inout) :: U !< Field to be integrated. real(R_P), intent(in) :: Dt !< Time step. real(R_P), intent(in) :: t !< Time. real(R_P), intent(out), optional :: new_Dt !< New adapted time step. real(R_P) :: Dt_ !< Time step, local variable. real(R_P) :: error !< Local truncation error estimation. integer(I_P) :: s !< First stages counter. integer(I_P) :: ss !< Second stages counter. Dt_ = Dt do ! compute stages do s=1, self%stages self%stage(s) = U do ss=1, s - 1 self%stage(s) = self%stage(s) + (self%stage(ss) * (Dt_ * self%alph(s, ss))) enddo self%stage(s) = self%stage(s)%t(t=t + self%gamm(s) * Dt_) enddo ! compute new time step self%U1 = U self%U2 = U do s=1, self%stages self%U1 = self%U1 + (self%stage(s) * (Dt_ * self%beta(s, 1))) self%U2 = self%U2 + (self%stage(s) * (Dt_ * self%beta(s, 2))) enddo error = self%U2.lterror.self%U1 if (error <= self%tolerance) exit call self%new_Dt(error=error, Dt=Dt_) enddo U = self%U1 if (present(new_Dt)) new_Dt = Dt_ endsubroutine integrate subroutine integrate_fast(self, U, Dt, t, new_Dt) !< Integrate field with explicit embedded Runge-Kutta scheme, fast mode. !< !< The time steps is adaptively resized using the local truncation error estimation by means of the embedded pairs of RK schemes. class(integrator_runge_kutta_emd), intent(inout) :: self !< Integrator. class(integrand_object), intent(inout) :: U !< Field to be integrated. real(R_P), intent(in) :: Dt !< Time step. real(R_P), intent(in) :: t !< Time. real(R_P), intent(out), optional :: new_Dt !< New adapted time step. real(R_P) :: Dt_ !< Time step, local variable. real(R_P) :: error !< Local truncation error estimation. integer(I_P) :: s !< First stages counter. integer(I_P) :: ss !< Second stages counter. Dt_ = Dt do ! compute stages do s=1, self%stages self%stage(s) = U do ss=1, s - 1 call self%buffer%multiply_fast(lhs=self%stage(ss), rhs=Dt_ * self%alph(s, ss)) call self%stage(s)%add_fast(lhs=self%stage(s), rhs=self%buffer) enddo call self%stage(s)%t_fast(t=t + self%gamm(s) * Dt_) enddo ! compute new time step self%U1 = U self%U2 = U do s=1, self%stages call self%buffer%multiply_fast(lhs=self%stage(s), rhs=Dt_ * self%beta(s, 1)) call self%U1%add_fast(lhs=self%U1, rhs=self%buffer) call self%buffer%multiply_fast(lhs=self%stage(s), rhs=Dt_ * self%beta(s, 2)) call self%U2%add_fast(lhs=self%U2, rhs=self%buffer) enddo error = self%U2.lterror.self%U1 if (error <= self%tolerance) exit call self%new_Dt(error=error, Dt=Dt_) enddo U = self%U1 if (present(new_Dt)) new_Dt = Dt_ endsubroutine integrate_fast elemental function is_supported(self, scheme) !< Return .true. if the integrator class support the given scheme. class(integrator_runge_kutta_emd), intent(in) :: self !< Integrator. character(*), intent(in) :: scheme !< Selected scheme. logical :: is_supported !< Inquire result. integer(I_P) :: s !< Counter. is_supported = .false. do s=lbound(supported_schemes_, dim=1), ubound(supported_schemes_, dim=1) if (trim(adjustl(scheme)) == trim(adjustl(supported_schemes_(s)))) then is_supported = .true. return endif enddo endfunction is_supported pure function supported_schemes(self) result(schemes) !< Return the list of supported schemes. class(integrator_runge_kutta_emd), intent(in) :: self !< Integrator. character(len=99), allocatable :: schemes(:) !< Queried scheme. allocate(schemes(lbound(supported_schemes_, dim=1):ubound(supported_schemes_, dim=1))) schemes = supported_schemes_ endfunction supported_schemes ! public methods elemental subroutine destroy(self) !< Destroy the integrator. class(integrator_runge_kutta_emd), intent(inout) :: self !< Integrator. call self%destroy_multistage self%tolerance = 0._R_P self%pp1_inv = 0._R_P if (allocated(self%alph)) deallocate(self%alph) if (allocated(self%beta)) deallocate(self%beta) if (allocated(self%gamm)) deallocate(self%gamm) endsubroutine destroy subroutine initialize(self, scheme, U, tolerance, stop_on_fail) !< Create the actual RK integrator: initialize the Butcher' table coefficients. class(integrator_runge_kutta_emd), intent(inout) :: self !< Integrator. character(*), intent(in) :: scheme !< Selected scheme. class(integrand_object), intent(in), optional :: U !< Integrand molding prototype. real(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. if (self%is_supported(scheme=scheme)) then call self%destroy self%description_ = trim(adjustl(scheme)) self%tolerance = 0.01_R_P ; if (present(tolerance)) self%tolerance = tolerance select case(trim(adjustl(scheme))) case('runge_kutta_emd_stages_2_order_2') ! do not use, seems to not work! self%stages = 2 allocate(self%beta(1:self%stages, 1:2 )) ; self%beta = 0._R_P allocate(self%alph(1:self%stages, 1:self%stages)) ; self%alph = 0._R_P allocate(self%gamm( 1:self%stages)) ; self%gamm = 0._R_P self%pp1_inv = 1._R_P/(2._R_P + 1._R_P) self%beta(1, 1) = 0.5_R_P ; self%beta(1, 2) = 5._R_P self%beta(2, 1) = 1._R_P ; self%beta(2, 2) = 0._R_P self%alph(2, 1) = 1._R_P self%gamm(2) = 1._R_P case('runge_kutta_emd_stages_6_order_5') self%stages = 6 allocate(self%beta(1:self%stages, 1:2 )) ; self%beta = 0._R_P allocate(self%alph(1:self%stages, 1:self%stages)) ; self%alph = 0._R_P allocate(self%gamm( 1:self%stages)) ; self%gamm = 0._R_P self%pp1_inv = 1._R_P/(5._R_P + 1._R_P) self%beta(1, 1) = 37._R_P/378._R_P ; self%beta(1, 2) = 2825._R_P/27648._R_P self%beta(2, 1) = 0._R_P ; self%beta(2, 2) = 0._R_P self%beta(3, 1) = 250._R_P/621._R_P ; self%beta(3, 2) = 18575._R_P/48384._R_P self%beta(4, 1) = 125._R_P/594._R_P ; self%beta(4, 2) = 13525._R_P/55296._R_P self%beta(5, 1) = 0._R_P ; self%beta(5, 2) = 277._R_P/14336._R_P self%beta(6, 1) = 512._R_P/1771._R_P ; self%beta(6, 2) = 1._R_P/4._R_P self%alph(2, 1)=1._R_P/5._R_P self%alph(3, 1)=3._R_P/40._R_P ;self%alph(3, 2)= 9._R_P/40._R_P self%alph(4, 1)=3._R_P/10._R_P ;self%alph(4, 2)= -9._R_P/10._R_P ;self%alph(4, 3)=6._R_P/5._R_P self%alph(5, 1)=-11._R_P/54._R_P ;self%alph(5, 2)= 5._R_P/2._R_P ;self%alph(5, 3)=-70._R_P/27._R_P self%alph(6, 1)=1631._R_P/55296._R_P;self%alph(6, 2)=175._R_P/512._R_P;self%alph(6, 3)=575._R_P/13824._R_P self%alph(5, 4)=35._R_P/27._R_P self%alph(6, 4)=44275._R_P/110592._R_P;self%alph(6, 5)=253._R_P/4096._R_P self%gamm(2) = 1._R_P/5._R_P self%gamm(3) = 3._R_P/10._R_P self%gamm(4) = 3._R_P/5._R_P self%gamm(5) = 1._R_P self%gamm(6) = 7._R_P/8._R_P case('runge_kutta_emd_stages_7_order_4') self%stages = 7 allocate(self%beta(1:self%stages, 1:2 )) ; self%beta = 0._R_P allocate(self%alph(1:self%stages, 1:self%stages)) ; self%alph = 0._R_P allocate(self%gamm( 1:self%stages)) ; self%gamm = 0._R_P self%pp1_inv = 1._R_P/(4._R_P + 1._R_P) self%beta(1, 1) = 35._R_P/384._R_P ; self%beta(1, 2) = 5179._R_P/57600._R_P self%beta(2, 1) = 0._R_P ; self%beta(2, 2) = 0._R_P self%beta(3, 1) = 500._R_P/1113._R_P ; self%beta(3, 2) = 7571._R_P/16695._R_P self%beta(4, 1) = 125._R_P/192._R_P ; self%beta(4, 2) = 393._R_P/640._R_P self%beta(5, 1) = -2187._R_P/6784._R_P ; self%beta(5, 2) = -92097._R_P/339200._R_P self%beta(6, 1) = 11._R_P/84._R_P ; self%beta(6, 2) = 187._R_P/2100._R_P self%beta(7, 1) = 0._R_P ; self%beta(7, 2) = 1._R_P/40._R_P self%alph(2, 1)=1._R_P/5._R_P self%alph(3, 1)=3._R_P/40._R_P ;self%alph(3, 2)= 9._R_P/40._R_P self%alph(4, 1)=44._R_P/45._R_P ;self%alph(4, 2)=-56._R_P/15._R_P ;self%alph(4, 3)=32._R_P/9._R_P self%alph(5, 1)=19372._R_P/6561._R_P;self%alph(5, 2)=-25360._R_P/2187._R_P;self%alph(5, 3)=64448._R_P/6561._R_P self%alph(6, 1)=9017._R_P/3168._R_P ;self%alph(6, 2)=-355._R_P/33._R_P ;self%alph(6, 3)=46732._R_P/5247._R_P self%alph(7, 1)=35._R_P/384._R_P ;self%alph(7, 2)= 0._R_P ;self%alph(7, 3)=500._R_P/1113._R_P self%alph(5, 4)=-212._R_P/729._R_P self%alph(6, 4)= 49._R_P/176._R_P ;self%alph(6, 5)=-5103._R_P/18656._R_P self%alph(7, 4)= 125._R_P/192._R_P;self%alph(7, 5)=-2187._R_P/6784._R_P ;self%alph(7, 6)=11._R_P/84._R_P self%gamm(2) = 1._R_P/5._R_P self%gamm(3) = 3._R_P/10._R_P self%gamm(4) = 4._R_P/5._R_P self%gamm(5) = 8._R_P/9._R_P self%gamm(6) = 1._R_P self%gamm(7) = 1._R_P case('runge_kutta_emd_stages_9_order_6') self%stages = 9 allocate(self%beta(1:self%stages, 1:2 )) ; self%beta = 0._R_P allocate(self%alph(1:self%stages, 1:self%stages)) ; self%alph = 0._R_P allocate(self%gamm( 1:self%stages)) ; self%gamm = 0._R_P self%pp1_inv = 1._R_P/(6._R_P + 1._R_P) self%beta(1, 1) = 17572349._R_P/289262523._R_P ; self%beta(1, 2) = 15231665._R_P/510830334._R_P self%beta(2, 1) = 0._R_P ; self%beta(2, 2) = 0._R_P self%beta(3, 1) = 57513011._R_P/201864250._R_P ; self%beta(3, 2) = 59452991._R_P/116050448._R_P self%beta(4, 1) = 15587306._R_P/354501571._R_P ; self%beta(4, 2) = -28398517._R_P/122437738._R_P self%beta(5, 1) = 71783021._R_P/234982865._R_P ; self%beta(5, 2) = 56673824._R_P/137010559._R_P self%beta(6, 1) = 29672000._R_P/180480167._R_P ; self%beta(6, 2) = 68003849._R_P/426673583._R_P self%beta(7, 1) = 65567621._R_P/127060952._R_P ; self%beta(7, 2) = 7097631._R_P/37564021._R_P self%beta(8, 1) = -79074570._R_P/210557597._R_P ; self%beta(8, 2) = -71226429._R_P/583093742._R_P self%beta(9, 1) = 0._R_P ; self%beta(9, 2) = 1._R_P/20._R_P self%alph(2, 1)=2._R_P/15._R_P self%alph(3, 1)=1._R_P/20._R_P ; self%alph(3, 2)=3._R_P/20._R_P self%alph(4, 1)=3._R_P/40._R_P ; self%alph(4, 2)=0._R_P self%alph(5, 1)=86727015._R_P/196851553._R_P ; self%alph(5, 2)=-60129073._R_P/52624712._R_P self%alph(6, 1)=-86860849._R_P/45628967._R_P ; self%alph(6, 2)=111022885._R_P/25716487._R_P self%alph(7, 1)=77759591._R_P/16096467._R_P ; self%alph(7, 2)=-49252809._R_P/6452555._R_P self%alph(8, 1)=237564263._R_P/39280295._R_P ; self%alph(8, 2)=-100523239._R_P/10677940._R_P self%alph(9, 1)=17572349._R_P/289262523._R_P ; self%alph(9, 2)=0._R_P self%alph(4, 3)=9._R_P/40._R_P self%alph(5, 3)=957436434._R_P/1378352377._R_P ; self%alph(5, 4)=83886832._R_P/147842441._R_P self%alph(6, 3)=108046682._R_P/101167669._R_P ; self%alph(6, 4)=-141756746._R_P/36005461._R_P self%alph(7, 3)=-381680111._R_P/51572984._R_P ; self%alph(7, 4)=879269579._R_P/66788831._R_P self%alph(8, 3)=-265574846._R_P/27330247._R_P ; self%alph(8, 4)=317978411._R_P/18988713._R_P self%alph(9, 3)=57513011._R_P/201864250._R_P ; self%alph(9, 4)=15587306._R_P/354501571._R_P self%alph(6, 5)=73139862._R_P/60170633._R_P self%alph(7, 5)=-90453121._R_P/33722162._R_P ; self%alph(7, 6)=111179552._R_P/157155827._R_P self%alph(8, 5)=-124494385._R_P/35453627._R_P ; self%alph(8, 6)=86822444._R_P/100138635._R_P self%alph(9, 5)=71783021._R_P/234982865._R_P ; self%alph(9, 6)=29672000._R_P/180480167._R_P self%alph(8, 7)=-12873523._R_P/724232625._R_P self%alph(9, 7)=65567621._R_P/127060952._R_P ;self%alph(9, 8)=-79074570._R_P/210557597._R_P self%gamm(2) = 2._R_P/15._R_P self%gamm(3) = 1._R_P/5._R_P self%gamm(4) = 3._R_P/10._R_P self%gamm(5) = 14._R_P/25._R_P self%gamm(6) = 19._R_P/25._R_P self%gamm(7) = 35226607._R_P/35688279._R_P self%gamm(8) = 1._R_P self%gamm(9) = 1._R_P case('runge_kutta_emd_stages_17_order_10') self%stages = 17 allocate(self%beta(1:self%stages, 1:2 )) ; self%beta = 0._R_P allocate(self%alph(1:self%stages, 1:self%stages)) ; self%alph = 0._R_P allocate(self%gamm( 1:self%stages)) ; self%gamm = 0._R_P self%pp1_inv = 1._R_P/(10._R_P + 1._R_P) self%beta(1 , 1) = 0.033333333333333333333_R_P; self%beta(1 , 2) = 0.033333333333333333333_R_P self%beta(2 , 1) = 0.025_R_P ; self%beta(2 , 2) = 1._R_P/36._R_P self%beta(3 , 1) = 0.033333333333333333333_R_P; self%beta(3 , 2) = 0.033333333333333333333_R_P self%beta(4 , 1) = 0._R_P ; self%beta(4 , 2) = 0._R_P self%beta(5 , 1) = 0.05_R_P ; self%beta(5 , 2) = 0.05_R_P self%beta(6 , 1) = 0._R_P ; self%beta(6 , 2) = 0._R_P self%beta(7 , 1) = 0.04_R_P ; self%beta(7 , 2) = 0.04_R_P self%beta(8 , 1) = 0._R_P ; self%beta(8 , 2) = 0._R_P self%beta(9 , 1) = 0.189237478148923490158_R_P; self%beta(9 , 2) = 0.189237478148923490158_R_P self%beta(10, 1) = 0.277429188517743176508_R_P; self%beta(10, 2) = 0.277429188517743176508_R_P self%beta(11, 1) = 0.277429188517743176508_R_P; self%beta(11, 2) = 0.277429188517743176508_R_P self%beta(12, 1) = 0.189237478148923490158_R_P; self%beta(12, 2) = 0.189237478148923490158_R_P self%beta(13, 1) =-0.04_R_P ; self%beta(13, 2) =-0.04_R_P self%beta(14, 1) =-0.05_R_P ; self%beta(14, 2) =-0.05_R_P self%beta(15, 1) =-0.033333333333333333333_R_P; self%beta(15, 2) =-0.033333333333333333333_R_P self%beta(16, 1) =-0.025_R_P ; self%beta(16, 2) =-1._R_P/36._R_P self%beta(17, 1) = 0.033333333333333333333_R_P; self%beta(17, 2) = 0.033333333333333333333_R_P self%alph(2 , 1)= 0.1_R_P self%alph(3 , 1)=-0.915176561375291440520_R_P; self%alph(3 , 2)=1.454534402178273228052_R_P self%alph(4 , 1)= 0.202259190301118170324_R_P; self%alph(4 , 2)=0._R_P self%alph(5 , 1)= 0.184024714708643575149_R_P; self%alph(5 , 2)=0._R_P self%alph(6 , 1)= 0.087900734020668133731_R_P; self%alph(6 , 2)=0._R_P self%alph(7 , 1)= 0.085970050490246030218_R_P; self%alph(7 , 2)=0._R_P self%alph(8 , 1)= 0.120930449125333720660_R_P; self%alph(8 , 2)=0._R_P self%alph(9 , 1)= 0.110854379580391483508_R_P; self%alph(9 , 2)=0._R_P self%alph(10, 1)= 0.112054414752879004829_R_P; self%alph(10, 2)=0._R_P self%alph(11, 1)= 0.113976783964185986138_R_P; self%alph(11, 2)=0._R_P self%alph(12, 1)= 0.079831452828019604635_R_P; self%alph(12, 2)=0._R_P self%alph(13, 1)= 0.985115610164857280120_R_P; self%alph(13, 2)=0._R_P self%alph(14, 1)= 0.895080295771632891049_R_P; self%alph(14, 2)=0._R_P self%alph(15, 1)=-0.915176561375291440520_R_P; self%alph(15, 2)=1.454534402178273228052_R_P self%alph(16, 1)= 0.1_R_P ; self%alph(16, 2)=0._R_P self%alph(17, 1)= 0.181781300700095283888_R_P; self%alph(17, 2)=0.675_R_P self%alph(4 , 3)= 0.606777570903354510974_R_P self%alph(5 , 3)= 0.197966831227192369068_R_P; self%alph(5 , 4)=-0.072954784731363262918_R_P self%alph(6 , 3)= 0._R_P ; self%alph(6 , 4)= 0.410459702520260645318_R_P self%alph(7 , 3)= 0._R_P ; self%alph(7 , 4)= 0.330885963040722183948_R_P self%alph(8 , 3)= 0._R_P ; self%alph(8 , 4)= 0._R_P self%alph(9 , 3)= 0._R_P ; self%alph(9 , 4)= 0._R_P self%alph(10, 3)= 0._R_P ; self%alph(10, 4)= 0._R_P self%alph(11, 3)= 0._R_P ; self%alph(11, 4)= 0._R_P self%alph(12, 3)= 0._R_P ; self%alph(12, 4)= 0._R_P self%alph(13, 3)= 0._R_P ; self%alph(13, 4)= 0.330885963040722183948_R_P self%alph(14, 3)= 0.197966831227192369068_R_P; self%alph(14, 4)=-0.072954784731363262918_R_P self%alph(15, 3)= 0._R_P ; self%alph(15, 4)= 0._R_P self%alph(16, 3)=-0.157178665799771163367_R_P; self%alph(16, 4)= 0._R_P self%alph(17, 3)= 0.342758159847189839942_R_P; self%alph(17, 4)= 0._R_P self%alph(6 , 5)= 0.482713753678866489204_R_P self%alph(7 , 5)= 0.489662957309450192844_R_P; self%alph(7 , 6)=-0.073185637507085073678_R_P self%alph(8 , 5)= 0.260124675758295622809_R_P; self%alph(8 , 6)= 0.032540262154909133015_R_P self%alph(9 , 5)= 0._R_P ; self%alph(9 , 6)=-0.060576148825500558762_R_P self%alph(10, 5)= 0._R_P ; self%alph(10, 6)=-0.144942775902865915672_R_P self%alph(11, 5)= 0._R_P ; self%alph(11, 6)=-0.076881336420335693858_R_P self%alph(12, 5)= 0._R_P ; self%alph(12, 6)=-0.052032968680060307651_R_P self%alph(13, 5)= 0.489662957309450192844_R_P; self%alph(13, 6)=-1.378964865748435675821_R_P self%alph(14, 5)= 0._R_P ; self%alph(14, 6)=-0.851236239662007619739_R_P self%alph(15, 5)=-0.777333643644968233538_R_P; self%alph(15, 6)= 0._R_P self%alph(16, 5)= 0._R_P ; self%alph(16, 6)= 0._R_P self%alph(17, 5)= 0.259111214548322744512_R_P; self%alph(17, 6)=-0.358278966717952089048_R_P self%alph(8 , 7)=-0.059578021181736100156_R_P self%alph(9 , 7)= 0.321763705601778390100_R_P; self%alph(9 , 8)=0.510485725608063031577_R_P self%alph(10, 7)=-0.333269719096256706589_R_P; self%alph(10, 8)=0.499269229556880061353_R_P self%alph(11, 7)= 0.239527360324390649107_R_P; self%alph(11, 8)=0.397774662368094639047_R_P self%alph(12, 7)=-0.057695414616854888173_R_P; self%alph(12, 8)=0.194781915712104164976_R_P self%alph(13, 7)=-0.861164195027635666673_R_P; self%alph(13, 8)=5.784288136375372200229_R_P self%alph(14, 7)= 0.398320112318533301719_R_P; self%alph(14, 8)=3.639372631810356060294_R_P self%alph(15, 7)=-0.091089566215517606959_R_P; self%alph(15, 8)=0._R_P self%alph(16, 7)= 0._R_P ; self%alph(16, 8)=0._R_P self%alph(17, 7)=-1.045948959408833060950_R_P; self%alph(17, 8)=0.930327845415626983292_R_P self%alph(10, 9)=0.509504608929686104236_R_P self%alph(11, 9)=0.010755895687360745555_R_P; self%alph(11, 10)=-0.327769124164018874147_R_P self%alph(12, 9)=0.145384923188325069727_R_P; self%alph(12, 10)=-0.078294271035167077755_R_P self%alph(13, 9)=3.288077619851035668904_R_P; self%alph(13, 10)=-2.386339050931363840134_R_P self%alph(14, 9)=1.548228770398303223653_R_P; self%alph(14, 10)=-2.122217147040537160260_R_P self%alph(15, 9)=0._R_P ; self%alph(15, 10)= 0._R_P self%alph(16, 9)=0._R_P ; self%alph(16, 10)= 0._R_P self%alph(17, 9)=1.779509594317081024461_R_P; self%alph(17, 10)= 0.1_R_P self%alph(12, 11)=-0.114503299361098912184_R_P self%alph(13, 11)=-3.254793424836439186545_R_P; self%alph(13, 12)=-2.163435416864229823539_R_P self%alph(14, 11)=-1.583503985453261727133_R_P; self%alph(14, 12)=-1.715616082859362649220_R_P self%alph(15, 11)= 0._R_P ; self%alph(15, 12)= 0._R_P self%alph(16, 11)= 0._R_P ; self%alph(16, 12)= 0._R_P self%alph(17, 11)=-0.282547569539044081612_R_P; self%alph(17, 12)=-0.159327350119972549169_R_P self%alph(14, 13)=-0.024403640575012745213_R_P self%alph(15, 13)= 0.091089566215517606959_R_P; self%alph(15, 14)= 0.777333643644968233538_R_P self%alph(16, 13)= 0._R_P ; self%alph(16, 14)= 0._R_P self%alph(17, 13)=-0.145515894647001510860_R_P; self%alph(17, 14)=-0.259111214548322744512_R_P self%alph(16, 15)= 0.157178665799771163367_R_P self%alph(17, 15)=-0.342758159847189839942_R_P; self%alph(17, 16)=-0.675_R_P self%gamm(2 ) = 0.1_R_P self%gamm(3 ) = 0.539357840802981787532_R_P self%gamm(4 ) = 0.809036761204472681298_R_P self%gamm(5 ) = 0.309036761204472681298_R_P self%gamm(6 ) = 0.981074190219795268254_R_P self%gamm(7 ) = 0.833333333333333333333_R_P self%gamm(8 ) = 0.354017365856802376329_R_P self%gamm(9 ) = 0.882527661964732346425_R_P self%gamm(10) = 0.642615758240322548157_R_P self%gamm(11) = 0.357384241759677451842_R_P self%gamm(12) = 0.117472338035267653574_R_P self%gamm(13) = 0.833333333333333333333_R_P self%gamm(14) = 0.309036761204472681298_R_P self%gamm(15) = 0.539357840802981787532_R_P self%gamm(16) = 0.1_R_P self%gamm(17) = 1._R_P endselect self%registers = self%stages if (present(U)) call self%allocate_integrand_members(U=U) else call self%trigger_error(error=ERROR_UNSUPPORTED_SCHEME, & error_message='"'//trim(adjustl(scheme))//'" unsupported scheme', & is_severe=stop_on_fail) endif endsubroutine initialize ! overridden public methods subroutine allocate_integrand_members(self, U) !< Allocate members of interpolator being of [[integrand_object]] class. !< !< @note It is assumed that the integrator has been properly initialized before calling this method. class(integrator_runge_kutta_emd), intent(inout) :: self !< Integrator. class(integrand_object), intent(in) :: U !< Integrand. integer(I_P) :: s !< Counter. if (self%is_multistage() .and. self%registers > 0) then if (allocated(self%stage)) deallocate(self%stage) allocate(self%stage(1:self%registers), mold=U) do s=1, self%registers self%stage(s) = U enddo endif if (self%has_fast_mode()) then if (allocated(self%buffer)) deallocate(self%buffer) allocate(self%buffer, mold=U) self%buffer = U endif if (allocated(self%U1)) deallocate(self%U1) allocate(self%U1, mold=U) self%U1 = U if (allocated(self%U2)) deallocate(self%U2) allocate(self%U2, mold=U) self%U2 = U endsubroutine allocate_integrand_members ! private methods elemental subroutine new_Dt(self, error, Dt) !< Compute new estimation of the time step Dt. !< !< The formula employed is: !< !< $$ Dt_{new} = 0.9 Dt_{old} \left( \frac{tolerance}{error} \right)^{\frac{1}{p+1}} $$ !< !< @note 0.9 is a safety factor. class(integrator_runge_kutta_emd), intent(in) :: self !< Integrator. real(R_P), intent(in) :: error !< Local truncation error estimation. real(R_P), intent(inout) :: Dt !< Time step. if (error>self%tolerance) Dt = 0.9_R_P * Dt * (self%tolerance/error) ** self%pp1_inv endsubroutine new_Dt endmodule foodie_integrator_runge_kutta_emd