current_calculations.f90

Go to the documentation of this file.

module current_calculations

  use itm_types
  use mod_helena_io

  implicit none

  real(r8), public :: j0      ! constant if alfa constant
  real(r8), public :: cscale  ! used for calculation of cur0

  private

  real(r8), allocatable :: j_tor_hat(:)
  real(r8), allocatable :: j_tor_hat_av(:), I_tor_hat(:)

  real(r8), allocatable, public :: j_tor_loc(:), j_tor_av(:), I_tor(:)

  public :: current_densities
  public :: current_calculations_constructor
  public :: current_calculations_destructor
  public :: total_current, si_currents
  public :: q_calculation

  contains

!--------------------------------------------------------------------------
  subroutine current_calculations_constructor(nr, np)
!--------------------------------------------------------------------------

    implicit none

    integer(itm_i4), intent(in) :: nr, np

    allocate(j_tor_hat((nr - 1) * (np - 1)))
    allocate(j_tor_hat_av(nr))
    allocate(i_tor_hat(nr))
    allocate(j_tor_loc((nr - 1) * (np - 1)))
    allocate(j_tor_av(nr))
    allocate(i_tor(nr))

!-- initialize I_tor_hat
    i_tor_hat = 0._r8

  end subroutine current_calculations_constructor

!--------------------------------------------------------------------------
  subroutine current_calculations_destructor
!--------------------------------------------------------------------------

    implicit none

    deallocate(j_tor_hat, j_tor_hat_av)
    deallocate(j_tor_loc, j_tor_av)
    deallocate(i_tor_hat, i_tor)

  end subroutine current_calculations_destructor

!--------------------------------------------------------------------------
  subroutine current_densities(a, xaxis, type)
!--------------------------------------------------------------------------
! Calculates averaged normalized current densities
! The averaging process depends on the type:
!   type = 'line  ' : line average in poloidal cross section
!        = 'area  ' : surface average over poloidal cross section between
!                     two flux surfaces
!        = 'volume' : volume average between two flux surfaces
!        = 'theta ' : average over surface area of a flux surface
!--------------------------------------------------------------------------
  
    use mod_dat, only : eps
    use mod_mesh, only : nr, np
    use math_functions, only : element_average, flux_surface_average
    use phys_functions, only : rh_side, identity
  
    implicit none
  
    real(r8), intent(in) :: a, xaxis
    character(len = *), intent(in) :: type
  
    real(r8) :: dl(np - 1), df(np - 1)
    integer(itm_i4) :: i, j, n

!-- initialize dl and df
    dl = 0._r8
    df = 0._r8

!-- current density calculation excluding magnetic axis (i = 1)
    n = 0
    do i = 2, nr
      call element_average(i, type, xaxis, a, 0._r8, rh_side, df, dl)
      do j = 1, np - 1
        n = n + 1
!-- normalized current density
        j_tor_hat(n) = -eps * df(j)
      end do
    end do
!-- flux surface averaged current density (including magnetic axis, i = 1)
    do i = 1, nr
      j_tor_hat_av(i) = -eps * flux_surface_average(i, xaxis, a, 0._r8, &
       type, rh_side)
    end do

    return
  end subroutine current_densities

!--------------------------------------------------------------------------
  subroutine total_current(a, xaxis, toroidal_current)
!--------------------------------------------------------------------------
! calculates the normalized current profile I_tor_hat(i) but returns the
! total plasma toroidal current in SI units in toroidal_current
! Attention: recalculates the current densities for averaging type 'area  '
!--------------------------------------------------------------------------

    use mod_dat, only : rvac, bvac, eps
    use mod_mesh, only : nr, np, ias
    use mod_helena_io, only : out_he
    use mod_output
    use mod_parameter
    use math_functions, only : element_average
    use phys_functions, only : identity

    implicit none

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

    real(r8), intent(in) :: a, xaxis
    real(r8), intent(out) :: toroidal_current

    real(r8) :: df(np - 1), da(np - 1)
    real(r8) :: da_sum, i_sum, factas
    integer(itm_i4) :: i, j

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

    call current_densities(a, xaxis, type)

    factas = 2._r8
    if (ias == 1) factas = 1._r8

    i_sum = 0._r8
    do i = 1, nr
      da_sum = 0._r8
      call element_average(i, type, xaxis, 1._r8, 0._r8, identity, df, da)
      do j = 1, np -1
        da_sum = da_sum + da(j)
      end do
      i_sum = i_sum + factas * j_tor_hat_av(i) * da_sum
      i_tor_hat(i) = i_sum
    end do

    toroidal_current = i_tor_hat(nr) * rvac * bvac * eps / (4.e-7_r8 * pi)

    if (standard_output) &
     write(out_he, 4) toroidal_current
 4  format(' total current (alfa = 1) : ', e12.4)
  
    return
  end subroutine total_current

!--------------------------------------------------------------------------
  subroutine si_currents
!--------------------------------------------------------------------------
! denormalizes the currents and current densities
!--------------------------------------------------------------------------

    use mod_parameter
    use mod_dat

    implicit none

    real(r8) :: j0, i0

!-- normalization constants
    j0 = bvac / (mu0 * eps * rvac * alfa)
    i0 = j0 * (eps * rvac)**2

!-- denormalization
    j_tor_loc = j0 * j_tor_hat
    j_tor_av = j0 * j_tor_hat_av

    i_tor = i_tor_hat * i0

  end subroutine si_currents

!--------------------------------------------------------------------------
  subroutine q_calculation(xaxis, cx, cy, a, q)
!--------------------------------------------------------------------------
! calculated the q profile
! requires flux surface aligned coordinates, i.e. prior call to remesh
!--------------------------------------------------------------------------

    use mod_parameter
    use mod_dat
    use mod_gauss
    use mod_interpolation
    use mod_map
    use mod_mesh

    implicit none

    real(r8), intent(in) :: xaxis, cx, cy
    real(r8), intent(in) :: a
    real(r8), intent(inout) :: q(nr)
    real(r8) :: factas
    real(r8) :: r, s
    real(r8) :: sumq
    real(r8) :: x, xr, xs, xrs, xrr, xss
    real(r8) :: y, yr, ys, yrs, yrr, yss
    real(r8) :: ps, psr, pss, psrs, psrr, psss
    real(r8) :: jacobian, bigr, dl, grad_psi2
    integer(itm_i4) :: i, j, ngs
    integer(itm_i4) :: n1, n2, n3, n4

! profiles yields the pressure normalized to eps * Bvac^2 / (mu_0 * alfa^2)
! and R Bphi normalized to Rvac * Bvac
    call profiles(p0, rbphi, dp0, drbphi, a)

    factas = 2._r8
    if (ias == 1) factas = 1._r8

!------------------------------------- nelm elements ------------------
    r = -1._r8
    do i = 1, nr - 1
      sumq = 0._r8
      do j = 1, np - 1
        n1 = (i - 1) * np + j
        n2 = n1 + 1
        n3 = n2 + np
        n4 = n1 + np
!------------------------------------- 4 point gaussian int. in s -----
        do ngs = 1, 4
          s = xgauss(ngs)
          call interpolation(2, xx(1, n1), xx(1, n2), xx(1, n3), &
           xx(1, n4), r, s, x, xr, xs, xrs, xrr, xss)
          call interpolation(2, yy(1, n1), yy(1, n2), yy(1, n3), &
           yy(1, n4), r, s, y, yr, ys, yrs, yrr, yss)
          call interpolation(2, psi(4 * (n1 - 1) + 1), &
           psi(4 * (n2 - 1) + 1), psi(4 * (n3 - 1) + 1), &
           psi(4 * (n4 - 1) + 1), r, s, ps, psr, pss, psrs, psrr, psss)
          jacobian = xr * ys - xs * yr
          bigr = 1._r8 + eps * x
          dl = sqrt(xs**2 + ys**2)
          grad_psi2= psr**2 * (xs**2 + ys**2) / jacobian**2
          sumq = sumq + wgauss(ngs) / (bigr * sqrt(grad_psi2)) * dl
        end do
      end do
      q(nr - i + 1) = 0.5_r8 * factas * sumq * rbphi(nr - i + 1)
    end do
    q(1) = rbphi(1) * pi / (2._r8 * sqrt(cx * cy) * (1. + eps * xaxis))
    do i = 1, nr
      q(i) = q(i) * alfa / pi
    end do

  end subroutine q_calculation

end module current_calculations