!!****************************************************************************** !! !! 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-2024 Grzegorz Kowal !! !! 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 . !! !!***************************************************************************** !! !! 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(t, dt, injected) real(kind=8), intent(in ) :: t, dt logical , intent( out) :: injected 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 real(kind=8), save :: xlo = 9.9d+99 real(kind=8), save :: xup = 9.9d+99 real(kind=8), save :: xth = 0.0d+00 real(kind=8), save :: xfc = 1.0d+00 real(kind=8), save :: ylo = 9.9d+99 real(kind=8), save :: yup = 9.9d+99 real(kind=8), save :: yth = 0.0d+00 real(kind=8), save :: yfc = 1.0d+00 #if NDIMS == 3 real(kind=8), save :: zlo = 9.9d+99 real(kind=8), save :: zup = 9.9d+99 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 ! 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 = 0.0d+00 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) use constants , only : pi2 use coordinates, only : xlen, ylen, zlen 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 */ integer :: kxs, kxm integer :: kys, kym #if NDIMS == 3 integer :: kzs, kzm #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) ! get the injection profile parameters ! call get_parameter("driving_x_lower", xlo) call get_parameter("driving_x_upper", xup) call get_parameter("driving_y_lower", ylo) call get_parameter("driving_y_upper", yup) #if NDIMS == 3 call get_parameter("driving_z_lower", zlo) call get_parameter("driving_z_upper", zup) #endif /* NDIMS == 3 */ xth = xup - xlo if (xth > 0.0d+00) then xfc = 2.5d-01 * pi2 / xth else xfc = 1.0d+00 end if yth = yup - ylo if (yth > 0.0d+00) then yfc = 2.5d-01 * pi2 / yth else yfc = 1.0d+00 end if #if NDIMS == 3 zth = zup - zlo if (zth > 0.0d+00) then zfc = 2.5d-01 * pi2 / zth else zfc = 1.0d+00 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)) 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 call get_parameter("driving_velocity", vscale) if (vscale < 0.0d+00) then if (verbose) then write(*,*) write(*,wfmt) "'driving_velocity' must be positive!" write(*,wfmt) "Resetting it to the default value: 0.0!" end if vscale = 0.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 call get_parameter("driving_velocity", vscale) if (vscale < 0.0d+00) then if (verbose) then write(*,*) write(*,wfmt) "'driving_velocity' must be positive!" write(*,wfmt) "Resetting it to the default value: 0.0!" end if vscale = 0.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 (kf * min(xlen, ylen, zlen) < 1.0d+00) then if (verbose) then write(*,*) write(*,wfmt) "'driving_wavenumber' is too small for the domain!" kf = 1.0d+00 / min(xlen, ylen, zlen) write(*,wfmt) "Resetting it to the minimum allowed value!" end if 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) ! get the wavenumber step and the maximum value due to the domain size ! kxs = max(1, int(1.0d+00 / xlen)) kys = max(1, int(1.0d+00 / ylen)) #if NDIMS == 3 kzs = max(1, int(1.0d+00 / zlen)) #endif /* NDIMS == 3 */ kxm = kmax - mod(kmax, kxs) kym = kmax - mod(kmax, kys) #if NDIMS == 3 kzm = kmax - mod(kmax, kzs) #endif /* NDIMS == 3 */ ! 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 = 0, kxm, kxs if (i > 0) then do j = - kym, kym, kys do k = - kzm, kzm, kzs 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 if end do do j = 0, kym, kys if (j > 0) then do k = - kzm, kzm, kzs 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 if end do do k = 0, kzm, kzs if (k > 0) then 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 if end do #else /* NDIMS == 3 */ do i = 0, kxm, kxs if (i > 0) then do j = - kym, kym, kys 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 if end do do j = 0, kym, kys if (j > 0) then 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 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 = 0, kxm, kxs if (i > 0) then do j = - kym, kym, kys do k = - kzm, kzm, kzs 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 if end do do j = 0, kym, kys if (j > 0) then do k = - kzm, kzm, kzs 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 if end do do k = 0, kzm, kzs if (k > 0) then 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 if end do #else /* NDIMS == 3 */ do i = 0, kxm, kxs if (i > 0) then do j = - kym, kym, kys 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 if end do do j = 0, kym, kys if (j > 0) then 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 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 if (vscale > 0.0d+00) then tscale = lscale / vscale else vscale = lscale / tscale end if 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 if (vscale > 0.0d+00) then tscale = lscale / vscale fpower = 5.0d-01 * (vscale**2 / lscale)**2 else tscale = sqrt(lscale / sqrt(2.0d+00 * fpower)) vscale = lscale / tscale end if 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 coordinates, only : periodic 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 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)) & call print_parameter(verbose, "profile Y transition between", ylo, yup) #if NDIMS == 3 if (.not. periodic(3)) & call print_parameter(verbose, "profile Z transition between", zlo, zup) #endif /* NDIMS == 3 */ end if !------------------------------------------------------------------------------- ! end subroutine print_forcing ! !=============================================================================== ! ! subroutine UPDATE_FORCING_EDDIES: ! -------------------------------- ! ! This subroutine adds the energy injection terms using simple eddy mixing. ! ! Arguments: ! ! t - time; ! dt - time step; ! injected - flag indicating whether the injection was performed; ! !=============================================================================== ! subroutine update_forcing_eddies(t, dt, injected) use coordinates, only : xmin, ymin, zmin, xlen, ylen, zlen use random , only : randuni, randsym implicit none real(kind=8), intent(in ) :: t, dt logical , intent( out) :: injected integer :: ni, n #if NDIMS == 3 real(kind=8) :: tmp #endif /* NDIMS == 3 */ real(kind=8), dimension(3) :: xp, ap !------------------------------------------------------------------------------- ! injected = .false. if (t < tinj_start .or. t > tinj_stop) then arms = 0.0d+00 dinj = 0.0d+00 rinj = 0.0d+00 return end if ! 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 injected = .true. end if !------------------------------------------------------------------------------- ! end subroutine update_forcing_eddies ! !=============================================================================== ! ! subroutine UPDATE_FORCING_OH: ! ---------------------------- ! ! This subroutine adds the energy injection terms using Ornstein–Uhlenbeck ! driving. ! ! Arguments: ! ! t - time; ! dt - time step; ! injected - flag indicating whether the injection was performed; ! !=============================================================================== ! subroutine update_forcing_oh(t, dt, injected) #if NDIMS == 3 use constants, only : pi2 #endif /* NDIMS == 3 */ use random , only : randuni, randnorz implicit none real(kind=8), intent(in ) :: t, dt logical , intent( out) :: injected 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 !------------------------------------------------------------------------------- ! injected = .false. if (t < tinj_start .or. t > tinj_stop) then arms = 0.0d+00 dinj = 0.0d+00 rinj = 0.0d+00 return end if ! 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 injected = .true. !------------------------------------------------------------------------------- ! end subroutine update_forcing_oh ! !=============================================================================== ! ! subroutine UPDATE_FORCING_ALVELIUS: ! ---------------------------------- ! ! This subroutine adds the energy injection terms using Alvelius ! driving. ! ! Arguments: ! ! t - time; ! dt - time step; ! injected - flag indicating whether the injection was performed; ! !=============================================================================== ! subroutine update_forcing_alvelius(t, dt, injected) use constants, only : pi2 use random , only : randuni implicit none real(kind=8), intent(in ) :: t, dt logical , intent( out) :: injected 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 */ !------------------------------------------------------------------------------- ! injected = .false. if (t < tinj_start .or. t > tinj_stop) then arms = 0.0d+00 dinj = 0.0d+00 rinj = 0.0d+00 return end if ! 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 injected = .true. !------------------------------------------------------------------------------- ! 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, data_blocks, get_dblocks implicit none real(kind=8), dimension(3), intent(in) :: xp, ap integer :: m, n type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! n = get_dblocks() !$omp parallel do default(shared) private(pdata) do m = 1, n pdata => data_blocks(m)%ptr call inject_eddy_block(pdata, xp(:), ap(:)) end do !$omp end parallel do !------------------------------------------------------------------------------- ! 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 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 */ #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ ! 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, data_blocks, get_dblocks implicit none real(kind=8), intent(in) :: dt integer :: m, n type(block_data), pointer :: pdata !------------------------------------------------------------------------------- ! n = get_dblocks() !$omp parallel do default(shared) private(pdata) do m = 1, n pdata => data_blocks(m)%ptr call inject_fmodes_block(pdata, dt) end do !$omp end parallel do !------------------------------------------------------------------------------- ! 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 coordinates, only : periodic 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, 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):: xk, yk, kx, ky, snkx, snky, cskx, csky #if NDIMS == 3 real(kind=8), dimension(nn):: zk, kz, snkz, cskz #endif /* NDIMS == 3 */ real(kind=8), save :: deinj real(kind=8), dimension(:,:,:,:), pointer, save :: acc real(kind=8), dimension(:,:,:) , pointer, save :: den integer, save :: nt !$ integer :: omp_get_thread_num !$omp threadprivate(first, nt, acc, den, deinj) character(len=*), parameter :: loc = 'FORCING:inject_fmodes_block()' !------------------------------------------------------------------------------- ! nt = 0 !$ nt = omp_get_thread_num() 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,nt) den(1:nn,1:nn,1:nn) => work(i+1:j,nt) #else /* NDIMS == 3 */ acc(1:2,1:nn,1:nn,1: 1) => work( 1:i,nt) den(1:nn,1:nn,1: 1) => work(i+1:j,nt) #endif /* NDIMS == 3 */ first = .false. end if #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ x(:) = pdata%meta%xmin + ax(pdata%meta%level,:) y(:) = pdata%meta%ymin + ay(pdata%meta%level,:) #if NDIMS == 3 z(:) = pdata%meta%zmin + az(pdata%meta%level,:) #else /* NDIMS == 3 */ z(:) = 0.0d+00 #endif /* NDIMS == 3 */ xk(:) = - pi2 * x(:) yk(:) = - pi2 * y(:) #if NDIMS == 3 zk(:) = - pi2 * z(:) #endif /* NDIMS == 3 */ dvol = advol(pdata%meta%level) if (work_in_use(nt)) & call print_message(loc, "Workspace is being used right now! " // & "Corruptions can occur!") work_in_use(nt) = .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) * xk(:) ky(:) = kmodes(l,2) * yk(:) #if NDIMS == 3 kz(:) = kmodes(l,3) * zk(:) #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 ! apply the injection profiles ! if (.not. periodic(1)) then do i = 1, nn if (abs(x(i)) > xup) then acc(1:NDIMS,i,:,:) = 0.0d+00 else if (abs(x(i)) >= xlo) then acc(1:NDIMS,i,:,:) = acc(1:NDIMS,i,:,:) & * cos(xfc * (abs(x(i)) - xlo))**2 end if end do end if if (.not. periodic(2)) then do j = 1, nn if (abs(y(j)) > yup) then acc(1:NDIMS,:,j,:) = 0.0d+00 else if (abs(y(j)) >= ylo) then acc(1:NDIMS,:,j,:) = acc(1:NDIMS,:,j,:) & * cos(yfc * (abs(y(j)) - ylo))**2 end if end do end if #if NDIMS == 3 if (.not. periodic(3)) then do k = 1, nn if (abs(z(k)) > zup) then acc(1:NDIMS,:,:,k) = 0.0d+00 else if (abs(z(k)) >= zlo) then acc(1:NDIMS,:,:,k) = acc(1:NDIMS,:,:,k) & * cos(zfc * (abs(z(k)) - zlo))**2 end if end do end if #endif /* NDIMS == 3 */ ! 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 deinj = sum(den(nb:ne,nb:ne,nb:ne)) * dvol #else /* NDIMS == 3 */ deinj = sum(den(nb:ne,nb:ne, 1 )) * dvol #endif /* NDIMS == 3 */ !$omp atomic update einj = einj + deinj ! add the energy increment ! if (ien > 0) then pdata%u(ien,:,:,:) = pdata%u(ien,:,:,:) + den(:,:,:) end if work_in_use(nt) = .false. 100 continue !------------------------------------------------------------------------------- ! end subroutine inject_fmodes_block ! !=============================================================================== ! ! subroutine GET_VCOEFS: ! --------------------- ! ! Subroutine determines velocity Fourier coefficients. ! ! !=============================================================================== ! subroutine get_vcoefs() use blocks , only : block_data, data_blocks, get_dblocks use helpers , only : print_message #ifdef MPI use mpitools, only : reduce_sum #endif /* MPI */ implicit none integer :: m, n, status type(block_data), pointer :: pdata complex(kind=8), dimension(:,:,:), allocatable :: vc character(len=*), parameter :: loc = 'FORCING:get_vcoefs()' !------------------------------------------------------------------------------- ! n = get_dblocks() allocate(vc(nmodes, NDIMS, n), stat=status) if (status /= 0) then call print_message(loc, & "Could not allocate memory for the Fourier coefficients!") go to 100 end if vc(:,:,:) = cmplx(0.0d+00, 0.0d+00, kind=8) !$omp parallel do default(shared) private(pdata) do m = 1, n pdata => data_blocks(m)%ptr call get_vcoefs_block(pdata, vc(:,:,m)) end do !$omp end parallel do vcoefs = sum(vc, 3) deallocate(vc, stat=status) if (status /= 0) & call print_message(loc, & "Could not release memory of the Fourier coefficients!") 100 continue #ifdef MPI 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; ! vc - an array for the velocity Fourier coefficients; ! !=============================================================================== ! subroutine get_vcoefs_block(pdata, vc) 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 complex(kind=8), dimension(:,:), intent(inout) :: vc integer :: i, j, k, 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 */ !------------------------------------------------------------------------------- ! #if NDIMS == 2 k = 1 #endif /* NDIMS == 2 */ ! 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 vc(l,1) = vc(l,1) + pdata%q(ivx,i,j,k) * cf vc(l,2) = vc(l,2) + pdata%q(ivy,i,j,k) * cf #if NDIMS == 3 vc(l,3) = vc(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