kt_grids_range.f90 Source File


Contents

Source Code


Source Code

!> Set up ranges of kx and ky for linear runs.
module kt_grids_range
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN

  implicit none

  private

  public :: init_kt_grids_range, range_get_sizes, range_get_grids
  public :: check_kt_grids_range
  public :: wnml_kt_grids_range
  public :: read_parameters_range
  public :: finish_parameters_range

  public :: kt_grids_range_config_type
  public :: set_kt_grids_range_config
  public :: get_kt_grids_range_config

  integer :: naky, ntheta0, nn0, n0_min, n0_max
  real :: aky_min, aky_max, theta0_min, theta0_max
  real :: akx_min, akx_max, rhostar_range
  character(20) :: kyspacing_option
  integer :: kyspacingopt_switch

  logical :: parameters_read = .false.
  logical :: initialized = .false.

  !Note if we ever want to offer different spacing for theta0 we could
  !reuse these flags (rename to spacingopt_...).
  integer, parameter :: kyspacingopt_linear=1, kyspacingopt_exp=2

  !> Used to represent the input configuration of kt_grids_range
  type, extends(abstract_config_type) :: kt_grids_range_config_type
     ! namelist : kt_grids_range_parameters
     ! indexed : false
     !> Max kx for periodic finite kx ballooning space runs with
     !> \(\hat{s}=0\).
     real :: akx_max = 0.0
     !> Min kx for periodic finite kx ballooning space runs with
     !> \(\hat{s}=0\).
     real :: akx_min = 0.0
     !> Upper limit of \(k_y \rho\) range. Should set to something other
     !> than zero.
     real :: aky_max = 0.0
     !> Lower limit of \(k_y \rho\) range. Should typically set to
     !> something other than zero.
     real :: aky_min = 0.0
     !> Sets the type of spacing between ky grid points, available options are :
     !>
     !> -  'default' : Same as 'linear'
     !> -  'exponential' : Evenly spaced in log(ky).
     !> -  'linear' : Evenly spaced in ky.
     !>
     character(len = 20) :: kyspacing_option = 'default'
     !> Maximum toroidal mode number. Can use instead of `aky_max`.
     integer :: n0_max = 0
     !> Minimum toroidal mode number. Can use instead of `aky_min`.
     integer :: n0_min = 0
     !> The number of 'actual' ky modes.
     integer :: naky = 1
     !> Number of toroidal modes, only used if `n0_min`>0. Overrides `naky`
     !> in kt_grids_range_parameters. Note if the number of modes isn't compatible
     !> with the requested min and max toroidal mode numbers then we just run with
     !> one mode, determined by `n0_max`.
     integer :: nn0 = 1
     !> Number of \(\theta_0\) (kx) modes
     integer :: ntheta0 = 1
     !> Used to convert `n0_min`, `n0_max` range into `aky_min`,
     !> `aky_max`, if `n0_min`>0. If `n0_min` is set,
     !> `aky_min=n0_min*drhodpsi*rhostar_range` and
     !> `aky_max=n0_max*drhodpsi*rhostar_range` where `drhodpsi` is
     !> calculated as a part of the geometry setup.
     real :: rhostar_range = 1.0e-4
     !> Upper limit of `theta_0` range
     real :: theta0_max = 0.0
     !> Lower limit of `theta_0` range
     real :: theta0_min = 0.0
   contains
     procedure, public :: read => read_kt_grids_range_config
     procedure, public :: write => write_kt_grids_range_config
     procedure, public :: reset => reset_kt_grids_range_config
     procedure, public :: broadcast => broadcast_kt_grids_range_config
     procedure, public, nopass :: get_default_name => get_default_name_kt_grids_range_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_kt_grids_range_config
  end type kt_grids_range_config_type

  type(kt_grids_range_config_type) :: kt_grids_range_config

contains

  !> FIXME : Add documentation
  subroutine read_parameters_range(kt_grids_range_config_in)
    use file_utils, only: error_unit
    use text_options, only: text_option, get_option_value
    implicit none
    type(kt_grids_range_config_type), intent(in), optional :: kt_grids_range_config_in
    integer :: ierr
    type (text_option), dimension(3), parameter :: kyspacingopts = &
         [ text_option('default', kyspacingopt_linear), &
         text_option('linear', kyspacingopt_linear), &
         text_option('exponential', kyspacingopt_exp) ]

    if (parameters_read) return
    parameters_read = .true.

    if (present(kt_grids_range_config_in)) kt_grids_range_config = kt_grids_range_config_in

    call kt_grids_range_config%init(name = 'kt_grids_range_parameters', requires_index = .false.)

    ! Copy out internal values into module level parameters
    associate(self => kt_grids_range_config)
#include "kt_grids_range_copy_out_auto_gen.inc"
    end associate

    ierr = error_unit()
    call get_option_value(kyspacing_option, kyspacingopts, kyspacingopt_switch,&
         ierr, "kyspacing_option in kt_grids_range_parameters",.true.)
  end subroutine read_parameters_range

  !> FIXME : Add documentation
  subroutine finish_parameters_range
    implicit none
    parameters_read = .false.
    initialized = .false.
    call kt_grids_range_config%reset()
  end subroutine finish_parameters_range

  !> FIXME : Add documentation
  subroutine init_kt_grids_range(kt_grids_range_config_in)
    !CMR, 14/10/2013:
    ! New namelist variables nn0, n0_min, n0_max, rhostar_range to set ky grid
    !                                             using toroidal mode numbers.
    ! Toroidal modenumbers are used if n0_min> 0 prescribed in input file.
    use theta_grid, only: drhodpsi
    use file_utils, only: error_unit
    implicit none
    type(kt_grids_range_config_type), intent(in), optional :: kt_grids_range_config_in
    integer :: ierr
    !! Temporary value of n0 as floating point
    real :: n0_tmp
    !! Temporary variable for swapping n0_{min,max}
    integer :: n0_swap

    if (initialized) return
    initialized = .false.

    call read_parameters_range(kt_grids_range_config_in)

    ierr = error_unit()

    !Override kyspacing_option in certain cases
    select case (kyspacingopt_switch)
    case (kyspacingopt_exp)
       if(aky_min<=0) then
          write(ierr,'("Cannot use kyspacing_option=",A," with aky_min<=0.0 --> setting to",A)') &
               "'exponential'","'linear'"
          kyspacingopt_switch=kyspacingopt_linear
       endif
    end select

    if (n0_min > 0) then
       !CMR if n0_min>0 then override aky inputs and use nn0, n0_min, n0_max to determine aky range

       !Important to only do the following check and fix if nn0 > 0 as the second argument of `mod`
       !must be non-zero. Failing to guard against this can lead to different outcomes with different
       !compilers
       if(nn0 > 1) then

          !If toroidal mode number range would lead to non-integer mode numbers then
          !we adjust the range to try to fix this.
          if(mod(n0_max-n0_min,nn0-1)/=0) then
             !Give a warning message that we're changing things
             write(ierr,'("Warning: toroidal mode number range setup would lead to non-integer")')
             write(ierr,'("         mode numbers --> Attempting to adjust range slightly.")')

             !n0_max should be n0_min + I*(nn0-1) or n0_min + (I+1)*(nn0-1)
             !where I is int(n0_max-n0_min/(nn0-1)) and we add 1 if the
             !remainder (n0_max-n0_min)/(nn0-1) - I is > 0.5
             !Note it's nn0-1 rather than nn0 as this is number of intervals

             !First calculate the floating step size
             n0_tmp = (n0_max-n0_min*1.0)/(nn0-1)

             !Now construct the new upper limit, int(n0_tmp) = I
             !nint(n0_tmp-int(n0_tmp)) should be either 0 or 1 depending
             !on if the remainder is < or > 0.5
             n0_max = n0_min + (int(n0_tmp)+nint(n0_tmp-int(n0_tmp)))*(nn0-1)

             !Double check it's all OK now, should always be fine but just
             !putting this here in case of logic error or strange cases.
             if(mod(n0_max-n0_min,nn0-1)/=0) then
                write(ierr,'("Error: Attempt to fix toroidal mode number range failed. Forcing nn0=1")')
                nn0=1
                n0_max = n0_min
             endif
          endif
       end if

       !If n0_min>n0_max swap values
       if(n0_min>n0_max) then
          write(ierr,'("Warning: Swapping max and min n0 values")')
          n0_swap = n0_min
          n0_min = n0_max
          n0_max = n0_swap
       endif

       !If n0_min == n0_max ensure nn0=1
       if(n0_min==n0_max) then
          if(nn0/=1) write(ierr,'("Warning: Forcing nn0=1 as n0_min==n0_max.")')
          nn0 = 1
       endif

       !If there's only one mode then force n0_max=n0_min
       if(nn0 == 1) then
          n0_max = n0_min
       endif

       !Set the upper and lower aky limits
       aky_max=n0_max*drhodpsi*rhostar_range
       aky_min=n0_min*drhodpsi*rhostar_range

       !Set the number of aky values
       naky=nn0
    endif

  end subroutine init_kt_grids_range

  !> FIXME : Add documentation
  subroutine wnml_kt_grids_range(unit)
    implicit none
    integer, intent(in) :: unit
    call kt_grids_range_config%write(unit)
  end subroutine wnml_kt_grids_range

  !> FIXME : Add documentation
  subroutine range_get_sizes (naky_x, ntheta0_x, nx, ny)
    implicit none
    integer, intent (out) :: naky_x, ntheta0_x, nx, ny
    naky_x = naky  ;  ntheta0_x = ntheta0
    nx = 0         ;  ny = 0
  end subroutine range_get_sizes

  !> FIXME : Add documentation
  subroutine range_get_grids (aky, theta0, akx, ikx)
    ! BD: Could add some logic here to set theta0 if akx is given?  When do we need what?
    use theta_grid, only: shat
    use mp, only: mp_abort
    implicit none
    real, dimension (:), intent (out) :: akx, aky
    real, dimension (:,:), intent (out) :: theta0
    !> Discrete kx wavenumber grid indices
    integer, dimension (:), intent (out) :: ikx

    real :: dkx, dky, dtheta0
    integer :: i, j

    if ( size(aky) /= naky) then
       call mp_abort('range_get_grids: size(aky) /= naky',.true.)
    endif

    if ( size(akx) /= ntheta0) then
       call mp_abort('range_get_grids: size(akx) /= ntheta0',.true.)
    endif

    dky = 0.0
    if (naky > 1)then
       select case (kyspacingopt_switch)
       case (kyspacingopt_linear)
          dky = (aky_max - aky_min)/real(naky - 1)
          aky = [ (aky_min + dky*real(i), i = 0,naky-1) ]
       case (kyspacingopt_exp)
          dky = (log(aky_max) - log(aky_min))/real(naky - 1)
          aky = [ (exp(log(aky_min) + dky*real(i)), i = 0,naky-1) ]
       end select
    else
       aky = [ (aky_min, i = 0,naky-1) ]
    endif

    ! set default theta0 to 0
    theta0=0.0

    !
    ! BD: Assumption here differs from convention that abs(shat) <= 1.e-5 triggers periodic bc
    !
    if (shat /= 0.0d0) then  ! ie assumes boundary_option .eq. 'linked'
       dtheta0 = 0.0
       if (ntheta0 > 1) dtheta0 = (theta0_max - theta0_min)/real(ntheta0 - 1)

       do j = 1, naky
          theta0(:,j) &
               = [ (theta0_min + dtheta0*real(i), i=0,ntheta0-1) ]
       end do

       !<DD>Adding support for ky=0, kx/=0
       if(aky(1)==0)then
          if(naky>1)then
             akx = theta0(:,2) * shat * aky(2)
          else
             dkx = 0.0
             if (ntheta0 > 1) dkx = (akx_max - akx_min)/real(ntheta0 - 1)
             akx = [ (akx_min + dkx*real(i), i = 0,ntheta0-1) ]
          end if
       else
          !This is the original behaviour
          akx = theta0(:,1) * shat * aky(1)
       endif
    else

       !CMR, 22/9/2010:  ie here assume boundary_option .eq. 'periodic'
       !new code for periodic finite kx ballooning space runs with shat=0
       dkx = 0.0
       if (ntheta0 > 1) dkx = (akx_max - akx_min)/real(ntheta0 - 1)
       akx = [ (akx_min + dkx*real(i), i = 0,ntheta0-1) ]
    endif

    do j = 1, ntheta0
       ikx(j) = j - 1
    end do
  end subroutine range_get_grids

  !> FIXME : Add documentation
  subroutine check_kt_grids_range(report_unit)
    use constants, only: twopi
    use theta_grid, only: shat
    implicit none
    integer, intent(in) :: report_unit
    real :: dtheta0
    integer :: i, j
    real, dimension(:), allocatable:: aky, akx
    real, dimension(:,:), allocatable:: theta0
    integer, dimension(:), allocatable :: ikx

    write (report_unit, *)
    write (report_unit, fmt="('A range of k_perps will be evolved.')")
    if (n0_min >0) write (report_unit, fmt="('ky set using toroidal mode numbers with n0_min=',i8/T34,'rhostar_range=',1pe12.4)") n0_min,rhostar_range
    write (report_unit, *)
    write (report_unit, fmt="('There are ',i3,' values of ky rho and ',i3,' values of theta_0/kx rho:')") naky, ntheta0
    write (report_unit, *)

    !<DD>Calculate the kt grids
    allocate(aky(naky),theta0(ntheta0,naky),akx(ntheta0), ikx(ntheta0))
    call range_get_grids(aky, theta0, akx, ikx)

    !Report grid values
    do j = 1, naky
       do i = 1, ntheta0
          write (report_unit, fmt="('ky rho = ',e11.4,' theta0 = ',e11.4,' kx rho = ',e11.4)") &
               aky(j),theta0(i,j),akx(i)
       end do
    end do
    deallocate(aky,theta0,akx)

    ! CMR, add some !!!error checking!!! for ballooning space runs for shat /= 0
    ! using flow shear: check that the constraints on theta0 grid are satisfied!

    if (shat /= 0) then
       !It would be nice to only write this information if g_exb*gexbfac/=0 but currently
       !dependencies prevent this.
       dtheta0 = 0.0    ;  if (ntheta0 > 1) dtheta0 = (theta0_max - theta0_min)/real(ntheta0 - 1)
       if (abs(mod(twopi-theta0_max+theta0_min,twopi)-dtheta0) > 1.0e-3*dtheta0) then
          write (report_unit, *)
          write (report_unit, fmt="('IF using perp ExB flow shear in BALLOONING SPACE there is an ERROR that will corrupt results.')")
          write (report_unit, fmt="('check_kt_grids_range: inappropriate theta0 grid')")
          write (report_unit, fmt="('In ballooning space with sheared flow, 2pi-theta0_max+theta0_min =',e11.4,' must be set equal to dtheta = ',e11.4)") twopi-theta0_max+theta0_min, dtheta0
       endif
    endif

  end subroutine check_kt_grids_range

#include "kt_grids_range_auto_gen.inc"
end module kt_grids_range