psi_error.f90

Go to the documentation of this file.

subroutine psi_error
!-----------------------------------------------------------------------
! subroutine to evaluate the error in psi and grad(psi) on the nodes
! input is a restored psi vector (not called by Helena)
!-----------------------------------------------------------------------

  use itm_types
  use mod_mesh
  use mod_helena_io, only : out_he
  use mod_interpolation
  use mod_output

  implicit none

  real(r8), dimension(10, 10) :: abspsi, abspsx, abspsy, abpsxy
  real(r8) :: errmps, errmpx, errmpy, errpxy, errps, errpsx, errpsy
  real(r8) :: s, r
  real(r8) :: x, xr, xs, xrs, xrr, xss
  real(r8) :: y, yr, ys, yrs, yrr, yss
  real(r8) :: eqpsi, psir, psis, psrs, psrr, psss
  real(r8) :: solpsi, spsix, spsiy, spsixy
  real(r8) :: ejac
  real(r8) :: psix, psiy
  real(r8) :: rmax, smax
  real(r8) :: rx, ry, sx, sy
  real(r8) :: er, es
  real(r8) :: rxy, sxy, psixy
  integer(itm_i4) :: i, j, k, l
  integer(itm_i4) :: node, n1, n2, n3, n4
  integer(itm_i4) :: imax, jmax

  do k = 1, 10
    do l = 1, 10
      abspsi(k, l) = 0._r8
      abspsx(k, l) = 0._r8
      abspsy(k, l) = 0._r8
      abpsxy(k, l) = 0._r8
    end do
  end do
  errmps = 0._r8
  errmpx = 0._r8
  errmpy = 0._r8
  errpxy = 0._r8
  errps = 0._r8
  errpsx = 0._r8
  errpsy = 0._r8
  do l = 1, 10
    s = -1._r8 + 2._r8 * (l - 1) / 10._r8
    do i = 1, nr - 1
      do j = 1, np - 1
        do k = 1, 10
          r = -1._r8 +  2._r8 * (k - 1) / 10._r8
          node = (i - 1) * np + j
          n1 = node
          n2 = node + 1
          n3 = node + 1 + np
          n4 = node + np
          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, &
           eqpsi, psir, psis, psrs, psrr, psss)
          call soloviev(x, y, solpsi, spsix, spsiy, spsixy)
          abspsi(k, l) = abspsi(k, l) + abs(solpsi - eqpsi)
          ejac = xr * ys - xs * yr
          psix = (ys * psir - yr * psis) / ejac
          psiy = (-xs * psir + xr * psis) / ejac
          abspsx(k, l) = abspsx(k, l) + abs(psix - spsix)
          abspsy(k, l) = abspsy(k, l) + abs(psiy - spsiy)
          if (abs(solpsi - eqpsi) > errmps) then
            errmps = abs(solpsi - eqpsi)
          end if
          if (abs(psix - spsix) > errmpx) then
            errmpx = abs(psix - spsix)
          end if
          if (abs(psiy - spsiy) > errmpy) then
            errmpy = abs(psiy - spsiy)
            rmax = r
            smax = s
            imax = i
            jmax = j
          end if
          ry = -xs / ejac
          rx = ys / ejac
          sy = xr / ejac
          sx = -yr / ejac
          er = xrr * ys + xr * yrs - xrs * yr - xs * yrr
          es = xrs * ys + xr * yss - xss * yr - xs * yrs
          rxy = -(er * ry + es * sy) * ys / ejac**2 + (yrs * ry &
           + yss * sy) / ejac
          sxy = (er * ry + es * sy) * yr / ejac**2 - (yrr * ry &
           + yrs * sy) / ejac
          psixy = psrr * rx * ry + psrs * (rx * sy + sx * ry) &
           + psir * rxy + psss * sx * sy + psis * sxy
          if (abs(psixy - spsixy) > errpxy) then
            errpxy = abs(psixy - spsixy)
          end if
          abpsxy(k, l) = abpsxy(k, l) + abs(psixy - spsixy)
        end do
      end do
    end do
    do k = 1, 10
      r = -1._r8 + 2._r8 * (k - 1) / 10._r8
      abspsi(k, l) = abspsi(k, l) / ((nr - 1) * (np - 1))
      abspsx(k, l) = abspsx(k, l) / ((nr - 1) * (np - 1))
      abspsy(k, l) = abspsy(k, l) / ((nr - 1) * (np - 1))
      errps = errps + abspsi(k, l)
      errpsx = errpsx + abspsx(k, l)
      errpsy = errpsy + abspsy(k, l)
    end do
  end do
  if (standard_output) then
    write(out_he, *) 'max err in psi :', errmps
    write(out_he, *) 'max err in psix:', errmpx
    write(out_he, *) 'max err in psiy:', errmpy
    write(out_he, *) ' at i,j : ', i, j
    write(out_he, *) '    r,s : ', rmax, smax
    write(out_he, *) 'max err in psixy :', errpxy
    write(out_he, *) 'average error in psi : ', errps / 100.
    write(out_he, *) 'average error in psix : ', errpsx / 100.
    write(out_he, *) 'average error in psiy : ', errpsy / 100.
  end if

  return
end subroutine psi_error