From a3e0c0ccd2ca5a0b95864a3cfcf797f14fa8f5d5 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 26 Oct 2021 12:29:16 -0300 Subject: [PATCH] INTEGRALS: Remove reconnection rate terms from INTEGRALS. Signed-off-by: Grzegorz Kowal --- sources/integrals.F90 | 639 +-------------------------------------- sources/user_problem.F90 | 4 +- 2 files changed, 7 insertions(+), 636 deletions(-) diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 57b5f42..f9f1bfa 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -57,7 +57,6 @@ module integrals integer(kind=4), save :: funit = 11 integer(kind=4), save :: sunit = 12 integer(kind=4), save :: eunit = 13 - integer(kind=4), save :: runit = 14 integer(kind=4), save :: iintd = 1 ! flag indicating if the files are actually written to disk @@ -294,55 +293,6 @@ module integrals ! efmt = "(i9," // trim(adjustl(stmp)) // "(1x,1es18.8e3))" -! depending on the append parameter, choose the right file initialization for -! the reconnection file -! - append = "off" - call get_parameter("reconnection_append", append) - select case(trim(append)) - case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - write(fname, "('reconnection.dat')") - inquire(file = fname, exist = flag) - case default - write(fname, "('reconnection_',i2.2,'.dat')") irun - flag = .false. - end select - -! check if the file exists; if not, create a new one, otherwise move to the end -! - if (flag .and. irun > 1) then -#ifdef __INTEL_COMPILER - open(newunit = runit, file = fname, form = 'formatted', status = 'old' & - , position = 'append', buffered = 'yes') -#else /* __INTEL_COMPILER */ - open(newunit = runit, file = fname, form = 'formatted', status = 'old' & - , position = 'append') -#endif /* __INTEL_COMPILER */ - write(runit,"('#')") - else -#ifdef __INTEL_COMPILER - open(newunit = runit, file = fname, form = 'formatted' & - , status = 'replace', buffered = 'yes') -#else /* __INTEL_COMPILER */ - open(newunit = runit, file = fname, form = 'formatted' & - , status = 'replace') -#endif /* __INTEL_COMPILER */ - end if - -! write the integral file header -! - write(runit,'("#",a8,27a25)') & - 'step', 'time', '|Bx| int', '|Bx| inf' & - , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' & - , '|Bx| y-dif', '|Bx| z-dif', '|Bx| Psi' & - , 'Vin lower' , 'Vin upper' & - , 'Emag', 'Emag inf' & - , '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' & - , 'Emag-Ekin', 'Emag-Eint', 'Emag-Psi' - write(runit,"('#')") - end if ! store #ifdef PROFILE @@ -396,7 +346,6 @@ module integrals close(funit) close(sunit) close(eunit) - close(runit) end if #ifdef PROFILE @@ -425,16 +374,9 @@ module integrals ! import external variables and subroutines ! use blocks , only : block_meta, block_data, list_data - use coordinates , only : ni => ncells, nn => bcells - use coordinates , only : nb, nbm, nbp, ne, nem, nep - use coordinates , only : adx, ady, adz, advol, voli - use coordinates , only : periodic - use coordinates , only : xmin, xmax - use coordinates , only : ymin, ymax, yarea -#if NDIMS == 3 - use coordinates , only : zmin, zmax -#endif /* NDIMS == 3 */ - use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp + use coordinates , only : ni => ncells, nb, ne + use coordinates , only : advol, voli + use equations , only : idn, ipr, ivx, ivz, ibx, iby, ibz, ibp use equations , only : ien, imx, imy, imz use equations , only : magnetized, adiabatic_index, csnd use equations , only : errors @@ -444,8 +386,6 @@ module integrals #ifdef MPI use mpitools , only : reduce_minimum, reduce_maximum, reduce_sum #endif /* MPI */ - use operators , only : gradient, curl - use sources , only : resistivity ! local variables are not implicit by default ! @@ -453,11 +393,7 @@ module integrals ! local variables ! - real(kind=8) :: dvol, dvolh, dyz, dxz -#if NDIMS == 3 - real(kind=8) :: dxy -#endif /* NDIMS == 3 */ - real(kind=8), dimension(3) :: dh + real(kind=8) :: dvol, dvolh ! local pointers ! @@ -465,7 +401,7 @@ module integrals ! local parameters ! - integer, parameter :: narr = 64 + integer, parameter :: narr = 16 ! local arrays ! @@ -475,11 +411,6 @@ module integrals #else /* NDIMS == 3 */ real(kind=8), dimension(ni,ni, 1) :: vel, mag, sqd, tmp #endif /* NDIMS == 3 */ -#if NDIMS == 3 - real(kind=8), dimension(3,nn,nn,nn) :: cur -#else /* NDIMS == 3 */ - real(kind=8), dimension(3,nn,nn, 1) :: cur -#endif /* NDIMS == 3 */ ! parameters ! @@ -532,14 +463,6 @@ module integrals ! dvol = advol(pdata%meta%level) dvolh = 0.5d+00 * dvol - dyz = ady(pdata%meta%level) * adz(pdata%meta%level) -#if NDIMS == 3 - dxy = adx(pdata%meta%level) * ady(pdata%meta%level) -#endif /* NDIMS == 3 */ - dxz = adx(pdata%meta%level) * adz(pdata%meta%level) - dh(1) = adx(pdata%meta%level) - dh(2) = ady(pdata%meta%level) - dh(3) = adz(pdata%meta%level) ! sum up density and momenta components ! @@ -687,546 +610,6 @@ module integrals avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol mnarr(7) = min(mnarr(7), minval(tmp(:,:,:))) mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:))) - -!! RECONNECTION RATE MEASUREMENTS -!! -! calculate current density (J = ∇xB) -! - call curl(dh(:), pdata%q(ibx:ibz,:,:,:), cur(1:3,:,:,:)) - -! the integral of |Bx| -! -#if NDIMS == 3 - inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol -#else /* NDIMS == 3 */ - inarr(11) = inarr(11) + sum(abs(pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol -#endif /* NDIMS == 3 */ - -! magnetic energy -! - inarr(22) = inarr(7) - -! fluxes across the X boundaries -! - if (.not. periodic(1)) then - -! lower X-boundary -! - if (pdata%meta%xmin <= xmin + dh(1)) then - -! advection of magnetic energy through the lower X boundary -! -#if NDIMS == 3 - inarr(24) = inarr(24) & - + sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne)**2,1) & - * pdata%q(ivx ,nbm:nb,nb:ne,nb:ne)) * dyz - inarr(27) = inarr(27) & - - sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne,nb:ne) & - * pdata%q(ibx:ibz,nbm:nb,nb:ne,nb:ne),1) & - * pdata%q(ibx ,nbm:nb,nb:ne,nb:ne)) * dyz -#else /* NDIMS == 3 */ - inarr(24) = inarr(24) & - + sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne, : )**2,1) & - * pdata%q(ivx ,nbm:nb,nb:ne, : )) * dyz - inarr(27) = inarr(27) & - - sum(sum(pdata%q(ivx:ivz,nbm:nb,nb:ne, : ) & - * pdata%q(ibx:ibz,nbm:nb,nb:ne, : ),1) & - * pdata%q(ibx ,nbm:nb,nb:ne, : )) * dyz -#endif /* NDIMS == 3 */ - - if (resistivity > 0.0d+00) then - -#if NDIMS == 3 -! diffusion of magnetic energy through the lower X boundary -! - inarr(30) = inarr(30) & - + sum(cur(2,nbm:nb,nb:ne,nb:ne) & - * pdata%q(ibz,nbm:nb,nb:ne,nb:ne) & - - cur(3,nbm:nb,nb:ne,nb:ne) & - * pdata%q(iby,nbm:nb,nb:ne,nb:ne)) * dyz -#else /* NDIMS == 3 */ -! diffusion of magnetic energy through the lower X boundary -! - inarr(30) = inarr(30) & - + sum(cur(2,nbm:nb,nb:ne, : ) & - * pdata%q(ibz,nbm:nb,nb:ne, : ) & - - cur(3,nbm:nb,nb:ne, : ) & - * pdata%q(iby,nbm:nb,nb:ne, : )) * dyz -#endif /* NDIMS == 3 */ - - end if ! resistivity - - end if ! xmin - -! upper X-boundary -! - if (pdata%meta%xmax >= xmax - dh(1)) then - -! advection of magnetic energy through the upper X boundary -! -#if NDIMS == 3 - inarr(24) = inarr(24) & - - sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne)**2,1) & - * pdata%q(ivx ,ne:nep,nb:ne,nb:ne)) * dyz - inarr(27) = inarr(27) & - + sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne,nb:ne) & - * pdata%q(ibx:ibz,ne:nep,nb:ne,nb:ne),1) & - * pdata%q(ibx ,ne:nep,nb:ne,nb:ne)) * dyz -#else /* NDIMS == 3 */ - inarr(24) = inarr(24) & - - sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne, : )**2,1) & - * pdata%q(ivx ,ne:nep,nb:ne, : )) * dyz - inarr(27) = inarr(27) & - + sum(sum(pdata%q(ivx:ivz,ne:nep,nb:ne, : ) & - * pdata%q(ibx:ibz,ne:nep,nb:ne, : ),1) & - * pdata%q(ibx ,ne:nep,nb:ne, : )) * dyz -#endif /* NDIMS == 3 */ - - if (resistivity > 0.0d+00) then - -#if NDIMS == 3 -! diffusion of magnetic energy through the upper X boundary -! - - inarr(30) = inarr(30) & - - sum(cur(2,ne:nep,nb:ne,nb:ne) & - * pdata%q(ibz,ne:nep,nb:ne,nb:ne) & - - cur(3,ne:nep,nb:ne,nb:ne) & - * pdata%q(iby,ne:nep,nb:ne,nb:ne)) * dyz -#else /* NDIMS == 3 */ -! diffusion of magnetic energy through the upper X boundary -! - inarr(30) = inarr(30) & - - sum(cur(2,ne:nep,nb:ne, : ) & - * pdata%q(ibz,ne:nep,nb:ne, : ) & - - cur(3,ne:nep,nb:ne, : ) & - * pdata%q(iby,ne:nep,nb:ne, : )) * dyz -#endif /* NDIMS == 3 */ - - end if ! resistivity - - end if ! xmax - - end if ! not periodic in X - -! fluxes across the Y boundaries -! -! lower Y-boundary -! - if (pdata%meta%ymin <= ymin + dh(2)) then - -! mean Bx at boundary -! -#if NDIMS == 3 - inarr(12) = inarr(12) & - + sum(abs(pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz -#else /* NDIMS == 3 */ - inarr(12) = inarr(12) & - + sum(abs(pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz -#endif /* NDIMS == 3 */ - -! advection of Bx along Y -! -#if NDIMS == 3 - inarr(13) = inarr(13) & - + sum(abs(pdata%q(ibx,nb:ne,nbm:nb,nb:ne)) & - * pdata%q(ivy,nb:ne,nbm:nb,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(13) = inarr(13) & - + sum(abs(pdata%q(ibx,nb:ne,nbm:nb, : )) & - * pdata%q(ivy,nb:ne,nbm:nb, : )) * dxz -#endif /* NDIMS == 3 */ - -! shear of By along X -! -#if NDIMS == 3 - inarr(15) = inarr(15) & - - sum(sign(pdata%q(iby,nb:ne,nbm:nb,nb:ne) & - * pdata%q(ivx,nb:ne,nbm:nb,nb:ne), & - pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz -#else /* NDIMS == 3 */ - inarr(15) = inarr(15) & - - sum(sign(pdata%q(iby,nb:ne,nbm:nb, : ) & - * pdata%q(ivx,nb:ne,nbm:nb, : ), & - pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz -#endif /* NDIMS == 3 */ - -! mean magnetic energy at the lower Y boundary -! -#if NDIMS == 3 - inarr(23) = inarr(23) & - + sum(pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne)**2) * dxz -#else /* NDIMS == 3 */ - inarr(23) = inarr(23) & - + sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2) * dxz -#endif /* NDIMS == 3 */ - -! advection of magnetic energy through the lower Y boundary -! -#if NDIMS == 3 - inarr(25) = inarr(25) & - + sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne)**2,1) & - * pdata%q(ivy ,nb:ne,nbm:nb,nb:ne)) * dxz - inarr(28) = inarr(28) & - - sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb,nb:ne) & - * pdata%q(ibx:ibz,nb:ne,nbm:nb,nb:ne),1) & - * pdata%q(iby ,nb:ne,nbm:nb,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(25) = inarr(25) & - + sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2,1) & - * pdata%q(ivy ,nb:ne,nbm:nb, : )) * dxz - inarr(28) = inarr(28) & - - sum(sum(pdata%q(ivx:ivz,nb:ne,nbm:nb, : ) & - * pdata%q(ibx:ibz,nb:ne,nbm:nb, : ),1) & - * pdata%q(iby ,nb:ne,nbm:nb, : )) * dxz -#endif /* NDIMS == 3 */ - - if (resistivity > 0.0d+00) then - -#if NDIMS == 3 -! diffusion of Bx through the lower Y boundary -! - inarr(17) = inarr(17) & - + sum(sign(cur(3,nb:ne,nbm:nb,nb:ne), & - pdata%q(ibx,nb:ne,nbm:nb,nb:ne))) * dxz - -! diffusion of magnetic energy through the lower Y boundary -! - inarr(31) = inarr(31) & - + sum(cur(3,nb:ne,nbm:nb,nb:ne) & - * pdata%q(ibx,nb:ne,nbm:nb,nb:ne) & - - cur(1,nb:ne,nbm:nb,nb:ne) & - * pdata%q(ibz,nb:ne,nbm:nb,nb:ne)) * dxz -#else /* NDIMS == 3 */ -! diffusion of Bx through the lower Y boundary -! - inarr(17) = inarr(17) & - + sum(sign(cur(3,nb:ne,nbm:nb, : ), & - pdata%q(ibx,nb:ne,nbm:nb, : ))) * dxz - -! diffusion of magnetic energy through the lower Y boundary -! - inarr(31) = inarr(31) & - + sum(cur(3,nb:ne,nbm:nb, : ) & - * pdata%q(ibx,nb:ne,nbm:nb, : ) & - - cur(1,nb:ne,nbm:nb, : ) & - * pdata%q(ibz,nb:ne,nbm:nb, : )) * dxz -#endif /* NDIMS == 3 */ - - end if ! resistivity - - end if ! ymin - -! upper Y-boundary -! - if (pdata%meta%ymax >= ymax - dh(2)) then - -! mean Bx at boundary -! -#if NDIMS == 3 - inarr(12) = inarr(12) & - + sum(abs(pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz -#else /* NDIMS == 3 */ - inarr(12) = inarr(12) & - + sum(abs(pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz -#endif /* NDIMS == 3 */ - -! advection of Bx along Y -! -#if NDIMS == 3 - inarr(13) = inarr(13) & - - sum(abs(pdata%q(ibx,nb:ne,ne:nep,nb:ne)) & - * pdata%q(ivy,nb:ne,ne:nep,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(13) = inarr(13) & - - sum(abs(pdata%q(ibx,nb:ne,ne:nep, : )) & - * pdata%q(ivy,nb:ne,ne:nep, : )) * dxz -#endif /* NDIMS == 3 */ - -! shear of By along X -! -#if NDIMS == 3 - inarr(15) = inarr(15) & - + sum(sign(pdata%q(iby,nb:ne,ne:nep,nb:ne) & - * pdata%q(ivx,nb:ne,ne:nep,nb:ne), & - pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz -#else /* NDIMS == 3 */ - inarr(15) = inarr(15) & - + sum(sign(pdata%q(iby,nb:ne,ne:nep, : ) & - * pdata%q(ivx,nb:ne,ne:nep, : ), & - pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz -#endif /* NDIMS == 3 */ - -! mean magnetic energy at the upper Y boundary -! -#if NDIMS == 3 - inarr(23) = inarr(23) & - + sum(pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne)**2) * dxz -#else /* NDIMS == 3 */ - inarr(23) = inarr(23) & - + sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2) * dxz -#endif /* NDIMS == 3 */ - -! advection of magnetic energy through the upper Y boundary -! -#if NDIMS == 3 - inarr(25) = inarr(25) & - - sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne)**2,1) & - * pdata%q(ivy ,nb:ne,ne:nep,nb:ne)) * dxz - inarr(28) = inarr(28) & - + sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep,nb:ne) & - * pdata%q(ibx:ibz,nb:ne,ne:nep,nb:ne),1) & - * pdata%q(iby ,nb:ne,ne:nep,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(25) = inarr(25) & - - sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2,1) & - * pdata%q(ivy ,nb:ne,ne:nep, : )) * dxz - inarr(28) = inarr(28) & - + sum(sum(pdata%q(ivx:ivz,nb:ne,ne:nep, : ) & - * pdata%q(ibx:ibz,nb:ne,ne:nep, : ),1) & - * pdata%q(iby ,nb:ne,ne:nep, : )) * dxz -#endif /* NDIMS == 3 */ - - if (resistivity > 0.0d+00) then - -#if NDIMS == 3 -! diffusion of Bx through the upper Y boundary -! - inarr(17) = inarr(17) & - - sum(sign(cur(3,nb:ne,ne:nep,nb:ne), & - pdata%q(ibx,nb:ne,ne:nep,nb:ne))) * dxz - -! diffusion of magnetic energy through the upper Y boundary -! - inarr(31) = inarr(31) & - - sum(cur(3,nb:ne,ne:nep,nb:ne) & - * pdata%q(ibx,nb:ne,ne:nep,nb:ne) & - - cur(1,nb:ne,ne:nep,nb:ne) & - * pdata%q(ibz,nb:ne,ne:nep,nb:ne)) * dxz -#else /* NDIMS == 3 */ -! diffusion of Bx through the upper Y boundary -! - inarr(17) = inarr(17) & - - sum(sign(cur(3,nb:ne,ne:nep, : ), & - pdata%q(ibx,nb:ne,ne:nep, : ))) * dxz - -! diffusion of magnetic energy through the upper Y boundary -! - inarr(31) = inarr(31) & - - sum(cur(3,nb:ne,ne:nep, : ) & - * pdata%q(ibx,nb:ne,ne:nep, : ) & - - cur(1,nb:ne,ne:nep, : ) & - * pdata%q(ibz,nb:ne,ne:nep, : )) * dxz -#endif /* NDIMS == 3 */ - - end if ! resistivity - - end if ! ymax - -#if NDIMS == 3 -! fluxes across the Z boundaries -! - if (.not. periodic(3)) then - -! lower Z-boundary -! - if (pdata%meta%zmin <= zmin + dh(3)) then - -! advection of Bx along Z -! - inarr(14) = inarr(14) & - + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbm:nb)) & - * pdata%q(ivz,nb:ne,nb:ne,nbm:nb)) * dxy - -! shear of Bz along X -! - inarr(16) = inarr(16) & - - sum(sign(pdata%q(ibz,nb:ne,nb:ne,nbm:nb) & - * pdata%q(ivx,nb:ne,nb:ne,nbm:nb), & - pdata%q(ibx,nb:ne,nb:ne,nbm:nb))) * dxy - -! advection of magnetic energy through the lower Z boundary -! - inarr(26) = inarr(26) & - + sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb)**2,1) & - * pdata%q(ivz ,nb:ne,nb:ne,nbm:nb)) * dxy - inarr(29) = inarr(29) & - - sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nbm:nb) & - * pdata%q(ibx:ibz,nb:ne,nb:ne,nbm:nb),1) & - * pdata%q(ibz ,nb:ne,nb:ne,nbm:nb)) * dxy - - if (resistivity > 0.0d+00) then - -! diffusion of Bx through the lower Z boundary -! - inarr(18) = inarr(18) & - + sum(sign(cur(2,nb:ne,nb:ne,nbm:nb), & - pdata%q(ibx,nb:ne,nb:ne,nbm:nb))) * dxy - -! diffusion of magnetic energy through the lower Z boundary -! - inarr(32) = inarr(32) & - + sum(cur(1,nb:ne,nb:ne,nbm:nb) & - * pdata%q(iby,nb:ne,nb:ne,nbm:nb) & - - cur(2,nb:ne,nb:ne,nbm:nb) & - * pdata%q(ibx,nb:ne,nb:ne,nbm:nb)) * dxy - - end if ! resistivity - - end if ! zmin - -! upper Z-boundary -! - if (pdata%meta%zmax >= zmax - dh(3)) then - -! advection of Bx along Z -! - inarr(14) = inarr(14) & - - sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:nep)) & - * pdata%q(ivz,nb:ne,nb:ne,ne:nep)) * dxy - -! shear of Bz along X -! - inarr(16) = inarr(16) & - + sum(sign(pdata%q(ibz,nb:ne,nb:ne,ne:nep) & - * pdata%q(ivx,nb:ne,nb:ne,ne:nep), & - pdata%q(ibx,nb:ne,nb:ne,ne:nep))) * dxy - -! advection of magnetic energy through the upper Z boundary -! - inarr(26) = inarr(26) & - - sum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep)**2,1) & - * pdata%q(ivz ,nb:ne,nb:ne,ne:nep)) * dxy - inarr(29) = inarr(29) & - + sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,ne:nep) & - * pdata%q(ibx:ibz,nb:ne,nb:ne,ne:nep),1) & - * pdata%q(ibz ,nb:ne,nb:ne,ne:nep)) * dxy - - if (resistivity > 0.0d+00) then - -! diffusion of Bx through the upper Z boundary -! - inarr(18) = inarr(18) & - - sum(sign(cur(2,nb:ne,nb:ne,ne:nep), & - pdata%q(ibx,nb:ne,nb:ne,ne:nep))) * dxy - -! diffusion of magnetic energy through the upper Z boundary -! - inarr(32) = inarr(32) & - - sum(cur(1,nb:ne,nb:ne,ne:nep) & - * pdata%q(iby,nb:ne,nb:ne,ne:nep) & - - cur(2,nb:ne,nb:ne,ne:nep) & - * pdata%q(ibx,nb:ne,nb:ne,ne:nep)) * dxy - - end if ! resistivity - - end if ! zmax - - end if ! not periodic in Z -#endif /* NDIMS == 3 */ - -! get the conversion between the magnetic energy and kinetic and -! internal energies -! -#if NDIMS == 3 - inarr(33) = inarr(33) & - + 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)) & - * cur(1,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(33) = inarr(33) & - + 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, : )) & - * cur(1,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - -#if NDIMS == 3 - inarr(33) = inarr(33) & - + 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)) & - * cur(2,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(33) = inarr(33) & - + 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, : )) & - * cur(2,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - -#if NDIMS == 3 - inarr(33) = inarr(33) & - + 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)) & - * cur(3,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(33) = inarr(33) & - + 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, : )) & - * cur(3,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - - if (resistivity > 0.0d+00) then -#if NDIMS == 3 - inarr(34) = inarr(34) - sum(cur(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol -#else /* NDIMS == 3 */ - inarr(34) = inarr(34) - sum(cur(1:3,nb:ne,nb:ne, : )**2) * dvol -#endif /* NDIMS == 3 */ - end if - -! calculate gradient (∇ψ) -! - call gradient(dh(:), pdata%q(ibp,:,:,:), cur(1:3,:,:,:)) - -#if NDIMS == 3 - inarr(35) = inarr(35) & - - sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) & - * cur(1:3,nb:ne,nb:ne,nb:ne)) * dvol -#else /* NDIMS == 3 */ - inarr(35) = inarr(35) & - - sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) & - * cur(1:3,nb:ne,nb:ne, : )) * dvol -#endif /* NDIMS == 3 */ - -! contribution to |Bx| from ∇ψ -! -#if NDIMS == 3 - inarr(19) = inarr(19) & - - sum(sign(cur(1,nb:ne,nb:ne,nb:ne), & - pdata%q(ibx,nb:ne,nb:ne,nb:ne))) * dvol -#else /* NDIMS == 3 */ - inarr(19) = inarr(19) & - - sum(sign(cur(1,nb:ne,nb:ne, : ), & - pdata%q(ibx,nb:ne,nb:ne, : ))) * dvol -#endif /* NDIMS == 3 */ - - end if ! magnetized - -! get the inflow speed -! - if (pdata%meta%ymin <= ymin + dh(2)) then -#if NDIMS == 3 - inarr(20) = inarr(20) + sum(pdata%q(ivy,nb:ne,nb,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(20) = inarr(20) + sum(pdata%q(ivy,nb:ne,nb, : )) * dxz -#endif /* NDIMS == 3 */ - end if - if (pdata%meta%ymax >= ymax - dh(2)) then -#if NDIMS == 3 - inarr(21) = inarr(21) - sum(pdata%q(ivy,nb:ne,ne,nb:ne)) * dxz -#else /* NDIMS == 3 */ - inarr(21) = inarr(21) - sum(pdata%q(ivy,nb:ne,ne, : )) * dxz -#endif /* NDIMS == 3 */ end if ! associate the pointer with the next block on the list @@ -1255,16 +638,6 @@ module integrals ! avarr(:) = avarr(:) * voli -! apply factors to the reconnection rate terms -! - inarr(12) = 2.50d-01 * inarr(12) / yarea - inarr(13:16) = 5.00d-01 * inarr(13:16) - inarr(17:18) = 5.00d-01 * resistivity * inarr(17:18) - inarr(23) = 1.25d-01 * inarr(23) / yarea - inarr(24:29) = 5.00d-01 * inarr(24:29) - inarr(30:32) = 5.00d-01 * resistivity * inarr(30:32) - inarr(34) = resistivity * inarr(34) - ! write down the integrals and statistics to appropriate files ! if (stored) then @@ -1278,12 +651,10 @@ module integrals , avarr(6), mnarr(6), mxarr(6) & , avarr(7), mnarr(7), mxarr(7) write(eunit,efmt) step, time, dth, dte, maxval(errors(:)), errors(:) - write(runit,"(i9, 26es25.15e3)") step, time, inarr(11:35) call flush_and_sync(funit) call flush_and_sync(sunit) call flush_and_sync(eunit) - call flush_and_sync(runit) end if #ifdef PROFILE diff --git a/sources/user_problem.F90 b/sources/user_problem.F90 index 00bc3fd..6bfbbe6 100644 --- a/sources/user_problem.F90 +++ b/sources/user_problem.F90 @@ -372,10 +372,10 @@ module user_problem call get_parameter("reconnection_append", append) select case(trim(append)) case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES") - write(fname, "('reconnection-new.dat')") + write(fname, "('reconnection.dat')") inquire(file = fname, exist = flag) case default - write(fname, "('reconnection-new_',i2.2,'.dat')") rcount + write(fname, "('reconnection_',i2.2,'.dat')") rcount flag = .false. end select