p_j_iteration.f90

Go to the documentation of this file.

subroutine p_j_iteration(equil_in, equil_out, first, a, cx, cy, psaxis, &
 xaxis, yaxis, nax, rax, sax)
!-----------------------------------------------------------------------------
! subroutine to solve the Grad-Shafranov equation for either p and j_tor
! or dp and j_tor input
! The following input combinations are allowed so far:
! 1.) input_type = "p and j_tor", radial_coordinate = "psi":
!     HELENA will solve the GS equation for dp/dpsi(psi) and j_tor(psi) after
!     deriving dp/dpsi(psi) from p(psi) and psi specified in equil_in;
!     whether psi is maintained as in equil_in or rather Ip depends on the
!     specification of either alfa or Ip
! 2.) input_type = "p and j_tor", radial_coordinate = "rho_vol"
!     HELENA will solve the GS equation for dp/dpsi(psi) and j_tor(psi) by
!     updating them on every iteration step. I.e. for a current solution with
!     psi and rho_vol, psi is transferred to the rho_vol grid of equil_in
!     and the input p and j_tor are interpolated and derived with respect to
!     this transferred psi grid and then evaluated at the current psi knots.
!     dp/dpsi(psi) and j_tor(psi) therefore change with every iteration step.
! 3.) input_type = "p' and j_tor", radial_coordinate = "psi"
!     HELENA will solve the GS equation for dp/dpsi(psi) and j_tor(psi) as
!     specified in equil_in.
! 4.) input_type = "p' and FF'", radial_coordinate = "psi"
!     HELENA solves the GS equation by a simple call to dp_fdf_iteration
!     (no current iteration needed)
! Attention: still assumes that equil_in%profiles_1d%psi be set on input for
!            all cases (will be lifted later)
! Subroutine assumes that equil_out be allocated on entry!
!-----------------------------------------------------------------------------

  use itm_types
  use mod_dat
  use mod_equilibrium
  use global_plasma_parameters
  use mod_helena_io, only : iu6, out_he
  use mod_interpolation
  use mod_mesh
  use mod_output
  use mod_parameter, only : mu0
  use mod_profiles
  use current_calculations
  use math_functions, only : integrate
  use euitm_schemas
  use plot_data
  use helena_spline

  implicit none

  type(type_equilibrium), intent(inout) :: equil_in, equil_out
  type(spline_coefficients) :: p_spline, coeff

  logical, intent(inout) :: first

  real(r8), intent(inout) :: a
  real(r8), intent(inout) :: xaxis, yaxis
  real(r8), intent(inout) :: cx, cy
  real(r8), intent(inout) :: rax, sax

  real(r8), intent(out) :: psaxis
  integer(itm_i4), intent(out) :: nax


  real(r8), parameter :: fdf_eps = 1.e-7_r8

  character(len = 2) :: file_no
  real(r8), dimension(nr) :: psi_hat
  real(r8), dimension(size(equil_in%profiles_1d%psi)) :: rho_sqr_in
  real(r8), dimension(size(equil_out%profiles_1d%psi)) :: rho_sqr_out
  real(r8), dimension(npts) :: psi_tmp
  real(r8), dimension(npts) :: fdf_old
  real(r8), dimension(3) :: abltg
  real(r8) :: pressure
  real(r8) :: diff_fdf
  real(r8) :: b_ref
  integer(itm_i4) :: i, j
  integer(itm_i4) :: meshno
  logical :: converged_current

!-- initializations
  converged_current = .false.
  meshno = 0

!-- prepare and check input for "p and j_tor" cases
  if (trim(input_type) == "p and j_tor") then
    if (.not. associated(equil_in%profiles_1d%pressure)) then
      write(iu6, *) 'ERROR in p_j_iteration: pressure not provided ', &
       'on input.'
      stop
    end if
    if (.not. associated(equil_in%profiles_1d%jphi)) then
      write(iu6, *) 'ERROR in p_j_iteration: jphi not provided on input.'
      stop
    end if

    if (trim(radial_coordinate) == "psi") then

      if (.not. associated(equil_in%profiles_1d%psi)) then
        write(iu6, *) 'ERROR in p_j_iteration: psi not provided ', &
         'on input.'
        stop
      end if

    else if (trim(radial_coordinate) == "rho_vol") then

      if (.not. associated(equil_in%profiles_1d%rho_vol)) then
        write(iu6, *) 'ERROR in p_j_iteration: rho_vol not provided ', &
         'on input.'
        stop
      end if

!-- avoid singular derivative at axis by interpolating rho**2
      rho_sqr_in = equil_in%profiles_1d%rho_vol**2
      call transfer_1d_profile(equil_in%profiles_1d%psi, &
       rho_sqr_in, equil_out%profiles_1d%psi, &
       rho_sqr_out, spline_type = 1, monotonic = .true.)
      equil_out%profiles_1d%rho_vol = rho_sqr_out**0.5

    else if (trim(radial_coordinate) == "rho_tor") then

      if (.not. associated(equil_in%profiles_1d%rho_tor)) then
        write(iu6, *) 'ERROR in p_j_iteration: rho_tor not provided ', &
         'on input.'
        stop
      end if

!-- avoid singular derivative at axis by interpolating rho**2
      rho_sqr_in = equil_in%profiles_1d%rho_tor**2
      call transfer_1d_profile(equil_in%profiles_1d%psi, &
       rho_sqr_in, equil_out%profiles_1d%psi, &
       rho_sqr_out, spline_type = 1, monotonic = .true.)
      equil_out%profiles_1d%rho_tor = rho_sqr_out**0.5

    end if

!-- prepare and check input for "p' and j_tor" cases
  else if (trim(input_type) == "p' and j_tor") then

    if (.not. associated(equil_in%profiles_1d%pprime)) then
      write(iu6, *) 'ERROR in p_j_iteration: pprime not provided ', &
       'on input.'
      stop
    end if
    if (.not. associated(equil_in%profiles_1d%jphi)) then
      write(iu6, *) 'ERROR in p_j_iteration: jphi not provided on input.'
      stop
    end if

    call transfer_1d_profile(equil_in%profiles_1d%psi, &
     equil_in%profiles_1d%pprime, equil_out%profiles_1d%psi, &
     equil_out%profiles_1d%pprime, spline_type = 1)

    call transfer_1d_profile(equil_in%profiles_1d%psi, &
     equil_in%profiles_1d%jphi, equil_out%profiles_1d%psi, &
     equil_out%profiles_1d%jphi, spline_type = 1)

  end if

!-- beginning of current loop
  do i = 1, nmesh
    if (verbosity > 3) then
      write(file_no, '(i2.2)') i
    end if

    if (verbosity > 1) then
      write(iu6, *) 'input_type ', input_type
      write(iu6, *) 'radial_coordinate ', radial_coordinate
    end if

    if (trim(input_type) == "p and j_tor") then 
      if (trim(radial_coordinate) == "rho_vol") then
!-- interpolate new psi (equil_out) on new rho_vol, to get new psi
!   on reference rho_vol (equil_in)
!-- avoid singular derivative at axis by interpolating rho**2
        rho_sqr_in = equil_in%profiles_1d%rho_vol**2
        rho_sqr_out = equil_out%profiles_1d%rho_vol**2
        call transfer_1d_profile(rho_sqr_out, &
         equil_out%profiles_1d%psi, rho_sqr_in, &
         psi_tmp, spline_type = 1, monotonic = .true.)

        if (verbosity > 2) then
          write(iu6, *) 'old cpsurfin = ', &
           equil_in%profiles_1d%psi(size(equil_in%profiles_1d%psi))
          write(iu6, *) 'new cpsurfin = ', psi_tmp(size(psi_tmp))
        end if

      else if (trim(radial_coordinate) == "rho_tor") then
!-- interpolate new psi (equil_out) on new rho_tor, to get new psi
!   on reference rho_tor (equil_in)
!-- avoid singular derivative at axis by interpolating rho**2
        rho_sqr_in = (equil_in%profiles_1d%rho_tor &
         / equil_in%profiles_1d%rho_tor(nr))**2
        rho_sqr_out = (equil_out%profiles_1d%rho_tor &
         / equil_out%profiles_1d%rho_tor(nr))**2
        call transfer_1d_profile(rho_sqr_out, &
         equil_out%profiles_1d%psi, rho_sqr_in, &
         psi_tmp, spline_type = 1, monotonic = .true.)

        if (verbosity > 2) then
          write(iu6, *) 'old cpsurfin = ', &
           equil_in%profiles_1d%psi(size(equil_in%profiles_1d%psi))
          write(iu6, *) 'new cpsurfin = ', psi_tmp(size(psi_tmp))
          write(iu6, *) 'old rho_tor(nr) = ', &
           equil_in%profiles_1d%rho_tor(size(equil_in%profiles_1d%rho_tor))
          write(iu6, *) 'new rho_tor(nr) = ', &
           equil_out%profiles_1d%rho_tor(size(equil_out%profiles_1d%rho_tor))
        end if


      else if (trim(radial_coordinate) == "psi") then
!-- maintain p(psi_bar) and j_tor(psi_bar) with psi_bar being the normalized
!   psi specified on input
        psi_tmp = (equil_in%profiles_1d%psi &
         / equil_in%profiles_1d%psi(size(equil_in%profiles_1d%psi))) &
         * equil_out%profiles_1d%psi(size(equil_out%profiles_1d%psi))
      end if

      if (i <= 50) then
!-- calculate dp/dpsi for new psi (equil_out) from reference
!   p given on reference rho_vol
        call transfer_1d_profile(psi_tmp, equil_in%profiles_1d%pressure, &
         equil_out%profiles_1d%psi, equil_out%profiles_1d%pressure, &
         equil_out%profiles_1d%pprime, spline_type = 3)
      else
!-- calculate dp/dpsi for new psi (equil_out) from reference
!   pprime given on reference rho_vol
        call transfer_1d_profile(equil_in%profiles_1d%psi, &
         equil_in%profiles_1d%pprime, equil_out%profiles_1d%psi, &
         equil_out%profiles_1d%pprime, spline_type = 3)
      end if
      if (xmgrace_output .and. verbosity > 3) then
       call plot_data_1d('line', equil_in%profiles_1d%psi, &
        equil_in%profiles_1d%pressure, size(equil_in%profiles_1d%psi), &
        'xmgr_p_in_' // file_no)
       call plot_data_1d('line', equil_out%profiles_1d%psi, &
        equil_out%profiles_1d%pressure, size(equil_out%profiles_1d%psi), &
        'xmgr_p_out_' // file_no)
       call plot_data_1d('line', equil_out%profiles_1d%psi, &
        equil_out%profiles_1d%pprime, size(equil_out%profiles_1d%psi), &
        'xmgr_pprime_out_' // file_no)
      end if

!-- set new reference profiles after 5 iterations
      if (i == 50) then
        equil_in%profiles_1d%psi = psi_tmp
        call transfer_1d_profile(equil_out%profiles_1d%psi, &
         equil_out%profiles_1d%pressure, psi_tmp, &
         equil_in%profiles_1d%pressure, equil_in%profiles_1d%pprime, &
         spline_type = 3)
      end if
!-- calculate jphi for new psi (equil_out) from reference
!   jphi given on reference grid (rho_vol or rho_tor)
      call transfer_1d_profile(psi_tmp, equil_in%profiles_1d%jphi, &
       equil_out%profiles_1d%psi, equil_out%profiles_1d%jphi, spline_type = 3)

!TODO: The following should be one contingent block (no call to
!      fill_equilibrium allowed in this block!
!-- calculate dpres on non-equidistant psivec grid with npts points
      call deallocate_spline_coefficients(dp_spline)
      call allocate_spline_coefficients(dp_spline, nr)

!-- splining has to on psi_hat to be able to transfer to psivec which
!   goes from 0 to 1 !
      psi_hat = equil_out%profiles_1d%psi / equil_out%profiles_1d%psi(nr)
      call spline(nr, psi_hat, equil_out%profiles_1d%pprime, 0._r8, &
       0._r8, 2, dp_spline)

!TODO: Is psi_hat really different from psivec? Isn't npts=nr?
      if (verbosity > 3) then
        write(iu6, *) 'npts = ', npts, ' nr = ', nr
      end if
!-- transfer equil_out%profiles_1d%pprime to psivec
      do j = 1, npts
        dpres(j) = spwert(nr, psivec(j), dp_spline, psi_hat, abltg, 0)
      end do
      call deallocate_spline_coefficients(dp_spline)

!-- integrate to get new pressure
      call allocate_spline_coefficients(dp_spline, npts)
!TODO: Is this extra splining necessary (not if psivec=psi_hat)?
      call spline(npts, psivec, dpres, 0._r8, 0._r8, 2, dp_spline)
      p_int = integrate(dp_spline, psivec, npts, .true.) &
       * sign(1._r8, cpsurfin)
      if (verbosity >= 2) &
       write(iu6, *) 'new pressure on axis ', p_int(1)
    end if

!-- update FF' for input_type /= "p' and FF'"
    if (trim(input_type) == "p' and j_tor" &
     .or. trim(input_type) == "p and j_tor" &
     .or. trim(input_type) == "p' and q") then

      if (i > 1) then
        fdf_old = fdf
      else
        fdf_old = 0._r8
      end if

!-- calculate FF' from new p' and (old) jtor
      call update_fdf(equil_out, xaxis)
!-- resulting FF' (fdf) already normalized to 1 on axis

      diff_fdf = 0._r8
      do j = 1, npts
        diff_fdf = diff_fdf + abs(fdf(j) - fdf_old(j)) / (fdf_eps &
         + abs(fdf(j) + fdf_old(j)) / 2._r8)
      end do

      if (verbosity > 1) &
       write(iu6, *) "residual for FF' = ", diff_fdf

!-- exit current loop if fdf converged
      if (diff_fdf < errcur) converged_current = .true.

!-- re-calculate gam_spline and gam_int after change of fdf (already
!   normalized)
      dgam = fdf
!-- spline profiles
      call spline(npts, psivec, dgam, 0._r8, 0._r8, 2, gam_spline)
!-- integrate profiles
      gam_int = integrate(gam_spline, psivec, npts, .true.) &
       * sign(1._r8, cpsurfin)

!-- adapt B for i > 1
      b_ref = b
      if (verbosity > 2) &
       write(iu6, *) 'old B = ', b_ref
      if (hbt) then
        b = 2._r8 * eps * equil_out%profiles_1d%pprime(1) &
         / (equil_out%profiles_1d%pprime(1) &
         + equil_out%profiles_1d%ffprime(1) / (mu0 * rvac**2))
      else
        b = mu0 * rvac**2 * equil_out%profiles_1d%pprime(1) &
         / equil_out%profiles_1d%ffprime(1)
      end if
      if (verbosity > 2) &
       write(iu6, *) 'new B = ', b

!-- re-normalize pressure profiles (gamma profiles normalized in update_fdf)
      pscale = dpres(1)
      if (verbosity > 3) &
       write(iu6, *) 'new pscale = ', pscale
      dpres = dpres / pscale
      p_int = p_int / pscale
      dp_spline%sp1 = dp_spline%sp1 / pscale
      dp_spline%sp2 = dp_spline%sp2 / pscale
      dp_spline%sp3 = dp_spline%sp3 / pscale
      dp_spline%sp4 = dp_spline%sp4 / pscale

    else

      converged_current = .true.

    end if

!-- re-install old mesh for inner iteration
    xx = xxold
    yy = yyold
    psi = psiold

!-- inner iteration with fixed dp and fdf
    call dp_fdf_iteration(first, a, psaxis, xaxis, yaxis, nax, rax, sax)

!-- remeshing needed for calculation of flux surface quantities
    call remesh(a, nr, np, meshno, cx, cy, xaxis, yaxis, nax, rax, sax)

!-- adapt alfa to new a
    call current_calculations_constructor(nr, np)
    call adapt_alfa(a, xaxis, cx, cy, npts)
    call current_calculations_destructor

!-- fill equilibrium with new converged solution
    if (trim(input_type) == "p and j_tor") then
      call fill_equilibrium(a, xaxis, yaxis, current_averaging, equil_out)
      if (trim(radial_coordinate) == "rho_tor") then
!-- calculate q and rho_tor
        call q_calculation(xaxis, cx, cy, a, equil_out%profiles_1d%q)

        call allocate_spline_coefficients(coeff, nr)
        call spline(nr, equil_out%profiles_1d%psi, equil_out%profiles_1d%q, &
         0._r8, 0._r8, 3, coeff)
        equil_out%profiles_1d%phi = integrate(coeff, &
         equil_out%profiles_1d%psi, nr, .false.)
        call deallocate_spline_coefficients(coeff)

        equil_out%profiles_1d%rho_tor = sqrt(equil_out%profiles_1d%phi / pi &
         / equil_out%global_param%toroid_field%b0)

      end if
    end if

    if (verbosity > 1) &
     write(iu6, *) 'current iteration ', i, ' completed'

!-- calculate W_MHD
!    rminor = equil_out%eqgeometry%a_minor
!    mhd_energy = w_mhd(equil_out, rminor, eps)
!    volume = plasma_volume(equil_out, rminor, eps)
!    if (verbosity >= 2) then
!      write(iu6, *) 'W_MHD = ', mhd_energy, ' J'
!      write(iu6, *) 'plasma volume = ', volume, ' m^3'
!    end if

    if (verbosity > 1) &
     write(iu6, *) "FF' converged? :", converged_current

    if (converged_current) exit

  end do
!-- end of loop over current

  return
end subroutine p_j_iteration