module mms_braginskii_slab_m !! MMS solutions and sources for braginskii model in slab geometry use precision_grillix_m, only : GP use constants_m, only : Pi, TWO_PI use slab_equilibrium_m, only : slab_t use equilibrium_m, only : equilibrium_t ! Parameters use params_brag_model_m, only : & rhos, delta, etapar0, eta_i0, heat_fac, tratio, beta, mass_ratio_ei, nu_e0, thermal_force_coeff use params_brag_pardiss_model_m, only : & chipar0_e, chipar0_i, free_streaming_fraction_e, free_streaming_fraction_i, & free_streaming_limiter_qfac, aspect_ratio_inv, & heatflux_model, landau_flutter_lhs_on use params_brag_numdiss_m, only : & perpdiss_coeff_ne, perpdiss_coeff_vort, perpdiss_coeff_upar, perpdiss_coeff_ohm, & perpdiss_coeff_te, perpdiss_coeff_ti, perpdiss_coeff_pe, perpdiss_coeff_pi, & pardiss_coeff_ne, pardiss_coeff_vort use params_brag_buffer_m, only : & buffer_coeff_ne, buffer_coeff_vort, buffer_coeff_upar, buffer_coeff_ohm, & buffer_coeff_te, buffer_coeff_ti, buffer_coeff_pe, buffer_coeff_pi use params_brag_switches_m implicit none public :: mms_sol_braginskii_slab_ne public :: mms_sol_braginskii_slab_pot public :: mms_sol_braginskii_slab_vort public :: mms_sol_braginskii_slab_vortbsq public :: mms_sol_braginskii_slab_upar public :: mms_sol_braginskii_slab_jpar public :: mms_sol_braginskii_slab_apar public :: mms_sol_braginskii_slab_te public :: mms_sol_braginskii_slab_ti public :: mms_sol_braginskii_slab_qe public :: mms_sol_braginskii_slab_qi public :: mms_source_braginskii_slab_continuity public :: mms_source_braginskii_slab_vorticity public :: mms_source_braginskii_slab_vorticitybsq public :: mms_source_braginskii_slab_parmomentum public :: mms_source_braginskii_slab_ohm public :: mms_source_braginskii_slab_etemp public :: mms_source_braginskii_slab_itemp real(GP), parameter, private :: twopi = TWO_PI !! Amplitudes of MMS solutions real(GP), private, parameter :: ampne = 1.0_GP real(GP), private, parameter :: amppot = 0.8_GP real(GP), private, parameter :: ampupar = 1.3_GP real(GP), private, parameter :: ampapar = 0.7_GP real(GP), private, parameter :: ampte = 1.2_GP real(GP), private, parameter :: ampti = 0.9_GP real(GP), private, parameter :: ampqe = 0.5_GP real(GP), private, parameter :: ampqi = 0.6_GP !! Offset of MMS solutions (for positively defined quantites, i.e. density and temperatures) real(GP), private, parameter :: offsetne = ampne + 1.5E-1_GP real(GP), private, parameter :: offsette = ampte + 1.0E-1_GP real(GP), private, parameter :: offsetti = ampti + 1.7E-1_GP !! Radial mode number of MMS Solutions real(GP), private, parameter :: kxne = 1.0_GP real(GP), private, parameter :: kxpot = 1.0_GP real(GP), private, parameter :: kxupar = 2.0_GP real(GP), private, parameter :: kxapar = 1.0_GP real(GP), private, parameter :: kxte = 2.0_GP real(GP), private, parameter :: kxti = 2.0_GP real(GP), private, parameter :: kxqe = 2.0_GP real(GP), private, parameter :: kxqi = 1.0_GP !! Poloidal mode number of MMS Solutions real(GP), private, parameter :: kyne = 2.0_GP real(GP), private, parameter :: kypot = 3.0_GP real(GP), private, parameter :: kyupar = 2.0_GP real(GP), private, parameter :: kyapar = 3.0_GP real(GP), private, parameter :: kyte = 3.0_GP real(GP), private, parameter :: kyti = 2.0_GP real(GP), private, parameter :: kyqe = 3.0_GP real(GP), private, parameter :: kyqi = 2.0_GP !! Poloidal phase shift of MMS Solutions real(GP), private, parameter :: phyne = 0.1_GP real(GP), private, parameter :: phypot = 0.7_GP real(GP), private, parameter :: phyupar = 1.2_GP real(GP), private, parameter :: phyapar = 0.5_GP real(GP), private, parameter :: phyte = 0.3_GP real(GP), private, parameter :: phyti = 0.4_GP real(GP), private, parameter :: phyqe = 0.8_GP real(GP), private, parameter :: phyqi = 0.6_GP !! Axial mode number of MMS Solutions real(GP), private, parameter :: kzne = 1.0_GP real(GP), private, parameter :: kzpot = 1.0_GP real(GP), private, parameter :: kzupar = 1.0_GP real(GP), private, parameter :: kzapar = 1.0_GP real(GP), private, parameter :: kzte = 1.0_GP real(GP), private, parameter :: kzti = 1.0_GP real(GP), private, parameter :: kzqe = 1.0_GP real(GP), private, parameter :: kzqi = 1.0_GP !! Axial phase shift of MMS Solutions real(GP), private, parameter :: phzne = 1.1_GP real(GP), private, parameter :: phzpot = 0.2_GP real(GP), private, parameter :: phzupar = 0.1_GP real(GP), private, parameter :: phzapar = 0.6_GP real(GP), private, parameter :: phzte = 0.5_GP real(GP), private, parameter :: phzti = 0.8_GP real(GP), private, parameter :: phzqe = 0.4_GP real(GP), private, parameter :: phzqi = 0.3_GP !! Angular time frequency of MMS Solutions real(GP), private, parameter :: omegane = 59.0_GP real(GP), private, parameter :: omegapot = 83.0_GP real(GP), private, parameter :: omegaupar = 99.0_GP real(GP), private, parameter :: omegaapar = 66.0_GP real(GP), private, parameter :: omegate = 79.0_GP real(GP), private, parameter :: omegati = 61.0_GP real(GP), private, parameter :: omegaqe = 71.0_GP real(GP), private, parameter :: omegaqi = 51.0_GP !! Temporal phase shift of MMS Solutions real(GP), private, parameter :: phtne = 0.4_GP real(GP), private, parameter :: phtpot = 1.1_GP real(GP), private, parameter :: phtupar = 0.5_GP real(GP), private, parameter :: phtapar = 0.3_GP real(GP), private, parameter :: phtte = 0.1_GP real(GP), private, parameter :: phtti = 0.6_GP real(GP), private, parameter :: phtqe = 0.7_GP real(GP), private, parameter :: phtqi = 0.2_GP contains real(GP) function mms_sol_braginskii_slab_ne(equi, x, y, z, t) !! MMS solution for electron density 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_ne_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_ne.txt' mms_sol_braginskii_slab_ne = res end function real(GP) function mms_sol_braginskii_slab_pot(equi, x, y, z, t) !! MMS solution for electrostatic potential 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_pot_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_pot.txt' mms_sol_braginskii_slab_pot = res end function real(GP) function mms_sol_braginskii_slab_vort(equi, x, y, z, t) !! MMS solution for generalised vorticity 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) :: xmin, xmax, ymin, ymax, res real(GP) :: swvortdia include 'mms_sol_slab_vort_head.txt' swvortdia = swb_vort_dia xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_vort.txt' mms_sol_braginskii_slab_vort = res end function real(GP) function mms_sol_braginskii_slab_vortbsq(equi, x, y, z, t) !! MMS solution for Boussinesq vorticity 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) :: xmin, xmax, ymin, ymax, res real(GP) :: swvortdia include 'mms_sol_slab_vortbsq_head.txt' swvortdia = swb_vort_dia xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_vortbsq.txt' mms_sol_braginskii_slab_vortbsq = res end function real(GP) function mms_sol_braginskii_slab_upar(equi, x, y, z, t) !! MMS solution for parallel ion velocity 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_upar_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_upar.txt' mms_sol_braginskii_slab_upar = res end function real(GP) function mms_sol_braginskii_slab_jpar(equi, x, y, z, t) !! MMS solution for parallel current 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_jpar_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_jpar.txt' mms_sol_braginskii_slab_jpar = res end function real(GP) function mms_sol_braginskii_slab_apar(equi, x, y, z, t) !! MMS solution for parallel electromagnetic potential 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_apar_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_apar.txt' mms_sol_braginskii_slab_apar = res end function real(GP) function mms_sol_braginskii_slab_te(equi, x, y, z, t) !! MMS solution for electron temperature 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_te_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_te.txt' mms_sol_braginskii_slab_te = res end function real(GP) function mms_sol_braginskii_slab_ti(equi, x, y, z, t) !! MMS solution for ion temperature 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 coordinatchipar0ie) real(GP), intent(in) :: t !! Time real(GP) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_ti_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_ti.txt' mms_sol_braginskii_slab_ti = res end function real(GP) function mms_sol_braginskii_slab_qe(equi, x, y, z, t) !! MMS solution for ion heat flux 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_qe_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_qe.txt' mms_sol_braginskii_slab_qe = res end function real(GP) function mms_sol_braginskii_slab_qi(equi, x, y, z, t) !! MMS solution for ion heat flux 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) :: xmin, xmax, ymin, ymax, res include 'mms_sol_slab_qi_head.txt' xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_sol_slab_qi.txt' mms_sol_braginskii_slab_qi = res end function real(GP) function mms_source_braginskii_slab_continuity(equi, x, y, z, t, src_ne) !! MMS source for continuity equation 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) :: src_ne !! Particle source real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: PerpdissCoeffNe, dissparallelne, bufferNe real(GP) :: srcne real(GP) :: swcontexb, & swcontcurvpot, swcontcurvte, swcontcurvne, & swcontdivjpar, swcontdivparflx, & swcontdissperp, swcontdissparallel, swcontsource, & swcontflutterparadv, swcontflutterdivjpar, swcontflutterdivupar include 'mms_source_slab_continuity_head.txt' PerpdissCoeffNe = perpdiss_coeff_ne dissparallelne = pardiss_coeff_ne bufferNe = buffer_coeff_ne srcne = src_ne swcontexb = swb_cont_exb swcontcurvpot = swb_cont_curv_pot swcontcurvte = swb_cont_curv_te swcontcurvne = swb_cont_curv_ne swcontdivjpar = swb_cont_divjpar swcontdivparflx = swb_cont_paradv ! swb_cont_divupar = swb_cont_paradv is currently assumed! swcontdissperp = swb_cont_dissperp swcontdissparallel = swb_cont_dissparallel swcontsource = swb_cont_source swcontflutterparadv = swb_cont_flutter_paradv swcontflutterdivjpar = swb_cont_flutter_divjpar swcontflutterdivupar = swb_cont_flutter_divupar xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_continuity.txt' mms_source_braginskii_slab_continuity = res end function real(GP) function mms_source_braginskii_slab_vorticity(equi, x, y, z, t, src_vort) !! MMS source for continuity equation 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) :: src_vort !! Vorticity source real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: etai0, heatfac, chipar0i, nuperpvort, dissparallelvort, buffervort real(GP) :: swvortdia, swvortexb, swvortparadv real(GP) :: swvortnecurvte, swvorttecurvne, swvortnecurvti, swvortticurvne real(GP) :: swvortdivjpar, swvortviscion, swvortdissperp, swvortdissparallel real(GP) :: swvortsource, srcvort, swvortflutterdivjpar, swvortflutterparadv real(GP) :: swfluttervisc, switempflutterpardiffgrad, swneoclvisc real(GP) :: aspectratioinv, freestreaminglimiterqfac, freestreamingfractioni real(GP) :: facheatfluxbraginskiilim, facheatfluxlandau include 'mms_source_slab_vort_head.txt' etai0 = eta_i0 heatfac = heat_fac chipar0i = chipar0_i nuperpvort = perpdiss_coeff_vort dissparallelvort = pardiss_coeff_vort buffervort = buffer_coeff_vort srcvort = src_vort swvortdia = swb_vort_dia swvortexb = swb_vort_exb swvortparadv = swb_vort_paradv swvortnecurvte = swb_vort_necurvte swvorttecurvne = swb_vort_tecurvne swvortnecurvti = swb_vort_necurvti swvortticurvne = swb_vort_ticurvne swvortdivjpar = swb_vort_divjpar swvortviscion = swb_vort_viscion swvortdissperp = swb_vort_dissperp swvortdissparallel = swb_vort_dissparallel swvortsource = swb_vort_source swvortflutterdivjpar = swb_vort_flutter_divjpar swvortflutterparadv = swb_vort_flutter_paradv switempflutterpardiffgrad = swb_itemp_flutter_pardiff_grad swfluttervisc = swb_flutter_visc swneoclvisc = swb_neocl_visc aspectratioinv = aspect_ratio_inv freestreaminglimiterqfac = free_streaming_limiter_qfac freestreamingfractioni = free_streaming_fraction_i if (heatflux_model == 'LANDAU') then facheatfluxbraginskiilim = 0.0_GP facheatfluxlandau = 1.0_GP else facheatfluxbraginskiilim = 1.0_GP facheatfluxlandau = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_vort.txt' mms_source_braginskii_slab_vorticity = res end function real(GP) function mms_source_braginskii_slab_vorticitybsq(equi, x, y, z, t, src_vort) !! MMS source for continuity equation 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) :: src_vort !! Vorticity source real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: etai0, heatfac, chipar0i, nuperpvort, dissparallelvort, buffervort real(GP) :: swvortdia, swvortexb, swvortparadv real(GP) :: swvortnecurvte, swvorttecurvne, swvortnecurvti, swvortticurvne real(GP) :: swvortdivjpar, swvortviscion, swvortdissperp, swvortdissparallel real(GP) :: swvortsource, srcvort, swvortflutterdivjpar, swvortflutterparadv real(GP) :: swfluttervisc, swneoclvisc, aspectratioinv, freestreaminglimiterqfac, freestreamingfractioni, & switempflutterpardiffgrad real(GP) :: facheatfluxbraginskiilim, facheatfluxlandau include 'mms_source_slab_vortbsq_head.txt' etai0 = eta_i0 heatfac = heat_fac chipar0i = chipar0_i nuperpvort = perpdiss_coeff_vort dissparallelvort = pardiss_coeff_vort buffervort = buffer_coeff_vort srcvort = src_vort swvortdia = swb_vort_dia swvortexb = swb_vort_exb swvortparadv = swb_vort_paradv swvortnecurvte = swb_vort_necurvte swvorttecurvne = swb_vort_tecurvne swvortnecurvti = swb_vort_necurvti swvortticurvne = swb_vort_ticurvne swvortdivjpar = swb_vort_divjpar swvortviscion = swb_vort_viscion swvortdissperp = swb_vort_dissperp swvortdissparallel = swb_vort_dissparallel swvortsource = swb_vort_source swvortflutterdivjpar = swb_vort_flutter_divjpar swvortflutterparadv = swb_vort_flutter_paradv switempflutterpardiffgrad = swb_itemp_flutter_pardiff_grad swfluttervisc = swb_flutter_visc swneoclvisc = swb_neocl_visc aspectratioinv = aspect_ratio_inv freestreaminglimiterqfac = free_streaming_limiter_qfac freestreamingfractioni = free_streaming_fraction_i if (heatflux_model == 'LANDAU') then facheatfluxbraginskiilim = 0.0_GP facheatfluxlandau = 1.0_GP else facheatfluxbraginskiilim = 1.0_GP facheatfluxlandau = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_vortbsq.txt' mms_source_braginskii_slab_vorticitybsq = res end function real(GP) function mms_source_braginskii_slab_parmomentum(equi, x, y, z, t, src_upar) !! MMS source for parallel momentum equation 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) :: src_upar !! Parallel momentum source real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: etai0, heatfac, chipar0i, nuperpupar, bufferupar, swuparexb, swuparparadv, swuparcurv real(GP) :: swupargradparn, swupargradparTe, swupargradparTi, swuparviscion, swupardissperp real(GP) :: swuparsource, srcupar real(GP) :: swneoclvisc, aspectratioinv, freestreaminglimiterqfac, freestreamingfractioni real(GP) :: swuparflutterparadv, swuparfluttergradne, swuparfluttergradte, swuparfluttergradti, & swfluttervisc, swuparfluttervisciongrad, switempflutterpardiffgrad real(GP) :: facheatfluxbraginskiilim, facheatfluxlandau include 'mms_source_slab_parmomentum_head.txt' etai0 = eta_i0 heatfac = heat_fac chipar0i = chipar0_i nuperpupar = perpdiss_coeff_upar bufferupar = buffer_coeff_upar srcupar = src_upar swuparexb = swb_upar_exb swuparparadv = swb_upar_paradv swuparcurv = swb_upar_curv swupargradparn = swb_upar_gradpar_n swupargradparTe = swb_upar_gradparTe swupargradparTi = swb_upar_gradparTi swupardissperp = swb_upar_dissperp swuparviscion = swb_upar_viscion swuparsource = swb_upar_source swuparflutterparadv = swb_upar_flutter_paradv swuparfluttergradne = swb_upar_flutter_gradne swuparfluttergradte = swb_upar_flutter_gradte swuparfluttergradti = swb_upar_flutter_gradti swuparfluttervisciongrad = swb_upar_flutter_viscion_grad switempflutterpardiffgrad = swb_itemp_flutter_pardiff_grad swfluttervisc = swb_flutter_visc swneoclvisc = swb_neocl_visc aspectratioinv = aspect_ratio_inv freestreaminglimiterqfac = free_streaming_limiter_qfac freestreamingfractioni = free_streaming_fraction_i if (heatflux_model == 'LANDAU') then facheatfluxbraginskiilim = 0.0_GP facheatfluxlandau = 1.0_GP else facheatfluxbraginskiilim = 1.0_GP facheatfluxlandau = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_parmomentum.txt' mms_source_braginskii_slab_parmomentum = res end function real(GP) function mms_source_braginskii_slab_ohm(equi, x, y, z, t) !! MMS source for parallel momentum equation 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) :: xmin, xmax, ymin, ymax, res real(GP) :: nuperpohm, bufferohm, mu real(GP) :: swohmexb, swohmparadv, swohmphysdiss, swohmgradparne, swohmgradparte, swohmgradparpot, swohmdissperp, & swohmflutterparadv, swohmfluttergradne, swohmfluttergradte, swohmfluttergradpot real(GP) :: thermalforcecoeff include 'mms_source_slab_ohm_head.txt' thermalforcecoeff = thermal_force_coeff nuperpohm = perpdiss_coeff_ohm bufferohm = buffer_coeff_ohm mu = mass_ratio_ei swohmexb = swb_ohm_exb swohmparadv = swb_ohm_paradv swohmphysdiss = swb_ohm_physdiss swohmgradparne = swb_ohm_gradpar_ne swohmgradparte = swb_ohm_gradpar_te swohmgradparpot = swb_ohm_gradpar_pot swohmdissperp = swb_ohm_dissperp swohmflutterparadv = swb_ohm_flutter_paradv swohmfluttergradne = swb_ohm_flutter_gradne swohmfluttergradte = swb_ohm_flutter_gradte swohmfluttergradpot= swb_ohm_flutter_gradpot xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_ohm.txt' mms_source_braginskii_slab_ohm = res end function real(GP) function mms_source_braginskii_slab_etemp(equi, x, y, z, t, src_te) !! MMS source for electron temperature equation 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) :: src_te !! Electron temperature source value real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: PerpdissCoeffNe, PerpdissCoeffTe, PerpdissCoeffPe, bufferNe, bufferTe, bufferPe, thermalforcecoeff, mu real(GP) :: srcte, chipar0e, nue0, freestreamingfractione, freestreaminglimiterqfac real(GP) :: swetempexb, swetempparadv , & swetempcurvpot, swetempcurvte, swetempcurvne, & swetempdissperp, swetempdivvpar, swetempdivjpar, & swetempflutterparadv, swetempflutterdivvpar, swetempflutterdivjpar, & swetempresist, swetempequipart, swetemppardiff, swetempsource, & swetempflutterpardiffdiv, swetempflutterpardiffgrad real(GP) :: facheatfluxbraginskiilim, facheatfluxlandau include 'mms_source_slab_etemp_head.txt' thermalforcecoeff = thermal_force_coeff PerpdissCoeffNe = perpdiss_coeff_ne PerpdissCoeffTe = perpdiss_coeff_te PerpdissCoeffPe = perpdiss_coeff_pe bufferNe = buffer_coeff_ne bufferTe = buffer_coeff_te bufferPe = buffer_coeff_pe chipar0e = chipar0_e nue0 = nu_e0 mu = mass_ratio_ei freestreamingfractione = free_streaming_fraction_e freestreaminglimiterqfac = free_streaming_limiter_qfac srcte = src_te swetempexb = swb_etemp_exb swetempparadv = swb_etemp_paradv swetempcurvpot = swb_etemp_curv_pot swetempcurvte = swb_etemp_curv_te swetempcurvne = swb_etemp_curv_ne swetempdivvpar = swb_etemp_divvpar swetempdivjpar = swb_etemp_divjpar swetempresist = swb_etemp_resist swetempequipart = swb_etemp_equipart swetemppardiff = swb_etemp_pardiff swetempdissperp= swb_etemp_dissperp swetempsource = swb_etemp_source swetempflutterparadv = swb_etemp_flutter_paradv swetempflutterdivvpar= swb_etemp_flutter_divvpar swetempflutterdivjpar= swb_etemp_flutter_divjpar swetempflutterpardiffdiv = swb_etemp_flutter_pardiff_div swetempflutterpardiffgrad= swb_etemp_flutter_pardiff_grad if (heatflux_model == 'LANDAU') then facheatfluxbraginskiilim = 0.0_GP facheatfluxlandau = 1.0_GP else facheatfluxbraginskiilim = 1.0_GP facheatfluxlandau = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_etemp.txt' mms_source_braginskii_slab_etemp = res end function real(GP) function mms_source_braginskii_slab_itemp(equi, x, y, z, t, src_ti) !! MMS source for ion temperature equation 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) :: src_ti !! Ion temperature source value real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: PerpdissCoeffNe, PerpdissCoeffTi, PerpdissCoeffPi, bufferNe, bufferTi, bufferPi, & etai0, heatfac, chipar0i, nue0, & freestreamingfractioni, freestreaminglimiterqfac, mu real(GP) :: srcti real(GP) :: switempexb, & switempcurvpot, switempcurvte, switempcurvti, switempcurvne, & switempdissperp, switempparadv, switemppardiff, switempdivupar, & switempflutterparadv, switempflutterdivupar, switempflutterdivjpar, & switempdivjpar, switempequipart, switempvischeat, switempsource, & switempflutterpardiffdiv, switempflutterpardiffgrad real(GP) :: swfluttervisc, swneoclvisc, AspectRatioInv real(GP) :: facheatfluxbraginskiilim, facheatfluxlandau include 'mms_source_slab_itemp_head.txt' PerpdissCoeffNe = perpdiss_coeff_ne PerpdissCoeffTi = perpdiss_coeff_ti PerpdissCoeffPi = perpdiss_coeff_pi bufferNe = buffer_coeff_ne bufferTi = buffer_coeff_ti bufferPi = buffer_coeff_pi etai0 = eta_i0 heatfac = heat_fac chipar0i = chipar0_i nue0 = nu_e0 mu = mass_ratio_ei freestreamingfractioni = free_streaming_fraction_i freestreaminglimiterqfac = free_streaming_limiter_qfac srcti = src_ti switempexb = swb_itemp_exb switempcurvpot = swb_itemp_curv_pot switempcurvte = swb_itemp_curv_te switempcurvti = swb_itemp_curv_ti switempcurvne = swb_itemp_curv_ne switempdissperp = swb_itemp_dissperp switempparadv = swb_itemp_paradv switemppardiff = swb_itemp_pardiff switempdivupar = swb_itemp_divupar switempdivjpar = swb_itemp_divjpar switempequipart = swb_itemp_equipart switempvischeat = swb_itemp_vischeat switempsource = swb_itemp_source switempflutterparadv = swb_itemp_flutter_paradv switempflutterdivupar= swb_itemp_flutter_divupar switempflutterdivjpar= swb_itemp_flutter_divjpar switempflutterpardiffdiv = swb_itemp_flutter_pardiff_div switempflutterpardiffgrad= swb_itemp_flutter_pardiff_grad swfluttervisc = swb_flutter_visc swneoclvisc = swb_neocl_visc AspectRatioInv = aspect_ratio_inv if (heatflux_model == 'LANDAU') then facheatfluxbraginskiilim = 0.0_GP facheatfluxlandau = 1.0_GP else facheatfluxbraginskiilim = 1.0_GP facheatfluxlandau = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_itemp.txt' mms_source_braginskii_slab_itemp = res end function real(GP) function mms_source_braginskii_slab_landau_e(equi, x, y, z, t, landaualpha, landaubeta) !! MMS source for elliptic Landau heat flux equation 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) :: landaualpha !! Numerical parameter alpha of current Lorentzian real(GP), intent(in) :: landaubeta !! Numerical parameter beta of current Lorentzian real(GP), parameter :: sqrteightbypi = sqrt(8.0_GP / PI) real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: mu, chipar0e real(GP) :: swetempflutterpardiffgrad real(GP) :: faclandauflutterlhs include 'mms_source_slab_landau_e_head.txt' mu = mass_ratio_ei chipar0e = chipar0_e swetempflutterpardiffgrad = swb_etemp_flutter_pardiff_grad if(landau_flutter_lhs_on) then faclandauflutterlhs = 1.0_GP else faclandauflutterlhs = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_landau_e.txt' mms_source_braginskii_slab_landau_e = res end function real(GP) function mms_source_braginskii_slab_landau_i(equi, x, y, z, t, landaualpha, landaubeta) !! MMS source for elliptic Landau heat flux equation 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) :: landaualpha !! Numerical parameter alpha of current Lorentzian real(GP), intent(in) :: landaubeta !! Numerical parameter beta of current Lorentzian real(GP), parameter :: sqrteightbypi = sqrt(8.0_GP / PI) real(GP) :: xmin, xmax, ymin, ymax, res real(GP) :: chipar0i real(GP) :: switempflutterpardiffgrad real(GP) :: faclandauflutterlhs include 'mms_source_slab_landau_i_head.txt' chipar0i = chipar0_i switempflutterpardiffgrad = swb_itemp_flutter_pardiff_grad if(landau_flutter_lhs_on) then faclandauflutterlhs = 1.0_GP else faclandauflutterlhs = 0.0_GP endif xmin = equi%xmin xmax = equi%xmax ymin = equi%ymin ymax = equi%ymax include 'mms_source_slab_landau_i.txt' mms_source_braginskii_slab_landau_i = res end function end module