gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2023-08-29 10:54:11 -03:00
commit fd3ff031a5
2 changed files with 124 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

54
sources/evolution.F90
View File

@ -1090,7 +1090,6 @@ module evolution
! references ! references
! !
use blocks , only : set_blocks_update
use coordinates, only : toplev use coordinates, only : toplev
use forcing , only : update_forcing, forcing_enabled use forcing , only : update_forcing, forcing_enabled
use mesh , only : update_mesh use mesh , only : update_mesh
@ -1117,10 +1116,6 @@ module evolution
! !
if (toplev > 1) then if (toplev > 1) then
! set all meta blocks to not be updated
!
call set_blocks_update(.false.)
! check refinement and refine ! check refinement and refine
! !
call update_mesh(status) call update_mesh(status)
@ -1131,10 +1126,6 @@ module evolution
call update_variables(time + dt, 0.0d+00, status) call update_variables(time + dt, 0.0d+00, status)
if (status /= 0) go to 100 if (status /= 0) go to 100
! set all meta blocks to be updated
!
call set_blocks_update(.true.)
end if ! toplev > 1 end if ! toplev > 1
! error entry point ! error entry point
@ -1736,9 +1727,6 @@ module evolution
integer :: l, n, status integer :: l, n, status
real(kind=8) :: t, ds 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 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)] != 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
! !
ds = f22 * dt ds = dt / 4.0d+00
t = time + 0.5d+00 * dt t = time + 5.0d-01 * dt
!$omp parallel do default(shared) private(pdata) !$omp parallel do default(shared) private(pdata)
do l = 1, n do l = 1, n
@ -1800,8 +1788,9 @@ module evolution
do l = 1, n do l = 1, n
pdata => data_blocks(l)%ptr pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = (3.0d+00 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)) / 4.0d+00
end do end do
!$omp end parallel 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)] != 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 t = time + dt
!$omp parallel do default(shared) private(pdata) !$omp parallel do default(shared) private(pdata)
@ -1831,8 +1820,9 @@ module evolution
do l = 1, n do l = 1, n
pdata => data_blocks(l)%ptr pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = ( pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:) + 2.0d+00 * (pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:))) / 3.0d+00
pdata%u => pdata%uu(:,:,:,:,1) pdata%u => pdata%uu(:,:,:,:,1)
end do end do
@ -1854,7 +1844,7 @@ module evolution
100 continue 100 continue
if (status /= 0) dt = 2.5d-01 * dt if (status /= 0) dt = dt / 4.0d+00
end do end do
@ -1904,9 +1894,6 @@ module evolution
integer :: l, n, status integer :: l, n, status
real(kind=8) :: t, ds 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 status = 1
@ -1997,9 +1984,9 @@ module evolution
do l = 1, n do l = 1, n
pdata => data_blocks(l)%ptr pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ b2 * pdata%uu(:,:,:,:,2) & + pdata%uu(:,:,:,:,2) &
+ b3 * dt * pdata%du(:,:,:,:) + 5.0d-01 * dt * pdata%du(:,:,:,:)) / 3.0d+00
end do end do
!$omp end parallel do !$omp end parallel do
@ -3036,9 +3023,6 @@ module evolution
integer :: i, l, n, nrej, status integer :: i, l, n, nrej, status
real(kind=8) :: tm, dtm, dh, fc 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()' character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()'
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
@ -3113,10 +3097,10 @@ module evolution
do l = 1, n do l = 1, n
pdata => data_blocks(l)%ptr pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,3) = ( pdata%uu(:,:,:,:,1) &
+ twothird * pdata%uu(:,:,:,:,2) + 2.0d+00 * pdata%uu(:,:,:,:,2)) / 3.0d+00
pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1) & pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ onethird * pdata%uu(:,:,:,:,2) + pdata%uu(:,:,:,:,2)) / 3.0d+00
end do end do
!$omp end parallel do !$omp end parallel do
@ -3224,8 +3208,8 @@ module evolution
call update_variables(tm, dtm, status) call update_variables(tm, dtm, status)
if (status /= 0) & if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // & call print_message(loc, "Possible bug: " // &
"the previous state found it unphysical.") "the previous state seems unphysical.")
end if end if
niterations = niterations + 1 niterations = niterations + 1

120
sources/statistics.F90
View File

@ -667,7 +667,7 @@ module statistics
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : ) tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol avarr(1) = avarr(1) + accsum(tmp(:,:,:)) * dvol
mnarr(1) = min(mnarr(1), minval(tmp(:,:,:))) mnarr(1) = min(mnarr(1), minval(tmp(:,:,:)))
mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:))) mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:)))
@ -679,7 +679,7 @@ module statistics
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : ) tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol avarr(2) = avarr(2) + accsum(tmp(:,:,:)) * dvol
mnarr(2) = min(mnarr(2), minval(tmp(:,:,:))) mnarr(2) = min(mnarr(2), minval(tmp(:,:,:)))
mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:))) mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:)))
end if end if
@ -691,7 +691,7 @@ module statistics
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : )**2, 1)) vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : )**2, 1))
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol avarr(3) = avarr(3) + accsum(vel(:,:,:)) * dvol
mnarr(3) = min(mnarr(3), minval(vel(:,:,:))) mnarr(3) = min(mnarr(3), minval(vel(:,:,:)))
mxarr(3) = max(mxarr(3), maxval(vel(:,:,:))) mxarr(3) = max(mxarr(3), maxval(vel(:,:,:)))
@ -704,7 +704,7 @@ module statistics
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : )**2, 1)) mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : )**2, 1))
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol avarr(4) = avarr(4) + accsum(mag(:,:,:)) * dvol
mnarr(4) = min(mnarr(4), minval(mag(:,:,:))) mnarr(4) = min(mnarr(4), minval(mag(:,:,:)))
mxarr(4) = max(mxarr(4), maxval(mag(:,:,:))) mxarr(4) = max(mxarr(4), maxval(mag(:,:,:)))
@ -713,7 +713,7 @@ module statistics
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne, : ) tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
avarr(5) = avarr(5) + sum(tmp(:,:,:)) * dvol avarr(5) = avarr(5) + accsum(tmp(:,:,:)) * dvol
mnarr(5) = min(mnarr(5), minval(tmp(:,:,:))) mnarr(5) = min(mnarr(5), minval(tmp(:,:,:)))
mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:))) mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:)))
end if end if
@ -738,7 +738,7 @@ module statistics
else else
tmp(:,:,:) = vel(:,:,:) / csnd tmp(:,:,:) = vel(:,:,:) / csnd
end if end if
avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol avarr(6) = avarr(6) + accsum(tmp(:,:,:)) * dvol
mnarr(6) = min(mnarr(6), minval(tmp(:,:,:))) mnarr(6) = min(mnarr(6), minval(tmp(:,:,:)))
mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:))) mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:)))
@ -746,7 +746,7 @@ module statistics
! !
if (magnetized) then if (magnetized) then
tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:)) 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(:,:,:))) mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
end if end if
@ -775,20 +775,20 @@ module statistics
! total mass ! total mass
! !
#if NDIMS == 3 #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 */ #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 */ #endif /* NDIMS == 3 */
! sum up kinetic energy ! sum up kinetic energy
! !
#if NDIMS == 3 #if NDIMS == 3
inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 & 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(ivy,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ivz,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 * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne, : )**2 & inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivy,nb:ne,nb:ne, : )**2 & + pdata%q(ivy,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivz,nb:ne,nb:ne, : )**2) & + pdata%q(ivz,nb:ne,nb:ne, : )**2) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvolh * pdata%q(idn,nb:ne,nb:ne, : )) * dvolh
@ -798,9 +798,9 @@ module statistics
! !
if (ipr > 0) then if (ipr > 0) then
#if NDIMS == 3 #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 */ #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 */ #endif /* NDIMS == 3 */
end if end if
@ -808,11 +808,11 @@ module statistics
! !
if (magnetized) then if (magnetized) then
#if NDIMS == 3 #if NDIMS == 3
inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 & 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(iby,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh + pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne, : )**2 & inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne, : )**2 &
+ pdata%q(iby,nb:ne,nb:ne, : )**2 & + pdata%q(iby,nb:ne,nb:ne, : )**2 &
+ pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh + pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -1036,13 +1036,13 @@ module statistics
! !
inarr(24) = inarr(24) & inarr(24) = inarr(24) &
#if NDIMS == 3 #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(ibz,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivz,nb:ne,nb:ne,nb:ne) & - pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
* pdata%q(iby,nb:ne,nb:ne,nb:ne)) & * pdata%q(iby,nb:ne,nb:ne,nb:ne)) &
* jc(1,nb:ne,nb:ne,nb:ne)) * dvol * jc(1,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */ #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(ibz,nb:ne,nb:ne, : ) &
- pdata%q(ivz,nb:ne,nb:ne, : ) & - pdata%q(ivz,nb:ne,nb:ne, : ) &
* pdata%q(iby,nb:ne,nb:ne, : )) & * pdata%q(iby,nb:ne,nb:ne, : )) &
@ -1050,13 +1050,13 @@ module statistics
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
inarr(24) = inarr(24) & inarr(24) = inarr(24) &
#if NDIMS == 3 #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(ibx,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivx,nb:ne,nb:ne,nb:ne) & - pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibz,nb:ne,nb:ne,nb:ne)) & * pdata%q(ibz,nb:ne,nb:ne,nb:ne)) &
* jc(2,nb:ne,nb:ne,nb:ne)) * dvol * jc(2,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */ #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(ibx,nb:ne,nb:ne, : ) &
- pdata%q(ivx,nb:ne,nb:ne, : ) & - pdata%q(ivx,nb:ne,nb:ne, : ) &
* pdata%q(ibz,nb:ne,nb:ne, : )) & * pdata%q(ibz,nb:ne,nb:ne, : )) &
@ -1064,13 +1064,13 @@ module statistics
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
inarr(24) = inarr(24) & inarr(24) = inarr(24) &
#if NDIMS == 3 #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(iby,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivy,nb:ne,nb:ne,nb:ne) & - pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibx,nb:ne,nb:ne,nb:ne)) & * pdata%q(ibx,nb:ne,nb:ne,nb:ne)) &
* jc(3,nb:ne,nb:ne,nb:ne)) * dvol * jc(3,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */ #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(iby,nb:ne,nb:ne, : ) &
- pdata%q(ivy,nb:ne,nb:ne, : ) & - pdata%q(ivy,nb:ne,nb:ne, : ) &
* pdata%q(ibx,nb:ne,nb:ne, : )) & * pdata%q(ibx,nb:ne,nb:ne, : )) &
@ -1082,9 +1082,9 @@ module statistics
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
inarr(26) = inarr(26) & inarr(26) = inarr(26) &
#if NDIMS == 3 #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 */ #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 */ #endif /* NDIMS == 3 */
end if ! resistive end if ! resistive
@ -1098,10 +1098,10 @@ module statistics
inarr(27) = inarr(27) & inarr(27) = inarr(27) &
#if NDIMS == 3 #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 * jc(2,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
- sum(jc(1,nb:ne,nb:ne, : ) & - accsum(jc(1,nb:ne,nb:ne, : ) &
* jc(2,nb:ne,nb:ne, : )) * dvol * jc(2,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -1111,11 +1111,11 @@ module statistics
inarr(28) = inarr(28) & inarr(28) = inarr(28) &
#if NDIMS == 3 #if NDIMS == 3
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & - accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & - accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : )) * dvol * jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end if ! magnetized end if ! magnetized
@ -1132,11 +1132,11 @@ module statistics
inarr(25) = inarr(25) & inarr(25) = inarr(25) &
#if NDIMS == 3 #if NDIMS == 3
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) & * jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) & * jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvol * pdata%q(idn,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
@ -1151,12 +1151,12 @@ module statistics
inarr(25) = inarr(25) & inarr(25) = inarr(25) &
#if NDIMS == 3 #if NDIMS == 3
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) & * jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) & * pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
* dvol / 3.0d+00 * dvol / 3.0d+00
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) & * jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) & * pdata%q(idn,nb:ne,nb:ne, : )) &
* dvol / 3.0d+00 * dvol / 3.0d+00
@ -1175,11 +1175,11 @@ module statistics
inarr(23) = inarr(23) & inarr(23) = inarr(23) &
#if NDIMS == 3 #if NDIMS == 3
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & - accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol
#else /* NDIMS == 3 */ #else /* NDIMS == 3 */
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & - accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : )) * dvol * jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
else else
@ -1282,6 +1282,52 @@ module statistics
!------------------------------------------------------------------------------- !-------------------------------------------------------------------------------
! !
end subroutine store_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
!=============================================================================== !===============================================================================
! !