diff --git a/sources/evolution.F90 b/sources/evolution.F90 index d45cc7d..6f02e67 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 @@ -1736,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 @@ -1780,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 @@ -1800,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 @@ -1811,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) @@ -1831,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 @@ -1854,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 @@ -1904,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 @@ -1997,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 @@ -3036,9 +3023,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 +3097,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 +3157,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 +3208,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 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 !=============================================================================== !