module mms_braginskii_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 ! 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 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_boundaries_parpen_m, only : & pen_epsinv use params_brag_switches_m implicit none public :: mms_sol_braginskii_circular_ne public :: mms_sol_braginskii_circular_pot public :: mms_sol_braginskii_circular_vort public :: mms_sol_braginskii_circular_upar public :: mms_sol_braginskii_circular_jpar public :: mms_sol_braginskii_circular_apar public :: mms_sol_braginskii_circular_te public :: mms_sol_braginskii_circular_ti public :: mms_source_braginskii_circular_continuity public :: mms_source_braginskii_circular_vorticity public :: mms_source_braginskii_circular_parmomentum public :: mms_source_braginskii_circular_ohm public :: mms_source_braginskii_circular_etemp public :: mms_source_braginskii_circular_itemp private :: cart_to_circ 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 !! 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.2E-1_GP !! Radial mode number of MMS Solutions real(GP), private, parameter :: krnne = 1.0_GP real(GP), private, parameter :: krnpot = 1.0_GP real(GP), private, parameter :: krnupar = 2.0_GP real(GP), private, parameter :: krnapar = 1.0_GP real(GP), private, parameter :: krnte = 2.0_GP real(GP), private, parameter :: krnti = 2.0_GP !! Poloidal mode number of MMS Solutions real(GP), private, parameter :: kpne = 2.0_GP real(GP), private, parameter :: kppot = 3.0_GP real(GP), private, parameter :: kpupar = 2.0_GP real(GP), private, parameter :: kpapar = 3.0_GP real(GP), private, parameter :: kpte = 3.0_GP real(GP), private, parameter :: kpti = 2.0_GP !! Poloidal phase shift of MMS Solutions real(GP), private, parameter :: phpne = 0.1_GP real(GP), private, parameter :: phppot = 0.7_GP real(GP), private, parameter :: phpupar = 1.2_GP real(GP), private, parameter :: phpapar = 0.5_GP real(GP), private, parameter :: phpte = 0.3_GP real(GP), private, parameter :: phpti = 0.4_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 !! 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 !! 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 !! 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 contains real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_ne_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_ne.txt' mms_sol_braginskii_circular_ne = res end function real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_pot_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_pot.txt' mms_sol_braginskii_circular_pot = res end function real(GP) function mms_sol_braginskii_circular_vort(equi, x, y, z, t) !! MMS solution for generalised vorticity !! TODO: Boussinesq version 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) :: r, p, q, rmin, rmax, res real(GP) :: swvortdia include 'mms_sol_circular_vort_head.txt' swvortdia = swb_vort_dia ! 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_vort.txt' mms_sol_braginskii_circular_vort = res end function real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_upar_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_upar.txt' mms_sol_braginskii_circular_upar = res end function real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_jpar_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_jpar.txt' mms_sol_braginskii_circular_jpar = res end function real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_apar_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_apar.txt' mms_sol_braginskii_circular_apar = res end function real(GP) function mms_sol_braginskii_circular_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) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_te_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_te.txt' mms_sol_braginskii_circular_te = res end function real(GP) function mms_sol_braginskii_circular_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 coordinate) real(GP), intent(in) :: t !! Time real(GP) :: r, p, q, rmin, rmax, res include 'mms_sol_circular_ti_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_ti.txt' mms_sol_braginskii_circular_ti = res end function real(GP) function mms_source_braginskii_circular_continuity(equi, x, y, z, t, chi, lognepen, 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: lognepen !! Value that shall be penalised to real(GP), intent(in) :: src_ne !! Particle source value real(GP) :: r, p, q, rmin, rmax, res real(GP) :: PerpdissCoeffNe, dissparallelne, bufferNe, epsinv real(GP) :: srcne real(GP) :: swcontexb, & swcontcurvpot, swcontcurvte, swcontcurvne, & swcontdivjpar, swcontdivparflx, & swcontdissperp, swcontdissparallel, & swcontsource, & swcontflutterparadv , swcontflutterdivjpar, swcontflutterdivupar include 'mms_source_circular_continuity_head.txt' PerpdissCoeffNe = perpdiss_coeff_ne dissparallelne = pardiss_coeff_ne bufferNe = buffer_coeff_ne epsinv = pen_epsinv 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 ! 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_continuity.txt' mms_source_braginskii_circular_continuity = res end function real(GP) function mms_source_braginskii_circular_vorticity(equi, x, y, z, t, chi, vortpen, 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: vortpen !! Value that shall be penalised to real(GP), intent(in) :: src_vort !! Vorticity source value real(GP) :: r, p, q, rmin, rmax, res real(GP) :: etai0, heatfac, nuperpvort, dissparallelvort, buffervort, epsinv, chipar0i 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, switempflutterpardiffgrad real(GP) :: aspectratioinv, freestreamingfractioni, freestreaminglimiterqfac include 'mms_source_circular_vort_head.txt' etai0 = eta_i0 heatfac = heat_fac nuperpvort = perpdiss_coeff_vort dissparallelvort = pardiss_coeff_vort buffervort = buffer_coeff_vort epsinv = pen_epsinv chipar0i = chipar0_i 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 ! 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_vort.txt' mms_source_braginskii_circular_vorticity = res end function real(GP) function mms_source_braginskii_circular_parmomentum(equi, x, y, z, t, chi, uparpen, 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: uparpen !! Value that shall be penalised to real(GP), intent(in) :: src_upar !! Parallel momentum source value real(GP) :: r, p, q, rmin, rmax, res real(GP) :: etai0, heatfac, chipar0i, nuperpupar, bufferupar, epsinv, 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 include 'mms_source_circular_parmomentum_head.txt' etai0 = eta_i0 heatfac = heat_fac chipar0i = chipar0_i nuperpupar = perpdiss_coeff_upar bufferupar = buffer_coeff_upar epsinv = pen_epsinv 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 swuparviscion = swb_upar_viscion swupardissperp = swb_upar_dissperp 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 ! 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_parmomentum.txt' mms_source_braginskii_circular_parmomentum = res end function real(GP) function mms_source_braginskii_circular_ohm(equi, x, y, z, t, chi, psiparpen) !! 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: psiparpen !! Value that shall be penalised to real(GP) :: r, p, q, rmin, rmax, res real(GP) :: nuperpohm, bufferohm, mu, epsinv real(GP) :: swohmexb, swohmparadv, swohmphysdiss, swohmgradparne, & swohmgradparte, swohmgradparpot, swohmdissperp, thermalforcecoeff, & swohmflutterparadv, swohmfluttergradne, swohmfluttergradte, swohmfluttergradpot include 'mms_source_circular_ohm_head.txt' thermalforcecoeff = thermal_force_coeff nuperpohm = perpdiss_coeff_ohm bufferohm = buffer_coeff_ohm mu = mass_ratio_ei epsinv = pen_epsinv 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 ! 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_ohm.txt' mms_source_braginskii_circular_ohm = res end function real(GP) function mms_source_braginskii_circular_etemp(equi, x, y, z, t, chi, logtepen, 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: logtepen !! Value that shall be penalised to real(GP), intent(in) :: src_te !! Electron temperature source value real(GP) :: r, p, q, rmin, rmax, res, mu, epsinv real(GP) :: PerpdissCoeffNe, PerpdissCoeffTe, PerpdissCoeffPe, bufferNe, bufferTe, bufferPe, & thermalforcecoeff, chipar0e, nue0, freestreamingfractione, freestreaminglimiterqfac real(GP) :: srcte real(GP) :: swetempexb, swetempparadv , & swetempcurvpot, swetempcurvte, swetempcurvne, & swetempdissperp, swetempdivvpar, swetempdivjpar, & swetempresist, swetempequipart, swetemppardiff, swetempsource, & swetempflutterparadv, swetempflutterdivvpar, swetempflutterdivjpar, & swetempflutterpardiffdiv, swetempflutterpardiffgrad include 'mms_source_circular_etemp_head.txt' thermalforcecoeff = thermal_force_coeff PerpdissCoeffNe = perpdiss_coeff_ne PerpdissCoeffTe = perpdiss_coeff_te PerpdissCoeffPe = perpdiss_coeff_pe chipar0e = chipar0_e nue0 = nu_e0 mu = mass_ratio_ei bufferNe = buffer_coeff_ne bufferTe = buffer_coeff_te bufferPe = buffer_coeff_pe freestreamingfractione = free_streaming_fraction_e freestreaminglimiterqfac = free_streaming_limiter_qfac epsinv = pen_epsinv 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 ! 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_etemp.txt' mms_source_braginskii_circular_etemp = res end function real(GP) function mms_source_braginskii_circular_itemp(equi, x, y, z, t, chi, logtipen, 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) :: chi !! Value of characteristc function of penalisation real(GP), intent(in) :: logtipen !! Value that shall be penalised to real(GP), intent(in) :: src_ti !! Ion temperature source value real(GP) :: r, p, q, rmin, rmax, res real(GP) :: PerpdissCoeffNe, PerpdissCoeffTi, PerpdissCoeffPi, bufferNe, bufferTi, bufferPi, & etai0, heatfac, chipar0i, nue0, mu, epsinv, freestreamingfractioni, freestreaminglimiterqfac real(GP) :: srcti real(GP) :: switempexb, & switempcurvpot, switempcurvte, switempcurvti, switempcurvne, & switempdissperp, switempparadv, switemppardiff, switempdivupar, & switempdivjpar, switempequipart, switempvischeat, switempsource, & switempflutterparadv, switempflutterdivupar, switempflutterdivjpar, & switempflutterpardiffdiv, switempflutterpardiffgrad real(GP) :: swfluttervisc, swneoclvisc, aspectratioinv include 'mms_source_circular_itemp_head.txt' etai0 = eta_i0 heatfac = heat_fac PerpdissCoeffNe = perpdiss_coeff_ne PerpdissCoeffTi = perpdiss_coeff_ti PerpdissCoeffPi = perpdiss_coeff_pi bufferNe = buffer_coeff_ne bufferTi = buffer_coeff_ti bufferPi = buffer_coeff_pi chipar0i = chipar0_i nue0 = nu_e0 mu = mass_ratio_ei freestreamingfractioni = free_streaming_fraction_i freestreaminglimiterqfac = free_streaming_limiter_qfac epsinv = pen_epsinv 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 ! 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_itemp.txt' mms_source_braginskii_circular_itemp = 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 ! Circular geometry is axisymmetric 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