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)