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

INTEGRALS: Add terms for reconnection speed based on magnetic energy.

Browse Source 
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 <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2020-12-25 11:21:54 -03:00
parent d34b780525
commit 4ece753176
1 changed files with 185 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

253
sources/integrals.F90
View File

@ -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)