!!******************************************************************************
!!
!! This file is part of the AMUN source code, a program to perform
!! Newtonian or relativistic magnetohydrodynamical simulations on uniform or
!! adaptive mesh.
!!
!! Copyright (C) 2017-2021 Grzegorz Kowal <grzegorz@amuncode.org>
!!
!! This program is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!*****************************************************************************
!!
!! module: FORCING
!!
!! This module provides energy injection, e.g. turbulence driving, supernova
!! explosions, etc.
!!
!!*****************************************************************************
!
module forcing
implicit none
! interfaces for procedure pointers
!
abstract interface
subroutine update_forcing_iface(dt)
real(kind=8), intent(in) :: dt
end subroutine
end interface
! pointers to the forcing methods
!
procedure(update_forcing_iface), pointer, save :: update_forcing => null()
! flag indicating if the energy injection is enabled and if it is done in
! the Fourier space
!
logical, save :: forcing_enabled = .false.
logical, save :: forcing_fourier = .false.
logical, save :: energy_profile = .false.
! implemented injection methods
!
integer, parameter :: injection_none = 0
integer, parameter :: injection_eddy = 1
integer, parameter :: injection_oh = 2
integer, parameter :: injection_alvelius = 3
! implemented driving force spectral profiles
!
integer, parameter :: profile_normal = 1
integer, parameter :: profile_parabolic = 2
integer, parameter :: profile_powerlaw = 3
! driving force method and spectral profile
!
integer, save :: injection_method = injection_none
integer, save :: injection_profile = profile_normal
! forcing parameters
!
real(kind=8), save :: power = 1.0d+00
real(kind=8), save :: rate = 1.0d+03
real(kind=8), save :: amp = 1.0d-03
real(kind=8), save :: del = 1.0d-01
! Fourier profile parameters
!
integer , save :: nmodes = 0
real(kind=8), save :: fpower = 1.0d+00
real(kind=8), save :: lscale = 5.0d-01
real(kind=8), save :: tscale = 1.0d+00
real(kind=8), save :: vscale = 5.0d-01
real(kind=8), save :: ascale = 5.0d-01
real(kind=8), save :: kf = 2.0d+00
real(kind=8), save :: kl = 0.0d+00
real(kind=8), save :: ku = 4.0d+00
real(kind=8), save :: kc = 5.0d-01
real(kind=8), save :: pidx = - 5.0d+00 / 3.0d+00
real(kind=8), save :: fmin = 1.0d-08
! module parameters
!
real(kind=8), save :: dinj = 0.0d+00
real(kind=8), save :: einj = 0.0d+00
real(kind=8), save :: rinj = 0.0d+00
real(kind=8), save :: arms = 0.0d+00
! bound for the driving mode wave number
!
integer, save :: kmax = 4
! arrays for driving mode positions and amplitudes
!
integer , dimension(:,: ) , allocatable :: kmodes
real(kind=8), dimension(:) , allocatable :: fmodes
! vector of driving directions in Alfvelius method
!
real(kind=8), dimension(:,:) , allocatable :: e1vecs, e2vecs
! array for driving mode coefficients
!
complex(kind=8), dimension(:,:), allocatable, target :: fcoefs
! velocity Fourier coefficients in Alfvelius method
!
complex(kind=8), dimension(:,:), allocatable :: vcoefs
! by default everything is private
!
private
! declare public subroutines
!
public :: initialize_forcing, finalize_forcing, print_forcing
public :: update_forcing
public :: forcing_enabled, einj, rinj, arms, nmodes, fcoefs
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
contains
!
!===============================================================================
!!
!!*** PUBLIC SUBROUTINES *****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INITIALIZE_FORCING:
! -----------------------------
!
! Subroutine initializes module FORCING.
!
! Arguments:
!
! verbose - flag determining if the subroutine should be verbose;
! status - return flag of the procedure execution status;
!
!===============================================================================
!
subroutine initialize_forcing(verbose, status)
#if NDIMS == 3
use constants , only : pi2
#endif /* NDIMS == 3 */
use helpers , only : print_message
use parameters, only : get_parameter
use random , only : randuni, randnorz
implicit none
logical, intent(in) :: verbose
integer, intent(out) :: status
character(len=32) :: wfmt = "(1x,'WARNING:',1x,a)"
character(len=64) :: injection = "none"
character(len=64) :: profile_type = "gauss"
character(len=64) :: profile_energy = "off"
integer :: i, j, l, k2
#if NDIMS == 3
integer :: k = 0
#endif /* NDIMS == 3 */
real(kind=8) :: kl2, ku2, kv2, kv
real(kind=8) :: fa, fi, uu
#if NDIMS == 3
real(kind=8) :: phi
#endif /* NDIMS == 3 */
real(kind=8), dimension(NDIMS) :: u, v
integer, dimension(:), allocatable :: kcount
character(len=*), parameter :: loc = 'FORCING::initialize_forcing()'
!-------------------------------------------------------------------------------
!
status = 0
! obtain the chosen injection method
!
call get_parameter("driving_method", injection)
! select the energy injection method
!
select case(trim(injection))
case("eddies", "eddy")
forcing_enabled = .true.
injection_method = injection_eddy
update_forcing => update_forcing_eddies
call get_parameter("injection_power", power)
call get_parameter("injection_rate" , rate )
call get_parameter("eddy_amplitude" , amp )
call get_parameter("eddy_width" , del )
case("Ornstein–Uhlenbeck", "ornstein–uhlenbeck", "oh")
forcing_enabled = .true.
forcing_fourier = .true.
injection_method = injection_oh
update_forcing => update_forcing_oh
! get the driving power
!
call get_parameter("driving_timescale", tscale)
if (tscale <= 0.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_timescale' must be larger than zero!"
write(*,wfmt) "Resetting it to the default value: 1.0!"
end if
tscale = 1.0d+00
end if
case("Alvelius", "alvelius", "al")
forcing_enabled = .true.
forcing_fourier = .true.
injection_method = injection_alvelius
update_forcing => update_forcing_alvelius
! get the driving power
!
call get_parameter("driving_power", fpower)
if (fpower <= 0.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_power' must be larger than zero!"
write(*,wfmt) "Resetting it to the default value: 1.0!"
end if
fpower = 1.0d+00
end if
case default
injection_method = injection_none
if (verbose .and. injection /= "none") then
write(*,*)
write(*,wfmt) "Unknown injection method! Forcing disabled!"
end if
end select
! initialization of the Fourier space driving
!
if (forcing_fourier) then
call get_parameter("driving_wavenumber" , kf)
call get_parameter("driving_profile_width", kc)
call get_parameter("driving_profile_lower", kl)
call get_parameter("driving_profile_upper", ku)
call get_parameter("amplitude_threshold" , fmin)
if (kf < 1.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_wavenumber' must be larger than 1.0!"
write(*,wfmt) "Resetting it to the default value: 2.0!"
end if
kf = 2.0d+00
end if
if (kc <= 0.0d+00) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_profile_width' must be larger than 0.0!"
write(*,wfmt) "Resetting it to the default value: 1.0!"
end if
kc = 1.0d+00
end if
if (kc > kf) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_profile_width' must be smaller than " // &
"'driving_wavenumber'!"
write(*,wfmt) "Resetting it to the half of 'driving_wavenumber'!"
end if
kc = 5.0d-01 * kf
end if
if (kl >= kf) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_profile_lower' must be smaller than " // &
"'driving_wavenumber'!"
write(*,wfmt) "Setting it to 'driving_wavenumber' " // &
"- 'driving_profile_width'!"
end if
kl = kf - kc
end if
if (ku <= kf) then
if (verbose) then
write(*,*)
write(*,wfmt) "'driving_profile_upper' must be larger than " // &
"'driving_wavenumber'!"
write(*,wfmt) "Setting it to 'driving_wavenumber' " // &
"+ 'driving_profile_width'!"
end if
ku = kf + kc
end if
fmin = max(1.0d-08, fmin)
kl2 = kl**2
ku2 = ku**2
! get the upper bound for the wave number
!
kmax = ceiling(ku)
! allocate arrays to count the modes with the same wave number
!
allocate(kcount(NDIMS * kmax**2), stat = status)
if (status == 0) then
! initialize mode counts
!
kcount(:) = 0
! count the modes with the same k, get the minimum sufficient number of
! driving modes
!
#if NDIMS == 3
do i = 1, kmax
do j = -kmax, kmax
do k = -kmax, kmax
k2 = i**2 + j**2 + k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
nmodes = nmodes + 1
kcount(k2) = kcount(k2) + 2
end if
end do
end do
end do
do j = 1, kmax
do k = -kmax, kmax
k2 = j**2 + k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
nmodes = nmodes + 1
kcount(k2) = kcount(k2) + 2
end if
end do
end do
do k = 1, kmax
k2 = k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
nmodes = nmodes + 1
kcount(k2) = kcount(k2) + 2
end if
end do
#else /* NDIMS == 3 */
do i = 1, kmax
do j = -kmax, kmax
k2 = i**2 + j**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
nmodes = nmodes + 1
kcount(k2) = kcount(k2) + 2
end if
end do
end do
do j = 1, kmax
k2 = j**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
nmodes = nmodes + 1
kcount(k2) = kcount(k2) + 2
end if
end do
#endif /* NDIMS == 3 */
! allocate arrays for driving modes (positions, amplitudes and coefficients)
!
allocate(e1vecs(nmodes,NDIMS), e2vecs(nmodes,NDIMS), &
kmodes(nmodes,NDIMS), fmodes(nmodes), &
fcoefs(nmodes,NDIMS), stat = status)
if (status == 0) then
! prepare wave vectors
!
l = 0
#if NDIMS == 3
do i = 1, kmax
do j = -kmax, kmax
do k = -kmax, kmax
k2 = i**2 + j**2 + k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
l = l + 1
kmodes(l,1) = i
kmodes(l,2) = j
kmodes(l,3) = k
end if
end do
end do
end do
do j = 1, kmax
do k = -kmax, kmax
k2 = j**2 + k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
l = l + 1
kmodes(l,1) = 0
kmodes(l,2) = j
kmodes(l,3) = k
end if
end do
end do
do k = 1, kmax
k2 = k**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
l = l + 1
kmodes(l,1) = 0
kmodes(l,2) = 0
kmodes(l,3) = k
end if
end do
#else /* NDIMS == 3 */
do i = 1, kmax
do j = -kmax, kmax
k2 = i**2 + j**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
l = l + 1
kmodes(l,1) = i
kmodes(l,2) = j
end if
end do
end do
do j = 1, kmax
k2 = j**2
kv2 = 1.0d+00 * k2
if (kl2 <= kv2 .and. kv2 <= ku2) then
l = l + 1
kmodes(l,1) = 0
kmodes(l,2) = j
end if
end do
#endif /* NDIMS == 3 */
! generate random vectors, perpendicular to the wave vector
!
do l = 1, nmodes
u(:) = kmodes(l,:)
#if NDIMS == 3
if (abs(kmodes(l,3)) < abs(kmodes(l,1))) then
v(:) = (/ kmodes(l,2), - kmodes(l,1), 0 /)
else
v(:) = (/ 0, - kmodes(l,3), kmodes(l,2) /)
end if
#else /* NDIMS == 3 */
v(:) = (/ kmodes(l,2), - kmodes(l,1) /)
#endif /* NDIMS == 3 */
uu = dot_product(v(:), v(:))
if (uu > 0.0d+00) then
e1vecs(l,:) = v(:)
#if NDIMS == 3
e2vecs(l,1) = u(2) * v(3) - u(3) * v(2)
e2vecs(l,2) = u(3) * v(1) - u(1) * v(3)
e2vecs(l,3) = u(1) * v(2) - u(2) * v(1)
#endif /* NDIMS == 3 */
e1vecs(l,:) = e1vecs(l,:) / sqrt(sum(e1vecs(l,:)**2))
#if NDIMS == 3
e2vecs(l,:) = e2vecs(l,:) / sqrt(sum(e2vecs(l,:)**2))
#else /* NDIMS == 3 */
e2vecs(l,:) = 0.0d+00
#endif /* NDIMS == 3 */
else
call print_message(loc, "v(:) is ZERO!")
end if
end do
! determine if the profile is set by energy or amplitude
!
call get_parameter("profile_by_energy", profile_energy)
select case(trim(profile_energy))
case ("on", "ON", "t", "T", "y", "Y", "true", "TRUE", "yes", "YES")
energy_profile = .true.
case default
energy_profile = .false.
kcount(:) = 1
end select
! initialize driving mode amplitudes depending on the chosen profile
!
call get_parameter("driving_profile", profile_type)
select case(trim(profile_type))
! --- power-law profile ---
!
case("power-law", "powerlaw", "band")
injection_profile = profile_powerlaw
! get the power index of the power-law distribution
!
call get_parameter("powerlaw_index", pidx)
! calculate amplitudes of driving modes
!
fi = 0.0d+00
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
kv = sqrt(1.0d+00 * k2)
fa = kv**pidx / kcount(k2)
fi = fi + fa
fmodes(l) = sqrt(fa)
end do
fmodes(:) = fmodes(:) / sqrt(fi)
! --- normal distribution profile ---
!
case("gauss", "gaussian", "normal")
injection_profile = profile_normal
! calculate amplitudes of driving modes
!
fi = 0.0d+00
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
kv = sqrt(1.0d+00 * k2)
fa = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2)
fi = fi + fa
fmodes(l) = sqrt(fa)
end do
fmodes(:) = fmodes(:) / sqrt(fi)
! --- parabolic profile ---
!
case("parabolic")
injection_profile = profile_parabolic
! correct the injection scale and width
!
kf = 0.5d+00 * (ku + kl)
kc = 0.5d+00 * (ku - kl)
! calculate amplitudes of driving modes
!
fi = 0.0d+00
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
kv = sqrt(1.0d+00 * k2)
fa = 4.0d+00 * (kv - kl) * (ku - kv) / (ku - kl)**2
fa = max(0.0d+00, fa)**2 / kcount(k2)
fi = fi + fa
fmodes(l) = sqrt(fa)
end do
fmodes(:) = fmodes(:) / sqrt(fi)
! --- normal distribution profile ---
!
case default
! warn if the specified profile is not implemented and the normal one is used
!
if (verbose .and. .not. (profile_type == "normal" .or. &
profile_type == "gauss" .or. &
profile_type == "gaussian")) then
write(*,*)
write(*,wfmt) "Unknown spectral profile of driving!"
write(*,wfmt) "Available profiles: " // &
"'power-law', 'normal', 'parabolic'."
write(*,wfmt) "Using default normal distribution profile " // &
"'gaussian'."
end if
injection_profile = profile_normal
! calculate amplitudes of driving modes
!
fi = 0.0d+00
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
kv = sqrt(1.0d+00 * k2)
fa = exp(- 0.5d+00 * ((kv - kf) / kc)**2) / kcount(k2)
fi = fi + fa
fmodes(l) = sqrt(fa)
end do
do l = 1, nmodes
k2 = sum(kmodes(l,:)**2)
fmodes(l) = fmodes(l) / sqrt(fi)
end do
end select
! === Ornstein–Uhlenbeck driving ===
!
if (injection_method == injection_oh) then
! get the characteristic driving length, velocity, acceleration scales, and
! determine the estimated driving power (since the main injection is through
! the forcing-velocity correlation, it may not represent the actual
! instantenous injection rate)
!
lscale = 1.0d+00 / kf
vscale = lscale / tscale
ascale = sqrt(2.0d+00) * lscale / tscale**2 / sqrt(tscale)
fpower = vscale**2 / tscale
! normalize the driving amplitude profile to the driving power
!
fmodes(:) = ascale * fmodes(:)
! generate the initial driving mode coefficients
!
do l = 1, nmodes
#if NDIMS == 3
phi = pi2 * randuni()
fcoefs(l,:) = fmodes(l) * (e1vecs(l,:) * cos(phi) &
+ e2vecs(l,:) * sin(phi)) * randnorz()
#else /* NDIMS == 3 */
fcoefs(l,:) = fmodes(l) * e1vecs(l,:) * randnorz()
#endif /* NDIMS == 3 */
end do ! l = 1, nmodes
! calculate the rms of acceleration
!
arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2 &
+ aimag(fcoefs(:,:))**2))
end if ! injection = 'OH'
! === Alvelius driving ===
!
if (injection_method == injection_alvelius) then
! get the characteristic driving length, velocity, time, and acceleration scales
!
lscale = 1.0d+00 / kf
tscale = sqrt(lscale / sqrt(2.0d+00 * fpower))
vscale = lscale / tscale
ascale = vscale / tscale
! normalize the driving amplitude profile to the driving power
!
fmodes(:) = sqrt(2.0d+00) * ascale * fmodes(:)
! reset the driving coefficients and calculate the rms of acceleration
!
fcoefs(:,:) = cmplx(0.0d+00, 0.0d+00, kind=8)
arms = 0.0d+00
! allocate vectors of velocity Fourier coefficients
!
allocate(vcoefs(nmodes,NDIMS), stat = status)
if (status /= 0) &
call print_message(loc, "Cannot allocate vcoefs(:,:)!")
end if ! injection = 'Alvelius'
end if ! allocation of driving mode arrays
end if ! allocation of kcount
! deallocate mode counter
!
if (allocated(kcount)) deallocate(kcount)
end if ! Fourier forcing
!-------------------------------------------------------------------------------
!
end subroutine initialize_forcing
!
!===============================================================================
!
! subroutine FINALIZE_FORCING:
! ---------------------------
!
! Subroutine releases memory used by the module variables.
!
! Arguments:
!
! status - return flag of the procedure execution status;
!
!===============================================================================
!
subroutine finalize_forcing(status)
implicit none
integer, intent(out) :: status
!-------------------------------------------------------------------------------
!
status = 0
if (allocated(kmodes)) deallocate(kmodes)
if (allocated(fmodes)) deallocate(fmodes)
if (allocated(fcoefs)) deallocate(fcoefs)
if (allocated(e1vecs)) deallocate(e1vecs)
if (allocated(e2vecs)) deallocate(e2vecs)
if (allocated(vcoefs)) deallocate(vcoefs)
!-------------------------------------------------------------------------------
!
end subroutine finalize_forcing
!
!===============================================================================
!
! subroutine PRINT_FORCING:
! ------------------------
!
! Subroutine prints module parameters.
!
! Arguments:
!
! verbose - a logical flag turning the information printing;
!
!===============================================================================
!
subroutine print_forcing(verbose)
use helpers, only : print_section, print_parameter
implicit none
logical, intent(in) :: verbose
character(len=64) :: injection = "none"
character(len=64) :: profile_type = "gaussian"
!-------------------------------------------------------------------------------
!
if (verbose) then
select case(injection_method)
case(injection_eddy)
injection = "random eddies"
case(injection_oh)
injection = "Ornstein–Uhlenbeck"
case(injection_alvelius)
injection = "Alvelius"
case default
injection = "none"
end select
call print_section(verbose, "Forcing")
call print_parameter(verbose, "method", injection)
select case(injection_method)
case(injection_eddy)
call print_parameter(verbose, "injection power", power)
call print_parameter(verbose, "injection rate" , rate)
call print_parameter(verbose, "eddy amplitude" , amp)
call print_parameter(verbose, "eddy width" , del)
end select
if (forcing_fourier) then
if (injection_method == injection_oh) then
call print_parameter(verbose, "power (estimated)", fpower)
end if
if (injection_method == injection_alvelius) then
call print_parameter(verbose, "power" , fpower)
end if
call print_parameter(verbose, "length scale" , lscale)
call print_parameter(verbose, "time scale" , tscale)
call print_parameter(verbose, "velocity scale" , vscale)
call print_parameter(verbose, "acceleration scale", ascale)
select case(injection_profile)
case(profile_powerlaw)
profile_type = "power-law"
case(profile_normal)
profile_type = "gaussian"
case(profile_parabolic)
profile_type = "parabolic"
case default
profile_type = "unknown"
end select
call print_parameter(verbose, "spectrum profile", profile_type)
if (energy_profile) then
call print_parameter(verbose, "profile of", "energy")
else
call print_parameter(verbose, "profile of", "amplitude")
end if
call print_parameter(verbose, "number of driving modes", nmodes)
call print_parameter(verbose, "driving wavenumber" , kf)
call print_parameter(verbose, "lower wavenumber limit" , kl)
call print_parameter(verbose, "upper wavenumber limit" , ku)
if (injection_profile == profile_powerlaw) then
call print_parameter(verbose, "power-law index", pidx)
end if
if (injection_profile == profile_normal) then
call print_parameter(verbose, "spectrum profile width", kc)
end if
call print_parameter(verbose, "amplitude threshold", fmin)
end if
end if
!-------------------------------------------------------------------------------
!
end subroutine print_forcing
!
!===============================================================================
!
! subroutine UPDATE_FORCING_EDDIES:
! --------------------------------
!
! Subroutine adds the energy injection terms.
!
! Arguments:
!
! dt - time increment;
!
!===============================================================================
!
subroutine update_forcing_eddies(dt)
use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen
use random , only : randuni, randsym
implicit none
real(kind=8), intent(in) :: dt
integer :: ni, n
#if NDIMS == 3
real(kind=8) :: tmp
#endif /* NDIMS == 3 */
real(kind=8), dimension(3) :: xp, ap
!-------------------------------------------------------------------------------
!
! calculate the number of eddies to be injected during this interval
!
dinj = dinj + rate * dt
ni = int(floor(dinj))
dinj = dinj - 1.0d+00 * ni
! inject the required number of eddies
!
if (ni > 0) then
do n = 1, ni
! get random position
!
xp(1) = xmin + xlen * randuni()
xp(2) = ymin + ylen * randuni()
xp(3) = zmin + zlen * randuni()
! get random orientation
!
#if NDIMS == 3
tmp = 0.0d+00
do while(tmp <= tiny(0.0d+00))
ap(1) = randsym()
ap(2) = randsym()
ap(3) = randsym()
tmp = sqrt(ap(1)**2 + ap(2)**2 + ap(3)**2)
end do
ap(:) = amp * (ap(:) / tmp) / del
#else /* NDIMS == 3 */
ap(:) = sign(1.0d+00, randsym()) * (/ 0.0d+00, 0.0d+00, amp / del /)
#endif /* NDIMS == 3 */
ap(:) = ap(:) * dt
! iterate over data blocks and add forcing components
!
call inject_eddy(xp(:), ap(:))
end do
end if
!-------------------------------------------------------------------------------
!
end subroutine update_forcing_eddies
!
!===============================================================================
!
! subroutine UPDATE_FORCING_OH:
! ----------------------------
!
! Subroutine adds the energy injection terms for Ornstein–Uhlenbeck driving.
!
! Arguments:
!
! dt - time increment;
!
!===============================================================================
!
subroutine update_forcing_oh(dt)
#if NDIMS == 3
use constants, only : pi2
#endif /* NDIMS == 3 */
use random , only : randuni, randnorz
implicit none
real(kind=8), intent(in) :: dt
integer :: l
real(kind=8) :: acoeff, dcoeff
#if NDIMS == 3
real(kind=8) :: phi
#endif /* NDIMS == 3 */
real(kind=8) :: dinj
complex(kind=8), dimension(NDIMS) :: finj
!-------------------------------------------------------------------------------
!
! calculate drift and diffusion coefficients
!
acoeff = exp(- dt / tscale)
dcoeff = sqrt(1.0d+00 - acoeff**2)
! iterate over driving modes
!
do l = 1, nmodes
! generate gaussian random vector along the driving modes
!
#if NDIMS == 3
phi = pi2 * randuni()
finj(:) = fmodes(l) * (e1vecs(l,:) * cos(phi) &
+ e2vecs(l,:) * sin(phi)) * randnorz()
#else /* NDIMS == 3 */
finj(:) = fmodes(l) * e1vecs(l,:) * randnorz()
#endif /* NDIMS == 3 */
! advance the driving force modes
!
fcoefs(l,:) = acoeff * fcoefs(l,:) + dcoeff * finj(:)
end do
! calculate the rms of acceleration
!
arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2 + aimag(fcoefs(:,:))**2))
! store previously injected energy
!
dinj = einj
! inject driving modes
!
call inject_fmodes(dt)
! calculate the energy injected during this step
!
dinj = einj - dinj
! calculate the injection rate
!
rinj = dinj / dt
!-------------------------------------------------------------------------------
!
end subroutine update_forcing_oh
!
!===============================================================================
!
! subroutine UPDATE_FORCING_ALVELIUS:
! ----------------------------------
!
! Subroutine adds the energy injection terms for Alvelius driving.
!
! Arguments:
!
! dt - time increment;
!
!===============================================================================
!
subroutine update_forcing_alvelius(dt)
use constants, only : pi2
use random , only : randuni
implicit none
real(kind=8), intent(in) :: dt
integer :: l
real(kind=8) :: th1, dinj, sqdt
#if NDIMS == 3
real(kind=8) :: th2, phi, psi, ga, gb
#endif /* NDIMS == 3 */
complex(kind=8) :: aran
#if NDIMS == 3
complex(kind=8) :: bran
complex(kind=8) :: xi1, xi2
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
! determine velocify coefficients in Fourier space
!
call get_vcoefs()
! square root of dt
!
sqdt = sqrt(dt)
! iterate over driving modes
!
do l = 1, nmodes
#if NDIMS == 3
! get xi1 and xi2
!
xi1 = dot_product(vcoefs(l,:), e1vecs(l,:))
xi2 = dot_product(vcoefs(l,:), e2vecs(l,:))
! get random angles phi and psi and determine th1 and th2 to reduce
! velocity-force correlation
!
phi = pi2 * randuni()
ga = sin(phi)
gb = cos(phi)
psi = pi2 * randuni()
th1 = - ga * aimag(xi1) + gb * (sin(psi) * real(xi2) &
- cos(psi) * aimag(xi2))
if (abs(th1) > epsilon(0.0d+00)) then
th1 = atan((ga * real(xi1) + gb * (sin(psi) * aimag(xi2) &
+ cos(psi) * real(xi2))) / th1)
end if
th2 = psi + th1
! calculate amplitudes
!
aran = cmplx(cos(th1), sin(th1), kind=8) * ga
bran = cmplx(cos(th2), sin(th2), kind=8) * gb
! calculate driving mode coefficient
!
fcoefs(l,:) = fmodes(l) * (aran * e1vecs(l,:) &
+ bran * e2vecs(l,:)) / sqdt
#else /* NDIMS == 3 */
! get random phase
!
th1 = pi2 * randuni()
aran = cmplx(cos(th1), sin(th1), kind=8)
! calculate driving mode coefficient
!
fcoefs(l,:) = fmodes(l) * aran * e1vecs(l,:) / sqdt
#endif /* NDIMS == 3 */
end do
! calculate the rms of acceleration
!
arms = sqrt(5.0d-01 * sum(real(fcoefs(:,:))**2 + aimag(fcoefs(:,:))**2))
! store previously injected energy
!
dinj = einj
! inject driving modes
!
call inject_fmodes(dt)
! calculate the energy injected during this step
!
dinj = einj - dinj
! calculate the injection rate
!
rinj = dinj / dt
!-------------------------------------------------------------------------------
!
end subroutine update_forcing_alvelius
!
!===============================================================================
!!
!!*** PRIVATE SUBROUTINES ****************************************************
!!
!===============================================================================
!
!===============================================================================
!
! subroutine INJECT_EDDY:
! ----------------------
!
! Subroutine injects one eddy to the domain at given position.
!
! Arguments:
!
! xp, ap - eddy position and amplitude vectors;
!
!===============================================================================
!
subroutine inject_eddy(xp, ap)
use blocks, only : block_data, list_data
implicit none
real(kind=8), dimension(3), intent(in) :: xp, ap
type(block_data), pointer :: pdata
!-------------------------------------------------------------------------------
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! inject eddy into the current block
!
call inject_eddy_block(pdata, xp(:), ap(:))
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!-------------------------------------------------------------------------------
!
end subroutine inject_eddy
!
!===============================================================================
!
! subroutine INJECT_EDDY_BLOCK:
! ----------------------------
!
! Subroutine injects one eddy to the local block.
!
! Arguments:
!
! xp, ap - eddy position and amplitude vectors;
!
!===============================================================================
!
subroutine inject_eddy_block(pdata, xp, ap)
use blocks , only : block_data
use coordinates, only : nm => bcells
use coordinates, only : ax, ay, xlen, ylen
#if NDIMS == 3
use coordinates, only : az, zlen
#endif /* NDIMS == 3 */
use coordinates, only : periodic
use equations , only : idn, imx, imy, imz, ien
implicit none
type(block_data), pointer , intent(inout) :: pdata
real(kind=8), dimension(3), intent(in) :: xp, ap
integer :: i, j, k = 1
real(kind=8) :: x2, y2, r2
#if NDIMS == 3
real(kind=8) :: z2
#endif /* NDIMS == 3 */
real(kind=8) :: fx, fy, fz, fp, e1, e2
real(kind=8), dimension(nm) :: x, y, z
!-------------------------------------------------------------------------------
!
! prepare block coordinates
!
if (periodic(1)) then
fx = 0.5d+00 * xlen
x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1)
x(:) = abs(x(:) / fx)
x(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(x(:))))) * fx / del
else
x(:) = (pdata%meta%xmin + ax(pdata%meta%level,:) - xp(1)) / del
end if
if (periodic(2)) then
fy = 0.5d+00 * ylen
y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2)
y(:) = abs(y(:) / fy)
y(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(y(:))))) * fy / del
else
y(:) = (pdata%meta%ymin + ay(pdata%meta%level,:) - xp(2)) / del
end if
#if NDIMS == 3
if (periodic(3)) then
fz = 0.5d+00 * zlen
z(:) = pdata%meta%zmin + az(pdata%meta%level,:) - xp(3)
z(:) = abs(z(:) / fz)
z(:) = (abs(1.0d+00 - abs(1.0d+00 - abs(z(:))))) * fz / del
else
z(:) = (pdata%meta%zmin + az(pdata%meta%level,:) - xp(3)) / del
end if
#else /* NDIMS == 3 */
z(:) = 0.0d+00
#endif /* NDIMS == 3 */
! iterate over the block coordinates
!
if (ien > 0) then
#if NDIMS == 3
do k = 1, nm
z2 = z(k)**2
#endif /* NDIMS == 3 */
do j = 1, nm
y2 = y(j)**2
do i = 1, nm
x2 = x(i)**2
#if NDIMS == 3
r2 = x2 + y2 + z2
#else /* NDIMS == 3 */
r2 = x2 + y2
#endif /* NDIMS == 3 */
fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2)
#if NDIMS == 3
fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp
fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp
fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp
#else /* NDIMS == 3 */
fx = (- ap(3) * y(j)) * fp
fy = (+ ap(3) * x(i)) * fp
fz = 0.0d+00
#endif /* NDIMS == 3 */
e1 = 0.5d+00 * sum(pdata%u(imx:imx,i,j,k)**2) / pdata%u(idn,i,j,k)
pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx
pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy
pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz
e2 = 0.5d+00 * sum(pdata%u(imx:imz,i,j,k)**2) / pdata%u(idn,i,j,k)
pdata%u(ien,i,j,k) = pdata%u(ien,i,j,k) + (e2 - e1)
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
else
#if NDIMS == 3
do k = 1, nm
z2 = z(k)**2
#endif /* NDIMS == 3 */
do j = 1, nm
y2 = y(j)**2
do i = 1, nm
x2 = x(i)**2
#if NDIMS == 3
r2 = x2 + y2 + z2
#else /* NDIMS == 3 */
r2 = x2 + y2
#endif /* NDIMS == 3 */
fp = pdata%u(idn,i,j,k) * exp(- 0.5d+00 * r2)
#if NDIMS == 3
fx = (- ap(3) * y(j) + ap(2) * z(k)) * fp
fy = (- ap(1) * z(k) + ap(3) * x(i)) * fp
fz = (- ap(2) * x(i) + ap(1) * y(j)) * fp
#else /* NDIMS == 3 */
fx = (- ap(3) * y(j)) * fp
fy = (+ ap(3) * x(i)) * fp
fz = 0.0d+00
#endif /* NDIMS == 3 */
pdata%u(imx,i,j,k) = pdata%u(imx,i,j,k) + fx
pdata%u(imy,i,j,k) = pdata%u(imy,i,j,k) + fy
pdata%u(imz,i,j,k) = pdata%u(imz,i,j,k) + fz
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
end if
!-------------------------------------------------------------------------------
!
end subroutine inject_eddy_block
!
!===============================================================================
!
! subroutine INJECT_FMODES:
! ------------------------
!
! Subroutine injects Fourier driving modes.
!
! Arguments:
!
! dt - the time increment;
!
!===============================================================================
!
subroutine inject_fmodes(dt)
use blocks, only : block_data, list_data
implicit none
real(kind=8), intent(in) :: dt
type(block_data), pointer :: pdata
!-------------------------------------------------------------------------------
!
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! inject eddy into the current block
!
call inject_fmodes_block(pdata, dt)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
!-------------------------------------------------------------------------------
!
end subroutine inject_fmodes
!
!===============================================================================
!
! subroutine INJECT_FMODES_BLOCK:
! ------------------------------
!
! Subroutine injects Fourier driving modes to the local block.
!
! Arguments:
!
! pdata - a pointer to the data block;
! dt - the time increment;
!
!===============================================================================
!
subroutine inject_fmodes_block(pdata, dt)
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
#endif /* NDIMS == 3 */
use equations , only : idn, imx, imy, imz, ien
use helpers , only : print_message
use workspace , only : resize_workspace, work, work_in_use
implicit none
type(block_data), pointer, intent(inout) :: pdata
real(kind=8) , intent(in) :: dt
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
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(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()'
!-------------------------------------------------------------------------------
!
if (first) then
i = 3 * nn**NDIMS
j = 4 * nn**NDIMS
call resize_workspace(j, status)
if (status /= 0) then
call print_message(loc, "Could not resize the workspace!")
go to 100
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 */
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
z(:) = - pi2 * (pdata%meta%zmin + az(pdata%meta%level,:))
#else /* NDIMS == 3 */
z(:) = 0.0d+00
#endif /* NDIMS == 3 */
dvol = advol(pdata%meta%level)
if (work_in_use) &
call print_message(loc, "Workspace is being used right now! " // &
"Corruptions can occur!")
work_in_use = .true.
acc(1:NDIMS,:,:,:) = 0.0d+00
! iterate over driving modes and compute the acceleration in the real space
!
do l = 1, nmodes
if (fmodes(l) > fmin) then
kx(:) = kmodes(l,1) * x(:)
ky(:) = kmodes(l,2) * y(:)
#if NDIMS == 3
kz(:) = kmodes(l,3) * z(:)
#endif /* NDIMS == 3 */
cskx(:) = cos(kx(:))
snkx(:) = sin(kx(:))
csky(:) = cos(ky(:))
snky(:) = sin(ky(:))
#if NDIMS == 3
cskz(:) = cos(kz(:))
snkz(:) = sin(kz(:))
#endif /* NDIMS == 3 */
#if NDIMS == 3
do k = 1, nn
#endif /* NDIMS == 3 */
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)
#if NDIMS == 3
tt = cs
cs = tt * cskz(k) - sn * snkz(k)
sn = tt * snkz(k) + sn * cskz(k)
#endif /* NDIMS == 3 */
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
end do
#endif /* NDIMS == 3 */
end if
end do ! l = 1, nmodes
! multiply the acceleration by density and time step
!
do n = 1, NDIMS
acc(n,:,:,:) = pdata%u(idn,:,:,:) * acc(n,:,:,:) * dt
end do
! calculate the kinetic energy before the injection
!
den(:,:,:) = sum(pdata%u(imx:imz,:,:,:)**2, 1)
! add the momentum increment
!
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,:,:,:) + acc(3,:,:,:)
#endif /* NDIMS == 3 */
! calculate the injected energy
!
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(den(nb:ne,nb:ne,nb:ne)) * dvol
#else /* NDIMS == 3 */
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,:,:,:) + den(:,:,:)
end if
work_in_use = .false.
100 continue
!-------------------------------------------------------------------------------
!
end subroutine inject_fmodes_block
!
!===============================================================================
!
! subroutine GET_VCOEFS:
! ---------------------
!
! Subroutine determines velocity Fourier coefficients.
!
!
!===============================================================================
!
subroutine get_vcoefs()
use blocks , only : block_data, list_data
#ifdef MPI
use mpitools, only : reduce_sum
#endif /* MPI */
implicit none
type(block_data), pointer :: pdata
!-------------------------------------------------------------------------------
!
! reset vcoefs
!
vcoefs(:,:) = cmplx(0.0d+00, 0.0d+00, kind=8)
! assign pdata with the first block on the data block list
!
pdata => list_data
! iterate over all data blocks
!
do while (associated(pdata))
! get contribution of velocity coefficients from the current block
!
call get_vcoefs_block(pdata)
! assign pdata to the next block
!
pdata => pdata%next
end do ! over data blocks
#ifdef MPI
! reduce velocity coefficients over all processes
!
call reduce_sum(vcoefs)
#endif /* MPI */
!-------------------------------------------------------------------------------
!
end subroutine get_vcoefs
!
!===============================================================================
!
! subroutine GET_VCOEFS_BLOCK:
! ---------------------------
!
! Subroutine determines Fourier coefficients of velocity.
!
! Arguments:
!
! pdata - a pointer to the data block;
!
!===============================================================================
!
subroutine get_vcoefs_block(pdata)
use blocks , only : block_data
use constants , only : pi2
use coordinates, only : nm => bcells, nb, ne
use coordinates, only : ax, ay, advol
#if NDIMS == 3
use coordinates, only : az
#endif /* NDIMS == 3 */
use equations , only : ivx, ivy
#if NDIMS == 3
use equations , only : ivz
#endif /* NDIMS == 3 */
implicit none
type(block_data), pointer, intent(inout) :: pdata
integer :: i, j, k = 1, l
real(kind=8) :: cs, sn, dvol
#if NDIMS == 3
real(kind=8) :: tt
#endif /* NDIMS == 3 */
complex(kind=8) :: cf
real(kind=8), dimension(nm):: x, y, z
real(kind=8), dimension(nm):: kx, ky, snkx, snky, cskx, csky
#if NDIMS == 3
real(kind=8), dimension(nm):: kz, snkz, cskz
#endif /* NDIMS == 3 */
!-------------------------------------------------------------------------------
!
! prepare block coordinates
!
x(:) = pi2 * (pdata%meta%xmin + ax(pdata%meta%level,:))
y(:) = pi2 * (pdata%meta%ymin + ay(pdata%meta%level,:))
#if NDIMS == 3
z(:) = pi2 * (pdata%meta%zmin + az(pdata%meta%level,:))
#else /* NDIMS == 3 */
z(:) = 0.0d+00
#endif /* NDIMS == 3 */
dvol = advol(pdata%meta%level)
! iterate over driving modes and compute the velocity Fourier coefficients
!
do l = 1, nmodes
kx(:) = kmodes(l,1) * x(:)
ky(:) = kmodes(l,2) * y(:)
#if NDIMS == 3
kz(:) = kmodes(l,3) * z(:)
#endif /* NDIMS == 3 */
cskx(:) = cos(kx(:))
snkx(:) = sin(kx(:))
csky(:) = cos(ky(:))
snky(:) = sin(ky(:))
#if NDIMS == 3
cskz(:) = cos(kz(:))
snkz(:) = sin(kz(:))
#endif /* NDIMS == 3 */
#if NDIMS == 3
do k = nb, ne
#endif /* NDIMS == 3 */
do j = nb, ne
do i = nb, ne
cs = cskx(i) * csky(j) - snkx(i) * snky(j)
sn = cskx(i) * snky(j) + snkx(i) * csky(j)
#if NDIMS == 3
tt = cs
cs = tt * cskz(k) - sn * snkz(k)
sn = tt * snkz(k) + sn * cskz(k)
#endif /* NDIMS == 3 */
cf = cmplx(cs, sn, kind=8) * dvol
vcoefs(l,1) = vcoefs(l,1) + pdata%q(ivx,i,j,k) * cf
vcoefs(l,2) = vcoefs(l,2) + pdata%q(ivy,i,j,k) * cf
#if NDIMS == 3
vcoefs(l,3) = vcoefs(l,3) + pdata%q(ivz,i,j,k) * cf
#endif /* NDIMS == 3 */
end do
end do
#if NDIMS == 3
end do
#endif /* NDIMS == 3 */
end do ! l = 1, nmodes
!-------------------------------------------------------------------------------
!
end subroutine get_vcoefs_block
!===============================================================================
!
end module forcing