dp_fdf_iteration.f90

Go to the documentation of this file.

subroutine dp_fdf_iteration(first, a, psaxis, xaxis, yaxis, nax, rax, sax)
!-----------------------------------------------------------------------------
! subroutine to solve the Grad-Shafranov equation for p' and FF' input
!-----------------------------------------------------------------------------

  use itm_types
  use mod_dat, only : errit, amix, niter, verbosity
  use mod_helena_io, only : iu6, out_he
  use mod_mesh, only : psi, psiold, nr, np
  use mod_output
  use mod_solver
  use solver
  use mod_interpolation
  use matrix_calculations

  implicit none

  logical, intent(inout) :: first

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

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

  real(r8) :: residual
  integer(itm_i4) :: i, j

!-- start of inner iteration (with fixed dp and fdf)
  do i = 1, niter

    psiold = psi

    !call system('date +%s%N')

!-- build matrix, no boundary conditions
    call build_matrix(a, first)

    !call system('date +%s%N')

!-- solve set of equations
    call gaussian_elimination(first)

    !call system('date +%s%N')

!-- matrix built and decomposed
    if (first) first = .false.

!-- find magnetic axis and psi on axis
    call find_axis(psaxis, xaxis, yaxis, nax, rax, sax)

    !call system('date +%s%N')

    a = a / (1._r8 - psaxis)

    if (verbosity > 1) &
     write(iu6, *) 'iteration ', i
    if (verbosity > 2) then
      write(iu6, 1) xaxis, yaxis
      write(iu6, *) 'new A = ', a, ' psaxis = ', psaxis
    end if

 1  format('magnetic axis : x   = ', e14.6, '  y   = ', e14.6)

!-- normalize flux to zero on axis
    call normalize_psi(psaxis)
    if (amix /= 0._r8) then
      psi = (1._r8 - amix) * psi + amix * psiold
      call find_axis(psaxis, xaxis, yaxis, nax, rax, sax)
      a = a / (1._r8 - psaxis)
      call normalize_psi(psaxis)
    end if

!-- calculate residual
    residual = 0._r8
    do j = 1, 4 * nr * np
      residual = residual + abs(psi(j) - psiold(j))
    end do
    residual = residual / dble(nr * np)
    write(iu6, *) 'helena: residual = ', residual
    if (residual < errit) exit

    if (verbosity > 0) then
      if (i == niter) then
         write(iu6, *) 'WARNING: inner loop not fully converged!'
         read(*,*) 
         stop 'not converging in fp_fdf_iteration'
      end if
    end if

    !call system('date +%s%N')

    write(*, *) 'iteration ', i, ' end'

  end do
!-- end of inner iteration

  if (standard_output) write(out_he, 2) i, residual

 2 format('    residual after ', i3, 'iterations = ', es10.3)

  return
end subroutine dp_fdf_iteration