! ============================================================================= !> This submodule defines internal kink modes in force-free magnetic fields. !! The geometry is cylindrical with parabolic density and velocity profiles, !! The geometry can be overridden in the parfile. !! !! This equilibrium is taken from section III.B in !! _Goedbloed, J. P. "The Spectral Web of stationary plasma equilibria. !! II. Internal modes." Physics of Plasmas 25.3 (2018): 032110_. !! @note Default values are given by !! !! - k2 = 1 !! - k3 = \( 0.16\alpha \) !! - cte_rho0 = 1 : used as prefactor in setting the density. !! - cte_v03 = 1 : used as prefactor in setting the z-component of velocity. !! - cte_p0 = 3 : used to set the pressure. !! - alpha = 5 / x_end : used in the Bessel functions. !! !! and can all be changed in the parfile. @endnote submodule (mod_equilibrium) smod_equil_internal_kink_instability use mod_equilibrium_params, only: cte_rho0, cte_v03, cte_p0, alpha implicit none real(dp) :: a0 contains module procedure internal_kink_eq call settings%grid%set_geometry("cylindrical") call settings%grid%set_grid_boundaries(0.0_dp, 1.0_dp) a0 = settings%grid%get_grid_end() if (settings%equilibrium%use_defaults) then ! LCOV_EXCL_START call settings%physics%enable_flow() cte_rho0 = 1.0_dp cte_v03 = 1.0_dp cte_p0 = 9.0_dp alpha = 5.0_dp / a0 k2 = 1.0_dp k3 = 0.16_dp * alpha end if ! LCOV_EXCL_STOP call background%set_density_funcs(rho0_func=rho0, drho0_func=drho0) call background%set_velocity_3_funcs(v03_func=v03, dv03_func=dv03) call background%set_temperature_funcs(T0_func=T0, dT0_func=dT0) call background%set_magnetic_2_funcs(B02_func=B02, dB02_func=dB02) call background%set_magnetic_3_funcs(B03_func=B03, dB03_func=dB03) end procedure internal_kink_eq real(dp) function rho0(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 rho0 = cte_rho0 * (1.0_dp - x**2 / a0**2) end function rho0 real(dp) function drho0(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 drho0 = -2.0_dp * cte_rho0 * x / a0 end function drho0 real(dp) function T0(r) real(dp), intent(in) :: r T0 = cte_p0 / rho0(r) end function T0 real(dp) function dT0(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 dT0 = 2.0_dp * x * cte_p0 / (a0**2 * cte_rho0 * (1.0_dp - x**2)**2) end function dT0 real(dp) function v03(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 v03 = cte_v03 * (1.0_dp - x**2 / a0**2) end function v03 real(dp) function dv03(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 dv03 = -2.0_dp * cte_v03 * x / a0 end function dv03 real(dp) function B02(r) real(dp), intent(in) :: r B02 = J1(r) end function B02 real(dp) function dB02(r) real(dp), intent(in) :: r dB02 = dJ1(r) end function dB02 real(dp) function B03(r) real(dp), intent(in) :: r B03 = J0(r) end function B03 real(dp) function dB03(r) real(dp), intent(in) :: r dB03 = dJ0(r) end function dB03 real(dp) function J0(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 J0 = bessel_jn(0, alpha * x) end function J0 real(dp) function J1(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 J1 = bessel_jn(1, alpha * x) end function J1 real(dp) function J2(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 J2 = bessel_jn(2, alpha * x) end function J2 real(dp) function dJ0(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 dJ0 = -alpha * J1(r) end function dJ0 real(dp) function dJ1(r) real(dp), intent(in) :: r real(dp) :: x x = r / a0 dJ1 = alpha * (0.5_dp * J0(r) - 0.5_dp * J2(r)) end function dJ1 end submodule smod_equil_internal_kink_instability