!!****************************************************************************** !! !! 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 !! !! 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 . !! !!****************************************************************************** !! !! 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) + accsum(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) + accsum(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) + accsum(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) + accsum(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) + accsum(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) + accsum(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) + accsum(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) + accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol #endif /* NDIMS == 3 */ ! sum up kinetic energy ! #if NDIMS == 3 inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 & + pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 & + pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) & * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh #else /* NDIMS == 3 */ inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne, : )**2 & + pdata%q(ivy,nb:ne,nb:ne, : )**2 & + pdata%q(ivz,nb:ne,nb:ne, : )**2) & * pdata%q(idn,nb:ne,nb:ne, : )) * dvolh #endif /* NDIMS == 3 */ ! sum up internal energy ! if (ipr > 0) then #if NDIMS == 3 inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ inarr(9) = inarr(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 inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 & + pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 & + pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh #else /* NDIMS == 3 */ inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne, : )**2 & + pdata%q(iby,nb:ne,nb:ne, : )**2 & + pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh #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 - 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 */ inarr(24) = inarr(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 */ inarr(24) = inarr(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 inarr(26) = inarr(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,:,:,:)) inarr(27) = inarr(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,:,:,:)) inarr(28) = inarr(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,:,:,:)) inarr(25) = inarr(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,:,:,:)) inarr(25) = inarr(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,:,:,:)) inarr(23) = inarr(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,:,:,:)) 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 ! !=============================================================================== ! ! 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