! ============================================================================= !> This submodule defines a magnetic flux tube embedded in a uniform magnetic !! environment. In this case the flux tube is under photospheric conditions !! \( c_{Ae} < c_s < c_{se} < c_A \) where the subscript e denotes the outer region. !! More specifically the equilibrium is defined as !! \( c_{Ae} = c_s/2, c_{se} = 3c_s/2, c_A = 2c_s \). The geometry can be overridden !! in the parfile, and is cylindrical by default for \( r \in [0, 10] \). !! !! This equilibrium is taken from chapter 6, fig. 6.5 in !! _Roberts, Bernard (2019). MHD Waves in the Solar Atmosphere. !! Cambridge University Press._ [DOI](https://doi.org/10.1017/9781108613774). !! @note For best results, it is recommended to enable mesh accumulation. @endnote !! @note Default values are given by !! !! - k2 = 0 !! - k3 = 2 !! - cte_rho0 = 1 : density value for the inner tube. !! - cte_p0 = 1 : pressure value for the inner tube. !! - r0 = 1 : radius of the inner tube. !! !! and can all be changed in the parfile. @endnote ! SUBMODULE: smod_equil_coronal_flux_tube submodule(mod_equilibrium) smod_equil_photospheric_flux_tube use mod_equilibrium_params, only: cte_rho0, cte_p0, r0 implicit none real(dp) :: rho_e, p_e, B_0, B_e contains module procedure photospheric_flux_tube_eq use mod_global_variables, only: dp_LIMIT real(dp), allocatable :: custom_grid(:) real(dp) :: width, a_l, a_r, pct1, pct2, pct3, dx, dx1, dx2, dx3 real(dp) :: gamma real(dp) :: x_start, x_end integer :: gridpts integer :: i, N1, N2, N3 if (settings%equilibrium%use_defaults) then ! LCOV_EXCL_START call settings%grid%set_geometry("cylindrical") x_start = 0.0_dp x_end = 10.0_dp call settings%grid%set_grid_boundaries(x_start, x_end) cte_rho0 = 1.0_dp cte_p0 = 1.0_dp r0 = 1.0_dp k2 = 0.0_dp k3 = 2.0_dp end if ! LCOV_EXCL_STOP gamma = settings%physics%get_gamma() gridpts = settings%grid%get_gridpts() x_start = settings%grid%get_grid_start() x_end = settings%grid%get_grid_end() ! width of transition region width = 0.1_dp ! start of transition region a_l = r0 - width / 2.0_dp ! end of transition region a_r = a_l + width pct1 = 0.6_dp ! percentage of points in inner tube pct2 = 0.3_dp ! percentage of points in transition region pct3 = 0.1_dp ! percentage of points in outer tube ! for a custom grid we have to manually avoid setting r0 = 0 if ( & settings%grid%get_geometry() == "cylindrical" & .and. .not. settings%grid%force_r0 & ) then x_start = 2.5e-2_dp end if ! sanity checks if (r0 > x_end) then call logger%error("equilibrium: inner cylinder radius r0 > x_end") else if (r0 < x_start) then call logger%error("equilibrium: inner cylinder radius r0 < x_start") end if N1 = int(pct1 * gridpts) dx1 = (a_l - x_start) / (N1 - 1) ! -1 since first point is x0 N2 = int(pct2 * gridpts) dx2 = (a_r - a_l) / N2 N3 = int(pct3 * gridpts) dx3 = (x_end - a_r) / N3 ! fill the grid allocate(custom_grid(gridpts)) custom_grid(1) = x_start do i = 2, gridpts if (i <= N1) then dx = dx1 else if (N1 < i .and. i <= N1 + N2) then dx = dx2 else dx = dx3 end if custom_grid(i) = custom_grid(i - 1) + dx end do call grid%set_custom_grid(custom_grid) deallocate(custom_grid) rho_e = 8.0_dp * (2.0_dp * gamma + 1.0_dp) * cte_rho0 / (gamma + 18.0_dp) p_e = 18.0_dp * (2.0_dp * gamma + 1.0_dp) * cte_p0 / (gamma + 18.0_dp) B_0 = 2.0_dp * sqrt(gamma * cte_p0) B_e = sqrt(2.0_dp * gamma * cte_p0 * (2.0_dp * gamma + 1.0_dp) / (gamma + 18.0_dp)) ! check pressure balance if (abs(cte_p0 + 0.5_dp * B_0**2 - p_e - 0.5_dp * B_e**2) > dp_LIMIT) then call logger%error("equilibrium: total pressure balance not satisfied") end if call background%set_density_funcs(rho0_func=rho0) call background%set_temperature_funcs(T0_func=T0) call background%set_magnetic_3_funcs(B03_func=B03) end procedure photospheric_flux_tube_eq real(dp) function rho0(r) real(dp), intent(in) :: r if (r > r0) then rho0 = rho_e else rho0 = cte_rho0 end if end function rho0 real(dp) function T0(r) real(dp), intent(in) :: r if (r > r0) then T0 = p_e / rho_e else T0 = cte_p0 / cte_rho0 end if end function T0 real(dp) function B03(r) real(dp), intent(in) :: r if (r > r0) then B03 = B_e else B03 = B_0 end if end function B03 end submodule smod_equil_photospheric_flux_tube