ballstab.f90 Source File


Source Code

Source Code

!> A small module used to calculate ideal ballooning stability.
!! Based on the original ball program in geo/ball.f90 of GS2
!! but tweaked to integrate closer to full GS2 runs etc.
module ballstab
  use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN
  implicit none


  public :: init_ballstab, finish_ballstab
  public :: run_stability_check, is_unstable
  public :: write_stability_ascii_1d, write_stability_ascii_2d

  public :: stability

  public :: ballstab_config_type
  public :: set_ballstab_config
  public :: get_ballstab_config
  logical :: make_salpha
  integer :: n_shat, n_beta
  real :: shat_min, shat_max, theta0
  real :: beta_mul, beta_div
  real, dimension(:), allocatable :: shat_arr, beta_arr, dbdrho_arr
  integer, dimension(:,:), allocatable :: stability
  real :: diff

  logical :: initialised = .false. 
  character(len=14) :: namelist_name="ballstab_knobs"

  !> Used to represent the input configuration of ballstab
  type, extends(abstract_config_type) :: ballstab_config_type
     ! namelist : ballstab_knobs
     ! indexed : false
     !> The minimum \(\beta^\prime\) in scans is
     !> `beta_prime_equilibrium/beta_div`
     real :: beta_div = 1.0
     !> The maximum \(\beta^\prime\) in scans is
     !> `beta_prime_equilibrium*beta_mul`
     real :: beta_mul = 1.0
     !> Controls decentring used in the numerical integration used for
     !> solving the ballooning equation. Recommended values are either
     !> 0 (default) or 1/3.
     real :: diff = 0.0
     !> Not currently used for anything.
     !> @todo Consider removing this variable
     logical :: make_salpha = .false.
     !> How many \(\beta^\prime\) values should be used in s-alpha
     !> type scans.
     integer :: n_beta = 1
     !> How many \(\hat{s}\) values should be used in s-alpha type
     !> scans.
     integer :: n_shat = 1
     !> The maximum value of \(\hat{s}\) to use in s-alpha type scans.
     !> @note This gets a smart default of the equilibrium value of \(\hat{s}\)
     real :: shat_max = 0.0
     !> The minimum value of \(\hat{s}\) to use in s-alpha type scans.
     !> @note This gets a smart default of the equilibrium value of \(\hat{s}\)
     real :: shat_min = 0.0
     !> The theta0 value to use in solving the ideal ballooning problem. We don't
     !> currently provide a means to study multiple theta0 in a single run. This
     !> may change in the future.
     real :: theta0 = 0.0
     procedure, public :: read => read_ballstab_config
     procedure, public :: write => write_ballstab_config
     procedure, public :: reset => reset_ballstab_config
     procedure, public :: broadcast => broadcast_ballstab_config
     procedure, public, nopass :: get_default_name => get_default_name_ballstab_config
     procedure, public, nopass :: get_default_requires_index => get_default_requires_index_ballstab_config
  end type ballstab_config_type

  type(ballstab_config_type) :: ballstab_config
  !> Initialise grids, other modules etc.
  subroutine init_ballstab(ballstab_config_in)
    use theta_grid, only: init_theta_grid
    implicit none
    type(ballstab_config_type), intent(in), optional :: ballstab_config_in
    !Other modules
    call init_theta_grid

    !This module
    call real_init_ballstab(ballstab_config_in)
  end subroutine init_ballstab

  !> Finalise this module
  subroutine finish_ballstab
    implicit none
    call dealloc_arrays

    initialised = .false.
    call ballstab_config%reset()
  end subroutine finish_ballstab

  !> Initialise this module
  subroutine real_init_ballstab(ballstab_config_in)
    use mp, only: proc0
    implicit none
    type(ballstab_config_type), intent(in), optional :: ballstab_config_in    

    if(initialised) return
    initialised = .true.

    !Get settings
    call read_parameters(ballstab_config_in)

    !Note only proc0 will be allowed to do any work
    if (.not.proc0) return
    !Verify parameters
    call verify_parameters

    !Create arrays
    call setup_arrays

  end subroutine real_init_ballstab

  !> Read parameters
  subroutine read_parameters(ballstab_config_in)
    use theta_grid, only: shat
    implicit none
    type(ballstab_config_type), intent(in), optional :: ballstab_config_in    
    logical :: exist

    if (present(ballstab_config_in)) ballstab_config = ballstab_config_in

    !Smart defaults
    if (.not.ballstab_config%is_initialised()) then    
       ballstab_config%shat_min = shat
       ballstab_config%shat_max = shat
    end if
    call ballstab_config%init(name = 'ballstab_knobs', requires_index = .false.)

    ! Copy out internal values into module level parameters
    beta_div = ballstab_config%beta_div
    beta_mul = ballstab_config%beta_mul
    diff = ballstab_config%diff
    make_salpha = ballstab_config%make_salpha
    n_beta = ballstab_config%n_beta
    n_shat = ballstab_config%n_shat
    shat_max = ballstab_config%shat_max
    shat_min = ballstab_config%shat_min
    theta0 = ballstab_config%theta0

    exist = ballstab_config%exist
  end subroutine read_parameters

  !> Check parameters are consistent etc.
  subroutine verify_parameters
    use geometry, only: bishop, local_eq
    use file_utils, only: error_unit
    implicit none
    integer :: ierr

    ierr = error_unit()

    !Disable scan in s_hat when incompatible bishop used
    if (n_shat > 1) then
       select case(bishop)
          !These are valid values
       case default
          if( (bishop /= 0) .and. local_eq) then
             !This is also ok
             !This is not ok
             n_shat = 1
             write(ierr,'("ERROR : Cannot scan in shat for current combination of bishop and equilbrium type")')
          end if
       end select
    end if

    !Disable scan in beta_prime when incompatible bishop used
    if (n_beta > 1) then
       select case(bishop)
          !These are valid values
       case default
          !This is not ok
          n_beta = 1
          write(ierr,'("ERROR : Cannot scan in beta_prime for current combination of bishop and equilbrium type")')
       end select
    end if
  end subroutine verify_parameters

  !> Allocate+populate arrays
  subroutine setup_arrays
    implicit none
    integer :: is, ib
    real :: ds, db, b_min, b_max, bp
    dbdrho_arr = 0. !This will be populated later

    stability = 0

    ds = 0
    if( ds = (shat_max-shat_min)/(n_shat-1)
    do is=1,n_shat

    db = 0
    call get_beta_prime(bp)
    b_max = bp*beta_mul
    b_min = bp/beta_div
    if( then
       db = (b_max-b_min)/(n_beta-1)
    end if
    do ib=1,n_beta
  end subroutine setup_arrays

  !> Deallocate arrays
  subroutine dealloc_arrays
    implicit none
    if(allocated(shat_arr)) deallocate(shat_arr)
    if(allocated(beta_arr)) deallocate(beta_arr)
    if(allocated(dbdrho_arr)) deallocate(dbdrho_arr)
    if(allocated(stability)) deallocate(stability)
  end subroutine dealloc_arrays

  !> Write a namelist based on current values
  subroutine wnml_ballstab(unit)
    implicit none
    integer, intent(in) :: unit

    write(unit,'("&",A)') trim(adjustl(namelist_name))

    write(unit,'(" make_salpha = ",L1)') make_salpha
    write(unit,'(" n_shat      = ",I0)') n_shat
    if( then
       write(unit,'(" shat_min    = ",F12.5)') shat_min
       write(unit,'(" shat_max    = ",F12.5)') shat_max
    write(unit,'(" n_beta      = ",I0)') n_beta

  end subroutine wnml_ballstab

  !> Check value consistency etc.
  subroutine check_ballstab(report_unit)
    implicit none
    integer, intent(in) :: report_unit
  end subroutine check_ballstab

  !> Gets the current value of shat
  subroutine get_shat(shat_out)
    use geometry, only: surf
    implicit none
    real, intent(out) :: shat_out
    shat_out = surf%shat
  end subroutine get_shat

  !> Sets the value of shat
  subroutine set_shat(shat_in)
    use geometry, only: surf
    implicit none
    real, intent(in) :: shat_in
    !Note for bishop = 6,7 or not in 2-9 and .not. local_eq cannot change shat
    !Should flag this somewhere (during init)
    surf%shat = shat_in
  end subroutine set_shat

  !> Gets the current value of beta_prime (or equivalent var)
  !! Really just gets the variable that we can use to control 
  !! the pressure gradient.
  subroutine get_beta_prime(bp_out)
    use geometry, only: bishop, p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult
    implicit none
    real, intent(out) :: bp_out

    select case (bishop)
       bp_out = p_prime_input
       bp_out = invLp_input
       bp_out = beta_prime_input
       bp_out = alpha_input
       bp_out = beta_prime_input
       bp_out = dp_mult
       bp_out = beta_prime_input
       bp_out = beta_prime_input
    case default
       bp_out = 0.0
       print*,"ERROR: For bishop not in 2-9 cannot modify the pressure gradient used"
    end select
  end subroutine get_beta_prime

  !> Sets the current value of beta_prime (or equivalent var)
  !! Really just sets the variable that we can use to control 
  !! the pressure gradient.
  subroutine set_beta_prime(bp_in)
    use geometry, only: bishop, p_prime_input, invLp_input, beta_prime_input, alpha_input, dp_mult
    implicit none
    real, intent(in) :: bp_in

    select case (bishop)
       p_prime_input = bp_in
       invLp_input = bp_in
       beta_prime_input = bp_in
       alpha_input = bp_in
       beta_prime_input = bp_in
       dp_mult = bp_in
       beta_prime_input = bp_in
       beta_prime_input = bp_in
    case default
       print*,"ERROR: For bishop not in 2-9 cannot modify the pressure gradient used"
    end select
  end subroutine set_beta_prime

  !> Runs the stability check scan for the module level shat/beta array values
  subroutine run_stability_check
    use file_utils, only: open_output_file, close_output_file
    implicit none
    integer :: is, ib

    !Exit if not initialised
    if(.not.initialised) return

    !Check stability over range
    do ib = 1, n_beta
       do is = 1, n_shat
          call check_stability(shat_arr(is), beta_arr(ib), &
               dbdrho_arr(ib), stability(is, ib), theta0)
       end do
    end do
  end subroutine run_stability_check

  !> Routine to write out stability data to 1D ascii file
  subroutine write_stability_ascii_1d
    use file_utils, only: open_output_file, close_output_file
    implicit none
    integer :: is, ib, unit

    call open_output_file(unit,".ballstab_1d")
    write(unit,'(4(A12," "))') "shear", "beta_prime", "dbetadrho", "stability"
    do ib=1,n_beta
       do is=1,n_shat
          write(unit,'(3(F12.4," "),I12)') shat_arr(is), beta_arr(ib), dbdrho_arr(ib), stability(is,ib)
    call close_output_file(unit)
  end subroutine write_stability_ascii_1d

  !> Routine to write out stability data to 2D ascii file + 1d axis data
  subroutine write_stability_ascii_2d
    use file_utils, only: open_output_file, close_output_file
    implicit none
    integer :: is, ib, unit

    !/2D -- data
    call open_output_file(unit,".ballstab_2d")
    do ib=1,n_beta
       do is=1,n_shat
          write(unit,'(I0," ")',advance="no") stability(is,ib)
       write(unit,'()') !Move to the next line

    !/2D -- axis
    call open_output_file(unit,".ballstab_shat")
    !Commented lines here can be used to give row output
    do is=1,n_shat
       !write(unit,'(F9.2," ")',advance="no") shat_arr(is)
       write(unit,'(F9.2)') shat_arr(is)
    call close_output_file(unit)

    call open_output_file(unit,".ballstab_bp")
    do ib=1,n_beta
       write(unit,'(F12.4)') beta_arr(ib)
    call close_output_file(unit)

    call open_output_file(unit,".ballstab_dbdr")
    do ib=1,n_beta
       write(unit,'(F12.4)') dbdrho_arr(ib)
    call close_output_file(unit)
  end subroutine write_stability_ascii_2d

  !> Test if given shat/beta is unstable
  subroutine check_stability(shat_in, beta_prime_in, dbetadrho_out, &
       iunstable, theta0, restore)
    !Note iunstable>0 means unstable - Could use a logical but size may be interesting
    use theta_grid_params, only: ntheta, nperiod
    use geometry, only: eikcoefs, eikcoefs_output_type
    use optionals, only: get_option_with_default
    implicit none
    real, intent(in) :: shat_in, beta_prime_in
    real, intent(out) :: dbetadrho_out
    integer, intent(out) :: iunstable
    real, intent(in), optional :: theta0
    logical, intent(in), optional :: restore
    type(eikcoefs_output_type) :: eikcoefs_results
    real :: shat_bak, beta_prime_bak
    integer :: ntgrid
    logical :: local_restore
    local_restore = get_option_with_default(restore, .false.)

    !Store initial state if we're going to restore
    if (local_restore) then
       call get_shat(shat_bak)
       call get_beta_prime(beta_prime_bak)
    end if

    !Setup geometry
    !NOTE: As how we can influence beta_prime (and shat) varies with bishop
    !we actually end up setting different parameters in different cases. The
    !physical parameter dbetadrho should really be the relevant parameter.
    call set_beta_prime(beta_prime_in)
    call set_shat(shat_in)
    call eikcoefs(ntheta, nperiod, eikcoefs_results) !Recalculate geometry terms

    !Here we pass out the value of dbetadrho, this should be bishop value independent
    dbetadrho_out = eikcoefs_results%dbetadrho

    !Note we don't use theta_grid:ntgrid as in the case that
    !gridgen resizes the theta grid used by GS2 we don't end up
    !resizing any of the geometry module arrays (that we use directly here)
    !Really we would like to use everything from theta_grid directly but we'd
    !need some way to force theta_grid to update it's internal arrays when we
    !change parameters

    !Test stability
    iunstable = is_unstable(ntgrid, eikcoefs_results%theta, eikcoefs_results%bmag, &
         eikcoefs_results%dbetadrho, eikcoefs_results%gradpar, eikcoefs_results%gds2, &
         eikcoefs_results%gds21, eikcoefs_results%gds22, eikcoefs_results%cvdrift, &
         eikcoefs_results%cvdrift0, theta0)

    !Restore geometry if requested
    if (local_restore) then
       call set_beta_prime(beta_prime_bak)
       call set_shat(shat_bak)
       call eikcoefs(ntheta, nperiod, eikcoefs_results)
    end if
  end subroutine check_stability

  !> Test stability of current system for passed theta0 and geometry -- return integer
  !> Note is_unstable>0 means unstable - Could use a logical but size
  !> may be interesting.
  pure integer function is_unstable(ntgrid, theta, bmag, dbetadrho, gradpar, &
       gds2, gds21, gds22, cvdrift, cvdrift0, theta0_in) &
    use optionals, only: get_option_with_default
    implicit none
    integer, intent(in) :: ntgrid
    real, dimension(-ntgrid:ntgrid), intent(in) :: theta, bmag, gradpar
    real, intent(in) :: dbetadrho
    real, dimension(-ntgrid:ntgrid), intent(in) :: gds2, gds21, gds22, cvdrift, cvdrift0
    real, intent(in), optional :: theta0_in
    real, dimension(-ntgrid:ntgrid) :: g, c, psi, c1, c2, ch, g1, g2, gh
    real, dimension(-ntgrid:ntgrid-1) :: delthet
    integer :: ig
    real :: one_m_diff, psi_prime, theta0

    theta0 = get_option_with_default(theta0_in, 0.0)
    is_unstable_out = 0
    one_m_diff = 1 - diff

    !Calculate geometry arrays
    delthet = theta(-ntgrid+1:) - theta(:ntgrid-1)

    g = (gds2 + 2 * theta0 * gds21 + theta0**2 * gds22) * gradpar / bmag
    c = -0.5 * dbetadrho * (cvdrift + theta0 * cvdrift0) / (bmag * gradpar)

    ! Get grid centre values -- note we don't set the value at -ntgrid
    do ig=-ntgrid+1, ntgrid

    !Calculate coefficients -- note we don't set the value at -ntgrid or ntgrid
    do ig=-ntgrid+1,ntgrid-1
       c1(ig) = -delthet(ig)*(one_m_diff*c(ig)+0.5*diff*ch(ig+1)) &
            -delthet(ig-1)*(one_m_diff*c(ig)+0.5*diff*ch(ig)) &

    !Calculate coefficients -- note we don't set the value at -ntgrid
    do ig=-ntgrid+1,ntgrid
       c2(ig) = -0.25*diff*ch(ig)*delthet(ig-1)
       g1(ig) = gh(ig)/delthet(ig-1)
       g2(ig) = 1.0/(0.25*diff*ch(ig)*delthet(ig-1)+gh(ig)/delthet(ig-1))

    !Set boundary values
    psi_prime=psi(-ntgrid+1)/g2(-ntgrid+1)!-g1(-ntgrid+1)*psi(-ntgrid) !Second term zero by defn

    !Integrate along in theta
    do ig=-ntgrid+1,ntgrid-1

    !Find how many times psi crosses axis || Actually how many times it hits zero
    !If we really do just want to count crossings then need to replace le with lt
    ! Could we replace this with
    ! is_unstable_out = count(psi(-ntgrid+1:ntgrid-1) * psi(-ntgrid+2:) <= 0)
    do ig=-ntgrid+1,ntgrid-1
       if(psi(ig)*psi(ig+1) .le. 0 ) is_unstable_out = is_unstable_out + 1

    !Probably want to think about writing out psi if debug
  end function is_unstable

  !> Set the module level config type
  !> Will abort if the module has already been initialised to avoid
  !> inconsistencies.
  subroutine set_ballstab_config(ballstab_config_in)
    use mp, only: mp_abort
    type(ballstab_config_type), intent(in), optional :: ballstab_config_in
    if (initialised) then
       call mp_abort("Trying to set ballstab_config when already initialised.", to_screen = .true.)
    end if
    if (present(ballstab_config_in)) then
       ballstab_config = ballstab_config_in
    end if
  end subroutine set_ballstab_config

  !> Get the module level config instance
  function get_ballstab_config()
    type(ballstab_config_type) :: get_ballstab_config
    get_ballstab_config = ballstab_config
  end function get_ballstab_config

  ! Following is for the config_type
  !> Reads in the ballstab_knobs namelist and populates the member variables
  subroutine read_ballstab_config(self)
    use file_utils, only: input_unit_exist, get_indexed_namelist_unit
    use mp, only: proc0
    implicit none
    class(ballstab_config_type), intent(in out) :: self
    logical :: exist
    integer :: in_file

    ! Note: When this routine is in the module where these variables live
    ! we shadow the module level variables here. This is intentional to provide
    ! isolation and ensure we can move this routine to another module easily.    
    integer :: n_beta, n_shat
    logical :: make_salpha
    real :: beta_div, beta_mul, diff, shat_max, shat_min, theta0

    namelist /ballstab_knobs/ beta_div, beta_mul, diff, make_salpha, n_beta, n_shat, shat_max, shat_min, theta0

    ! Only proc0 reads from file
    if (.not. proc0) return

    ! First set local variables from current values
    beta_div = self%beta_div
    beta_mul = self%beta_mul
    diff = self%diff
    make_salpha = self%make_salpha
    n_beta = self%n_beta
    n_shat = self%n_shat
    shat_max = self%shat_max
    shat_min = self%shat_min
    theta0 = self%theta0

    ! Now read in the main namelist
    in_file = input_unit_exist(self%get_name(), exist)
    if (exist) read(in_file, nml = ballstab_knobs)

    ! Now copy from local variables into type members
    self%beta_div = beta_div
    self%beta_mul = beta_mul
    self%diff = diff
    self%make_salpha = make_salpha
    self%n_beta = n_beta
    self%n_shat = n_shat
    self%shat_max = shat_max
    self%shat_min = shat_min
    self%theta0 = theta0

    self%exist = exist
  end subroutine read_ballstab_config

  !> Writes out a namelist representing the current state of the config object
  subroutine write_ballstab_config(self, unit)
    implicit none
    class(ballstab_config_type), intent(in) :: self
    integer, intent(in) , optional:: unit
    integer :: unit_internal

    unit_internal = 6 ! @todo -- get stdout from file_utils
    if (present(unit)) then
       unit_internal = unit

    call self%write_namelist_header(unit_internal)
    call self%write_key_val("beta_div", self%beta_div, unit_internal)
    call self%write_key_val("beta_mul", self%beta_mul, unit_internal)
    call self%write_key_val("diff", self%diff, unit_internal)
    call self%write_key_val("make_salpha", self%make_salpha, unit_internal)
    call self%write_key_val("n_beta", self%n_beta, unit_internal)
    call self%write_key_val("n_shat", self%n_shat, unit_internal)
    call self%write_key_val("shat_max", self%shat_max, unit_internal)
    call self%write_key_val("shat_min", self%shat_min, unit_internal)
    call self%write_key_val("theta0", self%theta0, unit_internal)
    call self%write_namelist_footer(unit_internal)
  end subroutine write_ballstab_config

  !> Resets the config object to the initial empty state
  subroutine reset_ballstab_config(self)
    class(ballstab_config_type), intent(in out) :: self
    type(ballstab_config_type) :: empty
    select type (self)
    type is (ballstab_config_type)
       self = empty
    end select
  end subroutine reset_ballstab_config

  !> Broadcasts all config parameters so object is populated identically on
  !! all processors
  subroutine broadcast_ballstab_config(self)
    use mp, only: broadcast
    implicit none
    class(ballstab_config_type), intent(in out) :: self
    call broadcast(self%beta_div)
    call broadcast(self%beta_mul)
    call broadcast(self%diff)
    call broadcast(self%make_salpha)
    call broadcast(self%n_beta)
    call broadcast(self%n_shat)
    call broadcast(self%shat_max)
    call broadcast(self%shat_min)
    call broadcast(self%theta0)

    call broadcast(self%exist)
  end subroutine broadcast_ballstab_config

  !> Gets the default name for this namelist
  function get_default_name_ballstab_config()
    implicit none
    character(len = CONFIG_MAX_NAME_LEN) :: get_default_name_ballstab_config
    get_default_name_ballstab_config = "ballstab_knobs"
  end function get_default_name_ballstab_config

  !> Gets the default requires index for this namelist
  function get_default_requires_index_ballstab_config()
    implicit none
    logical :: get_default_requires_index_ballstab_config
    get_default_requires_index_ballstab_config = .false.
  end function get_default_requires_index_ballstab_config
end module ballstab