!> FIXME : Add documentation module species use abstract_config, only: abstract_config_type, CONFIG_MAX_NAME_LEN implicit none private public :: init_species, finish_species, wnml_species, check_species public :: nspec, spec, tracer_species public :: has_ion_species, has_electron_species public :: has_hybrid_electron_species, is_electron_species, is_hybrid_electron_species public :: calculate_f0_arrays, set_current_f0_species, eval_f0 public :: nonmaxw_corr, f0_values, dlnf0drho, f0_maxwellian, f0_sdanalytic, set_overrides public :: species_config_type, species_element_config_type public :: set_species_config, get_species_config, get_species_element_config !> FIXME : Add documentation type :: specie !> Main physical properties of the species, see namelist documentation. real :: z, mass, dens, temp, tprim, fprim, vnewk !> Common combinations of species properties real :: stm, zstm, tz, smz, zt !> Provides offsets to parallel flow shear drive term real :: uprim, uprim2 !> Parameters associated with a few initialisation options, including Orszag-Tang. real :: dens0, u0, tpar0, tperp0 !> nu_h controls a hyperviscous term embedded in the collision operator real :: nu_h !> Artificial factor multiplying the Bessel function argument real :: bess_fac !> Parameters associated with analytic slowing-down distribution real :: vcrit, vcprim !> Flags to identify the species type and the background distribution treatment integer :: type, f0type end type specie integer, parameter :: ion_species = 1 integer, parameter :: electron_species = 2 ! for collision operator integer, parameter :: tracer_species = 3 ! for test particle diffusion studies integer, parameter :: hybrid_electron_species = 4 ! Adibatic passing, non-adiabatic trapped integer, dimension(*), parameter :: electron_like_species = & [electron_species, hybrid_electron_species] integer, parameter :: & !> Regular Maxwellian species f0_maxwellian = 1, & !> f0 determined from table in external input file f0_tabulated = 2, & !> Analytic Gaffey slowing-down distribution f0_sdanalytic = 3 !> Various quantities needed in the calculation of the !> slowing-down distribution in case electrons are adiabatic real:: me, ZI_fac !> Defines a working species so that eval_f0 only needs to take one !> argument (needed for genquad) integer:: is_global = -1 integer :: nspec type (specie), dimension (:), allocatable :: spec real, dimension (:,:), allocatable :: nonmaxw_corr, f0_values, dlnf0drho logical :: initialized = .false. !> Used to represent the input configuration of species type, extends(abstract_config_type) :: species_config_type ! namelist : species_knobs ! indexed : false !> Specifies the electron mass used in calculating the slowing !> down parameters in simulations without kinetic electrons. real :: me = -1.0 !> Number of kinetic species to evolve. integer :: nspec = 2 !> Specify the ion charge factor used in calculating the slowing !> down parameters. If not specified that will attempt to !> calculate from the kinetic ion species. real :: zi_fac = -1.0 contains procedure, public :: read => read_species_config procedure, public :: write => write_species_config procedure, public :: reset => reset_species_config procedure, public :: broadcast => broadcast_species_config procedure, public, nopass :: get_default_name => get_default_name_species_config procedure, public, nopass :: get_default_requires_index => get_default_requires_index_species_config end type species_config_type type(species_config_type) :: species_config !> Used to represent the input configuration of species type, extends(abstract_config_type) :: species_element_config_type ! namelist : species_parameters ! indexed : true !> Multiplies the argument of the Bessel function for this !> species. Useful for exploring the effect of the Bessel !> functions. real :: bess_fac = 1.0 !> Set the normalised density for this species real :: dens = 1.0 !> Parameter used for a few specific initial condition types. !> !> @todo Consider moving this parameter to [[init_g_knobs]]. real :: dens0 = 1.0 !> Select which type of background distribution should be assumed !> for this species. Must be one of: !> !> - 'default' : equivalent to maxwellian !> - 'maxwellian' : Standard maxwellian !> - 'tabulated' : A general tabulated form !> - 'sdanalytic' : A analytic slowing down form. !> !> See [G. Wilkie !> thesis](https://drum.lib.umd.edu/handle/1903/17302). character(len = 20) :: f0type = "default" !> Normalised inverse density gradient: \(-\frac{1}{n}\frac{dn}{d !> \rho_N}\) (Note here \(\rho_N\) is the normalised radius !> \(\frac{\rho}{L_\textrm{ref}}\). real :: fprim = 2.2 !> Normalised mass of this species real :: mass = 1.0 !> Set the species hyper collision frequency used in the pitch !> angle scattering collision operator if !> [[collision_knobs::hypermult]] is `true`. real :: nu_h = 0.0 !> Set the normalised species temperature real :: temp = 1.0 !> Parameter used for a few specific initial condition types. !> !> @todo Consider moving this parameter to [[init_g_knobs]]. real :: tpar0 = 0. !> Parameter used for a few specific initial condition types. !> !> @todo Consider moving this parameter to [[init_g_knobs]]. real :: tperp0 = 0. !> Normalised inverse temperature gradient: \(-\frac{1}{T}\frac{dT}{d !> \rho_N}\) (Note here \(\rho_N\) is the normalised radius !> \(\frac{\rho}{L_\textrm{ref}}\). real :: tprim = 6.9 !> Sets the characterisation of the species. Should be one of !> !> - 'ion' : Thermal ion species !> - 'default' : Same as 'ion' !> - 'electron' : Thermal electron species !> - 'e' : Same as 'electron' !> - 'trace' : Tracer species !> - 'hybrid_electron' : Adiabatic passing electrons, full !> trapped electron treatment. *NOTE: This option is currently !> considered experimental and is intended for use in !> electrostatic simulations.* !> !> This primarily controls the treatment in the collision !> operator, the detection of adiabatic species. Tracer species !> do not contribute to velocity space integrals. character(len = 20) :: type = "default" !> Parameter used for a few specific initial condition types. !> !> @todo Consider moving this parameter to [[init_g_knobs]]. real :: u0 = 1.0 !> Controls an energy independent part of a source term. !> !> @note Documentation needs improving, looks like it's an early !> attempt at adding in the parallel flow shear drive term. This !> effectively provides a constant offset to the parallel flow !> shear term, allowing one to add a species dependent shift. real :: uprim = 0.0 !> Controls an energy dependent part of a source term. !> !> @note Documentation needs improving, looks like it's an early !> attempt at adding in the parallel flow shear drive term. This !> effectively provides a energy dependent offset to the parallel !> flow shear term, allowing one to add a species and energy !> dependent shift. real :: uprim2 = 0.0 !> Part of the `sdanalytic` slowing down model. See [G. Wilkie !> thesis](https://drum.lib.umd.edu/handle/1903/17302). real :: vcprim = -999.9 !> Part of the `sdanalytic` slowing down model. See [G. Wilkie !> thesis](https://drum.lib.umd.edu/handle/1903/17302). real :: vcrit = -1.0 !> Sets the normalised species-species collisionality !> frequency. For species s, `vnewk = nu_ss Lref/vref` where Lref !> is the reference length, `vref` is the reference thermal speed !> (not the thermal speed of species s!), \(\nu_{ss}=\sqrt{2}\pi !> n_s Z_s^4 e^4 \log(\lambda)/\sqrt{m_s}T^{3/2}_s\) is a !> dimensional collision frequency, e is the proton charge, !> \(\log(\lambda)\) is the Coloumb logarithm, and !> \(Z_s,T_s,n_s,m_s\) are the (charge, temperature, density, !> mass) of species s. (The dimensional collision frequency given !> here is in Gaussian units. For SI units, include an extra !> factor \((4\pi\epsilon_0)^2\) in the denominator of the !> definition of \(\nu_{ss}\). real :: vnewk = 0.0 !> Normalised species charge real :: z = 1.0 contains procedure, public :: read => read_species_element_config procedure, public :: write => write_species_element_config procedure, public :: reset => reset_species_element_config procedure, public :: broadcast => broadcast_species_element_config procedure, public, nopass :: get_default_name => get_default_name_species_element_config procedure, public, nopass :: get_default_requires_index => get_default_requires_index_species_element_config end type species_element_config_type type(species_element_config_type), dimension(:), allocatable :: species_element_config contains !> FIXME : Add documentation subroutine check_species(report_unit,beta,tite,alne,dbetadrho_spec) implicit none integer, intent(in) :: report_unit real, intent(in) :: beta, tite real, intent(out) :: alne, dbetadrho_spec integer :: is real :: aln, alp, charge, ee, ne, ptot, zeff_calc write (report_unit, fmt="('Number of species: ',i3)") nspec zeff_calc = 0. charge = 0. aln = 0. alne = 0. alp = 0. ptot = 0. do is=1, nspec write (report_unit, *) write (report_unit, fmt="(' Species ',i3)") is if (spec(is)%type == ion_species) write (report_unit, fmt="(' Type: Ion')") if (spec(is)%type == electron_species) write (report_unit, fmt="(' Type: Electron')") if (spec(is)%type == tracer_species) write (report_unit, fmt="(' Type: Trace')") if (spec(is)%type == hybrid_electron_species) write (report_unit, fmt="(' Type: Hybrid-Electron')") write (report_unit, fmt="(' Charge: ',f7.3)") spec(is)%z write (report_unit, fmt="(' Mass: ',es11.4)") spec(is)%mass write (report_unit, fmt="(' Density: ',f7.3)") spec(is)%dens write (report_unit, fmt="(' Temperature: ',f7.3)") spec(is)%temp write (report_unit, fmt="(' Collisionality: ',es11.4)") spec(is)%vnewk write (report_unit, fmt="(' Normalized Inverse Gradient Scale Lengths:')") write (report_unit, fmt="(' Temperature: ',f7.3)") spec(is)%tprim write (report_unit, fmt="(' Density: ',f7.3)") spec(is)%fprim write (report_unit, fmt="(' Parallel v: ',f7.3)") spec(is)%uprim if (spec(is)%bess_fac.ne.1.0) & write (report_unit, fmt="(' Bessel function arg multiplied by: ',f7.3)") spec(is)%bess_fac ! write (report_unit, fmt="(' Ignore this:')") ! write (report_unit, fmt="(' D_0: ',es10.4)") spec(is)%dens0 if (spec(is)%type /= 2) then zeff_calc = zeff_calc + spec(is)%dens*spec(is)%z**2 charge = charge + spec(is)%dens*spec(is)%z aln = aln + spec(is)%dens*spec(is)%z*spec(is)%fprim else alne = alne + spec(is)%dens*spec(is)%z*spec(is)%fprim ne = spec(is)%dens ee = spec(is)%z end if alp = alp + spec(is)%dens * spec(is)%temp *(spec(is)%fprim + spec(is)%tprim) ptot = ptot + spec(is)%dens * spec(is)%temp end do if (.not. has_electron_species(spec)) then ptot = ptot + 1./tite ! electron contribution to pressure alp = alp + aln/tite ! assuming charge neutrality, electron contribution to alp end if alp = alp / ptot write (report_unit, *) write (report_unit, fmt="('------------------------------------------------------------')") write (report_unit, fmt="('Calculated Z_eff: ',f7.3)") zeff_calc if (has_electron_species(spec)) then if (abs(charge+ne*ee) > 1.e-2) then if (charge+ne*ee < 0.) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('You are neglecting an ion species.')") write (report_unit, fmt="('This species has a charge fraction of ',f7.3)") abs(charge+ne*ee) write (report_unit, & & fmt="('and a normalized inverse density gradient scale length of ',f7.3)") & (aln+alne)/(charge+ne*ee) write (report_unit, fmt="('################# WARNING #######################')") else write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('There is an excess ion charge fraction of ',f7.3)") abs(charge+ne*ee) write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") end if else if (abs(aln+alne) > 1.e-2) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('The density gradients are inconsistent'/' a/lni =',e12.4,' but alne =',e12.4)") aln, alne write (report_unit, fmt="('################# WARNING #######################')") end if end if else if (charge > 1.01) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('There is an excess ion charge fraction of ',f7.3)") charge-1. write (report_unit, fmt="('THIS IS PROBABLY AN ERROR.')") write (report_unit, fmt="('################# WARNING #######################')") end if if (charge < 0.99) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('You are neglecting an ion species.')") write (report_unit, fmt="('This species has a charge fraction of ',f7.3)") abs(charge-1.) write (report_unit, fmt="('################# WARNING #######################')") end if end if if (has_hybrid_electron_species(spec)) then write (report_unit, *) write (report_unit, fmt="('################# WARNING #######################')") write (report_unit, fmt="('You are using a hybrid electron species.')") write (report_unit, fmt="('This is currently considered experimental.')") write (report_unit, fmt="('################# WARNING #######################')") end if write (report_unit, *) write (report_unit, fmt="('------------------------------------------------------------')") write (report_unit, *) write (report_unit, fmt="('GS2 beta parameter = ',f9.4)") beta write (report_unit, fmt="('Total beta = ',f9.4)") beta*ptot write (report_unit, *) write (report_unit, fmt="('The total normalized inverse pressure gradient scale length is ',f10.4)") alp dbetadrho_spec = -beta*ptot*alp write (report_unit, fmt="('corresponding to d beta / d rho = ',f10.4)") dbetadrho_spec end subroutine check_species !> FIXME : Add documentation subroutine wnml_species(unit) implicit none integer, intent(in) :: unit integer :: i character (100) :: line write (unit, *) write (unit, fmt="(' &',a)") "species_knobs" write (unit, fmt="(' nspec = ',i2)") nspec write (unit, fmt="(' /')") do i=1,nspec write (unit, *) write (line, *) i write (unit, fmt="(' &',a)") & & trim("species_parameters_"//trim(adjustl(line))) write (unit, fmt="(' z = ',e13.6)") spec(i)%z write (unit, fmt="(' mass = ',e13.6)") spec(i)%mass write (unit, fmt="(' dens = ',e13.6)") spec(i)%dens write (unit, fmt="(' temp = ',e13.6)") spec(i)%temp write (unit, fmt="(' tprim = ',e13.6)") spec(i)%tprim write (unit, fmt="(' fprim = ',e13.6)") spec(i)%fprim write (unit, fmt="(' uprim = ',e13.6)") spec(i)%uprim write (unit, fmt="(' vcrit = ',e13.6)") spec(i)%vcrit write (unit, fmt="(' vcprim = ',e13.6)") spec(i)%vcprim if (spec(i)%uprim2 /= 0.) write (unit, fmt="(' uprim2 = ',e13.6)") spec(i)%uprim2 write (unit, fmt="(' vnewk = ',e13.6)") spec(i)%vnewk if (spec(i)%type == ion_species) & write (unit, fmt="(a)") ' type = "ion" /' if (spec(i)%type == electron_species) & write (unit, fmt="(a)") ' type = "electron" /' if (spec(i)%type == hybrid_electron_species) & write (unit, fmt="(a)") ' type = "hybrid_electron" /' if (spec(i)%type == tracer_species) & write (unit, fmt="(a)") ' type = "trace" /' if (spec(i)%f0type == f0_maxwellian) & write (unit, fmt="(a)") ' f0type = "maxwellian" /' if (spec(i)%f0type == f0_sdanalytic) & write (unit, fmt="(a)") ' f0type = "sdanalytic" /' if (spec(i)%f0type == f0_tabulated) & write (unit, fmt="(a)") ' f0type = "tabulated" /' write (unit, fmt="(' dens0 = ',e13.6)") spec(i)%dens0 if(spec(i)%bess_fac.ne.1.0) & write (unit, fmt="(' bess_fac = ',e13.6)") spec(i)%bess_fac end do end subroutine wnml_species !> FIXME : Add documentation subroutine init_species(species_config_in, species_elements_config_in) implicit none type(species_config_type), intent(in), optional :: species_config_in type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_elements_config_in if (initialized) return initialized = .true. call read_parameters(species_config_in, species_elements_config_in) end subroutine init_species !> FIXME : Add documentation subroutine read_parameters(species_config_in, species_elements_config_in) use file_utils, only: error_unit use text_options, only: text_option, get_option_value use mp, only: proc0, mp_abort implicit none type(species_config_type), intent(in), optional :: species_config_in type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_elements_config_in integer :: is, ierr character(len = 20) :: type, f0type type (text_option), dimension (*), parameter :: typeopts = [ & text_option('default', ion_species), & text_option('ion', ion_species), & text_option('electron', electron_species), & text_option('e', electron_species), & text_option('trace', tracer_species), & text_option('hybrid_electron', hybrid_electron_species) & ] type (text_option), dimension (*), parameter :: f0_opts = [ & text_option('default', f0_maxwellian), & text_option('maxwellian', f0_maxwellian), & text_option('tabulated', f0_tabulated), & text_option('sdanalytic', f0_sdanalytic) & ] if (present(species_config_in)) species_config = species_config_in call species_config%init(name = 'species_knobs', requires_index = .false.) ! Copy out internal values into module level parameters associate(self => species_config) #include "species_copy_out_auto_gen.inc" end associate ierr = error_unit() if (proc0 .and. nspec < 1) then write (unit=ierr, fmt="('Invalid nspec in species_knobs: ', i5)") nspec call mp_abort('Invalid nspec in species_knobs') end if allocate (spec(nspec)) if (present(species_elements_config_in)) species_element_config = species_elements_config_in if (.not.allocated(species_element_config)) allocate(species_element_config(nspec)) ! This might look odd given the above line _but_ we might have set species_element_config ! elsewhere and hence not allocated in the above, potentially giving inconsistencies if (size(species_element_config) /= nspec) then if (proc0) print*,"inconsistent number of config elements" end if do is = 1, nspec call species_element_config(is)%init(name = 'species_parameters', requires_index = .true., index = is) ! Copy out internal values into module level parameters associate(self => species_element_config(is), bess_fac => spec(is)%bess_fac, & dens => spec(is)%dens, dens0 => spec(is)%dens0, fprim => spec(is)%fprim, & mass => spec(is)%mass, nu_h => spec(is)%nu_h, temp => spec(is)%temp, & tpar0 => spec(is)%tpar0, tperp0 => spec(is)%tperp0, tprim => spec(is)%tprim, & u0 => spec(is)%u0, uprim => spec(is)%uprim, uprim2 => spec(is)%uprim2, & vcprim => spec(is)%vcprim, vcrit => spec(is)%vcrit, vnewk => spec(is)%vnewk, & z => spec(is)%z) #include "species_element_copy_out_auto_gen.inc" end associate call get_option_value (species_element_config(is)%type, typeopts, spec(is)%type, & ierr, "type in species_parameters_x",.true.) call get_option_value (species_element_config(is)%f0type, f0_opts, spec(is)%f0type, & ierr, "f0type in species_parameters_x",.true.) end do spec%stm = sqrt(spec%temp / spec%mass) spec%zstm = spec%z / sqrt(spec%temp * spec%mass) spec%tz = spec%temp / spec%z spec%zt = spec%z / spec%temp spec%smz = abs(sqrt(spec%temp * spec%mass) / spec%z) end subroutine read_parameters !> FIXME : Add documentation subroutine calculate_slowingdown_parameters(alpha_index, vc, vcprim) use constants, only: pi use mp, only: proc0 implicit none integer, intent(in):: alpha_index real, intent(out):: vc, vcprim integer:: electron_spec, main_ion_species, is, js real:: vinj, Te_prim, ZI_fac_local, ne_prim, ni_prim, ne, vte, charge_contribution electron_spec = -1 main_ion_species = -1 ! Find species index corresponding to electrons and "main ions" ! where we assume the first ion species corresponds to the main ions. ! Perhaps we're looking for the species with mass = 1 (i.e. the reference mass?) do is = 1, nspec if (is_electron_species(spec(is))) electron_spec = is if (main_ion_species < 1 .and. spec(is)%type == ion_species) main_ion_species = is end do if (electron_spec > 0) then ! If we found the electron index extract its parameters me = spec(electron_spec)%mass ne = spec(electron_spec)%dens ne_prim = spec(electron_spec)%fprim Te_prim = spec(electron_spec)%tprim vte = sqrt(spec(electron_spec)%temp/spec(electron_spec)%mass) else ! If we don't have an electron species try to work out the relevant parameters if (me < 0.0) then if (proc0) write(*,*) "WARNING: me not specified in species_knobs. Assuming reference species is deuterium, so me = 0.0002723085" me = 0.0002723085 end if ! Set electron density and gradient from quasineutrality ne = sum(spec%z * spec%dens) ne_prim = sum(spec%z * spec%dens * spec%fprim) / ne ! Assume Te' = Ti' of the main ions Te_prim = spec(main_ion_species)%tprim ! Compute electron thermal velocity using the main ion properties? This seems ! surprising. Whilst Te = Ti is perhaps reasonable, we specify the electron mass ! so why not use it here? Could we also use the tite input to scale the ion T? vte = sqrt(spec(main_ion_species)%temp/spec(main_ion_species)%mass) if (proc0) then write(*,*) "Since electrons are adiabatic, we need to improvise electron parameters to " write(*,*) "caluclate alpha species properties. Imposed by global quasineutrality:" write(*,*) " me = ", me write(*,*) " ne = ", ne write(*,*) " vte = ", vte write(*,*) " ne_prim = ", ne_prim write(*,*) " Te_prim = ", Te_prim end if end if ni_prim = 0.0 ! If ZI_fac is not specified, calculate it from existing ion species: if (ZI_fac < 0.0) then ZI_fac_local = 0.0 do js = 1, nspec if (spec(js)%type == ion_species) then charge_contribution = spec(alpha_index)%mass * spec(js)%dens * spec(js)%z**2 / (spec(js)%mass * ne) ZI_fac_local = ZI_fac_local + charge_contribution ni_prim = ni_prim + charge_contribution * spec(js)%fprim end if end do else ZI_fac_local = ZI_fac end if ni_prim = ni_prim / ZI_fac_local ! Calculate vc/vinj vinj = sqrt(spec(alpha_index)%temp / spec(alpha_index)%mass) vc = (0.25 * 3 * sqrt(pi) * ZI_fac_local * me / spec(alpha_index)%mass)**(1/3.0) * & (vte / vinj) vcprim = ni_prim - ne_prim + 1.5 * Te_prim end subroutine calculate_slowingdown_parameters !> FIXME : Add documentation subroutine calculate_f0_arrays(egrid) implicit none real, dimension (:,:), intent (in) :: egrid integer:: negrid,is,ie real:: v, f0, df0dE, f0prim negrid = size(egrid(:,1)) allocate(nonmaxw_corr(negrid,nspec)) allocate(f0_values(negrid,nspec)) allocate(dlnf0drho(negrid,nspec)) do is = 1, nspec is_global = is if (spec(is)%f0type == f0_tabulated) then call calculate_f0_arrays_tabulated(is, egrid) else do ie = 1, negrid v = sqrt(egrid(ie,is)) select case (spec(is)%f0type) case (f0_maxwellian) f0_values(ie,is) = eval_f0(v) nonmaxw_corr(ie,is) = 1.0 dlnf0drho(ie,is) = -spec(is)%fprim - (egrid(ie,is) - 1.5) * spec(is)%tprim case (f0_sdanalytic) call eval_f0_sdanalytic(is, v, f0, df0dE, f0prim) f0_values(ie, is) = f0 nonmaxw_corr(ie, is) = -spec(is)%temp * df0dE dlnf0drho(ie, is) = f0prim end select end do end if end do end subroutine calculate_f0_arrays !> FIXME : Add documentation subroutine set_current_f0_species(is) implicit none integer, intent(in):: is is_global = is end subroutine set_current_f0_species !> FIXME : Add documentation !> !> @note Relies upon is_global to determine species. This function only takes one argument so that !> it can be passed as an argument to genquad real function eval_f0(v) use constants, only: pi use mp, only: mp_abort implicit none real, intent(in):: v real::dummy1, dummy2 eval_f0 = 0.0 !Set a default, avoids warnings about maybe unset if (is_global <= 0) call mp_abort("is_global not specified when calling eval_f0") select case(spec(is_global)%f0type) case (f0_maxwellian) eval_f0 = exp(-v**2) / pi**1.5 case (f0_sdanalytic) call eval_f0_sdanalytic(is_global, v, eval_f0, dummy1, dummy2) case (f0_tabulated) call mp_abort("ERROR: eval_f0 cannot yet handle tabulated option.") end select end function eval_f0 !> FIXME : Add documentation subroutine eval_f0_sdanalytic(is, v, f0, df0dE_out, f0prim) use constants, only: pi implicit none integer,intent(in):: is real,intent(in):: v real,intent(out):: f0, df0dE_out, f0prim real:: A, df0dv, vcprim, vcrit ! Early exit -- if v > 1 then assume f0 has been cutoff, i.e. ! we are above the threshold for the Heaviside function to ! zero the distribution. Note for v==1 our f0prim etc. are ! don't account for the Heaviside. if (v > 1.0) then f0 = 0.0 f0prim = 0.0 df0dE_out = 1.0 return end if if (spec(is)%vcrit <= 0.0) then ! If vcrit hasn't been set work out vcrit and vcprim using ! eqn 1.17 and 1.25 of [G. Wilkie thesis](https://drum.lib.umd.edu/handle/1903/17302) ! Note this is independent of `v` so we could precalculate for the relevant species. ! This could be "trivially" cached by simply storing the calculated values in spec(is) call calculate_slowingdown_parameters(is, vcrit, vcprim) ! This line is a bit odd, it _looks_ like it is checking if spec(is)%vcprim ! has not been set, and if it _hasn't_ then force vcprim = spec(is)%vcprim. ! In fact as spec(is)%vcprim defaults to -999.9 currently this logic _doesn't_ ! trigger without the user manually setting vcprim more negative. I suspect ! this was intended to be if(was_set(spec(is)%vcprim)) vcprim = spec(is)%vcprim if (abs(spec(is)%vcprim + 999.0) < 1.0e-3) vcprim = spec(is)%vcprim else vcrit = spec(is)%vcrit ! If vcrit is specified, there is no way of knowing its radial ! derivative, so set dvc/drho = 0. As above, this doesn't actually trigger ! for the default value of spec(is)%vcprim. I suspect this is meant to be ! if(not_set(spec(is)%vcprim)) vcprim = 0. if (abs(spec(is)%vcprim + 999.0) < 1.0e-3) vcprim = 0.0 end if A = (3 / (4 * pi)) / log((vcrit**3 + 1) / (vcrit**3)) f0 = A / (vcrit**3 + v**3) df0dv = -f0 * 3 * v**2 / (vcrit**3 + v**3) df0dE_out = (0.5 / v) * df0dv / ( f0 * spec(is)%temp ) ! See eqn 1.22 of Wilkie's thesis with v_alpha -> 1 f0prim = -spec(is)%fprim - (-(vcrit**3 / (vcrit**3 + v**3)) + & (1 / ( log(1 + vcrit**(-3)) * (1 + vcrit**3)))) * vcprim end subroutine eval_f0_sdanalytic !> This subroutine calculates f0 on the grid from an external !! input file. The grid on the input file can differ from that !! of gs2. A cubic spline is used to interpolate between the two. !! The user can either specify df0/dE or it can be estimated internally. !! The radial derivative of f0 as a function of energy should also be given. subroutine calculate_f0_arrays_tabulated(is,egrid) ! ! The egrid is already given by get_legendre_grids. ! ! The input file can take two forms, as controlled by the ! control parameter mode, the first integer read from the file ! mode = 1: E points in first columns, F0(E) in second. To calculate ! generalized temperature, a cubic spline is used to ! interpolate and the slope is taken from that. ! = 2: First column is E, second is F0(E), the third is (1/F0)*dF0/dE ! = 3: E, F0(E), -(1/F0)*dF0/drho ! = 4: E, F0(E),dF0/dE, -(1/F0)*dF0/drho use mp, only: broadcast use constants, only: pi use file_utils, only: run_name use splines, only: spline, new_spline, delete_spline, splint, dsplint implicit none integer, intent(in) :: is real, dimension(:,:),intent(in):: egrid type(spline) :: spl integer:: f0in_unit, mode, ie, numdat, negrid logical:: acceptable_spline real:: moment0, sigma, sigma_fac, df0dE real, dimension(:), allocatable:: f0_values_dat, df0dE_dat, egrid_dat, df0drho_dat ! Open file and read column option open(newunit = f0in_unit,file=trim(run_name)//'.f0in',status='old',action='read') read(f0in_unit,*) mode read(f0in_unit,*) numdat allocate(f0_values_dat(numdat)) allocate(df0dE_dat(numdat)) allocate(df0drho_dat(numdat)) allocate(egrid_dat(numdat)) negrid = size(egrid(:,1)) sigma = 1.0 sigma_fac = 2.0 acceptable_spline = .false. if (mode == 1) then ! Read f0 values do ie = 1,numdat read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie) end do close(f0in_unit) ! Perform cubic spline to get F0 and its slope at grid points ! If F0 ends up being negative or zero anywhere on the grid, double the tension ! factor and try again do while (.not. acceptable_spline) ! Generate spline parameters spl = new_spline(egrid_dat, f0_values_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get f0 at grid points f0_values(ie,is) = splint(egrid(ie, is), spl) ! Calculate d/dE lnF0 to get generalised temperature df0dE = dsplint(egrid(ie,is), spl) / (f0_values(ie,is) * spec(is)%temp) nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE end do call delete_spline(spl) if (minval(f0_values(:,is)) > epsilon(0.0)) then acceptable_spline = .true. else sigma = sigma * sigma_fac end if end do dlnf0drho(:,is) = - spec(is)%fprim else if (mode == 2) then ! Read both f0 and df0/dE do ie = 1,numdat read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie) end do close(f0in_unit) do while (.not. acceptable_spline) ! Generate spline parameters for f0 spl = new_spline(egrid_dat, f0_values_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get F0 at grid points f0_values(ie,is) = splint(egrid(ie, is), spl) end do call delete_spline(spl) if (minval(f0_values(:,is)) > epsilon(0.0)) then acceptable_spline = .true. else sigma = sigma * sigma_fac cycle end if ! Generate spline parameters for df0/dE spl = new_spline(egrid_dat, df0dE_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get f0 at grid points df0dE = splint(egrid(ie,is), spl) / spec(is)%temp nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE end do call delete_spline(spl) end do dlnf0drho(:,is) = - spec(is)%fprim else if (mode == 3) then ! Read f0 and df0/drho do ie = 1,numdat read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0drho_dat(ie) end do close(f0in_unit) do while (.not. acceptable_spline) ! Generate spline parameters for f0 spl = new_spline(egrid_dat, f0_values_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get F0 at grid points f0_values(ie,is) = splint(egrid(ie,is), spl) ! Calculate d/dE lnF0 to get generalised temperature df0dE = dsplint(egrid(ie,is), spl) / (f0_values(ie,is) * spec(is)%temp) nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE end do call delete_spline(spl) if (minval(f0_values(:,is)) > epsilon(0.0)) then acceptable_spline = .true. else sigma = sigma * sigma_fac cycle end if ! Generate spline parameters for df0/drho spl = new_spline(egrid_dat, df0drho_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get f0prim at grid points dlnf0drho(ie,is) = splint(egrid(ie,is), spl) end do call delete_spline(spl) end do else if (mode == 4) then ! Read f0, df0/dE, df0/drho do ie = 1,numdat read(f0in_unit,*) egrid_dat(ie), f0_values_dat(ie), df0dE_dat(ie), df0drho_dat(ie) end do close(f0in_unit) do while (.not. acceptable_spline) ! Generate spline parameters for f0 spl = new_spline(egrid_dat, f0_values_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get F0 at grid points f0_values(ie,is) = splint(egrid(ie,is), spl) end do call delete_spline(spl) if (minval(f0_values(:,is)) > epsilon(0.0)) then acceptable_spline = .true. else sigma = sigma * sigma_fac cycle end if ! Generate spline parameters for df0/dE spl = new_spline(egrid_dat, df0dE_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get f0 at grid points df0dE = splint(egrid(ie,is), spl) / spec(is)%temp nonmaxw_corr(ie,is) = -spec(is)%temp*df0dE end do call delete_spline(spl) ! Generate spline parameters for df0/drho spl = new_spline(egrid_dat, df0drho_dat, tension = sigma) do ie = 1,negrid ! Interpolate to get f0 at grid points dlnf0drho(ie,is) = splint(egrid(ie,is), spl) end do call delete_spline(spl) end do dlnf0drho(:,is) = - spec(is)%fprim else write(*,*) "ERROR. First line in f0 input file should be mode" stop 1 end if ! Now calculate moments of F0 from trapezoidal rule. ! Should probably do this by quadrature, ! but le_grids is not defined yet moment0 = 0.0 do ie = 1,numdat-1 moment0 = moment0 + 0.5*(egrid_dat(ie+1)-egrid_dat(ie)) * & ( sqrt(egrid_dat(ie))*f0_values_dat(ie) + & sqrt(egrid_dat(ie+1))*f0_values_dat(ie+1) ) end do moment0 = moment0*4.0*pi deallocate(egrid_dat,f0_values_dat,df0dE_dat,df0drho_dat) end subroutine calculate_f0_arrays_tabulated !> Determine if any of the known species are electron like !> !> @note This should be effectively constant, could cache during init pure logical function has_electron_species (spec) implicit none type (specie), dimension (:), intent (in) :: spec has_electron_species = any(is_electron_species(spec)) end function has_electron_species !> Determine if any of the known species are hybrid electrons !> !> @note This should be effectively constant, could cache during init pure logical function has_hybrid_electron_species (spec) implicit none type (specie), dimension (:), intent (in) :: spec has_hybrid_electron_species = any(is_hybrid_electron_species(spec)) end function has_hybrid_electron_species !> Determine if any of the known species are ion like !> !> @note This should be effectively constant, could cache during init pure logical function has_ion_species (spec) implicit none type (specie), dimension (:), intent (in) :: spec has_ion_species = any(is_ion_species(spec)) end function has_ion_species !> Determine if the passed species instance is considered electron !> like. In other words if spec%type matches any of the recognised !> electron like types (electron, hybrid_electron) elemental logical function is_electron_species(spec) implicit none type(specie), intent(in) :: spec is_electron_species = any(spec%type == electron_like_species) end function is_electron_species !> Determine if the passed species instance corresponds to hybrid electrons elemental logical function is_hybrid_electron_species(spec) implicit none type(specie), intent(in) :: spec is_hybrid_electron_species = spec%type == hybrid_electron_species end function is_hybrid_electron_species !> Determine if the passed species instance corresponds to ions elemental logical function is_ion_species(spec) implicit none type(specie), intent(in) :: spec is_ion_species = spec%type == ion_species end function is_ion_species !> Determine if any of the known species are non-Maxwellian pure logical function has_nonmaxw_species (spec) implicit none type (specie), dimension (:), intent (in) :: spec has_nonmaxw_species = .NOT. all(spec%f0type == f0_maxwellian) end function has_nonmaxw_species !> FIXME : Add documentation subroutine finish_species implicit none if(allocated(spec)) deallocate (spec) ! if(allocated(f0_values)) deallocate (f0_values,nonmaxw_corr,dlnf0drho) call species_config%reset() if (allocated(species_element_config)) deallocate(species_element_config) initialized = .false. end subroutine finish_species !> FIXME : Add documentation subroutine set_overrides(prof_ov) use overrides, only: profiles_overrides_type use unit_tests, only: debug_message type(profiles_overrides_type), intent(in) :: prof_ov integer :: is integer, parameter :: verbosity=3 character(len=800) :: message call debug_message(verbosity, 'species::set_overrides starting') do is = 1,nspec message = '' write(message,*) 'species::set_overrides setting species ', is call debug_message(verbosity, trim(message)) if (prof_ov%override_tprim(is)) spec(is)%tprim = prof_ov%tprim(is) if (prof_ov%override_fprim(is)) spec(is)%fprim = prof_ov%fprim(is) if (prof_ov%override_temp(is)) spec(is)%temp = prof_ov%temp(is) if (prof_ov%override_dens(is)) spec(is)%dens = prof_ov%dens(is) if (prof_ov%override_vnewk(is)) spec(is)%vnewk = prof_ov%vnewk(is) end do call debug_message(verbosity,& 'species::set_overrides calling calculate_and_broadcast_species_properties') call calculate_and_broadcast_species_properties call debug_message(verbosity, 'species::set_overrides finishing') end subroutine set_overrides !> FIXME : Add documentation subroutine calculate_and_broadcast_species_properties use mp, only: proc0, broadcast integer :: is if (proc0) then do is = 1, nspec spec(is)%stm = sqrt(spec(is)%temp/spec(is)%mass) spec(is)%zstm = spec(is)%z/sqrt(spec(is)%temp*spec(is)%mass) spec(is)%tz = spec(is)%temp/spec(is)%z spec(is)%zt = spec(is)%z/spec(is)%temp spec(is)%smz = abs(sqrt(spec(is)%temp*spec(is)%mass)/spec(is)%z) end do end if call broadcast (nspec) do is = 1, nspec call broadcast (spec(is)%dens) call broadcast (spec(is)%temp) call broadcast (spec(is)%fprim) call broadcast (spec(is)%tprim) call broadcast (spec(is)%vnewk) call broadcast (spec(is)%stm) call broadcast (spec(is)%zstm) call broadcast (spec(is)%tz) call broadcast (spec(is)%zt) call broadcast (spec(is)%smz) end do end subroutine calculate_and_broadcast_species_properties !> Set the module level config types !> Will abort if the module has already been initialised to avoid !> inconsistencies. subroutine set_species_config(species_config_in, species_element_config_in) use mp, only: mp_abort type(species_config_type), intent(in), optional :: species_config_in type(species_element_config_type), intent(in), dimension(:), allocatable, optional :: species_element_config_in if (initialized) then call mp_abort("Trying to set species_config when already initialized.", to_screen = .true.) end if if (present(species_config_in)) then species_config = species_config_in end if if (present(species_element_config_in)) then species_element_config = species_element_config_in end if end subroutine set_species_config #include "species_auto_gen.inc" #include "species_element_auto_gen.inc" end module species