smod_equil_kelvin_helmholtz_cd.f08 Source File

! =============================================================================
!> This submodule defines an unperturbed magnetised jet model in cylindrical
!! geometry, giving rise to Kelvin-Helmholtz and current-driven instabilities.
!! The geometry is fixed for this problem; the cylinder wall is dependent
!! on the equilibrium parameters and is given by <tt>2rj</tt>.
!!
!! This equilibrium is taken from
!! _Baty, H., & Keppens, R. (2002). Interplay between Kelvin-Helmholtz and
!!  current-driven instabilities in jets. The Astrophysical Journal, 580(2), 800_.
!! @note Default values are given by
!!
!! - <tt>k2</tt> = -1
!! - <tt>V</tt> = 1.63 : amplitude of the velocity shear
!! - <tt>cte_p0</tt> = 1 : used to set the pressure.
!! - <tt>cte_rho0</tt> = 1 : used to set the density.
!! - <tt>Bz0</tt> = 0.25 : used to set Bz.
!! - <tt>rc</tt> = 0.5 : length for radial variation.
!! - <tt>rj</tt> = 1 : jet radius
!!
!! and can all be changed in the parfile. @endnote
!! @note The default setup is _HEL2_ in the original paper.
!!       For _HEL1_ you can set <tt>rc = 2</tt>. @endnote
submodule (mod_equilibrium) smod_equil_kelvin_helmholtz_cd
  implicit none

contains

  !> Sets the equilibrium.
  module subroutine kh_cd_instability_eq()
    use mod_equilibrium_params, only: V, cte_rho0, cte_p0, Bz0, rc, Bth0, rj

    real(dp)    :: r, a
    integer     :: i

    if (use_defaults) then ! LCOV_EXCL_START
      flow = .true.

      V = 1.63d0
      cte_rho0 = 1.0d0
      cte_p0 = 1.0d0
      Bz0 = 0.25d0
      rj = 1.0d0

      rc    = 0.5d0
      k2  = -1.0d0
    end if ! LCOV_EXCL_STOP

    Bth0  = 0.4d0 * (rc**2+rj**2) / (rj*rc)
    a = 0.1d0 * rj
    k3  = dpi / rj

    geometry = "cylindrical"
    x_start = 0.0d0
    x_end   = 2.0d0 * rj
    call initialise_grid()

    do i = 1, gauss_gridpts
      r = grid_gauss(i)

      rho_field % rho0(i) = cte_rho0
      v_field % v03(i)    = (V/2.0d0) * tanh((rj-r)/a)
      B_field % B02(i)    = Bth0 * r*rc / (rc**2 + r**2 )
      B_field % B03(i)    = Bz0
      B_field % B0(i)     = sqrt((B_field % B02(i))**2 + (B_field % B03(i))**2)
      T_field % T0(i)     = cte_p0 / (rho_field % rho0(i)) &
        - (Bth0**2/(2.0d0*(rho_field % rho0(i)))) * (1 - rc**4/(rc**2+r**2)**2)

      B_field % d_B02_dr(i) = Bth0 * rc * (rc**2-r**2) / (r**2+rc**2)**2
      v_field % d_v03_dr(i) = - (V/(2.0d0*a)) / cosh((rj-r)/a)**2
      T_field % d_T0_dr(i)  = - (2.0d0*Bth0**2/(rho_field % rho0(i))) &
        * rc**4*r / (r**2+rc**2)**3
    end do
  end subroutine kh_cd_instability_eq

end submodule smod_equil_kelvin_helmholtz_cd