check_solution.f90

Go to the documentation of this file.

subroutine check_solution(equilibrium_out)
!-----------------------------------------------------------------------
!    - checks the error in the grad-shafranov equation using the 
!      straight field line system (not called by Helena)
!-----------------------------------------------------------------------

  use itm_types
  use mod_map
  use mod_mesh, only : nr
  use mod_parameter
  use mod_profiles
  use mod_spline
  use mod_helena_io, only : out_he
  use mod_output
  use helena_spline

  use euitm_schemas

  implicit none

  type(type_equilibrium), intent(in) :: equilibrium_out

  type(spline_coefficients) :: gp_spline
  integer(itm_i4) :: i, j, im, jm
  real(r8) :: maxerr, err
  real(r8) :: dgp0, dgp1, dgp
  real(r8) :: sps2
  real(r8) :: dp_axis, dp_boundary
  real(r8) :: drbphi_axis, drbphi_boundary
  real(r8), dimension(npsi - 1) :: c1
  real(r8), dimension(3) :: dummy
  real(r8), dimension(nchi) :: dgt
  real(r8), dimension(npsi) :: gp
  real(r8), dimension(npsi - 1, nchi) :: error, error2
!  real :: ff(nchi), fp(nchi), zc(10)

!-- derivatives at plasma boundaries (all zero except for dp_boundary)
  dp_axis = 0._r8
  drbphi_axis = 0._r8
  drbphi_boundary = 0._r8
  dp_boundary = equilibrium_out%profiles_1d%pprime(nr) &
   * 2 * equilibrium_out%profiles_1d%psi(nr) * mu0 &
   / equilibrium_out%global_param%mag_axis%bphi**2

!TODO: is this really supposed to change npsi?
  npsi = nr

  do i = 1, nchi
    c1(:) = g12_hel(2 : npsi - 1, i)
    call spline(npsi - 1, cs(2), c1, 0._r8, 0._r8, 2, q_spline)
    g12_hel(1, i) = spwert(npsi - 1, 0._r8, q_spline, cs(2), dummy, 0)
  end do

  call spline(npsi, cs, p0, dp_axis, dp_boundary, 1, p_spline)
  call spline(npsi, cs, rbphi, drbphi_axis, drbphi_boundary, 1, rbphi_spline)

  maxerr = -1.e20_r8

  call allocate_spline_coefficients(gp_spline, npsi)

  do j = 1, nchi
    do i = 1, npsi
      gp(i) = g11_hel(i, j) * qs(i) / rbphi(i)
    end do
    dgp0 = (gp(2) - gp(1)) / (cs(2) - cs(1))
    dgp1 = (gp(npsi) - gp(npsi - 1)) / (cs(npsi) - cs(npsi - 1))
    call spline(npsi, cs, gp, dgp0, dgp1, 2, gp_spline)
!    call spline(npsi, cs, gp, 0._r8, 0._r8, 2, gp_spline)

    do i = 2, npsi - 1
      sps2 = 2._r8 * cs(i) * cpsurf
      call theta_derivative(g12_hel(i, 1), dgt, nchi)
      dgp = (g11_hel(i + 1, j) * qs(i + 1) / rbphi(i + 1) &
       - g11_hel(i - 1, j) * qs(i - 1) / rbphi(i - 1)) &
       / (cs(i + 1) - cs(i - 1))
      err = abs(rbphi(i) / qs(i) * dgp + dgt(j) * sps2 &
       + (rbphi(i) * rbphi_spline%sp2(i) + g33_hel(i, j) * p_spline%sp2(i)))
      if (err > maxerr) then
        maxerr = err
        im = i
        jm = j
      end if
!     if (standard_output) &
!      write(out_he, 21) i, j, cs(i), chi(j), err, log(err)
      error(i - 1, j) = log(err)
      error2(i - 1, j) = err
    end do
  end do
  if (standard_output) then
    write(out_he, *)
    write(out_he, 23) maxerr, im, jm
    write(out_he, *)
  end if

! if (standard_output) then
!   write(out_he, *) 
!   write(out_he, *) ' error in gs : harmonics'
!   write(out_he, *)
!   do j = 1, nchi
!     fp(j) = float(j - 1)
!   end do 
!   do i = 2, npsi
!     do j = 1, nchi
!       ff(j) = error2(i - 1, j)
!     end do
!     call rft2(ff, nchi, 1)
!     do m = 1, nchi / 2
!       ff(m) = 2 * ff(m) / real(nchi)
!       write(out_he, 22) i, m, cs(i), ff(m)
!     end do
!   end do
! end if

  call deallocate_spline_coefficients(gp_spline)

  return

! 21 format(2i3, 2f8.3, 3e12.4)
! 22 format(2i3, f8.3, 3e12.4)
 23 format(' max error in gs after mapping : ', e12.4, 2i5)

end subroutine check_solution