solver.f90

Go to the documentation of this file.

module solver

  use itm_types
  use mod_corners
  use mod_mesh
  use mod_helena_io, only : out_he, iu6
  use mod_output, only : standard_output

  implicit none

  real(r8), allocatable, private :: qq(:)
  real(r8), allocatable, private :: gpx(:, :, :), gpjac(:, :, :)
  integer(itm_i4), allocatable, private :: indices(:), inv_index(:)


  contains

  subroutine allocate_solver(nr, np)

    use itm_types

    implicit none

    integer(itm_i4) :: nr, np

    allocate(qq(4 * nr * np))
    allocate(gpx(4, 4, (nr - 1) * (np - 1)))
    allocate(gpjac(4, 4, (nr - 1) * (np - 1)))
    allocate(indices(nr * np))
    allocate(inv_index(nr * np))

!-- initialize
    qq = 0._r8
    gpx = 0._r8
    gpjac = 0._r8
    indices = 0
    inv_index = 0

    return
  end subroutine allocate_solver

  subroutine deallocate_solver

    implicit none

    deallocate(qq)
    deallocate(gpx, gpjac)
    deallocate(indices, inv_index)

  end subroutine deallocate_solver

  subroutine find_axis(psaxis, xaxis, yaxis, nax, rax, sax)
!-----------------------------------------------------------------------
! subroutine to localize the position of the magnetic axis ; the minimum
! of psi of all elements
!-----------------------------------------------------------------------
    
    use mod_interpolation
    use matrix_calculations
    
    implicit none
    
    real(r8), intent(out) :: psaxis
    real(r8), intent(inout) :: xaxis
    real(r8), intent(inout) :: yaxis
    real(r8), intent(inout) :: rax, sax
    integer(itm_i4), intent(out) :: nax
    
    real(r8), dimension(2) :: x
    real(r8) :: psmin
    real(r8) :: x_err, f_err
    real(r8), parameter :: tol_x=1.e-8_r8, tol_f=1.e-8_r8
    real(r8), parameter :: tol_rs=1.0001_r8
    real(r8) :: r, s, zpsi
    real(r8) :: psi_00, psi_00_min
    integer(itm_i4) :: i, n_start
    integer(itm_i4) :: naxis
    integer(itm_i4) :: n1, n2, n3, n4
    integer(itm_i4) :: nelm
    integer(itm_i4) :: neq2, ntrial, ifail
    integer(itm_i4) :: naxis_00
    logical :: minimum_found

    minimum_found = .false.

    psmin = 1.e20_r8
    psi_00_min = psmin
    naxis_00 = -1

    if (ias == 1) then
!----------------------------------- asymmetric part -----------      
      nelm = (nr - 1) * (np - 1)
      neq2 = 2
      ntrial = 150
      do i = nelm, nelm / 2, -1
        naxis = i
        ifail = 1
        x(1) = 0._r8
        x(2) = 0._r8
        call set_node_number(naxis, n1, n2, n3, n4)
        call interpolation(0, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) + 1), &
         psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1),  x(1), x(2), psi_00)
        if (psi_00 .lt. psi_00_min) then
          psi_00_min = psi_00
          naxis_00 = naxis
        end if
        call mnewtax(ntrial, naxis, x, neq2, tol_x, tol_f, x_err, f_err)
        if (x_err <= tol_x .or. f_err <= tol_f) then
          ifail = 0
        else
!         if (standard_output) &
!          write(out_he, *) ' accuracy not reached : ', x_err, f_err, naxis
        end if 
        r = x(1)
        s = x(2)
        if (ifail == 0 .and. abs(r) <= tol_rs .and. abs(s) <= tol_rs) then
          minimum_found = .true.
          call set_node_number(naxis, n1, n2, n3, n4)
          call interpolation(0, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) &
           + 1), psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, zpsi)
          if (zpsi < psmin) then
            psmin = zpsi
            nax = naxis
            rax = r
            sax = s
          end if
        end if
      end do
      if (naxis_00 .gt. 0) then
!      write(*,*) 'naxis_00, psi_00_min ', naxis_00, psi_00_min 
      end if
      if (.not. minimum_found) then
        write(iu6, *) 'fatal error: no axis found'
        stop
      end if
      call set_node_number(nax, n1, n2, n3, n4)
      call interpolation(0, xx(1, n1), xx(1, n2), xx(1, n3), xx(1, n4), &
       rax, sax, xaxis)
      call interpolation(0, yy(1, n1), yy(1, n2), yy(1, n3), yy(1, n4), &
       rax, sax, yaxis)
      psaxis = psmin
    else
!**************************************** symmetric part ***************
      n_start = 1
      call axis_symm(psaxis, xaxis, yaxis, nax, rax, sax, n_start)
      n_start = np - 1
      call axis_symm(psaxis, xaxis, yaxis, nax, rax, sax, n_start)
      yaxis = 0._r8
    end if
    
    return
  end subroutine find_axis   

  subroutine axis_symm(psaxis, xaxis, yaxis, nax, rax, sax, n_start)

    use math_functions

    implicit none

    real(r8), intent(inout) :: psaxis
    real(r8), intent(inout) :: xaxis
    real(r8), intent(inout) :: yaxis
    real(r8), intent(inout) :: rax, sax
    integer(itm_i4), intent(inout) :: nax
    integer(itm_i4), intent(in) :: n_start

    real(r8) :: r, dummy, psmin
    real(r8), parameter :: tol = 1._r8 + 1.e-8_r8 
    integer(itm_i4) :: n, n1, n2, ifail

    psmin = 1.e20_r8

    do n = n_start, (nr - 1) * (np - 1), np - 1
      if (n_start == 1) then
        n1 = nodeno(n, 1)
        n2 = nodeno(n, 4)
      else 
        n1 = nodeno(n, 2)
        n2 = nodeno(n, 3)
      end if
      if (psi(4 * (n1 - 1) + 2) * psi(4 * (n2 - 1) + 2) <= 0.) then
!------------------------------------- quad. eq for r value at minimum -
        call quad_equation(psi(4 * (n1 - 1) + 1), psi(4 * (n1 - 1) + 2), &
         psi(4 * (n2 - 1) + 1), psi(4 * (n2 - 1) + 2), tol, r, ifail) 
!-------- the sign of r changes for elements on the left  (see remesh) -
        call cubic_interp_1d(xx(1, n1), xx(2, n1), xx(1, n2), &
         xx(2, n2), r, xaxis, dummy)
        call cubic_interp_1d(psi(4 * (n1 - 1) + 1), psi(4 * (n1 - 1) + 2), &
         psi(4 * (n2 - 1) + 1), psi(4 * (n2 - 1) + 2),r, psaxis, dummy)
        if (psaxis < psmin) then
          psmin = psaxis
          nax = n
          rax = r
          if (n_start == 1) then
            sax = -1._r8
          else 
            sax = 1._r8
          end if
          if (standard_output) &
           write(out_he, *) 'local minimum at xaxis = ', xaxis
        end if
      end if
    end do
  
    return
  end subroutine axis_symm

  subroutine build_matrix(a, first)
!----------------------------------------------------------------------
! subroutine to calculate the matrix kk and the right hand side array qq
! no boundary conditions are yet used.
! number of rows and columns : 4 * nr * np
! nr : number of radial nodes
! np : number of poloidal nodes
! a  : the total amplitude of the rhs (hbt definition)
! b  : measure of the total pressure   (hbt definition)
! eps : the inverse aspect ratio
! hbt = .true. : hbt definition of gamma, .false. : ff' as gamma.
! gpx, gpjac, indices used as static variables !
!----------------------------------------------------------------------

    use itm_types
    use mod_corners
    use mod_dat
    use mod_gauss
    use mod_mesh
    use mod_solver
    use mod_interpolation
    use matrix_calculations
    use phys_functions

    implicit none

    real(r8), intent(in) :: a
    logical, intent(in) :: first

    real(r8) :: r, wr, s, ws, wrs
    real(r8) :: x, xr, xs, yr, ys, ps
    real(r8) :: xjac
    real(r8) :: arhs
    real(r8) :: sumq, sumk
    real(r8) :: vx, vy
    real(r8) :: psix, psiy
    integer(itm_i4) :: i, j, k, l, n, ngr, ngs
    integer(itm_i4) :: nelm
    integer(itm_i4) :: n1, n2, n3, n4
    integer(itm_i4) :: nrow, ncol, noff, nend
    logical, parameter :: invert_index = .false.

    nelm = (np - 1) * (nr - 1)

!-- initialize qq to zero
    qq = 0._r8

    if (first) then
      call set_index(indices, invert_index)
!------------------------------------- nelm elements ------------------
      do i = 1, 4 * nr * np
        do j = 1, kklda
           kkbig(j, i) = 0._r8
        end do
      end do
      do n = 1, nelm
        call set_node_number(n, n1, n2, n3, n4)
!------------------------------------- 4 point gaussian int. in r -----
        do ngr = 1, 4
          r = xgauss(ngr)
          wr = wgauss(ngr)
!------------------------------------- 4 point gaussian int. in s -----
          do ngs = 1, 4
            s = xgauss(ngs)
            ws = wgauss(ngs)
            wrs = wr * ws
            call interpolation(3, xx(1, n1), xx(1, n2), xx(1, n3), &
             xx(1, n4), r, s, x, xr, xs, yn1 = yy(1, n1), yn2 = yy(1, n2), & 
             yn3 = yy(1, n3), yn4 = yy(1, n4), pn1 = psi(4 * (n1- 1) + 1), &
             pn2 = psi(4 * (n2 - 1) + 1), pn3 = psi(4 * (n3 - 1) + 1), &
             pn4 = psi(4 * (n4 - 1) + 1), yr = yr, ys = ys, ps = ps)

            xjac = xr * ys - xs * yr
            gpx(ngr, ngs, n) = x
            gpjac(ngr, ngs, n) = xjac
!----------------------------------------- calculate right hand side ---
            arhs = rh_side(a, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, 0._r8)
!------------------------------------- 4 nodes per element ------------
            do i = 1, 4
!------------------------------------- 4 functions v per node ---------
              do j = 1, 4
                nrow = 4 * (indices(nodeno(n, i)) - 1) + j
                sumq = -arhs * h(j, i, ngr, ngs) * xjac
                vx = ys * hr(j, i, ngr, ngs) - yr * hs(j, i, ngr, ngs)
                vy = -xs * hr(j, i, ngr, ngs) + xr * hs(j, i, ngr, ngs)
                qq(nrow) = qq(nrow) - wrs * sumq
!------------------------------------- 4 nodes of function psi --------
                do k = 1, 4
!------------------------------------- 4 functions h in psi -----------
                  do l = 1, 4
                    ncol = 4 * (indices(nodeno(n, k)) - 1) + l
                    noff = nrow - ncol
                    if (noff >= 0) then
                      psix = ys * hr(l, k, ngr, ngs) - yr * hs(l, k, ngr, ngs)
                      psiy = -xs * hr(l, k, ngr, ngs) + xr * hs(l, k, ngr, ngs)
                      sumk = -1. / (1. + eps * x) * (psix * vx + psiy * vy) &
                       / xjac
                      kkbig(noff + 1, ncol) = kkbig(noff + 1, ncol) + wrs &
                       * sumk
                    end if
                  end do
                end do
              end do
            end do
          end do
        end do
      end do
!------------------------------------------- remove empty columns (bnd.)
      nend = np - 1
      if (ias == 0) nend = np
      do j = 1, nend
        kkbig(1, 4 * j - 3) = 1.e20_r8
        kkbig(1, 4 * j - 1) = 1.e20_r8
      end do
!------------------------------------------- if matrix exists then
    else
      do n = 1, nelm
        call set_node_number(n, n1, n2, n3, n4)
!------------------------------------- 4 point gaussian int. in r -----
        do ngr = 1, 4
          r = xgauss(ngr)
          wr = wgauss(ngr)
!------------------------------------- 4 point gaussian int. in s -----
          do ngs = 1, 4
            s = xgauss(ngs)
            ws = wgauss(ngs)
            x = gpx(ngr, ngs, n)
            xjac = gpjac(ngr, ngs, n)
            call interpolation(0, psi(4 * (n1 - 1) + 1), psi(4 * (n2 - 1) &
             + 1), psi(4 * (n3 - 1) + 1), psi(4 * (n4 - 1) + 1), r, s, ps)
!----------------------------------------- calculate right hand side ---
            arhs = rh_side(a, x, 0._r8, 0._r8, 0._r8, 0._r8, ps, 0._r8, 0._r8)
!------------------------------------- 4 nodes per element ------------
            do i = 1, 4
!------------------------------------- 4 functions v per node ---------
              do j = 1, 4
                nrow = 4 * (indices(nodeno(n, i)) - 1) + j
                sumq = -arhs * h(j, i, ngr, ngs) * xjac
                qq(nrow) = qq(nrow) - wr * ws * sumq
              end do
            end do
          end do
        end do
      end do
    end if

    return
  end subroutine build_matrix

  subroutine gaussian_elimination(first)
!-----------------------------------------------------------------------
! subroutine to solve the system of equations using gaussian elimination
!-----------------------------------------------------------------------

    use itm_types
    use mod_parameter
    use mod_mesh
    use mod_solver
    use matrix_calculations, only : set_index

    implicit none

    logical, intent(in) :: first
  
    integer(itm_i4) :: nd, nrow, nvar
    integer(itm_i4) :: i, j, k
    integer(itm_i4) :: info
    integer(itm_i4) :: ik, ik2, ikn
    integer(itm_i4) :: nbase, nbase2
    logical, parameter :: invert = .true.

    if (ias == 0) then
      nd = 4 * nr * np
    else
      nd = 4 * nr * (np - 1)
    end if
    nrow = 4 * nr * np
    nvar = 4 * np + 7

    if (first) then
!------------------------------inverse of indices in build_matrix to restore
      if (ias == 1) then
        inv_index = 0
        call set_index(inv_index, invert)
      end if

!---------------------------------------------------- essl version
!      call dpbf(kkbig, kklda, nd, 4 * np + 8)
!---------------------------------------------------- lapack version      
      call dpbtrf('l', nd, nvar, kkbig, kklda, info)        
    end if

    psi(1 : 4 * nr * np) = qq 

!-------------------------------------------------- essl version
!    call dpbs(kkbig, kklda, nd, 4 * np + 8, psi)
!-------------------------------------------------- lapack version
    call dpbtrs('l', nd, nvar, 1, kkbig, kklda, psi, nrow, info)

!-------------- restore to simple clockwise numbering
    if (ias == 1) then
      qq = psi(1 : 4 * nr * np)
      do i = 1, nr * (np - 1)
        do k = 1, 4
          ik = 4 * (i - 1) + k
          ikn = 4 * (inv_index(i) - 1) + k
          psi(ikn) = qq(ik)
        end do  
      end do
      do i = 1, nr
        do k = 1, 4
          ik = 4 * (i - 1) * np + k
          ik2 = 4 * (i - 1) * np + 4 * (np - 1) + k
          psi(ik2) = psi(ik)
        end do
      end do
    end if
!-------------------------------- fill in boundary conditions      
    do i = 1, nr * np
      psi(4 * (i - 1) + 1) = psi(4 * (i - 1) + 1) + 1._r8
    end do
    do j = 1, np
      psi(4 * j - 3) = 1._r8
      psi(4 * j - 1) = 0._r8
    end do
    if (ias == 1) then
      do i = 1, nr
        do k = 1, 4
          nbase = 4 * (i - 1) * np + k
          nbase2 = nbase + 4 * (np - 1)
          psi(nbase2) = psi(nbase)
        end do
      end do
    end if

    return
  end subroutine gaussian_elimination

  subroutine normalize_psi(psaxis)
!-----------------------------------------------------------------------
! subroutine to normalize psi to one on the boundary and zero on axis
!-----------------------------------------------------------------------
  
    use itm_types
    use mod_mesh
  
    implicit none
  
    real(r8), intent(in) :: psaxis
  
    integer(itm_i4) :: i, l
  
    do i = 1, nr * np
      psi(4 * (i - 1) + 1) = 1._r8 - (1._r8 - psi(4 * (i - 1) + 1)) &
       / (1._r8 - psaxis)
      do l = 2, 4
        psi(4 * (i - 1) + l) = psi(4 * (i - 1) + l) / (1._r8 - psaxis)
      end do
    end do
  
    return
  end subroutine normalize_psi

end module solver