diff --git a/papers/announcement/errors_analysis/errors_analysis.sh b/papers/announcement/errors_analysis/errors_analysis.sh index 923f9721..3e2c1f05 100755 --- a/papers/announcement/errors_analysis/errors_analysis.sh +++ b/papers/announcement/errors_analysis/errors_analysis.sh @@ -11,6 +11,7 @@ euler=0 lf=0 lf_raw=0 lmm_ssp=0 +lmm_ssp_vss=0 ls_rk=0 tvd_rk=0 @@ -20,18 +21,19 @@ function print_usage { echo "script automotion for FOODIE's errors analysis" echo "Usage: `basename $0` #solver" echo "Valid solver switches are:" - echo " -ab # => adams-bashforth" - echo " -abm # => adams-bashforth-moulton" - echo " -am # => adams-moulton" - echo " -bdf # => backward differentiation formula" - echo " -emd-rk # => embeded RK" - echo " -euler # => euler" - echo " -lf # => leapfrog" - echo " -lf-raw # => leapfrog raw filtered" - echo " -lmm-ssp # => linear multistep methods SSP" - echo " -ls-rk # => low storage RK" - echo " -tvd-rk # => TVD RK" - echo " -all # => all solvers!" + echo " -ab # => adams-bashforth" + echo " -abm # => adams-bashforth-moulton" + echo " -am # => adams-moulton" + echo " -bdf # => backward differentiation formula" + echo " -emd-rk # => embeded RK" + echo " -euler # => euler" + echo " -lf # => leapfrog" + echo " -lf-raw # => leapfrog raw filtered" + echo " -lmm-ssp # => linear multistep methods SSP" + echo " -lmm-ssp-vss # => linear multistep methods SSP with Variable Step Size" + echo " -ls-rk # => low storage RK" + echo " -tvd-rk # => TVD RK" + echo " -all # => all solvers!" } #parsing command line @@ -64,6 +66,9 @@ while [ $# -gt 0 ]; do "-lmm-ssp") lmm_ssp=1 ;; + "-lmm-ssp-vss") + lmm_ssp_vss=1 + ;; "-ls-rk") ls_rk=1 ;; @@ -149,6 +154,10 @@ fi if [ $lmm_ssp -eq 1 ] ; then ../tests/oscillation -r -tf 1e6 --solver lmm-ssp --ss 3 5 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 50.d0 10.d0 --errors_analysis > analysis-lmm-ssp.log fi +if [ $lmm_ssp_vss -eq 1 ] ; then + ../tests/oscillation -r -tf 1e6 --solver lmm-ssp-vss --ss 3 5 --order 2 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis > analysis-lmm-ssp-vss.log + ../tests/oscillation -r -tf 1e6 --solver lmm-ssp-vss --ss 3 5 --order 3 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis > analysis-lmm-ssp-vss.log +fi if [ $ls_rk -eq 1 ] ; then ../tests/oscillation -r -tf 1e6 --solver ls-runge-kutta --ss 1 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis #> analysis-leapfrog.log ../tests/oscillation -r -tf 1e6 --solver ls-runge-kutta --ss 5 7 -Dt 5000.d0 2500.d0 1250.d0 625.d0 320.d0 100.d0 --errors_analysis #> analysis-leapfrog.log diff --git a/src/lib/foodie.f90 b/src/lib/foodie.f90 index e000e961..6fd0110c 100644 --- a/src/lib/foodie.f90 +++ b/src/lib/foodie.f90 @@ -72,6 +72,7 @@ module foodie use foodie_integrator_euler_explicit, only : integrator_euler_explicit use foodie_integrator_leapfrog, only : integrator_leapfrog use foodie_integrator_lmm_ssp, only : integrator_lmm_ssp +use foodie_integrator_lmm_ssp_vss, only : integrator_lmm_ssp_vss use foodie_integrator_runge_kutta_emd, only : integrator_runge_kutta_emd use foodie_integrator_runge_kutta_low_storage, only : integrator_runge_kutta_ls use foodie_integrator_runge_kutta_tvd, only : integrator_runge_kutta_tvd @@ -87,13 +88,14 @@ module foodie public :: integrator_euler_explicit public :: integrator_leapfrog public :: integrator_lmm_ssp +public :: integrator_lmm_ssp_vss public :: integrator_runge_kutta_emd public :: integrator_runge_kutta_ls public :: foodie_integrator public :: integrator_runge_kutta_tvd contains - function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(integrator) + function foodie_integrator(scheme, steps, stages, order, tolerance, nu, alpha) result(integrator) !< Return a concrete instance of [[integrator]] given a scheme selection. !< !< This is the FOODIE integrators factory. @@ -102,6 +104,7 @@ function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(i character(*), intent(in) :: scheme !< Selected integrator given. integer(I_P), intent(in), optional :: steps !< Number of time steps used in multi-step schemes. integer(I_P), intent(in), optional :: stages !< Number of Runge-Kutta stages used in multi-stage schemes. + integer(I_P), intent(in), optional :: order !< Order of accuracy. real(R_P), intent(in), optional :: tolerance !< Tolerance on the local truncation error. real(R_P), intent(in), optional :: nu !< Williams-Robert-Asselin filter coefficient. real(R_P), intent(in), optional :: alpha !< Robert-Asselin filter coefficient. @@ -176,6 +179,18 @@ function foodie_integrator(scheme, steps, stages, tolerance, nu, alpha) result(i error_message='missing steps number for initializing integrator!', & is_severe=.true.) endif + case('lmm_ssp_vss') + allocate(integrator_lmm_ssp_vss :: integrator) + if ((.not.present(steps)).or.(.not.present(order))) then + call integrator%trigger_error(error=ERROR_MISSING_STEPS_NUMBER, & + error_message='missing steps number for initializing integrator!', & + is_severe=.true.) + else + select type(integrator) + type is(integrator_lmm_ssp_vss) + call integrator%init(steps=steps, order=order) + endselect + endif case('runge_kutta_emd') allocate(integrator_runge_kutta_emd :: integrator) if (present(stages)) then diff --git a/src/lib/foodie_error_codes.f90 b/src/lib/foodie_error_codes.f90 index fe484743..a18efb57 100644 --- a/src/lib/foodie_error_codes.f90 +++ b/src/lib/foodie_error_codes.f90 @@ -11,9 +11,11 @@ module foodie_error_codes public :: ERROR_BAD_STEPS_NUMBER public :: ERROR_MISSING_STAGES_NUMBER public :: ERROR_BAD_STAGES_NUMBER +public :: ERROR_INTEGRATOR_INIT_FAIL integer(I_P), parameter :: ERROR_MISSING_STEPS_NUMBER = 1 !< Error missing steps number (for FOODIE factory). integer(I_P), parameter :: ERROR_BAD_STEPS_NUMBER = 2 !< Error bad (unsupported) steps number. integer(I_P), parameter :: ERROR_MISSING_STAGES_NUMBER = 3 !< Error missing stages number (for FOODIE factory). integer(I_P), parameter :: ERROR_BAD_STAGES_NUMBER = 4 !< Error bad (unsupported) stages number. +integer(I_P), parameter :: ERROR_INTEGRATOR_INIT_FAIL = 5 !< Error integrator initialization failed. endmodule foodie_error_codes diff --git a/src/lib/foodie_integrator_lmm_ssp.f90 b/src/lib/foodie_integrator_lmm_ssp.f90 index 2f45c456..efc18b5c 100644 --- a/src/lib/foodie_integrator_lmm_ssp.f90 +++ b/src/lib/foodie_integrator_lmm_ssp.f90 @@ -115,39 +115,39 @@ subroutine init(self, steps) if (self%is_supported(steps)) then call self%destroy self%steps = steps - allocate(self%a(1:steps)) ; self%a = 0.0_R_P - allocate(self%b(1:steps)) ; self%b = 0.0_R_P + allocate(self%a(1:steps)) ; self%a = 0._R_P + allocate(self%b(1:steps)) ; self%b = 0._R_P select case(steps) case(3) ! LMM-SSP(3,2) - self%a(1) = 1.0_R_P/4.0_R_P - self%a(2) = 0.0_R_P - self%a(3) = 3.0_R_P/4.0_R_P + self%a(1) = 1._R_P/4._R_P + self%a(2) = 0._R_P + self%a(3) = 3._R_P/4._R_P - self%b(1) = 0.0_R_P - self%b(2) = 0.0_R_P - self%b(3) = 3.0_R_P/2.0_R_P + self%b(1) = 0._R_P + self%b(2) = 0._R_P + self%b(3) = 3._R_P/2._R_P case(4) ! LMM-SSP(4,3) - self%a(1) = 11.0_R_P/27.0_R_P - self%a(2) = 0.0_R_P - self%a(3) = 0.0_R_P - self%a(4) = 16.0_R_P/27.0_R_P - - self%b(1) = 12.0_R_P/27.0_R_P - self%b(2) = 0.0_R_P - self%b(3) = 0.0_R_P - self%b(4) = 16.0_R_P/9.0_R_P + self%a(1) = 11._R_P/27._R_P + self%a(2) = 0._R_P + self%a(3) = 0._R_P + self%a(4) = 16._R_P/27._R_P + + self%b(1) = 12._R_P/27._R_P + self%b(2) = 0._R_P + self%b(3) = 0._R_P + self%b(4) = 16._R_P/9._R_P case(5) ! LMM-SSP(5,3) - self%a(1) = 7.0_R_P/32.0_R_P - self%a(2) = 0.0_R_P - self%a(3) = 0.0_R_P - self%a(4) = 0.0_R_P - self%a(5) = 25.0_R_P/32.0_R_P - - self%b(1) = 5.0_R_P/16.0_R_P - self%b(2) = 0.0_R_P - self%b(3) = 0.0_R_P - self%b(4) = 0.0_R_P - self%b(5) = 25.0_R_P/16.0_R_P + self%a(1) = 7._R_P/32._R_P + self%a(2) = 0._R_P + self%a(3) = 0._R_P + self%a(4) = 0._R_P + self%a(5) = 25._R_P/32._R_P + + self%b(1) = 5._R_P/16._R_P + self%b(2) = 0._R_P + self%b(3) = 0._R_P + self%b(4) = 0._R_P + self%b(5) = 25._R_P/16._R_P endselect else call self%trigger_error(error=ERROR_BAD_STEPS_NUMBER, & @@ -170,7 +170,9 @@ subroutine integrate(self, U, previous, Dt, t, autoupdate) autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate U = U * 0._R_P do s=1, self%steps - U = U + previous(s) * self%a(s) + previous(s)%t(t=t(s)) * (Dt * self%b(s)) + if (self%a(s) /= 0._R_P) U = U + previous(s) * self%a(s) + if (self%b(s) /= 0._R_P) U = U + previous(s)%t(t=t(s)) * (Dt * self%b(s)) + ! U = U + previous(s) * self%a(s) + previous(s)%t(t=t(s)) * (Dt * self%b(s)) enddo if (autoupdate_) call self%update_previous(U=U, previous=previous) endsubroutine integrate diff --git a/src/lib/foodie_integrator_lmm_ssp_vss.f90 b/src/lib/foodie_integrator_lmm_ssp_vss.f90 new file mode 100644 index 00000000..533cd8b6 --- /dev/null +++ b/src/lib/foodie_integrator_lmm_ssp_vss.f90 @@ -0,0 +1,278 @@ +!< FOODIE integrator: provide an explicit class of Linear Multi-step Methods (LLM) with Strong Stability Preserving property and +!< variable stepsize (VSS), from 2nd to 3rd order accurate. + +module foodie_integrator_lmm_ssp_vss +!< FOODIE integrator: provide an explicit class of Linear Multi-step Methods (LLM) with Strong Stability Preserving property and +!< variable stepsize (VSS), from 2nd to 3rd order accurate. +!< +!< 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 LMM-SSP class scheme implemented is: +!< +!< $$ U^{n+N_s} = \frac{1}{\Omega_{N_s-1}^2} U^n + \frac{\Omega_{N_s-1}^2 - 1}{\Omega_{N_s-1}^2} U^{n+N_s-1} + +!< \frac{\Omega_{N_s-1} + 1}{\Omega_{N_s-1}} \Delta t^{n+N_s} R(U^{n+N_s-1}) $$ +!< +! integrate_order_2 !< Integrate integrand field. + contains + ! deferred methods + procedure, pass(self) :: description !< Return pretty-printed object description. + procedure, pass(lhs) :: integr_assign_integr !< Operator `=`. + ! public methods + procedure, pass(self) :: destroy !< Destroy the integrator. + procedure, pass(self) :: init !< Initialize (create) the integrator. + procedure, pass(self) :: integrate !< Integrate integrand field. + procedure, nopass :: is_supported !< Check if the queried number of steps is supported or not. + procedure, nopass :: min_steps !< Return the minimum number of steps supported. + procedure, nopass :: max_steps !< Return the maximum number of steps supported. + procedure, pass(self) :: update_previous !< Cyclic update previous time steps. + ! private methods + procedure, pass(self), private :: integrate_order_2 !< Integrate integrand field by 2nd order formula. +endtype integrator_lmm_ssp_vss + +abstract interface + subroutine integrate_interface(self, U, previous, Dt, t, autoupdate) + !< Integrate field with LMM-SSP class scheme. + import :: integrand, integrator_lmm_ssp_vss, R_P + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + class(integrand), intent(inout) :: U !< Field to be integrated. + class(integrand), intent(inout) :: previous(1:) !< Previous time steps solutions of integrand field. + real(R_P), intent(inout) :: Dt(:) !< Time steps. + real(R_P), intent(in) :: t(:) !< Times. + logical, optional, intent(in) :: autoupdate !< Perform cyclic autoupdate of previous time steps. + endsubroutine integrate_interface +endinterface + +contains + ! deferred methods + pure function description(self, prefix) result(desc) + !< Return a pretty-formatted object description. + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + character(*), intent(in), optional :: prefix !< Prefixing string. + character(len=:), allocatable :: desc !< Description. + character(len=:), allocatable :: prefix_ !< Prefixing string, local variable. + character(len=1), parameter :: NL=new_line('a') !< New line character. + + prefix_ = '' ; if (present(prefix)) prefix_ = prefix + desc = '' + desc = desc//prefix_//'Strong Stability preserving Linear-Multistep-Methods Variable Stepsize class'//NL + desc = desc//prefix_//' Supported steps numbers: ['//trim(adjustl(supported_steps))//']'//NL + desc = desc//prefix_//' Supported orders: ['//trim(adjustl(supported_orders))//']'//NL + desc = desc//prefix_//' Not all steps/order combinations are supported' + endfunction description + + pure subroutine integr_assign_integr(lhs, rhs) + !< Operator `=`. + class(integrator_lmm_ssp_vss), 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_lmm_ssp_vss) + lhs%steps = rhs%steps + lhs%order = rhs%order + if (associated(rhs%integrate_)) lhs%integrate_ => rhs%integrate_ + endselect + endsubroutine integr_assign_integr + + ! public methods + elemental subroutine destroy(self) + !< Destroy the integrator. + class(integrator_lmm_ssp_vss), intent(inout) :: self !< Integrator. + + call self%destroy_abstract + self%steps = 0 + self%order = 0 + self%integrate_ => integrate_order_2 + endsubroutine destroy + + subroutine init(self, steps, order, stop_on_fail) + !< Create the actual LMM-SSP-VSS integrator. + !< + !< @note If the integrator is initialized with a bad (unsupported) number of required time steps the initialization fails and + !< the integrator error status is updated consistently for external-provided errors handling. + class(integrator_lmm_ssp_vss), intent(inout) :: self !< Integrator. + integer(I_P), intent(in) :: steps !< Number of time steps used. + integer(I_P), intent(in) :: order !< Order of accuracy. + logical, intent(in), optional :: stop_on_fail !< Stop execution if initialization fail. + + if (self%is_supported(steps=steps, order=order)) then + call self%destroy + self%steps = steps + self%order = order + select case(order) + case(2) ! LMM-SSP-VSS(steps,2) + self%integrate_ => integrate_order_2 + case(3) ! LMM-SSP-VSS(steps,3) + self%integrate_ => integrate_order_3 + endselect + else + call self%trigger_error(error=ERROR_INTEGRATOR_INIT_FAIL, & + error_message='bad (unsupported) scheme'//new_line('a')//self%description(), & + is_severe=stop_on_fail) + endif + endsubroutine init + + subroutine integrate(self, U, previous, Dt, t, autoupdate) + !< Integrate field with LMM-SSP class scheme. + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + class(integrand), intent(inout) :: U !< Field to be integrated. + class(integrand), intent(inout) :: previous(1:) !< Previous time steps solutions of integrand field. + real(R_P), intent(inout) :: Dt(:) !< Time steps. + real(R_P), intent(in) :: t(:) !< Times. + logical, optional, intent(in) :: autoupdate !< Perform cyclic autoupdate of previous time steps. + + call self%integrate_(U=U, previous=previous, Dt=Dt, t=t, autoupdate=autoupdate) + endsubroutine integrate + + elemental function is_supported(steps, order) + !< Check if the queried number of steps is supported or not. + integer(I_P), intent(in) :: steps !< Number of time steps used. + integer(I_P), intent(in) :: order !< Order of accuracy. + logical :: is_supported !< Is true if the steps number is in *supported_steps*. + + is_supported = is_admissible(n=steps, adm_range=trim(supported_steps)) + if (is_supported) is_supported = is_admissible(n=order, adm_range=trim(supported_orders)) + if (is_supported) then + select case(order) + case(2) + is_supported = steps >= 2 + case(3) + is_supported = steps >= 4 + endselect + endif + endfunction is_supported + + pure function min_steps() + !< Return the minimum number of steps supported. + integer(I_P) :: min_steps !< Minimum number of steps supported. + + min_steps = supported_steps_(1) + endfunction min_steps + + pure function max_steps() + !< Return the maximum number of steps supported. + integer(I_P) :: max_steps !< Maximum number of steps supported. + + max_steps = supported_steps_(1) + endfunction max_steps + + subroutine update_previous(self, U, previous, Dt) + !< Cyclic update previous time steps. + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + class(integrand), intent(in) :: U !< Field to be integrated. + class(integrand), intent(inout) :: previous(1:) !< Previous time steps solutions of integrand field. + real(R_P), intent(inout) :: Dt(:) !< Time steps. + integer(I_P) :: s !< Steps counter. + + do s=1, self%steps - 1 + previous(s) = previous(s + 1) + Dt(s) = Dt(s + 1) + enddo + previous(self%steps) = U + endsubroutine update_previous + + ! private methods + subroutine integrate_order_2(self, U, previous, Dt, t, autoupdate) + !< Integrate field with LMM-SSP-VSS 2nd order class scheme. + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + class(integrand), intent(inout) :: U !< Field to be integrated. + class(integrand), intent(inout) :: previous(1:) !< Previous time steps solutions of integrand field. + real(R_P), intent(inout) :: Dt(:) !< Time steps. + real(R_P), intent(in) :: t(:) !< Times. + logical, optional, intent(in) :: autoupdate !< Perform cyclic autoupdate of previous time steps. + logical :: autoupdate_ !< Perform cyclic autoupdate of previous time steps, dummy var. + integer(I_P) :: s !< Steps counter. + real(R_P) :: omega_ !< Omega coefficient. + real(R_P) :: omega_sq !< Square of omega coefficient. + + autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate + omega_= omega(Dt=Dt, s=self%steps-1) + omega_sq = omega_ * omega_ + U = previous(1) * (1._R_P / omega_sq) + previous(self%steps) * ((omega_sq - 1._R_P) / omega_sq) + & + previous(self%steps)%t(t=t(self%steps)) * (Dt(self%steps) * (omega_ + 1._R_P) / omega_) + if (autoupdate_) call self%update_previous(U=U, previous=previous, Dt=Dt) + endsubroutine integrate_order_2 + + subroutine integrate_order_3(self, U, previous, Dt, t, autoupdate) + !< Integrate field with LMM-SSP-VSS 3rd order class scheme. + class(integrator_lmm_ssp_vss), intent(in) :: self !< Integrator. + class(integrand), intent(inout) :: U !< Field to be integrated. + class(integrand), intent(inout) :: previous(1:) !< Previous time steps solutions of integrand field. + real(R_P), intent(inout) :: Dt(:) !< Time steps. + real(R_P), intent(in) :: t(:) !< Times. + logical, optional, intent(in) :: autoupdate !< Perform cyclic autoupdate of previous time steps. + logical :: autoupdate_ !< Perform cyclic autoupdate of previous time steps, dummy var. + integer(I_P) :: s !< Steps counter. + real(R_P) :: omega_ !< Omega coefficient. + + autoupdate_ = .true. ; if (present(autoupdate)) autoupdate_ = autoupdate + omega_= omega(Dt=Dt, s=self%steps-1) + U = previous(1) * ((3._R_P * omega_ + 2._R_P) / omega_ ** 3 ) + & + previous(self%steps) * (((omega_ + 1._R_P) ** 2) * (omega_ - 2._R_P) / omega_ ** 3) + & + previous(1)%t(t=t(1)) * (Dt(self%steps) * (omega_ + 1._R_P) / omega_ ** 2) + & + previous(self%steps)%t(t=t(self%steps)) * (Dt(self%steps) * (omega_ + 1._R_P) ** 2 / omega_ ** 2) + if (autoupdate_) call self%update_previous(U=U, previous=previous, Dt=Dt) + endsubroutine integrate_order_3 + + ! private non TBP + pure function dt_ratio(Dt, s) result(ratio) + !< Return `Dt(n+s)/Dt(n+Ns)` ratio. + real(R_P), intent(in) :: Dt(:) !< Time steps. + integer(I_P), intent(in) :: s !< Step index. + real(R_P) :: ratio !< Time steps ratio. + + ratio = Dt(s)/Dt(ubound(Dt, dim=1)) + endfunction dt_ratio + + pure function omega(Dt, s) + !< Return `omega=sum(dt_ratio(i)), i=1, s`. + real(R_P), intent(in) :: Dt(:) !< Time steps. + integer(I_P), intent(in) :: s !< Step index. + real(R_P) :: omega !< Omega sum. + integer(I_P) :: i !< Counter. + + omega = 0._R_P + do i=1, s + omega = omega + dt_Ratio(Dt=Dt, s=i) + enddo + endfunction omega +endmodule foodie_integrator_lmm_ssp_vss diff --git a/src/tests/accuracy/oscillation/oscillation.f90 b/src/tests/accuracy/oscillation/oscillation.f90 index 33403a37..6862a16d 100644 --- a/src/tests/accuracy/oscillation/oscillation.f90 +++ b/src/tests/accuracy/oscillation/oscillation.f90 @@ -11,6 +11,7 @@ module oscillation_test_t integrator_euler_explicit, & integrator_leapfrog, & integrator_lmm_ssp, & + integrator_lmm_ssp_vss, & integrator_runge_kutta_emd, & integrator_runge_kutta_ls, & integrator_runge_kutta_tvd @@ -24,7 +25,7 @@ module oscillation_test_t integer, parameter :: space_dimension=2 !< Space dimensions. real(R_P), parameter :: initial_state(1:space_dimension)=[0._R_P,1._R_P] !< Initial state. -character(99), parameter :: solvers(1:12) = ["all ", & +character(99), parameter :: solvers(1:13) = ["all ", & "adams-bashforth ", & "adams-bashforth-moulton", & "adams-moulton ", & @@ -34,6 +35,7 @@ module oscillation_test_t "leapfrog ", & "leapfrog-raw ", & "lmm-ssp ", & + "lmm-ssp-vss ", & "ls-runge-kutta ", & "tvd-runge-kutta "] !< List of available solvers. @@ -54,6 +56,7 @@ module oscillation_test_t logical :: plots=.false. !< Flag for activating plots saving. logical :: results=.false. !< Flag for activating results saving. character(99) :: solver='adams-bashforth' !< Solver used. + integer(I_P) :: order=0 !< (Formal) Order of accuracy. integer(I_P), allocatable :: stages_steps(:) !< Number of stages/steps used. real(R_P), allocatable :: Dt(:) !< Time step(s) exercised. real(R_P), allocatable :: tolerance(:) !< Tolerance(s) exercised on local truncation error. @@ -108,6 +111,7 @@ subroutine set_cli() call cli%add(switch='--iterations', help='Number of iterations for implicit solvers', required=.false., act='store', def='5') call cli%add(switch='--frequency', switch_ab='-f', help='Oscillation frequency', required=.false., def='1e-4', act='store') call cli%add(switch='--ss', nargs='+', help='Stages/steps used', required=.false., def='-1', act='store') + call cli%add(switch='--order', switch_ab='-o', help='ODE solver order', required=.false., def='1', act='store') call cli%add(switch='--time_step', switch_ab='-Dt', nargs='+', help='Time step', required=.false., def='100.d0', act='store') call cli%add(switch='--tolerance', switch_ab='-tol', nargs='+', help='Error Tolerance', required=.false., def='0.001d0', & act='store') @@ -128,6 +132,7 @@ subroutine parse_cli() call self%cli%get(switch='--iterations', val=self%implicit_iterations, error=self%error) ; if (self%error/=0) stop call self%cli%get(switch='-f', val=self%frequency, error=self%error) ; if (self%error/=0) stop call self%cli%get_varying(switch='--ss', val=self%stages_steps, error=self%error) ; if (self%error/=0) stop + call self%cli%get(switch='-o', val=self%order, error=self%error) ; if (self%error/=0) stop call self%cli%get_varying(switch='-Dt', val=self%Dt, error=self%error) ; if (self%error/=0) stop call self%cli%get_varying(switch='-tol', val=self%tolerance, error=self%error) ; if (self%error/=0) stop call self%cli%get(switch='-tf', val=self%final_time, error=self%error) ; if (self%error/=0) stop @@ -163,7 +168,7 @@ function is_solver_valid() integer(I_P) :: s !< Counter. is_solver_valid = .false. - do s=1, ubound(solvers, dim=1) + do s=lbound(solvers, dim=1), ubound(solvers, dim=1) is_solver_valid = (trim(adjustl(self%solver))==trim(adjustl(solvers(s)))) if (is_solver_valid) exit enddo @@ -187,7 +192,7 @@ function list_solvers() result(list) integer(I_P) :: s !< Solvers counter. list = 'Valid solver names are:' // new_line('a') - do s=1, ubound(solvers, dim=1) + do s=lbound(solvers, dim=1), ubound(solvers, dim=1) list = list // ' + ' // trim(adjustl(solvers(s))) // new_line('a') enddo endfunction list_solvers @@ -198,16 +203,17 @@ subroutine test(self, solver) class(oscillation_test), intent(in) :: self !< Test. character(*), intent(in) :: solver !< Selected solver. ! FOODIE integrators - type(integrator_adams_bashforth) :: ab_integrator !< Adams-Bashforth integrator. - type(integrator_adams_bashforth_moulton) :: abm_integrator !< Adams-Bashforth-Moulton integrator. - type(integrator_adams_moulton) :: am_integrator !< Adams-Moulton integrator. - type(integrator_back_df) :: bdf_integrator !< BDF integrator. - type(integrator_runge_kutta_emd) :: emd_rk_integrator !< Runge-Kutta integrator. - type(integrator_euler_explicit) :: euler_integrator !< Euler integrator. - type(integrator_leapfrog) :: lf_integrator !< Leapfrog integrator. - type(integrator_lmm_ssp) :: lmm_ssp_integrator !< LMM SSP integrator. - type(integrator_runge_kutta_ls) :: ls_rk_integrator !< Low Storage Runge-Kutta integrator. - type(integrator_runge_kutta_tvd) :: tvd_rk_integrator !< TVD Runge-Kutta integrator. + type(integrator_adams_bashforth) :: ab_integrator !< Adams-Bashforth integrator. + type(integrator_adams_bashforth_moulton) :: abm_integrator !< Adams-Bashforth-Moulton integrator. + type(integrator_adams_moulton) :: am_integrator !< Adams-Moulton integrator. + type(integrator_back_df) :: bdf_integrator !< BDF integrator. + type(integrator_runge_kutta_emd) :: emd_rk_integrator !< Runge-Kutta integrator. + type(integrator_euler_explicit) :: euler_integrator !< Euler integrator. + type(integrator_leapfrog) :: lf_integrator !< Leapfrog integrator. + type(integrator_lmm_ssp) :: lmm_ssp_integrator !< LMM SSP integrator. + type(integrator_lmm_ssp_vss) :: lmm_ssp_vss_integrator !< LMM SSP VSS integrator. + type(integrator_runge_kutta_ls) :: ls_rk_integrator !< Low Storage Runge-Kutta integrator. + type(integrator_runge_kutta_tvd) :: tvd_rk_integrator !< TVD Runge-Kutta integrator. ! Auxiliary variables real(R_P), allocatable :: solution(:,:) !< Solution at each time step. real(R_P), allocatable :: error(:,:) !< Error (norm L2) with respect the exact solution. @@ -246,6 +252,8 @@ subroutine test(self, solver) stages_steps_range = [lf_integrator%min_steps(), lf_integrator%max_steps()] case('lmm-ssp') stages_steps_range = [lmm_ssp_integrator%min_steps(), lmm_ssp_integrator%max_steps()] + case('lmm-ssp-vss') + stages_steps_range = [lmm_ssp_vss_integrator%min_steps(), lmm_ssp_vss_integrator%max_steps()] case('ls-runge-kutta') stages_steps_range = [ls_rk_integrator%min_stages(), ls_rk_integrator%max_stages()] case('tvd-runge-kutta') @@ -284,6 +292,7 @@ subroutine test(self, solver) frequency=self%frequency, & final_time=self%final_time, & stages_steps=s, & + order=self%order, & iterations=self%implicit_iterations, & solution=solution, & error=error(:, t), & @@ -300,11 +309,11 @@ subroutine test(self, solver) order(:, t-1) = observed_order(error=error(:, t-1:t), Dt=Dt_mean(t-1:t)) print "(A,F10.2,A,F10.2)", " Observed order, O(x): ", order(1, t-1), ", O(y): " , order(2, t-1) endif - call save_results(results=self%results, & - plots=self%plots, & - output_cli=self%output_cli, & - solver=trim(adjustl(solver))//'-'//trim(str(s, .true.)), & - frequency=self%frequency, & + call save_results(results=self%results, & + plots=self%plots, & + output_cli=self%output_cli, & + solver=trim(adjustl(solver))//'-'//trim(str(s, .true.))//'-order_'//trim(str(self%order, .true.)), & + frequency=self%frequency, & solution=solution(:, 0:last_step)) endif enddo @@ -314,6 +323,7 @@ subroutine test(self, solver) frequency=self%frequency, & final_time=self%final_time, & stages_steps=s, & + order=self%order, & iterations=self%implicit_iterations, & solution=solution, & error=error(:, t), & @@ -329,11 +339,11 @@ subroutine test(self, solver) order(:, t-1) = observed_order(error=error(:, t-1:t), Dt=self%Dt(t-1:t)) print "(A,F10.2,A,F10.2)", " Observed order, O(x): ", order(1, t-1), ", O(y): " , order(2, t-1) endif - call save_results(results=self%results, & - plots=self%plots, & - output_cli=self%output_cli, & - solver=trim(adjustl(solver))//'-'//trim(str(s, .true.)), & - frequency=self%frequency, & + call save_results(results=self%results, & + plots=self%plots, & + output_cli=self%output_cli, & + solver=trim(adjustl(solver))//'-'//trim(str(s, .true.))//'-order_'//trim(str(self%order, .true.)), & + frequency=self%frequency, & solution=solution(:, 0:last_step)) endif enddo @@ -342,7 +352,7 @@ subroutine test(self, solver) endsubroutine test ! non type bound procedures - subroutine solve(solver, frequency, final_time, stages_steps, iterations, solution, error, last_step, Dt, tolerance) + subroutine solve(solver, frequency, final_time, stages_steps, order, iterations, solution, error, last_step, Dt, tolerance) !< Rune the solver selected. !< !< The actual solver is selected by means of the *solver* input string that must be a valid string as defined into *solvers* @@ -351,6 +361,7 @@ subroutine solve(solver, frequency, final_time, stages_steps, iterations, soluti real(R_P), intent(in) :: frequency !< Oscillation frequency. real(R_P), intent(in) :: final_time !< Final integration time. integer(I_P), intent(in) :: stages_steps !< Number of stages/steps used. + integer(I_P), intent(in) :: order !< (Formal) order of accuracy. integer(I_P), intent(in) :: iterations !< Number of fixed point iterations. real(R_P), allocatable, intent(out) :: solution(:,:) !< Solution at each time step, X-Y. real(R_P), intent(out) :: error(1:) !< Error (norm L2) with respect the exact solution. @@ -358,16 +369,17 @@ subroutine solve(solver, frequency, final_time, stages_steps, iterations, soluti real(R_P), optional, intent(in) :: Dt !< Time step. real(R_P), optional, intent(in) :: tolerance !< Local error tolerance. ! FOODIE integrators - type(integrator_adams_bashforth) :: ab_integrator !< Adams-Bashforth integrator. - type(integrator_adams_bashforth_moulton) :: abm_integrator !< Adams-Bashforth-Moulton integrator. - type(integrator_adams_moulton) :: am_integrator !< Adams-Moulton integrator. - type(integrator_back_df) :: bdf_integrator !< BDF integrator. - type(integrator_runge_kutta_emd) :: emd_rk_integrator !< Runge-Kutta integrator. - type(integrator_euler_explicit) :: euler_integrator !< Euler integrator. - type(integrator_leapfrog) :: lf_integrator !< Leapfrog integrator. - type(integrator_lmm_ssp) :: lmm_ssp_integrator !< LMM SSP integrator. - type(integrator_runge_kutta_ls) :: ls_rk_integrator !< Low Storage Runge-Kutta integrator. - type(integrator_runge_kutta_tvd) :: tvd_rk_integrator !< TVD Runge-Kutta integrator. + type(integrator_adams_bashforth) :: ab_integrator !< Adams-Bashforth integrator. + type(integrator_adams_bashforth_moulton) :: abm_integrator !< Adams-Bashforth-Moulton integrator. + type(integrator_adams_moulton) :: am_integrator !< Adams-Moulton integrator. + type(integrator_back_df) :: bdf_integrator !< BDF integrator. + type(integrator_runge_kutta_emd) :: emd_rk_integrator !< Runge-Kutta integrator. + type(integrator_euler_explicit) :: euler_integrator !< Euler integrator. + type(integrator_leapfrog) :: lf_integrator !< Leapfrog integrator. + type(integrator_lmm_ssp) :: lmm_ssp_integrator !< LMM SSP integrator. + type(integrator_lmm_ssp_vss) :: lmm_ssp_vss_integrator !< LMM SSP VSS integrator. + type(integrator_runge_kutta_ls) :: ls_rk_integrator !< Low Storage Runge-Kutta integrator. + type(integrator_runge_kutta_tvd) :: tvd_rk_integrator !< TVD Runge-Kutta integrator. ! Auxiliary variables integer(I_P), parameter :: max_rk_stages=5 !< Max RK stages used to init high order multi-step solver. type(oscillation) :: oscillator !< Oscillation field. @@ -376,6 +388,7 @@ subroutine solve(solver, frequency, final_time, stages_steps, iterations, soluti integer :: step_offset !< Time steps counter offset for slicing previous data array. logical :: adaptive !< Flag for tagging time step adaptive class of solvers. real(R_P) :: Dt_a !< Adaptive time step. + real(R_P), allocatable :: Dts(:) !< Time steps for variable stepsize methods. type(oscillation), allocatable :: rk_stage(:) !< Runge-Kutta stages. type(oscillation), allocatable :: previous(:) !< Previous time steps solutions. type(oscillation) :: filter !< Filter displacement. @@ -454,6 +467,13 @@ subroutine solve(solver, frequency, final_time, stages_steps, iterations, soluti call lmm_ssp_integrator%init(steps=stages_steps) if (allocated(previous)) deallocate(previous) ; allocate(previous(1:stages_steps)) + case("lmm-ssp-vss") + supported = lmm_ssp_vss_integrator%is_supported(steps=stages_steps, order=order) + multistep = .true. + call lmm_ssp_vss_integrator%init(steps=stages_steps, order=order) + if (allocated(Dts)) deallocate(Dts) ; allocate(Dts(1:stages_steps)) ; Dts = Dt + if (allocated(previous)) deallocate(previous) ; allocate(previous(1:stages_steps)) + case("ls-runge-kutta") supported = ls_rk_integrator%is_supported(stages_steps) call ls_rk_integrator%init(stages=stages_steps) @@ -546,6 +566,11 @@ subroutine solve(solver, frequency, final_time, stages_steps, iterations, soluti previous=previous, & Dt=Dt, & t=solution(0, step-step_offset:step-1)) + case("lmm-ssp-vss") + call lmm_ssp_vss_integrator%integrate(U=oscillator, & + previous=previous, & + Dt=Dts, & + t=solution(0, step-step_offset:step-1)) endselect endif else