From 86cf55e836224cc1b49a4960ee5e16ad8fb1bf39 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 12:32:49 -0300 Subject: [PATCH 1/8] EVOLUTION: Add status flag to update_variables(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 141 +++++++++++++++++++----------------------- 1 file changed, 64 insertions(+), 77 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 93e2f8c..8136418 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1098,7 +1098,8 @@ module evolution ! update primitive variables ! - call update_variables(time + dt, 0.0d+00) + call update_variables(time + dt, 0.0d+00, status) + if (status /= 0) go to 100 ! set all meta blocks to be updated ! @@ -1163,7 +1164,7 @@ module evolution type(block_data), pointer :: pdata - integer :: mlev + integer :: mlev, status real(kind=8) :: dx_min, fnorm, h0, h1 real(kind=8), dimension(3) :: d @@ -1264,7 +1265,7 @@ module evolution pdata => pdata%next end do - call update_variables(time + h0, h0) + call update_variables(time + h0, h0, status) pdata => list_data do while (associated(pdata)) @@ -1438,6 +1439,7 @@ module evolution type(block_data), pointer :: pdata + integer :: status real(kind=8) :: tm, dtm ! !------------------------------------------------------------------------------- @@ -1482,7 +1484,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -1519,6 +1521,7 @@ module evolution type(block_data), pointer :: pdata + integer :: status real(kind=8) :: tm, dtm ! !------------------------------------------------------------------------------- @@ -1553,7 +1556,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 2nd step: U(n+1) = 1/2 U(n) + 1/2 U(1) + 1/2 dt * F[U(1)] ! @@ -1595,7 +1598,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -1640,7 +1643,7 @@ module evolution type(block_data), pointer :: pdata logical :: test - integer :: n, l, nrej + integer :: n, l, nrej, status real(kind=8) :: tm, dtm, ds, umax, emax real(kind=8) :: fc, fcmn, fcmx @@ -1737,7 +1740,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! n = 1, stages - 1 @@ -1775,7 +1778,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) ! find umax ! @@ -1866,7 +1869,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -1904,6 +1907,7 @@ module evolution type(block_data), pointer :: pdata + integer :: status real(kind=8) :: ds real(kind=8) :: tm, dtm @@ -1943,7 +1947,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)] ! @@ -1971,7 +1975,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)] ! @@ -2013,7 +2017,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -2052,6 +2056,7 @@ module evolution type(block_data), pointer :: pdata + integer :: status real(kind=8) :: ds real(kind=8) :: tm, dtm ! @@ -2088,7 +2093,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 2nd step: U(2) = U(2) + 1/2 dt F[U(1)] ! @@ -2113,7 +2118,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 3rd step: U(3) = 2/3 U(n) + 1/3 (U(2) + 1/2 dt F[U(2)]) ! @@ -2140,7 +2145,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != the final step: U(n+1) = U(3) + 1/2 dt F[U(3)] ! @@ -2179,7 +2184,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -2228,7 +2233,7 @@ module evolution type(block_data), pointer :: pdata logical :: test - integer :: nrej, i + integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00, & @@ -2295,7 +2300,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! i = 1, 3 @@ -2312,7 +2317,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate ! @@ -2338,7 +2343,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation ! @@ -2403,7 +2408,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end if @@ -2460,6 +2465,7 @@ module evolution type(block_data), pointer :: pdata + integer :: status real(kind=8) :: ds real(kind=8) :: tm, dtm @@ -2507,7 +2513,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 2nd step: U(2) = U(1) + b1 dt F[U(1)] ! @@ -2533,7 +2539,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] ! @@ -2561,7 +2567,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] ! @@ -2591,7 +2597,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] ! @@ -2631,7 +2637,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -2679,7 +2685,7 @@ module evolution real(kind=8), save :: f1, f2, f3, f4 logical :: test - integer :: nrej, i + integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc ! !------------------------------------------------------------------------------- @@ -2766,7 +2772,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! i = 1, i1 @@ -2806,7 +2812,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! i = i2, i3 @@ -2827,7 +2833,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) != 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2+1...stages ! @@ -2855,7 +2861,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! i = i4, stages @@ -2921,7 +2927,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end if @@ -2977,7 +2983,7 @@ module evolution type(block_data), pointer :: pdata - integer :: n + integer :: n, status real(kind=8) :: tm, dtm real(kind=8), dimension(9) :: ds @@ -3030,7 +3036,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! n = 1, 5 @@ -3075,7 +3081,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! n = 6, 9 @@ -3117,7 +3123,7 @@ module evolution end do end if - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) #ifdef PROFILE call stop_timer(imu) @@ -3160,7 +3166,7 @@ module evolution type(block_data), pointer :: pdata logical :: test - integer :: i, l, nrej + integer :: i, l, nrej, status real(kind=8) :: tm, dtm, fc ! !------------------------------------------------------------------------------- @@ -3224,7 +3230,7 @@ module evolution end do dtm = c(2) * dt - call update_variables(time, dtm) + call update_variables(time, dtm, status) != remaining stages ! @@ -3260,7 +3266,7 @@ module evolution end do dtm = c(i+1) * dt - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) end do ! i = 2, stages @@ -3345,7 +3351,7 @@ module evolution pdata => pdata%next end do - call update_variables(tm, dtm) + call update_variables(tm, dtm, status) if (fsal) fsal_update = .false. end if @@ -3367,7 +3373,7 @@ module evolution end do tm = time + dt - call update_variables(tm, dt) + call update_variables(tm, dt, status) #ifdef PROFILE call stop_timer(imu) @@ -3585,49 +3591,38 @@ module evolution ! ! Arguments: ! -! tm - time at the moment of update; -! dtm - time step since the last update; +! tm - time at the moment of update; +! dtm - time step since the last update; +! status - status flag indicating if the update was successful; ! !=============================================================================== ! - subroutine update_variables(tm, dtm) + subroutine update_variables(tm, dtm, status) -! include external procedures -! - use blocks , only : set_neighbors_update - use boundaries , only : boundary_variables - use equations , only : update_primitive_variables - use equations , only : fix_unphysical_cells, correct_unphysical_states - use shapes , only : update_shapes + use blocks , only : block_meta, list_meta + use blocks , only : block_data, list_data + use blocks , only : set_neighbors_update + use boundaries, only : boundary_variables + use equations , only : update_primitive_variables + use equations , only : fix_unphysical_cells, correct_unphysical_states + use shapes , only : update_shapes -! include external variables -! - use blocks , only : block_meta, list_meta - use blocks , only : block_data, list_data - -! local variables are not implicit by default -! implicit none -! subroutine arguments -! - real(kind=8), intent(in) :: tm, dtm + real(kind=8), intent(in) :: tm, dtm + integer , intent(out) :: status -! local pointers -! type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable update -! call start_timer(imv) #endif /* PROFILE */ -! update primitive variables in the changed blocks -! + status = 0 + pdata => list_data do while (associated(pdata)) pmeta => pdata%meta @@ -3637,12 +3632,8 @@ module evolution pdata => pdata%next end do -! update boundaries -! call boundary_variables(tm, dtm) -! apply shapes in blocks which need it -! pdata => list_data do while (associated(pdata)) pmeta => pdata%meta @@ -3652,8 +3643,6 @@ module evolution pdata => pdata%next end do -! correct unphysical states if detected -! if (fix_unphysical_cells) then ! if an unphysical cell appeared in a block while updating its primitive @@ -3683,8 +3672,6 @@ module evolution end if #ifdef PROFILE -! stop accounting time for variable update -! call stop_timer(imv) #endif /* PROFILE */ From 17f867cb082271db7959ba983c6f78f4ab826d1f Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 12:51:23 -0300 Subject: [PATCH 2/8] EQUATIONS: Add status flag to update_primitive_variables(). Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 5 ++++- sources/evolution.F90 | 7 ++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index e8f1723..da50da6 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -1412,7 +1412,7 @@ module equations ! !=============================================================================== ! - subroutine update_primitive_variables(uu, qq) + subroutine update_primitive_variables(uu, qq, status) ! include external procedures and variables ! @@ -1427,6 +1427,7 @@ module equations ! real(kind=8), dimension(:,:,:,:), intent(inout) :: uu real(kind=8), dimension(:,:,:,:), intent(inout) :: qq + integer , intent(out) :: status ! temporary variables ! @@ -1434,6 +1435,8 @@ module equations ! !------------------------------------------------------------------------------- ! + status = 0 + ! update primitive variables ! #if NDIMS == 3 diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 8136418..878f117 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -3627,7 +3627,10 @@ module evolution do while (associated(pdata)) pmeta => pdata%meta - if (pmeta%update) call update_primitive_variables(pdata%u, pdata%q) + if (pmeta%update) then + call update_primitive_variables(pdata%u, pdata%q, status) + if (status /= 0) go to 100 + end if pdata => pdata%next end do @@ -3671,6 +3674,8 @@ module evolution end do end if + 100 continue + #ifdef PROFILE call stop_timer(imv) #endif /* PROFILE */ From ea49382c499c023386372b37a78802e13ce790e7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 13:21:32 -0300 Subject: [PATCH 3/8] EQUATIONS: Add status flag to all cons2prim_*() subroutines. We also check if the positive variables are indeed positive. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 283 ++++++++++++++++++------------------------ 1 file changed, 123 insertions(+), 160 deletions(-) diff --git a/sources/equations.F90 b/sources/equations.F90 index da50da6..9b87d91 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -75,9 +75,10 @@ module equations real(kind=8), dimension(:,:), intent(out) :: u logical , optional , intent(in) :: s end subroutine - subroutine cons2prim_iface(u, q) + subroutine cons2prim_iface(u, q, s) real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s end subroutine subroutine fluxspeed_iface(q, u, f, c) real(kind=8), dimension(:,:) , intent(in) :: q, u @@ -1414,47 +1415,35 @@ module equations ! subroutine update_primitive_variables(uu, qq, status) -! include external procedures and variables -! use coordinates, only : nn => bcells use coordinates, only : nb, ne, nbl, neu -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:,:,:), intent(inout) :: uu real(kind=8), dimension(:,:,:,:), intent(inout) :: qq integer , intent(out) :: status -! temporary variables -! integer :: i, j, k = 1 -! + !------------------------------------------------------------------------------- ! status = 0 -! update primitive variables -! #if NDIMS == 3 do k = nb, ne #endif /* NDIMS == 3 */ do j = nb, ne -! convert conserved variables to primitive ones -! - call cons2prim(uu(1:nv,nb:ne,j,k), qq(1:nv,nb:ne,j,k)) + call cons2prim(uu(1:nv,nb:ne,j,k), qq(1:nv,nb:ne,j,k), status) + + if (status /= 0) go to 100 end do ! j = nb, ne #if NDIMS == 3 end do ! k = nb, ne #endif /* NDIMS == 3 */ -! fill out the borders -! #if NDIMS == 3 do i = nbl, 1, -1 qq(1:nv, i,nb:ne,nb:ne) = qq(1:nv,nb,nb:ne,nb:ne) @@ -1489,6 +1478,8 @@ module equations end do #endif /* NDIMS == 3 */ + 100 continue + !------------------------------------------------------------------------------- ! end subroutine update_primitive_variables @@ -1793,54 +1784,51 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_hd_iso(u, q) + subroutine cons2prim_hd_iso(u, q, s) -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! integer :: i, p -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ -! iterate over all positions -! + s = 0 + do i = 1, size(u,2) - q(idn,i) = u(idn,i) - q(ivx,i) = u(imx,i) / u(idn,i) - q(ivy,i) = u(imy,i) / u(idn,i) - q(ivz,i) = u(imz,i) / u(idn,i) + if (u(idn,i) > 0.0d+00) then + q(idn,i) = u(idn,i) + q(ivx,i) = u(imx,i) / u(idn,i) + q(ivy,i) = u(imy,i) / u(idn,i) + q(ivz,i) = u(imz,i) / u(idn,i) + else + s = 1 + go to 100 + end if end do -! update primitive passive scalars -! if (ns > 0) then do p = isl, isu q(p,:) = u(p,:) / u(idn,:) end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ @@ -2219,59 +2207,59 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_hd_adi(u, q) + subroutine cons2prim_hd_adi(u, q, s) -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! integer :: i, p real(kind=8) :: ek, ei -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ -! iterate over all positions -! + s = 0 + do i = 1, size(u,2) - q(idn,i) = u(idn,i) - q(ivx,i) = u(imx,i) / u(idn,i) - q(ivy,i) = u(imy,i) / u(idn,i) - q(ivz,i) = u(imz,i) / u(idn,i) - ek = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i) & - + u(imz,i) * q(ivz,i)) - ei = u(ien,i) - ek - q(ipr,i) = gammam1 * ei - + if (u(idn,i) > 0.0d+00) then + q(idn,i) = u(idn,i) + q(ivx,i) = u(imx,i) / u(idn,i) + q(ivy,i) = u(imy,i) / u(idn,i) + q(ivz,i) = u(imz,i) / u(idn,i) + ek = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i)) + ei = u(ien,i) - ek + if (ei > 0.0d+00) then + q(ipr,i) = gammam1 * ei + else + s = 1 + go to 100 + end if + else + s = 1 + go to 100 + end if end do -! update primitive passive scalars -! if (ns > 0) then do p = isl, isu q(p,:) = u(p,:) / u(idn,:) end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ @@ -2695,58 +2683,55 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_mhd_iso(u, q) + subroutine cons2prim_mhd_iso(u, q, s) -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! integer :: i, p -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ -! iterate over all positions -! + s = 0 + do i = 1, size(u,2) - q(idn,i) = u(idn,i) - q(ivx,i) = u(imx,i) / u(idn,i) - q(ivy,i) = u(imy,i) / u(idn,i) - q(ivz,i) = u(imz,i) / u(idn,i) - q(ibx,i) = u(ibx,i) - q(iby,i) = u(iby,i) - q(ibz,i) = u(ibz,i) - q(ibp,i) = u(ibp,i) + if (u(idn,i) > 0.0d+00) then + q(idn,i) = u(idn,i) + q(ivx,i) = u(imx,i) / u(idn,i) + q(ivy,i) = u(imy,i) / u(idn,i) + q(ivz,i) = u(imz,i) / u(idn,i) + q(ibx,i) = u(ibx,i) + q(iby,i) = u(iby,i) + q(ibz,i) = u(ibz,i) + q(ibp,i) = u(ibp,i) + else + s = 1 + go to 100 + end if end do -! update primitive passive scalars -! if (ns > 0) then do p = isl, isu q(p,:) = u(p,:) / u(idn,:) end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ @@ -3311,64 +3296,65 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_mhd_adi(u, q) + subroutine cons2prim_mhd_adi(u, q, s) -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! integer :: i, p real(kind=8) :: ei, ek, em -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ -! iterate over all positions -! + s = 0 + do i = 1, size(u,2) - q(idn,i) = u(idn,i) - q(ivx,i) = u(imx,i) / u(idn,i) - q(ivy,i) = u(imy,i) / u(idn,i) - q(ivz,i) = u(imz,i) / u(idn,i) - q(ibx,i) = u(ibx,i) - q(iby,i) = u(iby,i) - q(ibz,i) = u(ibz,i) - q(ibp,i) = u(ibp,i) - ek = 0.5d+00 * (u(imx,i) * q(ivx,i) + u(imy,i) * q(ivy,i) & - + u(imz,i) * q(ivz,i)) - em = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i)) - ei = u(ien,i) - (ek + em) - q(ipr,i) = gammam1 * ei + if (u(idn,i) > 0.0d+00) then + q(idn,i) = u(idn,i) + q(ivx,i) = u(imx,i) / u(idn,i) + q(ivy,i) = u(imy,i) / u(idn,i) + q(ivz,i) = u(imz,i) / u(idn,i) + q(ibx,i) = u(ibx,i) + q(iby,i) = u(iby,i) + q(ibz,i) = u(ibz,i) + q(ibp,i) = u(ibp,i) + ek = 0.5d+00 * sum(u(imx:imz,i) * q(ivx:ivz,i)) + em = 0.5d+00 * sum(q(ibx:ibz,i) * q(ibx:ibz,i)) + ei = u(ien,i) - (ek + em) + if (ei > 0.0d+00) then + q(ipr,i) = gammam1 * ei + else + s = 1 + go to 100 + end if + else + s = 1 + go to 100 + end if end do -! update primitive passive scalars -! if (ns > 0) then do p = isl, isu q(p,:) = u(p,:) / u(idn,:) end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ @@ -3990,45 +3976,35 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_srhd_adi(u, q) + subroutine cons2prim_srhd_adi(u, q, s) -! include external procedures -! use iso_fortran_env, only : error_unit -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! logical :: info integer :: i, p real(kind=8) :: mm, bb, mb, en, dn real(kind=8) :: w , vv, vm, vs -! local parameters -! character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srhd_adi()' -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ -! iterate over all positions -! + s = 0 + do i = 1, size(u,2) ! prepare variables which do not change during the Newton-Ralphson iterations @@ -4065,9 +4041,8 @@ module equations write(error_unit,"(a,5(1x,1es24.16e3))") "U = ", u(:,i) write(error_unit,"(a,3(1x,1es24.16e3))") "D, |m|², E = ", dn, mm, en -! set pressure to zero so we can hopefully fix it later -! - q(ipr,i) = 0.0d+00 + s = 1 + go to 100 end if @@ -4081,9 +4056,9 @@ module equations end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ @@ -5129,43 +5104,35 @@ module equations ! ! u - the input array of conservative variables; ! q - the output array of primitive variables; +! s - the status flag; ! !=============================================================================== ! - subroutine cons2prim_srmhd_adi(u, q) + subroutine cons2prim_srmhd_adi(u, q, s) -! include external procedures -! use iso_fortran_env, only : error_unit -! local variables are not implicit by default -! implicit none -! input/output arguments -! real(kind=8), dimension(:,:), intent(in) :: u real(kind=8), dimension(:,:), intent(out) :: q + integer , intent(out) :: s -! local variables -! logical :: info integer :: i, p real(kind=8) :: mm, mb, bb, en, dn real(kind=8) :: w, wt, vv, vm, vs, fc -! local parameters -! character(len=*), parameter :: loc = 'EQUATIONS::cons2prim_srmhd_adi()' -! + !------------------------------------------------------------------------------- ! #ifdef PROFILE -! start accounting time for variable conversion -! call start_timer(imc) #endif /* PROFILE */ + s = 0 + ! iterate over all positions ! do i = 1, size(u,2) @@ -5219,9 +5186,8 @@ module equations write(error_unit,"(a,6(1x,1es24.16e3))") "D, |m|², m.B, |B|², E, W = "& , dn, mm, mb, bb, en, w -! set pressure to zero so we can hopefully fix it later -! - q(ipr,i) = 0.0d+00 + s = 1 + go to 100 end if ! p <= 0 @@ -5234,25 +5200,22 @@ module equations write(error_unit,"(a,5(1x,1es24.16e3))") "D, |m|², m.B, |B|², E = " & , dn, mm, mb, bb, en -! set pressure to zero so we can hopefully fix it later -! - q(ipr,i) = 0.0d+00 + s = 1 + go to 100 end if ! unphysical state end do -! update primitive passive scalars -! if (ns > 0) then do p = isl, isu q(p,:) = u(p,:) / u(idn,:) end do end if + 100 continue + #ifdef PROFILE -! stop accounting time for variable conversion -! call stop_timer(imc) #endif /* PROFILE */ From 89bc70667343307bd91a6b0849208d1e9e168ee5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 15:37:47 -0300 Subject: [PATCH 4/8] EVOLUTION: Treat any unphysical state properly in 3S*+ method. If any unphysical cell is detected, the integration is immediately reverted and resumed with a time step reduced by a factor or four. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 38 +++++++++++++++++++++++++++++--------- 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 878f117..758def9 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -3156,10 +3156,11 @@ module evolution ! subroutine evolve_3sstarplus() - use blocks , only : block_data, list_data - use boundaries , only : boundary_fluxes - use equations , only : errors, ibp, cmax - use sources , only : update_sources + use blocks , only : block_data, list_data + use boundaries , only : boundary_fluxes + use equations , only : errors, ibp, cmax + use iso_fortran_env, only : error_unit + use sources , only : update_sources implicit none @@ -3168,7 +3169,9 @@ module evolution logical :: test integer :: i, l, nrej, status real(kind=8) :: tm, dtm, fc -! + + character(len=*), parameter :: loc = 'EVOLUTION:evolve_3sstarplus()' + !------------------------------------------------------------------------------- ! #ifdef PROFILE @@ -3231,6 +3234,7 @@ module evolution dtm = c(2) * dt call update_variables(time, dtm, status) + if (status /= 0) go to 100 != remaining stages ! @@ -3267,6 +3271,7 @@ module evolution dtm = c(i+1) * dt call update_variables(tm, dtm, status) + if (status /= 0) go to 100 end do ! i = 2, stages @@ -3327,7 +3332,9 @@ module evolution fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi) dte = dt * fc - if (fc > fac .or. nrej >= mrej) then + 100 continue + + if ((fc > fac .or. nrej >= mrej) .and. status == 0) then test = .false. errs(3) = errs(2) @@ -3335,9 +3342,12 @@ module evolution if (fsal) fsal_update = .true. else - dt = dte - - 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 @@ -3352,6 +3362,11 @@ module evolution 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 if (fsal) fsal_update = .false. end if @@ -3374,6 +3389,11 @@ module evolution tm = time + dt call update_variables(tm, dt, status) + if (status /= 0) then + write(error_unit,"('[', a, ']: ', a)") trim(loc), & + "Possible bug: the new state has been accepted, " // & + "nevertheless the variable update found it unphysical." + end if #ifdef PROFILE call stop_timer(imu) From 380e4c8d7e532d0335b64565fcabb7255fea22ef Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 15:59:34 -0300 Subject: [PATCH 5/8] EVOLUTION: Synchronize status among processes in update_variables(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 758def9..6ddcec3 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -3625,6 +3625,9 @@ module evolution use boundaries, only : boundary_variables use equations , only : update_primitive_variables use equations , only : fix_unphysical_cells, correct_unphysical_states +#ifdef MPI + use mpitools , only : reduce_maximum +#endif /* MPI */ use shapes , only : update_shapes implicit none @@ -3655,6 +3658,12 @@ module evolution pdata => pdata%next end do + 100 continue +#ifdef MPI + call reduce_maximum(status) +#endif /* MPI */ + if (status /= 0) go to 200 + call boundary_variables(tm, dtm) pdata => list_data @@ -3694,7 +3703,7 @@ module evolution end do end if - 100 continue + 200 continue #ifdef PROFILE call stop_timer(imv) From a52da2830d675ca051d69a3f587b578f604ed3b4 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 16:05:22 -0300 Subject: [PATCH 6/8] =?UTF-8?q?EVOLUTION:=20Treat=20any=20unphysical=20sta?= =?UTF-8?q?te=20properly=20in=20SSPRK3(2)n=C2=B2=20method.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 6ddcec3..11516dc 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -2671,10 +2671,11 @@ module evolution ! subroutine evolve_ssprk32m() - use blocks , only : block_data, list_data - use boundaries , only : boundary_fluxes - use equations , only : errors, ibp, cmax - use sources , only : update_sources + use blocks , only : block_data, list_data + use boundaries , only : boundary_fluxes + use equations , only : errors, ibp, cmax + use iso_fortran_env, only : error_unit + use sources , only : update_sources implicit none @@ -2687,7 +2688,9 @@ module evolution logical :: test integer :: nrej, i, status real(kind=8) :: tm, dtm, dh, fc -! + + character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssprk32m()' + !------------------------------------------------------------------------------- ! #ifdef PROFILE @@ -2773,6 +2776,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 end do ! i = 1, i1 @@ -2813,6 +2817,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 end do ! i = i2, i3 @@ -2834,6 +2839,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 != 6th step: U(1) = U(1) + dt/r F[U(1)], for i = n*(n+1)/2+1...stages ! @@ -2862,6 +2868,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 end do ! i = i4, stages @@ -2904,16 +2911,20 @@ module evolution fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi) dte = dt * fc - if (fc > fac .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) - else - dt = dte - - 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 @@ -2928,6 +2939,11 @@ module evolution 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 From 6dcd4002190e839ca5e0a370a6ebd4b79931cee4 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 16:09:32 -0300 Subject: [PATCH 7/8] EVOLUTION: Treat any unphysical state properly in SSPRK3(2)4 method. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 36 ++++++++++++++++++++++++++---------- 1 file changed, 26 insertions(+), 10 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 11516dc..f506265 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -2222,11 +2222,12 @@ module evolution ! subroutine evolve_ssp324() - use blocks , only : block_data, list_data - use boundaries , only : boundary_fluxes - use coordinates, only : nb, ne - use equations , only : errors, ibp, cmax - 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 @@ -2238,7 +2239,9 @@ module evolution real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00, & twothird = 2.0d+00 / 3.0d+00 -! + + character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()' + !------------------------------------------------------------------------------- ! #ifdef PROFILE @@ -2301,6 +2304,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 end do ! i = 1, 3 @@ -2318,6 +2322,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 != 5th step: U(1) = U(1) + ½ dt F[U(1)] <- 3ʳᵈ-order candidate ! @@ -2344,6 +2349,7 @@ module evolution end do call update_variables(tm, dtm, status) + if (status /= 0) go to 100 != 6th step: U(2) = ½ (U(1) + U(2)) <- 2ⁿᵈ-order approximation ! @@ -2385,16 +2391,21 @@ module evolution fc = 1.0d+00 + chi * atan((fc - 1.0d+00) / chi) dte = dt * fc - if (fc > fac .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) else - dt = dte - - 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 @@ -2409,6 +2420,11 @@ module evolution 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 From a87bcce9078e7b4b6c7054aeaf04b0d5f4507f18 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 9 Nov 2021 16:16:26 -0300 Subject: [PATCH 8/8] EVOLUTION: Reorganize methods. Put first non-embedded and then embedded methods. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 1360 +++++++++++++++++++++-------------------- 1 file changed, 682 insertions(+), 678 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index f506265..43b5984 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -240,16 +240,6 @@ module evolution registers = 2 name_int = "2nd order Runge-Kutta" - case ("ssprk(m,2)", "SSPRK(m,2)") - - error_control = .true. - evolve => evolve_ssprk2_m - order = 2 - 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 - case ("rk3", "RK3") evolve => evolve_rk3 @@ -273,6 +263,24 @@ module evolution cfl = 2.65062919143939d+00 * cfl name_int = "3rd order SSPRK(5,3)" + case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") + + evolve => evolve_ssprk4_10 + order = 4 + registers = 2 + cfl = 6.0d+00 * cfl + name_int = "Optimal 4th order SSPRK(10,4)" + + case ("ssprk(m,2)", "SSPRK(m,2)") + + error_control = .true. + evolve => evolve_ssprk2_m + order = 2 + 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 + case ("SSPRK3(2)4", "SSP3(2)4", "ssprk3(2)4", "ssp3(2)4") error_control = .true. @@ -324,14 +332,6 @@ module evolution " significantly used less memory." end if - case ("rk4.10", "RK4.10", "ssprk(10,4)", "SSPRK(10,4)") - - evolve => evolve_ssprk4_10 - order = 4 - registers = 2 - cfl = 6.0d+00 * cfl - name_int = "Optimal 4th order SSPRK(10,4)" - case ("RK3(2)5", "rk3(2)5") ! References: @@ -1610,277 +1610,6 @@ module evolution ! !=============================================================================== ! -! subroutine EVOLVE_SSPRK2_M: -! -------------------------- -! -! Subroutine advances the solution by one time step using the 2nd order -! m-stage Strong Stability Preserving Runge-Kutta time integration method. -! Up to 9 stages are allowed, due to stability problems with more stages. -! -! References: -! -! [1] 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, -! 2003, vol. 18, no. 1, pp. 83-109 -! -!=============================================================================== -! - 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 - - 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 , save :: first = .true. - real(kind=8), save :: ft, fl, fr, gt, gl - - 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 - - do while(test) - - lerr(:,:,:,:,:) = 0.0d+00 - - ds = ft * dt - -!= 1st step: U(1) = U(n), U(2) = U(n) -! - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) - pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) - - pdata%u => pdata%uu(:,:,:,:,2) - - pdata => pdata%next - end do - -!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 -! - do n = 1, stages - 1 - - tm = time + n * ds - 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) = 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 => pdata%next - end do - - call update_variables(tm, dtm, status) - - end do ! n = 1, stages - 1 - -!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] -! - tm = time + dt - dtm = fl * dt - - 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) = 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 => pdata%next - end do - - call update_variables(tm, dtm, status) - -! find umax -! - 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 - -! get the maximum error of each variable -! - do l = 1, nf - errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) - end do - -#ifdef MPI - call reduce_maximum(umax) - call reduce_maximum(errors) -#endif /* MPI */ - -! calculate tolerance and time step -! - emax = maxval(errors) / (atol + rtol * umax) - - if (emax <= 1.0d+00 .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 - end if - - niterations = niterations + 1 - - end do - - if (allocated(lerr)) deallocate(lerr) - -!= final step: U(n+1) = U(1) -! - tm = time + dt - dtm = ft * dt - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) - - pdata%u => pdata%uu(:,:,:,:,1) - - 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 */ - -!------------------------------------------------------------------------------- -! - end subroutine evolve_ssprk2_m -! -!=============================================================================== -! ! subroutine EVOLVE_RK3: ! --------------------- ! @@ -2196,6 +1925,670 @@ module evolution ! !=============================================================================== ! +! subroutine EVOLVE_SSPRK35: +! ------------------------- +! +! Subroutine advances the solution by one time step using the 3rd order +! 5-stage Strong Stability Preserving Runge-Kutta time integration method. +! +! References: +! +! [1] Ruuth, S. J., +! "Global Optimization of Explicit Strong-Stability-Preserving +! Runge-Kutta methods", +! Mathematics of Computation, +! 2006, vol. 75, no. 253, pp. 183-207 +! +!=============================================================================== +! + subroutine evolve_ssprk35() + + use blocks , only : block_data, list_data + use boundaries, only : boundary_fluxes + use equations , only : ibp, cmax + use sources , only : update_sources + + implicit none + + type(block_data), pointer :: pdata + + integer :: status + real(kind=8) :: ds + real(kind=8) :: tm, dtm + + real(kind=8), parameter :: b1 = 3.77268915331368d-01 + real(kind=8), parameter :: b3 = 2.42995220537396d-01 + real(kind=8), parameter :: b4 = 2.38458932846290d-01 + real(kind=8), parameter :: b5 = 2.87632146308408d-01 + real(kind=8), parameter :: a31 = 3.55909775063327d-01 + real(kind=8), parameter :: a33 = 6.44090224936673d-01 + real(kind=8), parameter :: a41 = 3.67933791638137d-01 + real(kind=8), parameter :: a44 = 6.32066208361863d-01 + real(kind=8), parameter :: a53 = 2.37593836598569d-01 + real(kind=8), parameter :: a55 = 7.62406163401431d-01 +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE + call start_timer(imu) +#endif /* PROFILE */ + +!= 1st step: U(1) = U(n) + b1 dt F[U(n)] +! + ds = b1 * dt + tm = time + ds + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + +!= 2nd step: U(2) = U(1) + b1 dt F[U(1)] +! + tm = time + 2.0d+00 * ds + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + +!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] +! + ds = b3 * dt + tm = time + (2.0d+00 * a33 * b1 + b3) * 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) & + + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + +!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] +! + ds = b4 * dt + tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = a41 * pdata%uu(:,:,:,:,1) & + + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + + pdata%u => pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + +!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] +! + ds = b5 * dt + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2) & + + a55 * pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) + + 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 */ + +!------------------------------------------------------------------------------- +! + end subroutine evolve_ssprk35 +! +!=============================================================================== +! +! subroutine EVOLVE_SSPRK4_10: +! --------------------------- +! +! Subroutine advances the solution by one time step using the 4rd order +! 10-stage Strong Stability Preserving Runge-Kutta time integration method. +! +! References: +! +! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., +! "Strong stability preserving Runge-Kutta and multistep +! time discretizations", +! World Scientific Publishing, 2011 +! +!=============================================================================== +! + subroutine evolve_ssprk4_10() + + use blocks , only : block_data, list_data + use boundaries, only : boundary_fluxes + use equations , only : ibp, cmax + use sources , only : update_sources + + implicit none + + type(block_data), pointer :: pdata + + integer :: n, status + real(kind=8) :: tm, dtm + + real(kind=8), dimension(9) :: ds + + real(kind=8), dimension(9), parameter :: c = & + (/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00, & + 3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00 +! +!------------------------------------------------------------------------------- +! +#ifdef PROFILE + call start_timer(imu) +#endif /* PROFILE */ + + ds(:) = c(:) * dt + +!= 1st step: U(2) = U(1) +! + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + +!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 +! + do n = 1, 5 + + tm = time + ds(n) + dtm = ds(1) + + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + + end do ! n = 1, 5 + +!= 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) +! + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2) & + + 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01 + pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2) & + - 5.0d+00 * pdata%uu(:,:,:,:,1) + + pdata => pdata%next + end do + +!= 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9 +! +! integrate the intermediate steps +! + do n = 6, 9 + + tm = time + ds(n) + dtm = ds(1) + + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) + + pdata => pdata%next + end do + + call update_variables(tm, dtm, status) + + end do ! n = 6, 9 + +!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] +! + tm = time + dt + dtm = dt / 1.0d+01 + + 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() + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) & + + 6.0d-01 * pdata%uu(:,:,:,:,1) & + + dtm * pdata%du(:,:,:,:) + + 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 */ + +!------------------------------------------------------------------------------- +! + end subroutine evolve_ssprk4_10 +! +!=============================================================================== +! +! EMBEDDED METHODS +! +!=============================================================================== +! +! subroutine EVOLVE_SSPRK2_M: +! -------------------------- +! +! Subroutine advances the solution by one time step using the 2nd order +! m-stage Strong Stability Preserving Runge-Kutta time integration method. +! Up to 9 stages are allowed, due to stability problems with more stages. +! +! References: +! +! [1] 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, +! 2003, vol. 18, no. 1, pp. 83-109 +! +!=============================================================================== +! + 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 + + 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 , save :: first = .true. + real(kind=8), save :: ft, fl, fr, gt, gl + + 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 + + do while(test) + + lerr(:,:,:,:,:) = 0.0d+00 + + ds = ft * dt + +!= 1st step: U(1) = U(n), U(2) = U(n) +! + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + pdata%uu(:,:,:,:,3) = pdata%uu(:,:,:,:,1) + + pdata%u => pdata%uu(:,:,:,:,2) + + pdata => pdata%next + end do + +!= 2nd step: U(1) = U(1) + dt/(m-1) F[U(1)], for i = 1, ..., m-1 +! + do n = 1, stages - 1 + + tm = time + n * ds + 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) = 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 => pdata%next + end do + + call update_variables(tm, dtm, status) + + end do ! n = 1, stages - 1 + +!= the final step: U(1) = (m-1)/m U(1) + 1/m U(2) + dt/m F[U(1)] +! + tm = time + dt + dtm = fl * dt + + 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) = 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 => pdata%next + end do + + call update_variables(tm, dtm, status) + +! find umax +! + 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 + +! get the maximum error of each variable +! + do l = 1, nf + errors(l) = dt * maxval(abs(lerr(:,l,:,:,:))) + end do + +#ifdef MPI + call reduce_maximum(umax) + call reduce_maximum(errors) +#endif /* MPI */ + +! calculate tolerance and time step +! + emax = maxval(errors) / (atol + rtol * umax) + + if (emax <= 1.0d+00 .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 + end if + + niterations = niterations + 1 + + end do + + if (allocated(lerr)) deallocate(lerr) + +!= final step: U(n+1) = U(1) +! + tm = time + dt + dtm = ft * dt + + pdata => list_data + do while (associated(pdata)) + + pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) + + pdata%u => pdata%uu(:,:,:,:,1) + + 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 */ + +!------------------------------------------------------------------------------- +! + end subroutine evolve_ssprk2_m +! +!=============================================================================== +! ! subroutine EVOLVE_SSP324: ! -------------------------- ! @@ -2454,217 +2847,6 @@ module evolution ! !=============================================================================== ! -! subroutine EVOLVE_SSPRK35: -! ------------------------- -! -! Subroutine advances the solution by one time step using the 3rd order -! 5-stage Strong Stability Preserving Runge-Kutta time integration method. -! -! References: -! -! [1] Ruuth, S. J., -! "Global Optimization of Explicit Strong-Stability-Preserving -! Runge-Kutta methods", -! Mathematics of Computation, -! 2006, vol. 75, no. 253, pp. 183-207 -! -!=============================================================================== -! - subroutine evolve_ssprk35() - - use blocks , only : block_data, list_data - use boundaries, only : boundary_fluxes - use equations , only : ibp, cmax - use sources , only : update_sources - - implicit none - - type(block_data), pointer :: pdata - - integer :: status - real(kind=8) :: ds - real(kind=8) :: tm, dtm - - real(kind=8), parameter :: b1 = 3.77268915331368d-01 - real(kind=8), parameter :: b3 = 2.42995220537396d-01 - real(kind=8), parameter :: b4 = 2.38458932846290d-01 - real(kind=8), parameter :: b5 = 2.87632146308408d-01 - real(kind=8), parameter :: a31 = 3.55909775063327d-01 - real(kind=8), parameter :: a33 = 6.44090224936673d-01 - real(kind=8), parameter :: a41 = 3.67933791638137d-01 - real(kind=8), parameter :: a44 = 6.32066208361863d-01 - real(kind=8), parameter :: a53 = 2.37593836598569d-01 - real(kind=8), parameter :: a55 = 7.62406163401431d-01 -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE - call start_timer(imu) -#endif /* PROFILE */ - -!= 1st step: U(1) = U(n) + b1 dt F[U(n)] -! - ds = b1 * dt - tm = time + ds - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) - - pdata%u => pdata%uu(:,:,:,:,2) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - -!= 2nd step: U(2) = U(1) + b1 dt F[U(1)] -! - tm = time + 2.0d+00 * ds - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - -!= 3rd step: U(3) = a31 U(n) + a33 U(2) + b3 dt F[U(2)] -! - ds = b3 * dt - tm = time + (2.0d+00 * a33 * b1 + b3) * 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = a31 * pdata%uu(:,:,:,:,1) & - + a33 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - -!= 4th step: U(4) = a41 U(n) + a44 U(3) + b4 dt F[U(3)] -! - ds = b4 * dt - tm = time + ((2.0d+00 * b1 * a33 + b3) * a44 + b4) * 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = a41 * pdata%uu(:,:,:,:,1) & - + a44 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) - - pdata%u => pdata%uu(:,:,:,:,1) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - -!= the final step: U(n+1) = a53 U(2) + a55 U(4) + b5 dt F[U(4)] -! - ds = b5 * dt - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = a53 * pdata%uu(:,:,:,:,2) & - + a55 * pdata%uu(:,:,:,:,1) + ds * pdata%du(:,:,:,:) - - 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 */ - -!------------------------------------------------------------------------------- -! - end subroutine evolve_ssprk35 -! -!=============================================================================== -! ! subroutine EVOLVE_SSPRK32M: ! -------------------------- ! @@ -2989,184 +3171,6 @@ module evolution ! !=============================================================================== ! -! subroutine EVOLVE_SSPRK4_10: -! --------------------------- -! -! Subroutine advances the solution by one time step using the 4rd order -! 10-stage Strong Stability Preserving Runge-Kutta time integration method. -! -! References: -! -! [1] Gottlieb, S., Ketcheson, D., Shu C. - W., -! "Strong stability preserving Runge-Kutta and multistep -! time discretizations", -! World Scientific Publishing, 2011 -! -!=============================================================================== -! - subroutine evolve_ssprk4_10() - - use blocks , only : block_data, list_data - use boundaries, only : boundary_fluxes - use equations , only : ibp, cmax - use sources , only : update_sources - - implicit none - - type(block_data), pointer :: pdata - - integer :: n, status - real(kind=8) :: tm, dtm - - real(kind=8), dimension(9) :: ds - - real(kind=8), dimension(9), parameter :: c = & - (/ 1.0d+00, 2.0d+00, 3.0d+00, 4.0d+00, 2.0d+00, & - 3.0d+00, 4.0d+00, 5.0d+00, 6.0d+00 /) / 6.0d+00 -! -!------------------------------------------------------------------------------- -! -#ifdef PROFILE - call start_timer(imu) -#endif /* PROFILE */ - - ds(:) = c(:) * dt - -!= 1st step: U(2) = U(1) -! - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = pdata%uu(:,:,:,:,1) - - pdata => pdata%next - end do - -!= 2nd step: U(1) = [1 + dt/6 L] U(1), for i = 1, ..., 5 -! - do n = 1, 5 - - tm = time + ds(n) - dtm = ds(1) - - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - - end do ! n = 1, 5 - -!= 3rd step: U(2) = U(2)/25 + 9/25 U(1), U(1) = 15 U(2) - 5 U(1) -! - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,2) = (pdata%uu(:,:,:,:,2) & - + 9.0d+00 * pdata%uu(:,:,:,:,1)) / 2.5d+01 - pdata%uu(:,:,:,:,1) = 1.5d+01 * pdata%uu(:,:,:,:,2) & - - 5.0d+00 * pdata%uu(:,:,:,:,1) - - pdata => pdata%next - end do - -!= 4th step: U(1) = [1 + dt/6 L] U(1), for i = 6, ..., 9 -! -! integrate the intermediate steps -! - do n = 6, 9 - - tm = time + ds(n) - dtm = ds(1) - - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,1) + dtm * pdata%du(:,:,:,:) - - pdata => pdata%next - end do - - call update_variables(tm, dtm, status) - - end do ! n = 6, 9 - -!= the final step: U(n+1) = U(2) + 3/5 U(1) + 1/10 dt F[U(1)] -! - tm = time + dt - dtm = dt / 1.0d+01 - - 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() - - pdata => list_data - do while (associated(pdata)) - - pdata%uu(:,:,:,:,1) = pdata%uu(:,:,:,:,2) & - + 6.0d-01 * pdata%uu(:,:,:,:,1) & - + dtm * pdata%du(:,:,:,:) - - 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 */ - -!------------------------------------------------------------------------------- -! - end subroutine evolve_ssprk4_10 -! -!=============================================================================== -! ! subroutine EVOLVE_3SSTARPLUS: ! ---------------------------- !