module test_timestepping_m use funit implicit none contains @test subroutine multistep_storage !! Test of multistep_storage_m use precision_grillix_m, only : GP use helpers_m, only : almost_equal use multistep_m, only : multistep_storage_t implicit none logical :: fine integer :: iret type(multistep_storage_t) :: mtstore integer :: ndim, nstorage, nfields, it, tt, i, istorage, ifield real(GP), allocatable, dimension(:,:) :: y real(GP), allocatable, dimension(:, :, :) :: yfirst write(*,*)'' write(*,*)'' write(*,'(A80)')'test_multistep_storage ----------------------------------------------------' ! Set dimensions and initial state ndim = 5 nstorage = 4 nfields = 2 allocate(y(ndim, nfields)) allocate(yfirst(ndim, nstorage, nfields)) do ifield = 1, nfields do istorage = 1, nstorage tt = -istorage + 1 do i = 1, ndim yfirst(i, istorage, ifield) = some_yval(tt, i, ifield) enddo enddo enddo ! Initialise storage call mtstore%init_storage(ndim, nstorage, nfields, yfirst) ! Test metadata iret = mtstore%get_ndim() @assertEqual(iret, ndim) iret = mtstore%get_nstorage() @assertEqual(iret, nstorage) iret = mtstore%get_nfields() @assertEqual(iret, nfields) ! Something like a time-loop do it = 1, 9 ! Do something to y do ifield = 1, nfields do i = 1, ndim y(i, ifield) = some_yval(it, i, ifield) enddo enddo ! Shift storage call mtstore%shift_storage(y) do ifield = 1, nfields do istorage = 1, nstorage tt = it-(istorage-1) do i = 1, ndim fine = almost_equal(mtstore%vstore(i,istorage,ifield), some_yval(tt, i, ifield) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) enddo enddo enddo enddo write(*,'(A80)')'test_multistep_storage complete -------------------------------------------' write(*,*)'' write(*,*)'' contains real(GP) function some_yval(it, i, ifield) integer, intent(in) :: it !! Something like a time-stamp integer, intent(in) :: i !! Something like a grid-index integer, intent(in) :: ifield !! Field index (e.g. y or dy) some_yval = 1.0_GP * (it + 10*(i-1) + 100*(ifield-1)) end function end subroutine @test subroutine karniadakis !! Test of module karniadakis_m !! integrates !! - dy_1/dtau = y --> y_1(tau) = exp(tau) !! - dy_2/dtau = 2*tau*y --> y_2(tau) = exp(tau^2) !! until tau = 2 !! Analytic result of y_1(2) = exp(2) ~ 7.38905609893 !! Analytic result of y_2(2) = exp(4) ~ 54.5981500331 use precision_grillix_m, only : GP use helpers_m, only : almost_equal use multistep_m, only : karniadakis_t implicit none type(karniadakis_t), allocatable :: tstep_krk logical :: fine integer :: ndim, iorder, it, iac, iret real(GP) :: tau, dtau, ret real(GP), dimension(2) :: y, yn, dy real(GP), dimension(2, 5) :: y_prev, dy_prev real(GP), dimension(2,6,2) :: res real(GP), dimension(6) :: res_dtau_bdf write(*,*)'' write(*,*)'' write(*,'(A80)')'test_karniadakis ----------------------------------------------------------' dtau = 1.0E-2_GP ndim = 2 ! Final result of integration ----------------------------------------------------------------- ! Analytic result of res(1,:,:) is exp(2) ~ 7.38905609893 ! Analytic result of res(2,:,:) is exp(4) ~ 54.5981500331 res(1,1,1) = 7.31601785182993_GP res(2,1,1) = 50.8102141064147_GP res(1,2,1) = 7.38753899776555_GP res(2,2,1) = 54.3967993224238_GP res(1,3,1) = 7.38857042247586_GP res(2,3,1) = 54.5815321294632_GP res(1,4,1) = 7.38854881287465_GP res(2,4,1) = 54.5901909943513_GP res(1,5,1) = 7.38855669091906_GP res(2,5,1) = 54.5907706710100_GP res(1,6,1) = 7.38855971714400_GP res(2,6,1) = 54.5908348471513_GP res(1,1,2) = 7.31601785182993_GP res(2,1,2) = 50.8102141064147_GP res(1,2,2) = 7.38808450026534_GP res(2,2,2) = 54.4049585707350_GP res(1,3,2) = 7.38904521846053_GP res(2,3,2) = 54.5885021250675_GP res(1,4,2) = 7.38905598343100_GP res(2,4,2) = 54.5976714191383_GP res(1,5,2) = 7.38905609773340_GP res(2,5,2) = 54.5981247196009_GP res(1,6,2) = 7.38905609891841_GP res(2,6,2) = 54.5981486517818_GP res_dtau_bdf(1) = dtau res_dtau_bdf(2) = 2.0_GP/3.0_GP*dtau res_dtau_bdf(3) = 6.0_GP/11.0_GP*dtau res_dtau_bdf(4) = 12.0_GP/25.0_GP*dtau res_dtau_bdf(5) = 60.0_GP/137.0_GP*dtau res_dtau_bdf(6) = 60.0_GP/147.0_GP*dtau ! Analytic solutions and values for initial timesteps ----------------------------------------- do it = 1,5 tau = -dtau*it y_prev(1,it) = exp(tau) dy_prev(1,it) = y_prev(1,it) y_prev(2,it) = exp(tau**2) dy_prev(2,it) = 2*tau*y_prev(2,it) enddo ! Perform integration with all 6 orders (iorder loop) and with and without initial timesteps given (iac loop) do iac = 1, 2 do iorder = 1, 6 allocate(tstep_krk) if (iac == 1) then call tstep_krk%init(ndim, iorder, dtau) elseif (iac == 2) then call tstep_krk%init(ndim, iorder, dtau, y_prev(:,1:iorder-1), dy_prev(:,1:iorder-1)) endif tau = 0.0_GP y(1) = 1.0_GP y(2) = 1.0_GP ! Time integration loop --------------------------------------------------------------- do it = 1, 200 dy(1) = y(1) ! dy/dtau = y --> y(tau) = exp(tau) dy(2) = 2*tau*y(2) ! dy/dtau = 2*tau*y --> y(tau) = exp(tau^2) call tstep_krk%advance(dy, y, yn) call tstep_krk%shift_storage((/y, dy/)) y = yn tau = tau + dtau ! Check order of last time-step and dtau_bdf iret = tstep_krk%get_current_order() ret = tstep_krk%get_dtau_bdf() if (iac == 1) then @assertEqual(iret, min(it, iorder)) fine = almost_equal(ret, res_dtau_bdf(min(it, iorder)) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) elseif (iac == 2) then @assertEqual(iret, iorder) fine = almost_equal(ret, res_dtau_bdf(iorder) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) endif enddo ! Check result from integration ------------------------------------------------------- fine = almost_equal(y(1), res(1, iorder, iac) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) fine = almost_equal(y(2), res(2, iorder, iac) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) fine = almost_equal(tau, 2.0_GP , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) ! Check Metainformation --------------------------------------------------------------- ! Check metadata if (iac == 1) then iret = tstep_krk%get_ncount() @assertEqual(iret, 200) elseif (iac == 2) then iret = tstep_krk%get_ncount() @assertEqual(iret, 200) endif iret = tstep_krk%get_ndim() @assertEqual(iret, 2) ret = tstep_krk%get_dtau() fine = almost_equal(ret, dtau , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) ret = tstep_krk%get_dtau_bdf() fine = almost_equal(ret, res_dtau_bdf(iorder) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) deallocate(tstep_krk) enddo enddo write(*,'(A80)')'test_karniadakis complete -------------------------------------------------' write(*,*)'' write(*,*)'' end subroutine @test subroutine runge_kutta !! Test of module runge_kutta_m !! integrates !! - dy_1/dtau = y --> y_1(tau) = exp(tau) !! - dy_2/dtau = 2*tau*y --> y_2(tau) = exp(tau^2) !! until tau = 2 !! Analytic result of y_1(2) = exp(2) ~ 7.38905609893 !! Analytic result of y_2(2) = exp(4) ~ 54.5981500331 use precision_grillix_m, only : GP use helpers_mesh_sample_m, only : mesh_sample_circular use helpers_m, only : almost_equal use runge_kutta_m, only : runge_kutta_t implicit none type(runge_kutta_t), allocatable :: tsteprk logical :: fine integer :: iorder, iret, ndim, it real(GP) :: tau, dtau real(GP), dimension(2) :: y real(GP), dimension(2,4) :: res write(*,*)'' write(*,*)'' write(*,'(A80)')'test_runge_kutta ----------------------------------------------------------' dtau = 1.0E-2_GP ndim = 2 y(1) = 1.0_GP y(2) = 1.0_GP res(1,1) = 7.31601785182993_GP res(2,1) = 50.8102141064147_GP res(1,2) = 7.38881164097980_GP res(2,2) = 54.5696748189694_GP res(1,3) = 7.38905548808157_GP res(2,3) = 54.5979406960639_GP res(1,4) = 7.38905609770937_GP res(2,4) = 54.5981485215352_GP do iorder = 1, 4 allocate(tsteprk) call tsteprk%init(ndim, iorder) tau = 0.0_GP do it = 1, 200 call tsteprk%advance(dtau, tau, y, right_hand_side) tau = tau + dtau enddo fine = almost_equal(y(1), res(1,iorder) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) fine = almost_equal(y(2), res(2,iorder) , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) fine = almost_equal(2.0_GP, tau , 0.0_GP, 1.0E-10_GP) @assertTrue(fine) iret = tsteprk%get_ncount() @assertEqual(iret, 200) iret = tsteprk%get_ndim() @assertEqual(iret, 2) iret = tsteprk%get_order() @assertEqual(iret, iorder) deallocate(tsteprk) tau = 0.0_GP y(1) = 1.0_GP y(2) = 1.0_GP enddo write(*,'(A80)')'test_runge_kutta complete -------------------------------------------------' write(*,*)'' write(*,*)'' contains subroutine right_hand_side(ndim, tau, y, dy) !! Some sample right hand side integer, intent(in) :: ndim real(GP), intent(in) :: tau real(GP), dimension(ndim), intent(in) :: y real(GP), dimension(ndim), intent(out) :: dy dy(1) = y(1) dy(2) = 2*tau * y(2) end subroutine end subroutine end module