module mms_diffusion_circular_m !! MMS Solutions and Sources for diffusion model in circular geometry use precision_grillix_m, only : GP use error_handling_grillix_m, only: handle_error, error_info_t use status_codes_grillix_m, only : GRILLIX_ERR_OTHER, GRILLIX_ERR_NAMELIST use params_diffusion_m, only : diffcoeff_perp, diffcoeff_par, epsinv use constants_m, only : Pi, TWO_PI use screen_io_m, only : get_stdout use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t implicit none public :: mms_sol_diffusion_circular public :: mms_source_diffusion_circular private :: cart_to_circ real(GP), private :: amp = 1.0_GP !! Amplitude of MMS solution real(GP), private :: krn = 1.0_GP !! Radial mode number of MMS Solution real(GP), private :: kp = 2.0_GP !! Poloidal mode number of MMS Solution real(GP), private :: php = 0.2_GP !! Poloidal phase shift of MMS Solution real(GP), private :: kz = 1.0_GP !! Axial mode number of MMS Solution real(GP), private :: phz = 1.1_GP !! Axial phase shift of MMS Solution real(GP), private :: omega = 0.5_GP !! Angular time frequency of MMS Solution real(GP), private :: pht = 0.4_GP !! Temporal phase shift of MMS Solution contains real(GP) function mms_sol_diffusion_circular(equi, x, y, z, t) !! MMS solution class(equilibrium_t) :: equi !! Equilibrium real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time real(GP) :: res, r, p, q, rmin, rmax real(GP), parameter :: twopi = TWO_PI include 'mms_sol_circular_head.txt' ! Convert Cartesian to polar coordinates and safety factor call cart_to_circ(equi, x, y, rho = r, theta = p, rhomin = rmin, rhomax = rmax, qval = q) include 'mms_sol_circular.txt' mms_sol_diffusion_circular = res end function real(GP) function mms_source_diffusion_circular(equi, x, y, z, t, charfun, penval) !! Source corresponding to MMS solution class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(in) :: z !! Axial coordinate z (represented by toroidal coordinate) real(GP), intent(in) :: t !! Time real(GP), intent(in) :: charfun !! Value of characteristic function (penalisation) real(GP), intent(in) :: penval !! Value that solution shall be penalised to real(GP), parameter :: twopi = TWO_PI real(GP) :: res, r, p, q, rmin, rmax real(GP) :: chiperp, chipar, penepsinv include 'mms_source_circular_head.txt' chiperp = diffcoeff_perp chipar = diffcoeff_par penepsinv = epsinv ! Convert Cartesian to polar coordinates and safety factor call cart_to_circ(equi, x, y, rho = r, theta = p, rhomin = rmin, rhomax = rmax, qval = q) include 'mms_source_circular.txt' mms_source_diffusion_circular = res end function subroutine cart_to_circ(equi, x, y, rho, theta, rhomin, rhomax, qval) !! Converts Cartesian (x,y) coordinates to (rho, theta) coordinates and returns equilibrium specific quantities class(equilibrium_t), intent(inout) :: equi !! Equilibrium real(GP), intent(in) :: x !! x-coordinate real(GP), intent(in) :: y !! y-coordinate real(GP), intent(out) :: rho !! radial coordinate real(GP), intent(out) :: theta !! poloidal angle real(GP), intent(out) :: rhomin !! Inner limiting flux surface real(GP), intent(out) :: rhomax !! Outer limiting flux surface real(GP), intent(out) :: qval !! safety factor real(GP) :: phi select type(equi) type is(circular_t) phi = 0.0_GP ! Axisymmetry rho = equi%rho(x, y, phi) theta = equi%theta(x, y) rhomin = equi%rhomin rhomax = equi%rhomax qval = equi%qval(rho) class default call handle_error('Equilibrium is not circular', & GRILLIX_ERR_OTHER, __LINE__, __FILE__) end select end subroutine end module