From 4af2db266c2e6f9de24911d7590bd0069bf6568e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Tue, 29 Dec 2020 17:47:57 -0300 Subject: [PATCH] INTEGRALS: Add remaining terms of the magnetic energy evolution. Signed-off-by: Grzegorz Kowal --- sources/integrals.F90 | 470 ++++++++++++++++++++++++------------------ 1 file changed, 275 insertions(+), 195 deletions(-) diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 3e6bcf3..8ba8cc8 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -331,13 +331,15 @@ module integrals ! write the integral file header ! - write(runit,'("#",a8,20a25)') & + write(runit,'("#",a8,25a25)') & 'step', 'time', '|Bx| int', '|Bx| inf' & , '|Bx| y-adv', '|Bx| z-adv', '|Bx| y-shr', '|Bx| z-shr' & , '|Bx| y-dif', '|Bx| z-dif', 'Vin lower' , 'Vin upper' & , 'Emag', 'Emag inf' & , 'Emag x-adv', 'Emag y-adv', 'Emag z-adv' & - , 'Emag x-dif', 'Emag y-dif', 'Emag z-dif' + , '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' write(runit,"('#')") end if ! store @@ -423,7 +425,7 @@ module integrals ! use blocks , only : block_meta, block_data, list_data use coordinates , only : ni => ncells - use coordinates , only : nb, nbl, nbu, ne, nel, neu + 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 @@ -461,7 +463,7 @@ module integrals ! local parameters ! - integer, parameter :: narr = 32 + integer, parameter :: narr = 64 ! local arrays ! @@ -704,23 +706,21 @@ module integrals ! advection of magnetic energy through the lower X boundary ! #if NDIMS == 3 - inarr(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne,nb:ne)**2 & - + pdata%q(ibz,nbl:nb,nb:ne,nb:ne)**2) & - * pdata%q(ivx,nbl:nb,nb:ne,nb:ne) & - - (pdata%q(ivy,nbl:nb,nb:ne,nb:ne) & - * pdata%q(iby,nbl:nb,nb:ne,nb:ne) & - + pdata%q(ivz,nbl:nb,nb:ne,nb:ne) & - * pdata%q(ibz,nbl:nb,nb:ne,nb:ne)) & - * pdata%q(ibx,nbl:nb,nb:ne,nb:ne)) * dyz + inarr(23) = inarr(23) & + + 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(26) = inarr(26) & + - 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(23) = inarr(23) + sum((pdata%q(iby,nbl:nb,nb:ne, : )**2 & - + pdata%q(ibz,nbl:nb,nb:ne, : )**2) & - * pdata%q(ivx,nbl:nb,nb:ne, : ) & - - (pdata%q(ivy,nbl:nb,nb:ne, : ) & - * pdata%q(iby,nbl:nb,nb:ne, : ) & - + pdata%q(ivz,nbl:nb,nb:ne, : ) & - * pdata%q(ibz,nbl:nb,nb:ne, : )) & - * pdata%q(ibx,nbl:nb,nb:ne, : )) * dyz + inarr(23) = inarr(23) & + + sum(sum(pdata%q(ibx:ibz,nbm:nb,nb:ne, : )**2,1) & + * pdata%q(ivx ,nbm:nb,nb:ne, : )) * dyz + inarr(26) = inarr(26) & + - 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 @@ -728,40 +728,40 @@ module integrals #if NDIMS == 3 ! the Y component of current density ! - tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,nb ,nb:ne,nbu:neu) & - - pdata%q(ibx,nb ,nb:ne,nbl:nel)) / dh(3) & + tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,nb ,nb:ne,nbp:nep) & + - pdata%q(ibx,nb ,nb:ne,nbm:nem)) / dh(3) & - (pdata%q(ibz,nb ,nb:ne,nb :ne ) & - - pdata%q(ibz,nbl,nb:ne,nb :ne )) / dh(1) + - pdata%q(ibz,nbm,nb:ne,nb :ne )) / dh(1) ! the Z component of current density ! tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne ,nb:ne) & - - pdata%q(iby,nbl,nb :ne ,nb:ne)) / dh(1) & - - 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu,nb:ne) & - - pdata%q(ibx,nb ,nbl:nel,nb:ne)) / dh(2) + - pdata%q(iby,nbm,nb :ne ,nb:ne)) / dh(1) & + - 0.5d+00 * (pdata%q(ibx,nb ,nbp:nep,nb:ne) & + - pdata%q(ibx,nb ,nbm:nem,nb:ne)) / dh(2) ! diffusion of magnetic energy through the lower X boundary ! - inarr(26) = inarr(26) & - - sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne,nb:ne) & + inarr(29) = inarr(29) & + + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne,nb:ne) & - tmp(2,:,:) * pdata%q(iby,nb,nb:ne,nb:ne)) * dyz #else /* NDIMS == 3 */ ! the Y component of current density ! - tmp(1,:,:) = (pdata%q(ibz,nbl,nb:ne, : ) & - - pdata%q(ibz,nb ,nb:ne, : )) / dh(1) + tmp(1,:,:) = - (pdata%q(ibz,nb ,nb:ne, : ) & + - pdata%q(ibz,nbm,nb:ne, : )) / dh(1) ! the Z component of current density ! tmp(2,:,:) = (pdata%q(iby,nb ,nb :ne , : ) & - - pdata%q(iby,nbl,nb :ne , : )) / dh(1) & - - 0.5d+00 * (pdata%q(ibx,nb ,nbu:neu, : ) & - - pdata%q(ibx,nb ,nbl:nel, : )) / dh(2) + - pdata%q(iby,nbm,nb :ne , : )) / dh(1) & + - 0.5d+00 * (pdata%q(ibx,nb ,nbp:nep, : ) & + - pdata%q(ibx,nb ,nbm:nem, : )) / dh(2) ! diffusion of magnetic energy through the lower X boundary ! - inarr(26) = inarr(26) & - - sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) & + inarr(29) = inarr(29) & + + sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) & - tmp(2,:,:) * pdata%q(iby,nb,nb:ne, : )) * dyz #endif /* NDIMS == 3 */ @@ -776,23 +776,21 @@ module integrals ! advection of magnetic energy through the upper X boundary ! #if NDIMS == 3 - inarr(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne,nb:ne)**2 & - + pdata%q(ibz,ne:neu,nb:ne,nb:ne)**2) & - * pdata%q(ivx,ne:neu,nb:ne,nb:ne) & - - (pdata%q(ivy,ne:neu,nb:ne,nb:ne) & - * pdata%q(iby,ne:neu,nb:ne,nb:ne) & - + pdata%q(ivz,ne:neu,nb:ne,nb:ne) & - * pdata%q(ibz,ne:neu,nb:ne,nb:ne)) & - * pdata%q(ibx,ne:neu,nb:ne,nb:ne)) * dyz + inarr(23) = inarr(23) & + - 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(26) = inarr(26) & + + 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(23) = inarr(23) + sum((pdata%q(iby,ne:neu,nb:ne, : )**2 & - + pdata%q(ibz,ne:neu,nb:ne, : )**2) & - * pdata%q(ivx,ne:neu,nb:ne, : ) & - - (pdata%q(ivy,ne:neu,nb:ne, : ) & - * pdata%q(iby,ne:neu,nb:ne, : ) & - + pdata%q(ivz,ne:neu,nb:ne, : ) & - * pdata%q(ibz,ne:neu,nb:ne, : )) & - * pdata%q(ibx,ne:neu,nb:ne, : )) * dyz + inarr(23) = inarr(23) & + - sum(sum(pdata%q(ibx:ibz,ne:nep,nb:ne, : )**2,1) & + * pdata%q(ivx ,ne:nep,nb:ne, : )) * dyz + inarr(26) = inarr(26) & + + 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 @@ -800,40 +798,40 @@ module integrals #if NDIMS == 3 ! the Y component of current density ! - tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,ne ,nb:ne,nbu:neu) & - - pdata%q(ibx,ne ,nb:ne,nbl:nel)) / dh(3) & - - (pdata%q(ibz,neu,nb:ne,nb :ne ) & - - pdata%q(ibz,ne ,nb:ne,nb :ne )) / dh(1) + tmp(1,:,:) = 0.5d+00 * (pdata%q(ibx,ne ,nb:ne,nbp:nep) & + - pdata%q(ibx,ne ,nb:ne,nbm:nem)) / dh(3) & + - (pdata%q(ibz,nep,nb:ne,nb :ne ) & + - pdata%q(ibz,ne ,nb:ne,nb :ne )) / dh(1) ! the Z component of current density ! - tmp(2,:,:) = (pdata%q(iby,neu,nb :ne ,nb:ne) & + tmp(2,:,:) = (pdata%q(iby,nep,nb :ne ,nb:ne) & - pdata%q(iby,ne ,nb :ne ,nb:ne)) / dh(1) & - - 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu,nb:ne) & - - pdata%q(ibx,ne ,nbl:nel,nb:ne)) / dh(2) + - 0.5d+00 * (pdata%q(ibx,ne ,nbp:nep,nb:ne) & + - pdata%q(ibx,ne ,nbm:nem,nb:ne)) / dh(2) ! diffusion of magnetic energy through the upper X boundary ! - inarr(26) = inarr(26) & - + sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne,nb:ne) & + inarr(29) = inarr(29) & + - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne,nb:ne) & - tmp(2,:,:) * pdata%q(iby,ne,nb:ne,nb:ne)) * dyz #else /* NDIMS == 3 */ ! the Y component of current density ! - tmp(1,:,:) = (pdata%q(ibz,neu,nb:ne, : ) & + tmp(1,:,:) = (pdata%q(ibz,nep,nb:ne, : ) & - pdata%q(ibz,ne ,nb:ne, : )) / dh(1) ! the Z component of current density ! - tmp(2,:,:) = (pdata%q(iby,neu,nb :ne , : ) & + tmp(2,:,:) = (pdata%q(iby,nep,nb :ne , : ) & - pdata%q(iby,ne ,nb :ne , : )) / dh(1) & - - 0.5d+00 * (pdata%q(ibx,ne ,nbu:neu, : ) & - - pdata%q(ibx,ne ,nbl:nel, : )) / dh(2) + - 0.5d+00 * (pdata%q(ibx,ne ,nbp:nep, : ) & + - pdata%q(ibx,ne ,nbm:nem, : )) / dh(2) ! diffusion of magnetic energy through the upper X boundary ! - inarr(26) = inarr(26) & - + sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) & + inarr(29) = inarr(29) & + - sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) & - tmp(2,:,:) * pdata%q(iby,ne,nb:ne, : )) * dyz #endif /* NDIMS == 3 */ @@ -860,25 +858,25 @@ module integrals ! advection of Bx along Y ! #if NDIMS == 3 - inarr(13) = inarr(13) + sum(abs(pdata%q(ibx,nb:ne,nbl:nb,nb:ne)) & - * pdata%q(ivy,nb:ne,nbl:nb,nb:ne)) * dxz + 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,nbl:nb, : )) & - * pdata%q(ivy,nb:ne,nbl:nb, : )) * dxz + 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,nbl:nb,nb:ne) & - * pdata%q(ivx,nb:ne,nbl:nb,nb:ne), & - pdata%q(ibx,nb:ne,nbl:nb,nb:ne))) * dxz + - 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,nbl:nb, : ) & - * pdata%q(ivx,nb:ne,nbl:nb, : ), & - pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz + - 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 @@ -892,23 +890,21 @@ module integrals ! advection of magnetic energy through the lower Y boundary ! #if NDIMS == 3 - inarr(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb,nb:ne)**2 & - + pdata%q(ibz,nb:ne,nbl:nb,nb:ne)**2) & - * pdata%q(ivy,nb:ne,nbl:nb,nb:ne) & - - (pdata%q(ivx,nb:ne,nbl:nb,nb:ne) & - * pdata%q(ibx,nb:ne,nbl:nb,nb:ne) & - + pdata%q(ivz,nb:ne,nbl:nb,nb:ne) & - * pdata%q(ibz,nb:ne,nbl:nb,nb:ne)) & - * pdata%q(iby,nb:ne,nbl:nb,nb:ne)) * dxz + inarr(24) = inarr(24) & + + 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(27) = inarr(27) & + - 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(24) = inarr(24) + sum((pdata%q(ibx,nb:ne,nbl:nb, : )**2 & - + pdata%q(ibz,nb:ne,nbl:nb, : )**2) & - * pdata%q(ivy,nb:ne,nbl:nb, : ) & - - (pdata%q(ivx,nb:ne,nbl:nb, : ) & - * pdata%q(ibx,nb:ne,nbl:nb, : ) & - + pdata%q(ivz,nb:ne,nbl:nb, : ) & - * pdata%q(ibz,nb:ne,nbl:nb, : )) & - * pdata%q(iby,nb:ne,nbl:nb, : )) * dxz + inarr(24) = inarr(24) & + + sum(sum(pdata%q(ibx:ibz,nb:ne,nbm:nb, : )**2,1) & + * pdata%q(ivy ,nb:ne,nbm:nb, : )) * dxz + inarr(27) = inarr(27) & + - 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 @@ -916,17 +912,17 @@ module integrals #if NDIMS == 3 ! the Z component of current density ! - tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb ,nb:ne) & - - pdata%q(iby,nbl:nel,nb ,nb:ne)) / dh(1) & + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,nb ,nb:ne) & + - pdata%q(iby,nbm:nem,nb ,nb:ne)) / dh(1) & - (pdata%q(ibx,nb :ne ,nb ,nb:ne) & - - pdata%q(ibx,nb :ne ,nbl,nb:ne)) / dh(2) + - pdata%q(ibx,nb :ne ,nbm,nb:ne)) / dh(2) ! the X component of current density ! tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb ,nb :ne ) & - - pdata%q(ibz,nb:ne,nbl,nb :ne )) / dh(2) & - - 0.5d+00 * (pdata%q(iby,nb:ne,nb ,nbu:neu) & - - pdata%q(iby,nb:ne,nb ,nbl:nel)) / dh(3) + - pdata%q(ibz,nb:ne,nbm,nb :ne )) / dh(2) & + - 0.5d+00 * (pdata%q(iby,nb:ne,nb ,nbp:nep) & + - pdata%q(iby,nb:ne,nb ,nbm:nem)) / dh(3) ! diffusion of Bx through the lower Y boundary ! @@ -935,21 +931,21 @@ module integrals ! diffusion of magnetic energy through the lower Y boundary ! - inarr(27) = inarr(27) & - - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb,nb:ne) & + inarr(30) = inarr(30) & + + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb,nb:ne) & - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb,nb:ne)) * dxz #else /* NDIMS == 3 */ ! the Z component of current density ! - tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) & - - pdata%q(iby,nbl:nel,nb , : )) / dh(1) & + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,nb , : ) & + - pdata%q(iby,nbm:nem,nb , : )) / dh(1) & - (pdata%q(ibx,nb :ne ,nb , : ) & - - pdata%q(ibx,nb :ne ,nbl, : )) / dh(2) + - pdata%q(ibx,nb :ne ,nbm, : )) / dh(2) ! the X component of current density ! tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb , : ) & - - pdata%q(ibz,nb:ne,nbl, : )) / dh(2) + - pdata%q(ibz,nb:ne,nbm, : )) / dh(2) ! diffusion of Bx through the lower Y boundary ! @@ -958,8 +954,8 @@ module integrals ! diffusion of magnetic energy through the lower Y boundary ! - inarr(27) = inarr(27) & - - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) & + inarr(30) = inarr(30) & + + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) & - tmp(:,2,:) * pdata%q(ibz,nb:ne,nb, : )) * dxz #endif /* NDIMS == 3 */ @@ -982,25 +978,25 @@ module integrals ! advection of Bx along Y ! #if NDIMS == 3 - inarr(13) = inarr(13) - sum(abs(pdata%q(ibx,nb:ne,ne:neu,nb:ne)) & - * pdata%q(ivy,nb:ne,ne:neu,nb:ne)) * dxz + 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:neu, : )) & - * pdata%q(ivy,nb:ne,ne:neu, : )) * dxz + 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:neu,nb:ne) & - * pdata%q(ivx,nb:ne,ne:neu,nb:ne), & - pdata%q(ibx,nb:ne,ne:neu,nb:ne))) * dxz + + 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:neu, : ) & - * pdata%q(ivx,nb:ne,ne:neu, : ), & - pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz + + 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 @@ -1014,23 +1010,21 @@ module integrals ! advection of magnetic energy through the upper Y boundary ! #if NDIMS == 3 - inarr(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu,nb:ne)**2 & - + pdata%q(ibz,nb:ne,ne:neu,nb:ne)**2) & - * pdata%q(ivy,nb:ne,ne:neu,nb:ne) & - - (pdata%q(ivx,nb:ne,ne:neu,nb:ne) & - * pdata%q(ibx,nb:ne,ne:neu,nb:ne) & - + pdata%q(ivz,nb:ne,ne:neu,nb:ne) & - * pdata%q(ibz,nb:ne,ne:neu,nb:ne)) & - * pdata%q(iby,nb:ne,ne:neu,nb:ne)) * dxz + inarr(24) = inarr(24) & + - 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(27) = inarr(27) & + + 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(24) = inarr(24) - sum((pdata%q(ibx,nb:ne,ne:neu, : )**2 & - + pdata%q(ibz,nb:ne,ne:neu, : )**2) & - * pdata%q(ivy,nb:ne,ne:neu, : ) & - - (pdata%q(ivx,nb:ne,ne:neu, : ) & - * pdata%q(ibx,nb:ne,ne:neu, : ) & - + pdata%q(ivz,nb:ne,ne:neu, : ) & - * pdata%q(ibz,nb:ne,ne:neu, : )) & - * pdata%q(iby,nb:ne,ne:neu, : )) * dxz + inarr(24) = inarr(24) & + - sum(sum(pdata%q(ibx:ibz,nb:ne,ne:nep, : )**2,1) & + * pdata%q(ivy ,nb:ne,ne:nep, : )) * dxz + inarr(27) = inarr(27) & + + 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 @@ -1038,17 +1032,17 @@ module integrals #if NDIMS == 3 ! the Z component of current density ! - tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne ,nb:ne) & - - pdata%q(iby,nbl:nel,ne ,nb:ne)) / dh(1) & - - (pdata%q(ibx,nb :ne ,neu,nb:ne) & + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,ne ,nb:ne) & + - pdata%q(iby,nbm:nem,ne ,nb:ne)) / dh(1) & + - (pdata%q(ibx,nb :ne ,nep,nb:ne) & - pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2) ! the X component of current density - tmp(:,2,:) = 0.5d+00 * (pdata%q(ibx,nb :ne ,ne,nbu:neu) & - - pdata%q(ibx,nb :ne ,ne,nbl:nel)) / dh(3) & - - 0.5d+00 * (pdata%q(ibz,nbu:neu,ne,nb :ne ) & - - pdata%q(ibz,nbl:nel,ne,nb :ne )) / dh(1) + tmp(:,2,:) = 0.5d+00 * (pdata%q(ibx,nb :ne ,ne,nbp:nep) & + - pdata%q(ibx,nb :ne ,ne,nbm:nem)) / dh(3) & + - 0.5d+00 * (pdata%q(ibz,nbp:nep,ne,nb :ne ) & + - pdata%q(ibz,nbm:nem,ne,nb :ne )) / dh(1) ! diffusion of Bx through the upper Y boundary ! @@ -1057,21 +1051,21 @@ module integrals ! diffusion of magnetic energy through the upper Y boundary ! - inarr(27) = inarr(27) & - + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne,nb:ne) & + inarr(30) = inarr(30) & + - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne,nb:ne) & - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne,nb:ne)) * dxz #else /* NDIMS == 3 */ ! the Z component of current density ! - tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne , : ) & - - pdata%q(iby,nbl:nel,ne , : )) / dh(1) & - - (pdata%q(ibx,nb :ne ,neu, : ) & + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbp:nep,ne , : ) & + - pdata%q(iby,nbm:nem,ne , : )) / dh(1) & + - (pdata%q(ibx,nb :ne ,nep, : ) & - pdata%q(ibx,nb :ne ,ne , : )) / dh(2) ! the X component of current density ! - tmp(:,2,:) = 0.5d+00 * (pdata%q(ibz,nbl:nel,ne, : ) & - - pdata%q(ibz,nbu:neu,ne, : )) / dh(1) + tmp(:,2,:) = 0.5d+00 * (pdata%q(ibz,nbm:nem,ne, : ) & + - pdata%q(ibz,nbp:nep,ne, : )) / dh(1) ! diffusion of Bx through the upper Y boundary ! @@ -1080,8 +1074,8 @@ module integrals ! diffusion of magnetic energy through the upper Y boundary ! - inarr(27) = inarr(27) & - + sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) & + inarr(30) = inarr(30) & + - sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) & - tmp(:,2,:) * pdata%q(ibz,nb:ne,ne, : )) * dxz #endif /* NDIMS == 3 */ @@ -1100,41 +1094,40 @@ module integrals ! advection of Bx along Z ! - inarr(14) = inarr(14) + sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbl:nb)) & - * pdata%q(ivz,nb:ne,nb:ne,nbl:nb)) * dxy + 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,nbl:nb) & - * pdata%q(ivx,nb:ne,nb:ne,nbl:nb), & - pdata%q(ibx,nb:ne,nb:ne,nbl:nb))) * dxy + 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(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,nbl:nb)**2 & - + pdata%q(iby,nb:ne,nb:ne,nbl:nb)**2) & - * pdata%q(ivz,nb:ne,nb:ne,nbl:nb) & - - (pdata%q(ivx,nb:ne,nb:ne,nbl:nb) & - * pdata%q(ibx,nb:ne,nb:ne,nbl:nb) & - + pdata%q(ivy,nb:ne,nb:ne,nbl:nb) & - * pdata%q(iby,nb:ne,nb:ne,nbl:nb)) & - * pdata%q(ibz,nb:ne,nb:ne,nbl:nb)) * dxy + inarr(25) = inarr(25) & + + 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(28) = inarr(28) & + - 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 ! the X component of current density ! - tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,nb ) & - - pdata%q(ibz,nb:ne,nbl:nel,nb )) / dh(2) & + tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbp:nep,nb ) & + - pdata%q(ibz,nb:ne,nbm:nem,nb )) / dh(2) & - (pdata%q(iby,nb:ne,nb :ne ,nb ) & - - pdata%q(iby,nb:ne,nb :ne ,nbl)) / dh(3) + - pdata%q(iby,nb:ne,nb :ne ,nbm)) / dh(3) ! the Y component of current density ! tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nb ) & - - pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3) & - - 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,nb ) & - - pdata%q(ibz,nbl:nel,nb:ne,nb )) / dh(1) + - pdata%q(ibx,nb :ne ,nb:ne,nbm)) / dh(3) & + - 0.5d+00 * (pdata%q(ibz,nbp:nep,nb:ne,nb ) & + - pdata%q(ibz,nbm:nem,nb:ne,nb )) / dh(1) ! diffusion of Bx through the lower Z boundary ! @@ -1143,9 +1136,9 @@ module integrals ! diffusion of magnetic energy through the lower Z boundary ! - inarr(28) = inarr(28) & - - sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) & - - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy + inarr(31) = inarr(31) & + + sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) & + - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy end if ! resistivity @@ -1157,41 +1150,40 @@ module integrals ! advection of Bx along Z ! - inarr(14) = inarr(14) - sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:neu)) & - * pdata%q(ivz,nb:ne,nb:ne,ne:neu)) * dxy + 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:neu) & - * pdata%q(ivx,nb:ne,nb:ne,ne:neu), & - pdata%q(ibx,nb:ne,nb:ne,ne:neu))) * dxy + 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(25) = inarr(25) + sum((pdata%q(ibx,nb:ne,nb:ne,ne:neu)**2 & - + pdata%q(iby,nb:ne,nb:ne,ne:neu)**2) & - * pdata%q(ivz,nb:ne,nb:ne,ne:neu) & - - (pdata%q(ivx,nb:ne,nb:ne,ne:neu) & - * pdata%q(ibx,nb:ne,nb:ne,ne:neu) & - + pdata%q(ivy,nb:ne,nb:ne,ne:neu) & - * pdata%q(iby,nb:ne,nb:ne,ne:neu)) & - * pdata%q(ibz,nb:ne,nb:ne,ne:neu)) * dxy + inarr(25) = inarr(25) & + - 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(28) = inarr(28) & + + 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 ! the X component of current density ! - tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,ne ) & - - pdata%q(ibz,nb:ne,nbl:nel,ne )) / dh(2) & - - (pdata%q(iby,nb:ne,nb :ne ,neu ) & + tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbp:nep,ne ) & + - pdata%q(ibz,nb:ne,nbm:nem,ne )) / dh(2) & + - (pdata%q(iby,nb:ne,nb :ne ,nep) & - pdata%q(iby,nb:ne,nb :ne ,ne )) / dh(3) ! the Y component of current density ! - tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,neu) & + tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,nep) & - pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) & - - 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) & - - pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1) + - 0.5d+00 * (pdata%q(ibz,nbp:nep,nb:ne,ne ) & + - pdata%q(ibz,nbm:nem,nb:ne,ne )) / dh(1) ! diffusion of Bx through the upper Z boundary ! @@ -1200,9 +1192,9 @@ module integrals ! diffusion of magnetic energy through the upper Z boundary ! - inarr(28) = inarr(28) & - + sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) & - - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy + inarr(31) = inarr(31) & + - sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) & + - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy end if ! resistivity @@ -1211,6 +1203,92 @@ module integrals end if ! not periodic in Z #endif /* NDIMS == 3 */ +! get the conversion between the magnetic energy and kinetic and +! internal energies +! +#if NDIMS == 3 +! the X component of current density +! + tmp(:,:,:) = 5.0d-01 * (pdata%q(ibz,nb:ne,nbp:nep,nb :ne ) & + - pdata%q(ibz,nb:ne,nbm:nem,nb :ne )) / dh(2) & + - 5.0d-01 * (pdata%q(iby,nb:ne,nb :ne ,nbp:nep) & + - pdata%q(iby,nb:ne,nb :ne ,nbm:nem)) / dh(3) + + inarr(32) = inarr(32) & + + 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)) * tmp(:,:,:)) * dvol +#else /* NDIMS == 3 */ + tmp(:,:,:) = 5.0d-01 * (pdata%q(ibz,nb:ne,nbp:nep, : ) & + - pdata%q(ibz,nb:ne,nbm:nem, : )) / dh(2) + + inarr(32) = inarr(32) & + + 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, : )) * tmp(:,:,:)) * dvol +#endif /* NDIMS == 3 */ + + if (resistivity > 0.0d+00) & + inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol + +! the Y component of current density +! +#if NDIMS == 3 + tmp(:,:,:) = 5.0d-01 * (pdata%q(ibx,nb :ne ,nb:ne,nbp:nep) & + - pdata%q(ibx,nb :ne ,nb:ne,nbm:nem)) / dh(3) & + - 5.0d-01 * (pdata%q(ibz,nbp:nep,nb:ne,nb :ne ) & + - pdata%q(ibz,nbm:nem,nb:ne,nb :ne )) / dh(1) + + inarr(32) = inarr(32) & + + 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)) * tmp(:,:,:)) * dvol +#else /* NDIMS == 3 */ + tmp(:,:,:) = - 5.0d-01 * (pdata%q(ibz,nbp:nep,nb:ne, : ) & + - pdata%q(ibz,nbm:nem,nb:ne, : )) / dh(1) + + inarr(32) = inarr(32) & + + 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, : )) * tmp(:,:,:)) * dvol +#endif /* NDIMS == 3 */ + + if (resistivity > 0.0d+00) & + inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol + +! the Z component of current density +! +#if NDIMS == 3 + tmp(:,:,:) = (pdata%q(iby,nbp:nep,nb :ne ,nb:ne) & + - pdata%q(iby,nbm:nem,nb :ne ,nb:ne)) / dh(1) & + - (pdata%q(ibx,nb :ne ,nbp:nep,nb:ne) & + - pdata%q(ibx,nb :ne ,nbm:nem,nb:ne)) / dh(2) + + inarr(32) = inarr(32) & + + 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)) * tmp(:,:,:)) * dvol +#else /* NDIMS == 3 */ + tmp(:,:,:) = (pdata%q(iby,nbp:nep,nb :ne , : ) & + - pdata%q(iby,nbm:nem,nb :ne , : )) / dh(1) & + - (pdata%q(ibx,nb :ne ,nbp:nep, : ) & + - pdata%q(ibx,nb :ne ,nbm:nem, : )) / dh(2) + + inarr(32) = inarr(32) & + + 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, : )) * tmp(:,:,:)) * dvol +#endif /* NDIMS == 3 */ + + if (resistivity > 0.0d+00) & + inarr(33) = inarr(33) - sum(tmp(:,:,:)**2) * dvol + end if ! magnetized ! get the inflow speed @@ -1262,8 +1340,10 @@ module integrals inarr(13:16) = 5.0d-01 * inarr(13:16) inarr(17:18) = resistivity * inarr(17:18) inarr(22) = 2.5d-01 * inarr(22) / yarea - inarr(23:25) = 5.0d-01 * inarr(23:25) - inarr(26:28) = resistivity * inarr(26:28) + inarr(23:28) = 5.0d-01 * inarr(23:28) + inarr(29:31) = resistivity * inarr(29:31) + inarr(32) = 5.0d-01 * inarr(32) + inarr(33) = 2.5d-01 * resistivity * inarr(33) ! write down the integrals and statistics to appropriate files ! @@ -1278,7 +1358,7 @@ module integrals , avarr(6), mnarr(6), mxarr(6) & , avarr(7), mnarr(7), mxarr(7) write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:) - write(runit,"(i9, 19es25.15e3)") step, time, inarr(11:28) + write(runit,"(i9, 24es25.15e3)") step, time, inarr(11:33) call flush_and_sync(funit) call flush_and_sync(sunit)