gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

INTEGRALS: Implement remainig fluxes of magnetic energy.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2020-12-26 18:04:48 -03:00
parent a30bfa0dd8
commit ffb1a2b6c8
1 changed files with 338 additions and 138 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

476
sources/integrals.F90
View File

@ -331,13 +331,14 @@ module integrals
! write the integral file header ! write the integral file header
! !
write(runit,'("#",a8,17a25)') & write(runit,'("#",a8,20a25)') &
'step', 'time', '|Bx| int', '|Bx| inf' & 'step', 'time', '|Bx| int', '|Bx| inf' &
, '|Bx| y-adv', '|Bx| y-shr', '|Bx| y-dif' & , '|Bx| y-adv', '|Bx| y-shr', '|Bx| y-dif' &
, '|Bx| z-adv', '|Bx| z-shr', '|Bx| z-dif' & , '|Bx| z-adv', '|Bx| z-shr', '|Bx| z-dif' &
, 'Vin lower' , 'Vin upper' & , 'Vin lower' , 'Vin upper' &
, 'Emag', 'Emag(inf)' & , 'Emag', 'Emag(inf)' &
, 'Emag(x-adv)', 'Emag(y-adv)', 'Emag(z-adv)' , 'Emag(x-adv)', 'Emag(y-adv)', 'Emag(z-adv)' &
, 'Emag(x-dif)', 'Emag(y-dif)', 'Emag(z-dif)'
write(runit,"('#')") write(runit,"('#')")
end if ! store end if ! store
@ -426,10 +427,10 @@ module integrals
use coordinates , only : nb, nbl, nbu, ne, nel, neu use coordinates , only : nb, nbl, nbu, ne, nel, neu
use coordinates , only : adx, ady, adz, advol, voli use coordinates , only : adx, ady, adz, advol, voli
use coordinates , only : periodic use coordinates , only : periodic
use coordinates , only : xmin, xmax, xarea use coordinates , only : xmin, xmax
use coordinates , only : ymin, ymax, yarea use coordinates , only : ymin, ymax, yarea
#if NDIMS == 3 #if NDIMS == 3
use coordinates , only : zmin, zmax, zarea use coordinates , only : zmin, zmax
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp use equations , only : idn, ipr, ivx, ivy, ivz, ibx, iby, ibz, ibp
use equations , only : ien, imx, imy, imz use equations , only : ien, imx, imy, imz
@ -523,11 +524,11 @@ module integrals
! !
dvol = advol(pdata%meta%level) dvol = advol(pdata%meta%level)
dvolh = 0.5d+00 * dvol dvolh = 0.5d+00 * dvol
dyz = ady(pdata%meta%level) * adz(pdata%meta%level) / xarea dyz = ady(pdata%meta%level) * adz(pdata%meta%level)
#if NDIMS == 3 #if NDIMS == 3
dxy = adx(pdata%meta%level) * ady(pdata%meta%level) / zarea dxy = adx(pdata%meta%level) * ady(pdata%meta%level)
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
dxz = adx(pdata%meta%level) * adz(pdata%meta%level) / yarea dxz = adx(pdata%meta%level) * adz(pdata%meta%level)
dh(1) = adx(pdata%meta%level) dh(1) = adx(pdata%meta%level)
dh(2) = ady(pdata%meta%level) dh(2) = ady(pdata%meta%level)
dh(3) = adz(pdata%meta%level) dh(3) = adz(pdata%meta%level)
@ -693,6 +694,158 @@ module integrals
! !
inarr(21) = inarr(7) inarr(21) = 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(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
#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
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then
#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) &
- (pdata%q(ibz,nb ,nb:ne,nb :ne ) &
- pdata%q(ibz,nbl,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)
! diffusion of magnetic energy through the lower X boundary
!
inarr(26) = inarr(26) &
+ 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)
! 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)
! diffusion of magnetic energy through the lower X boundary
!
inarr(26) = inarr(26) &
+ sum(tmp(1,:,:) * pdata%q(ibz,nb,nb:ne, : ) &
- tmp(2,:,:) * pdata%q(iby,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(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
#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
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then
#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)
! the Z component of current density
!
tmp(2,:,:) = (pdata%q(iby,neu,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)
! diffusion of magnetic energy through the upper X boundary
!
inarr(26) = inarr(26) &
- 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, : ) &
- pdata%q(ibz,ne ,nb:ne, : )) / dh(1)
! the Z component of current density
!
tmp(2,:,:) = (pdata%q(iby,neu,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)
! diffusion of magnetic energy through the upper X boundary
!
inarr(26) = inarr(26) &
- sum(tmp(1,:,:) * pdata%q(ibz,ne,nb:ne, : ) &
- tmp(2,:,:) * pdata%q(iby,ne,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 ! lower Y-boundary
! !
if (pdata%meta%ymin <= ymin + dh(2)) then if (pdata%meta%ymin <= ymin + dh(2)) then
@ -733,60 +886,86 @@ module integrals
pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz pdata%q(ibx,nb:ne,nbl:nb, : ))) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! diffusion along Y ! mean magnetic energy at the lower Y boundary
! !
#if NDIMS == 3
inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,nb,nb:ne)**2) * dxz
#else /* NDIMS == 3 */
inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,nb, : )**2) * dxz
#endif /* NDIMS == 3 */
! 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
#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
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
#if NDIMS == 3 #if NDIMS == 3
! current density Z component ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb ,nb:ne) & 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(iby,nbl:nel,nb ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb ,nb:ne) & - (pdata%q(ibx,nb :ne ,nb ,nb:ne) &
- pdata%q(ibx,nb :ne ,nbl,nb:ne)) / dh(2) - pdata%q(ibx,nb :ne ,nbl,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)
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), & inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,nb ,nb:ne))) * dxz pdata%q(ibx,nb:ne,nb,nb:ne))) * dxz
! diffusion of magnetic energy through the lower Y boundary
!
inarr(27) = inarr(27) &
+ 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 */ #else /* NDIMS == 3 */
! current density Z component ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,nb , : ) &
- pdata%q(iby,nbl:nel,nb , : )) / dh(1) & - pdata%q(iby,nbl:nel,nb , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,nb , : ) & - (pdata%q(ibx,nb :ne ,nb , : ) &
- pdata%q(ibx,nb :ne ,nbl, : )) / dh(2) - pdata%q(ibx,nb :ne ,nbl, : )) / dh(2)
! the X component of current density
!
tmp(:,2,:) = (pdata%q(ibz,nb:ne,nb , : ) &
- pdata%q(ibz,nb:ne,nbl, : )) / dh(2)
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), & inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,nb , : ))) * dxz pdata%q(ibx,nb:ne,nb , : ))) * dxz
! diffusion of magnetic energy through the lower Y boundary
!
inarr(27) = inarr(27) &
+ sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,nb, : ) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,nb, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end if ! resistivity end if ! resistivity
! mean magnetic energy at boundary
!
#if NDIMS == 3
inarr(22) = inarr(22) &
+ 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb,nb:ne)**2 &
+ pdata%q(ibz,nb:ne,nb,nb:ne)**2)) * dxz
#else /* NDIMS == 3 */
inarr(22) = inarr(22) &
+ 2.5d-01 * sum((pdata%q(ibx,nb:ne,nb, : )**2 &
+ pdata%q(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 end if ! ymin
! upper Y-boundary ! upper Y-boundary
@ -829,107 +1008,91 @@ module integrals
pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz pdata%q(ibx,nb:ne,ne:neu, : ))) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
! diffusion along Y ! mean magnetic energy at the upper Y boundary
! !
#if NDIMS == 3
inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,ne,nb:ne)**2) * dxz
#else /* NDIMS == 3 */
inarr(22) = inarr(22) + sum(pdata%q(ibx:ibz,nb:ne,ne, : )**2) * dxz
#endif /* NDIMS == 3 */
! 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
#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
#endif /* NDIMS == 3 */
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
#if NDIMS == 3 #if NDIMS == 3
! current density Z component ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne ,nb:ne) & 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(iby,nbl:nel,ne ,nb:ne)) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu,nb:ne) & - (pdata%q(ibx,nb :ne ,neu,nb:ne) &
- pdata%q(ibx,nb :ne ,ne ,nb:ne)) / dh(2) - 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)
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), & inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,ne,nb:ne))) * dxz pdata%q(ibx,nb:ne,ne,nb:ne))) * dxz
! diffusion of magnetic energy through the upper Y boundary
!
inarr(27) = inarr(27) &
- 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 */ #else /* NDIMS == 3 */
! current density Z component ! the Z component of current density
! !
tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne , : ) & tmp(:,1,:) = 0.5d+00 * (pdata%q(iby,nbu:neu,ne , : ) &
- pdata%q(iby,nbl:nel,ne , : )) / dh(1) & - pdata%q(iby,nbl:nel,ne , : )) / dh(1) &
- (pdata%q(ibx,nb :ne ,neu, : ) & - (pdata%q(ibx,nb :ne ,neu, : ) &
- pdata%q(ibx,nb :ne ,ne , : )) / dh(2) - 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)
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), & inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%q(ibx,nb:ne,ne, : ))) * dxz pdata%q(ibx,nb:ne,ne, : ))) * dxz
! diffusion of magnetic energy through the upper Y boundary
!
inarr(27) = inarr(27) &
- sum(tmp(:,1,:) * pdata%q(ibx,nb:ne,ne, : ) &
- tmp(:,2,:) * pdata%q(ibz,nb:ne,ne, : )) * dxz
#endif /* NDIMS == 3 */ #endif /* NDIMS == 3 */
end if ! resistivity end if ! resistivity
! mean magnetic energy at boundary
!
#if NDIMS == 3
inarr(22) = inarr(22) &
+ 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne,nb:ne)**2 &
+ pdata%q(ibz,nb:ne,ne,nb:ne)**2)) * dxz
#else /* NDIMS == 3 */
inarr(22) = inarr(22) &
+ 2.5d-01 * sum((pdata%q(ibx,nb:ne,ne, : )**2 &
+ pdata%q(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 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 #if NDIMS == 3
inarr(23) = inarr(23) & ! fluxes across the Z boundaries
+ 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 if (.not. periodic(3)) then
! lower Z-boundary ! lower Z-boundary
@ -949,28 +1112,43 @@ module integrals
* pdata%q(ivx,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 pdata%q(ibx,nb:ne,nb:ne,nbl:nb))) * dxy
! diffusion along Z ! 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
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
! current density Y component ! the X component of current density
! !
tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,nb ) & tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,nb ) &
- pdata%q(ibz,nbl:nel,nb:ne,nb )) / dh(1) & - pdata%q(ibz,nb:ne,nbl:nel,nb )) / dh(2) &
- (pdata%q(ibx,nb :ne ,nb:ne,nb ) & - (pdata%q(iby,nb:ne,nb :ne ,nb ) &
- pdata%q(ibx,nb :ne ,nb:ne,nbl)) / dh(3) - pdata%q(iby,nb:ne,nb :ne ,nbl)) / dh(3)
inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), & ! 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)
inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,2), &
pdata%q(ibx,nb:ne,nb:ne,nb))) * dxy pdata%q(ibx,nb:ne,nb:ne,nb))) * dxy
end if ! resistivity ! diffusion of magnetic energy through the lower Z boundary
! advection of magnetic field along Z
! !
inarr(25) = inarr(25) & inarr(28) = inarr(28) &
+ 0.5d+00 * sum((pdata%q(ibx,nb:ne,nb:ne,nbl:nb)**2 & + sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,nb) &
+ pdata%q(iby,nb:ne,nb:ne,nbl:nb)**2) & - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,nb)) * dxy
* pdata%q(ivz,nb:ne,nb:ne,nbl:nb)) * dxy
end if ! resistivity
end if ! zmin end if ! zmin
@ -991,28 +1169,43 @@ module integrals
* pdata%q(ivx,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 pdata%q(ibx,nb:ne,nb:ne,ne:neu))) * dxy
! diffusion along Z ! 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
if (resistivity > 0.0d+00) then if (resistivity > 0.0d+00) then
! current density Y component ! the X component of current density
! !
tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nbu:neu,nb:ne,ne ) & tmp(:,:,1) = 0.5d+00 * (pdata%q(ibz,nb:ne,nbu:neu,ne ) &
- pdata%q(ibz,nbl:nel,nb:ne,ne )) / dh(1) & - pdata%q(ibz,nb:ne,nbl:nel,ne )) / dh(2) &
- (pdata%q(ibx,nb :ne ,nb:ne,neu) & - (pdata%q(iby,nb:ne,nb :ne ,neu ) &
- pdata%q(ibx,nb :ne ,nb:ne,ne )) / dh(3) - pdata%q(iby,nb:ne,nb :ne ,ne )) / dh(3)
inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,1), & ! the Y component of current density
!
tmp(:,:,2) = (pdata%q(ibx,nb :ne ,nb:ne,neu) &
- 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)
inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), &
pdata%q(ibx,nb:ne,nb:ne,ne))) * dxy pdata%q(ibx,nb:ne,nb:ne,ne))) * dxy
end if ! resistivity ! diffusion of magnetic energy through the upper Z boundary
! advection of magnetic field along Z
! !
inarr(25) = inarr(25) & inarr(28) = inarr(28) &
- 0.5d+00 * sum((pdata%q(ibx,nb:ne,nb:ne,ne:neu)**2 & - sum(tmp(:,:,1) * pdata%q(iby,nb:ne,nb:ne,ne) &
+ pdata%q(iby,nb:ne,nb:ne,ne:neu)**2) & - tmp(:,:,2) * pdata%q(ibx,nb:ne,nb:ne,ne)) * dxy
* pdata%q(ivz,nb:ne,nb:ne,ne:neu)) * dxy
end if ! resistivity
end if ! zmax end if ! zmax
@ -1064,6 +1257,13 @@ module integrals
! !
avarr(:) = avarr(:) * voli avarr(:) = avarr(:) * voli
! apply factors to the reconnection rate terms
!
inarr(12) = inarr(12) / yarea
inarr(22) = 2.5d-01 * inarr(22) / yarea
inarr(23:25) = 5.0d-01 * inarr(23:25)
inarr(26:28) = resistivity * inarr(26:28)
! write down the integrals and statistics to appropriate files ! write down the integrals and statistics to appropriate files
! !
if (stored) then if (stored) then
@ -1077,7 +1277,7 @@ module integrals
, avarr(6), mnarr(6), mxarr(6) & , avarr(6), mnarr(6), mxarr(6) &
, avarr(7), mnarr(7), mxarr(7) , avarr(7), mnarr(7), mxarr(7)
write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:) write(eunit,efmt) step, time, dtn, dte, maxval(errors(:)), errors(:)
write(runit,"(i9, 16es25.15e3)") step, time, inarr(11:25) write(runit,"(i9, 19es25.15e3)") step, time, inarr(11:28)
call flush_and_sync(funit) call flush_and_sync(funit)
call flush_and_sync(sunit) call flush_and_sync(sunit)