From 78e3b728895bdbf779693adb4ce38c2714d53762 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 31 Dec 2023 17:45:17 -0300 Subject: [PATCH] STATISTICS: Parallelize statistics calculation with OpenMP. Signed-off-by: Grzegorz Kowal --- sources/statistics.F90 | 391 +++++++++++++++++++++-------------------- 1 file changed, 196 insertions(+), 195 deletions(-) diff --git a/sources/statistics.F90 b/sources/statistics.F90 index aa2b849..f4a3c94 100644 --- a/sources/statistics.F90 +++ b/sources/statistics.F90 @@ -450,10 +450,10 @@ module statistics ! subroutine store_statistics() - use blocks , only : block_leaf, block_meta, block_data + use blocks , only : block_leaf, block_meta, block_data, data_blocks use blocks , only : list_leaf, list_data - use blocks , only : get_mblocks, get_nleafs - use coordinates, only : ni => ncells, nn => bcells, nb, ne, nbm, nep + use blocks , only : get_mblocks, get_nleafs, get_dblocks + use coordinates, only : nc => ncells, nn => bcells, nb, ne, nbm, nep use coordinates, only : advol, voli use coordinates, only : toplev #if NDIMS == 3 @@ -481,52 +481,41 @@ module statistics implicit none - logical, save :: first = .true. - integer :: n, l, u, nk, nl, nm, status -#if NDIMS == 3 - real(kind=8) :: xl, xu, yl, yu, zl, zu -#else /* NDIMS == 3 */ - real(kind=8) :: xl, xu, yl, yu -#endif /* NDIMS == 3 */ - real(kind=8) :: dvol, dvolh -#if NDIMS == 3 - real(kind=8) :: dyz, dxz, dxy -#else /* NDIMS == 3 */ - real(kind=8) :: dyz, dxz -#endif /* NDIMS == 3 */ + character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()' + integer , parameter :: narr = 32 + real(kind=8) , parameter :: eps = epsilon(1.0d+00) + real(kind=8) , parameter :: big = huge(1.0d+00) type(block_leaf), pointer :: pleaf type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata - real(kind=8), dimension(3) :: dh + integer(kind=4), dimension(toplev) :: ldist +#ifdef MPI + integer(kind=4), dimension(nprocs) :: cdist +#endif /* MPI */ + real(kind=8), dimension(narr) :: gint, gavg, gmin, gmax + real(kind=8), dimension(narr) :: lint, lavg, lmin, lmax + real(kind=8), dimension(3) :: dh + + integer :: nblks + integer :: n, l, u, nl, nm, status + real(kind=8) :: xl, xu, yl, yu, zl, zu + real(kind=8) :: dvol, dvolh + real(kind=8) :: dyz, dxz, dxy + + logical, save :: first = .true. + integer, save :: nt = 0 real(kind=8), dimension(:,:,:,:), pointer, save :: jc real(kind=8), dimension(:,:,:) , pointer, save :: qq real(kind=8), dimension(:,:,:) , pointer, save :: vel, mag, sqd, tmp - integer , parameter :: narr = 32 - - real(kind=8), dimension(narr) :: inarr, avarr, mnarr, mxarr - - real(kind=8) , parameter :: eps = epsilon(1.0d+00) - real(kind=8) , parameter :: big = huge(1.0d+00) - - integer(kind=4), dimension(toplev) :: ldist -#ifdef MPI - integer(kind=4), dimension(nprocs) :: cdist -#endif /* MPI */ - - integer, save :: nt +!$omp threadprivate(first, nt, jc, qq, vel, mag, sqd, tmp) !$ integer :: omp_get_thread_num -!$omp threadprivate(first, nt, jc, qq) - - character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()' !------------------------------------------------------------------------------- ! - nt = 0 -!$ nt = omp_get_thread_num() ! process and store the mesh statistics only on the master node ! if (master) then @@ -574,20 +563,34 @@ module statistics ! if (mod(step, iintd) > 0) return +! reset the integrals array +! + gint(:) = 0.0d+00 + gavg(:) = 0.0d+00 + gmin(:) = big + gmax(:) = - big + +! reset some statistics if they are not used +! + if (ipr < 1) then + gmin(2) = 0.0d+00 + gmax(2) = 0.0d+00 + end if + if (.not. magnetized) then + gmin(4) = 0.0d+00 + gmin(5) = 0.0d+00 + gmin(7) = 0.0d+00 + gmax(4) = 0.0d+00 + gmax(5) = 0.0d+00 + gmax(7) = 0.0d+00 + end if + + nblks = get_dblocks() + +!$omp parallel default(shared) private(pdata,pmeta,lint,lavg,lmin,lmax,dvol,dvolh,n,l,u,dh,xl,xu,yl,yu,zl,zu,dxy,dxz,dyz) if (first) then - n = 4 * nn**NDIMS + (nv + 1) * ni**(NDIMS - 1) + 4 * ni**NDIMS - - call resize_workspace(n, status) - if (status /= 0) then - call print_message(loc, "Could not resize the workspace!") - go to 100 - end if - -#if NDIMS == 3 - nk = ni -#else /* NDIMS == 3 */ - nk = 1 -#endif /* NDIMS == 3 */ +!$ nt = omp_get_thread_num() + n = 4 * nn**NDIMS + (nv + 1) * nc**(NDIMS - 1) + 4 * nc**NDIMS l = 1 u = l - 1 + 4 * nn**NDIMS @@ -597,64 +600,56 @@ module statistics jc(0:3,1:nn,1:nn,1: 1) => work(l:u,nt) #endif /* NDIMS == 3 */ l = u + 1 - u = l - 1 + (nv + 1) * ni**(NDIMS - 1) + u = l - 1 + (nv + 1) * nc**(NDIMS - 1) #if NDIMS == 3 - qq(0:nv,1:ni,1:ni) => work(l:u,nt) + qq(0:nv,1:nc,1:nc) => work(l:u,nt) #else /* NDIMS == 3 */ - qq(0:nv,1:ni,1: 1) => work(l:u,nt) + qq(0:nv,1:nc,1: 1) => work(l:u,nt) #endif /* NDIMS == 3 */ - n = ni**NDIMS + n = nc**NDIMS l = u + 1 u = u + n - vel(1:ni,1:ni,1:nk) => work(l:u,nt) +#if NDIMS == 3 + vel(1:nc,1:nc,1:nc) => work(l:u,nt) +#else /* NDIMS == 3 */ + vel(1:nc,1:nc,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ l = u + 1 u = u + n - mag(1:ni,1:ni,1:nk) => work(l:u,nt) +#if NDIMS == 3 + mag(1:nc,1:nc,1:nc) => work(l:u,nt) +#else /* NDIMS == 3 */ + mag(1:nc,1:nc,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ l = u + 1 u = u + n - sqd(1:ni,1:ni,1:nk) => work(l:u,nt) +#if NDIMS == 3 + sqd(1:nc,1:nc,1:nc) => work(l:u,nt) +#else /* NDIMS == 3 */ + sqd(1:nc,1:nc,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ l = u + 1 u = u + n - tmp(1:ni,1:ni,1:nk) => work(l:u,nt) +#if NDIMS == 3 + tmp(1:nc,1:nc,1:nc) => work(l:u,nt) +#else /* NDIMS == 3 */ + tmp(1:nc,1:nc,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ first = .false. end if -! reset the integrals array -! - inarr(:) = 0.0d+00 - avarr(:) = 0.0d+00 - mnarr(:) = big - mxarr(:) = - big - -! reset some statistics if they are not used -! - if (ipr < 1) then - mnarr(2) = 0.0d+00 - mxarr(2) = 0.0d+00 - end if - if (.not. magnetized) then - mnarr(4) = 0.0d+00 - mxarr(4) = 0.0d+00 - mnarr(5) = 0.0d+00 - mxarr(5) = 0.0d+00 - mnarr(7) = 0.0d+00 - mxarr(7) = 0.0d+00 - end if - if (work_in_use(nt)) & call print_message(loc, "Workspace is being used right now! " // & "Corruptions can occur!") work_in_use(nt) = .true. +!$omp barrier -! associate the pointer with the first block on the data block list -! - pdata => list_data +!$omp do reduction(+:gint) reduction(+:gavg) reduction(min:gmin) reduction(max:gmax) + do n = 1, nblks -! iterate over all data blocks on the list -! - do while(associated(pdata)) + pdata => data_blocks(n)%ptr pmeta => pdata%meta dvol = advol(pmeta%level) @@ -667,9 +662,10 @@ module statistics #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(1) = avarr(1) + accsum(tmp(:,:,:)) * dvol - mnarr(1) = min(mnarr(1), minval(tmp(:,:,:))) - mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:))) + + lavg(1) = accsum(tmp) * dvol + lmin(1) = minval(tmp) + lmax(1) = maxval(tmp) ! get average, minimum and maximum values of pressure ! @@ -679,9 +675,9 @@ module statistics #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(2) = avarr(2) + accsum(tmp(:,:,:)) * dvol - mnarr(2) = min(mnarr(2), minval(tmp(:,:,:))) - mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:))) + lavg(2) = accsum(tmp) * dvol + lmin(2) = minval(tmp) + lmax(2) = maxval(tmp) end if ! get average, minimum and maximum values of velocity amplitude @@ -691,9 +687,9 @@ 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) + accsum(vel(:,:,:)) * dvol - mnarr(3) = min(mnarr(3), minval(vel(:,:,:))) - mxarr(3) = max(mxarr(3), maxval(vel(:,:,:))) + lavg(3) = accsum(vel) * dvol + lmin(3) = minval(vel) + lmax(3) = maxval(vel) ! get average, minimum and maximum values of magnetic field amplitude, and ! divergence potential @@ -704,18 +700,18 @@ 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) + accsum(mag(:,:,:)) * dvol - mnarr(4) = min(mnarr(4), minval(mag(:,:,:))) - mxarr(4) = max(mxarr(4), maxval(mag(:,:,:))) + lavg(4) = accsum(mag) * dvol + lmin(4) = minval(mag) + lmax(4) = maxval(mag) #if NDIMS == 3 tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne,nb:ne) #else /* NDIMS == 3 */ tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne, : ) #endif /* NDIMS == 3 */ - avarr(5) = avarr(5) + accsum(tmp(:,:,:)) * dvol - mnarr(5) = min(mnarr(5), minval(tmp(:,:,:))) - mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:))) + lavg(5) = accsum(tmp) * dvol + lmin(5) = minval(tmp) + lmax(5) = maxval(tmp) end if ! calculate square root of density @@ -738,21 +734,23 @@ module statistics else tmp(:,:,:) = vel(:,:,:) / csnd end if - avarr(6) = avarr(6) + accsum(tmp(:,:,:)) * dvol - mnarr(6) = min(mnarr(6), minval(tmp(:,:,:))) - mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:))) + lavg(6) = accsum(tmp) * dvol + lmin(6) = minval(tmp) + lmax(6) = maxval(tmp) ! get average, minimum and maximum values of Alfvénic Mach number ! if (magnetized) then tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:)) - avarr(7) = avarr(7) + accsum(tmp(:,:,:)) * dvol - mnarr(7) = min(mnarr(7), minval(tmp(:,:,:))) - mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) + lavg(7) = accsum(tmp) * dvol + lmin(7) = minval(tmp) + lmax(7) = maxval(tmp) end if if (track_conservation) then + lint(:) = 0.0d+00 + xl = pmeta%xmin - adx(pmeta%level) xu = pmeta%xmax + adx(pmeta%level) yl = pmeta%ymin - ady(pmeta%level) @@ -775,32 +773,32 @@ module statistics ! total mass ! #if NDIMS == 3 - inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol + lint(1) = accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol + lint(1) = accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ ! sum up kinetic energy ! #if NDIMS == 3 - 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 + lint(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) + 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 + lint(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) + accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol + lint(9) = accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ - inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol + lint(9) = accsum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ end if @@ -808,13 +806,13 @@ module statistics ! if (magnetized) then #if NDIMS == 3 - 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 + lint(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) + 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 + lint(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) @@ -832,18 +830,18 @@ module statistics #endif /* NDIMS == 3 */ qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(2) = inarr(2) + sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz - inarr(6) = inarr(6) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + lint(2) = lint(2) + sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz + lint(6) = lint(6) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz if (ipr > 0) then - inarr(10) = inarr(10) + sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz + lint(10) = lint(10) + sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(14) = inarr(14) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + lint(14) = lint(14) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(17) = inarr(17) - sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz + lint(17) = lint(17) - sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz if (resistivity > 0.0d+00) then #if NDIMS == 3 @@ -853,7 +851,7 @@ module statistics qq(2,:,:) = 5.0d-01 * sum(jc(2,nbm:nb,nb:ne, : ), 1) qq(3,:,:) = 5.0d-01 * sum(jc(3,nbm:nb,nb:ne, : ), 1) #endif /* NDIMS == 3 */ - inarr(20) = inarr(20) & + lint(20) = lint(20) & + sum(qq(2,:,:) * qq(ibz,:,:) & - qq(3,:,:) * qq(iby,:,:)) * dyz end if ! resistivity @@ -868,18 +866,18 @@ module statistics #endif /* NDIMS == 3 */ qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(2) = inarr(2) - sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz - inarr(6) = inarr(6) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + lint(2) = lint(2) - sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz + lint(6) = lint(6) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz if (ipr > 0) then - inarr(10) = inarr(10) - sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz + lint(10) = lint(10) - sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(14) = inarr(14) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + lint(14) = lint(14) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(17) = inarr(17) + sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz + lint(17) = lint(17) + sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz if (resistivity > 0.0d+00) then #if NDIMS == 3 @@ -889,7 +887,7 @@ module statistics qq(2,:,:) = 5.0d-01 * sum(jc(2,ne:nep,nb:ne, : ), 1) qq(3,:,:) = 5.0d-01 * sum(jc(3,ne:nep,nb:ne, : ), 1) #endif /* NDIMS == 3 */ - inarr(20) = inarr(20) & + lint(20) = lint(20) & - sum(qq(2,:,:) * qq(ibz,:,:) & - qq(3,:,:) * qq(iby,:,:)) * dyz end if ! resistivity @@ -906,18 +904,18 @@ module statistics #endif /* NDIMS == 3 */ qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(3) = inarr(3) + sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz - inarr(7) = inarr(7) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + lint(3) = lint(3) + sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz + lint(7) = lint(7) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz if (ipr > 0) then - inarr(11) = inarr(11) + sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz + lint(11) = lint(11) + sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(15) = inarr(15) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + lint(15) = lint(15) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(18) = inarr(18) - sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz + lint(18) = lint(18) - sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz if (resistivity > 0.0d+00) then #if NDIMS == 3 @@ -927,7 +925,7 @@ module statistics qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nbm:nb, : ), 2) qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb, : ), 2) #endif /* NDIMS == 3 */ - inarr(21) = inarr(21) & + lint(21) = lint(21) & + sum(qq(3,:,:) * qq(ibx,:,:) & - qq(1,:,:) * qq(ibz,:,:)) * dxz end if ! resistivity @@ -942,18 +940,18 @@ module statistics #endif /* NDIMS == 3 */ qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(3) = inarr(3) - sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz - inarr(7) = inarr(7) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + lint(3) = lint(3) - sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz + lint(7) = lint(7) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz if (ipr > 0) then - inarr(11) = inarr(11) - sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz + lint(11) = lint(11) - sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(15) = inarr(15) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + lint(15) = lint(15) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(18) = inarr(18) + sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz + lint(18) = lint(18) + sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz if (resistivity > 0.0d+00) then #if NDIMS == 3 @@ -963,7 +961,7 @@ module statistics qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,ne:nep, : ), 2) qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep, : ), 2) #endif /* NDIMS == 3 */ - inarr(21) = inarr(21) & + lint(21) = lint(21) & - sum(qq(3,:,:) * qq(ibx,:,:) & - qq(1,:,:) * qq(ibz,:,:)) * dxz end if ! resistivity @@ -978,23 +976,23 @@ module statistics qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,nbm:nb), 4) qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(4) = inarr(4) + sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy - inarr(8) = inarr(8) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + lint(4) = lint(4) + sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy + lint(8) = lint(8) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy if (ipr > 0) then - inarr(12) = inarr(12) + sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy + lint(12) = lint(12) + sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(16) = inarr(16) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + lint(16) = lint(16) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(19) = inarr(19) - sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy + lint(19) = lint(19) - sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy if (resistivity > 0.0d+00) then qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nb:ne,nbm:nb), 3) qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,nbm:nb), 3) - inarr(22) = inarr(22) & + lint(22) = lint(22) & + sum(qq(1,:,:) * qq(iby,:,:) & - qq(2,:,:) * qq(ibx,:,:)) * dxy end if ! resistivity @@ -1005,23 +1003,23 @@ module statistics qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,ne:nep), 4) qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2, 1) - inarr(4) = inarr(4) - sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy - inarr(8) = inarr(8) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + lint(4) = lint(4) - sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy + lint(8) = lint(8) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy if (ipr > 0) then - inarr(12) = inarr(12) - sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy + lint(12) = lint(12) - sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy end if if (magnetized) then qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1) - inarr(16) = inarr(16) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + lint(16) = lint(16) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1) - inarr(19) = inarr(19) + sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy + lint(19) = lint(19) + sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy if (resistivity > 0.0d+00) then qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nb:ne,ne:nep), 3) qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,ne:nep), 3) - inarr(22) = inarr(22) & + lint(22) = lint(22) & - sum(qq(1,:,:) * qq(iby,:,:) & - qq(2,:,:) * qq(ibx,:,:)) * dxy end if ! resistivity @@ -1034,7 +1032,7 @@ module statistics ! conversion between kinetic and magnetic energies ! - inarr(24) = inarr(24) & + lint(24) = lint(24) & #if NDIMS == 3 - accsum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) & * pdata%q(ibz,nb:ne,nb:ne,nb:ne) & @@ -1048,7 +1046,7 @@ module statistics * pdata%q(iby,nb:ne,nb:ne, : )) & * jc(1,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ - inarr(24) = inarr(24) & + lint(24) = lint(24) & #if NDIMS == 3 - accsum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) & * pdata%q(ibx,nb:ne,nb:ne,nb:ne) & @@ -1062,7 +1060,7 @@ module statistics * pdata%q(ibz,nb:ne,nb:ne, : )) & * jc(2,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ - inarr(24) = inarr(24) & + lint(24) = lint(24) & #if NDIMS == 3 - accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) & * pdata%q(iby,nb:ne,nb:ne,nb:ne) & @@ -1080,7 +1078,7 @@ module statistics ! conversion between magnetic and internal energies ! if (resistivity > 0.0d+00) then - inarr(26) = inarr(26) & + lint(26) = lint(26) & #if NDIMS == 3 - accsum(sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2, 1)) * dvol #else /* NDIMS == 3 */ @@ -1096,7 +1094,7 @@ module statistics ! call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), jc(2,:,:,:)) - inarr(27) = inarr(27) & + lint(27) = lint(27) & #if NDIMS == 3 - accsum(jc(1,nb:ne,nb:ne,nb:ne) & * jc(2,nb:ne,nb:ne,nb:ne)) * dvol @@ -1109,7 +1107,7 @@ module statistics ! call gradient(dh(:), pdata%q(ibp,:,:,:), jc(1:3,:,:,:)) - inarr(28) = inarr(28) & + lint(28) = lint(28) & #if NDIMS == 3 - accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol @@ -1130,7 +1128,7 @@ module statistics call laplace(dh(:), pdata%q(ivy,:,:,:), jc(2,:,:,:)) call laplace(dh(:), pdata%q(ivz,:,:,:), jc(3,:,:,:)) - inarr(25) = inarr(25) & + lint(25) = lint(25) & #if NDIMS == 3 + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & * jc(1:3,nb:ne,nb:ne,nb:ne),1) & @@ -1149,7 +1147,7 @@ module statistics ! call gradient(dh(:), jc(0,:,:,:), jc(1:3,:,:,:)) - inarr(25) = inarr(25) & + lint(25) = lint(25) & #if NDIMS == 3 + accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & * jc(1:3,nb:ne,nb:ne,nb:ne),1) & @@ -1173,7 +1171,7 @@ module statistics ! call gradient(dh(:), pdata%q(ipr,:,:,:), jc(1:3,:,:,:)) - inarr(23) = inarr(23) & + lint(23) = lint(23) & #if NDIMS == 3 - accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & * jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol @@ -1188,7 +1186,7 @@ module statistics ! call gradient(dh(:), pdata%q(idn,:,:,:), jc(1:3,:,:,:)) - inarr(23) = inarr(23) & + lint(23) = lint(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 @@ -1204,71 +1202,76 @@ module statistics ! sum up the injected energy and injection rate ! if (forcing_enabled) then - inarr(29) = einj - inarr(30) = rinj + lint(29) = einj + lint(30) = rinj end if - pdata => pdata%next + gint = gint + lint + gavg = gavg + lavg + gmin = min(gmin, lmin) + gmax = max(gmax, lmax) - end do ! data blocks + end do +!$omp end do nowait work_in_use(nt) = .false. +!$omp end parallel #ifdef MPI ! sum the integral array from all processes ! - call reduce_sum(inarr(1:narr)) + call reduce_sum(gint(1:narr)) ! reduce average, minimum and maximum values ! - call reduce_sum(avarr(1:narr)) - call reduce_minimum(mnarr(1:narr)) - call reduce_maximum(mxarr(1:narr)) + call reduce_sum(gavg(1:narr)) + call reduce_minimum(gmin(1:narr)) + call reduce_maximum(gmax(1:narr)) #endif /* MPI */ if (track_conservation) then if (ipr > 0) then - inarr( 9:12) = inarr( 9:12) / (adiabatic_index - 1.0d+00) - inarr(10:12) = inarr(10:12) * adiabatic_index + gint( 9:12) = gint( 9:12) / (adiabatic_index - 1.0d+00) + gint(10:12) = gint(10:12) * adiabatic_index else - inarr(23) = csnd2 * inarr(23) + gint(23) = csnd2 * gint(23) end if if (viscosity > 0.0d+00) then - inarr(25) = viscosity * inarr(25) + gint(25) = viscosity * gint(25) end if ! resistivity if (magnetized) then if (resistivity > 0.0d+00) then - inarr(20:22) = resistivity * inarr(20:22) - inarr(26) = resistivity * inarr(26) + gint(20:22) = resistivity * gint(20:22) + gint(26) = resistivity * gint(26) end if ! resistivity end if ! magnetized end if ! track conservation ! normalize the averages by the volume of domain ! - avarr(:) = avarr(:) * voli + gavg(:) = gavg(:) * voli ! write down the integrals and statistics to appropriate files ! if (master) then if (track_conservation) then - write(cunit,"(31es25.15e3)") time, inarr(1:30) + write(cunit,"(31es25.15e3)") time, gint(1:30) call flush_and_sync(cunit) end if if (track_statistics) then write(sunit,"(i9,23(1x,1es18.8e3))") step, time, & - avarr(1), mnarr(1), mxarr(1), & - avarr(2), mnarr(2), mxarr(2), & - avarr(3), mnarr(3), mxarr(3), & - avarr(4), mnarr(4), mxarr(4), & - avarr(5), mnarr(5), mxarr(5), & - avarr(6), mnarr(6), mxarr(6), & - avarr(7), mnarr(7), mxarr(7) + gavg(1), gmin(1), gmax(1), & + gavg(2), gmin(2), gmax(2), & + gavg(3), gmin(3), gmax(3), & + gavg(4), gmin(4), gmax(4), & + gavg(5), gmin(5), gmax(5), & + gavg(6), gmin(6), gmax(6), & + gavg(7), gmin(7), gmax(7) call flush_and_sync(sunit) end if if (forcing_enabled) then - write(funit,"(4(1x,1es18.8e3))") time, inarr(29:30), arms + write(funit,"(4(1x,1es18.8e3))") time, gint(29:30), arms call flush_and_sync(funit) end if if (error_control) then @@ -1277,8 +1280,6 @@ module statistics end if end if - 100 continue - !------------------------------------------------------------------------------- ! end subroutine store_statistics