module mms_neutrals_circular_m !! MMS solutions and sources for braginskii model in circular geometry use precision_grillix_m, only : GP, MPI_GP use error_handling_grillix_m, only: handle_error use status_codes_grillix_m, only : GRILLIX_ERR_OTHER use constants_m, only : Pi, TWO_PI use screen_io_m, only : get_stdout use circular_equilibrium_m, only : circular_t use equilibrium_m, only : equilibrium_t use rate_coeff_neutrals_m, only: rate_coeff_neutrals_t use mms_braginskii_circular_m, only : mms_sol_braginskii_circular_ne, & mms_sol_braginskii_circular_te, & mms_sol_braginskii_circular_ti, & mms_sol_braginskii_circular_upar ! Parameters use params_brag_model_m, only : & tratio use params_neut_model_m, only : & s_cx, pardiss_coeff_dens, pardiss_coeff_parmom, pardiss_coeff_pressure use params_neut_switches_m implicit none public :: mms_sol_circular_neutrals_dens public :: mms_sol_circular_neutrals_parmom public :: mms_sol_circular_neutrals_pressure public :: mms_source_circular_neutrals_dens public :: mms_source_circular_neutrals_parmom public :: mms_source_circular_neutrals_pressure private :: cart_to_circ real(GP), parameter, private :: twopi = TWO_PI !! Amplitudes of MMS solutions real(GP), private, parameter :: ampneutralsdens = 0.8_GP real(GP), private, parameter :: ampneutralsparmom = 1.1_GP real(GP), private, parameter :: ampneutralspressure = 1.2_GP real(GP), private, parameter :: ampne = 1.0_GP real(GP), private, parameter :: ampti = 0.9_GP !! Offset of MMS solutions (for positively defined quantites, i.e. density and temperatures) real(GP), private, parameter :: offsetneutralsdens = ampneutralsdens + 0.3_GP real(GP), private, parameter :: offsetneutralspressure = ampneutralspressure + 1.4E-1_GP real(GP), private, parameter :: offsetne = ampne + 1.5E-1_GP real(GP), private, parameter :: offsetti = ampti + 1.2E-1_GP !! Radial mode number of MMS Solutions real(GP), private, parameter :: krneutralsdens = 2.0_GP real(GP), private, parameter :: krneutralsparmom = 1.0_GP real(GP), private, parameter :: krneutralspressure = 1.0_GP real(GP), private, parameter :: krne = 1.0_GP real(GP), private, parameter :: krti = 2.0_GP !! Poloidal mode number of MMS Solutions real(GP), private, parameter :: kpneutralsdens = 1.0_GP real(GP), private, parameter :: kpneutralsparmom = 1.0_GP real(GP), private, parameter :: kpneutralspressure = 2.0_GP real(GP), private, parameter :: kpne = 2.0_GP real(GP), private, parameter :: kpti = 2.0_GP !! Poloidal phase shift of MMS Solutions real(GP), private, parameter :: phpneutralsdens = 0.3_GP real(GP), private, parameter :: phpneutralsparmom = 0.5_GP real(GP), private, parameter :: phpneutralspressure = 0.4_GP real(GP), private, parameter :: phpne = 0.1_GP real(GP), private, parameter :: phpti = 0.4_GP !! Axial mode number of MMS Solutions real(GP), private, parameter :: kzneutralsdens = 1.0_GP real(GP), private, parameter :: kzneutralsparmom = 1.0_GP real(GP), private, parameter :: kzneutralspressure = 1.0_GP real(GP), private, parameter :: kzne = 1.0_GP real(GP), private, parameter :: kzti = 1.0_GP !! Axial phase shift of MMS Solutions real(GP), private, parameter :: phzneutralsdens = 0.2_GP real(GP), private, parameter :: phzneutralsparmom = 0.4_GP real(GP), private, parameter :: phzneutralspressure = 0.3_GP real(GP), private, parameter :: phzne = 1.1_GP real(GP), private, parameter :: phzti = 0.8_GP !! Angular time frequency of MMS Solutions real(GP), private, parameter :: omeganeutralsdens = 55.0_GP real(GP), private, parameter :: omeganeutralsparmom = 65.0_GP real(GP), private, parameter :: omeganeutralspressure = 73.0_GP real(GP), private, parameter :: omegane = 59.0_GP real(GP), private, parameter :: omegati = 61.0_GP !! Temporal phase shift of MMS Solutions real(GP), private, parameter :: phtneutralsdens = 0.15_GP real(GP), private, parameter :: phtneutralsparmom = 0.3_GP real(GP), private, parameter :: phtneutralspressure = 0.26_GP real(GP), private, parameter :: phtne = 0.4_GP real(GP), private, parameter :: phtti = 0.6_GP contains real(GP) function mms_sol_circular_neutrals_dens(equi, x, y, z, t) !! MMS solution for neutrals density class(equilibrium_t), intent(in) :: 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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_neutrals_dens_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_neutrals_dens.txt' mms_sol_circular_neutrals_dens = res end function real(GP) function mms_sol_circular_neutrals_parmom(equi, x, y, z, t) !! MMS solution for neutrals parallel momentum class(equilibrium_t), intent(in) :: 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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_neutrals_parmom_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_neutrals_parmom.txt' mms_sol_circular_neutrals_parmom = res end function real(GP) function mms_sol_circular_neutrals_pressure(equi, x, y, z, t) !! MMS solution for neutrals temperature class(equilibrium_t), intent(in) :: 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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_neutrals_pressure_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_neutrals_pressure.txt' mms_sol_circular_neutrals_pressure = res end function real(GP) function mms_source_circular_neutrals_dens(equi, x, y, z, t, & chi, epsinv, neutrals_dens_pen, & k_iz_o, k_rec_o) !! MMS source for neutral density equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_dens_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object real(GP) :: r, p, q, rmin, rmax, res real(GP) :: kiz, krec, kcx, scx, pardissdens real(GP) :: swdiffdenstor, swdiffdenspol, swparfluxdens real(GP) :: tiplasma, teplasma, neplasma, neutrals_dens real(GP) :: neutralsdenspen include 'mms_source_circular_neutrals_dens_head.txt' swdiffdenstor = swn_dens_diff_tor swdiffdenspol = swn_dens_diff_pol swparfluxdens = swn_dens_par_flux pardissdens = pardiss_coeff_dens neplasma = mms_sol_braginskii_circular_ne(equi, x, y, z, t) tiplasma = mms_sol_braginskii_circular_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_circular_te(equi, x, y, z, t) neutrals_dens = mms_sol_circular_neutrals_dens(equi, x, y, z, t) scx = s_cx kiz = swn_iz * k_iz_o%eval(teplasma, neplasma) krec = swn_rec * k_rec_o%eval(teplasma, neplasma) kcx = 2.93_GP * scx * sqrt(tiplasma * tratio) neutralsdenspen = neutrals_dens_pen ! 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_neutrals_dens.txt' mms_source_circular_neutrals_dens = res end function real(GP) function mms_source_circular_neutrals_parmom(equi, x, y, z, t, & chi, epsinv, neutrals_parmom_pen, k_iz_o, k_rec_o) !! MMS source for neutral parallel momentum equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_parmom_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object real(GP) :: r, p, q, rmin, rmax, res real(GP) :: kiz, krec, kcx, scx, pardissparmom real(GP) :: swdiffparmompol, swparfluxparmom, swparpressure, swsrcparmom real(GP) :: tiplasma, teplasma, neplasma, neutrals_dens, uparplasma real(GP) :: neutralsparmompen include 'mms_source_circular_neutrals_parmom_head.txt' swdiffparmompol = swn_parmom_diff_pol swparfluxparmom = swn_parmom_par_flux swparpressure = swn_parmom_par_pressure swsrcparmom = swn_src_parmom pardissparmom = pardiss_coeff_parmom tiplasma = mms_sol_braginskii_circular_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_circular_te(equi, x, y, z, t) neplasma = mms_sol_braginskii_circular_ne(equi, x, y, z, t) uparplasma = mms_sol_braginskii_circular_upar(equi, x, y, z, t) scx = s_cx kiz = swn_iz * k_iz_o%eval(teplasma, neplasma) krec = swn_rec * k_rec_o%eval(teplasma, neplasma) kcx = 2.93_GP * scx * sqrt(tiplasma * tratio) neutralsparmompen = neutrals_parmom_pen ! 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_neutrals_parmom.txt' mms_source_circular_neutrals_parmom = res end function real(GP) function mms_source_circular_neutrals_pressure(equi, x, y, z, t, & chi, epsinv, neutrals_pressure_pen, k_iz_o, k_rec_o) !! MMS source for neutral temperature equation class(equilibrium_t), intent(in) :: 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: epsinv !! Inverse of penalisation epsilon real(GP), intent(in) :: neutrals_pressure_pen !! Target value to be penalised towards type(rate_coeff_neutrals_t), intent(in) :: k_iz_o !! Ionization rate coefficient object type(rate_coeff_neutrals_t), intent(in) :: k_rec_o !! Recombination rate coefficient object real(GP) :: r, p, q, rmin, rmax, res real(GP) :: kiz, krec, kcx, scx, pardisspressure real(GP) :: tiplasma, teplasma, neplasma, neutrals_dens, uparplasma real(GP) :: neutralspressurepen real(GP) :: swpressurediffpol, swpressureparflux, swpressuredivvpar, & swpressurevischeat, swpressuresrc include 'mms_source_circular_neutrals_pressure_head.txt' tiplasma = mms_sol_braginskii_circular_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_circular_te(equi, x, y, z, t) neplasma = mms_sol_braginskii_circular_ne(equi, x, y, z, t) uparplasma = mms_sol_braginskii_circular_upar(equi, x, y, z, t) swpressurediffpol = swn_pressure_diff_pol swpressureparflux = swn_pressure_par_flux swpressuredivvpar = swn_pressure_div_vpar swpressurevischeat = swn_pressure_vischeat swpressuresrc = swn_src_pressure pardisspressure = pardiss_coeff_pressure scx = s_cx kiz = swn_iz * k_iz_o%eval(teplasma, neplasma) krec = swn_rec * k_rec_o%eval(teplasma, neplasma) kcx = 2.93_GP * scx * sqrt(tiplasma * tratio) neutralspressurepen = neutrals_pressure_pen ! 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_neutrals_pressure.txt' mms_source_circular_neutrals_pressure = 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(in) :: 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