math_functions.f90

Go to the documentation of this file.

module math_functions

  use itm_types
  use mod_parameter, only : twopi
  use mod_profiles

  implicit none

  private :: integrate_over_hermite_element, integrate_along_hermite_element
  public :: integrate, root, quad_equation
  public :: element_average, flux_surface_average
  public :: equidistant_profile

  contains

!---------------------------------------------------------------------------
  function integrate(f_spline, x, n, inv) result(g)
!---------------------------------------------------------------------------
! integrates the function f defined by f_spline(1 : n) on the grid points
! x(1 : n)
! returns: g(1:n) with
!          g(i) = integral of f from 1 to i for inv = .false and
!                 from n to i for inv = .true.
!---------------------------------------------------------------------------

    use itm_types
    use helena_spline

    implicit none

    type(spline_coefficients) :: f_spline
    real(r8), dimension(n), intent(in)  :: x
    real(r8), dimension(n) :: g
    integer(itm_i4), intent(in) :: n
    logical :: inv

    real(r8), dimension(3) :: ablt   
    integer(itm_i4) :: i

    g = 0._r8

    do i = 2, n
      if (inv) then
        g(n + 1 - i) = g(n + 2 - i) + (spwert(n, x(n + 2 - i), f_spline, x, &
         ablt, 0) + spwert(n, x(n + 1 - i), f_spline, x, ablt, 0)) &
         * (x(n + 1 - i) - x(n + 2 - i)) / 2._r8
      else
        g(i) = g(i - 1) + (spwert(n, x(i - 1), f_spline, x, ablt, 0) &
         + spwert(n, x(i), f_spline, x, ablt, 0)) * (x(i) - x(i - 1)) &
         / 2._r8
      end if
    end do

  end function integrate

!---------------------------------------------------------------------------
  function integrate_1d(f, x, x_low, x_high, sum) result(integral)
!---------------------------------------------------------------------------
! integrates the function f(1:n) given on grid points x(1:n) from
! x1 to x2
! It assumes that x be monotonically increasing.
! If sum == .true., the integration is done through a simple summation,
! else via linear quadrature.
! returns: integral
!---------------------------------------------------------------------------

    use itm_types

    implicit none

    real(r8), intent(in) :: f(:), x(:)
    real(r8), intent(in) :: x_low, x_high
    logical, intent(in) :: sum

    integer(itm_i4), parameter :: iu6 = 6

    real(r8) :: integral
    real(r8) :: f_x1, f_x2, x1, x2
    integer(itm_i4) :: i, i_x1, i_x2
    logical :: inv

    if (size(f) /= size(x)) then
      write(iu6, *) 'Error in integrate_1d: Array sizes do not agree!'
      stop
    end if

    if (x_low > x_high) then   ! inverted integral
      x1 = x_high
      x2 = x_low
      inv = .true.
    else
      x1 = x_low
      x2 = x_high
      inv = .false.
    end if

!-- find neighbouring nodes
    i_x1 = 1
    i_x2 = 1
    do i = 1, size(x)
      if (x1 >= x(i)) i_x1 = i
      if (x2 >= x(i)) i_x2 = i
    end do
!-- f(i) assigned to interval x(i - 1) .. x(i) for sum == .true.

    integral = 0._r8

    if (.not. x1 == x2) then
!-- linear interpolations for f(x1) and f(x2)
      if (i_x1 == size(x)) then
        f_x1 = f(size(x))
      else
        if (sum) then
          f_x1 = f(i_x1 + 1)
        else
          f_x1 = f(i_x1) + (x1 - x(i_x1)) * (f(i_x1 + 1) - f(i_x1)) &
           / (x(i_x1 + 1) - x(i_x1))
        end if
      end if
      if (i_x2 == size(x)) then
        f_x2 = f(size(x))
      else
        if (sum) then
          f_x2 = f(i_x2 + 1)
        else
          f_x2 = f(i_x2) + (x2 - x(i_x2)) * (f(i_x2 + 1) - f(i_x2)) &
           / (x(i_x2 + 1) - x(i_x2))
        end if
      end if
!-- linear quadrature or simple sum
      if (i_x1 == i_x2) then
        integral = integral + 0.5_r8 * (f_x1 + f_x2) * (x(i_x2) - x(i_x1))
      else
        if (sum) then
          integral = integral + f_x1 * (x(i_x1) - x1)
        else
          integral = integral + 0.5_r8 * (f_x1 + f(i_x1 + 1)) &
           * (x(i_x1 + 1) - x1)
        end if
        do i = i_x1 + 1, i_x2 - 1
          if (sum) then
            integral = integral + f(i) * (x(i) - x(i - 1))
          else
            integral = integral + 0.5_r8 * (f(i) + f(i + 1)) &
             * (x(i + 1) - x(i))
          end if
        end do
        if (sum) then
          integral = integral + f(i_x2) * (x(i_x2) - x(i_x2 - 1))
          integral = integral + f_x2 * (x2 - x(i_x2))
        else
          integral = integral + 0.5_r8 * (f(i_x2) + f_x2) &
           * (x2 - x(i_x2))
        end if
      end if
    end if

    if (inv) integral = -integral

  end function integrate_1d

!---------------------------------------------------------------------------
  subroutine equidistant_profile(f, x, f_equi, x_equi)
!---------------------------------------------------------------------------
! This subroutine takes a function f(x), defined on the array x (not
! necessarily equidistant), interpolates it and returns the
! function f_equi(x_equi) on an equidistant grid x_equi spanning the range
! of x
! requires that x be monotonically increasing
!---------------------------------------------------------------------------
    use mod_profiles
    use helena_spline

    implicit none

    real(r8), intent(in) :: f(:)
    real(r8), intent(inout) :: x(:)
    real(r8), intent(out) :: f_equi(:), x_equi(:)

    type(spline_coefficients) :: f_spline
    real(r8) :: abltg(3)
    integer(itm_i4) :: i

    if (size(f) /= size(x) .or. size(f_equi) /= size(x_equi)) then
      stop 'Error in equidistant_profile: Array sizes do not agree!'
    end if

    call allocate_spline_coefficients(f_spline, size(f))

    call spline(size(f), x, f, 0._r8, 0._r8, 2, f_spline)

    do i = 1, size(f_equi)
      x_equi(i) = x(1) + dble(i - 1) / dble(size(f_equi) - 1) &
       * (x(size(x)) - x(1))
      f_equi(i) = spwert(size(f), x_equi(i), f_spline, x, abltg, 0)
    end do
  
    call deallocate_spline_coefficients(f_spline)

  end subroutine equidistant_profile

!---------------------------------------------------------------------------
  function root(a, b, c, d, sgn) result(f_root)
!---------------------------------------------------------------------------
! this function gives better roots of quadratics by avoiding
! cancellation of smaller root
!---------------------------------------------------------------------------

    implicit none

    real(r8) :: a, b, c, d, sgn
    real(r8) :: f_root

    if (b * sgn >= 0.0) then
      f_root = -2.0_r8 * c / (b + sgn * sqrt(d))
    else
      f_root = (-b + sgn * sqrt(d)) / (2.0_r8 * a)
    end if

  end function root

!---------------------------------------------------------------------------
  function integrate_over_hermite_element(i, j, a, F_dia, func, volume) &
   result(f_integral)
!---------------------------------------------------------------------------
! This function integrates the quantity described by the function func over
! the Hermite element with index i,j.
! volume = .true.  : volume integral
!        = .false. : surface integral
!---------------------------------------------------------------------------
    use mod_dat, only : eps
    use mod_gauss
    use mod_interpolation
    use mod_mesh
    use matrix_calculations
    use mod_parameter, only : twopi

    implicit none

    integer(itm_i4), intent(in) :: i, j
    real(r8), intent(in) :: a, f_dia
    logical :: volume

    real(r8) :: f_integral

    real(r8) :: r, wr, s, ws
    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) :: xjac, major_radius
    integer(itm_i4) :: ngr, ngs
    integer(itm_i4) :: nelm
    integer(itm_i4) :: n1, n2, n3, n4

    real(r8), external:: func

    f_integral = 0._r8

    nelm = (i - 1) * (np - 1) + j      
    call set_node_number(nelm, n1, n2, n3, n4)
!---------------------------------- 4 point gaussian integration in r -----
    do ngr = 1, 4
      r = xgauss(ngr)
      wr = wgauss(ngr)
!---------------------------------- 4 point gaussian integration in s -----
      do ngs = 1, 4
        s = xgauss(ngs)
        ws = wgauss(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)
        xjac = xr * ys - xs * yr
        if (volume) then
          major_radius = 1._r8 + eps * x
          f_integral = f_integral + wr * ws * xjac * major_radius &
           * twopi * func(a, x, xr, xs, yr, ys, ps, psr, f_dia)
        else
          f_integral = f_integral + wr * ws * xjac &
           * func(a, x, xr, xs, yr, ys, ps, psr, f_dia)
        end if
      end do
    end do

  end function integrate_over_hermite_element

!---------------------------------------------------------------------------
  function integrate_along_hermite_element(i, j, a, F_dia, func, theta) &
   result(f_integral)
!---------------------------------------------------------------------------
! This function integrates the quantity described by the function func along
! the Hermite element with index i,j.
! theta = .true.  : surface integral (over flux surface)
!       = .false. : line integral (along flux surface)
!---------------------------------------------------------------------------
    use mod_dat, only : eps
    use mod_gauss
    use mod_interpolation
    use mod_mesh
    use matrix_calculations
    use mod_parameter, only : twopi

    implicit none

    integer(itm_i4), intent(in) :: i, j
    real(r8), intent(in) :: a, f_dia
    logical :: theta

    real(r8) :: f_integral

    real(r8) :: r, s, ws
    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) :: xjac, dl, major_radius
    integer(itm_i4) :: ngs
    integer(itm_i4) :: nelm
    integer(itm_i4) :: n1, n2, n3, n4

    real(r8), external:: func

    f_integral = 0._r8

    r = -1._r8
    nelm = (i - 1) * (np - 1) + j      
    call set_node_number(nelm, n1, n2, n3, n4)
!---------------------------------- 4 point gaussian integration in s -----
    do ngs = 1, 4
      s = xgauss(ngs)
      ws = wgauss(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)
      if (theta) then
        xjac = xr * ys - xs * yr
        major_radius = 1._r8 + eps * x
!        if (xjac /= 0._r8) then
          f_integral = f_integral + ws * xjac * major_radius / abs(psr) &
           * twopi * func(a, x, xr, xs, yr, ys, ps, psr, f_dia)
!        end if
      else
        dl = sqrt(xs**2 + ys**2)
        f_integral = f_integral + ws * dl &
         * func(a, x, xr, xs, yr, ys, ps, psr, f_dia)
      end if
    end do

  end function integrate_along_hermite_element

!---------------------------------------------------------------------------
  subroutine element_average(i, type, xaxis, a, F_dia, func, df, dl)
!---------------------------------------------------------------------------
! This subroutine averages a quantity given by function func over the
! Hermite elements with radial index i and returns the result in df plus
! the corresponding line/area/volume elements in dl.
! 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
!
! Note that df and dl are ordered from axis to separatrix, i.e., psi is
! not inverted!
! xaxis is needed to extend the average to the magnetic axis (psi = 0)
! The minus sign for the types 'area' and 'volume' account for the
! fact that the Jacobian determinant XJAC is negative due to the orientation
! of r (from separatrix towards axis)!
!---------------------------------------------------------------------------
    use mod_helena_io, only : iu6
    use mod_mesh, only : nr, np
    use phys_functions, only : identity

    implicit none

    character(len = 6), intent(in) :: type
    real(r8), intent(in) :: xaxis, a, f_dia
    integer(itm_i4), intent(in) :: i

    real(r8) :: dl(np - 1), df(np - 1)

    real(r8) :: one = 1._r8, zero = 0._r8
    integer(itm_i4) :: j

    real(r8), external :: func

!-- no flux surface average possible on axis since dl = 0
    if (i == 1) then
      dl = 0._r8
      do j = 1, np - 1
        df(j) = func(a, xaxis, zero, zero, zero, zero, zero, zero, f_dia)
      end do
    else
      do j = 1, np - 1
!-- elements in HELENA are ordered in inverse radial direction, i.e.,
!   psi(nr - 1) = 1, psi(1) = 0
        select case (trim(type))
          case ('line')
            dl(j) = integrate_along_hermite_element(nr + 1 - i, j, a, &
             f_dia, identity, theta = .false.)
            df(j) = integrate_along_hermite_element(nr + 1 - i, j, a, &
             f_dia, func, theta = .false.)
          case ('area')
            dl(j) = -integrate_over_hermite_element(nr + 1 - i, j, a, &
             f_dia, identity, volume = .false.)
            df(j) = -integrate_over_hermite_element(nr + 1 - i, j, a, &
             f_dia, func, volume = .false.)
          case ('volume')
            dl(j) = -integrate_over_hermite_element(nr + 1 - i, j, a, &
             f_dia, identity, volume = .true.)
            df(j) = -integrate_over_hermite_element(nr + 1 - i, j, a, &
             f_dia, func, volume = .true.)
          case ('theta')
            dl(j) = integrate_along_hermite_element(nr + 1 - i, j, a, &
             f_dia, identity, theta = .true.)
            df(j) = integrate_along_hermite_element(nr + 1 - i, j, a, &
             f_dia, func, theta = .true.)
          case default
            write(iu6, *) 'ERROR: wrong type of averaging parameter!'
            stop
        end select
        df(j) = df(j) / dl(j)
      end do
    end if

  end subroutine element_average

!---------------------------------------------------------------------------
  function flux_surface_average(i, xaxis, a, F_dia, type, func) &
   result(f_average)
!---------------------------------------------------------------------------
! This function averages a quantity given by function func over the flux
! surface with radial index i and returns the result in f_average.
! 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_mesh, only : np

    implicit none

    character(len = 6), intent(in) :: type
    real(r8), intent(in) :: xaxis, a, f_dia
    integer(itm_i4), intent(in) :: i

    real(r8) :: f_average

    real(r8) :: dl(np - 1), df(np - 1)
    real(r8) :: dl_total
    integer(itm_i4) :: j

    real(r8), external:: func

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

    dl_total = 0._r8
    f_average = 0._r8

    call element_average(i, type, xaxis, a, f_dia, func, df, dl)

    do j = 1, np - 1
      dl_total = dl_total + dl(j)
      f_average = f_average + df(j) * dl(j)
    end do

!-- average on axis = original value on axis
    if (i == 1) then
      f_average = df(1)
    else
!-- average over flux surface i
      f_average = f_average / dl_total
    end if

  end function flux_surface_average

!---------------------------------------------------------------------------
  subroutine quad_equation(p_m, p_mr, p_p, p_pr, tol, r, ifail)
!---------------------------------------------------------------------------
! this subroutine solves a quadratic equation
!---------------------------------------------------------------------------

    implicit none

    real(r8), intent(in) :: p_m, p_mr, p_p, p_pr, tol
    real(r8), intent(out) :: r
    integer(itm_i4), intent(out) :: ifail

    real(r8) :: a, b, c, d

    ifail = 0
    a = 3._r8 * (p_m + p_mr - p_p + p_pr ) / 4._r8
    b = (-p_mr + p_pr ) / 2._r8
    c = (-3._r8 * p_m - p_mr + 3._r8 * p_p - p_pr) / 4._r8
    d = b * b - 4._r8 * a * c

    if (d > 0.) then
      r  = root(a, b, c, d, 1._r8)
      if (abs(r) > tol) r = root(a, b, c, d, -1._r8)
    else
      r = 999._r8
      ifail = 1
    end if

  return
  end subroutine quad_equation

end module math_functions