From 4ece753176eb75db37c1a3188b5cb547dc8df5b4 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Fri, 25 Dec 2020 11:21:54 -0300 Subject: [PATCH] INTEGRALS: Add terms for reconnection speed based on magnetic energy. These terms help to measure how much of the magnetic energy is dissipated or converted to other energies in the system. Signed-off-by: Grzegorz Kowal --- sources/integrals.F90 | 253 ++++++++++++++++++++++++++++++------------ 1 file changed, 185 insertions(+), 68 deletions(-) diff --git a/sources/integrals.F90 b/sources/integrals.F90 index 1ec22d8..5a99fbf 100644 --- a/sources/integrals.F90 +++ b/sources/integrals.F90 @@ -331,11 +331,13 @@ module integrals ! write the integral file header ! - write(runit,'("#",a8,12a25)') & + write(runit,'("#",a8,17a25)') & 'step', 'time', '|Bx| int', '|Bx| inf' & , '|Bx| y-adv', '|Bx| y-shr', '|Bx| y-dif' & , '|Bx| z-adv', '|Bx| z-shr', '|Bx| z-dif' & - , 'Vin lower' , 'Vin upper' + , 'Vin lower' , 'Vin upper' & + , 'Emag', 'Emag(inf)' & + , 'Emag(x-adv)', 'Emag(y-adv)', 'Emag(z-adv)' write(runit,"('#')") end if ! store @@ -423,9 +425,8 @@ module integrals use coordinates , only : ni => ncells use coordinates , only : nb, nbl, nbu, ne, nel, neu use coordinates , only : adx, ady, adz, advol, voli -#if NDIMS == 3 use coordinates , only : periodic -#endif /* NDIMS == 3 */ + use coordinates , only : xmin, xmax, xarea use coordinates , only : ymin, ymax, yarea #if NDIMS == 3 use coordinates , only : zmin, zmax, zarea @@ -448,7 +449,7 @@ module integrals ! local variables ! - real(kind=8) :: dvol, dvolh, dxz + real(kind=8) :: dvol, dvolh, dyz, dxz #if NDIMS == 3 real(kind=8) :: dxy #endif /* NDIMS == 3 */ @@ -522,6 +523,7 @@ module integrals ! dvol = advol(pdata%meta%level) dvolh = 0.5d+00 * dvol + dyz = ady(pdata%meta%level) * adz(pdata%meta%level) / xarea #if NDIMS == 3 dxy = adx(pdata%meta%level) * ady(pdata%meta%level) / zarea #endif /* NDIMS == 3 */ @@ -687,6 +689,10 @@ module integrals inarr(11) = inarr(11) + sum(abs(pdata%u(ibx,nb:ne,nb:ne, : ))) * dvol #endif /* NDIMS == 3 */ +! magnetic energy +! + inarr(21) = inarr(7) + ! lower Y-boundary ! if (pdata%meta%ymin <= ymin + dh(2)) then @@ -705,26 +711,26 @@ module integrals ! #if NDIMS == 3 inarr(13) = inarr(13) & - + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nbl:nb,nb:ne)) & - * pdata%u(ivy,nb:ne,nbl:nb,nb:ne)) * dxz + + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nbl:nb,nb:ne)) & + * pdata%q(ivy,nb:ne,nbl:nb,nb:ne)) * dxz #else /* NDIMS == 3 */ inarr(13) = inarr(13) & - + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nbl:nb, : )) & - * pdata%u(ivy,nb:ne,nbl:nb, : )) * dxz + + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nbl:nb, : )) & + * pdata%q(ivy,nb:ne,nbl:nb, : )) * dxz #endif /* NDIMS == 3 */ ! shear of By along X ! #if NDIMS == 3 inarr(14) = inarr(14) & - - 0.5d+00 * sum(sign(pdata%u(iby,nb:ne,nbl:nb,nb:ne) & - * pdata%u(ivx,nb:ne,nbl:nb,nb:ne), & - pdata%u(ibx,nb:ne,nbl:nb,nb:ne))) * dxz + - 0.5d+00 * 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 #else /* NDIMS == 3 */ inarr(14) = inarr(14) & - - 0.5d+00 * sum(sign(pdata%u(iby,nb:ne,nbl:nb, : ) & - * pdata%u(ivx,nb:ne,nbl:nb, : ), & - pdata%u(ibx,nb:ne,nbl:nb, : ))) * dxz + - 0.5d+00 * sum(sign(pdata%q(iby,nb:ne,nbl:nb, : ) & + * pdata%q(ivx,nb:ne,nbl:nb, : ), & + pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz #endif /* NDIMS == 3 */ ! diffusion along Y @@ -734,27 +740,53 @@ module integrals #if NDIMS == 3 ! current density Z component ! - tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,nb ,nb:ne) & - - pdata%u(iby,nbl:nel,nb ,nb:ne)) / dh(1) & - - (pdata%u(ibx,nb :ne ,nb ,nb:ne) & - - pdata%u(ibx,nb :ne ,nbl,nb:ne)) / dh(2) + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb ,nb:ne) & + - pdata%q(iby,nbl:nel,nb ,nb:ne)) / dh(1) & + - (pdata%q(ibx,nb :ne ,nb ,nb:ne) & + - pdata%quibx,nb :ne ,nbl,nb:ne)) / dh(2) inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), & - pdata%u(ibx,nb:ne,nb ,nb:ne))) * dxz + pdata%q(ibx,nb:ne,nb ,nb:ne))) * dxz #else /* NDIMS == 3 */ ! current density Z component ! - tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,nb , : ) & - - pdata%u(iby,nbl:nel,nb , : )) / dh(1) & - - (pdata%u(ibx,nb :ne ,nb , : ) & - - pdata%u(ibx,nb :ne ,nbl, : )) / dh(2) + tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) & + - pdata%q(iby,nbl:nel,nb , : )) / dh(1) & + - (pdata%q(ibx,nb :ne ,nb , : ) & + - pdata%q(ibx,nb :ne ,nbl, : )) / dh(2) inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), & - pdata%u(ibx,nb:ne,nb , : ))) * dxz + pdata%q(ibx,nb:ne,nb , : ))) * dxz #endif /* NDIMS == 3 */ end if ! resistivity +! mean magnetic energy at boundary +! +#if NDIMS == 3 + inarr(22) = inarr(22) & + + 2.5d-01 * sum((pdata%u(ibx,nb:ne,nb,nb:ne)**2 & + + pdata%u(ibz,nb:ne,nb,nb:ne)**2)) * dxz +#else /* NDIMS == 3 */ + inarr(22) = inarr(22) & + + 2.5d-01 * sum((pdata%u(ibx,nb:ne,nb, : )**2 & + + pdata%u(ibz,nb:ne,nb, : )**2)) * dxz +#endif /* NDIMS == 3 */ + +! advection of magnetic field along Y +! +#if NDIMS == 3 + inarr(24) = inarr(24) & + + 5.0d-01 * 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)) * dxz +#else /* NDIMS == 3 */ + inarr(24) = inarr(24) & + + 5.0d-01 * sum((pdata%q(ibx,nb:ne,nbl:nb, : )**2 & + + pdata%q(ibz,nb:ne,nbl:nb, : )**2) & + * pdata%q(ivy,nb:ne,nbl:nb, : )) * dxz +#endif /* NDIMS == 3 */ + end if ! ymin ! upper Y-boundary @@ -765,36 +797,36 @@ module integrals ! #if NDIMS == 3 inarr(12) = inarr(12) & - + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,ne,nb:ne))) * dxz + + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,ne,nb:ne))) * dxz #else /* NDIMS == 3 */ inarr(12) = inarr(12) & - + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,ne, : ))) * dxz + + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,ne, : ))) * dxz #endif /* NDIMS == 3 */ ! advection of Bx along Y ! #if NDIMS == 3 inarr(13) = inarr(13) & - - 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,ne:neu,nb:ne)) & - * pdata%u(ivy,nb:ne,ne:neu,nb:ne)) * dxz + - 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,ne:neu,nb:ne)) & + * pdata%q(ivy,nb:ne,ne:neu,nb:ne)) * dxz #else /* NDIMS == 3 */ inarr(13) = inarr(13) & - - 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,ne:neu, : )) & - * pdata%u(ivy,nb:ne,ne:neu, : )) * dxz + - 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,ne:neu, : )) & + * pdata%q(ivy,nb:ne,ne:neu, : )) * dxz #endif /* NDIMS == 3 */ ! shear of By along X ! #if NDIMS == 3 inarr(14) = inarr(14) & - + 0.5d+00 * sum(sign(pdata%u(iby,nb:ne,ne:neu,nb:ne) & - * pdata%u(ivx,nb:ne,ne:neu,nb:ne), & - pdata%u(ibx,nb:ne,ne:neu,nb:ne))) * dxz + + 0.5d+00 * 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 #else /* NDIMS == 3 */ inarr(14) = inarr(14) & - + 0.5d+00 * sum(sign(pdata%u(iby,nb:ne,ne:neu, : ) & - * pdata%u(ivx,nb:ne,ne:neu, : ), & - pdata%u(ibx,nb:ne,ne:neu, : ))) * dxz + + 0.5d+00 * sum(sign(pdata%q(iby,nb:ne,ne:neu, : ) & + * pdata%q(ivx,nb:ne,ne:neu, : ), & + pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz #endif /* NDIMS == 3 */ ! diffusion along Y @@ -804,29 +836,99 @@ module integrals #if NDIMS == 3 ! current density Z component ! - tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,ne ,nb:ne) & - - pdata%u(iby,nbl:nel,ne ,nb:ne)) / dh(1) & - - (pdata%u(ibx,nb :ne ,neu,nb:ne) & - - pdata%u(ibx,nb :ne ,ne ,nb:ne)) / dh(2) + 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) & + - pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2) inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), & - pdata%u(ibx,nb:ne,ne,nb:ne))) * dxz + pdata%q(ibx,nb:ne,ne,nb:ne))) * dxz #else /* NDIMS == 3 */ ! current density Z component ! - tmp(:,1,:) = 0.5d+00 * (pdata%u(iby,nbu:neu,ne , : ) & - - pdata%u(iby,nbl:nel,ne , : )) / dh(1) & - - (pdata%u(ibx,nb :ne ,neu, : ) & - - pdata%u(ibx,nb :ne ,ne , : )) / dh(2) + 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, : ) & + - pdata%q(ibx,nb :ne ,ne , : )) / dh(2) inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), & - pdata%u(ibx,nb:ne,ne, : ))) * dxz + pdata%q(ibx,nb:ne,ne, : ))) * dxz #endif /* NDIMS == 3 */ end if ! resistivity +! mean magnetic energy at boundary +! +#if NDIMS == 3 + inarr(22) = inarr(22) & + + 2.5d-01 * sum((pdata%u(ibx,nb:ne,ne,nb:ne)**2 & + + pdata%u(ibz,nb:ne,ne,nb:ne)**2)) * dxz +#else /* NDIMS == 3 */ + inarr(22) = inarr(22) & + + 2.5d-01 * sum((pdata%u(ibx,nb:ne,ne, : )**2 & + + pdata%u(ibz,nb:ne,ne, : )**2)) * dxz +#endif /* NDIMS == 3 */ + +! advection of magnetic field along Y +! +#if NDIMS == 3 + inarr(24) = inarr(24) & + - 5.0d-01 * 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)) * dxz +#else /* NDIMS == 3 */ + inarr(24) = inarr(24) & + - 5.0d-01 * sum((pdata%q(ibx,nb:ne,ne:neu, : )**2 & + + pdata%q(ibz,nb:ne,ne:neu, : )**2) & + * pdata%q(ivy,nb:ne,ne:neu, : )) * dxz +#endif /* NDIMS == 3 */ + end if ! ymax + if (.not. periodic(1)) then + +! lower X-boundary +! + if (pdata%meta%xmin <= xmin + dh(1)) then + +! advection of magnetic field along X +! +#if NDIMS == 3 + inarr(23) = inarr(23) & + + 5.0d-01 * 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)) * dyz +#else /* NDIMS == 3 */ + inarr(23) = inarr(23) & + + 5.0d-01 * sum((pdata%q(iby,nbl:nb,nb:ne, : )**2 & + + pdata%q(ibz,nbl:nb,nb:ne, : )**2) & + * pdata%q(ivx,nbl:nb,nb:ne, : )) * dyz +#endif /* NDIMS == 3 */ + + end if ! xmin + +! upper X-boundary +! + if (pdata%meta%xmax >= xmax - dh(1)) then + +! advection of magnetic field along X +! +#if NDIMS == 3 + inarr(23) = inarr(23) & + - 5.0d-01 * 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)) * dyz +#else /* NDIMS == 3 */ + inarr(23) = inarr(23) & + - 5.0d-01 * sum((pdata%q(iby,ne:neu,nb:ne, : )**2 & + + pdata%q(ibz,ne:neu,nb:ne, : )**2) & + * pdata%q(ivx,ne:neu,nb:ne, : )) * dyz +#endif /* NDIMS == 3 */ + + end if ! xmax + + end if ! not periodic in X + #if NDIMS == 3 if (.not. periodic(3)) then @@ -837,15 +939,15 @@ module integrals ! advection of Bx along Z ! inarr(16) = inarr(16) & - + 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb:ne,nbl:nb)) & - * pdata%u(ivz,nb:ne,nb:ne,nbl:nb)) * dxy + + 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nb:ne,nbl:nb)) & + * pdata%q(ivz,nb:ne,nb:ne,nbl:nb)) * dxy ! shear of Bz along X ! inarr(17) = inarr(17) & - - 0.5d+00 * sum(sign(pdata%u(ibz,nb:ne,nb:ne,nbl:nb) & - * pdata%u(ivx,nb:ne,nb:ne,nbl:nb), & - pdata%u(ibx,nb:ne,nb:ne,nbl:nb))) * dxy + - 0.5d+00 * 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 ! diffusion along Z ! @@ -853,16 +955,23 @@ module integrals ! current density Y component ! - tmp(:,:,1) = 0.5d+00 * (pdata%u(ibz,nbu:neu,nb:ne,nb ) & - - pdata%u(ibz,nbl:nel,nb:ne,nb )) / dh(1) & - - (pdata%u(ibx,nb :ne ,nb:ne,nb ) & - - pdata%u(ibx,nb :ne ,nb:ne,nbl)) / dh(3) + tmp(:,:,1) = 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,nb ) & + - pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3) inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), & - pdata%u(ibx,nb:ne,nb:ne,nb))) * dxy + pdata%q(ibx,nb:ne,nb:ne,nb))) * dxy end if ! resistivity +! advection of magnetic field along Z +! + inarr(25) = inarr(25) & + + 0.5d+00 * 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)) * dxy + end if ! zmin ! upper Z-boundary @@ -872,15 +981,15 @@ module integrals ! advection of Bx along Z ! inarr(16) = inarr(16) & - - 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb:ne,ne:neu)) & - * pdata%u(ivz,nb:ne,nb:ne,ne:neu)) * dxy + - 0.5d+00 * sum(abs(pdata%q(ibx,nb:ne,nb:ne,ne:neu)) & + * pdata%q(ivz,nb:ne,nb:ne,ne:neu)) * dxy ! shear of Bz along X ! inarr(17) = inarr(17) & - + 0.5d+00 * sum(sign(pdata%u(ibz,nb:ne,nb:ne,ne:neu) & - * pdata%u(ivx,nb:ne,nb:ne,ne:neu), & - pdata%u(ibx,nb:ne,nb:ne,ne:neu))) * dxy + + 0.5d+00 * 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 ! diffusion along Z ! @@ -888,17 +997,25 @@ module integrals ! current density Y component ! - tmp(:,:,1) = 0.5d+00 * (pdata%u(ibz,nbu:neu,nb:ne,ne ) & - - pdata%u(ibz,nbl:nel,nb:ne,ne )) / dh(1) & - - (pdata%u(ibx,nb :ne ,nb:ne,neu) & - - pdata%u(ibx,nb :ne ,nb:ne,ne )) / dh(3) + tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) & + - pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1) & + - (pdata%q(ibx,nb :ne ,nb:ne,neu) & + - pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,1), & - pdata%u(ibx,nb:ne,nb:ne,ne))) * dxy + pdata%q(ibx,nb:ne,nb:ne,ne))) * dxy end if ! resistivity +! advection of magnetic field along Z +! + inarr(25) = inarr(25) & + - 0.5d+00 * 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)) * dxy + end if ! zmax + end if ! not periodic in Z #endif /* NDIMS == 3 */ @@ -960,7 +1077,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, 11es25.15e3)") step, time, inarr(11:20) + write(runit,"(i9, 16es25.15e3)") step, time, inarr(11:25) call flush_and_sync(funit) call flush_and_sync(sunit)