1289 lines
44 KiB
Fortran
1289 lines
44 KiB
Fortran
!!******************************************************************************
|
|
!!
|
|
!! This file is part of the AMUN source code, a program to perform
|
|
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
|
|
!! adaptive mesh.
|
|
!!
|
|
!! Copyright (C) 2008-2023 Grzegorz Kowal <grzegorz@amuncode.org>
|
|
!!
|
|
!! This program is free software: you can redistribute it and/or modify
|
|
!! it under the terms of the GNU General Public License as published by
|
|
!! the Free Software Foundation, either version 3 of the License, or
|
|
!! (at your option) any later version.
|
|
!!
|
|
!! This program is distributed in the hope that it will be useful,
|
|
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
!! GNU General Public License for more details.
|
|
!!
|
|
!! You should have received a copy of the GNU General Public License
|
|
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
!!
|
|
!!******************************************************************************
|
|
!!
|
|
!! module: STATISTICS
|
|
!!
|
|
!! This module provides subroutines to calculate and store mesh and variable
|
|
!! statistics.
|
|
!!
|
|
!!******************************************************************************
|
|
!
|
|
module statistics
|
|
|
|
implicit none
|
|
|
|
! MODULE PARAMETERS:
|
|
! =================
|
|
!
|
|
! munit - a file handler for the mesh statistics file;
|
|
! cunit - a file handler for the conserved variable integrals file;
|
|
! sunit - a file handler for the variable statistics file;
|
|
! funit - a file handler for the forcing statistics file;
|
|
! eunit - a file handler for the integration errors file;
|
|
!
|
|
integer(kind=4), save :: munit = 10
|
|
integer(kind=4), save :: cunit = 11
|
|
integer(kind=4), save :: sunit = 12
|
|
integer(kind=4), save :: funit = 13
|
|
integer(kind=4), save :: eunit = 14
|
|
|
|
! iintd - the number of steps between subsequent intervals;
|
|
!
|
|
integer(kind=4), save :: iintd = 1
|
|
|
|
! the mesh coverage and efficiency factors
|
|
!
|
|
real(kind=8), save :: fcv = 1.0d+00, fef = 1.0d+00
|
|
|
|
! the format strings
|
|
!
|
|
character(len=32), save :: mfmt1, mfmt2
|
|
#ifdef MPI
|
|
character(len=32), save :: mfmt3
|
|
#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
|
|
public :: store_statistics
|
|
|
|
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
|
!
|
|
contains
|
|
!
|
|
!===============================================================================
|
|
!!
|
|
!!*** PUBLIC SUBROUTINES *****************************************************
|
|
!!
|
|
!===============================================================================
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine INITIALIZE_STATISTICS:
|
|
! -------------------------------
|
|
!
|
|
! Subroutine initializes module INTEGRALS by preparing the statistics files.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! verbose - the verbose flag;
|
|
! nrun - the number of resumed run;
|
|
! status - the subroutine call status;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine initialize_statistics(verbose, nrun, status)
|
|
|
|
use coordinates, only : toplev, ncells, bcells, nghosts, domain_base_dims
|
|
use equations , only : pvars, nf
|
|
use evolution , only : error_control
|
|
use forcing , only : forcing_enabled
|
|
use mpitools , only : master
|
|
#ifdef MPI
|
|
use mpitools , only : nprocs
|
|
#endif /* MPI */
|
|
use parameters , only : get_parameter
|
|
|
|
implicit none
|
|
|
|
logical, intent(in) :: verbose
|
|
integer, intent(in) :: nrun
|
|
integer, intent(out) :: status
|
|
|
|
character(len=32) :: fname, stmp
|
|
logical :: append, flag
|
|
integer :: l
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
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
|
|
|
|
! process parameters
|
|
!
|
|
call get_parameter("statistics_interval", iintd)
|
|
iintd = max(1, iintd)
|
|
|
|
stmp = "off"
|
|
call get_parameter("statistics_append", stmp)
|
|
select case(trim(stmp))
|
|
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
|
|
append = .true.
|
|
case default
|
|
append = .false.
|
|
end select
|
|
|
|
! prepare the mesh statistics file
|
|
!
|
|
if (append) then
|
|
write(fname,"('mesh-statistics.dat')")
|
|
inquire(file=fname, exist=flag)
|
|
flag = flag .and. nrun > 1
|
|
else
|
|
write(fname,"('mesh-statistics_',i2.2,'.dat')") nrun
|
|
flag = .false.
|
|
end if
|
|
|
|
if (flag) then
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=munit, file=fname, form='formatted', status='old', &
|
|
position='append', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=munit, file=fname, form='formatted', status='old', &
|
|
position='append')
|
|
#endif /* __INTEL_COMPILER */
|
|
write(munit,"('#')")
|
|
else
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=munit, file=fname, form='formatted', &
|
|
status='replace', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=munit, file=fname, form='formatted', status='replace')
|
|
#endif /* __INTEL_COMPILER */
|
|
end if
|
|
|
|
! write the header of the mesh statistics file
|
|
!
|
|
write(munit,"('#',5x,'step',4x,'time',11x,'coverage',7x,'efficiency',"// &
|
|
"5x,'blocks',5x,'leafs',1x,'|',1x," // &
|
|
"'block distribution per level')", advance="no")
|
|
#ifdef MPI
|
|
write(stmp,"(a,i0,a)") "(", max(1, 10 * toplev - 28), "x,'|',1x,a)"
|
|
write(munit,stmp) "block distribution per process"
|
|
#else /* MPI */
|
|
write(munit,"('')")
|
|
#endif /* MPI */
|
|
write(munit,"('#',75x,'|')", advance="no")
|
|
do l = 1, toplev
|
|
write(munit,"(1x,i9)", advance="no") l
|
|
end do
|
|
#ifdef MPI
|
|
write(stmp,"(a,i0,a)") "(", max(1, 10 * (3 - toplev)), "x,'|')"
|
|
write(munit,stmp, advance="no")
|
|
do l = 0, nprocs - 1
|
|
write(munit,"(1x,i9)", advance="no") l
|
|
end do
|
|
#endif /* MPI */
|
|
write(munit,"('')")
|
|
|
|
! calculate the coverage and efficiency factors
|
|
!
|
|
fcv = 1.0d+00 / (product(domain_base_dims) * 2**(NDIMS*(toplev - 1)))
|
|
fef = 1.0d+00
|
|
do l = 1, NDIMS
|
|
fef = fef * (domain_base_dims(l) * ncells * 2**(toplev - 1) &
|
|
+ 2 * nghosts) / bcells
|
|
end do
|
|
|
|
! prepare the format strings for mesh statistics
|
|
!
|
|
mfmt1 = "(1x,i9,3(1x,1es14.6e3),2(1x,i9))"
|
|
write (mfmt2,"(a,i0,a)") '(2x,', toplev, '(1x,i9))'
|
|
#ifdef MPI
|
|
write(stmp,"(a,i0,a)") "(", max(2, 10 * (3 - toplev) + 1), "x,"
|
|
write(mfmt3,"(a,i0,a)") trim(stmp), nprocs, '(1x,i9))'
|
|
#endif /* MPI */
|
|
|
|
! prepare the file for tracking conservation of variables
|
|
!
|
|
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
|
|
#ifdef __INTEL_COMPILER
|
|
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')
|
|
#endif /* __INTEL_COMPILER */
|
|
write(cunit,"('#')")
|
|
else
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=cunit, file=fname, form='formatted', &
|
|
status='replace', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=cunit, file=fname, form='formatted', status='replace')
|
|
#endif /* __INTEL_COMPILER */
|
|
end if
|
|
|
|
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 (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
|
|
#ifdef __INTEL_COMPILER
|
|
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')
|
|
#endif /* __INTEL_COMPILER */
|
|
write(sunit,"('#')")
|
|
else
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=sunit, file=fname, form='formatted', &
|
|
status='replace', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=sunit, file=fname, form='formatted', status='replace')
|
|
#endif /* __INTEL_COMPILER */
|
|
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)' &
|
|
, 'mean(Malf)', 'min(Malf)', 'max(Malf)'
|
|
write(sunit,"('#')")
|
|
end if
|
|
|
|
! prepare the forcing statistics file
|
|
!
|
|
if (forcing_enabled) then
|
|
if (append) then
|
|
write(fname, "('forcing-statistics.dat')")
|
|
inquire(file=fname, exist=flag)
|
|
flag = flag .and. nrun > 1
|
|
else
|
|
write(fname, "('forcing-statistics_',i2.2,'.dat')") nrun
|
|
flag = .false.
|
|
end if
|
|
|
|
if (flag) then
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=funit, file=fname, form='formatted', status='old', &
|
|
position='append', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=funit, file=fname, form='formatted', status='old', &
|
|
position='append')
|
|
#endif /* __INTEL_COMPILER */
|
|
write(cunit,"('#')")
|
|
else
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=funit, file=fname, form='formatted', &
|
|
status='replace', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=funit, file=fname, form='formatted', status='replace')
|
|
#endif /* __INTEL_COMPILER */
|
|
end if
|
|
|
|
write(funit,"('#',a18,4(1x,a18))") 'time', 'einj', 'erat', 'arms'
|
|
write(funit,"('#')")
|
|
end if
|
|
|
|
! prepare the integration errors file
|
|
!
|
|
if (error_control) then
|
|
|
|
if (append) then
|
|
write(fname, "('integration-errors.dat')")
|
|
inquire(file=fname, exist=flag)
|
|
flag = flag .and. nrun > 1
|
|
else
|
|
write(fname, "('integration-errors_',i2.2,'.dat')") nrun
|
|
flag = .false.
|
|
end if
|
|
|
|
! check if the file exists; if not, create a new one, otherwise move to the end
|
|
!
|
|
if (flag) then
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=eunit, file=fname, form='formatted', status='old', &
|
|
position='append', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=eunit, file=fname, form='formatted', status='old', &
|
|
position='append')
|
|
#endif /* __INTEL_COMPILER */
|
|
write(cunit,"('#')")
|
|
else
|
|
#ifdef __INTEL_COMPILER
|
|
open(newunit=eunit, file=fname, form='formatted', &
|
|
status='replace', buffered='yes')
|
|
#else /* __INTEL_COMPILER */
|
|
open(newunit=eunit, file=fname, form='formatted', status='replace')
|
|
#endif /* __INTEL_COMPILER */
|
|
end if
|
|
|
|
write(stmp,"(i9)") nf + 4
|
|
efmt = "('#',a8," // trim(adjustl(stmp)) // "(1x,a18))"
|
|
write(eunit,efmt) 'step', 'time', 'dt_cfl', 'dt_err', &
|
|
'max_err', pvars(1:nf)
|
|
write(eunit,"('#')")
|
|
|
|
efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))"
|
|
|
|
end if ! error_control
|
|
|
|
end if ! master
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine initialize_statistics
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine FINALIZE_STATISTICS:
|
|
! -----------------------------
|
|
!
|
|
! Subroutine finalizes module STATISTICS.
|
|
!
|
|
! Arguments:
|
|
!
|
|
! status - an integer flag for error return value;
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine finalize_statistics(status)
|
|
|
|
use evolution, only : error_control
|
|
use forcing , only : forcing_enabled
|
|
use mpitools , only : master
|
|
|
|
implicit none
|
|
|
|
integer, intent(out) :: status
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
status = 0
|
|
|
|
if (master) then
|
|
close(munit)
|
|
if (track_conservation) close(cunit)
|
|
if (track_statistics) close(sunit)
|
|
if (forcing_enabled) close(funit)
|
|
if (error_control) close(eunit)
|
|
end if
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine finalize_statistics
|
|
!
|
|
!===============================================================================
|
|
!
|
|
! subroutine STORE_STATISTICS:
|
|
! --------------------------
|
|
!
|
|
! Subroutine calculates the processes the variable integrals and statistics,
|
|
! mesh statistics and integration errors, collects them from all processes
|
|
! and stores in the corresponding files.
|
|
!
|
|
!===============================================================================
|
|
!
|
|
subroutine store_statistics()
|
|
|
|
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, nn => bcells, nb, ne, nbm, nep
|
|
use coordinates, only : advol, voli
|
|
use coordinates, only : toplev
|
|
#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, dth, dte
|
|
use forcing , only : forcing_enabled, einj, rinj, arms
|
|
use helpers , only : print_message, flush_and_sync
|
|
use mesh , only : changed
|
|
use mpitools , only : master
|
|
#ifdef MPI
|
|
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(3) :: dh
|
|
|
|
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
|
|
!$ 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
|
|
if (changed) then
|
|
|
|
! get the new numbers of meta blocks and leafs
|
|
!
|
|
nm = get_mblocks()
|
|
nl = get_nleafs()
|
|
|
|
! determine the block distribution per level and per process
|
|
!
|
|
ldist(:) = 0
|
|
#ifdef MPI
|
|
cdist(:) = 0
|
|
#endif /* MPI */
|
|
pleaf => list_leaf
|
|
do while(associated(pleaf))
|
|
l = pleaf%meta%level
|
|
ldist(l) = ldist(l) + 1
|
|
#ifdef MPI
|
|
n = pleaf%meta%process+1
|
|
cdist(n) = cdist(n) + 1
|
|
#endif /* MPI */
|
|
pleaf => pleaf%next
|
|
end do
|
|
|
|
! store the block statistics
|
|
!
|
|
write (munit,mfmt1,advance='no') step, time, fcv * nl, fef / nl, nm, nl
|
|
write (munit,mfmt2,advance='no') ldist(:)
|
|
#ifdef MPI
|
|
write (munit,mfmt3,advance='no') cdist(:)
|
|
#endif /* MPI */
|
|
write(munit,"('')")
|
|
|
|
call flush_and_sync(munit)
|
|
|
|
changed = .false.
|
|
|
|
end if ! number of blocks or leafs changed
|
|
end if ! master
|
|
|
|
! return if the storage interval was not reached
|
|
!
|
|
if (mod(step, iintd) > 0) return
|
|
|
|
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 */
|
|
|
|
l = 1
|
|
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 = u + 1
|
|
u = u + n
|
|
mag(1:ni,1:ni,1:nk) => work(l:u,nt)
|
|
l = u + 1
|
|
u = u + n
|
|
sqd(1:ni,1:ni,1:nk) => work(l:u,nt)
|
|
l = u + 1
|
|
u = u + n
|
|
tmp(1:ni,1:ni,1:nk) => work(l:u,nt)
|
|
|
|
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.
|
|
|
|
! associate the pointer with the first block on the data block list
|
|
!
|
|
pdata => list_data
|
|
|
|
! iterate over all data blocks on the list
|
|
!
|
|
do while(associated(pdata))
|
|
pmeta => pdata%meta
|
|
|
|
dvol = advol(pmeta%level)
|
|
dvolh = 0.5d+00 * dvol
|
|
|
|
! get average, minimum and maximum values of density
|
|
!
|
|
#if NDIMS == 3
|
|
tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne,nb:ne)
|
|
#else /* NDIMS == 3 */
|
|
tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : )
|
|
#endif /* NDIMS == 3 */
|
|
avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol
|
|
mnarr(1) = min(mnarr(1), minval(tmp(:,:,:)))
|
|
mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:)))
|
|
|
|
! get average, minimum and maximum values of pressure
|
|
!
|
|
if (ipr > 0) then
|
|
#if NDIMS == 3
|
|
tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne,nb:ne)
|
|
#else /* NDIMS == 3 */
|
|
tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : )
|
|
#endif /* NDIMS == 3 */
|
|
avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol
|
|
mnarr(2) = min(mnarr(2), minval(tmp(:,:,:)))
|
|
mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:)))
|
|
end if
|
|
|
|
! get average, minimum and maximum values of velocity amplitude
|
|
!
|
|
#if NDIMS == 3
|
|
vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne)**2, 1))
|
|
#else /* NDIMS == 3 */
|
|
vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : )**2, 1))
|
|
#endif /* NDIMS == 3 */
|
|
avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol
|
|
mnarr(3) = min(mnarr(3), minval(vel(:,:,:)))
|
|
mxarr(3) = max(mxarr(3), maxval(vel(:,:,:)))
|
|
|
|
! get average, minimum and maximum values of magnetic field amplitude, and
|
|
! divergence potential
|
|
!
|
|
if (magnetized) then
|
|
#if NDIMS == 3
|
|
mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne)**2, 1))
|
|
#else /* NDIMS == 3 */
|
|
mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : )**2, 1))
|
|
#endif /* NDIMS == 3 */
|
|
avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol
|
|
mnarr(4) = min(mnarr(4), minval(mag(:,:,:)))
|
|
mxarr(4) = max(mxarr(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) + sum(tmp(:,:,:)) * dvol
|
|
mnarr(5) = min(mnarr(5), minval(tmp(:,:,:)))
|
|
mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:)))
|
|
end if
|
|
|
|
! calculate square root of density
|
|
!
|
|
#if NDIMS == 3
|
|
sqd(:,:,:) = sqrt(pdata%q(idn,nb:ne,nb:ne,nb:ne))
|
|
#else /* NDIMS == 3 */
|
|
sqd(:,:,:) = sqrt(pdata%q(idn,nb:ne,nb:ne, : ))
|
|
#endif /* NDIMS == 3 */
|
|
|
|
! get average, minimum and maximum values of sonic Mach number
|
|
!
|
|
if (ipr > 0) then
|
|
tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / &
|
|
#if NDIMS == 3
|
|
sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne,nb:ne))
|
|
#else /* NDIMS == 3 */
|
|
sqrt(adiabatic_index * pdata%q(ipr,nb:ne,nb:ne, : ))
|
|
#endif /* NDIMS == 3 */
|
|
else
|
|
tmp(:,:,:) = vel(:,:,:) / csnd
|
|
end if
|
|
avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol
|
|
mnarr(6) = min(mnarr(6), minval(tmp(:,:,:)))
|
|
mxarr(6) = max(mxarr(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) + sum(tmp(:,:,:)) * dvol
|
|
mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
|
|
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
|
|
end if
|
|
|
|
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)
|
|
|
|
! 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 */
|
|
|
|
! calculate current density (J = ∇xB)
|
|
!
|
|
call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(1: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), 3)
|
|
#else /* NDIMS == 3 */
|
|
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nbm:nb, : ), 3)
|
|
#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), 2)
|
|
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb,nb:ne), 2)
|
|
#else /* NDIMS == 3 */
|
|
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) &
|
|
+ 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), 3)
|
|
#else /* NDIMS == 3 */
|
|
qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,ne:nep, : ), 3)
|
|
#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), 2)
|
|
qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep,nb:ne), 2)
|
|
#else /* NDIMS == 3 */
|
|
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) &
|
|
- 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), 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
|
|
|
|
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), 3)
|
|
qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,nbm:nb), 3)
|
|
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), 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
|
|
|
|
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), 3)
|
|
qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,ne:nep), 3)
|
|
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
|
|
|
|
work_in_use(nt) = .false.
|
|
|
|
#ifdef MPI
|
|
! sum the integral array from all processes
|
|
!
|
|
call reduce_sum(inarr(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))
|
|
#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
|
|
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
|
|
!
|
|
avarr(:) = avarr(:) * 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)
|
|
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(29:30), arms
|
|
call flush_and_sync(funit)
|
|
end if
|
|
if (error_control) then
|
|
write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:)
|
|
call flush_and_sync(eunit)
|
|
end if
|
|
end if
|
|
|
|
100 continue
|
|
|
|
!-------------------------------------------------------------------------------
|
|
!
|
|
end subroutine store_statistics
|
|
|
|
!===============================================================================
|
|
!
|
|
end module
|