From a9cf78983376d1f683797459ee1e22704314b629 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 22:16:40 -0300 Subject: [PATCH] EVOLUTION: Rewrite SSPRK2(1)m method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 258 ++++++++++++++++++++---------------------- 1 file changed, 120 insertions(+), 138 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 43b5984..e806006 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -279,7 +279,30 @@ module evolution registers = 3 stages = max(2, min(9, stages)) cfl = (stages - 1) * cfl - write(name_int, "('Optimal 2nd order SSPRK(',i0,',2) with error control')") stages + write(name_int, "('2ⁿᵈ-order ',i0,'-step embedded " // & + "SSPRK(',i0,',2) method')") stages, stages + + allocate(c(stages), dl(stages), bt(stages), bh(stages)) + + do n = 1, stages + c(n) = n - 1.0d+00 + end do + c (:) = c(:) / (stages - 1) + bt(:) = 1.0d+00 / (stages - 1) + bt(stages) = 1.0d+00 / stages + bh(:) = 1.0d+00 / stages + bh(1) = (stages + 1.0d+00) / stages**2 + bh(stages) = (stages - 1.0d+00) / stages**2 + dl(1) = (stages - 1.0d+00) / stages + dl(2) = 1.0d+00 - dl(1) + + betas(:) = [ 5.8d-01, -2.1d-01, 1.0d-01 ] + + call get_parameter("beta(1)", betas(1)) + call get_parameter("beta(2)", betas(2)) + call get_parameter("beta(3)", betas(3)) + + betas(:) = - betas(:) / 2.0d+00 case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4") @@ -2327,7 +2350,11 @@ module evolution ! ! References: ! -! [1] Gottlieb, S. and Gottlieb, L.-A., J. +! [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. and Gottlieb, L.-A., J. ! "Strong Stability Preserving Properties of Runge-Kutta Time ! Discretization Methods for Linear Constant Coefficient Operators", ! Journal of Scientific Computing, @@ -2337,71 +2364,43 @@ module evolution ! subroutine evolve_ssprk2_m() - use blocks , only : block_data, list_data, get_dblocks - 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 sources , only : update_sources + use blocks , only : block_data, list_data + use boundaries , only : boundary_fluxes + use coordinates , only : nb, ne + use equations , only : errors, ibp, cmax + use iso_fortran_env, only : error_unit + use sources , only : update_sources implicit none type(block_data), pointer :: pdata - logical :: test - integer :: n, l, nrej, status - real(kind=8) :: tm, dtm, ds, umax, emax - real(kind=8) :: fc, fcmn, fcmx + logical :: test + integer :: nrej, i, status + real(kind=8) :: tm, dtm, dh, fc - logical , save :: first = .true. - real(kind=8), save :: ft, fl, fr, gt, gl + character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk2_m()' - real(kind=8), dimension(:,:,:,:,:), allocatable :: lerr - - real(kind=8), parameter :: k1 = -2.9d-01, k2 = 1.05d-01, k3 = -5.0d-02 -! !------------------------------------------------------------------------------- ! #ifdef PROFILE call start_timer(imu) #endif /* PROFILE */ - if (first) then - - ft = 1.0d+00 / (stages - 1) - fl = 1.0d+00 / stages - fr = 1.0d+00 - fl - - gt = fl - ft - gl = fl - - 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 2nd order solution, U(1) +! pdata%uu(:,:,:,:,3) - the new 1st order solution, U(2) +! do while(test) - lerr(:,:,:,:,:) = 0.0d+00 - - ds = ft * dt - -!= 1st step: U(1) = U(n), U(2) = U(n) +!= preparatory step: U(1) = U(n), U(2) = 0, U(3) = U(n) ! pdata => list_data do while (associated(pdata)) @@ -2414,12 +2413,12 @@ module evolution pdata => pdata%next end do -!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 +!= 1st step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 ! - do n = 1, stages - 1 + do i = 1, stages - 1 - tm = time + n * ds - dtm = ds + dtm = c(i) * dt + tm = time + dtm pdata => list_data do while (associated(pdata)) @@ -2432,30 +2431,27 @@ 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,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gt * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) & + + bt(i) * dt * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) & + + bh(i) * dt * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 - end do ! n = 1, stages - 1 + end do ! i = 1, stages - 1 -!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] +!= 2nd step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] ! - tm = time + dt - dtm = fl * dt + i = stages + dtm = c(i) * dt + tm = time + dtm pdata => list_data do while (associated(pdata)) @@ -2468,93 +2464,93 @@ module evolution call boundary_fluxes() - l = 1 pdata => list_data do while (associated(pdata)) - pdata%uu(:,:,:,:,2) = fr * pdata%uu(:,:,:,:,2) & - + fl * (pdata%uu(:,:,:,:,3) & - + dt * pdata%du(:,:,:,:)) - -#if NDIMS == 3 - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - lerr(l,:,:,:,:) = lerr(l,:,:,:,:) + gl * pdata%du(:,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - l = l + 1 + pdata%uu(:,:,:,:,2) = dl(1) * pdata%uu(:,:,:,:,2) & + + dl(2) * pdata%uu(:,:,:,:,1) & + + bt(i) * dt * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,3) & + + bh(i) * dt * pdata%du(:,:,:,:) pdata => pdata%next end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 -! find umax +!= 3rd step: decay the magnetic divergence potential of both solutions ! - 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 => pdata%next - end do + if (ibp > 0) then + adecay(:) = exp(aglm(:) * cmax * dt) -! get the maximum error of each variable + pdata => list_data + do while (associated(pdata)) + + pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,2) + pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) & + * pdata%uu(ibp,:,:,:,3) + + pdata => pdata%next + end do + end if + +! update the vector of errors ! - do l = 1, nf - errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) - end do + call update_errors(2, 3) -#ifdef MPI - call reduce_maximum(umax) - call reduce_maximum(errors) -#endif /* MPI */ - -! calculate tolerance and time step +! calculate the tolerance and estimate the next time step due to the error ! - emax = maxval(errors) / (atol + rtol * umax) + 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 (emax <= 1.0d+00 .or. nrej >= mrej) then + 100 continue + + if ((fc > fac .or. nrej >= mrej) .and. status == 0) 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 + if (status == 0) then + dt = dte + nrej = nrej + 1 ! rejection count in the current step + else + dt = 2.5d-01 * dt + end if 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, status) + if (status /= 0) then + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "Possible bug: the revert to the previous state " //& + "found it unphysical." + end if + 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 = ft * dt - pdata => list_data do while (associated(pdata)) @@ -2565,20 +2561,6 @@ module evolution pdata => pdata%next end do - 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 - - call update_variables(tm, dtm, status) - #ifdef PROFILE call stop_timer(imu) #endif /* PROFILE */