Merge branch 'master' into binaries
This commit is contained in:
commit
f79846d1bc
@ -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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user