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

Merge branch 'master' into flux-tubes

Browse Source
This commit is contained in:
Grzegorz Kowal 2023-08-29 10:53:47 -03:00
commit b624595682
2 changed files with 124 additions and 94 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
sources/evolution.F90
View File

@ -1090,7 +1090,6 @@ module evolution
! references
!
use blocks , only : set_blocks_update
use coordinates, only : toplev
use forcing , only : update_forcing, forcing_enabled
use mesh , only : update_mesh
@ -1117,10 +1116,6 @@ module evolution
!
if (toplev > 1) then
! set all meta blocks to not be updated
!
call set_blocks_update(.false.)
! check refinement and refine
!
call update_mesh(status)
@ -1131,10 +1126,6 @@ module evolution
call update_variables(time + dt, 0.0d+00, status)
if (status /= 0) go to 100
! set all meta blocks to be updated
!
call set_blocks_update(.true.)
end if ! toplev > 1
! error entry point
@ -1736,9 +1727,6 @@ module evolution
integer :: l, n, status
real(kind=8) :: t, ds
real(kind=8), parameter :: f21 = 3.0d+00 / 4.0d+00, f22 = 1.0d+00 / 4.0d+00
real(kind=8), parameter :: f31 = 1.0d+00 / 3.0d+00, f32 = 2.0d+00 / 3.0d+00
!-------------------------------------------------------------------------------
!
status = 1
@ -1780,8 +1768,8 @@ module evolution
!= 2nd step: U(2) = 3/4 U(n) + 1/4 U(1) + 1/4 dt F[U(1)]
!
ds = f22 * dt
t = time + 0.5d+00 * dt
ds = dt / 4.0d+00
t = time + 5.0d-01 * dt
!$omp parallel do default(shared) private(pdata)
do l = 1, n
@ -1800,8 +1788,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f21 * pdata%uu(:,:,:,:,1) &
+ f22 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = (3.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:)) / 4.0d+00
end do
!$omp end parallel do
@ -1811,7 +1800,7 @@ module evolution
!= 3rd step: U(n+1) = 1/3 U(n) + 2/3 U(2) + 2/3 dt F[U(2)]
!
ds = f32 * dt
ds = 2.0d+00 * dt / 3.0d+00
t = time + dt
!$omp parallel do default(shared) private(pdata)
@ -1831,8 +1820,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = f31 * pdata%uu(:,:,:,:,1) &
+ f32 * pdata%uu(:,:,:,:,2) + ds * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = ( pdata%uu(:,:,:,:,1) &
+ 2.0d+00 * (pdata%uu(:,:,:,:,2) &
+ ds * pdata%du(:,:,:,:))) / 3.0d+00
pdata%u => pdata%uu(:,:,:,:,1)
end do
@ -1854,7 +1844,7 @@ module evolution
100 continue
if (status /= 0) dt = 2.5d-01 * dt
if (status /= 0) dt = dt / 4.0d+00
end do
@ -1904,9 +1894,6 @@ module evolution
integer :: l, n, status
real(kind=8) :: t, ds
real(kind=8), parameter :: b1 = 2.0d+00 / 3.0d+00, b2 = 1.0d+00 / 3.0d+00, &
b3 = 1.0d+00 / 6.0d+00
!-------------------------------------------------------------------------------
!
status = 1
@ -1997,9 +1984,9 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,2) = b1 * pdata%uu(:,:,:,:,1) &
+ b2 * pdata%uu(:,:,:,:,2) &
+ b3 * dt * pdata%du(:,:,:,:)
pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2) &
+ 5.0d-01 * dt * pdata%du(:,:,:,:)) / 3.0d+00
end do
!$omp end parallel do
@ -3036,9 +3023,6 @@ module evolution
integer :: i, l, n, nrej, status
real(kind=8) :: tm, dtm, dh, fc
real(kind=8), parameter :: onethird = 1.0d+00 / 3.0d+00, &
twothird = 2.0d+00 / 3.0d+00
character(len=*), parameter :: loc = 'EVOLUTION:evolve_ssp324()'
!-------------------------------------------------------------------------------
@ -3113,10 +3097,10 @@ module evolution
do l = 1, n
pdata => data_blocks(l)%ptr
pdata%uu(:,:,:,:,3) = onethird * pdata%uu(:,:,:,:,1) &
+ twothird * pdata%uu(:,:,:,:,2)
pdata%uu(:,:,:,:,2) = twothird * pdata%uu(:,:,:,:,1) &
+ onethird * pdata%uu(:,:,:,:,2)
pdata%uu(:,:,:,:,3) = ( pdata%uu(:,:,:,:,1) &
+ 2.0d+00 * pdata%uu(:,:,:,:,2)) / 3.0d+00
pdata%uu(:,:,:,:,2) = (2.0d+00 * pdata%uu(:,:,:,:,1) &
+ pdata%uu(:,:,:,:,2)) / 3.0d+00
end do
!$omp end parallel do
@ -3173,9 +3157,9 @@ module evolution
pdata => data_blocks(l)%ptr
pdata%uu(ibp,:,:,:,2) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,2)
* pdata%uu(ibp,:,:,:,2)
pdata%uu(ibp,:,:,:,3) = adecay(pdata%meta%level) &
* pdata%uu(ibp,:,:,:,3)
* pdata%uu(ibp,:,:,:,3)
end do
!$omp end parallel do
end if
@ -3224,8 +3208,8 @@ module evolution
call update_variables(tm, dtm, status)
if (status /= 0) &
call print_message(loc, "Possible bug: the revert to " // &
"the previous state found it unphysical.")
call print_message(loc, "Possible bug: " // &
"the previous state seems unphysical.")
end if
niterations = niterations + 1

160
sources/statistics.F90
View File

@ -667,7 +667,7 @@ module statistics
#else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(idn,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */
avarr(1) = avarr(1) + sum(tmp(:,:,:)) * dvol
avarr(1) = avarr(1) + accsum(tmp(:,:,:)) * dvol
mnarr(1) = min(mnarr(1), minval(tmp(:,:,:)))
mxarr(1) = max(mxarr(1), maxval(tmp(:,:,:)))
@ -679,7 +679,7 @@ module statistics
#else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(ipr,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */
avarr(2) = avarr(2) + sum(tmp(:,:,:)) * dvol
avarr(2) = avarr(2) + accsum(tmp(:,:,:)) * dvol
mnarr(2) = min(mnarr(2), minval(tmp(:,:,:)))
mxarr(2) = max(mxarr(2), maxval(tmp(:,:,:)))
end if
@ -691,7 +691,7 @@ module statistics
#else /* NDIMS == 3 */
vel(:,:,:) = sqrt(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : )**2, 1))
#endif /* NDIMS == 3 */
avarr(3) = avarr(3) + sum(vel(:,:,:)) * dvol
avarr(3) = avarr(3) + accsum(vel(:,:,:)) * dvol
mnarr(3) = min(mnarr(3), minval(vel(:,:,:)))
mxarr(3) = max(mxarr(3), maxval(vel(:,:,:)))
@ -704,7 +704,7 @@ module statistics
#else /* NDIMS == 3 */
mag(:,:,:) = sqrt(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : )**2, 1))
#endif /* NDIMS == 3 */
avarr(4) = avarr(4) + sum(mag(:,:,:)) * dvol
avarr(4) = avarr(4) + accsum(mag(:,:,:)) * dvol
mnarr(4) = min(mnarr(4), minval(mag(:,:,:)))
mxarr(4) = max(mxarr(4), maxval(mag(:,:,:)))
@ -713,7 +713,7 @@ module statistics
#else /* NDIMS == 3 */
tmp(:,:,:) = pdata%q(ibp,nb:ne,nb:ne, : )
#endif /* NDIMS == 3 */
avarr(5) = avarr(5) + sum(tmp(:,:,:)) * dvol
avarr(5) = avarr(5) + accsum(tmp(:,:,:)) * dvol
mnarr(5) = min(mnarr(5), minval(tmp(:,:,:)))
mxarr(5) = max(mxarr(5), maxval(tmp(:,:,:)))
end if
@ -738,7 +738,7 @@ module statistics
else
tmp(:,:,:) = vel(:,:,:) / csnd
end if
avarr(6) = avarr(6) + sum(tmp(:,:,:)) * dvol
avarr(6) = avarr(6) + accsum(tmp(:,:,:)) * dvol
mnarr(6) = min(mnarr(6), minval(tmp(:,:,:)))
mxarr(6) = max(mxarr(6), maxval(tmp(:,:,:)))
@ -746,7 +746,7 @@ module statistics
!
if (magnetized) then
tmp(:,:,:) = sqd(:,:,:) * vel(:,:,:) / max(eps, mag(:,:,:))
avarr(7) = avarr(7) + sum(tmp(:,:,:)) * dvol
avarr(7) = avarr(7) + accsum(tmp(:,:,:)) * dvol
mnarr(7) = min(mnarr(7), minval(tmp(:,:,:)))
mxarr(7) = max(mxarr(7), maxval(tmp(:,:,:)))
end if
@ -775,32 +775,32 @@ module statistics
! total mass
!
#if NDIMS == 3
inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
inarr(1) = inarr(1) + sum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol
inarr(1) = inarr(1) + accsum(pdata%q(idn,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
! sum up kinetic energy
!
#if NDIMS == 3
inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh
inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ivy,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ivz,nb:ne,nb:ne,nb:ne)**2) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvolh
#else /* NDIMS == 3 */
inarr(5) = inarr(5) + sum((pdata%q(ivx,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivy,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivz,nb:ne,nb:ne, : )**2) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvolh
inarr(5) = inarr(5) + accsum((pdata%q(ivx,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivy,nb:ne,nb:ne, : )**2 &
+ pdata%q(ivz,nb:ne,nb:ne, : )**2) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvolh
#endif /* NDIMS == 3 */
! sum up internal energy
!
if (ipr > 0) then
#if NDIMS == 3
inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol
inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
inarr(9) = inarr(9) + sum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol
inarr(9) = inarr(9) + accsum(pdata%q(ipr,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
end if
@ -808,13 +808,13 @@ module statistics
!
if (magnetized) then
#if NDIMS == 3
inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh
inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(iby,nb:ne,nb:ne,nb:ne)**2 &
+ pdata%q(ibz,nb:ne,nb:ne,nb:ne)**2) * dvolh
#else /* NDIMS == 3 */
inarr(13) = inarr(13) + sum(pdata%q(ibx,nb:ne,nb:ne, : )**2 &
+ pdata%q(iby,nb:ne,nb:ne, : )**2 &
+ pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh
inarr(13) = inarr(13) + accsum(pdata%q(ibx,nb:ne,nb:ne, : )**2 &
+ pdata%q(iby,nb:ne,nb:ne, : )**2 &
+ pdata%q(ibz,nb:ne,nb:ne, : )**2) * dvolh
#endif /* NDIMS == 3 */
! calculate current density (J = ∇xB)
@ -1036,13 +1036,13 @@ module statistics
!
inarr(24) = inarr(24) &
#if NDIMS == 3
- sum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
- accsum((pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibz,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
* pdata%q(iby,nb:ne,nb:ne,nb:ne)) &
* jc(1,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
- sum((pdata%q(ivy,nb:ne,nb:ne, : ) &
- accsum((pdata%q(ivy,nb:ne,nb:ne, : ) &
* pdata%q(ibz,nb:ne,nb:ne, : ) &
- pdata%q(ivz,nb:ne,nb:ne, : ) &
* pdata%q(iby,nb:ne,nb:ne, : )) &
@ -1050,13 +1050,13 @@ module statistics
#endif /* NDIMS == 3 */
inarr(24) = inarr(24) &
#if NDIMS == 3
- sum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
- accsum((pdata%q(ivz,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibx,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibz,nb:ne,nb:ne,nb:ne)) &
* jc(2,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
- sum((pdata%q(ivz,nb:ne,nb:ne, : ) &
- accsum((pdata%q(ivz,nb:ne,nb:ne, : ) &
* pdata%q(ibx,nb:ne,nb:ne, : ) &
- pdata%q(ivx,nb:ne,nb:ne, : ) &
* pdata%q(ibz,nb:ne,nb:ne, : )) &
@ -1064,13 +1064,13 @@ module statistics
#endif /* NDIMS == 3 */
inarr(24) = inarr(24) &
#if NDIMS == 3
- sum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
- accsum((pdata%q(ivx,nb:ne,nb:ne,nb:ne) &
* pdata%q(iby,nb:ne,nb:ne,nb:ne) &
- pdata%q(ivy,nb:ne,nb:ne,nb:ne) &
* pdata%q(ibx,nb:ne,nb:ne,nb:ne)) &
* jc(3,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
- sum((pdata%q(ivx,nb:ne,nb:ne, : ) &
- accsum((pdata%q(ivx,nb:ne,nb:ne, : ) &
* pdata%q(iby,nb:ne,nb:ne, : ) &
- pdata%q(ivy,nb:ne,nb:ne, : ) &
* pdata%q(ibx,nb:ne,nb:ne, : )) &
@ -1082,9 +1082,9 @@ module statistics
if (resistivity > 0.0d+00) then
inarr(26) = inarr(26) &
#if NDIMS == 3
- sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2) * dvol
- accsum(sum(jc(1:3,nb:ne,nb:ne,nb:ne)**2, 1)) * dvol
#else /* NDIMS == 3 */
- sum(jc(1:3,nb:ne,nb:ne, : )**2) * dvol
- accsum(sum(jc(1:3,nb:ne,nb:ne, : )**2, 1)) * dvol
#endif /* NDIMS == 3 */
end if ! resistive
@ -1098,10 +1098,10 @@ module statistics
inarr(27) = inarr(27) &
#if NDIMS == 3
- sum(jc(1,nb:ne,nb:ne,nb:ne) &
- accsum(jc(1,nb:ne,nb:ne,nb:ne) &
* jc(2,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
- sum(jc(1,nb:ne,nb:ne, : ) &
- accsum(jc(1,nb:ne,nb:ne, : ) &
* jc(2,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
@ -1111,11 +1111,11 @@ module statistics
inarr(28) = inarr(28) &
#if NDIMS == 3
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol
- accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol
#else /* NDIMS == 3 */
- sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : )) * dvol
- accsum(sum(pdata%q(ibx:ibz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol
#endif /* NDIMS == 3 */
end if ! magnetized
@ -1132,13 +1132,13 @@ module statistics
inarr(25) = inarr(25) &
#if NDIMS == 3
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
+ accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvol
+ accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) * dvol
#endif /* NDIMS == 3 */
! calculate divergence (∇·v)
@ -1151,15 +1151,15 @@ module statistics
inarr(25) = inarr(25) &
#if NDIMS == 3
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
* dvol / 3.0d+00
+ accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne),1) &
* pdata%q(idn,nb:ne,nb:ne,nb:ne)) &
* dvol / 3.0d+00
#else /* NDIMS == 3 */
+ sum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) &
* dvol / 3.0d+00
+ accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ),1) &
* pdata%q(idn,nb:ne,nb:ne, : )) &
* dvol / 3.0d+00
#endif /* NDIMS == 3 */
end if
@ -1175,11 +1175,11 @@ module statistics
inarr(23) = inarr(23) &
#if NDIMS == 3
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne)) * dvol
- accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne,nb:ne) &
* jc(1:3,nb:ne,nb:ne,nb:ne), 1)) * dvol
#else /* NDIMS == 3 */
- sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : )) * dvol
- accsum(sum(pdata%q(ivx:ivz,nb:ne,nb:ne, : ) &
* jc(1:3,nb:ne,nb:ne, : ), 1)) * dvol
#endif /* NDIMS == 3 */
else
@ -1282,6 +1282,52 @@ module statistics
!-------------------------------------------------------------------------------
!
end subroutine store_statistics
!
!===============================================================================
!
! function ACCSUM:
! ---------------
!
! Function performs a more accurate Kahan summation of the input array
! elements.
!
!===============================================================================
!
real(kind=8) function accsum(q) result(s)
implicit none
real(kind=8), dimension(:,:,:), intent(in) :: q
integer :: i, j, k
real(kind=8) :: c, y, t
!-------------------------------------------------------------------------------
!
s = 0.0d+00
c = 0.0d+00
do k = 1, size(q, 3)
do j = 1, size(q, 2)
do i = 1, size(q, 1)
t = s + q(i,j,k)
if (abs(s) >= abs(q(i,j,k))) then
c = c + (s - t) + q(i,j,k)
else
c = c + (q(i,j,k) - t) + s
end if
s = t
end do
end do
end do
s = s + c
return
!-------------------------------------------------------------------------------
!
end function accsum
!===============================================================================
!