gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
amun-code/sources/statistics.F90
Grzegorz Kowal 29f1e94e9c HELPERS: Add functions to enable/disable I/O flush and sync. 
This change introduces the logical flag 'noflush' to control the flushing
and synchronization of I/O buffers. Subsequently, subroutines
'enable_flush_and_sync' and 'disable_flush_and_sync' were added to enable
and disable the flushing and synchronization, respectively.

The 'flush_and_sync' subroutine was modified to respect the noflush flag.

Additionally, this change updates 'statistics.F90' to read
the 'enable_io_flush' parameter and call the appropriate subroutine
('enable_flush_and_sync' or 'disable_flush_and_sync') based on its value.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
2024-07-19 10:00:18 -03:00

1354 lines
45 KiB
Fortran
Raw Blame History

!!******************************************************************************
!!
!! 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-2024 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 helpers , only : enable_flush_and_sync, disable_flush_and_sync
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 = "off"
call get_parameter("enable_io_flush", stmp)
select case(trim(stmp))
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
call enable_flush_and_sync
case default
call disable_flush_and_sync
end select
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
call get_parameter("statistics_interval", iintd)
iintd = max(1, iintd)
! only master process writes the files to the disk
!
if (master) then
! process parameters
!
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, data_blocks
use blocks , only : list_leaf, list_data
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
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
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
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
!$omp threadprivate(first, nt, jc, qq, vel, mag, sqd, tmp)
!$ integer :: 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
! 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
!$ 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
#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) * nc**(NDIMS - 1)
#if NDIMS == 3
qq(0:nv,1:nc,1:nc) => work(l:u,nt)
#else /* NDIMS == 3 */
qq(0:nv,1:nc,1: 1) => work(l:u,nt)
#endif /* NDIMS == 3 */
n = nc**NDIMS
l = u + 1
u = u + n
#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
#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
#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
#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
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
!$omp do reduction(+:gint) reduction(+:gavg) reduction(min:gmin) reduction(max:gmax)
do n = 1, nblks
pdata => data_blocks(n)%ptr
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 */
lavg(1) = accsum(tmp) * dvol
lmin(1) = minval(tmp)
lmax(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 */
lavg(2) = accsum(tmp) * dvol
lmin(2) = minval(tmp)
lmax(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 */
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
!
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 */
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 */
lavg(5) = accsum(tmp) * dvol
lmin(5) = minval(tmp)
lmax(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
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(:,:,:))
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)
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
lint(1) = accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
lint(1) = accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
! sum up kinetic energy
!
#if NDIMS == 3
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 */
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
lint(9) = accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
lint(9) = accsum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
end if
! sum up magnetic energy
!
if (magnetized) then
#if NDIMS == 3
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 */
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)
!
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)
lint(2) = lint(2) + sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz
lint(6) = lint(6) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz
if (ipr > 0) then
lint(10) = lint(10) + sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(14) = lint(14) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
lint(17) = lint(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 */
lint(20) = lint(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)
lint(2) = lint(2) - sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz
lint(6) = lint(6) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz
if (ipr > 0) then
lint(10) = lint(10) - sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(14) = lint(14) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
lint(17) = lint(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 */
lint(20) = lint(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)
lint(3) = lint(3) + sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz
lint(7) = lint(7) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz
if (ipr > 0) then
lint(11) = lint(11) + sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(15) = lint(15) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
lint(18) = lint(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 */
lint(21) = lint(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)
lint(3) = lint(3) - sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz
lint(7) = lint(7) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz
if (ipr > 0) then
lint(11) = lint(11) - sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(15) = lint(15) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
lint(18) = lint(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 */
lint(21) = lint(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)
lint(4) = lint(4) + sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy
lint(8) = lint(8) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy
if (ipr > 0) then
lint(12) = lint(12) + sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(16) = lint(16) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
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)
lint(22) = lint(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)
lint(4) = lint(4) - sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy
lint(8) = lint(8) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy
if (ipr > 0) then
lint(12) = lint(12) - sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy
end if
if (magnetized) then
qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2, 1)
lint(16) = lint(16) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy
qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:), 1)
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)
lint(22) = lint(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
!
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) &
- 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 */
- accsum((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 */
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) &
- 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 */
- accsum((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 */
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) &
- 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 */
- accsum((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
lint(26) = lint(26) &
#if NDIMS == 3
- accsum(sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2, 1)) * dvol
#else /* NDIMS == 3 */
- accsum(sum(jc(1:3,nb:ne,nb:ne, : )**2, 1)) * 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,:,:,:))
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
#else /* NDIMS == 3 */
- accsum(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,:,:,:))
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
#else /* NDIMS == 3 */
- accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ), 1)) * 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,:,:,:))
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) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
+ accsum(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,:,:,:))
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) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
* dvol / 3.0d+00
#else /* NDIMS == 3 */
+ accsum(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,:,:,:))
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
#else /* NDIMS == 3 */
- accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol
#endif /* NDIMS == 3 */
else
! calculate gradient (∇ρ)
!
call gradient(dh(:), pdata%q(idn,:,:,:), jc(1:3,:,:,:))
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
#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
gint = gint + lint
gavg = gavg + lavg
gmin = min(gmin, lmin)
gmax = max(gmax, lmax)
end do
!$omp end do nowait
work_in_use(nt) = .false.
!$omp end parallel
! sum up the injected energy and injection rate
!
if (forcing_enabled) then
gint(29) = gint(29) + einj
gint(30) = gint(30) + rinj
end if
#ifdef MPI
! sum the integral array from all processes
!
call reduce_sum(gint(1:narr))
! reduce average, minimum and maximum values
!
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
gint( 9:12) = gint( 9:12) / (adiabatic_index - 1.0d+00)
gint(10:12) = gint(10:12) * adiabatic_index
else
gint(23) = csnd2 * gint(23)
end if
if (viscosity > 0.0d+00) then
gint(25) = viscosity * gint(25)
end if ! resistivity
if (magnetized) then
if (resistivity > 0.0d+00) then
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
!
gavg(:) = gavg(:) * voli
! calculate pressure for the isothermal case
!
if (ipr <= 0) then
gavg(2) = csnd2 * gavg(1)
gmin(2) = csnd2 * gmin(1)
gmax(2) = csnd2 * gmax(1)
end if
! write down the integrals and statistics to appropriate files
!
if (master) then
if (track_conservation) then
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, &
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, gint(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
!-------------------------------------------------------------------------------
!
end subroutine store_statistics
!
!===============================================================================
!
! function ACCSUM:
! ---------------
!
! Function performs a more accurate Kahan summation of the input array
! elements.
!
!===============================================================================
!
real(kind=8) function accsum(q) result(s)
implicit none
real(kind=8), dimension(:,:,:), intent(in) :: q
integer :: i, j, k
real(kind=8) :: c, y, t
!-------------------------------------------------------------------------------
!
s = 0.0d+00
c = 0.0d+00
do k = 1, size(q, 3)
do j = 1, size(q, 2)
do i = 1, size(q, 1)
t = s + q(i,j,k)
if (abs(s) >= abs(q(i,j,k))) then
c = c + (s - t) + q(i,j,k)
else
c = c + (q(i,j,k) - t) + s
end if
s = t
end do
end do
end do
s = s + c
return
!-------------------------------------------------------------------------------
!
end function accsum
!===============================================================================
!
end module
Reference in New Issue View Git Blame Copy Permalink