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

INTEGRALS: Add Bx infinity to reconnection measurements.

Browse Source 
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2019-03-08 07:37:37 -03:00
parent ed51af4dbd
commit 39d28df53d
1 changed files with 46 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
sources/integrals.F90
View File

@ -269,9 +269,9 @@ module integrals
! write the integral file header
!
write(runit,'("#",a8,11a18)') &
'step', 'time' &
, '|Bx| int', '|Bx| y-adv', '|Bx| y-shr', '|Bx| y-dif' &
write(runit,'("#",a8,12a18)') &
'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'
write(runit,"('#')")
@ -622,14 +622,24 @@ module integrals
!
if (pdata%meta%ymin == ymin) then
! advection of Bx along Y
! mean Bx at boundary
!
#if NDIMS == 3
inarr(12) = inarr(12) &
+ 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb,nb:ne))) * dxz
#else /* NDIMS == 3 */
inarr(12) = inarr(12) &
+ 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nb, : ))) * 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,nbl:nb,nb:ne)) &
* pdata%u(ivy,nb:ne,nbl:nb,nb:ne)) * dxz
#else /* NDIMS == 3 */
inarr(12) = inarr(12) &
inarr(13) = inarr(13) &
+ 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,nbl:nb, : )) &
* pdata%u(ivy,nb:ne,nbl:nb, : )) * dxz
#endif /* NDIMS == 3 */
@ -637,12 +647,12 @@ module integrals
! shear of By along X
!
#if NDIMS == 3
inarr(13) = inarr(13) &
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
#else /* NDIMS == 3 */
inarr(13) = inarr(13) &
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
@ -660,7 +670,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,nb ,nb:ne) &
- pdata%u(ibx,nb :ne ,nbl,nb:ne)) / dh(2)
inarr(14) = inarr(14) - resistivity * sum(sign(tmp(:,1,:), &
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
pdata%u(ibx,nb:ne,nb ,nb:ne))) * dxz
#else /* NDIMS == 3 */
! current density Z component
@ -670,7 +680,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,nb , : ) &
- pdata%u(ibx,nb :ne ,nbl, : )) / dh(2)
inarr(14) = inarr(14) - resistivity * sum(sign(tmp(:,1,:), &
inarr(15) = inarr(15) - resistivity * sum(sign(tmp(:,1,:), &
pdata%u(ibx,nb:ne,nb , : ))) * dxz
#endif /* NDIMS == 3 */
@ -682,14 +692,24 @@ module integrals
!
if (pdata%meta%ymax == ymax) then
! advection of Bx along Y
! mean Bx at boundary
!
#if NDIMS == 3
inarr(12) = inarr(12) &
+ 0.5d+00 * sum(abs(pdata%u(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
#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
#else /* NDIMS == 3 */
inarr(12) = inarr(12) &
inarr(13) = inarr(13) &
- 0.5d+00 * sum(abs(pdata%u(ibx,nb:ne,ne:neu, : )) &
* pdata%u(ivy,nb:ne,ne:neu, : )) * dxz
#endif /* NDIMS == 3 */
@ -697,12 +717,12 @@ module integrals
! shear of By along X
!
#if NDIMS == 3
inarr(13) = inarr(13) &
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
#else /* NDIMS == 3 */
inarr(13) = inarr(13) &
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
@ -720,7 +740,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,neu,nb:ne) &
- pdata%u(ibx,nb :ne ,ne ,nb:ne)) / dh(2)
inarr(14) = inarr(14) + resistivity * sum(sign(tmp(:,1,:), &
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%u(ibx,nb:ne,ne,nb:ne))) * dxz
#else /* NDIMS == 3 */
! current density Z component
@ -730,7 +750,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,neu, : ) &
- pdata%u(ibx,nb :ne ,ne , : )) / dh(2)
inarr(14) = inarr(14) + resistivity * sum(sign(tmp(:,1,:), &
inarr(15) = inarr(15) + resistivity * sum(sign(tmp(:,1,:), &
pdata%u(ibx,nb:ne,ne, : ))) * dxz
#endif /* NDIMS == 3 */
@ -747,13 +767,13 @@ module integrals
! advection of Bx along Z
!
inarr(15) = inarr(15) &
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
! shear of Bz along X
!
inarr(16) = inarr(16) &
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
@ -769,7 +789,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,nb:ne,nb ) &
- pdata%u(ibx,nb :ne ,nb:ne,nbl)) / dh(3)
inarr(17) = inarr(17) + resistivity * sum(sign(tmp(:,:,1), &
inarr(18) = inarr(18) + resistivity * sum(sign(tmp(:,:,1), &
pdata%u(ibx,nb:ne,nb:ne,nb))) * dxy
end if ! resistivity
@ -782,13 +802,13 @@ module integrals
! advection of Bx along Z
!
inarr(15) = inarr(15) &
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
! shear of Bz along X
!
inarr(16) = inarr(16) &
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
@ -804,7 +824,7 @@ module integrals
- (pdata%u(ibx,nb :ne ,nb:ne,neu) &
- pdata%u(ibx,nb :ne ,nb:ne,ne )) / dh(3)
inarr(17) = inarr(17) - resistivity * sum(sign(tmp(:,:,1), &
inarr(18) = inarr(18) - resistivity * sum(sign(tmp(:,:,1), &
pdata%u(ibx,nb:ne,nb:ne,ne))) * dxy
end if ! resistivity
@ -819,16 +839,16 @@ module integrals
!
if (pdata%meta%ymin == ymin) then
#if NDIMS == 3
inarr(18) = inarr(18) + sum(pdata%q(ivy,nb:ne,nb,nb:ne)) * dxz
inarr(19) = inarr(19) + sum(pdata%q(ivy,nb:ne,nb,nb:ne)) * dxz
#else /* NDIMS == 3 */
inarr(18) = inarr(18) + sum(pdata%q(ivy,nb:ne,nb, : )) * dxz
inarr(19) = inarr(19) + sum(pdata%q(ivy,nb:ne,nb, : )) * dxz
#endif /* NDIMS == 3 */
end if
if (pdata%meta%ymax == ymax) then
#if NDIMS == 3
inarr(19) = inarr(19) - sum(pdata%q(ivy,nb:ne,ne,nb:ne)) * dxz
inarr(20) = inarr(20) - sum(pdata%q(ivy,nb:ne,ne,nb:ne)) * dxz
#else /* NDIMS == 3 */
inarr(19) = inarr(19) - sum(pdata%q(ivy,nb:ne,ne, : )) * dxz
inarr(20) = inarr(20) - sum(pdata%q(ivy,nb:ne,ne, : )) * dxz
#endif /* NDIMS == 3 */
end if
@ -870,7 +890,7 @@ module integrals
, avarr(5), mnarr(5), mxarr(5) &
, avarr(6), mnarr(6), mxarr(6) &
, avarr(7), mnarr(7), mxarr(7)
write(runit,"(i9, 10es18.8e3)") step, time, inarr(11:19)
write(runit,"(i9, 11es18.8e3)") step, time, inarr(11:20)
end if
#ifdef PROFILE