leq.f90 Source File

#include "unused_dummy.inc"
!> Analytic local equilibria via generalised forms of the Miller model
!>
!> Justin Ball's MIT Masters thesis and/or ["GS2 analytic geometry
!> specification"](https://bitbucket.org/gyrokinetics/wikifiles/raw/master/JRB/GS2GeoDoc.pdf)
!> provide a good overview of how this module works.
!>
!> There are four equilibrium models in this module:
!>
!> 1. generalised Miller
!> 2. global MHD
!> 3. generalised elongation
!> 4. Fourier
!>
!> These can be chosen via the `geoType` input parameter (`0-4`,
!> respectively)
!>
!> All of these models allow you to specify the amplitudes of two
!> shaping modes, along with their phases and radial derivatives. The
!> generalised elongation and Fourier models also allow you to specify
!> the exact mode numbers.
module leq
  use geo_utils, only: abstract_geo_type, flux_surface_type, geo_input_type, &
       geo_type_miller, geo_type_global, geo_type_generalized_elongation, geo_type_fourier_series, &
       geo_type_miller_extended_harmonic

  implicit none

  private

  public :: leq_type

  type, extends(abstract_geo_type) :: leq_type
     !> Used to translate definition of the poloidal angle such that the location
     !> of the maximium magnetic field is at pi (see section 3.2.3 of Ball MIT
     !> Masters thesis or section 2.4 of "GS2 analytic geometry specification")
     real :: thetaShift
     type(flux_surface_type) :: surf
   contains
     procedure :: finish_setup, finalise, read_and_set
     procedure :: rpos, zpos, rcenter, diameter, btori, dbtori, qfun
     procedure :: pfun, dpfun, betafun, psi, rfun, eqitem
     procedure :: geo_Miller, geo_global, geo_generalizedElongation, geo_FourierSeries
     procedure :: geo_MillerExtendedHarmonic
  end type leq_type

contains

  !> Set module-level parameters and initialise the leq module
  subroutine read_and_set(self, inputs)
    use mp, only: mp_abort
    use constants, only: twopi, pi
    implicit none
    class(leq_type), intent(in out) :: self
    type(geo_input_type), intent(in) :: inputs
    real, dimension(3) :: dr
    integer :: i, j
    real :: r, t
    self%type_name = 'Local-Analytic'
    if (inputs%surf%geoType > 4 .or. inputs%surf%geoType < 0) &
       call mp_abort("ERROR: invalid analytic geometry specification", .true.)

    self%surf = inputs%surf ;  self%thetaShift = inputs%thShift
    self%beta_0 = 1.
    self%nr = 3 ; self%nt = self%surf%nt + 1
    call self%alloc_arrays()
    dr = [-self%surf%dr, 0., self%surf%dr]
    ! Why do we set pressure to -dr here?
    self%pressure = -dr
    self%eqpsi = 1 + dr

    do j = 1, self%nt
       t = (j - 1) * twopi / real(self%nt - 1) - pi
       self%eqth(:, j) = t
       self%eqpsi_2d(:, j) = self%eqpsi
       do i = 1, self%nr
          r = self%surf%r + dr(i)
          self%R_psi(i, j) = self%Rpos(r, t)
          self%Z_psi(i, j) = self%Zpos(r, t)
       end do
    end do

    self%psi_0 = 0 ; self%psi_a = 1.0
    self%B_T = -1.0 ; self%R_mag = self%surf%Rmaj ; self%Z_mag = 0.0
    self%aminor = 1.0
    self%has_full_theta_range = .true.
  end subroutine read_and_set

  !> Finalise leq module
  pure subroutine finalise(self)
    implicit none
    class(leq_type), intent(in out) :: self
    self%initialised = .false.
    call self%dealloc_arrays
  end subroutine finalise

  !> Finish the leq setup
  pure subroutine finish_setup(self)
    class(leq_type), intent(in out) :: self
    call self%calculate_gradients
  end subroutine finish_setup

  !> Calculates `fstar` which is `f` interpolated at the location `theta_in`
  pure real function eqitem(self, r_in, theta_in, f, char) result(fstar)
    use geo_utils, only: mod2pi, find_lower_bound, bilinear_interpolate
    implicit none
    class(leq_type), intent(in) :: self
    real, intent(in) :: r_in !Ignored
    !> Desired theta location
    real, intent(in) :: theta_in
    !> Array to interpolate
    real, dimension(:, :), intent(in) :: f
    !> FIXME: Unused
    character(1), intent(in) :: char
    !> Found r and theta indices
    integer :: jstar, istar
    real :: thet
    !> Internal theta grid
    real, dimension(self%nt) :: mtheta
    real, dimension(self%nr) :: dim_1
    UNUSED_DUMMY(r_in); UNUSED_DUMMY(char)
    dim_1 = [1, 2, 3] ; mtheta = self%eqth(1, :)
    thet = mod2pi(theta_in)
    istar = 2 !find_lower_bound(2.5, dim_1, 0)  !Make 2 < r < 3 here so istar == 2
    jstar = find_lower_bound(thet, mtheta, self%nt - 1)
    fstar = bilinear_interpolate(f, dim_1, mtheta, istar, jstar, dim_1(istar), thet)
  end function eqitem

  !> Return the major radius at \((r, \theta)\)
  pure real function Rpos (self, r, theta)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real :: Zpos

    ! Choose one of four different analytic geometry specifications
    ! (see section 2.1 of "GS2 analytic geometry specification") ! JRB
    select case (self%surf%geoType)
       case (geo_type_miller)
          call self%geo_Miller(r, theta, Rpos, Zpos)
       case (geo_type_global)
          call self%geo_global(r, theta, Rpos, Zpos)
       case (geo_type_generalized_elongation)
          call self%geo_generalizedElongation(r, theta, Rpos, Zpos)
       case (geo_type_fourier_series)
          call self%geo_FourierSeries(r, theta, Rpos, Zpos)
       case (geo_type_miller_extended_harmonic)
          call self%geo_MillerExtendedHarmonic(r, theta, Rpos, Zpos)
    end select
  end function Rpos

  !> Return the axial coordinate at \((r, \theta)\)
  pure real function Zpos (self, r, theta)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real :: Rpos

    ! Choose one of four different analytic geometry specifications
    ! (see section 2.1 of "GS2 analytic geometry specification") ! JRB
    select case (self%surf%geoType)
       case (geo_type_miller)
          call self%geo_Miller(r, theta, Rpos, Zpos)
       case (geo_type_global)
          call self%geo_global(r, theta, Rpos, Zpos)
       case (geo_type_generalized_elongation)
          call self%geo_generalizedElongation(r, theta, Rpos, Zpos)
       case (geo_type_fourier_series)
          call self%geo_FourierSeries(r, theta, Rpos, Zpos)
       case (geo_type_miller_extended_harmonic)
          call self%geo_MillerExtendedHarmonic(r, theta, Rpos, Zpos)
    end select
  end function Zpos

  !> Calculate the cylindrical \((R, Z)$$ at $$(r, \theta)\) for the
  !> generalised Miller model
  !>
  !> See section 3.2.1 of Ball MIT Masters thesis or section 2.1.1 of
  !> "GS2 analytic geometry specification"
  pure subroutine geo_Miller(self, r, theta, Rpos, Zpos)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent(in) :: r, theta
    real, intent(out) :: Rpos, Zpos
    real :: dr, thAdj
    real :: Rcirc, Rcircp, Relong, Relongp, RelongTilt, RelongTiltp, Rtri, Rtrip, RtriTilt, RtriTiltp, Rfinal, Rfinalp
    real :: Zcirc, Zcircp, Zelong, Zelongp, ZelongTilt, ZelongTiltp, Ztri, Ztrip, ZtriTilt, ZtriTiltp, Zfinal, Zfinalp

    dr = r - self%surf%r

    thAdj = theta + self%thetaShift

    Rcirc = self%surf%r*cos(self%surf%thn - self%surf%thm - thAdj)
    Zcirc = -self%surf%r*sin(self%surf%thn - self%surf%thm - thAdj)
    Rcircp = cos(self%surf%thn - self%surf%thm - thAdj)
    Zcircp = -sin(self%surf%thn - self%surf%thm - thAdj)

    Relong = Rcirc
    Zelong = -self%surf%r*(self%surf%delm - 1)*sin(self%surf%thn - self%surf%thm - thAdj) + Zcirc
    Relongp = Rcircp
    Zelongp = -((self%surf%delm - 1) + self%surf%r*self%surf%delmp)*sin(self%surf%thn - self%surf%thm - thAdj) + Zcircp

    RelongTilt = cos(self%surf%thn - self%surf%thm)*Relong - sin(self%surf%thn - self%surf%thm)*Zelong
    ZelongTilt = sin(self%surf%thn - self%surf%thm)*Relong + cos(self%surf%thn - self%surf%thm)*Zelong
    RelongTiltp = cos(self%surf%thn - self%surf%thm)*Relongp - sin(self%surf%thn - self%surf%thm)*Zelongp
    ZelongTiltp = sin(self%surf%thn - self%surf%thm)*Relongp + cos(self%surf%thn - self%surf%thm)*Zelongp

    Rtri = -self%surf%r*cos(thAdj) + self%surf%r*cos(sin(thAdj)*self%surf%deln + thAdj) + RelongTilt
    Ztri = ZelongTilt
    Rtrip = -cos(thAdj) + cos(sin(thAdj)*self%surf%deln + thAdj) &
            - self%surf%r*sin(sin(thAdj)*self%surf%deln + thAdj)*sin(thAdj)*self%surf%delnp + RelongTiltp
    Ztrip = ZelongTiltp

    RtriTilt = cos(self%surf%thn)*Rtri + sin(self%surf%thn)*Ztri
    ZtriTilt = -sin(self%surf%thn)*Rtri + cos(self%surf%thn)*Ztri
    RtriTiltp = cos(self%surf%thn)*Rtrip + sin(self%surf%thn)*Ztrip
    ZtriTiltp = -sin(self%surf%thn)*Rtrip + cos(self%surf%thn)*Ztrip

    Rfinal = RtriTilt + self%surf%Rmaj
    Zfinal = ZtriTilt
    Rfinalp = RtriTiltp + self%surf%sHorz
    Zfinalp = ZtriTiltp + self%surf%sVert

    Rpos = Rfinal + dr*Rfinalp
    Zpos = Zfinal + dr*Zfinalp

  end subroutine geo_Miller

  !> Calculate the cylindrical \((R, Z)\) at \((r, \theta)\) for the
  !> global MHD model
  !>
  !> See section 2.1.2 of "GS2 analytic geometry specification"
  pure subroutine geo_global(self, r, theta, Rpos, Zpos)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real, intent (out) :: Rpos, Zpos
    real :: rho0, thAdj, rCyl, Rc, Zc
    real :: Cm, Cn, psiN, xAng, yAng, rAng

    rho0 = r/self%surf%aSurf

    thAdj = theta + self%thetaShift

    Cm=(self%surf%delm**2-1)/(self%surf%delm**2+1)
    Cn=(self%surf%deln**2-1)/(self%surf%deln**3+1)

    psiN=rho0**2+Cm*rho0**2+Cn*rho0**3

    if (abs(Cn*cos(3*(thAdj+self%surf%thn)))<10e-6) then
      rCyl=self%surf%aSurf*sqrt(psiN/(1+Cm*cos(2*(thAdj+self%surf%thm))))
    else
      xAng=2*(1+Cm*cos(2*(thAdj+self%surf%thm)))**3-27*Cn**2*psiN*cos(3*(thAdj+self%surf%thn))**2
      yAng=3*Cn*sqrt(3*psiN)*cos(3*(thAdj+self%surf%thn))*sqrt(xAng+2*(1+Cm*cos(2*(thAdj+self%surf%thm)))**3)

      rAng=(1/3.0)*atan2(yAng,xAng)

      rCyl=self%surf%aSurf*(1+Cm*cos(2*(thAdj+self%surf%thm)))/(3*Cn*cos(3*(thAdj+self%surf%thn))) &
           *(cos(rAng)+sqrt(3.0)*sin(rAng)-1)
    end if

    Rc=self%surf%sHorz-(self%surf%sHorz-self%surf%Rmaj)*rho0**2
    Zc=self%surf%sVert-(self%surf%sVert)*rho0**2

    Rpos=Rc+rCyl*cos(thAdj)
    Zpos=Zc+rCyl*sin(thAdj)

  end subroutine geo_global

  !> Calculate the cylindrical \((R, Z)\) at \((r, \theta)\) for the
  !> generalised elongation model
  !>
  !> See section 5.1.2 of Ball Oxford PhD thesis or section 2.1.3 of
  !> "GS2 analytic geometry specification"
  pure subroutine geo_generalizedElongation(self, r, theta, Rpos, Zpos)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real, intent (out) :: Rpos, Zpos
    real :: dr, thAdj, rCyl, rCylp, Rfinal, Zfinal, Rfinalp, Zfinalp

    dr = r - self%surf%r

    thAdj = theta + self%thetaShift

    rCyl=self%surf%r * (1 &
        + (-1+self%surf%delm/sqrt(1+(self%surf%delm**2-1)*cos(self%surf%mMode*(thAdj+self%surf%thm)/2.0)**2)) &
        + (-1+self%surf%deln/sqrt(1+(self%surf%deln**2-1)*cos(self%surf%nMode*(thAdj+self%surf%thn)/2.0)**2)))
    rCylp=1 &
     +(-1+(1+(self%surf%delm**2-1)*cos(self%surf%mMode*(thAdj+self%surf%thm)/2.0)**2)**(-0.5)*(self%surf%delm+self%surf%r*self%surf%delmp &
      *(1-(self%surf%delm*cos(self%surf%mMode*(thAdj+self%surf%thm)/2.0))**2/(1+(self%surf%delm**2-1)*cos(self%surf%mMode*(thAdj+self%surf%thm)/2.0)**2)))) &
     +(-1+(1+(self%surf%deln**2-1)*cos(self%surf%nMode*(thAdj+self%surf%thn)/2)**2)**(-0.5)*(self%surf%deln+self%surf%r*self%surf%delnp &
      *(1-(self%surf%deln*cos(self%surf%nMode*(thAdj+self%surf%thn)/2.0))**2/(1+(self%surf%deln**2-1)*cos(self%surf%nMode*(thAdj+self%surf%thn)/2.0)**2))))

    Rfinal=rCyl*cos(thAdj)+self%surf%Rmaj
    Zfinal=rCyl*sin(thAdj)
    Rfinalp=rCylp*cos(thAdj)+self%surf%sHorz
    Zfinalp=rCylp*sin(thAdj)+self%surf%sVert

    Rpos=Rfinal+dr*Rfinalp
    Zpos=Zfinal+dr*Zfinalp

  end subroutine geo_generalizedElongation

  !> Calculate the cylindrical \((R, Z)\) at \((r, \theta)\) for the
  !> Fourier model
  !>
  !> See section 5.1.1 of Ball Oxford PhD thesis or section 2.1.4 of
  !> "GS2 analytic geometry specification"
  pure subroutine geo_FourierSeries(self, r, theta, Rpos, Zpos)
    implicit none
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real, intent (out) :: Rpos, Zpos
    real :: dr, thAdj, rCyl, rCylp, Rfinal, Zfinal, Rfinalp, Zfinalp

    dr = r - self%surf%r

    thAdj = theta + self%thetaShift

    rCyl=self%surf%r * (1+((self%surf%delm-1)/2) * (1-cos(self%surf%mMode*(thAdj+self%surf%thm))) + ((self%surf%deln-1)/2) * (1-cos(self%surf%nMode*(thAdj+self%surf%thn))))
    rCylp=1 + (((self%surf%delm-1)/2) + (self%surf%r*self%surf%delmp/2)) * (1-cos(self%surf%mMode*(thAdj+self%surf%thm))) + (((self%surf%deln-1)/2) + (self%surf%r*self%surf%delnp/2)) * (1-cos(self%surf%nMode*(thAdj+self%surf%thn)))

    Rfinal=rCyl*cos(thAdj)+self%surf%Rmaj
    Zfinal=rCyl*sin(thAdj)
    Rfinalp=rCylp*cos(thAdj)+self%surf%sHorz
    Zfinalp=rCylp*sin(thAdj)+self%surf%sVert

    Rpos=Rfinal+dr*Rfinalp
    Zpos=Zfinal+dr*Zfinalp

  end subroutine geo_FourierSeries

  !> Calculate the cylindrical \((R, Z)\) at \((r, \theta)\) for the
  !> Miller Extended Harmonic model
  !>
  !> See [PPCF 63 (2021) 012001 (5pp)](https://doi.org/10.1088/1361-6587/abc63b)
  pure subroutine geo_MillerExtendedHarmonic(self, r, theta, Rpos, Zpos)
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    real, intent (out) :: Rpos, Zpos

    integer :: n
    real :: theta_R, dthetaR_dr, dr, dZdr, dRdr

    theta_R = theta &
         + sum([(self%surf%c_mxh(n) * cos((n - 1) * theta) &
         + self%surf%s_mxh(n) * sin((n - 1) * theta), &
                n = 1, self%surf%n_mxh)])

    dthetaR_dr = sum([(&
         (self%surf%dc_mxh_dr(n) * cos((n - 1) * theta)) &
         + (self%surf%ds_mxh_dr(n) * sin((n - 1) * theta)), &
         n = 1, self%surf%n_mxh)])

    dr = r - self%surf%r

    dRdr = self%surf%sHorz + cos(theta_R) - self%surf%r * sin(theta_R) * dthetaR_dr
    dZdr = self%surf%sVert + self%surf%delm * sin(theta) * (1 + self%surf%delmp)

    Rpos = self%surf%Rmaj + (self%surf%r * cos(theta_R)) + (dr * dRdr)
    Zpos = (self%surf%r * self%surf%delm * sin(theta)) + (dr * dZdr)

  end subroutine geo_MillerExtendedHarmonic

  !> FIXME : Add documentation
  pure real function psi (self, r, theta)
    use warning_helpers, only: is_zero, exactly_equal
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    if (is_zero(r)) then
       psi = self%psi_0
    else if (exactly_equal(r, 1.) .and. is_zero(theta)) then
       psi = self%psi_a
    else
       psi = r - self%surf%r
    end if
  end function psi

  !> FIXME : Add documentation | Should this be equivalent to psi as in others?
  pure real function rfun (self, r, theta)
    class(leq_type), intent(in) :: self
    real, intent (in) :: r, theta
    UNUSED_DUMMY(self); UNUSED_DUMMY(theta)    
    rfun = r
  end function rfun

  !> FIXME : Add documentation
  pure real function diameter (self, rp)
    class(leq_type), intent(in) :: self
    real, intent(in) :: rp
    UNUSED_DUMMY(self)
    diameter = 2. * rp
  end function diameter

  !> Return the major radius of the centre of the flux surface
  !> @Note rcenter for other geo types returns the major radius of the LCFS
  !> which is named R_geo in geometry. We do store R_geo here so could return that?
  pure real function rcenter(self, rp)
    class(leq_type), intent(in) :: self
    real, intent(in) :: rp
    UNUSED_DUMMY(rp)    
    ! If using the global geometry, Rmaj is the center of the shaped flux surface
    ! (i.e. the surface at rhoc=aSurf), so here we calculate the major radial
    ! location of the center of the flux surface of interest ! JRB
    rcenter = self%surf%Rmaj
    if (self%surf%geoType == geo_type_global) then
       rcenter = rcenter + (1 - (self%surf%r / self%surf%aSurf)**2) * self%surf%sHorz
    end if
  end function rcenter

  !> FIXME : Add documentation
  pure real function dbtori (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(self); UNUSED_DUMMY(pbar)
    dbtori = 0
  end function dbtori

  !> FIXME : Add documentation
  pure real function btori (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(pbar)    
    btori = self%surf%R_geo
  end function btori

  !> FIXME : Add documentation
  pure real function qfun (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(pbar)    
    qfun = self%surf%q
  end function qfun

  !> FIXME : Add documentation
  pure real function pfun (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(pbar)    
    pfun = 0.5 * self%beta_0
  end function pfun

  !> FIXME : Add documentation
  pure real function dpfun (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(self); UNUSED_DUMMY(pbar)    
    dpfun = 0
  end function dpfun

  !> FIXME : Add documentation
  pure real function betafun (self, pbar)
    class(leq_type), intent(in) :: self
    real, intent(in) :: pbar
    UNUSED_DUMMY(pbar)
    betafun = self%beta_0
  end function betafun
end module leq