!> 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, finish_kt_grids_range, get_sizes_and_grids_range public :: check_kt_grids_range, wnml_kt_grids_range public :: kt_grids_range_config_type, set_kt_grids_range_config, get_kt_grids_range_config integer :: naky, ntheta0, n0_min, n0_max, kyspacingopt_switch real :: akx_min, akx_max, aky_min, aky_max, theta0_min, theta0_max character(20) :: kyspacing_option logical :: parameters_read = .false., 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 ky modes. integer :: naky = 1 !> Number of \(\theta_0\) (kx) modes integer :: ntheta0 = 1 !> 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_kt_grids_range implicit none parameters_read = .false. ; initialized = .false. call kt_grids_range_config%reset() end subroutine finish_kt_grids_range !> FIXME : Add documentation subroutine init_kt_grids_range(rhostar, kt_grids_range_config_in) use theta_grid, only: drhodpsi use file_utils, only: error_unit use mp, only: mp_abort, proc0 implicit none real, intent(in) :: rhostar 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 = .true. 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 if (proc0) write(ierr,'("Cannot use kyspacing_option=",A," with aky_min<=0.0 --> setting to",A)') & "'exponential'","'linear'" kyspacingopt_switch=kyspacingopt_linear end if end select if (n0_min > 0) then !Important to only do the following check and fix if naky > 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 (naky > 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, naky - 1) /= 0) then !Give a warning message that we're changing things if (proc0) then write(ierr,'("Warning: toroidal mode number range setup would lead to non-integer")') write(ierr,'(" mode numbers --> Attempting to adjust range slightly.")') end if !n0_max should be n0_min + I*(naky-1) or n0_min + (I+1)*(naky-1) !where I is int(n0_max-n0_min/(naky-1)) and we add 1 if the !remainder (n0_max-n0_min)/(naky-1) - I is > 0.5 !Note it's naky-1 rather than naky as this is number of intervals !First calculate the floating step size n0_tmp = real(n0_max - n0_min) / (naky - 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))) * (naky-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, naky-1) /= 0) & call mp_abort("Attempt to fix toroidal mode number range failed - aborting.", .true.) end if end if !If n0_min>n0_max swap values if (n0_min > n0_max) then if (proc0) write(ierr,'("Warning: Swapping max and min n0 values")') n0_swap = n0_min n0_min = n0_max n0_max = n0_swap end if !If n0_min == n0_max ensure naky=1 if (n0_min == n0_max) then if (naky /= 1 .and. proc0) write(ierr,'("Warning: Forcing naky=1 as n0_min==n0_max.")') naky = 1 end if !If there's only one mode then force n0_max=n0_min if (naky == 1) n0_max = n0_min !Set the upper and lower aky limits aky_max = n0_max * drhodpsi * rhostar aky_min = n0_min * drhodpsi * rhostar end if 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 get_sizes_and_grids_range (aky_out, theta0_out, akx_out, ikx_out, & naky_x, ntheta0_x, nx, ny) use theta_grid, only: shat use warning_helpers, only: is_zero implicit none real, dimension(:), allocatable, intent(out) :: aky_out, akx_out real, dimension(:,:), allocatable, intent(out) :: theta0_out integer, dimension(:), allocatable, intent(out) :: ikx_out integer, intent (out) :: naky_x, ntheta0_x, nx, ny real :: dkx, dky, dtheta0 integer :: i, j naky_x = naky ; ntheta0_x = ntheta0 nx = 0 ; ny = 0 allocate(aky_out(naky), akx_out(ntheta0), ikx_out(ntheta0), theta0_out(ntheta0, naky)) dky = 0.0 if (naky > 1)then select case (kyspacingopt_switch) case (kyspacingopt_linear) dky = (aky_max - aky_min)/real(naky - 1) aky_out = [ (aky_min + dky*real(i), i = 0,naky-1) ] case (kyspacingopt_exp) dky = (log(aky_max) - log(aky_min))/real(naky - 1) aky_out = [ (exp(log(aky_min) + dky*real(i)), i = 0,naky-1) ] end select else aky_out = [ (aky_min, i = 0,naky-1) ] endif ! set default theta0 to 0 theta0_out = 0.0 ! ! BD: Assumption here differs from convention that abs(shat) <= 1.e-5 triggers periodic bc ! if (is_zero(shat)) then !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_out = [ (akx_min + dkx*real(i), i = 0,ntheta0-1) ] else ! 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_out(:,j) & = [ (theta0_min + dtheta0*real(i), i=0,ntheta0-1) ] end do ! Adding support for ky=0, kx/=0 if(is_zero(aky_out(1)))then if(naky>1)then akx_out = theta0_out(:,2) * shat * aky_out(2) else dkx = 0.0 if (ntheta0 > 1) dkx = (akx_max - akx_min)/real(ntheta0 - 1) akx_out = [ (akx_min + dkx*real(i), i = 0,ntheta0-1) ] end if else !This is the original behaviour akx_out = theta0_out(:,1) * shat * aky_out(1) endif end if do j = 1, ntheta0 ikx_out(j) = j - 1 end do end subroutine get_sizes_and_grids_range !> FIXME : Add documentation subroutine check_kt_grids_range(rhostar, report_unit) use constants, only: twopi use theta_grid, only: shat use warning_helpers, only: is_not_zero implicit none real, intent(in) :: rhostar integer, intent(in) :: report_unit real :: dtheta0 integer :: i, j integer :: naky, ntheta0, nx, ny ! Shadow module level vars 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=',1pe12.4)") n0_min,rhostar 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, *) call get_sizes_and_grids_range(aky, theta0, akx, ikx, naky, ntheta0, nx, ny) !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 (is_not_zero(shat)) 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