module mms_neutrals_slab_m !! MMS solutions and sources for Neutrals model in slab geometry use precision_grillix_m, only : GP, MPI_GP use constants_m, only : Pi, TWO_PI use slab_equilibrium_m, only : slab_t use equilibrium_m, only : equilibrium_t use rate_coeff_neutrals_m, only: rate_coeff_neutrals_t use mms_braginskii_slab_m, only : mms_sol_braginskii_slab_ne, & mms_sol_braginskii_slab_te, & mms_sol_braginskii_slab_ti, & mms_sol_braginskii_slab_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_slab_neutrals_dens public :: mms_sol_slab_neutrals_parmom public :: mms_sol_slab_neutrals_pressure public :: mms_source_slab_neutrals_dens public :: mms_source_slab_neutrals_parmom public :: mms_source_slab_neutrals_pressure 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 + 1.1E-1_GP real(GP), private, parameter :: offsetneutralspressure = ampneutralspressure + 0.9E-1_GP real(GP), private, parameter :: offsetne = ampne + 1.5E-1_GP real(GP), private, parameter :: offsetti = ampti + 1.7E-1_GP !! Radial mode number of MMS Solutions real(GP), private, parameter :: kxneutralsdens = 2.0_GP real(GP), private, parameter :: kxneutralsparmom = 3.0_GP real(GP), private, parameter :: kxneutralspressure = 1.0_GP real(GP), private, parameter :: kxne = 1.0_GP real(GP), private, parameter :: kxti = 2.0_GP !! Poloidal mode number of MMS Solutions real(GP), private, parameter :: kyneutralsdens = 1.0_GP real(GP), private, parameter :: kyneutralsparmom = 2.0_GP real(GP), private, parameter :: kyneutralspressure = 3.0_GP real(GP), private, parameter :: kyne = 2.0_GP real(GP), private, parameter :: kyti = 2.0_GP !! Poloidal phase shift of MMS Solutions real(GP), private, parameter :: phyneutralsdens = 0.3_GP real(GP), private, parameter :: phyneutralsparmom = 0.5_GP real(GP), private, parameter :: phyneutralspressure = 0.4_GP real(GP), private, parameter :: phyne = 0.1_GP real(GP), private, parameter :: phyti = 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 = 1.5_GP real(GP), private, parameter :: phzneutralspressure = 0.7_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 = 72.0_GP real(GP), private, parameter :: omeganeutralspressure = 66.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.6_GP real(GP), private, parameter :: phtneutralspressure = 0.42_GP real(GP), private, parameter :: phtne = 0.4_GP real(GP), private, parameter :: phtti = 0.6_GP contains real(GP) function mms_sol_slab_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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_neutrals_dens_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_neutrals_dens.txt' mms_sol_slab_neutrals_dens = res end function real(GP) function mms_sol_slab_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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_neutrals_parmom_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_neutrals_parmom.txt' mms_sol_slab_neutrals_parmom = res end function real(GP) function mms_sol_slab_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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_neutrals_pressure_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_neutrals_pressure.txt' mms_sol_slab_neutrals_pressure = res end function real(GP) function mms_source_slab_neutrals_dens(equi, x, y, z, t, & k_iz_o, k_rec_o) !! MMS source for neutrals 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 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) :: xmin, xmax, ymin, ymax, res real(GP) :: kiz, krec, kcx, scx, pardissdens real(GP) :: swdiffdenstor, swdiffdenspol, swparfluxdens real(GP) :: tiplasma, teplasma, neplasma, neutrals_dens include 'mms_source_slab_neutrals_dens_head.txt' tiplasma = mms_sol_braginskii_slab_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_slab_te(equi, x, y, z, t) neplasma = mms_sol_braginskii_slab_ne(equi, x, y, z, t) neutrals_dens = mms_sol_slab_neutrals_dens(equi, x, y, z, t) swdiffdenstor = swn_dens_diff_tor swdiffdenspol = swn_dens_diff_pol swparfluxdens = swn_dens_par_flux pardissdens = pardiss_coeff_dens 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) xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_neutrals_dens.txt' mms_source_slab_neutrals_dens = res end function real(GP) function mms_source_slab_neutrals_parmom(equi, x, y, z, t, & k_iz_o, k_rec_o) !! MMS source for neutrals 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 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) :: xmin, xmax, ymin, ymax, res real(GP) :: kiz, krec, kcx, scx, pardissparmom real(GP) :: swdiffparmompol, swparfluxparmom, swparpressure, swsrcparmom real(GP) :: tiplasma, teplasma, neplasma, uparplasma include 'mms_source_slab_neutrals_parmom_head.txt' tiplasma = mms_sol_braginskii_slab_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_slab_te(equi, x, y, z, t) neplasma = mms_sol_braginskii_slab_ne(equi, x, y, z, t) uparplasma = mms_sol_braginskii_slab_upar(equi, x, y, z, t) swdiffparmompol = swn_parmom_diff_pol swparfluxparmom = swn_parmom_par_flux swparpressure = swn_parmom_par_pressure swsrcparmom = swn_src_parmom pardissparmom = pardiss_coeff_parmom 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) xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_neutrals_parmom.txt' mms_source_slab_neutrals_parmom = res end function real(GP) function mms_source_slab_neutrals_pressure(equi, x, y, z, t, & k_iz_o, k_rec_o) !! MMS source for neutrals 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 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) :: xmin, xmax, ymin, ymax, res real(GP) :: tiplasma, teplasma, neplasma, uparplasma real(GP) :: kiz, krec, kcx, scx, pardisspressure real(GP) :: swpressurediffpol, swpressureparflux, swpressuredivvpar, & swpressurevischeat, swpressuresrc include 'mms_source_slab_neutrals_pressure_head.txt' tiplasma = mms_sol_braginskii_slab_ti(equi, x, y, z, t) teplasma = mms_sol_braginskii_slab_te(equi, x, y, z, t) neplasma = mms_sol_braginskii_slab_ne(equi, x, y, z, t) uparplasma = mms_sol_braginskii_slab_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) xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_neutrals_pressure.txt' mms_source_slab_neutrals_pressure = res end function end module