From 64c58bb4bdf7e86deac9d9c10cb15a3b7317851d Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:26:19 -0300 Subject: [PATCH 01/10] EQUATIONS: Make compiler happy with imported variables. Signed-off-by: Grzegorz Kowal --- sources/equations.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sources/equations.F90 b/sources/equations.F90 index ce7bef0..ba96e78 100644 --- a/sources/equations.F90 +++ b/sources/equations.F90 @@ -1391,7 +1391,11 @@ module equations use blocks , only : block_data use coordinates, only : nn => bcells use helpers , only : print_message +#if NDIMS == 3 use coordinates, only : ax, ay, az +#else /* NDIMS == 3 */ + use coordinates, only : ax, ay +#endif /* NDIMS == 3 */ implicit none From 336a04377bef02c5bf3c4741323720b33e9c48e7 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:29:58 -0300 Subject: [PATCH 02/10] BOUNDARIES: Remove unused variables. Signed-off-by: Grzegorz Kowal --- sources/boundaries.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sources/boundaries.F90 b/sources/boundaries.F90 index 44a49e0..f698f9d 100644 --- a/sources/boundaries.F90 +++ b/sources/boundaries.F90 @@ -316,7 +316,7 @@ module boundaries #ifdef MPI use blocks , only : block_info, pointer_info #endif /* MPI */ - use blocks , only : ndims, nsides + use blocks , only : nsides use coordinates, only : minlev, maxlev use coordinates, only : nh => ncells_half use coordinates, only : nb, ne, nbm, nbp, nem, nep @@ -862,7 +862,6 @@ module boundaries ! subroutine boundary_variables(t, dt, status) - use blocks , only : ndims use coordinates, only : minlev, maxlev implicit none From 5952bcdb34d42223221c31916a881550a3068cbd Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:31:30 -0300 Subject: [PATCH 03/10] CMAKE: Turn off OpenMP by default. Signed-off-by: Grzegorz Kowal --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 62cc7cc..f77a884 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -66,7 +66,7 @@ if(ENABLE_MPI) endif() endif() -option(ENABLE_OPENMP "Enable OpenMP support." ON) +option(ENABLE_OPENMP "Enable OpenMP support." OFF) if(ENABLE_OPENMP) find_package(OpenMP COMPONENTS Fortran) if(OpenMP_Fortran_FOUND) From 94ad4878502dcbd23945dbd03790504989b63c79 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:32:27 -0300 Subject: [PATCH 04/10] CMAKE: Disable FMA optimizations by default. Signed-off-by: Grzegorz Kowal --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f77a884..e7e936a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -76,7 +76,7 @@ if(ENABLE_OPENMP) endif() endif() -option(DISABLE_FMA "Disable FMA operations for slightly slower, but symmetric floating-point arithmetics." OFF) +option(DISABLE_FMA "Disable FMA operations for slightly slower, but symmetric floating-point arithmetics." ON) if(DISABLE_FMA) if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") target_compile_options(amun.x PRIVATE "-mno-fma;-mno-fma4") From a70a74d2313a2c3eac27ee7b69843196498f93be Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 11:35:03 -0300 Subject: [PATCH 05/10] Revert "CMAKE: Disable FMA optimizations by default." This reverts commit 94ad4878502dcbd23945dbd03790504989b63c79. It breaks CI on ARM64 machine. --- CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e7e936a..f77a884 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -76,7 +76,7 @@ if(ENABLE_OPENMP) endif() endif() -option(DISABLE_FMA "Disable FMA operations for slightly slower, but symmetric floating-point arithmetics." ON) +option(DISABLE_FMA "Disable FMA operations for slightly slower, but symmetric floating-point arithmetics." OFF) if(DISABLE_FMA) if(CMAKE_Fortran_COMPILER_ID MATCHES "GNU") target_compile_options(amun.x PRIVATE "-mno-fma;-mno-fma4") From 296da6a63e0ff4e4c0f076a8591a0377d4b9a1c5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 17 Oct 2022 12:37:20 -0300 Subject: [PATCH 06/10] STATISTICS: Introduce flags to enable/disable statistics collection. Signed-off-by: Grzegorz Kowal --- sources/statistics.F90 | 154 ++++++++++++++++++++++++----------------- 1 file changed, 92 insertions(+), 62 deletions(-) diff --git a/sources/statistics.F90 b/sources/statistics.F90 index 2112dcf..7e00578 100644 --- a/sources/statistics.F90 +++ b/sources/statistics.F90 @@ -63,6 +63,11 @@ module statistics #endif /* MPI */ character(len=32), save :: efmt +! flag for enabling statistics +! + logical, save :: track_conservation = .true. + logical, save :: track_statistics = .true. + private public :: initialize_statistics, finalize_statistics @@ -119,6 +124,24 @@ module statistics ! status = 0 + stmp = "on" + call get_parameter("enable_conservation", stmp) + select case(trim(stmp)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + track_conservation = .true. + case default + track_conservation = .false. + end select + + stmp = "on" + call get_parameter("enable_statistics", stmp) + select case(trim(stmp)) + case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") + track_statistics = .true. + case default + track_statistics = .false. + end select + ! only master process writes the files to the disk ! if (master) then @@ -208,78 +231,82 @@ module statistics write(mfmt3,"(a,i0,a)") trim(stmp), nprocs, '(1x,i9))' #endif /* MPI */ -! prepare the conserved variable integrals file +! prepare the file for tracking conservation of variables ! - if (append) then - write(fname, "('variable-integrals.dat')") - inquire(file=fname, exist=flag) - flag = flag .and. nrun > 1 - else - write(fname, "('variable-integrals_',i2.2,'.dat')") nrun - flag = .false. - end if + if (track_conservation) then + if (append) then + write(fname, "('variable-conservation.dat')") + inquire(file=fname, exist=flag) + flag = flag .and. nrun > 1 + else + write(fname, "('variable-conservation_',i2.2,'.dat')") nrun + flag = .false. + end if - if (flag) then + if (flag) then #ifdef __INTEL_COMPILER - open(newunit=cunit, file=fname, form='formatted', status='old', & - position='append', buffered='yes') + open(newunit=cunit, file=fname, form='formatted', status='old', & + position='append', buffered='yes') #else /* __INTEL_COMPILER */ - open(newunit=cunit, file=fname, form='formatted', status='old', & - position='append') + open(newunit=cunit, file=fname, form='formatted', status='old', & + position='append') #endif /* __INTEL_COMPILER */ write(cunit,"('#')") - else + else #ifdef __INTEL_COMPILER - open(newunit=cunit, file=fname, form='formatted', & - status='replace', buffered='yes') + open(newunit=cunit, file=fname, form='formatted', & + status='replace', buffered='yes') #else /* __INTEL_COMPILER */ - open(newunit=cunit, file=fname, form='formatted', status='replace') + open(newunit=cunit, file=fname, form='formatted', status='replace') #endif /* __INTEL_COMPILER */ - end if + end if - write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt' & - , 'mass', 'momx', 'momy', 'momz' & - , 'ener', 'ekin', 'emag', 'eint', 'entr' - write(cunit,"('#')") + write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt', 'mass', & + 'momx', 'momy', 'momz', 'ener', & + 'ekin', 'emag', 'eint', 'entr' + write(cunit,"('#')") + end if ! prepare the variable statistics file ! - if (append) then - write(fname, "('variable-statistics.dat')") - inquire(file=fname, exist=flag) - flag = flag .and. nrun > 1 - else - write(fname, "('variable-statistics_',i2.2,'.dat')") nrun - flag = .false. - end if + if (track_statistics) then + if (append) then + write(fname, "('variable-statistics.dat')") + inquire(file=fname, exist=flag) + flag = flag .and. nrun > 1 + else + write(fname, "('variable-statistics_',i2.2,'.dat')") nrun + flag = .false. + end if - if (flag) then + if (flag) then #ifdef __INTEL_COMPILER - open(newunit=sunit, file=fname, form='formatted', status='old', & - position='append', buffered='yes') + open(newunit=sunit, file=fname, form='formatted', status='old', & + position='append', buffered='yes') #else /* __INTEL_COMPILER */ - open(newunit=sunit, file=fname, form='formatted', status='old', & - position='append') + open(newunit=sunit, file=fname, form='formatted', status='old', & + position='append') #endif /* __INTEL_COMPILER */ - write(sunit,"('#')") - else + write(sunit,"('#')") + else #ifdef __INTEL_COMPILER - open(newunit=sunit, file=fname, form='formatted', & - status='replace', buffered='yes') + open(newunit=sunit, file=fname, form='formatted', & + status='replace', buffered='yes') #else /* __INTEL_COMPILER */ - open(newunit=sunit, file=fname, form='formatted', status='replace') + open(newunit=sunit, file=fname, form='formatted', status='replace') #endif /* __INTEL_COMPILER */ - end if + end if - write(sunit,"('#',a8,23(1x,a18))") 'step', 'time' & - , 'mean(dens)', 'min(dens)', 'max(dens)' & - , 'mean(pres)', 'min(pres)', 'max(pres)' & - , 'mean(velo)', 'min(velo)', 'max(velo)' & - , 'mean(magn)', 'min(magn)', 'max(magn)' & - , 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)' & - , 'mean(Mson)', 'min(Mson)', 'max(Mson)' & + write(sunit,"('#',a8,23(1x,a18))") 'step', 'time' & + , 'mean(dens)', 'min(dens)', 'max(dens)' & + , 'mean(pres)', 'min(pres)', 'max(pres)' & + , 'mean(velo)', 'min(velo)', 'max(velo)' & + , 'mean(magn)', 'min(magn)', 'max(magn)' & + , 'mean(bpsi)', 'min(bpsi)', 'max(bpsi)' & + , 'mean(Mson)', 'min(Mson)', 'max(Mson)' & , 'mean(Malf)', 'min(Malf)', 'max(Malf)' - write(sunit,"('#')") + write(sunit,"('#')") + end if ! prepare the forcing statistics file ! @@ -802,18 +829,21 @@ module statistics ! write down the integrals and statistics to appropriate files ! if (master) then - write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9) - write(sunit,"(i9,23(1x,1es18.8e3))") step, time & - , avarr(1), mnarr(1), mxarr(1) & - , avarr(2), mnarr(2), mxarr(2) & - , avarr(3), mnarr(3), mxarr(3) & - , avarr(4), mnarr(4), mxarr(4) & - , avarr(5), mnarr(5), mxarr(5) & - , avarr(6), mnarr(6), mxarr(6) & - , avarr(7), mnarr(7), mxarr(7) - - call flush_and_sync(cunit) - call flush_and_sync(sunit) + if (track_conservation) then + write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9) + call flush_and_sync(cunit) + end if + if (track_statistics) then + write(sunit,"(i9,23(1x,1es18.8e3))") step, time, & + avarr(1), mnarr(1), mxarr(1), & + avarr(2), mnarr(2), mxarr(2), & + avarr(3), mnarr(3), mxarr(3), & + avarr(4), mnarr(4), mxarr(4), & + avarr(5), mnarr(5), mxarr(5), & + avarr(6), mnarr(6), mxarr(6), & + avarr(7), mnarr(7), mxarr(7) + call flush_and_sync(sunit) + end if if (forcing_enabled) then write(funit,"(4(1x,1es18.8e3))") time, inarr(11:12), arms From 41bfe4576419335d6ebd7dabb90ade696cddf48c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 23 Oct 2022 15:55:08 -0300 Subject: [PATCH 07/10] PYTHON: Correct calculation of electric field. Signed-off-by: Grzegorz Kowal --- python/amunpy/src/amunpy/amun.py | 115 +++++++++++++++---------------- 1 file changed, 55 insertions(+), 60 deletions(-) diff --git a/python/amunpy/src/amunpy/amun.py b/python/amunpy/src/amunpy/amun.py index 962dcb4..9e612ad 100644 --- a/python/amunpy/src/amunpy/amun.py +++ b/python/amunpy/src/amunpy/amun.py @@ -316,144 +316,139 @@ class Amun: bz = self.__read_binary_data__('magz', chunk_number) vy = self.__read_binary_data__('vely', chunk_number) vz = self.__read_binary_data__('velz', chunk_number) - dset = vy * bz - vz * by + dset = vz * by - vy * bz if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(bz, -1, axis=iy) - numpy.roll(bz, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(by, 1, axis=1) - numpy.roll(by, -1, axis=1)) - tmp += (numpy.roll(bz, -1, axis=2) - numpy.roll(bz, 1, axis=2)) - else: - tmp = (numpy.roll(bz, -1, axis=1) - numpy.roll(bz, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(by, 1, axis=iz) - numpy.roll(by, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'eley': bx = self.__read_binary_data__('magx', chunk_number) bz = self.__read_binary_data__('magz', chunk_number) vx = self.__read_binary_data__('velx', chunk_number) vz = self.__read_binary_data__('velz', chunk_number) - dset = vz * bx - vx * bz + dset = vx * bz - vz * bx if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(bz, 1, axis=ix) - numpy.roll(bz, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(bx, -1, axis=1) - numpy.roll(bx, 1, axis=1)) - tmp += (numpy.roll(bz, 1, axis=3) - numpy.roll(bz, -1, axis=3)) - else: - tmp = (numpy.roll(bz, 1, axis=2) - numpy.roll(bz, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(bx, -1, axis=iz) - numpy.roll(bx, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'elez': bx = self.__read_binary_data__('magx', chunk_number) by = self.__read_binary_data__('magy', chunk_number) vx = self.__read_binary_data__('velx', chunk_number) vy = self.__read_binary_data__('vely', chunk_number) - dset = vx * by - vy * bx + dset = vy * bx - vx * by if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(bx, 1, axis=2) - numpy.roll(bx, -1, axis=2)) - tmp += (numpy.roll(by, -1, axis=3) - numpy.roll(by, 1, axis=3)) - else: - tmp = (numpy.roll(bx, 1, axis=1) - numpy.roll(bx, -1, axis=1)) - tmp += (numpy.roll(by, -1, axis=2) - numpy.roll(by, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(by, -1, axis=ix) - numpy.roll(by, 1, axis=ix)) + tmp += (numpy.roll(bx, 1, axis=iy) - numpy.roll(bx, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dset -= 0.5 * self.attributes['resistivity'] * tmp + dset += (self.attributes['resistivity'] / 2) * tmp elif dataset == 'elec': b1 = self.__read_binary_data__('magy', chunk_number) b2 = self.__read_binary_data__('magz', chunk_number) v1 = self.__read_binary_data__('vely', chunk_number) v2 = self.__read_binary_data__('velz', chunk_number) - dtmp = v1 * b2 - v2 * b1 + dtmp = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) - else: - tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp dset = dtmp**2 b1 = self.__read_binary_data__('magx', chunk_number) v1 = self.__read_binary_data__('velx', chunk_number) - dtmp = v2 * b1 - v1 * b2 + dtmp = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1)) - tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3)) - else: - tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp dset += dtmp**2 b2 = self.__read_binary_data__('magy', chunk_number) v2 = self.__read_binary_data__('vely', chunk_number) - dtmp = v1 * b2 - v2 * b1 + dtmp = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2)) - tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3)) - else: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix)) + tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - dtmp -= 0.5 * self.attributes['resistivity'] * tmp + dtmp += (self.attributes['resistivity'] / 2) * tmp + dset += dtmp**2 - dset = numpy.sqrt(dset) + dset = numpy.sqrt(dset) elif dataset == 'evec': b1 = self.__read_binary_data__('magy', chunk_number) b2 = self.__read_binary_data__('magz', chunk_number) v1 = self.__read_binary_data__('vely', chunk_number) v2 = self.__read_binary_data__('velz', chunk_number) - wx = v1 * b2 - v2 * b1 + wx = v2 * b1 - v1 * b2 if self.attributes['resistivity'] > 0: + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=iy) - numpy.roll(b2, 1, axis=iy)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) - else: - tmp = (numpy.roll(b2, -1, axis=1) - numpy.roll(b2, 1, axis=1)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, 1, axis=iz) - numpy.roll(b1, -1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wx -= 0.5 * self.attributes['resistivity'] * tmp + wx += (self.attributes['resistivity'] / 2) * tmp b1 = self.__read_binary_data__('magx', chunk_number) v1 = self.__read_binary_data__('velx', chunk_number) - wy = v2 * b1 - v1 * b2 + wy = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: + ix = self.attributes['ndims'] + tmp = (numpy.roll(b2, 1, axis=ix) - numpy.roll(b2, -1, axis=ix)) if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, -1, axis=1) - numpy.roll(b1, 1, axis=1)) - tmp += (numpy.roll(b2, 1, axis=3) - numpy.roll(b2, -1, axis=3)) - else: - tmp = (numpy.roll(b2, 1, axis=2) - numpy.roll(b2, -1, axis=2)) + iz = self.attributes['ndims'] - 2 + tmp += (numpy.roll(b1, -1, axis=iz) - numpy.roll(b1, 1, axis=iz)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wy -= 0.5 * self.attributes['resistivity'] * tmp + wy += (self.attributes['resistivity'] / 2) * tmp b2 = self.__read_binary_data__('magy', chunk_number) v2 = self.__read_binary_data__('vely', chunk_number) wz = v1 * b2 - v2 * b1 if self.attributes['resistivity'] > 0: - if self.attributes['ndims'] == 3: - tmp = (numpy.roll(b1, 1, axis=2) - numpy.roll(b1, -1, axis=2)) - tmp += (numpy.roll(b2, -1, axis=3) - numpy.roll(b2, 1, axis=3)) - else: - tmp = (numpy.roll(b1, 1, axis=1) - numpy.roll(b1, -1, axis=1)) - tmp += (numpy.roll(b2, -1, axis=2) - numpy.roll(b2, 1, axis=2)) + ix = self.attributes['ndims'] + iy = self.attributes['ndims'] - 1 + tmp = (numpy.roll(b2, -1, axis=ix) - numpy.roll(b2, 1, axis=ix)) + tmp += (numpy.roll(b1, 1, axis=iy) - numpy.roll(b1, -1, axis=iy)) for p in range(self.chunks[chunk_number]['dblocks']): tmp[p,...] /= self.cell_size[self.chunks[chunk_number]['levels'][p]] - wz -= 0.5 * self.attributes['resistivity'] * tmp + wz += (self.attributes['resistivity'] / 2) * tmp dset = [wx, wy, wz] elif dataset == 'eint': From f539dbe3bb24cdfcaf9c7376cc01f9d439b80da2 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 23 Oct 2022 15:56:43 -0300 Subject: [PATCH 08/10] PYTHON: Increase amunpy version to 0.9.9. Signed-off-by: Grzegorz Kowal --- python/amunpy/setup.py | 2 +- python/amunpy/src/amunpy/__init__.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/python/amunpy/setup.py b/python/amunpy/setup.py index bb1e62d..7aab18d 100644 --- a/python/amunpy/setup.py +++ b/python/amunpy/setup.py @@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh: setuptools.setup( name="amunpy", - version="0.9.8", + version="0.9.9", author="Grzegorz Kowal", author_email="grzegorz@amuncode.org", description="Python Interface for the AMUN code's snapshots", diff --git a/python/amunpy/src/amunpy/__init__.py b/python/amunpy/src/amunpy/__init__.py index 3ead146..18e18d5 100644 --- a/python/amunpy/src/amunpy/__init__.py +++ b/python/amunpy/src/amunpy/__init__.py @@ -21,6 +21,6 @@ __all__ = [ 'AmunXML', 'AmunH5', 'WriteVTK', \ __author__ = "Grzegorz Kowal" __copyright__ = "Copyright 2018-2022 Grzegorz Kowal " -__version__ = "0.9.8" +__version__ = "0.9.9" __maintainer__ = "Grzegorz Kowal" __email__ = "grzegorz@amuncode.org" From fbfc1983e6a192a4349b5d5d9b96f9d029382c3c Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Sun, 23 Oct 2022 22:42:17 -0300 Subject: [PATCH 09/10] STATISTICS: Implement mass and energy conservation terms calculation. Signed-off-by: Grzegorz Kowal --- sources/statistics.F90 | 658 +++++++++++++++++++++++++++++++++-------- 1 file changed, 540 insertions(+), 118 deletions(-) diff --git a/sources/statistics.F90 b/sources/statistics.F90 index 7e00578..90e4993 100644 --- a/sources/statistics.F90 +++ b/sources/statistics.F90 @@ -261,9 +261,16 @@ module statistics #endif /* __INTEL_COMPILER */ end if - write(cunit,"('#',a8,11(1x,a18))") 'step', 'time', 'dt', 'mass', & - 'momx', 'momy', 'momz', 'ener', & - 'ekin', 'emag', 'eint', 'entr' + 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 @@ -420,10 +427,10 @@ module statistics if (master) then close(munit) - close(cunit) - close(sunit) - if (forcing_enabled) close(funit) - if (error_control) close(eunit) + if (track_conservation) close(cunit) + if (track_statistics) close(sunit) + if (forcing_enabled) close(funit) + if (error_control) close(eunit) end if !------------------------------------------------------------------------------- @@ -443,15 +450,21 @@ module statistics ! subroutine store_statistics() - use blocks , only : block_leaf, block_data + use blocks , only : block_leaf, block_meta, block_data use blocks , only : list_leaf, list_data use blocks , only : get_mblocks, get_nleafs - use coordinates, only : ni => ncells, nb, ne + use coordinates, only : ni => ncells, nn => bcells, nb, ne, nbm, nep use coordinates, only : advol, voli use coordinates, only : toplev - use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp +#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 : ien, imx, imy, imz - use equations , only : magnetized, adiabatic_index, csnd, csnd2 + use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2 use equations , only : errors use evolution , only : error_control use evolution , only : step, time, dt, dth, dte @@ -463,20 +476,37 @@ module statistics use mpitools , only : nprocs use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum #endif /* MPI */ + use operators , only : curl, divergence, gradient, laplace + use sources , only : viscosity, resistivity use workspace , only : resize_workspace, work, work_in_use implicit none logical, save :: first = .true. integer :: n, l, u, nk, nl, nm, status +#if NDIMS == 3 + real(kind=8) :: xl, xu, yl, yu, zl, zu +#else /* NDIMS == 3 */ + real(kind=8) :: xl, xu, yl, yu +#endif /* NDIMS == 3 */ real(kind=8) :: dvol, dvolh +#if NDIMS == 3 + real(kind=8) :: dyz, dxz, dxy +#else /* NDIMS == 3 */ + real(kind=8) :: dyz, dxz +#endif /* NDIMS == 3 */ type(block_leaf), pointer :: pleaf + type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata - real(kind=8), dimension(:,:,:), pointer, save :: vel, mag, sqd, tmp + real(kind=8), dimension(3) :: dh - integer , parameter :: narr = 16 + real(kind=8), dimension(:,:,:,:), pointer, save :: jc + real(kind=8), dimension(:,:,:) , pointer, save :: qq + real(kind=8), dimension(:,:,:) , pointer, save :: vel, mag, sqd, tmp + + integer , parameter :: narr = 32 real(kind=8), dimension(narr) :: inarr, avarr, mnarr, mxarr @@ -490,7 +520,7 @@ module statistics integer, save :: nt !$ integer :: omp_get_thread_num -!$omp threadprivate(first, nt) +!$omp threadprivate(first, nt, jc, qq) character(len=*), parameter :: loc = 'INTEGRALS:store_statistics()' @@ -546,7 +576,7 @@ module statistics if (mod(step, iintd) > 0) return if (first) then - n = 4 * ni**NDIMS + n = 4 * nn**NDIMS + (nv + 1) * ni**(NDIMS - 1) + 4 * ni**NDIMS call resize_workspace(n, status) if (status /= 0) then @@ -560,17 +590,31 @@ module statistics nk = 1 #endif /* NDIMS == 3 */ - n = ni**NDIMS l = 1 - u = n + u = l - 1 + 4 * nn**NDIMS +#if NDIMS == 3 + jc(0:3,1:nn,1:nn,1:nn) => work(l:u,nt) +#else /* NDIMS == 3 */ + jc(0:3,1:nn,1:nn,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ + l = u + 1 + u = l - 1 + (nv + 1) * ni**(NDIMS - 1) +#if NDIMS == 3 + qq(0:nv,1:ni,1:ni) => work(l:u,nt) +#else /* NDIMS == 3 */ + qq(0:nv,1:ni,1: 1) => work(l:u,nt) +#endif /* NDIMS == 3 */ + n = ni**NDIMS + l = u + 1 + u = u + n vel(1:ni,1:ni,1:nk) => work(l:u,nt) - l = l + n + l = u + 1 u = u + n mag(1:ni,1:ni,1:nk) => work(l:u,nt) - l = l + n + l = u + 1 u = u + n sqd(1:ni,1:ni,1:nk) => work(l:u,nt) - l = l + n + l = u + 1 u = u + n tmp(1:ni,1:ni,1:nk) => work(l:u,nt) @@ -612,97 +656,11 @@ module statistics ! iterate over all data blocks on the list ! do while(associated(pdata)) + pmeta => pdata%meta -! obtain the volume elements for the current block -! - dvol = advol(pdata%meta%level) + dvol = advol(pmeta%level) dvolh = 0.5d+00 * dvol -! sum up density and momenta components -! -#if NDIMS == 3 - inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvol - inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne,nb:ne)) * dvol - inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne,nb:ne)) * dvol - inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(1) = inarr(1) + sum(pdata%u(idn,nb:ne,nb:ne, : )) * dvol - inarr(2) = inarr(2) + sum(pdata%u(imx,nb:ne,nb:ne, : )) * dvol - inarr(3) = inarr(3) + sum(pdata%u(imy,nb:ne,nb:ne, : )) * dvol - inarr(4) = inarr(4) + sum(pdata%u(imz,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - -! sum up total energy and entropy -! - if (ien > 0) then -#if NDIMS == 3 - inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(5) = inarr(5) + sum(pdata%u(ien,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - -#if NDIMS == 3 - tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne,nb:ne) & - / pdata%q(idn,nb:ne,nb:ne,nb:ne)**adiabatic_index) & - * pdata%q(idn,nb:ne,nb:ne,nb:ne) / (adiabatic_index - 1.0d+00) -#else /* NDIMS == 3 */ - tmp(:,:,:) = - log(pdata%q(ipr,nb:ne,nb:ne, : ) & - / pdata%q(idn,nb:ne,nb:ne, : )**adiabatic_index) & - * pdata%q(idn,nb:ne,nb:ne, : ) / (adiabatic_index - 1.0d+00) -#endif /* NDIMS == 3 */ - else -#if NDIMS == 3 - tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne,nb:ne)) & - * pdata%q(idn,nb:ne,nb:ne,nb:ne) -#else /* NDIMS == 3 */ - tmp(:,:,:) = csnd2 * log(pdata%q(idn,nb:ne,nb:ne, : )) & - * pdata%q(idn,nb:ne,nb:ne, : ) -#endif /* NDIMS == 3 */ - end if - inarr(9) = inarr(9) + sum(tmp(:,:,:)) * dvol - -! sum up kinetic energy -! -#if NDIMS == 3 - inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne,nb:ne)**2 & - + pdata%u(imy,nb:ne,nb:ne,nb:ne)**2 & - + pdata%u(imz,nb:ne,nb:ne,nb:ne)**2) & - / pdata%u(idn,nb:ne,nb:ne,nb:ne)) * dvolh -#else /* NDIMS == 3 */ - inarr(6) = inarr(6) + sum((pdata%u(imx,nb:ne,nb:ne, : )**2 & - + pdata%u(imy,nb:ne,nb:ne, : )**2 & - + pdata%u(imz,nb:ne,nb:ne, : )**2) & - / pdata%u(idn,nb:ne,nb:ne, : )) * dvolh -#endif /* NDIMS == 3 */ - -! sum up magnetic energy -! - if (magnetized) then -#if NDIMS == 3 - inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne,nb:ne)**2 & - + pdata%u(iby,nb:ne,nb:ne,nb:ne)**2 & - + pdata%u(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh -#else /* NDIMS == 3 */ - inarr(7) = inarr(7) + sum(pdata%u(ibx,nb:ne,nb:ne, : )**2 & - + pdata%u(iby,nb:ne,nb:ne, : )**2 & - + pdata%u(ibz,nb:ne,nb:ne, : )**2) * dvolh -#endif /* NDIMS == 3 */ - if (ien <= 0) then -#if NDIMS == 3 - inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne,nb:ne)**2) * dvolh -#else /* NDIMS == 3 */ - inarr(9) = inarr(9) + sum(pdata%u(ibp,nb:ne,nb:ne, : )**2) * dvolh -#endif /* NDIMS == 3 */ - end if - end if - -! sum up the injected energy and injection rate -! - if (forcing_enabled) then - inarr(11) = einj - inarr(12) = rinj - end if - ! get average, minimum and maximum values of density ! #if NDIMS == 3 @@ -794,8 +752,462 @@ module statistics mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) end if -! associate the pointer with the next block on the list + if (track_conservation) then + + xl = pmeta%xmin - adx(pmeta%level) + xu = pmeta%xmax + adx(pmeta%level) + yl = pmeta%ymin - ady(pmeta%level) + yu = pmeta%ymax + ady(pmeta%level) +#if NDIMS == 3 + zl = pmeta%zmin - adz(pmeta%level) + zu = pmeta%zmax + adz(pmeta%level) +#endif /* NDIMS == 3 */ + + dyz = ady(pmeta%level) * adz(pmeta%level) + dxz = adx(pmeta%level) * adz(pmeta%level) +#if NDIMS == 3 + dxy = adx(pmeta%level) * ady(pmeta%level) +#endif /* NDIMS == 3 */ + + dh(1) = adx(pmeta%level) + dh(2) = ady(pmeta%level) + dh(3) = adz(pmeta%level) + +! calculate current density (J = ∇xB) ! + call curl(dh(:), pdata%q(ibx:ibz,:,:,:), jc(1:3,:,:,:)) + +! total mass +! +#if NDIMS == 3 + inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + +! sum up kinetic energy +! +#if NDIMS == 3 + inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh +#else /* NDIMS == 3 */ + inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne, : )**2 & + + pdata%q(ivy,nb:ne,nb:ne, : )**2 & + + pdata%q(ivz,nb:ne,nb:ne, : )**2) & + * pdata%q(idn,nb:ne,nb:ne, : )) * dvolh +#endif /* NDIMS == 3 */ + +! sum up internal energy +! + if (ipr > 0) then +#if NDIMS == 3 + inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + end if + +! sum up magnetic energy +! + if (magnetized) then +#if NDIMS == 3 + inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 & + + pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh +#else /* NDIMS == 3 */ + inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne, : )**2 & + + pdata%q(iby,nb:ne,nb:ne, : )**2 & + + pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh +#endif /* NDIMS == 3 */ + end if + + if (.not. periodic(1)) then + if (xl < xmin) then +#if NDIMS == 3 + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nbm:nb,nb:ne,nb:ne), 2) +#else /* NDIMS == 3 */ + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nbm:nb,nb:ne, : ), 2) +#endif /* NDIMS == 3 */ + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(2) = inarr(2) + sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz + inarr(6) = inarr(6) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + + if (ipr > 0) then + inarr(10) = inarr(10) + sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(14) = inarr(14) + sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(17) = inarr(17) - sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz + + if (resistivity > 0.0d+00) then +#if NDIMS == 3 + qq(1,:,:) = 5.0d-01 * sum(jc(2,nbm:nb,nb:ne,nb:ne), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nbm:nb,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + qq(2,:,:) = 5.0d-01 * sum(jc(2,nbm:nb,nb:ne, : ), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nbm:nb,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + inarr(20) = inarr(20) & + + sum(qq(2,:,:) * qq(ibz,:,:) & + - qq(3,:,:) * qq(iby,:,:)) * dyz + end if ! resistivity + end if ! magnetized + end if + + if (xu > xmax) then +#if NDIMS == 3 + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,ne:nep,nb:ne,nb:ne), 2) +#else /* NDIMS == 3 */ + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,ne:nep,nb:ne, : ), 2) +#endif /* NDIMS == 3 */ + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(2) = inarr(2) - sum(qq(idn,:,:) * qq(ivx,:,:)) * dyz + inarr(6) = inarr(6) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + + if (ipr > 0) then + inarr(10) = inarr(10) - sum(qq(ipr,:,:) * qq(ivx,:,:)) * dyz + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(14) = inarr(14) - sum(qq( 0 ,:,:) * qq(ivx,:,:)) * dyz + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(17) = inarr(17) + sum(qq( 0 ,:,:) * qq(ibx,:,:)) * dyz + + if (resistivity > 0.0d+00) then +#if NDIMS == 3 + qq(1,:,:) = 5.0d-01 * sum(jc(2,ne:nep,nb:ne,nb:ne), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,ne:nep,nb:ne,nb:ne), 1) +#else /* NDIMS == 3 */ + qq(2,:,:) = 5.0d-01 * sum(jc(2,ne:nep,nb:ne, : ), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,ne:nep,nb:ne, : ), 1) +#endif /* NDIMS == 3 */ + inarr(20) = inarr(20) & + - sum(qq(2,:,:) * qq(ibz,:,:) & + - qq(3,:,:) * qq(iby,:,:)) * dyz + end if ! resistivity + end if ! magnetized + end if + end if + + if (.not. periodic(2)) then + if (yl < ymin) then +#if NDIMS == 3 + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nbm:nb,nb:ne), 2) +#else /* NDIMS == 3 */ + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nbm:nb, : ), 2) +#endif /* NDIMS == 3 */ + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(3) = inarr(3) + sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz + inarr(7) = inarr(7) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + + if (ipr > 0) then + inarr(11) = inarr(11) + sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(15) = inarr(15) + sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(18) = inarr(18) - sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz + + if (resistivity > 0.0d+00) then +#if NDIMS == 3 + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nbm:nb,nb:ne), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb,nb:ne), 1) +#else /* NDIMS == 3 */ + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nbm:nb, : ), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,nbm:nb, : ), 1) +#endif /* NDIMS == 3 */ + inarr(21) = inarr(21) & + + sum(qq(3,:,:) * qq(ibx,:,:) & + - qq(1,:,:) * qq(ibz,:,:)) * dxz + end if ! resistivity + end if ! magnetized + end if + + if (yu > ymax) then +#if NDIMS == 3 + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,ne:nep,nb:ne), 2) +#else /* NDIMS == 3 */ + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,ne:nep, : ), 2) +#endif /* NDIMS == 3 */ + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(3) = inarr(3) - sum(qq(idn,:,:) * qq(ivy,:,:)) * dxz + inarr(7) = inarr(7) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + + if (ipr > 0) then + inarr(11) = inarr(11) - sum(qq(ipr,:,:) * qq(ivy,:,:)) * dxz + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(15) = inarr(15) - sum(qq( 0 ,:,:) * qq(ivy,:,:)) * dxz + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(18) = inarr(18) + sum(qq( 0 ,:,:) * qq(iby,:,:)) * dxz + + if (resistivity > 0.0d+00) then +#if NDIMS == 3 + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,ne:nep,nb:ne), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep,nb:ne), 1) +#else /* NDIMS == 3 */ + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,ne:nep, : ), 1) + qq(3,:,:) = 5.0d-01 * sum(jc(3,nb:ne,ne:nep, : ), 1) +#endif /* NDIMS == 3 */ + inarr(21) = inarr(21) & + - sum(qq(3,:,:) * qq(ibx,:,:) & + - qq(1,:,:) * qq(ibz,:,:)) * dxz + end if ! resistivity + end if ! magnetized + end if + end if + +#if NDIMS == 3 + if (.not. periodic(3)) then + if (zl < zmin) then + + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,nbm:nb), 2) + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(4) = inarr(4) + sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy + inarr(8) = inarr(8) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + + if (ipr > 0) then + inarr(12) = inarr(12) + sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(16) = inarr(16) + sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(19) = inarr(19) - sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy + + if (resistivity > 0.0d+00) then + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nb:ne,nbm:nb), 1) + qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,nbm:nb), 1) + inarr(22) = inarr(22) & + + sum(qq(1,:,:) * qq(iby,:,:) & + - qq(2,:,:) * qq(ibx,:,:)) * dxy + end if ! resistivity + end if ! magnetized + end if + + if (zu > zmax) then + qq(1:nv,:,:) = 5.0d-01 * sum(pdata%q(1:nv,nb:ne,nb:ne,ne:nep), 1) + qq( 0 ,:,:) = 5.0d-01 * qq(idn,:,:) * sum(qq(ivx:ivz,:,:)**2,1) + + inarr(4) = inarr(4) - sum(qq(idn,:,:) * qq(ivz,:,:)) * dxy + inarr(8) = inarr(8) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + + if (ipr > 0) then + inarr(12) = inarr(12) - sum(qq(ipr,:,:) * qq(ivz,:,:)) * dxy + end if + if (magnetized) then + qq( 0 ,:,:) = sum(qq(ibx:ibz,:,:)**2,1) + inarr(16) = inarr(16) - sum(qq( 0 ,:,:) * qq(ivz,:,:)) * dxy + + qq( 0 ,:,:) = sum(qq(ivx:ivz,:,:) * qq(ibx:ibz,:,:),1) + inarr(19) = inarr(19) + sum(qq( 0 ,:,:) * qq(ibz,:,:)) * dxy + + if (resistivity > 0.0d+00) then + qq(1,:,:) = 5.0d-01 * sum(jc(1,nb:ne,nb:ne,ne:nep), 1) + qq(2,:,:) = 5.0d-01 * sum(jc(2,nb:ne,nb:ne,ne:nep), 1) + inarr(22) = inarr(22) & + - sum(qq(1,:,:) * qq(iby,:,:) & + - qq(2,:,:) * qq(ibx,:,:)) * dxy + end if ! resistivity + end if ! magnetized + end if + end if +#endif /* NDIMS == 3 */ + + if (magnetized) then + +! conversion between kinetic and magnetic energies +! + inarr(24) = inarr(24) & +#if NDIMS == 3 + - sum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibz,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivz,nb:ne,nb:ne,nb:ne) & + * pdata%q(iby,nb:ne,nb:ne,nb:ne)) & + * jc(1,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum((pdata%q(ivy,nb:ne,nb:ne, : ) & + * pdata%q(ibz,nb:ne,nb:ne, : ) & + - pdata%q(ivz,nb:ne,nb:ne, : ) & + * pdata%q(iby,nb:ne,nb:ne, : )) & + * jc(1,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + inarr(24) = inarr(24) & +#if NDIMS == 3 + - sum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibx,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivx,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibz,nb:ne,nb:ne,nb:ne)) & + * jc(2,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum((pdata%q(ivz,nb:ne,nb:ne, : ) & + * pdata%q(ibx,nb:ne,nb:ne, : ) & + - pdata%q(ivx,nb:ne,nb:ne, : ) & + * pdata%q(ibz,nb:ne,nb:ne, : )) & + * jc(2,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + inarr(24) = inarr(24) & +#if NDIMS == 3 + - sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) & + * pdata%q(iby,nb:ne,nb:ne,nb:ne) & + - pdata%q(ivy,nb:ne,nb:ne,nb:ne) & + * pdata%q(ibx,nb:ne,nb:ne,nb:ne)) & + * jc(3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum((pdata%q(ivx,nb:ne,nb:ne, : ) & + * pdata%q(iby,nb:ne,nb:ne, : ) & + - pdata%q(ivy,nb:ne,nb:ne, : ) & + * pdata%q(ibx,nb:ne,nb:ne, : )) & + * jc(3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + +! conversion between magnetic and internal energies +! + if (resistivity > 0.0d+00) then + inarr(26) = inarr(26) & +#if NDIMS == 3 + - sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol +#else /* NDIMS == 3 */ + - sum(jc(1:3,nb:ne,nb:ne, : )**2) * dvol +#endif /* NDIMS == 3 */ + end if ! resistive + +! calculate product (v·B) +! + jc(1,:,:,:) = sum(pdata%q(ivx:ivz,:,:,:) * pdata%q(ibx:ibz,:,:,:), 1) + +! calculate divergence (∇·B) +! + call divergence(dh(:), pdata%q(ibx:ibz,:,:,:), jc(2,:,:,:)) + + inarr(27) = inarr(27) & +#if NDIMS == 3 + - sum(jc(1,nb:ne,nb:ne,nb:ne) & + * jc(2,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum(jc(1,nb:ne,nb:ne, : ) & + * jc(2,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + +! calculate gradient (∇ψ) +! + call gradient(dh(:), pdata%q(ibp,:,:,:), jc(1:3,:,:,:)) + + inarr(28) = inarr(28) & +#if NDIMS == 3 + - sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + + end if ! magnetized + +! viscous heating +! + if (viscosity > 0.0d+00) then + +! calculate Laplace (∇²v) +! + call laplace(dh(:), pdata%q(ivx,:,:,:), jc(1,:,:,:)) + call laplace(dh(:), pdata%q(ivy,:,:,:), jc(2,:,:,:)) + call laplace(dh(:), pdata%q(ivz,:,:,:), jc(3,:,:,:)) + + inarr(25) = inarr(25) & +#if NDIMS == 3 + + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne),1) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ),1) & + * pdata%q(idn,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + +! calculate divergence (∇·v) +! + call divergence(dh(:), pdata%q(ivx:ivz,:,:,:), jc(0,:,:,:)) + +! calculate gradient (∇p) +! + call gradient(dh(:), jc(0,:,:,:), jc(1:3,:,:,:)) + + inarr(25) = inarr(25) & +#if NDIMS == 3 + + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne),1) & + * pdata%q(idn,nb:ne,nb:ne,nb:ne)) & + * dvol / 3.0d+00 +#else /* NDIMS == 3 */ + + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : ),1) & + * pdata%q(idn,nb:ne,nb:ne, : )) & + * dvol / 3.0d+00 +#endif /* NDIMS == 3 */ + + end if + +! kinetic energy due to pressure gradient and the change of pressure +! due to the compression +! + if (ipr > 0) then + +! calculate gradient (∇p) +! + call gradient(dh(:), pdata%q(ipr,:,:,:), jc(1:3,:,:,:)) + + inarr(23) = inarr(23) & +#if NDIMS == 3 + - sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + + else + +! calculate gradient (∇ρ) +! + call gradient(dh(:), pdata%q(idn,:,:,:), jc(1:3,:,:,:)) + + inarr(23) = inarr(23) & +#if NDIMS == 3 + - sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) & + * jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol +#else /* NDIMS == 3 */ + - sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) & + * jc(1:3,nb:ne,nb:ne, : )) * dvol +#endif /* NDIMS == 3 */ + + end if + + end if ! track conservation + +! sum up the injected energy and injection rate +! + if (forcing_enabled) then + inarr(29) = einj + inarr(30) = rinj + end if + pdata => pdata%next end do ! data blocks @@ -814,13 +1226,23 @@ module statistics call reduce_maximum(mxarr(1:narr)) #endif /* MPI */ -! calculate the internal energy, or update the entropy -! - if (ien > 0) then - inarr(8) = inarr(5) - inarr(6) - inarr(7) - else - inarr(9) = inarr(9) + inarr(6) + inarr(7) - end if + if (track_conservation) then + if (ipr > 0) then + inarr( 9:12) = inarr( 9:12) / (adiabatic_index - 1.0d+00) + inarr(10:12) = inarr(10:12) * adiabatic_index + else + inarr(23) = csnd2 * inarr(23) + end if + if (viscosity > 0.0d+00) then + inarr(25) = viscosity * inarr(25) + end if ! resistivity + if (magnetized) then + if (resistivity > 0.0d+00) then + inarr(20:22) = resistivity * inarr(20:22) + inarr(26) = resistivity * inarr(26) + end if ! resistivity + end if ! magnetized + end if ! track conservation ! normalize the averages by the volume of domain ! @@ -830,7 +1252,7 @@ module statistics ! if (master) then if (track_conservation) then - write(cunit,"(i9,11(1x,1es18.8e3))") step, time, dt, inarr(1:9) + write(cunit,"(31es25.15e3)") time, inarr(1:30) call flush_and_sync(cunit) end if if (track_statistics) then @@ -846,7 +1268,7 @@ module statistics end if if (forcing_enabled) then - write(funit,"(4(1x,1es18.8e3))") time, inarr(11:12), arms + write(funit,"(4(1x,1es18.8e3))") time, inarr(29:30), arms call flush_and_sync(funit) end if if (error_control) then From 5c8d13f385ee31327df8290136dba9cf605abd7b Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Mon, 24 Oct 2022 08:54:03 -0300 Subject: [PATCH 10/10] STATISTICS: Remove unused variables. Signed-off-by: Grzegorz Kowal --- sources/statistics.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sources/statistics.F90 b/sources/statistics.F90 index 90e4993..e325a96 100644 --- a/sources/statistics.F90 +++ b/sources/statistics.F90 @@ -463,11 +463,10 @@ module statistics #endif /* NDIMS == 3 */ use coordinates, only : adx, ady, adz use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp - use equations , only : ien, imx, imy, imz use equations , only : nv, magnetized, adiabatic_index, csnd, csnd2 use equations , only : errors use evolution , only : error_control - use evolution , only : step, time, dt, dth, dte + use evolution , only : step, time, dth, dte use forcing , only : forcing_enabled, einj, rinj, arms use helpers , only : print_message, flush_and_sync use mesh , only : changed