Merge branch 'master' into gem-reconnection-challenge
This commit is contained in:
commit
be01bc13fe
@ -66,7 +66,7 @@ if(ENABLE_MPI)
|
||||
endif()
|
||||
endif()
|
||||
|
||||
option(ENABLE_OPENMP "Enable OpenMP support." ON)
|
||||
option(ENABLE_OPENMP "Enable OpenMP support." OFF)
|
||||
if(ENABLE_OPENMP)
|
||||
find_package(OpenMP COMPONENTS Fortran)
|
||||
if(OpenMP_Fortran_FOUND)
|
||||
|
||||
@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
|
||||
|
||||
setuptools.setup(
|
||||
name="amunpy",
|
||||
version="0.9.8",
|
||||
version="0.9.9",
|
||||
author="Grzegorz Kowal",
|
||||
author_email="grzegorz@amuncode.org",
|
||||
description="Python Interface for the AMUN code's snapshots",
|
||||
|
||||
@ -21,6 +21,6 @@ __all__ = [ 'AmunXML', 'AmunH5', 'WriteVTK', \
|
||||
|
||||
__author__ = "Grzegorz Kowal"
|
||||
__copyright__ = "Copyright 2018-2022 Grzegorz Kowal <grzegorz@amuncode.org>"
|
||||
__version__ = "0.9.8"
|
||||
__version__ = "0.9.9"
|
||||
__maintainer__ = "Grzegorz Kowal"
|
||||
__email__ = "grzegorz@amuncode.org"
|
||||
|
||||
@ -316,144 +316,139 @@ class Amun:
|
||||
bz = self.__read_binary_data__('magz', chunk_number)
|
||||
vy = self.__read_binary_data__('vely', chunk_number)
|
||||
vz = self.__read_binary_data__('velz', chunk_number)
|
||||
dset = vy * bz - vz * by
|
||||
dset = vz * by - vy * bz
|
||||
if self.attributes['resistivity'] > 0:
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(bz, -1, axis=iy) - numpy.roll(bz, 1, axis=iy))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(by, 1, axis=1) - numpy.roll(by, -1, axis=1))
|
||||
tmp += (numpy.roll(bz, -1, axis=2) - numpy.roll(bz, 1, axis=2))
|
||||
else:
|
||||
tmp = (numpy.roll(bz, -1, axis=1) - numpy.roll(bz, 1, axis=1))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(by, 1, axis=iz) - numpy.roll(by, -1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dset -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dset += (self.attributes['resistivity'] / 2) * tmp
|
||||
elif dataset == 'eley':
|
||||
bx = self.__read_binary_data__('magx', chunk_number)
|
||||
bz = self.__read_binary_data__('magz', chunk_number)
|
||||
vx = self.__read_binary_data__('velx', chunk_number)
|
||||
vz = self.__read_binary_data__('velz', chunk_number)
|
||||
dset = vz * bx - vx * bz
|
||||
dset = vx * bz - vz * bx
|
||||
if self.attributes['resistivity'] > 0:
|
||||
ix = self.attributes['ndims']
|
||||
tmp = (numpy.roll(bz, 1, axis=ix) - numpy.roll(bz, -1, axis=ix))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(bx, -1, axis=1) - numpy.roll(bx, 1, axis=1))
|
||||
tmp += (numpy.roll(bz, 1, axis=3) - numpy.roll(bz, -1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(bz, 1, axis=2) - numpy.roll(bz, -1, axis=2))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(bx, -1, axis=iz) - numpy.roll(bx, 1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dset -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dset += (self.attributes['resistivity'] / 2) * tmp
|
||||
elif dataset == 'elez':
|
||||
bx = self.__read_binary_data__('magx', chunk_number)
|
||||
by = self.__read_binary_data__('magy', chunk_number)
|
||||
vx = self.__read_binary_data__('velx', chunk_number)
|
||||
vy = self.__read_binary_data__('vely', chunk_number)
|
||||
dset = vx * by - vy * bx
|
||||
dset = vy * bx - vx * by
|
||||
if self.attributes['resistivity'] > 0:
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(bx, 1, axis=2) - numpy.roll(bx, -1, axis=2))
|
||||
tmp += (numpy.roll(by, -1, axis=3) - numpy.roll(by, 1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(bx, 1, axis=1) - numpy.roll(bx, -1, axis=1))
|
||||
tmp += (numpy.roll(by, -1, axis=2) - numpy.roll(by, 1, axis=2))
|
||||
ix = self.attributes['ndims']
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(by, -1, axis=ix) - numpy.roll(by, 1, axis=ix))
|
||||
tmp += (numpy.roll(bx, 1, axis=iy) - numpy.roll(bx, -1, axis=iy))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dset -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dset += (self.attributes['resistivity'] / 2) * tmp
|
||||
elif dataset == 'elec':
|
||||
b1 = self.__read_binary_data__('magy', chunk_number)
|
||||
b2 = self.__read_binary_data__('magz', chunk_number)
|
||||
v1 = self.__read_binary_data__('vely', chunk_number)
|
||||
v2 = self.__read_binary_data__('velz', chunk_number)
|
||||
dtmp = v1 * b2 - v2 * b1
|
||||
dtmp = v2 * b1 - v1 * b2
|
||||
if self.attributes['resistivity'] > 0:
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
|
||||
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
|
||||
else:
|
||||
tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dtmp -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dtmp += (self.attributes['resistivity'] / 2) * tmp
|
||||
dset = dtmp**2
|
||||
|
||||
b1 = self.__read_binary_data__('magx', chunk_number)
|
||||
v1 = self.__read_binary_data__('velx', chunk_number)
|
||||
dtmp = v2 * b1 - v1 * b2
|
||||
dtmp = v1 * b2 - v2 * b1
|
||||
if self.attributes['resistivity'] > 0:
|
||||
ix = self.attributes['ndims']
|
||||
tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1))
|
||||
tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dtmp -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dtmp += (self.attributes['resistivity'] / 2) * tmp
|
||||
dset += dtmp**2
|
||||
|
||||
b2 = self.__read_binary_data__('magy', chunk_number)
|
||||
v2 = self.__read_binary_data__('vely', chunk_number)
|
||||
dtmp = v1 * b2 - v2 * b1
|
||||
dtmp = v2 * b1 - v1 * b2
|
||||
if self.attributes['resistivity'] > 0:
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2))
|
||||
tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
|
||||
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
|
||||
ix = self.attributes['ndims']
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix))
|
||||
tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
dtmp -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
dtmp += (self.attributes['resistivity'] / 2) * tmp
|
||||
|
||||
dset += dtmp**2
|
||||
dset = numpy.sqrt(dset)
|
||||
dset = numpy.sqrt(dset)
|
||||
elif dataset == 'evec':
|
||||
b1 = self.__read_binary_data__('magy', chunk_number)
|
||||
b2 = self.__read_binary_data__('magz', chunk_number)
|
||||
v1 = self.__read_binary_data__('vely', chunk_number)
|
||||
v2 = self.__read_binary_data__('velz', chunk_number)
|
||||
wx = v1 * b2 - v2 * b1
|
||||
wx = v2 * b1 - v1 * b2
|
||||
if self.attributes['resistivity'] > 0:
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
|
||||
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
|
||||
else:
|
||||
tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
wx -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
wx += (self.attributes['resistivity'] / 2) * tmp
|
||||
|
||||
b1 = self.__read_binary_data__('magx', chunk_number)
|
||||
v1 = self.__read_binary_data__('velx', chunk_number)
|
||||
wy = v2 * b1 - v1 * b2
|
||||
wy = v1 * b2 - v2 * b1
|
||||
if self.attributes['resistivity'] > 0:
|
||||
ix = self.attributes['ndims']
|
||||
tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix))
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1))
|
||||
tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2))
|
||||
iz = self.attributes['ndims'] - 2
|
||||
tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
wy -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
wy += (self.attributes['resistivity'] / 2) * tmp
|
||||
|
||||
b2 = self.__read_binary_data__('magy', chunk_number)
|
||||
v2 = self.__read_binary_data__('vely', chunk_number)
|
||||
wz = v1 * b2 - v2 * b1
|
||||
if self.attributes['resistivity'] > 0:
|
||||
if self.attributes['ndims'] == 3:
|
||||
tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2))
|
||||
tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3))
|
||||
else:
|
||||
tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1))
|
||||
tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2))
|
||||
ix = self.attributes['ndims']
|
||||
iy = self.attributes['ndims'] - 1
|
||||
tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix))
|
||||
tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy))
|
||||
|
||||
for p in range(self.chunks[chunk_number]['dblocks']):
|
||||
tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]]
|
||||
wz -= 0.5 * self.attributes['resistivity'] * tmp
|
||||
wz += (self.attributes['resistivity'] / 2) * tmp
|
||||
|
||||
dset = [wx, wy, wz]
|
||||
elif dataset == 'eint':
|
||||
|
||||
@ -316,7 +316,7 @@ module boundaries
|
||||
#ifdef MPI
|
||||
use blocks , only : block_info, pointer_info
|
||||
#endif /* MPI */
|
||||
use blocks , only : ndims, nsides
|
||||
use blocks , only : nsides
|
||||
use coordinates, only : minlev, maxlev
|
||||
use coordinates, only : nh => ncells_half
|
||||
use coordinates, only : nb, ne, nbm, nbp, nem, nep
|
||||
@ -862,7 +862,6 @@ module boundaries
|
||||
!
|
||||
subroutine boundary_variables(t, dt, status)
|
||||
|
||||
use blocks , only : ndims
|
||||
use coordinates, only : minlev, maxlev
|
||||
|
||||
implicit none
|
||||
|
||||
@ -1391,7 +1391,11 @@ module equations
|
||||
use blocks , only : block_data
|
||||
use coordinates, only : nn => bcells
|
||||
use helpers , only : print_message
|
||||
#if NDIMS == 3
|
||||
use coordinates, only : ax, ay, az
|
||||
#else /* NDIMS == 3 */
|
||||
use coordinates, only : ax, ay
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
implicit none
|
||||
|
||||
|
||||
@ -63,6 +63,11 @@ module statistics
|
||||
#endif /* MPI */
|
||||
character(len=32), save :: efmt
|
||||
|
||||
! flag for enabling statistics
|
||||
!
|
||||
logical, save :: track_conservation = .true.
|
||||
logical, save :: track_statistics = .true.
|
||||
|
||||
private
|
||||
|
||||
public :: initialize_statistics, finalize_statistics
|
||||
@ -119,6 +124,24 @@ module statistics
|
||||
!
|
||||
status = 0
|
||||
|
||||
stmp = "on"
|
||||
call get_parameter("enable_conservation", stmp)
|
||||
select case(trim(stmp))
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
track_conservation = .true.
|
||||
case default
|
||||
track_conservation = .false.
|
||||
end select
|
||||
|
||||
stmp = "on"
|
||||
call get_parameter("enable_statistics", stmp)
|
||||
select case(trim(stmp))
|
||||
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
||||
track_statistics = .true.
|
||||
case default
|
||||
track_statistics = .false.
|
||||
end select
|
||||
|
||||
! only master process writes the files to the disk
|
||||
!
|
||||
if (master) then
|
||||
@ -208,78 +231,89 @@ module statistics
|
||||
write(mfmt3,"(a,i0,a)") trim(stmp), nprocs, '(1x,i9))'
|
||||
#endif /* MPI */
|
||||
|
||||
! prepare the conserved variable integrals file
|
||||
! prepare the file for tracking conservation of variables
|
||||
!
|
||||
if (append) then
|
||||
write(fname, "('variable-integrals.dat')")
|
||||
inquire(file=fname, exist=flag)
|
||||
flag = flag .and. nrun > 1
|
||||
else
|
||||
write(fname, "('variable-integrals_',i2.2,'.dat')") nrun
|
||||
flag = .false.
|
||||
end if
|
||||
if (track_conservation) then
|
||||
if (append) then
|
||||
write(fname, "('variable-conservation.dat')")
|
||||
inquire(file=fname, exist=flag)
|
||||
flag = flag .and. nrun > 1
|
||||
else
|
||||
write(fname, "('variable-conservation_',i2.2,'.dat')") nrun
|
||||
flag = .false.
|
||||
end if
|
||||
|
||||
if (flag) then
|
||||
if (flag) then
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=cunit, file=fname, form='formatted', status='old', &
|
||||
position='append', buffered='yes')
|
||||
open(newunit=cunit, file=fname, form='formatted', status='old', &
|
||||
position='append', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=cunit, file=fname, form='formatted', status='old', &
|
||||
position='append')
|
||||
open(newunit=cunit, file=fname, form='formatted', status='old', &
|
||||
position='append')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
write(cunit,"('#')")
|
||||
else
|
||||
else
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=cunit, file=fname, form='formatted', &
|
||||
status='replace', buffered='yes')
|
||||
open(newunit=cunit, file=fname, form='formatted', &
|
||||
status='replace', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=cunit, file=fname, form='formatted', status='replace')
|
||||
open(newunit=cunit, file=fname, form='formatted', status='replace')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
end if
|
||||
end if
|
||||
|
||||
write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt' &
|
||||
, 'mass', 'momx', 'momy', 'momz' &
|
||||
, 'ener', 'ekin', 'emag', 'eint', 'entr'
|
||||
write(cunit,"('#')")
|
||||
write(cunit,"('#',a24,31a25)") 'Time', &
|
||||
'Mass', 'Mass x-adv', 'Mass y-adv', 'Mass z-adv', &
|
||||
'Ekin', 'Ekin x-adv', 'Ekin y-adv', 'Ekin z-adv', &
|
||||
'Eint', 'Eint x-adv', 'Eint y-adv', 'Eint z-adv', &
|
||||
'Emag', 'Emag x-adv', 'Emag y-adv', 'Emag z-adv', &
|
||||
'Emag x-v.b', 'Emag y-v.b', 'Emag z-v.b', &
|
||||
'Emag x-dif', 'Emag y-dif', 'Emag z-dif', &
|
||||
'Ekin-Eint' , 'Ekin-Emag' , 'Ekin-heat' , &
|
||||
'Emag-heat' , 'Emag-divB' , 'Emag-Psi' , &
|
||||
'Einj' , 'Einj-rate'
|
||||
write(cunit,"('#')")
|
||||
end if
|
||||
|
||||
! prepare the variable statistics file
|
||||
!
|
||||
if (append) then
|
||||
write(fname, "('variable-statistics.dat')")
|
||||
inquire(file=fname, exist=flag)
|
||||
flag = flag .and. nrun > 1
|
||||
else
|
||||
write(fname, "('variable-statistics_',i2.2,'.dat')") nrun
|
||||
flag = .false.
|
||||
end if
|
||||
if (track_statistics) then
|
||||
if (append) then
|
||||
write(fname, "('variable-statistics.dat')")
|
||||
inquire(file=fname, exist=flag)
|
||||
flag = flag .and. nrun > 1
|
||||
else
|
||||
write(fname, "('variable-statistics_',i2.2,'.dat')") nrun
|
||||
flag = .false.
|
||||
end if
|
||||
|
||||
if (flag) then
|
||||
if (flag) then
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=sunit, file=fname, form='formatted', status='old', &
|
||||
position='append', buffered='yes')
|
||||
open(newunit=sunit, file=fname, form='formatted', status='old', &
|
||||
position='append', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=sunit, file=fname, form='formatted', status='old', &
|
||||
position='append')
|
||||
open(newunit=sunit, file=fname, form='formatted', status='old', &
|
||||
position='append')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
write(sunit,"('#')")
|
||||
else
|
||||
write(sunit,"('#')")
|
||||
else
|
||||
#ifdef __INTEL_COMPILER
|
||||
open(newunit=sunit, file=fname, form='formatted', &
|
||||
status='replace', buffered='yes')
|
||||
open(newunit=sunit, file=fname, form='formatted', &
|
||||
status='replace', buffered='yes')
|
||||
#else /* __INTEL_COMPILER */
|
||||
open(newunit=sunit, file=fname, form='formatted', status='replace')
|
||||
open(newunit=sunit, file=fname, form='formatted', status='replace')
|
||||
#endif /* __INTEL_COMPILER */
|
||||
end if
|
||||
end if
|
||||
|
||||
write(sunit,"('#',a8,23(1x,a18))") 'step', 'time' &
|
||||
, 'mean(dens)', 'min(dens)', 'max(dens)' &
|
||||
, 'mean(pres)', 'min(pres)', 'max(pres)' &
|
||||
, 'mean(velo)', 'min(velo)', 'max(velo)' &
|
||||
, 'mean(magn)', 'min(magn)', 'max(magn)' &
|
||||
, 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)' &
|
||||
, 'mean(Mson)', 'min(Mson)', 'max(Mson)' &
|
||||
write(sunit,"('#',a8,23(1x,a18))") 'step', 'time' &
|
||||
, 'mean(dens)', 'min(dens)', 'max(dens)' &
|
||||
, 'mean(pres)', 'min(pres)', 'max(pres)' &
|
||||
, 'mean(velo)', 'min(velo)', 'max(velo)' &
|
||||
, 'mean(magn)', 'min(magn)', 'max(magn)' &
|
||||
, 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)' &
|
||||
, 'mean(Mson)', 'min(Mson)', 'max(Mson)' &
|
||||
, 'mean(Malf)', 'min(Malf)', 'max(Malf)'
|
||||
write(sunit,"('#')")
|
||||
write(sunit,"('#')")
|
||||
end if
|
||||
|
||||
! prepare the forcing statistics file
|
||||
!
|
||||
@ -393,10 +427,10 @@ module statistics
|
||||
|
||||
if (master) then
|
||||
close(munit)
|
||||
close(cunit)
|
||||
close(sunit)
|
||||
if (forcing_enabled) close(funit)
|
||||
if (error_control) close(eunit)
|
||||
if (track_conservation) close(cunit)
|
||||
if (track_statistics) close(sunit)
|
||||
if (forcing_enabled) close(funit)
|
||||
if (error_control) close(eunit)
|
||||
end if
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
@ -416,18 +450,23 @@ module statistics
|
||||
!
|
||||
subroutine store_statistics()
|
||||
|
||||
use blocks , only : block_leaf, block_data
|
||||
use blocks , only : block_leaf, block_meta, block_data
|
||||
use blocks , only : list_leaf, list_data
|
||||
use blocks , only : get_mblocks, get_nleafs
|
||||
use coordinates, only : ni => ncells, nb, ne
|
||||
use coordinates, only : ni => ncells, nn => bcells, nb, ne, nbm, nep
|
||||
use coordinates, only : advol, voli
|
||||
use coordinates, only : toplev
|
||||
use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp
|
||||
use equations , only : ien, imx, imy, imz
|
||||
use equations , only : magnetized, adiabatic_index, csnd, csnd2
|
||||
#if NDIMS == 3
|
||||
use coordinates, only : periodic, xmin, xmax, ymin, ymax, zmin, zmax
|
||||
#else /* NDIMS == 3 */
|
||||
use coordinates, only : periodic, xmin, xmax, ymin, ymax
|
||||
#endif /* NDIMS == 3 */
|
||||
use coordinates, only : adx, ady, adz
|
||||
use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
|
||||
use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2
|
||||
use equations , only : errors
|
||||
use evolution , only : error_control
|
||||
use evolution , only : step, time, dt, dth, dte
|
||||
use evolution , only : step, time, dth, dte
|
||||
use forcing , only : forcing_enabled, einj, rinj, arms
|
||||
use helpers , only : print_message, flush_and_sync
|
||||
use mesh , only : changed
|
||||
@ -436,20 +475,37 @@ module statistics
|
||||
use mpitools , only : nprocs
|
||||
use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum
|
||||
#endif /* MPI */
|
||||
use operators , only : curl, divergence, gradient, laplace
|
||||
use sources , only : viscosity, resistivity
|
||||
use workspace , only : resize_workspace, work, work_in_use
|
||||
|
||||
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 */
|
||||
|
||||
type(block_leaf), pointer :: pleaf
|
||||
type(block_meta), pointer :: pmeta
|
||||
type(block_data), pointer :: pdata
|
||||
|
||||
real(kind=8), dimension(:,:,:), pointer, save :: vel, mag, sqd, tmp
|
||||
real(kind=8), dimension(3) :: dh
|
||||
|
||||
integer , parameter :: narr = 16
|
||||
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
|
||||
|
||||
@ -463,7 +519,7 @@ module statistics
|
||||
|
||||
integer, save :: nt
|
||||
!$ integer :: omp_get_thread_num
|
||||
!$omp threadprivate(first, nt)
|
||||
!$omp threadprivate(first, nt, jc, qq)
|
||||
|
||||
character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()'
|
||||
|
||||
@ -519,7 +575,7 @@ module statistics
|
||||
if (mod(step, iintd) > 0) return
|
||||
|
||||
if (first) then
|
||||
n = 4 * ni**NDIMS
|
||||
n = 4 * nn**NDIMS + (nv + 1) * ni**(NDIMS - 1) + 4 * ni**NDIMS
|
||||
|
||||
call resize_workspace(n, status)
|
||||
if (status /= 0) then
|
||||
@ -533,17 +589,31 @@ module statistics
|
||||
nk = 1
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
n = ni**NDIMS
|
||||
l = 1
|
||||
u = n
|
||||
u = l - 1 + 4 * nn**NDIMS
|
||||
#if NDIMS == 3
|
||||
jc(0:3,1:nn,1:nn,1:nn) => work(l:u,nt)
|
||||
#else /* NDIMS == 3 */
|
||||
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)
|
||||
#if NDIMS == 3
|
||||
qq(0:nv,1:ni,1:ni) => work(l:u,nt)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(0:nv,1:ni,1: 1) => work(l:u,nt)
|
||||
#endif /* NDIMS == 3 */
|
||||
n = ni**NDIMS
|
||||
l = u + 1
|
||||
u = u + n
|
||||
vel(1:ni,1:ni,1:nk) => work(l:u,nt)
|
||||
l = l + n
|
||||
l = u + 1
|
||||
u = u + n
|
||||
mag(1:ni,1:ni,1:nk) => work(l:u,nt)
|
||||
l = l + n
|
||||
l = u + 1
|
||||
u = u + n
|
||||
sqd(1:ni,1:ni,1:nk) => work(l:u,nt)
|
||||
l = l + n
|
||||
l = u + 1
|
||||
u = u + n
|
||||
tmp(1:ni,1:ni,1:nk) => work(l:u,nt)
|
||||
|
||||
@ -585,97 +655,11 @@ module statistics
|
||||
! iterate over all data blocks on the list
|
||||
!
|
||||
do while(associated(pdata))
|
||||
pmeta => pdata%meta
|
||||
|
||||
! obtain the volume elements for the current block
|
||||
!
|
||||
dvol = advol(pdata%meta%level)
|
||||
dvol = advol(pmeta%level)
|
||||
dvolh = 0.5d+00 * dvol
|
||||
|
||||
! sum up density and momenta components
|
||||
!
|
||||
#if NDIMS == 3
|
||||
inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvol
|
||||
inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne,nb:ne)) * dvol
|
||||
inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne,nb:ne)) * dvol
|
||||
inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne,nb:ne)) * dvol
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne, : )) * dvol
|
||||
inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne, : )) * dvol
|
||||
inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne, : )) * dvol
|
||||
inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! sum up total energy and entropy
|
||||
!
|
||||
if (ien > 0) then
|
||||
#if NDIMS == 3
|
||||
inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne,nb:ne)) * dvol
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
#if NDIMS == 3
|
||||
tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne,nb:ne) &
|
||||
/ pdata%q(idn,nb:ne,nb:ne,nb:ne)**adiabatic_index) &
|
||||
* pdata%q(idn,nb:ne,nb:ne,nb:ne) / (adiabatic_index - 1.0d+00)
|
||||
#else /* NDIMS == 3 */
|
||||
tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne, : ) &
|
||||
/ pdata%q(idn,nb:ne,nb:ne, : )**adiabatic_index) &
|
||||
* pdata%q(idn,nb:ne,nb:ne, : ) / (adiabatic_index - 1.0d+00)
|
||||
#endif /* NDIMS == 3 */
|
||||
else
|
||||
#if NDIMS == 3
|
||||
tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
|
||||
* pdata%q(idn,nb:ne,nb:ne,nb:ne)
|
||||
#else /* NDIMS == 3 */
|
||||
tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne, : )) &
|
||||
* pdata%q(idn,nb:ne,nb:ne, : )
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
inarr(9) = inarr(9) + sum(tmp(:,:,:)) * dvol
|
||||
|
||||
! sum up kinetic energy
|
||||
!
|
||||
#if NDIMS == 3
|
||||
inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne,nb:ne)**2 &
|
||||
+ pdata%u(imy,nb:ne,nb:ne,nb:ne)**2 &
|
||||
+ pdata%u(imz,nb:ne,nb:ne,nb:ne)**2) &
|
||||
/ pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvolh
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne, : )**2 &
|
||||
+ pdata%u(imy,nb:ne,nb:ne, : )**2 &
|
||||
+ pdata%u(imz,nb:ne,nb:ne, : )**2) &
|
||||
/ pdata%u(idn,nb:ne,nb:ne, : )) * dvolh
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! sum up magnetic energy
|
||||
!
|
||||
if (magnetized) then
|
||||
#if NDIMS == 3
|
||||
inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne,nb:ne)**2 &
|
||||
+ pdata%u(iby,nb:ne,nb:ne,nb:ne)**2 &
|
||||
+ pdata%u(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne, : )**2 &
|
||||
+ pdata%u(iby,nb:ne,nb:ne, : )**2 &
|
||||
+ pdata%u(ibz,nb:ne,nb:ne, : )**2) * dvolh
|
||||
#endif /* NDIMS == 3 */
|
||||
if (ien <= 0) then
|
||||
#if NDIMS == 3
|
||||
inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne,nb:ne)**2) * dvolh
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne, : )**2) * dvolh
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
end if
|
||||
|
||||
! sum up the injected energy and injection rate
|
||||
!
|
||||
if (forcing_enabled) then
|
||||
inarr(11) = einj
|
||||
inarr(12) = rinj
|
||||
end if
|
||||
|
||||
! get average, minimum and maximum values of density
|
||||
!
|
||||
#if NDIMS == 3
|
||||
@ -767,8 +751,462 @@ module statistics
|
||||
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
|
||||
end if
|
||||
|
||||
! associate the pointer with the next block on the list
|
||||
if (track_conservation) then
|
||||
|
||||
xl = pmeta%xmin - adx(pmeta%level)
|
||||
xu = pmeta%xmax + adx(pmeta%level)
|
||||
yl = pmeta%ymin - ady(pmeta%level)
|
||||
yu = pmeta%ymax + ady(pmeta%level)
|
||||
#if NDIMS == 3
|
||||
zl = pmeta%zmin - adz(pmeta%level)
|
||||
zu = pmeta%zmax + adz(pmeta%level)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
dyz = ady(pmeta%level) * adz(pmeta%level)
|
||||
dxz = adx(pmeta%level) * adz(pmeta%level)
|
||||
#if NDIMS == 3
|
||||
dxy = adx(pmeta%level) * ady(pmeta%level)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
dh(1) = adx(pmeta%level)
|
||||
dh(2) = ady(pmeta%level)
|
||||
dh(3) = adz(pmeta%level)
|
||||
|
||||
! calculate current density (J = ∇xB)
|
||||
!
|
||||
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(1:3,:,:,:))
|
||||
|
||||
! total mass
|
||||
!
|
||||
#if NDIMS == 3
|
||||
inarr(1) = inarr(1) + sum(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
|
||||
#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
|
||||
#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
|
||||
#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
|
||||
#else /* NDIMS == 3 */
|
||||
inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
|
||||
! sum up magnetic energy
|
||||
!
|
||||
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
|
||||
#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
|
||||
#endif /* NDIMS == 3 */
|
||||
end if
|
||||
|
||||
if (.not. periodic(1)) then
|
||||
if (xl < xmin) then
|
||||
#if NDIMS == 3
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nbm:nb,nb:ne,nb:ne), 2)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nbm:nb,nb:ne, : ), 2)
|
||||
#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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(10) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(17) = inarr(17) - sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz
|
||||
|
||||
if (resistivity > 0.0d+00) then
|
||||
#if NDIMS == 3
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(2,nbm:nb,nb:ne,nb:ne), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,nbm:nb,nb:ne,nb:ne), 1)
|
||||
#else /* NDIMS == 3 */
|
||||
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) &
|
||||
+ sum(qq(2,:,:) * qq(ibz,:,:) &
|
||||
- qq(3,:,:) * qq(iby,:,:)) * dyz
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
|
||||
if (xu > xmax) then
|
||||
#if NDIMS == 3
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,ne:nep,nb:ne,nb:ne), 2)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,ne:nep,nb:ne, : ), 2)
|
||||
#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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(10) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(17) = inarr(17) + sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz
|
||||
|
||||
if (resistivity > 0.0d+00) then
|
||||
#if NDIMS == 3
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(2,ne:nep,nb:ne,nb:ne), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,ne:nep,nb:ne,nb:ne), 1)
|
||||
#else /* NDIMS == 3 */
|
||||
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) &
|
||||
- sum(qq(2,:,:) * qq(ibz,:,:) &
|
||||
- qq(3,:,:) * qq(iby,:,:)) * dyz
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
end if
|
||||
|
||||
if (.not. periodic(2)) then
|
||||
if (yl < ymin) then
|
||||
#if NDIMS == 3
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nbm:nb,nb:ne), 2)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nbm:nb, : ), 2)
|
||||
#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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(11) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(18) = inarr(18) - sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz
|
||||
|
||||
if (resistivity > 0.0d+00) then
|
||||
#if NDIMS == 3
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nbm:nb,nb:ne), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb,nb:ne), 1)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nbm:nb, : ), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb, : ), 1)
|
||||
#endif /* NDIMS == 3 */
|
||||
inarr(21) = inarr(21) &
|
||||
+ sum(qq(3,:,:) * qq(ibx,:,:) &
|
||||
- qq(1,:,:) * qq(ibz,:,:)) * dxz
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
|
||||
if (yu > ymax) then
|
||||
#if NDIMS == 3
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,ne:nep,nb:ne), 2)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,ne:nep, : ), 2)
|
||||
#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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(11) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(18) = inarr(18) + sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz
|
||||
|
||||
if (resistivity > 0.0d+00) then
|
||||
#if NDIMS == 3
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,ne:nep,nb:ne), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep,nb:ne), 1)
|
||||
#else /* NDIMS == 3 */
|
||||
qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,ne:nep, : ), 1)
|
||||
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep, : ), 1)
|
||||
#endif /* NDIMS == 3 */
|
||||
inarr(21) = inarr(21) &
|
||||
- sum(qq(3,:,:) * qq(ibx,:,:) &
|
||||
- qq(1,:,:) * qq(ibz,:,:)) * dxz
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
end if
|
||||
|
||||
#if NDIMS == 3
|
||||
if (.not. periodic(3)) then
|
||||
if (zl < zmin) then
|
||||
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,nbm:nb), 2)
|
||||
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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(12) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(19) = inarr(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), 1)
|
||||
qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,nbm:nb), 1)
|
||||
inarr(22) = inarr(22) &
|
||||
+ sum(qq(1,:,:) * qq(iby,:,:) &
|
||||
- qq(2,:,:) * qq(ibx,:,:)) * dxy
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
|
||||
if (zu > zmax) then
|
||||
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,ne:nep), 1)
|
||||
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
|
||||
|
||||
if (ipr > 0) then
|
||||
inarr(12) = inarr(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
|
||||
|
||||
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1)
|
||||
inarr(19) = inarr(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), 1)
|
||||
qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,ne:nep), 1)
|
||||
inarr(22) = inarr(22) &
|
||||
- sum(qq(1,:,:) * qq(iby,:,:) &
|
||||
- qq(2,:,:) * qq(ibx,:,:)) * dxy
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if
|
||||
end if
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
if (magnetized) then
|
||||
|
||||
! conversion between kinetic and magnetic energies
|
||||
!
|
||||
inarr(24) = inarr(24) &
|
||||
#if NDIMS == 3
|
||||
- sum((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, : ) &
|
||||
* pdata%q(ibz,nb:ne,nb:ne, : ) &
|
||||
- pdata%q(ivz,nb:ne,nb:ne, : ) &
|
||||
* pdata%q(iby,nb:ne,nb:ne, : )) &
|
||||
* jc(1,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
inarr(24) = inarr(24) &
|
||||
#if NDIMS == 3
|
||||
- sum((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, : ) &
|
||||
* pdata%q(ibx,nb:ne,nb:ne, : ) &
|
||||
- pdata%q(ivx,nb:ne,nb:ne, : ) &
|
||||
* pdata%q(ibz,nb:ne,nb:ne, : )) &
|
||||
* jc(2,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
inarr(24) = inarr(24) &
|
||||
#if NDIMS == 3
|
||||
- sum((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, : ) &
|
||||
* pdata%q(iby,nb:ne,nb:ne, : ) &
|
||||
- pdata%q(ivy,nb:ne,nb:ne, : ) &
|
||||
* pdata%q(ibx,nb:ne,nb:ne, : )) &
|
||||
* jc(3,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! conversion between magnetic and internal energies
|
||||
!
|
||||
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
|
||||
#else /* NDIMS == 3 */
|
||||
- sum(jc(1:3,nb:ne,nb:ne, : )**2) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
end if ! resistive
|
||||
|
||||
! calculate product (v·B)
|
||||
!
|
||||
jc(1,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * pdata%q(ibx:ibz,:,:,:), 1)
|
||||
|
||||
! calculate divergence (∇·B)
|
||||
!
|
||||
call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), jc(2,:,:,:))
|
||||
|
||||
inarr(27) = inarr(27) &
|
||||
#if NDIMS == 3
|
||||
- sum(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, : ) &
|
||||
* jc(2,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! calculate gradient (∇ψ)
|
||||
!
|
||||
call gradient(dh(:), pdata%q(ibp,:,:,:), jc(1:3,:,:,:))
|
||||
|
||||
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
|
||||
#else /* NDIMS == 3 */
|
||||
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
|
||||
* jc(1:3,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
end if ! magnetized
|
||||
|
||||
! viscous heating
|
||||
!
|
||||
if (viscosity > 0.0d+00) then
|
||||
|
||||
! calculate Laplace (∇²v)
|
||||
!
|
||||
call laplace(dh(:), pdata%q(ivx,:,:,:), jc(1,:,:,:))
|
||||
call laplace(dh(:), pdata%q(ivy,:,:,:), jc(2,:,:,:))
|
||||
call laplace(dh(:), pdata%q(ivz,:,:,:), jc(3,:,:,:))
|
||||
|
||||
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
|
||||
#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
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! calculate divergence (∇·v)
|
||||
!
|
||||
call divergence(dh(:), pdata%q(ivx:ivz,:,:,:), jc(0,:,:,:))
|
||||
|
||||
! calculate gradient (∇p)
|
||||
!
|
||||
call gradient(dh(:), jc(0,:,:,:), jc(1:3,:,:,:))
|
||||
|
||||
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
|
||||
#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
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
end if
|
||||
|
||||
! kinetic energy due to pressure gradient and the change of pressure
|
||||
! due to the compression
|
||||
!
|
||||
if (ipr > 0) then
|
||||
|
||||
! calculate gradient (∇p)
|
||||
!
|
||||
call gradient(dh(:), pdata%q(ipr,:,:,:), jc(1:3,:,:,:))
|
||||
|
||||
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
|
||||
#else /* NDIMS == 3 */
|
||||
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
|
||||
* jc(1:3,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
else
|
||||
|
||||
! calculate gradient (∇ρ)
|
||||
!
|
||||
call gradient(dh(:), pdata%q(idn,:,:,:), jc(1:3,:,:,:))
|
||||
|
||||
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
|
||||
#else /* NDIMS == 3 */
|
||||
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
|
||||
* jc(1:3,nb:ne,nb:ne, : )) * dvol
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
end if
|
||||
|
||||
end if ! track conservation
|
||||
|
||||
! sum up the injected energy and injection rate
|
||||
!
|
||||
if (forcing_enabled) then
|
||||
inarr(29) = einj
|
||||
inarr(30) = rinj
|
||||
end if
|
||||
|
||||
pdata => pdata%next
|
||||
|
||||
end do ! data blocks
|
||||
@ -787,13 +1225,23 @@ module statistics
|
||||
call reduce_maximum(mxarr(1:narr))
|
||||
#endif /* MPI */
|
||||
|
||||
! calculate the internal energy, or update the entropy
|
||||
!
|
||||
if (ien > 0) then
|
||||
inarr(8) = inarr(5) - inarr(6) - inarr(7)
|
||||
else
|
||||
inarr(9) = inarr(9) + inarr(6) + inarr(7)
|
||||
end if
|
||||
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
|
||||
else
|
||||
inarr(23) = csnd2 * inarr(23)
|
||||
end if
|
||||
if (viscosity > 0.0d+00) then
|
||||
inarr(25) = viscosity * inarr(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)
|
||||
end if ! resistivity
|
||||
end if ! magnetized
|
||||
end if ! track conservation
|
||||
|
||||
! normalize the averages by the volume of domain
|
||||
!
|
||||
@ -802,21 +1250,24 @@ module statistics
|
||||
! write down the integrals and statistics to appropriate files
|
||||
!
|
||||
if (master) then
|
||||
write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9)
|
||||
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)
|
||||
|
||||
call flush_and_sync(cunit)
|
||||
call flush_and_sync(sunit)
|
||||
if (track_conservation) then
|
||||
write(cunit,"(31es25.15e3)") time, inarr(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)
|
||||
call flush_and_sync(sunit)
|
||||
end if
|
||||
|
||||
if (forcing_enabled) then
|
||||
write(funit,"(4(1x,1es18.8e3))") time, inarr(11:12), arms
|
||||
write(funit,"(4(1x,1es18.8e3))") time, inarr(29:30), arms
|
||||
call flush_and_sync(funit)
|
||||
end if
|
||||
if (error_control) then
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user