From 39c3cdf7a46a35640974bea03d4910da7c18a6b6 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 28 Aug 2023 17:53:59 -0300 Subject: [PATCH 1/4] EVOLUTION: Improve the conservation in the SSPRK3(2)4 method. The 4th step produced errors large enought to give significant conservation errors. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index d45cc7d..ca7cdff 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -3036,9 +3036,6 @@ module evolution integer :: i, l, n, nrej, status real(kind=8) :: tm, dtm, dh, fc - 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()' !------------------------------------------------------------------------------- @@ -3113,10 +3110,10 @@ module evolution do l = 1, n pdata => data_blocks(l)%ptr - pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1) & - + twothird * pdata%uu(:,:,:,:,2) - pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1) & - + onethird * pdata%uu(:,:,:,:,2) + pdata%uu(:,:,:,:,3) = ( pdata%uu(:,:,:,:,1) & + + 2.0d+00 * pdata%uu(:,:,:,:,2)) / 3.0d+00 + pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) & + + pdata%uu(:,:,:,:,2)) / 3.0d+00 end do !$omp end parallel do @@ -3173,9 +3170,9 @@ module evolution pdata => data_blocks(l)%ptr pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) & - * pdata%uu(ibp,:,:,:,2) + * pdata%uu(ibp,:,:,:,2) pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) & - * pdata%uu(ibp,:,:,:,3) + * pdata%uu(ibp,:,:,:,3) end do !$omp end parallel do end if @@ -3224,8 +3221,8 @@ module evolution call update_variables(tm, dtm, status) if (status /= 0) & - call print_message(loc, "Possible bug: the revert to " // & - "the previous state found it unphysical.") + call print_message(loc, "Possible bug: " // & + "the previous state seems unphysical.") end if niterations = niterations + 1 From 8b36ce58e18bdc5af35c03571e77c5566efb2ae0 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 29 Aug 2023 10:15:08 -0300 Subject: [PATCH 2/4] EVOLUTION: Update all blocks after the mesh update. It seems that the selective block update, i.e., update of the blocks which have been (de)refined only, is not perfect. It violates the numerical conservation somehow, and requires more investigation. In order to fix this, perform the update of all data blocks, no matter if they were refined or not. This makes some blocks to be updated twice, but resolves the problem of variable conservation. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 9 --------- 1 file changed, 9 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index ca7cdff..5a1113b 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1090,7 +1090,6 @@ module evolution ! references ! - use blocks , only : set_blocks_update use coordinates, only : toplev use forcing , only : update_forcing, forcing_enabled use mesh , only : update_mesh @@ -1117,10 +1116,6 @@ module evolution ! if (toplev > 1) then -! set all meta blocks to not be updated -! - call set_blocks_update(.false.) - ! check refinement and refine ! call update_mesh(status) @@ -1131,10 +1126,6 @@ module evolution call update_variables(time + dt, 0.0d+00, status) if (status /= 0) go to 100 -! set all meta blocks to be updated -! - call set_blocks_update(.true.) - end if ! toplev > 1 ! error entry point From 91d3583a78066e57ecd8dd50afc557b87a48fc50 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 29 Aug 2023 10:32:52 -0300 Subject: [PATCH 3/4] STATISTICS: Use Kahan summation for statistics integration. Signed-off-by: Grzegorz Kowal --- sources/statistics.F90 | 160 ++++++++++++++++++++++++++--------------- 1 file changed, 103 insertions(+), 57 deletions(-) diff --git a/sources/statistics.F90 b/sources/statistics.F90 index 577d904..722a38f 100644 --- a/sources/statistics.F90 +++ b/sources/statistics.F90 @@ -667,7 +667,7 @@ module statistics #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol + avarr(1) = avarr(1) + accsum(tmp(:,:,:)) * dvol mnarr(1) = min(mnarr(1), minval(tmp(:,:,:))) mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:))) @@ -679,7 +679,7 @@ module statistics #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol + avarr(2) = avarr(2) + accsum(tmp(:,:,:)) * dvol mnarr(2) = min(mnarr(2), minval(tmp(:,:,:))) mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:))) end if @@ -691,7 +691,7 @@ module statistics #else /* NDIMS == 3 */ vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : )**2, 1)) #endif /* NDIMS == 3 */ - avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol + avarr(3) = avarr(3) + accsum(vel(:,:,:)) * dvol mnarr(3) = min(mnarr(3), minval(vel(:,:,:))) mxarr(3) = max(mxarr(3), maxval(vel(:,:,:))) @@ -704,7 +704,7 @@ module statistics #else /* NDIMS == 3 */ mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : )**2, 1)) #endif /* NDIMS == 3 */ - avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol + avarr(4) = avarr(4) + accsum(mag(:,:,:)) * dvol mnarr(4) = min(mnarr(4), minval(mag(:,:,:))) mxarr(4) = max(mxarr(4), maxval(mag(:,:,:))) @@ -713,7 +713,7 @@ module statistics #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(5) = avarr(5) + sum(tmp(:,:,:)) * dvol + avarr(5) = avarr(5) + accsum(tmp(:,:,:)) * dvol mnarr(5) = min(mnarr(5), minval(tmp(:,:,:))) mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:))) end if @@ -738,7 +738,7 @@ module statistics else tmp(:,:,:) = vel(:,:,:) / csnd end if - avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol + avarr(6) = avarr(6) + accsum(tmp(:,:,:)) * dvol mnarr(6) = min(mnarr(6), minval(tmp(:,:,:))) mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:))) @@ -746,7 +746,7 @@ module statistics ! if (magnetized) then tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:)) - avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol + avarr(7) = avarr(7) + accsum(tmp(:,:,:)) * dvol mnarr(7) = min(mnarr(7), minval(tmp(:,:,:))) mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) end if @@ -775,32 +775,32 @@ module statistics ! total mass ! #if NDIMS == 3 - inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol + inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol + inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ ! sum up kinetic energy ! #if NDIMS == 3 - inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 & - + pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 & - + pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) & - * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh + inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh #else /* NDIMS == 3 */ - inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne, : )**2 & - + pdata%q(ivy,nb:ne,nb:ne, : )**2 & - + pdata%q(ivz,nb:ne,nb:ne, : )**2) & - * pdata%q(idn,nb:ne,nb:ne, : )) * dvolh + inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne, : )**2 & + + pdata%q(ivy,nb:ne,nb:ne, : )**2 & + + pdata%q(ivz,nb:ne,nb:ne, : )**2) & + * pdata%q(idn,nb:ne,nb:ne, : )) * dvolh #endif /* NDIMS == 3 */ ! sum up internal energy ! if (ipr > 0) then #if NDIMS == 3 - inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol + inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol + inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ end if @@ -808,13 +808,13 @@ module statistics ! if (magnetized) then #if NDIMS == 3 - inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 & - + pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 & - + pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh + inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh #else /* NDIMS == 3 */ - inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne, : )**2 & - + pdata%q(iby,nb:ne,nb:ne, : )**2 & - + pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh + inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne, : )**2 & + + pdata%q(iby,nb:ne,nb:ne, : )**2 & + + pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh #endif /* NDIMS == 3 */ ! calculate current density (J = ∇xB) @@ -1036,13 +1036,13 @@ module statistics ! inarr(24) = inarr(24) & #if NDIMS == 3 - - sum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) & + - accsum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) & * pdata%q(ibz,nb:ne,nb:ne,nb:ne) & - pdata%q(ivz,nb:ne,nb:ne,nb:ne) & * pdata%q(iby,nb:ne,nb:ne,nb:ne)) & * jc(1,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - - sum((pdata%q(ivy,nb:ne,nb:ne, : ) & + - accsum((pdata%q(ivy,nb:ne,nb:ne, : ) & * pdata%q(ibz,nb:ne,nb:ne, : ) & - pdata%q(ivz,nb:ne,nb:ne, : ) & * pdata%q(iby,nb:ne,nb:ne, : )) & @@ -1050,13 +1050,13 @@ module statistics #endif /* NDIMS == 3 */ inarr(24) = inarr(24) & #if NDIMS == 3 - - sum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) & + - accsum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) & * pdata%q(ibx,nb:ne,nb:ne,nb:ne) & - pdata%q(ivx,nb:ne,nb:ne,nb:ne) & * pdata%q(ibz,nb:ne,nb:ne,nb:ne)) & * jc(2,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - - sum((pdata%q(ivz,nb:ne,nb:ne, : ) & + - accsum((pdata%q(ivz,nb:ne,nb:ne, : ) & * pdata%q(ibx,nb:ne,nb:ne, : ) & - pdata%q(ivx,nb:ne,nb:ne, : ) & * pdata%q(ibz,nb:ne,nb:ne, : )) & @@ -1064,13 +1064,13 @@ module statistics #endif /* NDIMS == 3 */ inarr(24) = inarr(24) & #if NDIMS == 3 - - sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) & + - accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) & * pdata%q(iby,nb:ne,nb:ne,nb:ne) & - pdata%q(ivy,nb:ne,nb:ne,nb:ne) & * pdata%q(ibx,nb:ne,nb:ne,nb:ne)) & * jc(3,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - - sum((pdata%q(ivx,nb:ne,nb:ne, : ) & + - accsum((pdata%q(ivx,nb:ne,nb:ne, : ) & * pdata%q(iby,nb:ne,nb:ne, : ) & - pdata%q(ivy,nb:ne,nb:ne, : ) & * pdata%q(ibx,nb:ne,nb:ne, : )) & @@ -1082,9 +1082,9 @@ module statistics if (resistivity > 0.0d+00) then inarr(26) = inarr(26) & #if NDIMS == 3 - - sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol + - accsum(sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2, 1)) * dvol #else /* NDIMS == 3 */ - - sum(jc(1:3,nb:ne,nb:ne, : )**2) * dvol + - accsum(sum(jc(1:3,nb:ne,nb:ne, : )**2, 1)) * dvol #endif /* NDIMS == 3 */ end if ! resistive @@ -1098,10 +1098,10 @@ module statistics inarr(27) = inarr(27) & #if NDIMS == 3 - - sum(jc(1,nb:ne,nb:ne,nb:ne) & + - accsum(jc(1,nb:ne,nb:ne,nb:ne) & * jc(2,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - - sum(jc(1,nb:ne,nb:ne, : ) & + - accsum(jc(1,nb:ne,nb:ne, : ) & * jc(2,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ @@ -1111,11 +1111,11 @@ module statistics inarr(28) = inarr(28) & #if NDIMS == 3 - - sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & - * jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol + - accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol #else /* NDIMS == 3 */ - - sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & - * jc(1:3,nb:ne,nb:ne, : )) * dvol + - accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol #endif /* NDIMS == 3 */ end if ! magnetized @@ -1132,13 +1132,13 @@ module statistics inarr(25) = inarr(25) & #if NDIMS == 3 - + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & - * jc(1:3,nb:ne,nb:ne,nb:ne),1) & - * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol + + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne),1) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & - * jc(1:3,nb:ne,nb:ne, : ),1) & - * pdata%q(idn,nb:ne,nb:ne, : )) * dvol + + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ),1) & + * pdata%q(idn,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ ! calculate divergence (∇·v) @@ -1151,15 +1151,15 @@ module statistics inarr(25) = inarr(25) & #if NDIMS == 3 - + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & - * jc(1:3,nb:ne,nb:ne,nb:ne),1) & - * pdata%q(idn,nb:ne,nb:ne,nb:ne)) & - * dvol / 3.0d+00 + + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne),1) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) & + * dvol / 3.0d+00 #else /* NDIMS == 3 */ - + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & - * jc(1:3,nb:ne,nb:ne, : ),1) & - * pdata%q(idn,nb:ne,nb:ne, : )) & - * dvol / 3.0d+00 + + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ),1) & + * pdata%q(idn,nb:ne,nb:ne, : )) & + * dvol / 3.0d+00 #endif /* NDIMS == 3 */ end if @@ -1175,11 +1175,11 @@ module statistics inarr(23) = inarr(23) & #if NDIMS == 3 - - sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & - * jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol + - accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol #else /* NDIMS == 3 */ - - sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & - * jc(1:3,nb:ne,nb:ne, : )) * dvol + - accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol #endif /* NDIMS == 3 */ else @@ -1282,6 +1282,52 @@ module statistics !------------------------------------------------------------------------------- ! end subroutine store_statistics +! +!=============================================================================== +! +! function ACCSUM: +! --------------- +! +! Function performs a more accurate Kahan summation of the input array +! elements. +! +!=============================================================================== +! + real(kind=8) function accsum(q) result(s) + + implicit none + + real(kind=8), dimension(:,:,:), intent(in) :: q + + integer :: i, j, k + real(kind=8) :: c, y, t + +!------------------------------------------------------------------------------- +! + s = 0.0d+00 + c = 0.0d+00 + + do k = 1, size(q, 3) + do j = 1, size(q, 2) + do i = 1, size(q, 1) + t = s + q(i,j,k) + if (abs(s) >= abs(q(i,j,k))) then + c = c + (s - t) + q(i,j,k) + else + c = c + (q(i,j,k) - t) + s + end if + s = t + end do + end do + end do + + s = s + c + + return + +!------------------------------------------------------------------------------- +! + end function accsum !=============================================================================== ! From dfe7310408e1ba15f7fc1ecc074f065ba04d4761 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 29 Aug 2023 10:43:42 -0300 Subject: [PATCH 4/4] EVOLUTION: Improve the conservation in the SSPRK3 and SSPRK(4,3) methods. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 30 +++++++++++++----------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 5a1113b..6f02e67 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -1727,9 +1727,6 @@ module evolution integer :: l, n, status real(kind=8) :: t, ds - real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00 - real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00 - !------------------------------------------------------------------------------- ! status = 1 @@ -1771,8 +1768,8 @@ module evolution != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)] ! - ds = f22 * dt - t = time + 0.5d+00 * dt + ds = dt / 4.0d+00 + t = time + 5.0d-01 * dt !$omp parallel do default(shared) private(pdata) do l = 1, n @@ -1791,8 +1788,9 @@ module evolution do l = 1, n pdata => data_blocks(l)%ptr - pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) & - + f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = (3.0d+00 * pdata%uu(:,:,:,:,1) & + + pdata%uu(:,:,:,:,2) & + + ds * pdata%du(:,:,:,:)) / 4.0d+00 end do !$omp end parallel do @@ -1802,7 +1800,7 @@ module evolution != 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)] ! - ds = f32 * dt + ds = 2.0d+00 * dt / 3.0d+00 t = time + dt !$omp parallel do default(shared) private(pdata) @@ -1822,8 +1820,9 @@ module evolution do l = 1, n pdata => data_blocks(l)%ptr - pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) & - + f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = ( pdata%uu(:,:,:,:,1) & + + 2.0d+00 * (pdata%uu(:,:,:,:,2) & + + ds * pdata%du(:,:,:,:))) / 3.0d+00 pdata%u => pdata%uu(:,:,:,:,1) end do @@ -1845,7 +1844,7 @@ module evolution 100 continue - if (status /= 0) dt = 2.5d-01 * dt + if (status /= 0) dt = dt / 4.0d+00 end do @@ -1895,9 +1894,6 @@ module evolution integer :: l, n, status real(kind=8) :: t, ds - real(kind=8), parameter :: b1 = 2.0d+00 / 3.0d+00, b2 = 1.0d+00 / 3.0d+00, & - b3 = 1.0d+00 / 6.0d+00 - !------------------------------------------------------------------------------- ! status = 1 @@ -1988,9 +1984,9 @@ module evolution do l = 1, n pdata => data_blocks(l)%ptr - pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) & - + b2 * pdata%uu(:,:,:,:,2) & - + b3 * dt * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) & + + pdata%uu(:,:,:,:,2) & + + 5.0d-01 * dt * pdata%du(:,:,:,:)) / 3.0d+00 end do !$omp end parallel do