global_plasma_parameters.f90

Go to the documentation of this file.

module global_plasma_parameters

  use itm_types
  use mod_mesh, only : nr, np, ias
  use mod_parameter
  use euitm_schemas

  implicit none

  public :: w_mhd, plasma_volume, cross_section, surface_area

  contains

!-----------------------------------------------------------------------------
  function w_mhd(equilibrium, rminor, eps) result(f_w_mhd)
!-----------------------------------------------------------------------------
! This function calculates W_MHD, i.e. the volume integral of the pressure
! times 1.5 in SI units.
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : identity

    implicit none

    type (type_equilibrium), intent(in) :: equilibrium
    real(r8) :: rminor, eps
    real(r8) :: f_w_mhd

    character(len = 6), parameter :: type = 'volume'

    real(r8) :: df(np - 1), dv(np - 1), dv_flux_tube(nr)
    real(r8) :: xaxis, factas
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- normalized magnetic axis
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor

!-- initialize dV and df
    dv = 0._r8
    df = 0._r8

!-- calculate volume of flux tubes
    do i = 1, nr
      dv_flux_tube(i) = 0._r8
      call element_average(i, type, xaxis, 1._r8, 0._r8, identity, df, dv)
      do j = 1, np -1
        dv_flux_tube(i) = dv_flux_tube(i) + dv(j)
      end do
    end do
    dv_flux_tube = dv_flux_tube * factas  ! symmetry factor
    dv_flux_tube = dv_flux_tube * rminor**3 / eps  ! denormalize

!-- interpolate pressure p for centers of flux tubes
!   p((x(i) + x(i + 1))/2) = p(x(i)) + (x(i + 1) - x(i)) / 2 &
!    * (0.75 * dp(x(i)) + 0.25 * dp(x(i + 1))) 
!   dV_flux_tube(i) = volume of flux tube between flux surfaces i - 1 and i
!   i.e. dV_flux_tube(1) = 0 (magnetic axis)
    f_w_mhd = 0._r8
    do i = 2, nr
      f_w_mhd = f_w_mhd + 1.5_r8 * dv_flux_tube(i) &
       * (equilibrium%profiles_1d%pressure(i - 1) &
       + 0.5_r8 * (equilibrium%profiles_1d%psi(i) &
       - equilibrium%profiles_1d%psi(i - 1)) * (0.75_r8 &
       * equilibrium%profiles_1d%pprime(i - 1) + 0.25_r8 &
       * equilibrium%profiles_1d%pprime(i)))
    end do

  end function w_mhd

!-----------------------------------------------------------------------------
  function plasma_volume(equilibrium, rminor, eps) &
   result(f_plasma_volume)
!-----------------------------------------------------------------------------
! This function calculates the volume enclosed by the flux surfaces
! (-> equilibrium) and returns the plasma volume in SI units.
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : identity

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: rminor, eps
    real(r8) :: f_plasma_volume

    character(len = 6), parameter :: type = 'volume'

    real(r8) :: df(np - 1), dv(np - 1), dv_flux_tube(nr)
    real(r8) :: xaxis, factas
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor

!-- initialize dV and df
    dv = 0._r8
    df = 0._r8

!-- calculate volume of flux tubes
    do i = 1, nr
      dv_flux_tube(i) = 0._r8
      call element_average(i, type, xaxis, 1._r8, 0._r8, identity, df, dv)
      do j = 1, np -1
        dv_flux_tube(i) = dv_flux_tube(i) + dv(j)
      end do
    end do
    dv_flux_tube = dv_flux_tube * factas  ! symmetry factor
    dv_flux_tube = dv_flux_tube * rminor**3 / eps  ! denormalize

    f_plasma_volume = 0._r8
    do i = 1, nr
      f_plasma_volume = f_plasma_volume + dv_flux_tube(i)
      equilibrium%profiles_1d%volume(i) = f_plasma_volume
    end do

  end function plasma_volume

!-----------------------------------------------------------------------------
  function cross_section(equilibrium, rminor, eps) &
   result(f_cross_section)
!-----------------------------------------------------------------------------
! This function calculates the area enclosed by the flux surfaces
! (-> equilibrium) and returns the total area of the plasma
! cross-section in SI units.
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : identity

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: rminor, eps
    real(r8) :: f_cross_section

    character(len = 6), parameter :: type = 'area  '

    real(r8) :: df(np - 1), da(np - 1), da_flux_tube(nr)
    real(r8) :: xaxis, factas
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor

!-- initialize dA and df
    da = 0._r8
    df = 0._r8

!-- calculate area of flux tubes
    do i = 1, nr
      da_flux_tube(i) = 0._r8
      call element_average(i, type, xaxis, 1._r8, 0._r8, identity, df, da)
      do j = 1, np -1
        da_flux_tube(i) = da_flux_tube(i) + da(j)
      end do
    end do
    da_flux_tube = da_flux_tube * factas     ! symmetry factor
    da_flux_tube = da_flux_tube * rminor**2  ! denormalize

    f_cross_section = 0._r8
    do i = 1, nr
      f_cross_section = f_cross_section + da_flux_tube(i)
      equilibrium%profiles_1d%area(i) = f_cross_section
    end do

  end function cross_section

!-----------------------------------------------------------------------------
  subroutine surface_area(equilibrium, rminor, eps)
!-----------------------------------------------------------------------------
! This subroutine calculates the surface area of all flux surfaces
! (-> equilibrium) in SI units.
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : r_major

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: rminor, eps

    character(len = 6), parameter :: type = 'line  '

    real(r8) :: df(np - 1), da(np - 1), da_flux_tube(nr)
    real(r8) :: xaxis, factas
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor

!-- initialize dA and df
    da = 0._r8
    df = 0._r8

!-- calculate surface area of flux surfaces
    do i = 1, nr
      da_flux_tube(i) = 0._r8
      call element_average(i, type, xaxis, 1._r8, 0._r8, r_major, df, da)
      do j = 1, np -1
        da_flux_tube(i) = da_flux_tube(i) + da(j) * df(j)
      end do
    end do
    da_flux_tube = da_flux_tube * factas             ! symmetry factor
    da_flux_tube = da_flux_tube * rminor &
     * equilibrium%eqgeometry%geom_axis%r * twopi    ! denormalize

    equilibrium%profiles_1d%surface = da_flux_tube

  end subroutine surface_area

!-----------------------------------------------------------------------------
  function beta_poloidal(a, equilibrium) &
   result(f_beta_poloidal)
!-----------------------------------------------------------------------------
! This function calculates the poloidal beta for all flux surfaces
! (-> equilibrium) and returns the poloidal beta of the equilibrium
! dependencies: profiles_1d%area
!               p_spline (profiles has to be called beforehand)
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : identity, p
    use mod_parameter, only : mu0
    use mod_dat

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: rminor
    real(r8) :: f_beta_poloidal

    character(len = 6), parameter :: type = 'area  '

    real(r8) :: df(np - 1), da(np - 1), dl(np - 1)
    real(r8), intent(in) :: a
    real(r8) :: factas, xaxis
    real(r8) :: lp, b_a
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor
!-- rminor = a_minor
    rminor = equilibrium%eqgeometry%a_minor

!-- initialize dA and df
    da = 0._r8
    df = 0._r8
    dl = 0._r8
    f_beta_poloidal = 0._r8

!-- calculate poloidal beta
    do i = 1, nr
      call element_average(i, type, xaxis, a, 0._r8, p, df, da)
      da = da * factas     ! symmetry factor
      da = da * rminor**2  ! denormalize
      do j = 1, np -1
        f_beta_poloidal = f_beta_poloidal + df(j) * da(j)
      end do
!-- average on axis = original value on axis
      if (i == 1) then
        equilibrium%profiles_1d%beta_pol(i) = df(1)
      else
!-- average over cross-section up to flux surface i
        equilibrium%profiles_1d%beta_pol(i) = f_beta_poloidal &
         / equilibrium%profiles_1d%area(i)
      end if
    end do

!-- calculate B_a
    call element_average(nr, 'line  ', xaxis, 1._r8, 0._r8, identity, df, dl)
    lp = 0._r8
    do j = 1, np -1
      lp = lp + dl(j)
    end do
    lp = lp * factas * rminor
    b_a = mu0 * equilibrium%global_param%i_plasma / lp

    equilibrium%profiles_1d%beta_pol = equilibrium%profiles_1d%beta_pol &
     * 2._r8 * mu0 / b_a**2

    f_beta_poloidal = equilibrium%profiles_1d%beta_pol(nr)

  end function beta_poloidal

!-----------------------------------------------------------------------------
  function beta_toroidal(a, equilibrium) &
   result(f_beta_toroidal)
!-----------------------------------------------------------------------------
! This function returns the toroidal beta for the last closed flux surface
! dependencies: global_param%volume, global_param%toroid_field%b0
!               p_spline (profiles has to be called beforehand)
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : p
    use mod_parameter, only : mu0
    use mod_dat

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: rminor
    real(r8) :: f_beta_toroidal

    character(len = 6), parameter :: type = 'volume'

    real(r8) :: df(np - 1), dv(np - 1)
    real(r8), intent(in) :: a
    real(r8) :: factas, xaxis
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor
!-- rminor = a_minor
    rminor = equilibrium%eqgeometry%a_minor

!-- initialize dA and df
    dv = 0._r8
    df = 0._r8
    f_beta_toroidal = 0._r8

!-- calculate toroidal beta (no contribution from axis
!   since dV = 0)
    do i = 2, nr
      call element_average(i, type, xaxis, a, 0._r8, p, df, dv)
      dv = dv * factas           ! symmetry factor
      dv = dv * rminor**3 / eps  ! denormalize
      do j = 1, np -1
        f_beta_toroidal = f_beta_toroidal + df(j) * dv(j)
      end do
    end do
!-- volume average
    f_beta_toroidal = f_beta_toroidal / equilibrium%global_param%volume
!-- normalization
    f_beta_toroidal = f_beta_toroidal * 2._r8 * mu0 &
     / equilibrium%global_param%toroid_field%b0**2

  end function beta_toroidal

!-----------------------------------------------------------------------------
  function internal_inductance(a, equilibrium) &
   result(f_internal_inductance)
!-----------------------------------------------------------------------------
! This function calculates the normalized internal inductance li for all flux
! surfaces (-> equilibrium) and returns the normalized internal inductance
! of the equilibrium
! dependencies: profiles_1d%area
!               profiles_1d%jphi
!               profiles_1d%psi
!-----------------------------------------------------------------------------

    use math_functions, only : element_average
    use phys_functions, only : bp2
    use mod_parameter, only : mu0
    use mod_dat
    use mod_profiles
    use helena_spline

    implicit none

    type (type_equilibrium), intent(inout) :: equilibrium
    real(r8) :: f_internal_inductance

    character(len = 6), parameter :: type = 'volume'

    type(spline_coefficients) :: coeff 
    real(r8) :: abltg(3)
    real(r8) :: df(np - 1), dv(np - 1), dl(np - 1)
    real(r8) :: i_tor(nr)
    real(r8), intent(in) :: a
    real(r8) :: factas, xaxis
    real(r8) :: rminor, r_average
    integer(itm_i4) :: i, j

    if (ias == 1) then
      factas = 1._r8
    else
      factas = 2._r8
    end if

!-- xaxis = (Rmag - Rvac) / a
    xaxis = (equilibrium%global_param%mag_axis%position%r &
     - equilibrium%eqgeometry%geom_axis%r) &
     / equilibrium%eqgeometry%a_minor
    rminor = equilibrium%eqgeometry%a_minor

!-- initialize dV and df
    dv = 0._r8
    df = 0._r8
    i_tor = 0._r8
    f_internal_inductance = 0._r8

!-- calculate internal inductance
    do i = 2, nr
      call element_average(i, type, xaxis, a, 0._r8, bp2, df, dv)
      dv = dv * factas           ! symmetry factor
      dv = dv * rminor**3 / eps  ! denormalize
      do j = 1, np -1
        f_internal_inductance = f_internal_inductance + df(j) * dv(j)
      end do
      i_tor(i) = i_tor(i - 1) + (equilibrium%profiles_1d%area(i) &
       - equilibrium%profiles_1d%area(i - 1)) &
       * equilibrium%profiles_1d%jphi(i)
      r_average = equilibrium%profiles_1d%volume(i) / (twopi &
       * equilibrium%profiles_1d%area(i))
      equilibrium%profiles_1d%li(i) = f_internal_inductance / i_tor(i)**2 &
       / (equilibrium%eqgeometry%a_minor &
       * equilibrium%eqgeometry%geom_axis%r &
       / (equilibrium%profiles_1d%psi(nr) &
       - equilibrium%profiles_1d%psi(1)))**2 / mu0 
!-- normalize internal inductance
      equilibrium%profiles_1d%li(i) = 2._r8 * equilibrium%profiles_1d%li(i) &
       / (mu0 * r_average)
    end do

!-- extrapolate internal inductance onto axis (over volume)
    call allocate_spline_coefficients(coeff, nr - 1)
    call spline(nr - 1, equilibrium%profiles_1d%volume(2 : nr), &
     equilibrium%profiles_1d%li(2 : nr), 0._r8, 0._r8, 2, coeff)
    equilibrium%profiles_1d%li(1) = spwert(nr - 1, 0._r8, coeff, &
     equilibrium%profiles_1d%volume(2 : nr), abltg, 0)
    call deallocate_spline_coefficients(coeff)

    f_internal_inductance = equilibrium%profiles_1d%li(nr)

  end function internal_inductance

end module global_plasma_parameters