From 0c8d8fac9c499374c77787023b4d10b583dcb244 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sat, 16 Oct 2021 19:03:25 -0300 Subject: [PATCH] EVOLUTION: Rewrite SSPRK4(2)M integration method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 330 +++++++++++++++++++----------------------- 1 file changed, 146 insertions(+), 184 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index f3339de..20397f2 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -276,11 +276,11 @@ module evolution betas(:) = - betas(:) / 3.0d+00 - case ("rk3.m", "ssprk(m,3)", "SSPRK(m,3)") + case ("SSPRK3(2)m", "SSP3(2)m", "ssprk3(2)m", "ssp3(2)m", "SSPRK(m,3)", "ssprk(m,3)") error_control = .true. - evolve => evolve_ssprk3_m - registers = 3 + evolve => evolve_ssprk32m + registers = 4 n = 2 do while(n**2 <= stages) n = n + 1 @@ -288,7 +288,25 @@ module evolution n = n - 1 stages = max(4, n**2) cfl = (n - 1) * n * cfl - write(name_int, "('Optimal 3rd order SSPRK(',i0,',3) with error control')") stages + write(name_int, "('3ʳᵈ-order order ',i0,'-step embedded SSPRK3(2)n² method')") stages + + betas(:) = [ 5.5d-01, -2.7d-01, 5.0d-02 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 3.0d+00 + + if (stages == 4 .and. verbose) then + write(*,*) + write(*,*) "ADVICE:" + write(*,*) "The method 'SSPRK3(2)m' is primarily designed for n²," // & + " where n > 2. Since you intend to use 4 stages with" // & + " this method, it is highly advised to use method" // & + " 'SSPRK3(2)4', since it slightly better optimized and" // & + " significantly used less memory." + end if case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") @@ -302,11 +320,13 @@ module evolution if (verbose) then write(*,*) write(*,"(1x,a)") "ERROR!" - write(*,"(1x,a)") "The selected time advance method is not " // & + write(*,"(1x,a)") "The selected time advance method is not " // & "implemented: " // trim(integration) - write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // & - " 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // & + write(*,"(1x,a)") "Available methods: 'euler' 'rk2', 'rk3'," // & + " 'ssprk(m,2)', 'ssprk(4,3)', 'ssprk(5,3)'," // & " 'ssprk(m,3)', 'ssprk(10,4)'." + write(*,"(1x,a)") "Available embedded methods: 'ssprk3(2)4'," // & + " 'ssprk3(2)m'." end if status = 1 @@ -1943,30 +1963,31 @@ module evolution ! !=============================================================================== ! -! subroutine EVOLVE_SSPRK3_M: +! subroutine EVOLVE_SSPRK32M: ! -------------------------- ! -! Subroutine advances the solution by one time step using the 3rd order -! m-stage Strong Stability Preserving Runge-Kutta time integration method. +! Subroutine advances the solution by one time step using the 3ʳᵈ-order +! m-stage (m=n²) embedded Strong Stability Preserving Runge-Kutta time +! integration method SSP3(2)m with the error control. ! ! References: ! -! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., -! "Strong stability preserving Runge-Kutta and multistep -! time discretizations", -! World Scientific Publishing, 2011 +! [1] Conde, S., Fekete, I., Shadid, J. N., +! "Embedded error estimation and adaptive step-size control for +! optimal explicit strong stability preserving Runge–Kutta methods" +! 2018, arXiv:1806.08693v1 +! [2] Gottlieb, S., Ketcheson, D., Shu, C.-W., +! "Strong stability preserving Runge-Kutta and multistep time +! discretizations", +! 2011, World Scientific Publishing ! !=============================================================================== ! - subroutine evolve_ssprk3_m() + subroutine evolve_ssprk32m() - use blocks , only : block_data, list_data, get_dblocks + use blocks , only : block_data, list_data use boundaries , only : boundary_fluxes - use coordinates, only : nc => ncells, nb, ne - use equations , only : errors, nf, ibp, cmax -#ifdef MPI - use mpitools , only : reduce_maximum -#endif /* MPI */ + use equations , only : errors, ibp, cmax use sources , only : update_sources implicit none @@ -1974,18 +1995,12 @@ module evolution type(block_data), pointer :: pdata logical , save :: first = .true. - integer , save :: n, n1, n2, n3, n4 - real(kind=8), save :: r, f1, f2, g1, g2 + integer , save :: i1, i2, i3, i4, n, m + real(kind=8), save :: f1, f2, f3, f4 logical :: test - integer :: i, l, nrej - real(kind=8) :: tm, dtm, ds, umax, emax - real(kind=8) :: fc, fcmn, fcmx - - real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr - - real(kind=8), parameter :: k1 = -5.8d-01 / 3.0d+00, k2 = 7.0d-02, & - k3 = -1.0d-01 / 3.0d+00 + integer :: nrej, i + real(kind=8) :: tm, dtm, dh, fc ! !------------------------------------------------------------------------------- ! @@ -1995,45 +2010,45 @@ module evolution if (first) then - n = int(sqrt(1.0d+00 * stages)) - n1 = (n - 1) * (n - 2) / 2 - n2 = n1 + 1 - n3 = n * (n + 1) / 2 - 1 - n4 = n * (n + 1) / 2 + 1 + n = 2 + do while(n**2 < stages) + n = n + 1 + end do + m = stages - n + i1 = (n - 1) * (n - 2) / 2 + i2 = i1 + 1 + i3 = n * (n + 1) / 2 + i4 = i3 + 1 - r = (1.0d+00 * (2 * n - 1)) - f1 = (1.0d+00 * n ) / r - f2 = (1.0d+00 * (n - 1)) / r - r = 1.0d+00 * (stages - n) - - g1 = 1.0d+00 / ( stages - n) - 1.0d+00 / stages - g2 = 1.0d+00 / (2 * stages - n) - 1.0d+00 / stages + fc = real(2 * n - 1, kind=8) + f1 = real(n - 1, kind=8) / fc + f2 = real(n , kind=8) / fc + f4 = real(m , kind=8) / stages + f3 = 1.0d+00 - f4 first = .false. end if - fc = fac - fcmn = facmin - fcmx = facmax - - l = get_dblocks() -#if NDIMS == 3 - allocate(lerr(l,nf,nc,nc,nc)) -#else /* NDIMS == 3 */ - allocate(lerr(l,nf,nc,nc, 1)) -#endif /* NDIMS == 3 */ - test = .true. nrej = 0 +! at the entry point we assume the previous solution of conserved variables U(n) +! is stored in pdata%u(:,:,:,:,1) and the primitive variables are already +! updated from this array; +! +! pdata%uu(:,:,:,:,1) - the previous solution, U(n) +! pdata%uu(:,:,:,:,2) - the new 3rd order solution, U(1) +! pdata%uu(:,:,:,:,3) - the intermediate solution, U(2) +! pdata%uu(:,:,:,:,4) - the new 2nd order solution, U(3) +! do while(test) - lerr(:,:,:,:,:) = 0.0d+00 +! initiate the fractional time step +! + dh = dt / m - ds = dt / r - -!= 1st step: U(1) = U(n) +!= preparation step: U(1) = U(n) ! pdata => list_data do while (associated(pdata)) @@ -2045,12 +2060,12 @@ module evolution pdata => pdata%next end do -!= 2ns step: U(1) = U(1) + dt/r F[U(1)], for i = 1, ..., (n-1)*(n-2)/2 +!= 1st step: U(1) = U(1) + dt/r F[U(1)], for i = 1...(n-1)*(n-2)/2 ! - do i = 1, n1 + do i = 1, i1 - tm = time + i * ds - dtm = ds + tm = time + (i - 1) * dh + dtm = dh pdata => list_data do while (associated(pdata)) @@ -2063,29 +2078,19 @@ module evolution call boundary_fluxes() - l = 1 pdata => list_data do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) - end do ! n = 1, n1 + end do ! i = 1, i1 -!= 3rd step: U(2) = U(1) -! -! iterate over all data blocks +!= 2nd step: U(2) = U(1) ! pdata => list_data do while (associated(pdata)) @@ -2093,15 +2098,14 @@ module evolution pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,2) pdata => pdata%next - end do -!= 4th step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1, ..., n*(n+1)/2-1 +!= 3rd step: U(1) = U(1) + dt/r F[U(1)], for i = (n-1)*(n-2)/2+1...n*(n+1)/2 ! - do i = n2, n3 + do i = i2, i3 - tm = time + i * ds - dtm = ds + tm = time + (i - 1) * dh + dtm = dh pdata => list_data do while (associated(pdata)) @@ -2114,68 +2118,43 @@ module evolution call boundary_fluxes() - l = 1 pdata => list_data do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) - end do ! i = n2, n3 + end do ! i = i2, i3 -!= 5th step: U(1) = n * U(2) + (n - 1) * ( U(1) + dt/r F[U(1)] ) / (2 * n - 1) +!= 4th step: U(3) = (1 - r/s) * U(n) + r/s * U(1), +! U(1) = [(n - 1) * U(1) + n * U(2)] / (2*n - 1), +! U(3) = U(3) - r/s * U(1) ! - tm = time + dt - dtm = ds - pdata => list_data do while (associated(pdata)) - call update_increment(pdata) - call update_sources(pdata, tm, dtm, pdata%du(:,:,:,:)) - - pdata => pdata%next - end do - - call boundary_fluxes() - - l = 1 - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,3) & - + f2 * (pdata%uu(:,:,:,:,2) & - + ds * pdata%du(:,:,:,:)) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g2 * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,4) = f3 * pdata%uu(:,:,:,:,1) & + + f4 * pdata%uu(:,:,:,:,2) + pdata%uu(:,:,:,:,2) = f1 * pdata%uu(:,:,:,:,2) & + + f2 * pdata%uu(:,:,:,:,3) + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) & + - f4 * pdata%uu(:,:,:,:,2) pdata => pdata%next end do call update_variables(tm, dtm) -!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2, ..., m +!= 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2+1...stages ! - do i = n4, stages + do i = i4, stages - tm = time + (i - n) * ds - dtm = ds + tm = time + (i - 1 - n) * dh + dtm = dh pdata => list_data do while (associated(pdata)) @@ -2188,93 +2167,90 @@ module evolution call boundary_fluxes() - l = 1 pdata => list_data do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + g1 * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + dh * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm) - end do ! i = n4, stages + end do ! i = i4, stages -! find umax +!= 7th step: U(3) = U(3) + r/s * U(1) ! - umax = 0.0d+00 pdata => list_data do while (associated(pdata)) -#if NDIMS == 3 - umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,1))), & - maxval(abs(pdata%uu(:,nb:ne,nb:ne,nb:ne,2)))) -#else /* NDIMS == 3 */ - umax = max(umax, maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,1))), & - maxval(abs(pdata%uu(:,nb:ne,nb:ne, : ,2)))) -#endif /* NDIMS == 3 */ + + pdata%uu(:,:,:,:,4) = pdata%uu(:,:,:,:,4) + f4 * pdata%uu(:,:,:,:,2) + pdata => pdata%next end do -! get the maximum error of each variable +!= 8th step: decay the magnetic divergence potential of both solutions ! - do l = 1, nf - errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) - end do + if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) -#ifdef MPI - call reduce_maximum(umax) - call reduce_maximum(errors) -#endif /* MPI */ + pdata => list_data + do while (associated(pdata)) -! calculate tolerance and time step + pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,2) + pdata%uu(ibp,:,:,:,4) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,4) + + pdata => pdata%next + end do + end if + +! update the vector of errors ! - emax = maxval(errors) / (atol + rtol * umax) + call update_errors(2, 4) - if (emax <= 1.0d+00 .or. nrej >= mrej) then +! calculate the tolerance and estimate the next time step due to the error +! + errtol = maxval(errors) + errs(1) = errtol + fc = product(errs(:)**betas(:)) + fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi) + dte = dt * fc + + if (fc > fac .or. nrej >= mrej) then test = .false. errs(3) = errs(2) errs(2) = errs(1) - errs(1) = emax - errtol = emax - - dte = dt * min(fcmx, max(fcmn, & - fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) - - fcmx = facmax else - errs(1) = emax - - dte = dt * min(fcmx, max(fcmn, & - fc * errs(1)**k1 * errs(2)**k2 * errs(3)**k3)) dt = dte - fcmx = fac - nrej = nrej + 1 ! rejection count in the current step nrejections = nrejections + 1 + +! since the solution was rejected, we have to revert the primitive variables +! to the previous state +! + pdata => list_data + do while (associated(pdata)) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + + call update_variables(tm, dtm) + end if niterations = niterations + 1 end do - if (allocated(lerr)) deallocate(lerr) - -!= final step: U(n+1) = U(1) +!= final step: U(n+1) = U(1) - update the accepted solution ! - tm = time + dt - dtm = dt / r - pdata => list_data do while (associated(pdata)) @@ -2285,27 +2261,13 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) - - if (ibp > 0) then - adecay(:) = exp(aglm(:) * cmax * dt) - - pdata => list_data - do while (associated(pdata)) - - pdata%u(ibp,:,:,:) = adecay(pdata%meta%level) * pdata%u(ibp,:,:,:) - - pdata => pdata%next - end do - end if - #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */ !------------------------------------------------------------------------------- ! - end subroutine evolve_ssprk3_m + end subroutine evolve_ssprk32m ! !=============================================================================== !