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

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2021-11-12 14:14:21 -03:00
commit aec1d3a0e1
1 changed files with 82 additions and 42 deletions

124
sources/forcing.F90

@ -1525,47 +1525,80 @@ module forcing
!
subroutine inject_fmodes_block(pdata, dt)
! include external variables
!
use blocks , only : block_data
use constants , only : pi2
use coordinates, only : nm => bcells, nb, ne
use coordinates, only : ax, ay, advol
use blocks , only : block_data
use constants , only : pi2
use coordinates , only : nn => bcells, nb, ne
use coordinates , only : ax, ay, advol
#if NDIMS == 3
use coordinates, only : az
use coordinates , only : az
#endif /* NDIMS == 3 */
use equations , only : idn, imx, imy, imz, ien
use equations , only : idn, imx, imy, imz, ien
use iso_fortran_env, only : error_unit
use mesh , only : work, nwork, work_in_use
! local variables are not implicit by default
!
implicit none
! subroutine arguments
!
type(block_data), pointer, intent(inout) :: pdata
real(kind=8) , intent(in) :: dt
! local variables
!
integer :: i, j, k = 1, l, n
logical, save :: first = .true.
integer :: i, j, k = 1, l, n, status
real(kind=8) :: cs, sn
#if NDIMS == 3
real(kind=8) :: tt
#endif /* NDIMS == 3 */
real(kind=8) :: dvol
! local arrays
!
real(kind=8), dimension(nm):: x, y, z
real(kind=8), dimension(nm):: kx, ky, snkx, snky, cskx, csky
real(kind=8), dimension(nn):: x, y, z
real(kind=8), dimension(nn):: kx, ky, snkx, snky, cskx, csky
#if NDIMS == 3
real(kind=8), dimension(nm):: kz, snkz, cskz
real(kind=8), dimension(nn):: kz, snkz, cskz
#endif /* NDIMS == 3 */
!
real(kind=8), dimension(:,:,:,:), pointer, save :: acc
real(kind=8), dimension(:,:,:) , pointer, save :: den
character(len=*), parameter :: loc = 'FORCING:inject_fmodes_block()'
!-------------------------------------------------------------------------------
!
! prepare block coordinates
!
if (first) then
status = 0
if (work_in_use) then
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), &
"Workspace is already occupied!"
status = 1
go to 100
else
i = 3 * nn**NDIMS
j = 4 * nn**NDIMS
if (j > nwork) then
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), &
"Workspace has too be increased!"
nwork = j
deallocate(work)
allocate(work(nwork), stat=status)
if (status /= 0) then
write(error_unit,"('[',a,']: ',a)") trim(loc), &
"Cannot increase the workspace's size!"
go to 100
end if
end if
#if NDIMS == 3
acc(1:3,1:nn,1:nn,1:nn) => work( 1:i)
den(1:nn,1:nn,1:nn) => work(i+1:j)
#else /* NDIMS == 3 */
acc(1:2,1:nn,1:nn,1: 1) => work(1:i)
den(1:nn,1:nn,1: 1) => work(i+1:j)
#endif /* NDIMS == 3 */
end if
first = .false.
end if
x(:) = - pi2 * (pdata%meta%xmin + ax(pdata%meta%level,:))
y(:) = - pi2 * (pdata%meta%ymin + ay(pdata%meta%level,:))
#if NDIMS == 3
@ -1575,9 +1608,14 @@ module forcing
#endif /* NDIMS == 3 */
dvol = advol(pdata%meta%level)
! reset the acceleration
!
pdata%du(1:NDIMS,:,:,:) = 0.0d+00
if (work_in_use) then
write(error_unit,"('[',a,']: ',a,3i4,a)") trim(loc), &
"Workspace is already occupied! Corruptions can occur!"
end if
work_in_use = .true.
acc(1:NDIMS,:,:,:) = 0.0d+00
! iterate over driving modes and compute the acceleration in the real space
!
@ -1600,10 +1638,10 @@ module forcing
#endif /* NDIMS == 3 */
#if NDIMS == 3
do k = 1, nm
do k = 1, nn
#endif /* NDIMS == 3 */
do j = 1, nm
do i = 1, nm
do j = 1, nn
do i = 1, nn
cs = cskx(i) * csky(j) - snkx(i) * snky(j)
sn = cskx(i) * snky(j) + snkx(i) * csky(j)
@ -1613,9 +1651,9 @@ module forcing
sn = tt * snkz(k) + sn * cskz(k)
#endif /* NDIMS == 3 */
pdata%du(1:NDIMS,i,j,k) = pdata%du(1:NDIMS,i,j,k) &
+ (real(fcoefs(l,:)) * cs &
+ aimag(fcoefs(l,:)) * sn)
acc(1:NDIMS,i,j,k) = acc(1:NDIMS,i,j,k) &
+ (real(fcoefs(l,:)) * cs &
+ aimag(fcoefs(l,:)) * sn)
end do
end do
#if NDIMS == 3
@ -1629,40 +1667,42 @@ module forcing
! multiply the acceleration by density and time step
!
do n = 1, NDIMS
pdata%du(n,:,:,:) = pdata%u(idn,:,:,:) * pdata%du(n,:,:,:) * dt
acc(n,:,:,:) = pdata%u(idn,:,:,:) * acc(n,:,:,:) * dt
end do
! calculate the kinetic energy before the injection
!
pdata%du(4,:,:,:) = sum(pdata%u(imx:imz,:,:,:)**2, 1)
den(:,:,:) = sum(pdata%u(imx:imz,:,:,:)**2, 1)
! add the momentum increment
!
pdata%u(imx,:,:,:) = pdata%u(imx,:,:,:) + pdata%du(1,:,:,:)
pdata%u(imy,:,:,:) = pdata%u(imy,:,:,:) + pdata%du(2,:,:,:)
pdata%u(imx,:,:,:) = pdata%u(imx,:,:,:) + acc(1,:,:,:)
pdata%u(imy,:,:,:) = pdata%u(imy,:,:,:) + acc(2,:,:,:)
#if NDIMS == 3
pdata%u(imz,:,:,:) = pdata%u(imz,:,:,:) + pdata%du(3,:,:,:)
pdata%u(imz,:,:,:) = pdata%u(imz,:,:,:) + acc(3,:,:,:)
#endif /* NDIMS == 3 */
! calculate the injected energy
!
pdata%du(4,:,:,:) = 5.0d-01 * (sum(pdata%u(imx:imz,:,:,:)**2, 1) &
- pdata%du(4,:,:,:)) / pdata%u(idn,:,:,:)
den(:,:,:) = 5.0d-01 * (sum(pdata%u(imx:imz,:,:,:)**2, 1) &
- den(:,:,:)) / pdata%u(idn,:,:,:)
! calculate the injected energy
!
#if NDIMS == 3
einj = einj + sum(pdata%du(4,nb:ne,nb:ne,nb:ne)) * dvol
einj = einj + sum(den(nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
einj = einj + sum(pdata%du(4,nb:ne,nb:ne, 1 )) * dvol
einj = einj + sum(den(nb:ne,nb:ne, 1 )) * dvol
#endif /* NDIMS == 3 */
! add the energy increment
!
if (ien > 0) then
pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + pdata%du(4,:,:,:)
pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:)
end if
100 continue
!-------------------------------------------------------------------------------
!
end subroutine inject_fmodes_block