Skip to content

Commit bcd7e27

Browse files
szaghiszaghi
committed
add ms-rk-ssp test, seems to work
1 parent 7c9ca89 commit bcd7e27

File tree

2 files changed

+134
-23
lines changed

2 files changed

+134
-23
lines changed

‎src/lib/foodie_integrator_ms_runge_kutta_ssp.f90‎

Lines changed: 74 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,8 @@ module foodie_integrator_ms_runge_kutta_ssp
4444
public :: integrator_ms_runge_kutta_ssp
4545

4646
character(len=99), parameter :: class_name_='ms_runge_kutta_ssp' !< Name of the class of schemes.
47-
character(len=99), parameter :: supported_schemes_(1:2)=[trim(class_name_)//'_steps_2_stages_2_order_3', &
47+
character(len=99), parameter :: supported_schemes_(1:3)=[trim(class_name_)//'_steps_2_stages_2_order_3', &
48+
trim(class_name_)//'_steps_3_stages_2_order_3', &
4849
trim(class_name_)//'_steps_4_stages_5_order_8'] !< List of supported
4950
!< schemes.
5051

@@ -54,8 +55,8 @@ module foodie_integrator_ms_runge_kutta_ssp
5455
!<
5556
!< @note The integrator must be created or initialized (initialize the *A,Ahat,B,Bhat,D,T* coefficients) before used.
5657
private
57-
integer(I_P) :: steps=0 !< Number of time steps.
58-
integer(I_P) :: stages=0 !< Number of stages.
58+
integer(I_P), public :: steps=0 !< Number of time steps.
59+
integer(I_P), public :: stages=0 !< Number of stages.
5960
real(R_P), allocatable :: A(:,:) !< *A* coefficients.
6061
real(R_P), allocatable :: Ahat(:,:) !< *Ahat* coefficients.
6162
real(R_P), allocatable :: B(:) !< *B* coefficients.
@@ -175,7 +176,57 @@ subroutine initialize(self, scheme, stop_on_fail)
175176
call self%destroy
176177
select case(trim(adjustl(scheme)))
177178
case('ms_runge_kutta_ssp_steps_2_stages_2_order_3')
178-
error stop 'error: "ms_runge_kutta_ssp_steps_2_stages_2_order_3" to be implemented!'
179+
self%steps = 2
180+
self%stages = 2
181+
182+
allocate(self%A(1:self%stages, 1:self%stages)) ; self%A = 0._R_P
183+
self%A(2, 1) = 0.910683602522959_R_P
184+
185+
allocate(self%Ahat(1:self%stages, 1:self%steps)) ; self%Ahat = 0._R_P
186+
self%Ahat(2, 1) = -1.11985107194706e-19_R_P
187+
188+
allocate(self%B(1:self%stages)) ; self%B = 0._R_P
189+
self%B(1) = 0.535898384862245_R_P
190+
self%B(2) = 0.803847577293368_R_P
191+
192+
allocate(self%Bhat(1:self%steps)) ; self%Bhat = 0._R_P
193+
self%Bhat(1) = 0.267949192431123_R_P
194+
195+
allocate(self%D(1:self%stages, 1:self%steps)) ; self%D = 0._R_P
196+
self%D(2, 1) = 1._R_P / 3._R_P
197+
self%D(2, 2) = 2._R_P / 3._R_P
198+
199+
allocate(self%T(1:self%steps)) ; self%T = 0._R_P
200+
self%T(1) = 0.607695154586736_R_P
201+
self%T(2) = 0.392304845413264_R_P
202+
203+
case('ms_runge_kutta_ssp_steps_3_stages_2_order_3')
204+
self%steps = 3
205+
self%stages = 2
206+
207+
allocate(self%A(1:self%stages, 1:self%stages)) ; self%A = 0._R_P
208+
self%A(2, 1) = 0.731058363135786_R_P
209+
210+
allocate(self%Ahat(1:self%stages, 1:self%steps)) ; self%Ahat = 0._R_P
211+
self%Ahat(2, 1) = 0.127467809251820_R_P
212+
213+
allocate(self%B(1:self%stages)) ; self%B = 0._R_P
214+
self%B(1) = 0.618048297723782_R_P
215+
self%B(2) = 0.759677988437936_R_P
216+
217+
allocate(self%Bhat(1:self%steps)) ; self%Bhat = 0._R_P
218+
self%Bhat(1) = 0.246670340394148_R_P
219+
220+
allocate(self%D(1:self%stages, 1:self%steps)) ; self%D = 0._R_P
221+
self%D(2, 1) = 0.186433848116852_R_P
222+
self%D(2, 2) = 1.80945758995975e-24_R_P
223+
self%D(2, 3) = 0.813566151883148_R_P
224+
225+
allocate(self%T(1:self%steps)) ; self%T = 0._R_P
226+
self%T(1) = 0.312198313277933_R_P
227+
self%T(2) = 2.58493941422821e-24_R_P
228+
self%T(3) = 0.687801686722067_R_P
229+
179230
case('ms_runge_kutta_ssp_steps_4_stages_5_order_8')
180231
self%steps = 4
181232
self%stages = 5
@@ -223,25 +274,25 @@ subroutine initialize(self, scheme, stop_on_fail)
223274
self%Bhat(3) = 0.541410372692778_R_P
224275

225276
allocate(self%D(1:self%stages, 1:self%steps)) ; self%D = 0._R_P
226-
self%D(1, 1) = 0.0273928990974108_R_P
227-
self%D(2, 1) = 0.283607987548794_R_P
228-
self%D(3, 1) = 0.0642241421937960_R_P
229-
self%D(4, 1) = 0.194681462814288_R_P
230-
231-
self%D(1, 2) = 0.554039201229659_R_P
232-
self%D(2, 2) = 0.0914454235177934_R_P
233-
self%D(3, 2) = 0.198601843371796_R_P
234-
self%D(4, 2) = 0.293086617882372_R_P
235-
236-
self%D(1, 3) = 0.348927848402249_R_P
237-
self%D(2, 3) = 0.437897855084625_R_P
238-
self%D(3, 3) = 0.662266498804591_R_P
239-
self%D(4, 3) = 0.158740367819382_R_P
240-
241-
self%D(1, 4) = 0.0696400512706807_R_P
242-
self%D(2, 4) = 0.187048733848788_R_P
243-
self%D(3, 4) = 0.0749075156298171_R_P
244-
self%D(4, 4) = 0.353491551483958_R_P
277+
self%D(2, 1) = 0.0273928990974108_R_P
278+
self%D(3, 1) = 0.283607987548794_R_P
279+
self%D(4, 1) = 0.0642241421937960_R_P
280+
self%D(5, 1) = 0.194681462814288_R_P
281+
282+
self%D(2, 2) = 0.554039201229659_R_P
283+
self%D(3, 2) = 0.0914454235177934_R_P
284+
self%D(4, 2) = 0.198601843371796_R_P
285+
self%D(5, 2) = 0.293086617882372_R_P
286+
287+
self%D(2, 3) = 0.348927848402249_R_P
288+
self%D(3, 3) = 0.437897855084625_R_P
289+
self%D(4, 3) = 0.662266498804591_R_P
290+
self%D(5, 3) = 0.158740367819382_R_P
291+
292+
self%D(2, 4) = 0.0696400512706807_R_P
293+
self%D(3, 4) = 0.187048733848788_R_P
294+
self%D(4, 4) = 0.0749075156298171_R_P
295+
self%D(5, 4) = 0.353491551483958_R_P
245296

246297
allocate(self%T(1:self%steps)) ; self%T = 0._R_P
247298
self%T(1) = 0.0273988434707855_R_P

‎src/tests/accuracy/oscillation/oscillation.f90‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@ module oscillation_test_t
1515
integrator_leapfrog, &
1616
integrator_lmm_ssp, &
1717
integrator_lmm_ssp_vss, &
18+
integrator_ms_runge_kutta_ssp, &
1819
integrator_object, &
1920
integrator_runge_kutta_emd, &
2021
integrator_runge_kutta_ls, &
@@ -265,6 +266,7 @@ subroutine initialize_test
265266
type(integrator_leapfrog) :: int_leapfrog !< Leapfrog integrator.
266267
type(integrator_lmm_ssp) :: int_lmm_ssp !< LMM SSP integrator.
267268
type(integrator_lmm_ssp_vss) :: int_lmm_ssp_vss !< LMM SSP VSS integrator.
269+
type(integrator_ms_runge_kutta_ssp) :: int_ms_runge_kutta_ssp !< Multistep Runge Kutta SSP integrator.
268270
type(integrator_runge_kutta_ls) :: int_runge_kutta_ls !< Runge Kutta low storage integrator.
269271
type(integrator_runge_kutta_ssp) :: int_runge_kutta_ssp !< Runge Kutta SSP integrator.
270272

@@ -284,6 +286,8 @@ subroutine initialize_test
284286
integrate => integrate_lmm_ssp_vss
285287
elseif (index(trim(adjustl(scheme)), trim(int_lmm_ssp%class_name())) > 0) then
286288
integrate => integrate_lmm_ssp
289+
elseif (index(trim(adjustl(scheme)), trim(int_ms_runge_kutta_ssp%class_name())) > 0) then
290+
integrate => integrate_ms_runge_kutta_ssp
287291
elseif (index(trim(adjustl(scheme)), trim(int_runge_kutta_emd%class_name())) > 0) then
288292
integrate => integrate_runge_kutta_emd
289293
elseif (index(trim(adjustl(scheme)), trim(int_runge_kutta_ls%class_name())) > 0) then
@@ -751,6 +755,62 @@ subroutine integrate_lmm_ssp_vss(scheme, frequency, final_time, solution, error,
751755
error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
752756
endsubroutine integrate_lmm_ssp_vss
753757

758+
subroutine integrate_ms_runge_kutta_ssp(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, tolerance)
759+
!< Integrate domain by means of the Multistep Runge-Kutta SSP scheme.
760+
character(*), intent(in) :: scheme !< Selected scheme.
761+
real(R_P), intent(in) :: frequency !< Oscillation frequency.
762+
real(R_P), intent(in) :: final_time !< Final integration time.
763+
real(R_P), allocatable, intent(out) :: solution(:,:) !< Solution at each time step, X-Y.
764+
real(R_P), intent(out) :: error(1:) !< Error (norm L2) with respect the exact solution.
765+
integer(I_P), intent(out) :: last_step !< Last time step computed.
766+
integer(I_P), intent(in), optional :: iterations !< Number of fixed point iterations.
767+
real(R_P), intent(in), optional :: Dt !< Time step.
768+
real(R_P), intent(in), optional :: tolerance !< Local error tolerance.
769+
type(integrator_ms_runge_kutta_ssp) :: integrator !< The integrator.
770+
type(integrator_runge_kutta_ssp) :: integrator_rk !< RK integrator for starting non self-starting integrators.
771+
type(oscillator) :: domain !< Oscillation field.
772+
type(oscillator), allocatable :: rk_stage(:) !< Runge-Kutta stages for the RK initial integrator.
773+
type(oscillator), allocatable :: stage(:) !< Runge-Kutta stages.
774+
type(oscillator), allocatable :: previous(:) !< Previous time steps solutions.
775+
integer :: step !< Time steps counter.
776+
777+
call domain%init(initial_state=initial_state, frequency=frequency)
778+
779+
if (allocated(solution)) deallocate(solution) ; allocate(solution(0:space_dimension, 0:int(final_time/Dt)))
780+
solution = 0.0_R_P
781+
solution(1:, 0) = domain%output()
782+
783+
call integrator%initialize(scheme=scheme)
784+
if (allocated(previous)) deallocate(previous) ; allocate(previous(1:integrator%steps))
785+
if (allocated(stage)) deallocate(stage) ; allocate(stage(1:integrator%stages))
786+
787+
call integrator_rk%initialize(scheme='runge_kutta_ssp_stages_5_order_4')
788+
if (allocated(rk_stage)) deallocate(rk_stage) ; allocate(rk_stage(1:integrator_rk%stages))
789+
790+
step = 0
791+
do while(solution(0, step) < final_time .and. step < ubound(solution, dim=2))
792+
step = step + 1
793+
794+
if (integrator%steps >= step) then
795+
call integrator_rk%integrate(U=domain, stage=rk_stage, Dt=Dt, t=solution(0, step))
796+
previous(step) = domain
797+
else
798+
call integrator%integrate(U=domain, &
799+
previous=previous, &
800+
stage=stage, &
801+
Dt=Dt, &
802+
t=solution(0, step-integrator%steps:step-1))
803+
endif
804+
805+
solution(0, step) = step * Dt
806+
807+
solution(1:, step) = domain%output()
808+
enddo
809+
last_step = step
810+
811+
error = error_L2(frequency=frequency, solution=solution(:, 0:last_step))
812+
endsubroutine integrate_ms_runge_kutta_ssp
813+
754814
subroutine integrate_runge_kutta_emd(scheme, frequency, final_time, solution, error, last_step, iterations, Dt, tolerance)
755815
!< Integrate domain by means of the Runge Kutta embedded scheme.
756816
character(*), intent(in) :: scheme !< Selected scheme.

0 commit comments

Comments
 (0)