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

FORCING: Allow to specify the period of turbulence injection.

Browse Source 
The period is controlled by two parameters:

'injection_time_start' - the beginning of the injection period,
'injection_time_stop'  - the end of the injection period.

Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
Grzegorz Kowal 2023-02-27 10:04:06 -03:00
parent 5b26a7f726
commit 99d5196736
2 changed files with 41 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
sources/evolution.F90
View File

@ -1109,7 +1109,7 @@ module evolution
! add forcing contribution
!
if (forcing_enabled) call update_forcing(dt)
if (forcing_enabled) call update_forcing(time, dt)
! check if we need to perform the refinement step
!

51
sources/forcing.F90
View File

@ -35,8 +35,8 @@ module forcing
! interfaces for procedure pointers
!
abstract interface
subroutine update_forcing_iface(dt)
real(kind=8), intent(in) :: dt
subroutine update_forcing_iface(t, dt)
real(kind=8), intent(in) :: t, dt
end subroutine
end interface
@ -89,6 +89,8 @@ module forcing
real(kind=8), save :: zth = 0.0d+00
real(kind=8), save :: zfc = 1.0d+00
#endif /* NDIMS == 3 */
real(kind=8) :: tinj_start = 0.00000d+00
real(kind=8) :: tinj_stop = 9.99999d+99
! Fourier profile parameters
!
@ -244,6 +246,21 @@ module forcing
end if
#endif /* NDIMS == 3 */
! get the injection period
!
call get_parameter("injection_time_start", tinj_start)
call get_parameter("injection_time_stop" , tinj_stop )
if (tinj_start >= tinj_stop) then
if (verbose) then
write(*,*)
write(*,wfmt) "Injection cannot stop before it starts!"
write(*,wfmt) "Please, modify parameter 'injection_time_start'" // &
" or 'injection_time_stop'."
write(*,wfmt) "Resetting 'injection_time_start' to 0.0!"
tinj_start = 0.0d+00
end if
end if
! select the energy injection method
!
select case(trim(injection))
@ -929,6 +946,9 @@ module forcing
call print_parameter(verbose, "amplitude threshold", fmin)
end if
call print_parameter(verbose, "injection starts at", tinj_start)
call print_parameter(verbose, "injection stops at", tinj_stop)
if (.not. periodic(1)) &
call print_parameter(verbose, "profile X transition between", xlo, xup)
if (.not. periodic(2)) &
@ -953,18 +973,19 @@ module forcing
!
! Arguments:
!
! dt - time increment;
! t - time;
! dt - time step;
!
!===============================================================================
!
subroutine update_forcing_eddies(dt)
subroutine update_forcing_eddies(t, dt)
use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen
use random , only : randuni, randsym
implicit none
real(kind=8), intent(in) :: dt
real(kind=8), intent(in) :: t, dt
integer :: ni, n
#if NDIMS == 3
@ -974,6 +995,8 @@ module forcing
!-------------------------------------------------------------------------------
!
if (t < tinj_start .or. t > tinj_stop) return
! calculate the number of eddies to be injected during this interval
!
dinj = dinj + rate * dt
@ -1029,11 +1052,12 @@ module forcing
!
! Arguments:
!
! dt - time increment;
! t - time;
! dt - time step;
!
!===============================================================================
!
subroutine update_forcing_oh(dt)
subroutine update_forcing_oh(t, dt)
#if NDIMS == 3
use constants, only : pi2
@ -1042,7 +1066,7 @@ module forcing
implicit none
real(kind=8), intent(in) :: dt
real(kind=8), intent(in) :: t, dt
integer :: l
real(kind=8) :: acoeff, dcoeff
@ -1055,6 +1079,8 @@ module forcing
!-------------------------------------------------------------------------------
!
if (t < tinj_start .or. t > tinj_stop) return
! calculate drift and diffusion coefficients
!
acoeff = exp(- dt / tscale)
@ -1113,18 +1139,19 @@ module forcing
!
! Arguments:
!
! dt - time increment;
! t - time;
! dt - time step;
!
!===============================================================================
!
subroutine update_forcing_alvelius(dt)
subroutine update_forcing_alvelius(t, dt)
use constants, only : pi2
use random , only : randuni
implicit none
real(kind=8), intent(in) :: dt
real(kind=8), intent(in) :: t, dt
integer :: l
real(kind=8) :: th1, dinj, sqdt
@ -1139,6 +1166,8 @@ module forcing
!-------------------------------------------------------------------------------
!
if (t < tinj_start .or. t > tinj_stop) return
! determine velocify coefficients in Fourier space
!
call get_vcoefs()