diff --git a/Project.toml b/Project.toml index a33a8e3f4..76563d054 100644 --- a/Project.toml +++ b/Project.toml @@ -43,6 +43,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [compat] julia = "1.7.0" diff --git a/docs/src/utils.md b/docs/src/utils.md new file mode 100644 index 000000000..5394ba87d --- /dev/null +++ b/docs/src/utils.md @@ -0,0 +1,6 @@ +`utils` +=============== + +```@autodocs +Modules = [moment_kinetics.utils] +``` diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml new file mode 100644 index 000000000..e9d44cc5b --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -0,0 +1,94 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +T_e = 1.0 +T_wall = 1.0 +rhostar = 1.0 +initial_density1 = 0.5 +initial_temperature1 = 1.0 +initial_density2 = 0.5 +initial_temperature2 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.001 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +z_IC_option2 = "sinusoid" +z_IC_density_amplitude2 = 0.001 +z_IC_density_phase2 = 0.0 +z_IC_upar_amplitude2 = 0.0 +z_IC_upar_phase2 = 0.0 +z_IC_temperature_amplitude2 = 0.0 +z_IC_temperature_phase2 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +krook_collisions_option = "reference_parameters" +nuii_krook = 1.0 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 17 +z_nelement = 16 +z_nelement_local = 16 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 17 +vpa_nelement = 16 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 1 +vperp_nelement = 1 +vperp_L = 6.0 +vperp_bc = "periodic" +vperp_discretization = "chebyshev_pseudospectral" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = -1.0 +z_dissipation_coefficient = -1.0 +r_dissipation_coefficient = -1.0 diff --git a/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml new file mode 100644 index 000000000..ae02bfd9c --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -0,0 +1,94 @@ +n_ion_species = 1 +n_neutral_species = 0 +electron_physics = "boltzmann_electron_response" +#electron_physics = "boltzmann_electron_response_with_simple_sheath" +evolve_moments_density = false +evolve_moments_parallel_flow = false +evolve_moments_parallel_pressure = false +evolve_moments_conservation = false +force_Er_zero_at_wall = false #true +Er_constant = 0.0 +T_e = 1.0 +T_wall = 1.0 +rhostar = 1.0 +initial_density1 = 0.5 +initial_temperature1 = 1.0 +initial_density2 = 0.5 +initial_temperature2 = 1.0 +z_IC_option1 = "sinusoid" +z_IC_density_amplitude1 = 0.001 +z_IC_density_phase1 = 0.0 +z_IC_upar_amplitude1 = 0.0 +z_IC_upar_phase1 = 0.0 +z_IC_temperature_amplitude1 = 0.0 +z_IC_temperature_phase1 = 0.0 +z_IC_option2 = "sinusoid" +z_IC_density_amplitude2 = 0.001 +z_IC_density_phase2 = 0.0 +z_IC_upar_amplitude2 = 0.0 +z_IC_upar_phase2 = 0.0 +z_IC_temperature_amplitude2 = 0.0 +z_IC_temperature_phase2 = 0.0 +charge_exchange_frequency = 0.0 +ionization_frequency = 0.0 +krook_collisions_option = "reference_parameters" +nuii_krook = 1.0 +nstep = 2000 +dt = 0.0005 +nwrite = 200 +nwrite_dfns = 200 +use_semi_lagrange = false +n_rk_stages = 4 +split_operators = false +z_ngrid = 17 +z_nelement = 16 +z_nelement_local = 16 +z_bc = "wall" +z_discretization = "chebyshev_pseudospectral" +r_ngrid = 1 +r_nelement = 1 +r_nelement_local = 1 +r_bc = "periodic" +r_discretization = "chebyshev_pseudospectral" +vpa_ngrid = 17 +vpa_nelement = 8 +vpa_L = 12.0 +vpa_bc = "zero" +vpa_discretization = "chebyshev_pseudospectral" +vperp_ngrid = 17 +vperp_nelement = 8 +vperp_L = 6.0 +vperp_bc = "periodic" +vperp_discretization = "chebyshev_pseudospectral" + +vz_ngrid = 17 +vz_nelement = 4 +vz_L = 12.0 +vz_bc = "periodic" +vz_discretization = "chebyshev_pseudospectral" + +vr_ngrid = 17 +vr_nelement = 4 +vr_L = 12.0 +vr_bc = "periodic" +vr_discretization = "chebyshev_pseudospectral" + +vzeta_ngrid = 17 +vzeta_nelement = 4 +vzeta_L = 12.0 +vzeta_bc = "periodic" +vzeta_discretization = "chebyshev_pseudospectral" + +[manufactured_solns] + use_for_advance=true + use_for_init=true + # constant to be used to control Ez divergence in MMS tests + epsilon_offset=0.1 + # bool to control if dfni is a function of vpa or vpabar in MMS test + use_vpabar_in_mms_dfni=true + alpha_switch=1.0 + type="default" +[numerical_dissipation] +vpa_dissipation_coefficient = -1.0 +z_dissipation_coefficient = -1.0 +r_dissipation_coefficient = -1.0 diff --git a/src/constants.jl b/src/constants.jl new file mode 100644 index 000000000..7f5214fc9 --- /dev/null +++ b/src/constants.jl @@ -0,0 +1,33 @@ +""" +Some physical constants +""" +module constants + +export epsilon0, mu0 +export electron_mass +export proton_charge, proton_mass +export deuteron_mass +export amu + +# https://physics.nist.gov/cgi-bin/cuu/Value?ep0 +const epsilon0 = 8.8541878128e-12 # F m^-1 + +# https://physics.nist.gov/cgi-bin/cuu/Value?mu0 +const mu0 = 1.25663706212e-6 # N A^-2 + +# https://physics.nist.gov/cgi-bin/cuu/Value?me +const electron_mass = 9.109383701e-31 # kg + +# https://physics.nist.gov/cgi-bin/cuu/Value?e +const proton_charge = 1.602176634e-19 # C + +# https://physics.nist.gov/cgi-bin/cuu/Value?mp +const proton_mass = 1.67262192369e-27 # kg + +# https://physics.nist.gov/cgi-bin/cuu/Value?md +const deuteron_mass = 3.3435837724e-27 + +# https://physics.nist.gov/cgi-bin/cuu/Value?ukg +const amu = 1.66053906660e-27 + +end diff --git a/src/file_io.jl b/src/file_io.jl index d77ae35dd..1a53c3cca 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -61,6 +61,8 @@ struct io_moments_info{Tfile, Ttime, Tphi, Tmomi, Tmomn} parallel_flow::Tmomi # handle for the charged species parallel pressure parallel_pressure::Tmomi + # handle for the charged species perpendicular pressure + perpendicular_pressure::Tmomi # handle for the charged species parallel heat flux parallel_heat_flux::Tmomi # handle for the charged species thermal speed @@ -591,6 +593,13 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species, parallel_io=parallel_io, description="charged species parallel pressure", units="n_ref*T_ref") + + # io_pperp is the handle for the ion parallel pressure + io_pperp = create_dynamic_variable!(dynamic, "perpendicular_pressure", mk_float, z, r; + n_ion_species=n_ion_species, + parallel_io=parallel_io, + description="charged species perpendicular pressure", + units="n_ref*T_ref") # io_qpar is the handle for the ion parallel heat flux io_qpar = create_dynamic_variable!(dynamic, "parallel_heat_flux", mk_float, z, r; @@ -642,7 +651,7 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species, units="c_ref") return io_moments_info(fid, io_time, io_phi, io_Er, io_Ez, io_density, io_upar, - io_ppar, io_qpar, io_vth, io_density_neutral, io_uz_neutral, + io_ppar, io_pperp, io_qpar, io_vth, io_density_neutral, io_uz_neutral, io_pz_neutral, io_qz_neutral, io_thermal_speed_neutral, parallel_io) end @@ -780,7 +789,8 @@ function reopen_moments_io(file_info) end return io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"), getvar("Ez"), getvar("density"), getvar("parallel_flow"), - getvar("parallel_pressure"), getvar("parallel_heat_flux"), + getvar("parallel_pressure"), getvar("perpendicular_pressure"), + getvar("parallel_heat_flux"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), getvar("qz_neutral"), getvar("thermal_speed_neutral"), @@ -865,6 +875,7 @@ function reopen_dfns_io(file_info) io_moments = io_moments_info(fid, getvar("time"), getvar("phi"), getvar("Er"), getvar("Ez"), getvar("density"), getvar("parallel_flow"), getvar("parallel_pressure"), + getvar("perpendicular_pressure"), getvar("parallel_heat_flux"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), @@ -920,6 +931,8 @@ function write_moments_data_to_binary(moments, fields, t, n_ion_species, n_ion_species) append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar, t_idx, z, r, n_ion_species) + append_to_dynamic_var(io_moments.perpendicular_pressure, moments.charged.pperp, t_idx, + z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_heat_flux, moments.charged.qpar, t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth, t_idx, z, r, @@ -1008,6 +1021,8 @@ end t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_pressure, moments.charged.ppar.data, t_idx, z, r, n_ion_species) + append_to_dynamic_var(io_moments.perpendicular_pressure.data, moments.charged.pperp, + t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.parallel_heat_flux, moments.charged.qpar.data, t_idx, z, r, n_ion_species) append_to_dynamic_var(io_moments.thermal_speed, moments.charged.vth.data, @@ -1271,7 +1286,7 @@ all the arrays have the same length, with an entry for each call to `debug_dump( function debug_dump end function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coordinate, vperp::coordinate, z::coordinate, r::coordinate, t::mk_float; - ff=nothing, dens=nothing, upar=nothing, ppar=nothing, qpar=nothing, + ff=nothing, dens=nothing, upar=nothing, ppar=nothing, pperp=nothing, qpar=nothing, vth=nothing, ff_neutral=nothing, dens_neutral=nothing, uz_neutral=nothing, #ur_neutral=nothing, uzeta_neutral=nothing, @@ -1370,6 +1385,11 @@ function debug_dump(vz::coordinate, vr::coordinate, vzeta::coordinate, vpa::coor else debug_output_file.moments.parallel_pressure[:,:,:,debug_output_counter[]] = ppar end + if pperp === nothing + debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = 0.0 + else + debug_output_file.moments.perpendicular_pressure[:,:,:,debug_output_counter[]] = pperp + end if qpar === nothing debug_output_file.moments.parallel_heat_flux[:,:,:,debug_output_counter[]] = 0.0 else @@ -1446,6 +1466,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, density = nothing upar = nothing ppar = nothing + pperp = nothing pdf_neutral = nothing density_neutral = nothing else @@ -1453,6 +1474,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, density = fvec.density upar = fvec.upar ppar = fvec.ppar + pperp = fvec.pperp pdf_neutral = fvec.pdf_neutral density_neutral = fvec.density_neutral end @@ -1466,7 +1488,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, Ez = fields.Ez end return debug_dump(vz, vr, vzeta, vpa, vperp, z, r, t; ff=pdf, dens=density, upar=upar, - ppar=ppar, ff_neutral=pdf_neutral, dens_neutral=density_neutral, + ppar=ppar, pperp=pperp, ff_neutral=pdf_neutral, dens_neutral=density_neutral, phi=phi, Er=Er, Ez=Ez, t, istage=istage, label=label) end diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl index f4fe358ec..56fc0444d 100644 --- a/src/initial_conditions.jl +++ b/src/initial_conditions.jl @@ -31,7 +31,7 @@ using ..velocity_moments: create_moments_charged, create_moments_neutral, update using ..velocity_moments: moments_charged_substruct, moments_neutral_substruct using ..velocity_moments: update_neutral_density!, update_neutral_pz!, update_neutral_pr!, update_neutral_pzeta! using ..velocity_moments: update_neutral_uz!, update_neutral_ur!, update_neutral_uzeta!, update_neutral_qz! -using ..velocity_moments: update_ppar!, update_upar!, update_density!, reset_moments_status! +using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, update_vth!, reset_moments_status! using ..manufactured_solns: manufactured_solutions @@ -157,8 +157,9 @@ function init_pdf_and_moments!(pdf, moments, boundary_distributions, geometry, init_upar!(moments.charged.upar, z, r, species.charged, n_ion_species) # initialise the parallel thermal speed profile init_vth!(moments.charged.vth, z, r, species.charged, n_ion_species) + # initialise pressures assuming isotropic distribution @. moments.charged.ppar = 0.5 * moments.charged.dens * moments.charged.vth^2 - + @. moments.charged.pperp = moments.charged.ppar if n_neutral_species > 0 #neutral particles init_density!(moments.neutral.dens, z, r, species.neutral, n_neutral_species) @@ -898,7 +899,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, manufactured_solns_input) manufactured_solns_list = manufactured_solutions(manufactured_solns_input, r.L, z.L, r.bc, z.bc, geometry, composition, - species, r.n) + species, r.n, vperp.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func dfnn_func = manufactured_solns_list.dfnn_func @@ -924,17 +925,13 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, moments.charged.dens, moments.charged.upar, pdf.charged.norm, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar) + update_pperp!(moments.charged.pperp, pdf.charged.norm, vpa, vperp, z, r, composition) update_qpar!(moments.charged.qpar, moments.charged.qpar_updated, moments.charged.dens, moments.charged.upar, moments.charged.vth, pdf.charged.norm, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar, moments.evolve_ppar) - # convert from particle particle flux to parallel flow - begin_s_r_z_region() - @loop_s_r_z is ir iz begin - # update the thermal speed - moments.charged.vth[iz,ir,is] = sqrt(2.0*moments.charged.ppar[iz,ir,is]/moments.charged.dens[iz,ir,is]) - end + update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition) if n_neutral_species > 0 begin_sn_r_z_region() diff --git a/src/input_structs.jl b/src/input_structs.jl index a08168e58..be8a19969 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -14,6 +14,8 @@ export collisions_input export io_input export pp_input export geometry_input +export set_defaults_and_check_top_level!, set_defaults_and_check_section!, + Dict_to_NamedTuple using ..type_definitions: mk_float, mk_int @@ -59,6 +61,7 @@ mutable struct advance_info ionization_collisions::Bool ionization_collisions_1V::Bool ionization_source::Bool + krook_collisions::Bool numerical_dissipation::Bool source_terms::Bool continuity::Bool @@ -260,13 +263,7 @@ mutable struct species_composition phi_wall::mk_float # constant for testing nonzero Er Er_constant::mk_float - # constant controlling divergence at wall boundaries in MMS test - epsilon_offset::mk_float - # logical controlling whether or not dfni(vpabar,z,r) or dfni(vpa,z,r) in MMS test - use_vpabar_in_mms_dfni::Bool - # associated float controlling form of assumed potential in MMS test - alpha_switch::mk_float - # ratio of the neutral particle mass to the ion mass + # ratio of the neutral particle mass to the ion mass mn_over_mi::mk_float # ratio of the electron particle mass to the ion mass me_over_mi::mk_float @@ -307,6 +304,10 @@ mutable struct collisions_input ionization::mk_float # if constant_ionization_rate = true, use an ionization term that is constant in z constant_ionization_rate::Bool + # Coulomb collision rate at the reference density and temperature + krook_collision_frequency_prefactor::mk_float + # Setting to switch between different options for Krook collision operator + krook_collisions_option::String end """ @@ -506,4 +507,112 @@ struct pp_input diagnostics_chodura_r::Bool end +import Base: get +""" +Utility method for converting a string to an Enum when getting from a Dict, based on the +type of the default value +""" +function get(d::Dict, key, default::Enum) + valstring = get(d, key, nothing) + if valstring == nothing + return default + # instances(typeof(default)) gets the possible values of the Enum. Then convert to + # Symbol, then to String. + elseif valstring ∈ (split(s, ".")[end] for s ∈ String.(Symbol.(instances(typeof(default))))) + return eval(Symbol(valstring)) + else + error("Expected a $(typeof(default)), but '$valstring' is not in " + * "$(instances(typeof(default)))") + end +end + +""" +Set the defaults for options in the top level of the input, and check that there are not +any unexpected options (i.e. options that have no default). + +Modifies the options[section_name]::Dict by adding defaults for any values that are not +already present. + +Ignores any sections, as these will be checked separately. +""" +function set_defaults_and_check_top_level!(options::AbstractDict; kwargs...) + DictType = typeof(options) + + # Check for any unexpected values in the options - all options that are set should be + # present in the kwargs of this function call + options_keys_symbols = keys(kwargs) + options_keys = (String(k) for k ∈ options_keys_symbols) + for (key, value) in options + # Ignore any ssections when checking + if !(isa(value, AbstractDict) || key ∈ options_keys) + error("Unexpected option '$key=$value' in top-level options") + end + end + + # Set default values if a key was not set explicitly + explicit_keys = keys(options) + for (key_sym, value) ∈ kwargs + key = String(key_sym) + if !(key ∈ explicit_keys) + options[key] = value + end + end + + return options +end + +""" +Set the defaults for options in a section, and check that there are not any unexpected +options (i.e. options that have no default). + +Modifies the options[section_name]::Dict by adding defaults for any values that are not +already present. +""" +function set_defaults_and_check_section!(options::AbstractDict, section_name; + kwargs...) + DictType = typeof(options) + + if !(section_name ∈ keys(options)) + # If section is not present, create it + options[section_name] = DictType() + end + + if !isa(options[section_name], AbstractDict) + error("Expected '$section_name' to be a section in the input file, but it has a " + * "value '$(options[section_name])'") + end + + section = options[section_name] + + # Check for any unexpected values in the section - all options that are set should be + # present in the kwargs of this function call + section_keys_symbols = keys(kwargs) + section_keys = (String(k) for k ∈ section_keys_symbols) + for (key, value) in section + if !(key ∈ section_keys) + error("Unexpected option '$key=$value' in section '$section_name'") + end + end + + # Set default values if a key was not set explicitly + explicit_keys = keys(section) + for (key_sym, value) ∈ kwargs + key = String(key_sym) + if !(key ∈ explicit_keys) + section[key] = value + end + end + + return section +end + +""" +Convert a Dict whose keys are String or Symbol to a NamedTuple + +Useful as NamedTuple is immutable, so option values cannot be accidentally changed. +""" +function Dict_to_NamedTuple(d) + return NamedTuple(Symbol(k)=>v for (k,v) ∈ d) +end + end diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl new file mode 100644 index 000000000..fe6477098 --- /dev/null +++ b/src/krook_collisions.jl @@ -0,0 +1,161 @@ +""" +""" +module krook_collisions + +export setup_krook_collisions, get_collision_frequency, krook_collisions! + +using ..constants: epsilon0, proton_charge +using ..looping + +""" +Calculate normalized collision frequency at reference parameters for Coulomb collisions. + +Currently valid only for hydrogenic ions (Z=1) +""" +function setup_krook_collisions(reference_params) + Nref = reference_params.Nref + Tref = reference_params.Tref + mref = reference_params.mref + timeref = reference_params.timeref + cref = reference_params.cref + + Nref_per_cm3 = Nref * 1.0e-6 + + # Coulomb logarithm at reference parameters for same-species ion-ion collisions, using + # NRL formulary. Formula given for n in units of cm^-3 and T in units of eV. + logLambda_ii = 23.0 - log(sqrt(2.0*Nref_per_cm3) / Tref^1.5) + + # Collision frequency, using \hat{\nu} from Appendix, p. 277 of Helander "Collisional + # Transport in Magnetized Plasmas" (2002). + nu_ii0_per_s = Nref * proton_charge^4 * logLambda_ii / + (4.0 * π * epsilon0^2 * mref^2 * cref^3) # s^-1 + nu_ii0 = nu_ii0_per_s * timeref + + return nu_ii0 +end + +""" + get_collision_frequency(collisions, n, vth) + +Calculate the collision frequency, depending on the settings/parameters in `collisions`, +for the given density `n` and thermal speed `vth`. + +`n` and `vth` may be scalars or arrays, but should have shapes that can be broadcasted +together. +""" +function get_collision_frequency(collisions, n, vth) + if collisions.krook_collisions_option == "reference_parameters" + return @. collisions.krook_collision_frequency_prefactor * n * vth^(-3) + elseif collisions.krook_collisions_option == "manual" + # Include 0.0*n so that the result gets promoted to an array if n is an array, + # which hopefully means this function will have a fixed return type given the + # types of the arguments (we don't want to be 'type unstable' for array inputs by + # returning a scalar from this branch but an array from the "reference_parameters" + # branch). + return @. collisions.krook_collision_frequency_prefactor + 0.0 * n + elseif collisions.krook_collisions_option == "none" + # Include 0.0*n so that the result gets promoted to an array if n is an array, + # which hopefully means this function will have a fixed return type given the + # types of the arguments (we don't want to be 'type unstable' for array inputs by + # returning a scalar from this branch but an array from the "reference_parameters" + # branch). + return @. 0.0 * n + else + error("Unrecognised option " + * "krook_collisions_option=$(collisions.krook_collisions_option)") + end +end + +""" +Add collision operator + +Currently Krook collisions +""" +function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, vperp, vpa, dt) + begin_s_r_z_region() + + if vperp.n > 1 && (moments.evolve_density || moments.evolve_upar || moments.evolve_ppar) + error("Krook collisions not implemented for 2V moment-kinetic cases yet") + end + + # Note: do not need 1/sqrt(pi) for the 'Maxwellian' term because the pdf is already + # normalized by sqrt(pi) (see velocity_moments.integrate_over_vspace). + if moments.evolve_ppar && moments.evolve_upar + # Compared to evolve_upar version, grid is already normalized by vth, and multiply + # through by vth, remembering pdf is already multiplied by vth + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + nu_ii = get_collision_frequency(collisions, n, vth) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - exp(-vpa.grid[ivpa]^2 - vperp.grid[ivperp]^2)) + end + end + elseif moments.evolve_ppar + # Compared to full-f collision operater, multiply through by vth, remembering pdf + # is already multiplied by vth, and grid is already normalized by vth + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + nu_ii = get_collision_frequency(collisions, n, vth) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is])/vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + elseif moments.evolve_upar + # Compared to evolve_density version, grid is already shifted by upar + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + nu_ii = get_collision_frequency(collisions, n, vth) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - 1.0 / vth * exp(-(vpa.grid[ivpa] / vth)^2 + - (vperp.grid[ivperp] / vth)^2)) + end + end + elseif moments.evolve_density + # Compared to full-f collision operater, divide through by density, remembering + # that pdf is already normalized by density + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + nu_ii = get_collision_frequency(collisions, n, vth) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - 1.0 / vth + * exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is]) / vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + else + @loop_s_r_z is ir iz begin + n = fvec_in.density[iz,ir,is] + vth = moments.charged.vth[iz,ir,is] + if vperp.n == 1 + vth_prefactor = 1.0 / vth + else + vth_prefactor = 1.0 / vth^3 + end + nu_ii = get_collision_frequency(collisions, n, vth) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * + (fvec_in.pdf[ivpa,ivperp,iz,ir,is] + - n * vth_prefactor + * exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is])/vth)^2 + - (vperp.grid[ivperp]/vth)^2)) + end + end + end + + return nothing +end + +end # krook_collisions diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 1c4fd8600..b4925486d 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -10,7 +10,7 @@ julia --project run_makie_post_processing.jl dir1 [dir2 [dir3 ...]] """ module makie_post_processing -export makie_post_process, generate_example_input_file, +export makie_post_process, generate_example_input_file, get_variable, setup_makie_post_processing_input!, get_run_info, irregular_heatmap, irregular_heatmap!, postproc_load_variable, positive_or_nan @@ -18,11 +18,12 @@ using ..analysis: analyze_fields_data, check_Chodura_condition, get_r_perturbati get_Fourier_modes_2D, get_Fourier_modes_1D, steady_state_residuals using ..array_allocation: allocate_float using ..coordinates: define_coordinate -using ..input_structs: grid_input, advection_input +using ..input_structs: grid_input, advection_input, set_defaults_and_check_top_level!, + set_defaults_and_check_section!, Dict_to_NamedTuple +using ..krook_collisions: get_collision_frequency using ..looping: all_dimensions, ion_dimensions, neutral_dimensions using ..manufactured_solns: manufactured_solutions, manufactured_electric_fields -using ..moment_kinetics_input: mk_input, set_defaults_and_check_top_level!, - set_defaults_and_check_section!, Dict_to_NamedTuple +using ..moment_kinetics_input: mk_input using ..load_data: open_readonly_output_file, get_group, load_block_data, load_coordinate_data, load_distributed_charged_pdf_slice, load_distributed_neutral_pdf_slice, load_input, load_mk_options, @@ -68,7 +69,8 @@ const input_dict_dfns = OrderedDict{String,Any}() const em_variables = ("phi", "Er", "Ez") const ion_moment_variables = ("density", "parallel_flow", "parallel_pressure", - "thermal_speed", "temperature", "parallel_heat_flux") + "thermal_speed", "temperature", "parallel_heat_flux", + "collision_frequency") const neutral_moment_variables = ("density_neutral", "uz_neutral", "pz_neutral", "thermal_speed_neutral", "temperature_neutral", "qz_neutral") @@ -313,7 +315,12 @@ function makie_post_process(run_dir::Union{String,Tuple}, sound_wave_plots(run_info; plot_prefix=plot_prefix) - manufactured_solutions_analysis(run_info; plot_prefix=plot_prefix) + if all(ri === nothing for ri ∈ run_info_dfns) + nvperp = nothing + else + nvperp = maximum(ri.vperp.n_global for ri ∈ run_info_dfns if ri !== nothing) + end + manufactured_solutions_analysis(run_info; plot_prefix=plot_prefix, nvperp=nvperp) manufactured_solutions_analysis_dfns(run_info_dfns; plot_prefix=plot_prefix) return nothing @@ -673,6 +680,7 @@ function _setup_single_input!(this_input_dict::OrderedDict{String,Any}, (o=>section_defaults[String(o)] for o ∈ (:it0, :ir0, :iz0, :ivperp0, :ivpa0, :ivzeta0, :ivr0, :ivz0))..., colormap=this_input_dict["colormap"], animation_ext=this_input_dict["animation_ext"], + show_element_boundaries=this_input_dict["show_element_boundaries"], ) sort!(this_input_dict["manufactured_solns"]) @@ -778,10 +786,9 @@ function get_run_info(run_dir, restart_index=nothing; itime_min=1, itime_max=-1, # obtain input options from moment_kinetics_input.jl # and check input to catch errors - io_input, evolve_moments, t_input, z_input, r_input, vpa_input, vperp_input, - gyrophase_input, vz_input, vr_input, vzeta_input, composition, species, - collisions, geometry, drive_input, num_diss_params, manufactured_solns_input = - mk_input(input) + io_input, evolve_moments, t_input, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, + composition, species, collisions, geometry, drive_input, num_diss_params, + manufactured_solns_input = mk_input(input) n_ion_species, n_neutral_species = load_species_data(file_final_restart) evolve_density, evolve_upar, evolve_ppar = load_mk_options(file_final_restart) @@ -810,7 +817,7 @@ function get_run_info(run_dir, restart_index=nothing; itime_min=1, itime_max=-1, dummy_comm = MPI.COMM_NULL dummy_input = grid_input("dummy", 1, 1, 1, 1, 0, 1.0, "chebyshev_pseudospectral", "", "periodic", - dummy_adv_input, dummy_comm) + dummy_adv_input, dummy_comm, "uniform") vzeta, vzeta_spectral = define_coordinate(dummy_input) vzeta_chunk_size = 1 vr, vr_spectral = define_coordinate(dummy_input) @@ -1177,6 +1184,43 @@ function postproc_load_variable(run_info, variable_name; it=nothing, is=nothing, return result end +""" + get_variable(run_info::Tuple, variable_name; kwargs...) + get_variable(run_info, variable_name; kwargs...) + +Get an array (or Tuple of arrays, if `run_info` is a Tuple) of the data for +`variable_name` from `run_info`. + +Some derived variables need to be calculated from the saved output, not just loaded from +file (with `postproc_load_variable`). This function takes care of that calculation, and +handles the case where `run_info` is a Tuple (which `postproc_load_data` does not handle). + +`kwargs...` are passed through to `postproc_load_variable()`. +""" +function get_variable end + +function get_variable(run_info::Tuple, variable_name; kwargs...) + return Tuple(get_variable(ri, variable_name; kwargs...) for ri ∈ run_info) +end + +function get_variable(run_info, variable_name; kwargs...) + if variable_name == "temperature" + vth = postproc_load_variable(run_info, "thermal_speed") + variable = vth.^2 + elseif variable_name == "collision_frequency" + n = postproc_load_variable(run_info, "density") + vth = postproc_load_variable(run_info, "thermal_speed") + variable = get_collision_frequency(run_info.collisions, n, vth) + elseif variable_name == "temperature_neutral" + vth = postproc_load_variable(run_info, "thermal_speed_neutral") + variable = vth.^2 + else + variable = postproc_load_variable(run_info, variable_name) + end + + return variable +end + const chunk_size_1d = 10000 const chunk_size_2d = 100 struct VariableCache{T1,T2,N} @@ -1286,23 +1330,17 @@ function plots_for_variable(run_info, variable_name; plot_prefix, is_1D=false, if is_1D && variable_name == "Er" return nothing + elseif variable_name == "collision_frequency" && + all(ri.collisions.krook_collisions_option == "none" for ri ∈ run_info) + # No Krook collisions active, so do not make plots. + return nothing end println("Making plots for $variable_name") flush(stdout) - if variable_name == "temperature" - vth = Tuple(postproc_load_variable(ri, "thermal_speed") - for ri ∈ run_info) - variable = Tuple(v.^2 for v ∈ vth) - elseif variable_name == "temperature_neutral" - vth = Tuple(postproc_load_variable(ri, "thermal_speed_neutral") - for ri ∈ run_info) - variable = Tuple(v.^2 for v ∈ vth) - else - variable = Tuple(postproc_load_variable(ri, variable_name) - for ri ∈ run_info) - end + variable = get_variable(run_info, variable_name) + if variable_name ∈ em_variables species_indices = (nothing,) elseif variable_name ∈ neutral_moment_variables || @@ -1647,10 +1685,10 @@ for dim ∈ one_dimension_combinations end function $function_name(run_info, var_name; is=1, data=nothing, - input=nothing, ax=nothing, label=nothing, - outfile=nothing, it=nothing, ir=nothing, iz=nothing, - ivperp=nothing, ivpa=nothing, ivzeta=nothing, - ivr=nothing, ivz=nothing, kwargs...) + input=nothing, fig=nothing, ax=nothing, + label=nothing, outfile=nothing, it=nothing, + ir=nothing, iz=nothing, ivperp=nothing, ivpa=nothing, + ivzeta=nothing, ivr=nothing, ivz=nothing, kwargs...) if input === nothing if run_info.dfns input = input_dict_dfns[var_name] @@ -1686,7 +1724,7 @@ for dim ∈ one_dimension_combinations if $idim !== nothing x = x[$idim] end - fig = plot_1d(x, data; label=label, ax=ax, kwargs...) + plot_1d(x, data; label=label, ax=ax, kwargs...) if input.show_element_boundaries && Symbol($dim_str) != :t && ax_was_nothing element_boundary_positions = run_info.$dim.grid[begin:run_info.$dim.ngrid-1:end] @@ -5433,11 +5471,12 @@ Returns `variable`, `variable_sym`. """ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; it=nothing, ir=nothing, iz=nothing, ivperp=nothing, ivpa=nothing, ivzeta=nothing, - ivr=nothing, ivz=nothing) + ivr=nothing, ivz=nothing, nvperp) variable_name = Symbol(variable_name) func_name_lookup = (phi=:phi_func, Er=:Er_func, Ez=:Ez_func, density=:densi_func, + parallel_flow=:upari_func, parallel_pressure=:ppari_func, density_neutral=:densn_func, f=:dfni_func, f_neutral=:dfnn_func) nt = run_info.nt @@ -5467,11 +5506,13 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; run_info.z.bc, run_info.composition, run_info.r.n, run_info.manufactured_solns_input, run_info.species) - elseif variable_name ∈ (:density, :density_neutral, :f, :f_neutral) + elseif variable_name ∈ (:density, :parallel_flow, :parallel_pressure, + :density_neutral, :f, :f_neutral) manufactured_funcs = manufactured_solutions(run_info.manufactured_solns_input, Lr_in, run_info.z.L, run_info.r.bc, run_info.z.bc, run_info.geometry, - run_info.composition, run_info.species, run_info.r.n) + run_info.composition, run_info.species, run_info.r.n, + nvperp) end variable_func = manufactured_funcs[func_name_lookup[variable_name]] @@ -5558,7 +5599,8 @@ the label for the error (the difference between the computed and manufactured so If `io` is passed then error norms will be written to that file. """ function compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_sym_label, - norm_label, variable_name; io=nothing, input=nothing) + norm_label, variable_name; io=nothing, + input=nothing, nvperp) println("Doing MMS analysis and making plots for $variable_name") flush(stdout) @@ -5568,7 +5610,7 @@ function compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_ end field, field_sym = - manufactured_solutions_get_field_and_field_sym(run_info, variable_name) + manufactured_solutions_get_field_and_field_sym(run_info, variable_name; nvperp=nvperp) error = field .- field_sym nt = run_info.nt @@ -5780,13 +5822,19 @@ plots/animations, `field_sym_label` is the label for the manufactured solution, omitting `:t`. """ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, neutrals) nt = run_info.nt time = run_info.time - all_plot_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ plot_dims) - all_animate_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ animate_dims) + if neutrals + all_dims_no_t = (:r, :z, :vzeta, :vr, :vz) + else + all_dims_no_t = (:r, :z, :vperp, :vpa) + end + all_dims = tuple(:t, all_dims_no_t...) + all_plot_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ all_dims) + all_animate_slices = Tuple(Symbol(:i, d)=>input[Symbol(:i, d, :0)] for d ∈ all_dims_no_t) # Options to produce either regular or log-scale plots epsilon = 1.0e-30 # minimum data value to include in log plots @@ -5800,7 +5848,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label slices = (k=>v for (k, v) ∈ all_plot_slices if k != Symbol(:i, dim)) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; yscale=yscale, get_legend_place=:below) @@ -5825,7 +5873,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label if k ∉ (Symbol(:i, dim1), Symbol(:i, dim2))) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -5849,7 +5897,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label slices = (k=>v for (k, v) ∈ all_animate_slices if k != Symbol(:i, dim)) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; yscale=yscale, get_legend_place=:below) @@ -5877,7 +5925,7 @@ function _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label if k ∉ (Symbol(:i, dim1), Symbol(:i, dim2))) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name; slices...) + run_info, variable_name; nvperp=run_info.vperp.n, slices...) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -5984,7 +6032,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for r_chunk ∈ r_chunks, z_chunk ∈ z_chunks f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=it, ir=r_chunk, iz=z_chunk) + run_info, variable_name; nvperp=run_info.vperp.n, it=it, + ir=r_chunk, iz=z_chunk) dummy += sum(@. (f - f_sym)^2) #dummy_N += sum(f_sym.^2) end @@ -6004,8 +6053,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for (iz, z_label) ∈ ((1, "wall-"), (z.n, "wall+")) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, ir=input.ir0, iz=iz, - ivperp=input.ivperp0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, + ir=input.ir0, iz=iz, ivperp=input.ivperp0) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; get_legend_place=:below) @@ -6024,7 +6073,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1D f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ivperp=input.ivperp0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ivperp=input.ivperp0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6045,7 +6095,8 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1V f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6077,7 +6128,7 @@ function compare_charged_pdf_symbolic_test(run_info, plot_prefix; io=nothing, end plot_dims = tuple(:t, animate_dims...) _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, false) return field_norm end @@ -6166,7 +6217,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for r_chunk ∈ r_chunks, z_chunk ∈ z_chunks f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=it, ir=r_chunk, iz=z_chunk) + run_info, variable_name; nvperp=run_info.vperp.n, it=it, + ir=r_chunk, iz=z_chunk) dummy += sum(@. (f - f_sym)^2) #dummy_N += sum(f_sym.^2) end @@ -6186,8 +6238,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, for (iz, z_label) ∈ ((1, "wall-"), (z.n, "wall+")) f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, ir=input.ir0, iz=iz, - ivzeta=input.ivzeta0, ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, + ir=input.ir0, iz=iz, ivzeta=input.ivzeta0, ivr=input.ivr0) error = f .- f_sym fig, ax, legend_place = get_1d_ax(2; get_legend_place=:below) @@ -6206,8 +6258,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1D f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ivzeta=input.ivzeta0, - ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ivzeta=input.ivzeta0, ivr=input.ivr0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6228,8 +6280,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, if !is_1V f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0, - ivzeta=input.ivzeta0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0, ivzeta=input.ivzeta0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6251,8 +6303,8 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, f, f_sym = manufactured_solutions_get_field_and_field_sym( - run_info, variable_name, it=input.it0, iz=iz, ir=input.ir0, - ivr=input.ivr0) + run_info, variable_name; nvperp=run_info.vperp.n, it=input.it0, iz=iz, + ir=input.ir0, ivr=input.ivr0) error = f .- f_sym fig, ax, colorbar_place = get_2d_ax(3) @@ -6284,7 +6336,7 @@ function compare_neutral_pdf_symbolic_test(run_info, plot_prefix; io=nothing, end plot_dims = tuple(:t, animate_dims...) _MMS_pdf_plots(run_info, input, variable_name, plot_prefix, field_label, - field_sym_label, norm_label, plot_dims, animate_dims) + field_sym_label, norm_label, plot_dims, animate_dims, true) return field_norm end @@ -6311,7 +6363,7 @@ greater than one will result in an error. """ function manufactured_solutions_analysis end -function manufactured_solutions_analysis(run_info::Tuple; plot_prefix) +function manufactured_solutions_analysis(run_info::Tuple; plot_prefix, nvperp) if !any(ri !== nothing && ri.manufactured_solns_input.use_for_advance && ri.manufactured_solns_input.use_for_init for ri ∈ run_info) # No manufactured solutions tests @@ -6330,18 +6382,24 @@ function manufactured_solutions_analysis(run_info::Tuple; plot_prefix) return nothing end try - return manufactured_solutions_analysis(run_info[1]; plot_prefix=plot_prefix) + return manufactured_solutions_analysis(run_info[1]; plot_prefix=plot_prefix, + nvperp=nvperp) catch e println("Error in manufactured_solutions_analysis(). Error was ", e) end end -function manufactured_solutions_analysis(run_info; plot_prefix) +function manufactured_solutions_analysis(run_info; plot_prefix, nvperp) manufactured_solns_input = run_info.manufactured_solns_input if !(manufactured_solns_input.use_for_advance && manufactured_solns_input.use_for_init) return nothing end + if nvperp === nothing + error("No `nvperp` found - must have distributions function outputs to plot MMS " + * "tests") + end + input = Dict_to_NamedTuple(input_dict["manufactured_solns"]) open(run_info.run_prefix * "MMS_errors.txt", "w") do io @@ -6353,14 +6411,20 @@ function manufactured_solutions_analysis(run_info; plot_prefix) ("Er", L"\tilde{E}_r", L"\tilde{E}_r^{sym}", L"\varepsilon(\tilde{E}_r)"), ("Ez", L"\tilde{E}_z", L"\tilde{E}_z^{sym}", L"\varepsilon(\tilde{E}_z)"), ("density", L"\tilde{n}_i", L"\tilde{n}_i^{sym}", L"\varepsilon(\tilde{n}_i)"), + ("parallel_flow", L"\tilde{u}_{i,\parallel}", L"\tilde{u}_{i,\parallel}^{sym}", L"\varepsilon(\tilde{u}_{i,\parallel})"), + ("parallel_pressure", L"\tilde{p}_{i,\parallel}", L"\tilde{p}_{i,\parallel}^{sym}", L"\varepsilon(\tilde{p}_{i,\parallel})"), ("density_neutral", L"\tilde{n}_n", L"\tilde{n}_n^{sym}", L"\varepsilon(\tilde{n}_n)")) if contains(variable_name, "neutral") && run_info.n_neutral_species == 0 continue end + if contains(variable_name, "Er") && run_info.r.n_global == 1 + continue + end compare_moment_symbolic_test(run_info, plot_prefix, field_label, field_sym_label, - norm_label, variable_name; io=io, input=input) + norm_label, variable_name; io=io, input=input, + nvperp=nvperp) end end diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 8f93faf2b..df267400e 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -24,10 +24,12 @@ using IfElse if dfni_vpa_power_opt == "2" pvpa = 2 nconst = 0.25 + pconst = 3.0/4.0 fluxconst = 0.5 elseif dfni_vpa_power_opt == "4" pvpa = 4 - nconst = (3.0/8.0) + nconst = 3.0/8.0 + pconst = 15.0/8.0 fluxconst = 1.0 end @@ -199,9 +201,74 @@ using IfElse * "manufactured_solns:type=$(manufactured_solns_input.type)") end end - - function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition, - manufactured_solns_input) + + # ion mean parallel flow symbolic function + function upari_sym(Lr,Lz,r_bc,z_bc,composition,geometry,nr,manufactured_solns_input,species) + if z_bc == "periodic" + upari = 0.0 #not supported + elseif z_bc == "wall" + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species) + rhostar = geometry.rhostar + bzed = geometry.bzed + epsilon = manufactured_solns_input.epsilon_offset + alpha = manufactured_solns_input.alpha_switch + upari = ( (fluxconst/(sqrt(pi)*densi))*((z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) + - (0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha)) + + alpha*(rhostar/(2.0*bzed))*Er ) + end + return upari + end + + # ion parallel pressure symbolic function + function ppari_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + # normalisation factor due to strange pressure normalisation convention in master + norm_fac = 0.5 + if z_bc == "periodic" + ppari = 0.0 # not supported + elseif z_bc == "wall" + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + epsilon = manufactured_solns_input.epsilon_offset + alpha = manufactured_solns_input.alpha_switch + ppari = ( pconst*((0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) + + (z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha)) + + (z/Lz + 0.5)*(0.5 - z/Lz)*nzero_sym(Lr,Lz,r_bc,z_bc,alpha) + - (2.0/(pi*densi))*((z/Lz + 0.5)*nplus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha) + - (0.5 - z/Lz)*nminus_sym(Lr,Lz,r_bc,z_bc,epsilon,alpha))^2 ) + end + return ppari*norm_fac + end + + # ion perpendicular pressure symbolic function + function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + normfac = 0.5 # if pressure normalised to 0.5* nref * Tref = mref cref^2 + #normfac = 1.0 # if pressure normalised to nref*Tref + if nvperp > 1 + pperpi = densi # simple vperp^2 dependence of dfni + else + pperpi = 0.0 # marginalised model has nvperp = 1, vperp[1] = 0 + end + return pperpi*normfac + end + + # ion thermal speed symbolic function + function vthi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + ppari = ppari_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + pperpi = pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species,nvperp) + isotropic_pressure = (1.0/3.0)*(ppari + 2.0*pperpi) + normfac = 2.0 # if pressure normalised to 0.5* nref * Tref = mref cref^2 + #normfac = 1.0 # if pressure normalised to nref*Tref + if nvperp > 1 + vthi = sqrt(normfac*isotropic_pressure/densi) # thermal speed definition of 2V model + else + vthi = sqrt(normfac*ppari/densi) # thermal speed definition of 1V model + end + return vthi + end + + function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input) if z_bc == "periodic" jpari_into_LHS_wall_sym = 0.0 elseif z_bc == "wall" @@ -273,7 +340,7 @@ using IfElse # get N_e factor for boltzmann response if composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr == 1 # so 1D MMS test with 3V neutrals where ion current can be calculated prior to knowing Er - jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, composition, + jpari_into_LHS_wall = jpari_into_LHS_wall_sym(Lr, Lz, r_bc, z_bc, manufactured_solns_input) N_e = -2.0*sqrt(pi*composition.me_over_mi)*exp(-composition.phi_wall/composition.T_e)*jpari_into_LHS_wall elseif composition.electron_physics == boltzmann_electron_response_with_simple_sheath && nr > 1 @@ -308,7 +375,7 @@ using IfElse end function manufactured_solutions(manufactured_solns_input, Lr, Lz, r_bc, z_bc, - geometry, composition, species, nr) + geometry, composition, species, nr, nvperp) charged_species = species.charged[1] if composition.n_neutral_species > 0 neutral_species = species.neutral[1] @@ -318,6 +385,14 @@ using IfElse densi = densi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, charged_species) + upari = upari_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, + charged_species) + ppari = ppari_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species) + pperpi = pperpi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species, nvperp) + vthi = vthi_sym(Lr, Lz, r_bc, z_bc, composition, manufactured_solns_input, + charged_species, nvperp) dfni = dfni_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, charged_species) @@ -329,6 +404,10 @@ using IfElse #build julia functions from these symbolic expressions # cf. https://docs.juliahub.com/Symbolics/eABRO/3.4.0/tutorials/symbolic_functions/ densi_func = build_function(densi, z, r, t, expression=Val{false}) + upari_func = build_function(upari, z, r, t, expression=Val{false}) + ppari_func = build_function(ppari, z, r, t, expression=Val{false}) + pperpi_func = build_function(pperpi, z, r, t, expression=Val{false}) + vthi_func = build_function(vthi, z, r, t, expression=Val{false}) densn_func = build_function(densn, z, r, t, expression=Val{false}) dfni_func = build_function(dfni, vpa, vperp, z, r, t, expression=Val{false}) dfnn_func = build_function(dfnn, vz, vr, vzeta, z, r, t, expression=Val{false}) @@ -339,7 +418,10 @@ using IfElse # densn_func(zval, rval, tval) # dfnn_func(vzval, vrval, vzetapval, zval, rval, tval) - manufactured_solns_list = (densi_func = densi_func, densn_func = densn_func, dfni_func = dfni_func, dfnn_func = dfnn_func) + manufactured_solns_list = (densi_func = densi_func, densn_func = densn_func, + dfni_func = dfni_func, dfnn_func = dfnn_func, + upari_func = upari_func, ppari_func = ppari_func, + pperpi_func = pperpi_func, vthi_func = vthi_func) return manufactured_solns_list end @@ -374,6 +456,9 @@ using IfElse # ion manufactured solutions densi = densi_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, manufactured_solns_input, charged_species) + upari = upari_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, geometry, r_coord.n, manufactured_solns_input, charged_species) + vthi = vthi_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, manufactured_solns_input, + charged_species, vperp_coord.n) dfni = dfni_sym(r_coord.L, z_coord.L, r_coord.bc, z_coord.bc, composition, geometry, r_coord.n, manufactured_solns_input, charged_species) #dfni in vr vz vzeta coordinates @@ -432,11 +517,28 @@ using IfElse Si += - num_diss_params.vpa_dissipation_coefficient*Dvpa(Dvpa(dfni)) end if num_diss_params.r_dissipation_coefficient > 0.0 && include_num_diss_in_MMS - Si += - num_diss_params.r_dissipation_coefficient*Dr(Dr(dfni)) + Si += - rfac*num_diss_params.r_dissipation_coefficient*Dr(Dr(dfni)) end if num_diss_params.z_dissipation_coefficient > 0.0 && include_num_diss_in_MMS Si += - num_diss_params.z_dissipation_coefficient*Dz(Dz(dfni)) end + nu_krook = collisions.krook_collision_frequency_prefactor + if nu_krook > 0.0 + Ti_over_Tref = vthi^2 + if collisions.krook_collisions_option == "manual" + nuii_krook = nu_krook + else # default option + nuii_krook = nu_krook * densi * Ti_over_Tref^(-1.5) + end + if vperp_coord.n > 1 + pvth = 3 + else + pvth = 1 + end + FMaxwellian = (densi/vthi^pvth)*exp( -( ( vpa-upari)^2 + vperp^2 )/vthi^2) + Si += -nuii_krook*(FMaxwellian - dfni) + end + Source_i = expand_derivatives(Si) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 880284a0f..0611eba69 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -11,6 +11,7 @@ using MPI # be defined include("../machines/shared/machine_setup.jl") # Included so Documenter.jl can find its docs include("command_line_options.jl") +include("constants.jl") include("debugging.jl") include("type_definitions.jl") include("communication.jl") @@ -26,6 +27,7 @@ include("quadrature.jl") include("hermite_spline_interpolation.jl") include("derivatives.jl") include("input_structs.jl") +include("reference_parameters.jl") include("coordinates.jl") include("file_io.jl") include("velocity_moments.jl") @@ -45,6 +47,7 @@ include("neutral_z_advection.jl") include("neutral_vz_advection.jl") include("charge_exchange.jl") include("ionization.jl") +include("krook_collisions.jl") include("continuity.jl") include("energy_equation.jl") include("force_balance.jl") @@ -61,6 +64,7 @@ include("plot_sequence.jl") include("makie_post_processing.jl") include("time_advance.jl") +include("utils.jl") using TimerOutputs using Dates using Glob @@ -71,7 +75,6 @@ using .file_io: write_moments_data_to_binary, write_dfns_data_to_binary using .command_line_options: get_options using .communication using .communication: _block_synchronize -using .coordinates: define_coordinate using .debugging using .input_structs using .initial_conditions: allocate_pdf_and_moments, init_pdf_and_moments!, @@ -266,28 +269,11 @@ function setup_moment_kinetics(input_dict::Dict; input = mk_input(input_dict; save_inputs_to_txt=true, ignore_MPI=false) # obtain input options from moment_kinetics_input.jl # and check input to catch errors - io_input, evolve_moments, - t_input, z_input, r_input, - vpa_input, vperp_input, gyrophase_input, - vz_input, vr_input, vzeta_input, - composition, species, collisions, - geometry, drive_input, num_diss_params, manufactured_solns_input = input - # initialize z grid and write grid point locations to file - z, z_spectral = define_coordinate(z_input, io_input.parallel_io) - # initialize r grid and write grid point locations to file - r, r_spectral = define_coordinate(r_input, io_input.parallel_io) - # initialize vpa grid and write grid point locations to file - vpa, vpa_spectral = define_coordinate(vpa_input, io_input.parallel_io) - # initialize vperp grid and write grid point locations to file - vperp, vperp_spectral = define_coordinate(vperp_input, io_input.parallel_io) - # initialize gyrophase grid and write grid point locations to file - gyrophase, gyrophase_spectral = define_coordinate(gyrophase_input, io_input.parallel_io) - # initialize vz grid and write grid point locations to file - vz, vz_spectral = define_coordinate(vz_input, io_input.parallel_io) - # initialize vr grid and write grid point locations to file - vr, vr_spectral = define_coordinate(vr_input, io_input.parallel_io) - # initialize vr grid and write grid point locations to file - vzeta, vzeta_spectral = define_coordinate(vzeta_input, io_input.parallel_io) + io_input, evolve_moments, t_input, z, z_spectral, r, r_spectral, vpa, vpa_spectral, + vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, + vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, + drive_input, num_diss_params, manufactured_solns_input = input + # Create loop range variables for shared-memory-parallel loops if debug_loop_type === nothing # Non-debug case used for all simulations diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index a7e4ab048..72e23d6f5 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -10,10 +10,13 @@ export read_input_file using ..type_definitions: mk_float, mk_int using ..array_allocation: allocate_float using ..communication +using ..coordinates: define_coordinate using ..file_io: io_has_parallel, input_option_error, open_ascii_output_file +using ..krook_collisions: setup_krook_collisions using ..finite_differences: fd_check_option using ..input_structs using ..numerical_dissipation: setup_numerical_dissipation +using ..reference_parameters using MPI using TOML @@ -35,114 +38,6 @@ function read_input_file(input_filename::String) return input end -import Base: get -""" -Utility method for converting a string to an Enum when getting from a Dict, based on the -type of the default value -""" -function get(d::Dict, key, default::Enum) - valstring = get(d, key, nothing) - if valstring == nothing - return default - # instances(typeof(default)) gets the possible values of the Enum. Then convert to - # Symbol, then to String. - elseif valstring ∈ (split(s, ".")[end] for s ∈ String.(Symbol.(instances(typeof(default))))) - return eval(Symbol(valstring)) - else - error("Expected a $(typeof(default)), but '$valstring' is not in " - * "$(instances(typeof(default)))") - end -end - -""" -Set the defaults for options in the top level of the input, and check that there are not -any unexpected options (i.e. options that have no default). - -Modifies the options[section_name]::Dict by adding defaults for any values that are not -already present. - -Ignores any sections, as these will be checked separately. -""" -function set_defaults_and_check_top_level!(options::AbstractDict; kwargs...) - DictType = typeof(options) - - # Check for any unexpected values in the options - all options that are set should be - # present in the kwargs of this function call - options_keys_symbols = keys(kwargs) - options_keys = (String(k) for k ∈ options_keys_symbols) - for (key, value) in options - # Ignore any ssections when checking - if !(isa(value, AbstractDict) || key ∈ options_keys) - error("Unexpected option '$key=$value' in top-level options") - end - end - - # Set default values if a key was not set explicitly - explicit_keys = keys(options) - for (key_sym, value) ∈ kwargs - key = String(key_sym) - if !(key ∈ explicit_keys) - options[key] = value - end - end - - return options -end - -""" -Set the defaults for options in a section, and check that there are not any unexpected -options (i.e. options that have no default). - -Modifies the options[section_name]::Dict by adding defaults for any values that are not -already present. -""" -function set_defaults_and_check_section!(options::AbstractDict, section_name; - kwargs...) - DictType = typeof(options) - - if !(section_name ∈ keys(options)) - # If section is not present, create it - options[section_name] = DictType() - end - - if !isa(options[section_name], AbstractDict) - error("Expected '$section_name' to be a section in the input file, but it has a " - * "value '$(options[section_name])'") - end - - section = options[section_name] - - # Check for any unexpected values in the section - all options that are set should be - # present in the kwargs of this function call - section_keys_symbols = keys(kwargs) - section_keys = (String(k) for k ∈ section_keys_symbols) - for (key, value) in section - if !(key ∈ section_keys) - error("Unexpected option '$key=$value' in section '$section_name'") - end - end - - # Set default values if a key was not set explicitly - explicit_keys = keys(section) - for (key_sym, value) ∈ kwargs - key = String(key_sym) - if !(key ∈ explicit_keys) - section[key] = value - end - end - - return section -end - -""" -Convert a Dict whose keys are String or Symbol to a NamedTuple - -Useful as NamedTuple is immutable, so option values cannot be accidentally changed. -""" -function Dict_to_NamedTuple(d) - return NamedTuple(Symbol(k)=>v for (k,v) ∈ d) -end - """ Process user-supplied inputs @@ -196,28 +91,17 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # constant to be used to test nonzero Er in wall boundary condition composition.Er_constant = get(scan_input, "Er_constant", 0.0) - # Get reference parameters for normalizations - reference_parameter_section = set_defaults_and_check_section!( - scan_input, "reference_params"; - Bref=1.0, - Lref=10.0, - Nref=1.0e19, - Tref=100.0, - ) - reference_parameters = Dict_to_NamedTuple(reference_parameter_section) + # Reference parameters that define the conversion between physical quantities and + # normalised values used in the code. + reference_params = setup_reference_parameters(scan_input) - elementary_charge = 1.602176634e-19 # C - mi = 3.3435837724e-27 # kg - cref = sqrt(2.0 * elementary_charge*reference_parameters.Tref / mi) # m/s - Omegaref = elementary_charge * reference_parameters.Bref / mi - ## set geometry_input geometry.Bzed = get(scan_input, "Bzed", 1.0) geometry.Bmag = get(scan_input, "Bmag", 1.0) geometry.bzed = geometry.Bzed/geometry.Bmag geometry.bzeta = sqrt(1.0 - geometry.bzed^2.0) geometry.Bzeta = geometry.Bmag*geometry.bzeta - geometry.rhostar = get(scan_input, "rhostar", cref/reference_parameters.Lref/Omegaref) + geometry.rhostar = get(scan_input, "rhostar", get_default_rhostar(reference_params)) #println("Info: Bzed is ",geometry.Bzed) #println("Info: Bmag is ",geometry.Bmag) #println("Info: rhostar is ",geometry.rhostar) @@ -284,6 +168,19 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.charge_exchange = get(scan_input, "charge_exchange_frequency", 2.0*sqrt(species.charged[1].initial_temperature)) collisions.ionization = get(scan_input, "ionization_frequency", collisions.charge_exchange) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) + collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "none") + nuii_krook_default = setup_krook_collisions(reference_params) + if collisions.krook_collisions_option == "reference_parameters" + collisions.krook_collision_frequency_prefactor = nuii_krook_default + elseif collisions.krook_collisions_option == "manual" # get the frequency from the input file + collisions.krook_collision_frequency_prefactor = get(scan_input, "nuii_krook", nuii_krook_default) + elseif collisions.krook_collisions_option == "none" + # By default, no krook collisions included + collisions.krook_collision_frequency_prefactor = -1.0 + else + error("Invalid option " + * "krook_collisions_option=$(collisions.krook_collisions_option) passed") + end # parameters related to the time stepping nstep = get(scan_input, "nstep", 5) @@ -475,7 +372,8 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) # need grid and MPI information to determine these values # MRH just put dummy values now if ignore_MPI - irank_z = nrank_z = irank_r = nrank_r = -1 + irank_z = irank_r = 0 + nrank_z = nrank_r = 1 comm_sub_z = comm_sub_r = MPI.COMM_NULL else irank_z, nrank_z, comm_sub_z, irank_r, nrank_r, comm_sub_r = setup_distributed_memory_MPI(z.nelement_global,z.nelement_local,r.nelement_global,r.nelement_local) @@ -595,6 +493,23 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) io_immutable = io_input(; output_dir=output_dir, run_name=run_name, Dict(Symbol(k)=>v for (k,v) in io_settings)...) + # initialize z grid and write grid point locations to file + z, z_spectral = define_coordinate(z_immutable, io_immutable.parallel_io) + # initialize r grid and write grid point locations to file + r, r_spectral = define_coordinate(r_immutable, io_immutable.parallel_io) + # initialize vpa grid and write grid point locations to file + vpa, vpa_spectral = define_coordinate(vpa_immutable, io_immutable.parallel_io) + # initialize vperp grid and write grid point locations to file + vperp, vperp_spectral = define_coordinate(vperp_immutable, io_immutable.parallel_io) + # initialize gyrophase grid and write grid point locations to file + gyrophase, gyrophase_spectral = define_coordinate(gyrophase_immutable, io_immutable.parallel_io) + # initialize vz grid and write grid point locations to file + vz, vz_spectral = define_coordinate(vz_immutable, io_immutable.parallel_io) + # initialize vr grid and write grid point locations to file + vr, vr_spectral = define_coordinate(vr_immutable, io_immutable.parallel_io) + # initialize vr grid and write grid point locations to file + vzeta, vzeta_spectral = define_coordinate(vzeta_immutable, io_immutable.parallel_io) + if global_rank[] == 0 && save_inputs_to_txt # Make file to log some information about inputs into. # check to see if output_dir exists in the current directory @@ -611,9 +526,10 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) save_inputs_to_txt) # return immutable structs for z, vpa, species and composition - all_inputs = (io_immutable, evolve_moments, t_input, - z_immutable, r_immutable, vpa_immutable, vperp_immutable, gyrophase_immutable, vz_immutable, vr_immutable, vzeta_immutable, - composition, species_immutable, collisions, geometry, drive_immutable, + all_inputs = (io_immutable, evolve_moments, t_input, z, z_spectral, r, r_spectral, + vpa, vpa_spectral, vperp, vperp_spectral, gyrophase, gyrophase_spectral, + vz, vz_spectral, vr, vr_spectral, vzeta, vzeta_spectral, composition, + species_immutable, collisions, geometry, drive_immutable, num_diss_params, manufactured_solns_input) println(io, "\nAll inputs returned from mk_input():") println(io, all_inputs) @@ -940,20 +856,13 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) phi_wall = 0.0 # constant to test nonzero Er Er_constant = 0.0 - # constant to control Ez divergence - epsilon_offset = 0.001 - # bool to control functional form of dfni in MMS test - use_vpabar_in_mms_dfni = true - # float to control form of MMS density/potential/Er/Ez - alpha_switch = 1.0 # ratio of the neutral particle mass to the ion particle mass mn_over_mi = 1.0 # ratio of the electron particle mass to the ion particle mass me_over_mi = 1.0/1836.0 composition = species_composition(n_species, n_ion_species, n_neutral_species, electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant, - epsilon_offset, use_vpabar_in_mms_dfni, alpha_switch, mn_over_mi, me_over_mi, - allocate_float(n_species)) + mn_over_mi, me_over_mi, allocate_float(n_species)) species_charged = Array{species_parameters_mutable,1}(undef,n_ion_species) species_neutral = Array{species_parameters_mutable,1}(undef,n_neutral_species) @@ -1051,7 +960,9 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) # ionization collision frequency ionization = 0.0 constant_ionization_rate = false - collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate) + krook_collision_frequency_prefactor = -1.0 + collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, + krook_collision_frequency_prefactor,"none") Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength @@ -1218,4 +1129,13 @@ function check_input_initialization(composition, species, io) end end +""" + function get_default_rhostar(reference_params) + +Calculate the normalised ion gyroradius at reference parameters +""" +function get_default_rhostar(reference_params) + return reference_params.cref / reference_params.Omegaref / reference_params.Lref +end + end diff --git a/src/moment_kinetics_structs.jl b/src/moment_kinetics_structs.jl index fa63cd9d0..949250317 100644 --- a/src/moment_kinetics_structs.jl +++ b/src/moment_kinetics_structs.jl @@ -16,6 +16,7 @@ struct scratch_pdf{n_distribution_charged, n_moment, n_distribution_neutral,n_mo density::MPISharedArray{mk_float, n_moment} upar::MPISharedArray{mk_float, n_moment} ppar::MPISharedArray{mk_float, n_moment} + pperp::MPISharedArray{mk_float, n_moment} temp_z_s::MPISharedArray{mk_float, n_moment} # neutral particles pdf_neutral::MPISharedArray{mk_float, n_distribution_neutral} diff --git a/src/plot_MMS_sequence.jl b/src/plot_MMS_sequence.jl index a95fc20d5..4d48e98a8 100644 --- a/src/plot_MMS_sequence.jl +++ b/src/plot_MMS_sequence.jl @@ -63,11 +63,10 @@ function get_MMS_error_data(path_list,scan_type,scan_name) scan_input = load_input(fid) # get run-time input/composition/geometry/collisions/species info for convenience - #run_name_internal, output_dir, evolve_moments, - # t_input, z_input, r_input, - # vpa_input, vperp_input, gyrophase_input, - # vz_input, vr_input, vzeta_input, - # composition, species, collisions, geometry, drive_input = mk_input(scan_input) + #io_input, evolve_moments, t_input, z, z_spectral, r, r_spectral, vpa, vpa_spectral, + # vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, + # vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, + # drive_input, num_diss_params, manufactured_solns_input = mk_input(scan_input) z_nelement, r_nelement, vpa_nelement, vperp_nelement, vz_nelement, vr_nelement, vzeta_nelement = get_coords_nelement(scan_input) if scan_type == "vpa_nelement" diff --git a/src/post_processing.jl b/src/post_processing.jl index a664199b0..b8ec9f5d9 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -223,9 +223,10 @@ function allocate_global_zr_charged_moments(nz_global,nr_global,n_ion_species,nt density = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_flow = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_pressure = allocate_float(nz_global,nr_global,n_ion_species,ntime) + perpendicular_pressure = allocate_float(nz_global,nr_global,n_ion_species,ntime) parallel_heat_flux = allocate_float(nz_global,nr_global,n_ion_species,ntime) thermal_speed = allocate_float(nz_global,nr_global,n_ion_species,ntime) - return density, parallel_flow, parallel_pressure, parallel_heat_flux, thermal_speed + return density, parallel_flow, parallel_pressure, perpendicular_pressure, parallel_heat_flux, thermal_speed end function allocate_global_zr_charged_dfns(nvpa_global, nvperp_global, nz_global, nr_global, @@ -311,8 +312,7 @@ function get_geometry_and_composition(scan_input,n_ion_species,n_neutral_species me_over_mi = 1.0/1836.0 composition = species_composition(n_species, n_ion_species, n_neutral_species, electron_physics, use_test_neutral_wall_pdf, T_e, T_wall, phi_wall, Er_constant, - epsilon_offset, use_vpabar_in_mms_dfni, alpha_switch, mn_over_mi, me_over_mi, - allocate_float(n_species)) + mn_over_mi, me_over_mi, allocate_float(n_species)) return geometry, composition end @@ -489,7 +489,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) Tuple(this_z.n_global for this_z ∈ z), Tuple(this_r.n_global for this_r ∈ r), ntime) - density, parallel_flow, parallel_pressure, parallel_heat_flux, thermal_speed = + density, parallel_flow, parallel_pressure, perpendicular_pressure, parallel_heat_flux, thermal_speed = get_tuple_of_return_values(allocate_global_zr_charged_moments, Tuple(this_z.n_global for this_z ∈ z), Tuple(this_r.n_global for this_r ∈ r), @@ -532,6 +532,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) "parallel_pressure", run_names, "moments", nblocks, Tuple(this_z.n for this_z ∈ z), Tuple(this_r.n for this_r ∈ r), iskip) + get_tuple_of_return_values(read_distributed_zr_data!, perpendicular_pressure, + "perpendicular_pressure", run_names, "moments", nblocks, + Tuple(this_z.n for this_z ∈ z), + Tuple(this_r.n for this_r ∈ r), iskip) get_tuple_of_return_values(read_distributed_zr_data!, parallel_heat_flux, "parallel_heat_flux", run_names, "moments", nblocks, Tuple(this_z.n for this_z ∈ z), @@ -1013,6 +1017,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) density = density[1] parallel_flow = parallel_flow[1] parallel_pressure = parallel_pressure[1] + perpendicular_pressure = perpendicular_pressure[1] parallel_heat_flux = parallel_heat_flux[1] thermal_speed = thermal_speed[1] time = time[1] @@ -1047,12 +1052,9 @@ function analyze_and_plot_data(prefix...; run_index=nothing) input = mk_input(scan_input) # obtain input options from moment_kinetics_input.jl # and check input to catch errors - io_input, evolve_moments, - t_input, z_input, r_input, - vpa_input, vperp_input, gyrophase_input, - vz_input, vr_input, vzeta_input, - composition, species, collisions, - geometry, drive_input, num_diss_params, manufactured_solns_input = input + io_input, evolve_moments, t_input, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, + composition, species, collisions, geometry, drive_input, num_diss_params, + manufactured_solns_input = input if !is_1D1V # make plots and animations of the phi, Ez and Er @@ -1099,9 +1101,13 @@ function analyze_and_plot_data(prefix...; run_index=nothing) manufactured_solns_list = manufactured_solutions(manufactured_solns_input, Lr_in, z_global.L, r_global.bc, z_global.bc, geometry, - composition, species, r_global.n) + composition, species, r_global.n, vperp.n) dfni_func = manufactured_solns_list.dfni_func densi_func = manufactured_solns_list.densi_func + upari_func = manufactured_solns_list.upari_func + ppari_func = manufactured_solns_list.ppari_func + pperpi_func = manufactured_solns_list.pperpi_func + vthi_func = manufactured_solns_list.vthi_func dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func manufactured_E_fields = @@ -1138,16 +1144,32 @@ function analyze_and_plot_data(prefix...; run_index=nothing) # ion test density_sym = copy(density[:,:,:,:]) + upar_sym = copy(density[:,:,:,:]) + ppar_sym = copy(density[:,:,:,:]) + pperp_sym = copy(density[:,:,:,:]) + vthi_sym = copy(density[:,:,:,:]) is = 1 for it in 1:ntime for ir in 1:r_global.n for iz in 1:z_global.n density_sym[iz,ir,is,it] = densi_func(z_global.grid[iz],r_global.grid[ir],time[it]) + upar_sym[iz,ir,is,it] = upari_func(z_global.grid[iz],r_global.grid[ir],time[it]) + ppar_sym[iz,ir,is,it] = ppari_func(z_global.grid[iz],r_global.grid[ir],time[it]) + pperp_sym[iz,ir,is,it] = pperpi_func(z_global.grid[iz],r_global.grid[ir],time[it]) + vthi_sym[iz,ir,is,it] = vthi_func(z_global.grid[iz],r_global.grid[ir],time[it]) end end end compare_moments_symbolic_test(run_name_label,density,density_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, L"\widetilde{n}_i",L"\widetilde{n}_i^{sym}",L"\varepsilon(\widetilde{n}_i)","dens") + compare_moments_symbolic_test(run_name_label,parallel_flow,upar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{u}_{\|\|i}",L"\widetilde{u}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{u}_{\|\|i})","upar") + compare_moments_symbolic_test(run_name_label,parallel_pressure,ppar_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{p}_{\|\|i}",L"\widetilde{p}_{\|\|i}^{sym}",L"\varepsilon(\widetilde{p}_{\|\|i})","ppar") + compare_moments_symbolic_test(run_name_label,perpendicular_pressure,pperp_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{p}_{\perp i}",L"\widetilde{p}_{\perp i}^{sym}",L"\varepsilon(\widetilde{p}_{\perp i})","pperp") + compare_moments_symbolic_test(run_name_label,thermal_speed,vthi_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + L"\widetilde{v}_{th,i}",L"\widetilde{v}_{th,i}^{sym}",L"\varepsilon(\widetilde{v}_{th,i})","vthi") compare_charged_pdf_symbolic_test(run_name_label,manufactured_solns_list,"ion", L"\widetilde{f}_i",L"\widetilde{f}^{sym}_i",L"\varepsilon(\widetilde{f}_i)","pdf") diff --git a/src/reference_parameters.jl b/src/reference_parameters.jl new file mode 100644 index 000000000..efdbb4ded --- /dev/null +++ b/src/reference_parameters.jl @@ -0,0 +1,36 @@ +""" +Reference parameters + +Reference parameters are not needed or used by the main part of the code, but define the +physical units of the simulation, and are needed for a few specific steps during setup +(e.g. calculation of normalised collision frequency). +""" +module reference_parameters + +export setup_reference_parameters + +using ..constants +using ..input_structs + +""" +""" +function setup_reference_parameters(input_dict) + # Get reference parameters for normalizations + reference_parameter_section = copy(set_defaults_and_check_section!( + input_dict, "reference_params"; + Bref=1.0, + Lref=10.0, + Nref=1.0e19, + Tref=100.0, + mref=deuteron_mass, + )) + reference_parameter_section["cref"] = sqrt(2.0 * proton_charge * reference_parameter_section["Tref"] / (reference_parameter_section["mref"])) + reference_parameter_section["timeref"] = reference_parameter_section["Lref"] / reference_parameter_section["cref"] + reference_parameter_section["Omegaref"] = proton_charge * reference_parameter_section["Bref"] / reference_parameter_section["mref"] + + reference_params = Dict_to_NamedTuple(reference_parameter_section) + + return reference_params +end + +end diff --git a/src/time_advance.jl b/src/time_advance.jl index 350d80f25..e7b276367 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -14,9 +14,8 @@ using ..debugging using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_data_to_binary, debug_dump using ..looping using ..moment_kinetics_structs: scratch_pdf -using ..moment_kinetics_structs: scratch_pdf using ..velocity_moments: update_moments!, update_moments_neutral!, reset_moments_status! -using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_qpar! +using ..velocity_moments: update_density!, update_upar!, update_ppar!, update_pperp!, update_qpar!, update_vth! using ..velocity_moments: update_neutral_density!, update_neutral_qz! using ..velocity_moments: update_neutral_uzeta!, update_neutral_uz!, update_neutral_ur! using ..velocity_moments: update_neutral_pzeta!, update_neutral_pz!, update_neutral_pr! @@ -40,6 +39,7 @@ using ..vperp_advection: update_speed_vperp!, vperp_advection! using ..vpa_advection: update_speed_vpa!, vpa_advection! using ..charge_exchange: charge_exchange_collisions_1V!, charge_exchange_collisions_3V! using ..ionization: ionization_collisions_1V!, ionization_collisions_3V!, constant_ionization_source! +using ..krook_collisions: krook_collisions! using ..numerical_dissipation: vpa_boundary_buffer_decay!, vpa_boundary_buffer_diffusion!, vpa_dissipation!, z_dissipation!, r_dissipation!, @@ -404,6 +404,7 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss advance_ionization = false advance_ionization_1V = false advance_ionization_source = false + advance_krook_collisions = false advance_numerical_dissipation = false advance_sources = false advance_continuity = false @@ -473,6 +474,9 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss if collisions.ionization > 0.0 && collisions.constant_ionization_rate advance_ionization_source = true end + if collisions.krook_collision_frequency_prefactor > 0.0 + advance_krook_collisions = true + end advance_numerical_dissipation = true # if evolving the density, must advance the continuity equation, # in addition to including sources arising from the use of a modified distribution @@ -521,12 +525,12 @@ function setup_advance_flags(moments, composition, t_input, collisions, num_diss advance_neutral_z_advection, advance_neutral_r_advection, advance_neutral_vz_advection, advance_cx, advance_cx_1V, advance_ionization, advance_ionization_1V, - advance_ionization_source, advance_numerical_dissipation, - advance_sources, advance_continuity, advance_force_balance, - advance_energy, advance_neutral_sources, - advance_neutral_continuity, advance_neutral_force_balance, - advance_neutral_energy, rk_coefs, manufactured_solns_test, - r_diffusion, vpa_diffusion, vz_diffusion) + advance_ionization_source, advance_krook_collisions, + advance_numerical_dissipation, advance_sources, + advance_continuity, advance_force_balance, advance_energy, + advance_neutral_sources, advance_neutral_continuity, + advance_neutral_force_balance, advance_neutral_energy, rk_coefs, + manufactured_solns_test, r_diffusion, vpa_diffusion, vz_diffusion) end function setup_dummy_and_buffer_arrays(nr,nz,nvpa,nvperp,nvz,nvr,nvzeta,nspecies_ion,nspecies_neutral) @@ -651,6 +655,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag density_array = allocate_shared_float(moment_dims...) upar_array = allocate_shared_float(moment_dims...) ppar_array = allocate_shared_float(moment_dims...) + pperp_array = allocate_shared_float(moment_dims...) temp_z_s_array = allocate_shared_float(moment_dims...) pdf_neutral_array = allocate_shared_float(pdf_neutral_dims...) @@ -660,7 +665,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag scratch[istage] = scratch_pdf(pdf_array, density_array, upar_array, - ppar_array, temp_z_s_array, + ppar_array, pperp_array, temp_z_s_array, pdf_neutral_array, density_neutral_array, uz_neutral_array, pz_neutral_array) @serial_region begin @@ -668,6 +673,7 @@ function setup_scratch_arrays(moments, pdf_charged_in, pdf_neutral_in, n_rk_stag scratch[istage].density .= moments.charged.dens scratch[istage].upar .= moments.charged.upar scratch[istage].ppar .= moments.charged.ppar + scratch[istage].pperp .= moments.charged.pperp scratch[istage].pdf_neutral .= pdf_neutral_in scratch[istage].density_neutral .= moments.neutral.dens @@ -1064,6 +1070,14 @@ function time_advance_split_operators!(pdf, scratch, t, t_input, vpa, z, advance.ionization_collisions = false end end + if collisions.krook_collision_frequency_prefactor > 0.0 + advance.krook_collisions = true + time_advance_no_splitting!(pdf, scratch, t, t_input, z, vpa, + z_spectral, vpa_spectral, moments, fields, z_advect, vpa_advect, + z_SL, vpa_SL, composition, collisions, sources, num_diss_params, + advance, istep) + advance.krook_collisions = false + end # and add the source terms associated with redefining g = pdf/density or pdf*vth/density # to the kinetic equation if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar @@ -1215,7 +1229,9 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v ## # update the charged particle distribution and moments ## - + # MRH here we seem to have duplicate arrays for storing n, u||, p||, etc, but not for vth + # MRH 'scratch' is for the multiple stages of time advanced quantities, but 'moments' can be updated directly at each stage + # MRH in the standard drift-kinetic model. Consider taking moment quantities out of scratch for clarity. @loop_s_r_z_vperp_vpa is ir iz ivperp ivpa begin new_scratch.pdf[ivpa,ivperp,iz,ir,is] = rk_coefs[1]*pdf.charged.norm[ivpa,ivperp,iz,ir,is] + rk_coefs[2]*old_scratch.pdf[ivpa,ivperp,iz,ir,is] + rk_coefs[3]*new_scratch.pdf[ivpa,ivperp,iz,ir,is] end @@ -1246,15 +1262,12 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v vpa) end end - # update remaining velocity moments that are calculable from the evolved pdf update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composition) # update the thermal speed begin_s_r_z_region() try #below block causes DomainError if ppar < 0 or density, so exit cleanly if possible - @loop_s_r_z is ir iz begin - moments.charged.vth[iz,ir,is] = sqrt(2.0*new_scratch.ppar[iz,ir,is]/new_scratch.density[iz,ir,is]) - end + update_vth!(moments.charged.vth, new_scratch.ppar, new_scratch.pperp, new_scratch.density, vperp, z, r, composition) catch e if global_size[] > 1 println("ERROR: error calculating vth in time_advance.jl") @@ -1418,7 +1431,8 @@ function update_derived_moments!(new_scratch, moments, vpa, vperp, z, r, composi update_ppar!(new_scratch.ppar, moments.charged.ppar_updated, new_scratch.density, new_scratch.upar, new_scratch.pdf, vpa, vperp, z, r, composition, moments.evolve_density, moments.evolve_upar) - end + end + update_pperp!(new_scratch.pperp, new_scratch.pdf, vpa, vperp, z, r, composition) end """ @@ -1463,6 +1477,7 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, first_scratch.density[iz,ir,is] = moments.charged.dens[iz,ir,is] first_scratch.upar[iz,ir,is] = moments.charged.upar[iz,ir,is] first_scratch.ppar[iz,ir,is] = moments.charged.ppar[iz,ir,is] + first_scratch.pperp[iz,ir,is] = moments.charged.pperp[iz,ir,is] end if composition.n_neutral_species > 0 @@ -1513,6 +1528,7 @@ function ssp_rk!(pdf, scratch, t, t_input, vz, vr, vzeta, vpa, vperp, gyrophase, moments.charged.dens[iz,ir,is] = final_scratch.density[iz,ir,is] moments.charged.upar[iz,ir,is] = final_scratch.upar[iz,ir,is] moments.charged.ppar[iz,ir,is] = final_scratch.ppar[iz,ir,is] + moments.charged.pperp[iz,ir,is] = final_scratch.pperp[iz,ir,is] end if composition.n_neutral_species > 0 # No need to synchronize here as we only change neutral quantities and previous @@ -1677,6 +1693,12 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, collisions, dt) end + # Add Krook collision operator for ions + if advance.krook_collisions + krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions, + vperp, vpa, dt) + end + # add numerical dissipation if advance.numerical_dissipation vpa_dissipation!(fvec_out.pdf, fvec_in.pdf, vpa, vpa_spectral, dt, diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 000000000..ad8e5f19e --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,99 @@ +""" +Utility functions +""" +module utils + +export get_unnormalized_parameters, print_unnormalized_parameters + +using ..constants +using ..moment_kinetics_input: mk_input +using ..reference_parameters + +using OrderedCollections +using TOML +using Unitful + +Unitful.@unit eV "eV" "electron volt" proton_charge*Unitful.J true + +function __init__() + Unitful.register(utils) +end + +""" + get_unnormalized_parameters(input::Dict) + get_unnormalized_parameters(input_filename::String) + +Get many parameters for the simulation setup given by `input` or in the file +`input_filename`, in SI units and eV, returned as an OrderedDict. +""" +function get_unnormalized_parameters end +function get_unnormalized_parameters(input::Dict) + io_input, evolve_moments, t_input, z, z_spectral, r, r_spectral, vpa, vpa_spectral, + vperp, vperp_spectral, gyrophase, gyrophase_spectral, vz, vz_spectral, vr, + vr_spectral, vzeta, vzeta_spectral, composition, species, collisions, geometry, + drive_input, external_source_settings, num_diss_params, manufactured_solns_input = + mk_input(input) + + reference_params = setup_reference_parameters(input) + + Nnorm = reference_params.Nref * Unitful.m^(-3) + Tnorm = reference_params.Tref * eV + Lnorm = reference_params.Lref * Unitful.m + Bnorm = reference_params.Bref * Unitful.T + cnorm = reference_params.cref * Unitful.m / Unitful.s + timenorm = reference_params.timeref * Unitful.s + + # Assume single ion species so normalised ion mass is always 1 + mi = reference_params.mref * Unitful.kg + + parameters = OrderedDict{String,Any}() + parameters["run_name"] = run_name + + parameters["Nnorm"] = Nnorm + parameters["Tnorm"] = Tnorm + parameters["Lnorm"] = Lnorm + + parameters["Lz"] = Lnorm * z_input.L + + parameters["cs0"] = cnorm + + dt = t_input.dt * timenorm + parameters["dt"] = dt + parameters["output time step"] = dt * t_input.nwrite + parameters["total simulated time"] = dt * t_input.nstep + + parameters["T_e"] = Tnorm * composition.T_e + parameters["T_wall"] = Tnorm * composition.T_wall + + parameters["CX_rate_coefficient"] = collisions.charge_exchange / Nnorm / timenorm + parameters["ionization_rate_coefficient"] = collisions.ionization / Nnorm / timenorm + parameters["coulomb_collision_frequency0"] = + collisions.coulomb_collision_frequency_prefactor / timenorm + + return parameters +end +function get_unnormalized_parameters(input_filename::String, args...; kwargs...) + return get_unnormalized_parameters(TOML.parsefile(input_filename), args...; + kwargs...) +end + +""" + print_unnormalized_parameters(input) + +Print many parameters for the simulation setup given by `input` (a Dict of parameters or +a String giving a filename), in SI units and eV. +""" +function print_unnormalized_parameters(args...; kwargs...) + + parameters = get_unnormalized_parameters(args...; kwargs...) + + println("Dimensional parameters for '$(parameters["run_name"])'") + + for (k,v) ∈ parameters + println("$k = $v") + end + + return nothing +end + +end #utils diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 8db08f9c3..60ef582e4 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -10,7 +10,9 @@ export update_moments! export update_density! export update_upar! export update_ppar! +export update_pperp! export update_qpar! +export update_vth! export reset_moments_status! export moments_chrg_substruct, moments_ntrl_substruct export update_neutral_density! @@ -22,6 +24,13 @@ export update_neutral_pr! export update_neutral_pzeta! export update_neutral_qz! +# for testing +export get_density +export get_upar +export get_ppar +export get_pperp +export get_pressure + using ..type_definitions: mk_float using ..array_allocation: allocate_shared_float, allocate_bool, allocate_float using ..calculus: integral @@ -58,6 +67,8 @@ struct moments_charged_substruct # sets/uses the value for the same subset of species. This means ppar_update does # not need to be a shared memory array. ppar_updated::Vector{Bool} + # this is the perpendicular pressure + pperp::MPISharedArray{mk_float,3} # this is the parallel heat flux qpar::MPISharedArray{mk_float,3} # flag that keeps track of whether or not qpar needs updating before use @@ -182,6 +193,8 @@ function create_moments_charged(nz, nr, n_species, evolve_density, evolve_upar, # allocate array of Bools that indicate if the parallel pressure is updated for each species parallel_pressure_updated = allocate_bool(n_species) parallel_pressure_updated .= false + # allocate array used for the perpendicular pressure + perpendicular_pressure = allocate_shared_float(nz, nr, n_species) # allocate array used for the parallel flow parallel_heat_flux = allocate_shared_float(nz, nr, n_species) # allocate array of Bools that indicate if the parallel flow is updated for each species @@ -247,7 +260,7 @@ function create_moments_charged(nz, nr, n_species, evolve_density, evolve_upar, # return struct containing arrays needed to update moments return moments_charged_substruct(density, density_updated, parallel_flow, - parallel_flow_updated, parallel_pressure, parallel_pressure_updated, + parallel_flow_updated, parallel_pressure, parallel_pressure_updated,perpendicular_pressure, parallel_heat_flux, parallel_heat_flux_updated, thermal_speed, v_norm_fac, ddens_dz_upwind, d2dens_dz2, dupar_dz, dupar_dz_upwind, d2upar_dz2, dppar_dz, dppar_dz_upwind, d2ppar_dz2, dqpar_dz, dvth_dz) @@ -381,10 +394,7 @@ function update_moments!(moments, ff, vpa, vperp, z, r, composition) moments.evolve_upar) moments.charged.ppar_updated[is] = true end - @loop_r_z ir iz begin - moments.charged.vth[iz,ir,is] = - sqrt(2*moments.charged.ppar[iz,ir,is]/moments.charged.dens[iz,ir,is]) - end + @views update_pperp_species!(moments.charged.pperp[:,:,is], ff[:,:,:,:,is], vpa, vperp, z, r) if moments.charged.qpar_updated[is] == false @views update_qpar_species!(moments.charged.qpar[:,:,is], moments.charged.dens[:,:,is], @@ -395,6 +405,7 @@ function update_moments!(moments, ff, vpa, vperp, z, r, composition) moments.charged.qpar_updated[is] = true end end + update_vth!(moments.charged.vth, moments.charged.ppar, moments.charged.pperp, moments.charged.dens, vperp, z, r, composition) return nothing end @@ -431,13 +442,16 @@ function update_density_species!(dens, ff, vpa, vperp, z, r) @loop_r_z ir iz begin # When evolve_density = false, the evolved pdf is the 'true' pdf, and the vpa # coordinate is (dz/dt) / c_s. - # Integrating calculates n_s / N_e = (1/√π)∫d(vpa/c_s) (√π f_s c_s / N_e) - dens[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) + dens[iz,ir] = get_density(@view(ff[:,:,iz,ir]), vpa, vperp) end return nothing end +function get_density(ff, vpa, vperp) + # Integrating calculates n_s / N_e = (0/√π)∫d(vpa/c_s) (√π f_s c_s / N_e) + return integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 0, vperp.wgts) +end + """ NB: if this function is called and if upar_updated is false, then the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density @@ -475,10 +489,10 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # Integrating calculates # (upar_s / vth_s) = (1/√π)∫d(vpa/vth_s) * (vpa/vth_s) * (√π f_s vth_s / n_s) # so convert from upar_s / vth_s to upar_s / c_s + # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0 @loop_r_z ir iz begin vth = sqrt(2.0*ppar[iz,ir]/density[iz,ir]) - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) * vth + upar[iz,ir] = vth*get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0) end elseif evolve_density # corresponds to case where only the density is evolved separately from the @@ -486,9 +500,9 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # (dz/dt) / c_s. # Integrating calculates # (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / n_s) + # we set the input density to get_upar = 1.0 as the normalised distribution has density of 1.0 @loop_r_z ir iz begin - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) + upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, 1.0) end else # When evolve_density = false, the evolved pdf is the 'true' pdf, @@ -496,14 +510,21 @@ function update_upar_species!(upar, density, ppar, ff, vpa, vperp, z, r, evolve_ # Integrating calculates # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e) @loop_r_z ir iz begin - upar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), - vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) / - density[iz,ir] + upar[iz,ir] = get_upar(@view(ff[:,:,iz,ir]), vpa, vperp, density[iz,ir]) end end return nothing end +function get_upar(ff, vpa, vperp, density) + # Integrating calculates + # (n_s / N_e) * (upar_s / c_s) = (1/√π)∫d(vpa/c_s) * (vpa/c_s) * (√π f_s c_s / N_e) + # so we divide by the density of f_s + upar = integrate_over_vspace(@view(ff[:,:]), vpa.grid, 1, vpa.wgts, vperp.grid, 0, vperp.wgts) + upar /= density + return upar +end + """ NB: if this function is called and if ppar_updated is false, then the incoming pdf is the un-normalized pdf that satisfies int dv pdf = density @@ -540,40 +561,102 @@ function update_ppar_species!(ppar, density, upar, ff, vpa, vperp, z, r, evolve_ if evolve_upar # this is the case where the parallel flow and density are evolved separately # from the normalized pdf, g_s = (√π f_s c_s / n_s); the vpa coordinate is - # ((dz/dt) - upar_s) / c_s> - # Integrating calculates (p_parallel/m_s n_s c_s^2) = (1/√π)∫d((vpa-upar_s)/c_s) (1/2)*((vpa-upar_s)/c_s)^2 * (√π f_s c_s / n_s) - # so convert from p_s / m_s n_s c_s^2 to ppar_s = p_s / m_s N_e c_s^2 + # ((dz/dt) - upar_s) / c_s> and so we set upar = 0 in the call to get_ppar + # because the mean flow of the normalised ff is zero @loop_r_z ir iz begin - ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) * - density[iz,ir] + ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, 0.0) end elseif evolve_density # corresponds to case where only the density is evolved separately from the # normalised pdf, given by g_s = (√π f_s c_s / n_s); the vpa coordinate is # (dz/dt) / c_s. - # Integrating calculates - # (p_parallel/m_s n_s c_s^2) + (upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / n_s) - # so subtract off the mean kinetic energy and multiply by density to get the - # internal energy density (aka pressure) @loop_r_z ir iz begin - ppar[iz,ir] = (integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) - - upar[iz,ir]^2) * density[iz,ir] + ppar[iz,ir] = density[iz,ir]*get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir]) end else # When evolve_density = false, the evolved pdf is the 'true' pdf, # and the vpa coordinate is (dz/dt) / c_s. - # Integrating calculates - # (p_parallel/m_s N_e c_s^2) + (n_s/N_e)*(upar_s/c_s)^2 = (1/√π)∫d(vpa/c_s) (vpa/c_s)^2 * (√π f_s c_s / N_e) - # so subtract off the mean kinetic energy density to get the internal energy - # density (aka pressure) @loop_r_z ir iz begin - ppar[iz,ir] = integrate_over_vspace(@view(ff[:,:,iz,ir]), vpa.grid, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) - - density[iz,ir]*upar[iz,ir]^2 + ppar[iz,ir] = get_ppar(@view(ff[:,:,iz,ir]), vpa, vperp, upar[iz,ir]) end end return nothing end +function get_ppar(ff, vpa, vperp, upar) + # Integrating calculates + # (p_parallel/m_s N_e c_s^2) = (1/√π)∫d(vpa/c_s) ((vpa-upar)/c_s)^2 * (√π f_s c_s / N_e) + # the internal energy density (aka pressure of f_s) + + # modify input vpa.grid to account for the mean flow + @. vpa.scratch = vpa.grid - upar + norm_fac = 1.0 # normalise to m_s N_e c_s^2 + #norm_fac = 2.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s + return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.scratch, 2, vpa.wgts, vperp.grid, 0, vperp.wgts) +end + +function update_pperp!(pperp, pdf, vpa, vperp, z, r, composition) + @boundscheck composition.n_ion_species == size(pperp,3) || throw(BoundsError(pperp)) + @boundscheck r.n == size(pperp,2) || throw(BoundsError(pperp)) + @boundscheck z.n == size(pperp,1) || throw(BoundsError(pperp)) + + begin_s_r_z_region() + + @loop_s is begin + @views update_pperp_species!(pperp[:,:,is], pdf[:,:,:,:,is], vpa, vperp, z, r) + end +end + +""" +calculate the updated perpendicular pressure (pperp) for a given species +""" +function update_pperp_species!(pperp, ff, vpa, vperp, z, r) + @boundscheck vpa.n == size(ff, 1) || throw(BoundsError(ff)) + @boundscheck vperp.n == size(ff, 2) || throw(BoundsError(ff)) + @boundscheck z.n == size(ff, 3) || throw(BoundsError(ff)) + @boundscheck r.n == size(ff, 4) || throw(BoundsError(ff)) + @boundscheck z.n == size(pperp, 1) || throw(BoundsError(pperp)) + @boundscheck r.n == size(pperp, 2) || throw(BoundsError(pperp)) + @loop_r_z ir iz begin + pperp[iz,ir] = get_pperp(@view(ff[:,:,iz,ir]), vpa, vperp) + end + return nothing +end + +function get_pperp(ff, vpa, vperp) + norm_fac = 0.5 # normalise to m_s N_e c_s^2 + #norm_fac = 1.0 # normalise to 0.5 m_s N_e c_s^2 = N_e T_s + return norm_fac*integrate_over_vspace(@view(ff[:,:]), vpa.grid, 0, vpa.wgts, vperp.grid, 2, vperp.wgts) +end + +function update_vth!(vth, ppar, pperp, dens, vperp, z, r, composition) + @boundscheck composition.n_ion_species == size(vth,3) || throw(BoundsError(vth)) + @boundscheck r.n == size(vth,2) || throw(BoundsError(vth)) + @boundscheck z.n == size(vth,1) || throw(BoundsError(vth)) + + begin_s_r_z_region() + normfac = 2.0 # if ppar normalised to 2*nref Tref = mref cref^2 + #normfac = 1.0 # if ppar normalised to nref Tref = 0.5 * mref cref^2 + if vperp.n > 1 #2V definition + @loop_s_r_z is ir iz begin + piso = get_pressure(ppar[iz,ir,is],pperp[iz,ir,is]) + vth[iz,ir,is] = sqrt(normfac*piso/dens[iz,ir,is]) + end + else #1V definition + @loop_s_r_z is ir iz begin + vth[iz,ir,is] = sqrt(normfac*ppar[iz,ir,is]/dens[iz,ir,is]) + end + end +end + +""" +compute the isotropic pressure from the already computed ppar and pperp +""" +function get_pressure(ppar::mk_float,pperp::mk_float) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres +end + """ NB: the incoming pdf is the normalized pdf """ diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl new file mode 100644 index 000000000..62e59c892 --- /dev/null +++ b/test/Krook_collisions_tests.jl @@ -0,0 +1,443 @@ +module KrookCollisionsTests + +# Test for Krook collision operator, based on NonlinearSoundWave test + +include("setup.jl") + +using Base.Filesystem: tempname +using MPI + +using moment_kinetics.coordinates: define_coordinate +using moment_kinetics.input_structs: grid_input, advection_input +using moment_kinetics.load_data: open_readonly_output_file, load_coordinate_data, + load_species_data, load_fields_data, + load_charged_particle_moments_data, load_pdf_data, + load_neutral_particle_moments_data, + load_neutral_pdf_data, load_time_data, load_species_data +using moment_kinetics.interpolation: interpolate_to_grid_z, interpolate_to_grid_vpa +using moment_kinetics.type_definitions: mk_float + +# Create a temporary directory for test output +test_output_directory = tempname() +mkpath(test_output_directory) + +# Useful parameters +const z_L = 1.0 # always 1 in normalized units? +const vpa_L = 8.0 + +# Use very small number of points in vpa_expected to reduce the amount of entries we +# need to store. First and last entries are within the grid (rather than at the ends) in +# order to get non-zero values. +# Note: in the arrays of numbers for expected data, space-separated entries have to stay +# on the same line. +const expected = + ( + z=[z for z in range(-0.5 * z_L, 0.5 * z_L, length=11)], + vpa=[vpa for vpa in range(-0.2 * vpa_L, 0.2 * vpa_L, length=3)], + phi=[-1.386282080324426 -1.2382381134436695; -1.2115129555832849 -1.1306145497660034; + -0.8609860698164498 -0.8726509017405432; -0.5494724768120176 -0.5904511161957423; + -0.35345976364887166 -0.37557956583926283; -0.28768207245186167 -0.2919214243915014; + -0.353459763648872 -0.3755795658392631; -0.5494724768120175 -0.5904511161957432; + -0.8609860698164502 -0.8726509017405427; -1.2115129555832849 -1.1306145497660032; + -1.3862820803244258 -1.2382381134436695], + n_charged=[0.2500030702177186 0.28989452952580286; 0.2977473631375158 0.32283464906590775; + 0.42274585818529853 0.41784232850006636; 0.5772542465450629 0.5540775204094593; + 0.7022542481909738 0.6868914177534788; 0.7499999999999394 0.7468272160708606; + 0.7022542481909738 0.6868914177534787; 0.577254246545063 0.554077520409459; + 0.42274585818529864 0.4178423285000665; 0.2977473631375159 0.3228346490659078; + 0.2500030702177185 0.2898945295258028], + n_neutral=[0.7499999999999382 0.7736770941648626; 0.7022542481909748 0.7056867075427516; + 0.5772542465450632 0.5582975660019874; 0.4227458581852985 0.4096913953484598; + 0.29774736313751604 0.3053964124252619; 0.2500030702177186 0.2681998023548167; + 0.29774736313751604 0.3053964124252619; 0.42274585818529836 0.4096913953484599; + 0.5772542465450631 0.5582975660019875; 0.7022542481909745 0.7056867075427524; + 0.7499999999999383 0.7736770941648626], + upar_charged=[-2.7135787559953277e-17 -1.6845791254993525e-16; -9.321028970172899e-18 -0.18245939812953485; + -2.8374879811351724e-18 -0.19666454846377826; 1.2124327390522635e-17 -0.11128043369942339; + 3.6525788403693063e-17 -0.03317985705380149; -2.0930856430671915e-17 4.720175801869314e-17; + 8.753545920086251e-18 0.033179857053801595; 1.1293771270243255e-17 0.11128043369942343; + 1.3739171132886587e-17 0.19666454846377784; -6.840453743089351e-18 0.18245939812953468; + -2.7135787559953277e-17 -1.9129596434811267e-16], + upar_neutral=[6.5569385065066925e-18 8.08747058038406e-18; 1.1054500872839027e-17 -0.03620988455458174; + -3.241833393685864e-17 -0.009156078199383568; -3.617637280460899e-17 0.05452623197292568; + 4.417578961284041e-17 0.07607875911384775; 4.9354467746194965e-17 1.635044638743921e-16; + 6.573091229872379e-18 -0.0760787591138477; 2.989662686945165e-17 -0.05452623197292564; + -3.1951996361666834e-17 0.009156078199383685; -4.395464518158184e-18 0.03620988455458165; + 6.5569385065066925e-18 1.8232586069007834e-18], + ppar_charged=[0.18749999999999992 0.23302732230115558; 0.20909325514551116 0.21936799130257528; + 0.24403180771238264 0.20856296024163393; 0.24403180771238278 0.2154266357557397; + 0.2090932551455113 0.2206183912107678; 0.1875 0.21979739387340663; + 0.20909325514551128 0.22061839121076784; 0.2440318077123828 0.21542663575573945; + 0.24403180771238256 0.20856296024163395; 0.20909325514551116 0.2193679913025754; + 0.18749999999999992 0.23302732230115553], + ppar_neutral=[0.18750000000000003 0.2480292382671593; 0.20909325514551122 0.24401255100297964; + 0.24403180771238286 0.22861763406831279; 0.24403180771238278 0.2058922545451891; + 0.20909325514551144 0.1926313699453636; 0.18749999999999992 0.19090651730415983; + 0.20909325514551141 0.19263136994536365; 0.2440318077123828 0.20589225454518903; + 0.24403180771238286 0.2286176340683127; 0.20909325514551114 0.24401255100297964; + 0.18750000000000006 0.24802923826715936], + f_charged=[0.0370462360994826 0.04059927063892091 0.0428431419871786 0.030398267195914062 0.01236045902698859 0.006338529470383425 0.012360459026988587 0.030398267195914028 0.04284314198717859 0.0405992706389209 0.0370462360994826; + 0.20411991941198782 0.25123395823993105 0.3934413727192304 0.6277900619432855 0.9100364506661008 1.0606601717796504 0.910036450666101 0.6277900619432859 0.39344137271923046 0.25123395823993094 0.20411991941198776; + 0.0370462360994826 0.04059927063892091 0.0428431419871786 0.030398267195914062 0.01236045902698859 0.006338529470383425 0.012360459026988587 0.030398267195914028 0.04284314198717859 0.0405992706389209 0.0370462360994826;;; + 0.0538996852594264 0.06066864433237418 0.03746866696438989 0.014783440166032301 0.010917691665145668 0.018422971878502774 0.027170953411068444 0.027269146560166702 0.026567569739750264 0.035612674100528624 0.05389968525942639; + 0.2118369019176154 0.24917436308523389 0.37345448114678914 0.5972219245577428 0.8859681860177208 1.0485988935814787 0.8859681860177204 0.5972219245577435 0.37345448114678825 0.24917436308523389 0.2118369019176155; + 0.05389968525942635 0.03561267410052869 0.02656756973975021 0.02726914656016675 0.027170953411068514 0.018422971878502753 0.01091769166514568 0.014783440166032254 0.037468666964389795 0.060668644332374164 0.05389968525942635], + f_neutral=[0.0063385294703834595 0.012360459026988546 0.030398267195914108 0.04284314198717859 0.040599270638920985 0.03704623609948259 0.040599270638920965 0.0428431419871786 0.030398267195914094 0.012360459026988546 0.006338529470383456; + 1.0606601717796493 0.9100364506661016 0.6277900619432857 0.3934413727192303 0.2512339582399308 0.20411991941198754 0.2512339582399307 0.3934413727192301 0.6277900619432853 0.9100364506661016 1.0606601717796487; + 0.0063385294703834595 0.012360459026988546 0.030398267195914108 0.04284314198717859 0.040599270638920985 0.03704623609948259 0.040599270638920965 0.0428431419871786 0.030398267195914094 0.012360459026988546 0.006338529470383456;;; + 0.0242848886629411 0.04071460358290305 0.04191389118981371 0.03638215764882266 0.03692283098105331 0.04164449216999481 0.03671950776850948 0.01928119243573099 0.008423252360063483 0.010011392733734206 0.02428488866294109; + 1.0530033604430462 0.9036809869030653 0.6251085339983469 0.3955308968816375 0.25710352416286547 0.21137159186144025 0.25710352416286547 0.3955308968816377 0.6251085339983473 0.9036809869030653 1.0530033604430464; + 0.024284888662941113 0.010011392733734206 0.008423252360063494 0.019281192435730943 0.036719507768509525 0.041644492169994836 0.03692283098105331 0.03638215764882269 0.04191389118981368 0.04071460358290303 0.024284888662941134]) + +# default inputs for tests +test_input_full_f = Dict("n_ion_species" => 1, + "n_neutral_species" => 1, + "boltzmann_electron_response" => true, + "run_name" => "full_f", + "base_directory" => test_output_directory, + "evolve_moments_density" => false, + "evolve_moments_parallel_flow" => false, + "evolve_moments_parallel_pressure" => false, + "evolve_moments_conservation" => true, + "krook_collisions_option" => "reference_parameters", + "T_e" => 1.0, + "initial_density1" => 0.5, + "initial_temperature1" => 1.0, + "initial_density2" => 0.5, + "initial_temperature2" => 1.0, + "z_IC_option1" => "sinusoid", + "z_IC_density_amplitude1" => 0.5, + "z_IC_density_phase1" => 0.0, + "z_IC_upar_amplitude1" => 0.0, + "z_IC_upar_phase1" => 0.0, + "z_IC_temperature_amplitude1" => 0.5, + "z_IC_temperature_phase1" => mk_float(π), + "z_IC_option2" => "sinusoid", + "z_IC_density_amplitude2" => 0.5, + "z_IC_density_phase2" => mk_float(π), + "z_IC_upar_amplitude2" => 0.0, + "z_IC_upar_phase2" => 0.0, + "z_IC_temperature_amplitude2" => 0.5, + "z_IC_temperature_phase2" => 0.0, + "charge_exchange_frequency" => 2*π*0.1, + "ionization_frequency" => 0.0, + "nstep" => 100, + "dt" => 0.001, + "nwrite" => 100, + "nwrite_dfns" => 100, + "use_semi_lagrange" => false, + "n_rk_stages" => 4, + "split_operators" => false, + "r_ngrid" => 1, + "r_nelement" => 1, + "r_bc" => "periodic", + "r_discretization" => "chebyshev_pseudospectral", + "z_ngrid" => 9, + "z_nelement" => 4, + "z_bc" => "periodic", + "z_discretization" => "chebyshev_pseudospectral", + "vpa_ngrid" => 17, + "vpa_nelement" => 8, + "vpa_L" => vpa_L, + "vpa_bc" => "periodic", + "vpa_discretization" => "chebyshev_pseudospectral", + "vz_ngrid" => 17, + "vz_nelement" => 8, + "vz_L" => vpa_L, + "vz_bc" => "periodic", + "vz_discretization" => "chebyshev_pseudospectral") + +test_input_split_1_moment = + merge(test_input_full_f, + Dict("run_name" => "split_1_moment", + "evolve_moments_density" => true)) + +test_input_split_2_moments = + merge(test_input_split_1_moment, + Dict("run_name" => "split_2_moments", + "evolve_moments_parallel_flow" => true)) + +test_input_split_3_moments = + merge(test_input_split_2_moments, + Dict("run_name" => "split_3_moments", + "evolve_moments_parallel_pressure" => true, + "vpa_L" => 12.0, "vz_L" => 12.0)) + + +""" +Run a sound-wave test for a single set of parameters +""" +# Note 'name' should not be shared by any two tests in this file +function run_test(test_input, rtol, atol; args...) + # by passing keyword arguments to run_test, args becomes a Dict which can be used to + # update the default inputs + + # Convert keyword arguments to a unique name + name = test_input["run_name"] + if length(args) > 0 + name = string(name, "_", (string(k, "-", v, "_") for (k, v) in args)...) + + # Remove trailing "_" + name = chop(name) + end + + # Provide some progress info + println(" - testing ", name) + + # Convert dict from symbol keys to String keys + modified_inputs = Dict(String(k) => v for (k, v) in args) + + # Update default inputs with values to be changed + input = merge(test_input, modified_inputs) + + input["run_name"] = name + + # Suppress console output while running + quietoutput() do + # run simulation + run_moment_kinetics(input) + end + + phi = nothing + n_charged = nothing + upar_charged = nothing + ppar_charged = nothing + f_charged = nothing + n_neutral = nothing + upar_neutral = nothing + ppar_neutral = nothing + f_neutral = nothing + z, z_spectral = nothing, nothing + vpa, vpa_spectral = nothing, nothing + + if global_rank[] == 0 + quietoutput() do + + # Load and analyse output + ######################### + + path = joinpath(realpath(input["base_directory"]), name, name) + + # open the netcdf file containing moments data and give it the handle 'fid' + fid = open_readonly_output_file(path, "moments") + + # load species, time coordinate data + n_ion_species, n_neutral_species = load_species_data(fid) + ntime, time = load_time_data(fid) + n_ion_species, n_neutral_species = load_species_data(fid) + + # load fields data + phi_zrt, Er_zrt, Ez_zrt = load_fields_data(fid) + + # load velocity moments data + n_charged_zrst, upar_charged_zrst, ppar_charged_zrst, qpar_charged_zrst, v_t_charged_zrst = load_charged_particle_moments_data(fid) + n_neutral_zrst, upar_neutral_zrst, ppar_neutral_zrst, qpar_neutral_zrst, v_t_neutral_zrst = load_neutral_particle_moments_data(fid) + z, z_spectral = load_coordinate_data(fid, "z") + + close(fid) + + # open the netcdf file containing pdf data + fid = open_readonly_output_file(path, "dfns") + + # load particle distribution function (pdf) data + f_charged_vpavperpzrst = load_pdf_data(fid) + f_neutral_vzvrvzetazrst = load_neutral_pdf_data(fid) + vpa, vpa_spectral = load_coordinate_data(fid, "vpa") + + close(fid) + + phi = phi_zrt[:,1,:] + n_charged = n_charged_zrst[:,1,:,:] + upar_charged = upar_charged_zrst[:,1,:,:] + ppar_charged = ppar_charged_zrst[:,1,:,:] + qpar_charged = qpar_charged_zrst[:,1,:,:] + v_t_charged = v_t_charged_zrst[:,1,:,:] + f_charged = f_charged_vpavperpzrst[:,1,:,1,:,:] + n_neutral = n_neutral_zrst[:,1,:,:] + upar_neutral = upar_neutral_zrst[:,1,:,:] + ppar_neutral = ppar_neutral_zrst[:,1,:,:] + qpar_neutral = qpar_neutral_zrst[:,1,:,:] + v_t_neutral = v_t_neutral_zrst[:,1,:,:] + f_neutral = f_neutral_vzvrvzetazrst[:,1,1,:,1,:,:] + + # Unnormalize f + if input["evolve_moments_density"] + for it ∈ 1:length(time), is ∈ 1:n_ion_species, iz ∈ 1:z.n + f_charged[:,iz,is,it] .*= n_charged[iz,is,it] + end + for it ∈ 1:length(time), isn ∈ 1:n_neutral_species, iz ∈ 1:z.n + f_neutral[:,iz,isn,it] .*= n_neutral[iz,isn,it] + end + end + if input["evolve_moments_parallel_pressure"] + for it ∈ 1:length(time), is ∈ 1:n_ion_species, iz ∈ 1:z.n + f_charged[:,iz,is,it] ./= v_t_charged[iz,is,it] + end + for it ∈ 1:length(time), isn ∈ 1:n_neutral_species, iz ∈ 1:z.n + f_neutral[:,iz,isn,it] ./= v_t_neutral[iz,isn,it] + end + end + end + + # Test against values interpolated onto 'expected' grid which is fairly coarse no we + # do not have to save too much data in this file + + # Use commented-out lines to get the test data to put in `expected` + #newgrid_phi = cat(interpolate_to_grid_z(expected.z, phi[:, 1], z, z_spectral), + # interpolate_to_grid_z(expected.z, phi[:, 2], z, z_spectral); + # dims=2) + #println("phi ", size(newgrid_phi)) + #println(newgrid_phi) + #println() + #newgrid_n_charged = cat(interpolate_to_grid_z(expected.z, n_charged[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, n_charged[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("n_charged ", size(newgrid_n_charged)) + #println(newgrid_n_charged) + #println() + #newgrid_n_neutral = cat(interpolate_to_grid_z(expected.z, n_neutral[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, n_neutral[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("n_neutral ", size(newgrid_n_neutral)) + #println(newgrid_n_neutral) + #println() + #newgrid_upar_charged = cat(interpolate_to_grid_z(expected.z, upar_charged[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, upar_charged[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("upar_charged ", size(newgrid_upar_charged)) + #println(newgrid_upar_charged) + #println() + #newgrid_upar_neutral = cat(interpolate_to_grid_z(expected.z, upar_neutral[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, upar_neutral[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("upar_neutral ", size(newgrid_upar_neutral)) + #println(newgrid_upar_neutral) + #println() + #newgrid_ppar_charged = cat(interpolate_to_grid_z(expected.z, ppar_charged[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, ppar_charged[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("ppar_charged ", size(newgrid_ppar_charged)) + #println(newgrid_ppar_charged) + #println() + #newgrid_ppar_neutral = cat(interpolate_to_grid_z(expected.z, ppar_neutral[:, :, 1], z, z_spectral)[:,1], + # interpolate_to_grid_z(expected.z, ppar_neutral[:, :, 2], z, z_spectral)[:,1]; + # dims=2) + #println("ppar_neutral ", size(newgrid_ppar_neutral)) + #println(newgrid_ppar_neutral) + #println() + #newgrid_f_charged = cat(interpolate_to_grid_vpa(expected.vpa, interpolate_to_grid_z(expected.z, f_charged[:, :, :, 1], z, z_spectral), vpa, vpa_spectral)[:,:,1], + # interpolate_to_grid_vpa(expected.vpa, interpolate_to_grid_z(expected.z, f_charged[:, :, :, 2], z, z_spectral), vpa, vpa_spectral)[:,:,1]; + # dims=4) + #println("f_charged ", size(newgrid_f_charged)) + #println(newgrid_f_charged) + #println() + #newgrid_f_neutral = cat(interpolate_to_grid_vpa(expected.vpa, interpolate_to_grid_z(expected.z, f_neutral[:, :, :, 1], z, z_spectral), vpa, vpa_spectral)[:,:,1], + # interpolate_to_grid_vpa(expected.vpa, interpolate_to_grid_z(expected.z, f_neutral[:, :, :, 2], z, z_spectral), vpa, vpa_spectral)[:,:,1]; + # dims=4) + #println("f_neutral ", size(newgrid_f_neutral)) + #println(newgrid_f_neutral) + #println() + function test_values(tind) + @testset "tind=$tind" begin + newgrid_phi = interpolate_to_grid_z(expected.z, phi[:, tind], z, z_spectral) + @test isapprox(expected.phi[:, tind], newgrid_phi, rtol=rtol) + + # Check charged particle moments and f + ###################################### + + newgrid_n_charged = interpolate_to_grid_z(expected.z, n_charged[:, :, tind], z, z_spectral) + @test isapprox(expected.n_charged[:, tind], newgrid_n_charged[:,1], rtol=rtol) + + newgrid_upar_charged = interpolate_to_grid_z(expected.z, upar_charged[:, :, tind], z, z_spectral) + @test isapprox(expected.upar_charged[:, tind], newgrid_upar_charged[:,1], rtol=rtol, atol=atol) + + newgrid_ppar_charged = interpolate_to_grid_z(expected.z, ppar_charged[:, :, tind], z, z_spectral) + @test isapprox(expected.ppar_charged[:, tind], newgrid_ppar_charged[:,1], rtol=rtol) + + newgrid_vth_charged = @. sqrt(2.0*newgrid_ppar_charged/newgrid_n_charged) + newgrid_f_charged = interpolate_to_grid_z(expected.z, f_charged[:, :, :, tind], z, z_spectral) + temp = newgrid_f_charged + newgrid_f_charged = fill(NaN, length(expected.vpa), + size(newgrid_f_charged, 2), + size(newgrid_f_charged, 3), + size(newgrid_f_charged, 4)) + for iz ∈ 1:length(expected.z) + wpa = copy(expected.vpa) + if input["evolve_moments_parallel_flow"] + wpa .-= newgrid_upar_charged[iz,1] + end + if input["evolve_moments_parallel_pressure"] + wpa ./= newgrid_vth_charged[iz,1] + end + newgrid_f_charged[:,iz,1] = interpolate_to_grid_vpa(wpa, temp[:,iz,1], vpa, vpa_spectral) + end + @test isapprox(expected.f_charged[:, :, tind], newgrid_f_charged[:,:,1], rtol=rtol) + + # Check neutral particle moments and f + ###################################### + + newgrid_n_neutral = interpolate_to_grid_z(expected.z, n_neutral[:, :, tind], z, z_spectral) + @test isapprox(expected.n_neutral[:, tind], newgrid_n_neutral[:,:,1], rtol=rtol) + + newgrid_upar_neutral = interpolate_to_grid_z(expected.z, upar_neutral[:, :, tind], z, z_spectral) + @test isapprox(expected.upar_neutral[:, tind], newgrid_upar_neutral[:,:,1], rtol=rtol, atol=atol) + + newgrid_ppar_neutral = interpolate_to_grid_z(expected.z, ppar_neutral[:, :, tind], z, z_spectral) + @test isapprox(expected.ppar_neutral[:, tind], newgrid_ppar_neutral[:,:,1], rtol=rtol) + + newgrid_vth_neutral = @. sqrt(2.0*newgrid_ppar_neutral/newgrid_n_neutral) + newgrid_f_neutral = interpolate_to_grid_z(expected.z, f_neutral[:, :, :, tind], z, z_spectral) + temp = newgrid_f_neutral + newgrid_f_neutral = fill(NaN, length(expected.vpa), + size(newgrid_f_neutral, 2), + size(newgrid_f_neutral, 3), + size(newgrid_f_neutral, 4)) + for iz ∈ 1:length(expected.z) + wpa = copy(expected.vpa) + if input["evolve_moments_parallel_flow"] + wpa .-= newgrid_upar_neutral[iz,1] + end + if input["evolve_moments_parallel_pressure"] + wpa ./= newgrid_vth_neutral[iz,1] + end + newgrid_f_neutral[:,iz,1] = interpolate_to_grid_vpa(wpa, temp[:,iz,1], vpa, vpa_spectral) + end + @test isapprox(expected.f_neutral[:, :, tind], newgrid_f_neutral[:,:,1], rtol=rtol) + end + end + + # Test initial values + test_values(1) + + # Test final values + test_values(2) + end +end + + +function runtests() + @testset "Krook collisions" verbose=use_verbose begin + println("Krook collisions tests") + + # Benchmark data is taken from this run (full-f with no splitting) + @testset "full-f" begin + run_test(test_input_full_f, 1.e-10, 3.e-16) + end + @testset "split 1" begin + run_test(test_input_split_1_moment, 1.e-3, 1.e-15) + end + @testset "split 2" begin + run_test(test_input_split_2_moments, 1.e-3, 1.e-15) + end + @testset "split 3" begin + run_test(test_input_split_3_moments, 1.e-3, 1.e-15) + end + end +end + +end # KrookCollisionsTests + + +using .KrookCollisionsTests + +KrookCollisionsTests.runtests() diff --git a/test/nonlinear_sound_wave_tests.jl b/test/nonlinear_sound_wave_tests.jl index 703c71b21..471e5c1b2 100644 --- a/test/nonlinear_sound_wave_tests.jl +++ b/test/nonlinear_sound_wave_tests.jl @@ -8,10 +8,11 @@ using TimerOutputs using moment_kinetics.coordinates: define_coordinate using moment_kinetics.input_structs: grid_input, advection_input -using moment_kinetics.load_data: open_readonly_output_file, load_coordinate_data, load_species_data, - load_fields_data, load_charged_particle_moments_data, load_pdf_data, - load_neutral_particle_moments_data, load_neutral_pdf_data, load_time_data, - load_species_data +using moment_kinetics.load_data: open_readonly_output_file, load_coordinate_data, + load_species_data, load_fields_data, + load_charged_particle_moments_data, load_pdf_data, + load_neutral_particle_moments_data, + load_neutral_pdf_data, load_time_data, load_species_data using moment_kinetics.interpolation: interpolate_to_grid_z, interpolate_to_grid_vpa using moment_kinetics.type_definitions: mk_float @@ -51,76 +52,77 @@ const expected = [z for z in range(-0.5 * z_L, 0.5 * z_L, length=11)], [vpa for vpa in range(-0.2 * vpa_L, 0.2 * vpa_L, length=3)], # Expected phi: - [-1.3862820803244256 -1.2383045504753758; -1.211510602668698 -1.1306858553168957; - -0.860938418534854 -0.8726873701297669; -0.5495322983936358 -0.5902920278548919; - -0.3534144494723056 -0.3756206847277757; -0.2876820724518619 -0.2919957813382994; - -0.35341444947230544 -0.37562068472777577; -0.5495322983936355 -0.5902920278548919; - -0.8609384185348539 -0.8726873701297669; -1.2115106026686981 -1.130685855316896; - -1.3862820803244256 -1.2383045504753758], + [-1.386282080324426 -1.2382641646968997; -1.2115129555832849 -1.130635565831393; + -0.8609860698164498 -0.872637046489647; -0.5494724768120176 -0.5903597644911374; + -0.35345976364887166 -0.37552658974835484; -0.28768207245186167 -0.2921445534164449; + -0.353459763648872 -0.3755265897483555; -0.5494724768120175 -0.5903597644911376; + -0.8609860698164502 -0.8726370464896476; -1.2115129555832849 -1.1306355658313922; + -1.3862820803244258 -1.2382641646968997], # Expected n_charged: - [0.2500030702177184 0.2898752704335188; 0.2977471383217195 0.3227988953227183; - 0.42274614626845974 0.417842206578383; 0.5772539714051019 0.5541536351162784; - 0.702254450621661 0.686868306132489; 0.7499999999999392 0.7467716863438243; - 0.702254450621661 0.6868683061324891; 0.577253971405102 0.5541536351162781; - 0.42274614626845963 0.41784220657838306; 0.29774713832171945 0.3227988953227184; - 0.25000307021771856 0.2898752704335188], + [0.2500030702177186 0.2898869775083742; 0.2977473631375158 0.3228278662412625; + 0.42274585818529853 0.417848119539277; 0.5772542465450629 0.5541281150892785; + 0.7022542481909738 0.6869277664245242; 0.7499999999999394 0.7466605958319346; + 0.7022542481909738 0.6869277664245237; 0.577254246545063 0.5541281150892783; + 0.42274585818529864 0.41784811953927686; 0.2977473631375159 0.32282786624126253; + 0.2500030702177185 0.2898869775083743], # Expected n_neutral: - [0.7499999999999392 0.7737996616909211; 0.702254450621661 0.7056147469533546; - 0.5772539714051019 0.5583199869826109; 0.42274614626845974 0.4096819689829928; - 0.29774713832171956 0.30537567265010457; 0.2500030702177185 0.26816544412496246; - 0.2977471383217197 0.30537567265010435; 0.4227461462684595 0.4096819689829924; - 0.5772539714051017 0.5583199869826102; 0.7022544506216611 0.7056147469533546; - 0.7499999999999394 0.7737996616909211], + [0.7499999999999382 0.7736769553678673; 0.7022542481909748 0.7056866352169496; + 0.5772542465450632 0.5582977481633454; 0.4227458581852985 0.40969188756651037; + 0.29774736313751604 0.30539644783353687; 0.2500030702177186 0.268198658560817; + 0.29774736313751604 0.305396447833537; 0.42274585818529836 0.4096918875665103; + 0.5772542465450631 0.5582977481633457; 0.7022542481909745 0.7056866352169494; + 0.7499999999999383 0.7736769553678673], # Expected upar_charged: - [1.1971912119126474e-17 -2.0968470015869656e-16; - -5.818706134342973e-17 -0.18232119131671534; - 9.895531571141618e-17 -0.1967239995126128; - -8.38774587436108e-18 -0.11125439467488389; - -1.7717259792356293e-17 -0.033747618153236424; - -3.143783114880992e-17 9.84455572616838e-17; - -2.499642278253839e-17 0.03374761815323648; - -2.9272523290371316e-17 0.1112543946748839; - 3.346728365577734e-17 0.19672399951261313; -8.193702949354942e-17 0.18232119131671523; - 1.1971912119126474e-17 -2.0187639060277426e-16], + [-2.7135787559953277e-17 -6.299214622140781e-17; + -9.321028970172899e-18 -0.1823721921091055; + -2.8374879811351724e-18 -0.19657035490893093; + 1.2124327390522635e-17 -0.11139486685283827; + 3.6525788403693063e-17 -0.033691837771623996; + -2.0930856430671915e-17 4.84147091991613e-17; + 8.753545920086251e-18 0.033691837771624024; + 1.1293771270243255e-17 0.11139486685283813; + 1.3739171132886587e-17 0.19657035490893102; + -6.840453743089351e-18 0.18237219210910513; + -2.7135787559953277e-17 -4.656897959900552e-17], # Expected upar_neutral: - [-3.143783114880993e-17 -3.003240017784847e-17; - -1.7717259792356296e-17 -0.03618184473095593; - -8.38774587436108e-18 -0.009226347310827927; - 9.895531571141618e-17 0.054572308562824384; -5.818706134342965e-17 0.0761077077764258; - 1.1971912119126477e-17 2.1033522146218786e-16; - 4.747047511913174e-17 -0.07610770777642595; - 6.558391323223502e-18 -0.05457230856282445; - 2.2213638810498713e-17 0.009226347310827925; - -2.7413075842225616e-17 0.036181844730955835; - -3.143783114880993e-17 -2.810175715538666e-17], + [6.5569385065066925e-18 7.469475342067322e-17; + 1.1054500872839027e-17 -0.036209130454625794; + -3.241833393685864e-17 -0.00915544640981337; + -3.617637280460899e-17 0.05452268209340691; 4.417578961284041e-17 0.07606644718003618; + 4.9354467746194965e-17 4.452343983947504e-17; + 6.573091229872379e-18 -0.07606644718003616; + 2.989662686945165e-17 -0.05452268209340687; + -3.1951996361666834e-17 0.009155446409813412; + -4.395464518158184e-18 0.03620913045462582; + 6.5569385065066925e-18 7.150569974151354e-17], # Expected ppar_charged: - [0.18749999999999997 0.23278940073547755; 0.20909100943488423 0.21912527958959363; - 0.24403280042122125 0.20817795270356831; 0.2440328004212212 0.21516422119834766; - 0.20909100943488412 0.2208129180125869; 0.18750000000000003 0.2213757117801786; - 0.2090910094348841 0.22081291801258685; 0.244032800421221 0.21516422119834752; - 0.24403280042122122 0.2081779527035683; 0.2090910094348842 0.21912527958959355; - 0.18749999999999992 0.23278940073547752], + [0.18749999999999992 0.2328164829490338; 0.20909325514551116 0.21912575009260987; + 0.24403180771238264 0.20822611102296495; 0.24403180771238278 0.21506741942934832; + 0.2090932551455113 0.22097085045011763; 0.1875 0.22119050467096843; + 0.20909325514551128 0.2209708504501176; 0.2440318077123828 0.2150674194293483; + 0.24403180771238256 0.20822611102296476; 0.20909325514551116 0.21912575009260982; + 0.18749999999999992 0.2328164829490338], # Expected ppar_neutral: - [0.18750000000000003 0.2482565420443467; 0.20909100943488412 0.24384477300624322; - 0.2440328004212212 0.228696297221881; 0.24403280042122125 0.20584407392704468; - 0.20909100943488423 0.19265693636741701; 0.18749999999999994 0.19084861718939317; - 0.2090910094348842 0.19265693636741688; 0.24403280042122116 0.20584407392704462; - 0.2440328004212211 0.228696297221881; 0.2090910094348841 0.2438447730062432; - 0.18750000000000003 0.24825654204434672], + [0.18750000000000003 0.2480244331470989; 0.20909325514551122 0.2440075646485762; + 0.24403180771238286 0.22861256884534023; 0.24403180771238278 0.20588932618946498; + 0.20909325514551144 0.19263633346696638; 0.18749999999999992 0.19091848744561835; + 0.20909325514551141 0.19263633346696654; 0.2440318077123828 0.20588932618946482; + 0.24403180771238286 0.22861256884534029; 0.20909325514551114 0.24400756464857642; + 0.18750000000000006 0.24802443314709893], # Expected f_charged: - [0.03704623609948259 0.04056128509273146 0.04289169811317835 0.030368915327672292 0.01235362235033934 0.0063385294703834204 0.012353622350339327 0.030368915327672247 0.04289169811317828 0.04056128509273145 0.0370462360994826; - 0.20411991941198782 0.251156132910555 0.3935556226209418 0.6276758497903185 0.9100827333021343 1.06066017177965 0.9100827333021342 0.6276758497903192 0.3935556226209421 0.25115613291055494 0.2041199194119877; - 0.03704623609948259 0.04056128509273146 0.04289169811317835 0.030368915327672292 0.01235362235033934 0.0063385294703834204 0.012353622350339327 0.030368915327672247 0.04289169811317828 0.04056128509273145 0.0370462360994826;;; - 0.05394101807537287 0.060554592436498814 0.036833910331906125 0.013633122209675089 0.010808508772375046 0.019345197472213343 0.027958746006592806 0.027628128813266543 0.026659296935378614 0.035613334811632306 0.05394101807537285; - 0.21177593449262566 0.24890460398430211 0.3731260689313831 0.5960741409510352 0.8872166610615642 1.0533874354116926 0.8872166610615648 0.5960741409510364 0.37312606893138234 0.24890460398430203 0.21177593449262566; - 0.053941018075372965 0.03561333481163235 0.026659296935378617 0.027628128813266564 0.02795874600659282 0.019345197472213388 0.010808508772375068 0.013633122209675051 0.03683391033190611 0.06055459243649881 0.05394101807537297], + [0.0370462360994826 0.04059927063892091 0.0428431419871786 0.030398267195914062 0.01236045902698859 0.006338529470383425 0.012360459026988587 0.030398267195914028 0.04284314198717859 0.0405992706389209 0.0370462360994826; + 0.20411991941198782 0.25123395823993105 0.3934413727192304 0.6277900619432855 0.9100364506661008 1.0606601717796504 0.910036450666101 0.6277900619432859 0.39344137271923046 0.25123395823993094 0.20411991941198776; + 0.0370462360994826 0.04059927063892091 0.0428431419871786 0.030398267195914062 0.01236045902698859 0.006338529470383425 0.012360459026988587 0.030398267195914028 0.04284314198717859 0.0405992706389209 0.0370462360994826;;; + 0.05392403019146985 0.06057819609646438 0.03676744157455075 0.013740507879552622 0.010777319583092297 0.019330359159894384 0.027982173790396116 0.027603104735767332 0.02667986700464528 0.035654512254837005 0.05392403019146984; + 0.21177720235387912 0.24902901234066305 0.3729377138332225 0.596281539172339 0.8870867512643452 1.0533860567375264 0.887086751264345 0.5962815391723388 0.3729377138332225 0.24902901234066285 0.21177720235387912; + 0.053924030191469796 0.035654512254837074 0.02667986700464531 0.02760310473576733 0.02798217379039615 0.019330359159894287 0.010777319583092311 0.013740507879552624 0.03676744157455069 0.060578196096464365 0.05392403019146979], # Expected f_neutral: - [0.006338529470383422 0.012353622350339336 0.030368915327672292 0.04289169811317835 0.04056128509273145 0.03704623609948259 0.04056128509273144 0.04289169811317833 0.03036891532767228 0.012353622350339327 0.00633852947038342; - 1.0606601717796502 0.9100827333021341 0.6276758497903185 0.3935556226209418 0.25115613291055494 0.2041199194119877 0.25115613291055505 0.3935556226209418 0.6276758497903185 0.9100827333021344 1.0606601717796504; - 0.006338529470383422 0.012353622350339336 0.030368915327672292 0.04289169811317835 0.04056128509273145 0.03704623609948259 0.04056128509273144 0.04289169811317833 0.03036891532767228 0.012353622350339327 0.00633852947038342;;; - 0.02430266037826662 0.040681127671090396 0.04194789083035397 0.036333231317680646 0.03689201427762576 0.04167458210401646 0.03666641377414817 0.019366777795459547 0.00833763923889437 0.009995913879120842 0.024302660378266627; - 1.0530041566990045 0.9037244245778953 0.6249909697155338 0.39563761233111516 0.2570339174716525 0.21139265577993924 0.25703391747165233 0.3956376123311148 0.6249909697155333 0.903724424577895 1.053004156699005; - 0.0243026603782666 0.009995913879120834 0.00833763923889439 0.01936677779545955 0.036666413774148206 0.041674582104016505 0.036892014277625805 0.036333231317680584 0.04194789083035393 0.04068112767109046 0.024302660378266613]) + [0.0063385294703834595 0.012360459026988546 0.030398267195914108 0.04284314198717859 0.040599270638920985 0.03704623609948259 0.040599270638920965 0.0428431419871786 0.030398267195914094 0.012360459026988546 0.006338529470383456; + 1.0606601717796493 0.9100364506661016 0.6277900619432857 0.3934413727192303 0.2512339582399308 0.20411991941198754 0.2512339582399307 0.3934413727192301 0.6277900619432853 0.9100364506661016 1.0606601717796487; + 0.0063385294703834595 0.012360459026988546 0.030398267195914108 0.04284314198717859 0.040599270638920985 0.03704623609948259 0.040599270638920965 0.0428431419871786 0.030398267195914094 0.012360459026988546 0.006338529470383456;;; + 0.024285034070612683 0.04071236753946936 0.04190483876050118 0.036374533667106086 0.0369234055803037 0.04165072188572372 0.03672486160719089 0.019283695804388743 0.008424202166370513 0.010011778724856858 0.02428503407061268; + 1.05300288883775 0.9036794996386066 0.6251037063201776 0.39552476644559265 0.25711045639921726 0.2113940344541052 0.25711045639921726 0.39552476644559253 0.6251037063201775 0.9036794996386066 1.0530028888377503; + 0.024285034070612672 0.01001177872485688 0.00842420216637052 0.019283695804388705 0.03672486160719087 0.04165072188572364 0.0369234055803037 0.036374533667106045 0.041904838760501203 0.040712367539469434 0.02428503407061266]) # default inputs for tests test_input_finite_difference = Dict("n_ion_species" => 1, @@ -192,13 +194,14 @@ test_input_finite_difference_split_2_moments = test_input_finite_difference_split_3_moments = merge(test_input_finite_difference_split_2_moments, Dict("run_name" => "finite_difference_split_3_moments", - "evolve_moments_parallel_pressure" => true)) + "evolve_moments_parallel_pressure" => true, + "vpa_L" => 12.0, "vz_L" => 12.0)) test_input_chebyshev = merge(test_input_finite_difference, Dict("run_name" => "chebyshev_pseudospectral", "z_discretization" => "chebyshev_pseudospectral", "z_ngrid" => 9, - "z_nelement" => 2, + "z_nelement" => 4, "vpa_discretization" => "chebyshev_pseudospectral", "vpa_ngrid" => 17, "vpa_nelement" => 8, @@ -209,8 +212,7 @@ test_input_chebyshev = merge(test_input_finite_difference, test_input_chebyshev_split_1_moment = merge(test_input_chebyshev, Dict("run_name" => "chebyshev_pseudospectral_split_1_moment", - "evolve_moments_density" => true, - "z_nelement" => 4)) + "evolve_moments_density" => true)) test_input_chebyshev_split_2_moments = merge(test_input_chebyshev_split_1_moment, @@ -220,7 +222,8 @@ test_input_chebyshev_split_2_moments = test_input_chebyshev_split_3_moments = merge(test_input_chebyshev_split_2_moments, Dict("run_name" => "chebyshev_pseudospectral_split_3_moments", - "evolve_moments_parallel_pressure" => true)) + "evolve_moments_parallel_pressure" => true, + "vpa_L" => 12.0, "vz_L" => 12.0)) # Not actually used in the tests, but needed for first argument of run_moment_kinetics @@ -259,20 +262,23 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) input["run_name"] = name # Suppress console output while running - phi = undef - n_charged = undef - upar_charged = undef - ppar_charged = undef - f_charged = undef - n_neutral = undef - upar_neutral = undef - ppar_neutral = undef - f_neutral = undef quietoutput() do # run simulation run_moment_kinetics(to, input) end + phi = nothing + n_charged = nothing + upar_charged = nothing + ppar_charged = nothing + f_charged = nothing + n_neutral = nothing + upar_neutral = nothing + ppar_neutral = nothing + f_neutral = nothing + z, z_spectral = nothing, nothing + vpa, vpa_spectral = nothing, nothing + if global_rank[] == 0 quietoutput() do @@ -305,6 +311,7 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) # load particle distribution function (pdf) data f_charged_vpavperpzrst = load_pdf_data(fid) f_neutral_vzvrvzetazrst = load_neutral_pdf_data(fid) + vpa, vpa_spectral = load_coordinate_data(fid, "vpa") close(fid) @@ -341,27 +348,6 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) end end - # Create coordinates - # - # create the 'input' struct containing input info needed to create a coordinate - # adv_input not actually used in this test so given values unimportant - adv_input = advection_input("default", 1.0, 0.0, 0.0) - nrank_per_block = 0 # dummy value - irank = 0 # dummy value - comm = MPI.COMM_NULL # dummy value - element_spacing_option = "uniform" - input = grid_input("coord", test_input["z_ngrid"], test_input["z_nelement"], - test_input["z_nelement"], nrank_per_block, irank, - z_L, test_input["z_discretization"], "", - "periodic", #test_input["z_bc"], - adv_input,comm, element_spacing_option) - z, z_spectral = define_coordinate(input) - input = grid_input("coord", test_input["vpa_ngrid"], test_input["vpa_nelement"], - test_input["vpa_nelement"], nrank_per_block, irank, - vpa_L, test_input["vpa_discretization"], "", - test_input["vpa_bc"], adv_input, comm, element_spacing_option) - vpa, vpa_spectral = define_coordinate(input) - # Test against values interpolated onto 'expected' grid which is fairly coarse no we # do not have to save too much data in this file @@ -437,8 +423,23 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) newgrid_ppar_charged = interpolate_to_grid_z(expected.z, ppar_charged[:, :, tind], z, z_spectral) @test isapprox(expected.ppar_charged[:, tind], newgrid_ppar_charged[:,1], rtol=rtol) + newgrid_vth_charged = @. sqrt(2.0*newgrid_ppar_charged/newgrid_n_charged) newgrid_f_charged = interpolate_to_grid_z(expected.z, f_charged[:, :, :, tind], z, z_spectral) - newgrid_f_charged = interpolate_to_grid_vpa(expected.vpa, newgrid_f_charged, vpa, vpa_spectral) + temp = newgrid_f_charged + newgrid_f_charged = fill(NaN, length(expected.vpa), + size(newgrid_f_charged, 2), + size(newgrid_f_charged, 3), + size(newgrid_f_charged, 4)) + for iz ∈ 1:length(expected.z) + wpa = copy(expected.vpa) + if input["evolve_moments_parallel_flow"] + wpa .-= newgrid_upar_charged[iz,1] + end + if input["evolve_moments_parallel_pressure"] + wpa ./= newgrid_vth_charged[iz,1] + end + newgrid_f_charged[:,iz,1] = interpolate_to_grid_vpa(wpa, temp[:,iz,1], vpa, vpa_spectral) + end @test isapprox(expected.f_charged[:, :, tind], newgrid_f_charged[:,:,1], rtol=rtol) # Check neutral particle moments and f @@ -453,8 +454,23 @@ function run_test(test_input, rtol, atol, upar_rtol=nothing; args...) newgrid_ppar_neutral = interpolate_to_grid_z(expected.z, ppar_neutral[:, :, tind], z, z_spectral) @test isapprox(expected.ppar_neutral[:, tind], newgrid_ppar_neutral[:,:,1], rtol=rtol) + newgrid_vth_neutral = @. sqrt(2.0*newgrid_ppar_neutral/newgrid_n_neutral) newgrid_f_neutral = interpolate_to_grid_z(expected.z, f_neutral[:, :, :, tind], z, z_spectral) - newgrid_f_neutral = interpolate_to_grid_vpa(expected.vpa, newgrid_f_neutral, vpa, vpa_spectral) + temp = newgrid_f_neutral + newgrid_f_neutral = fill(NaN, length(expected.vpa), + size(newgrid_f_neutral, 2), + size(newgrid_f_neutral, 3), + size(newgrid_f_neutral, 4)) + for iz ∈ 1:length(expected.z) + wpa = copy(expected.vpa) + if input["evolve_moments_parallel_flow"] + wpa .-= newgrid_upar_neutral[iz,1] + end + if input["evolve_moments_parallel_pressure"] + wpa ./= newgrid_vth_neutral[iz,1] + end + newgrid_f_neutral[:,iz,1] = interpolate_to_grid_vpa(wpa, temp[:,iz,1], vpa, vpa_spectral) + end @test isapprox(expected.f_neutral[:, :, tind], newgrid_f_neutral[:,:,1], rtol=rtol) end end @@ -479,10 +495,10 @@ function runtests() @testset "FD split 1" begin run_test(test_input_finite_difference_split_1_moment, 1.e-3, 1.e-11) end - @testset_skip "grids need shift/scale for collisions" "FD split 2" begin + @testset "FD split 2" begin run_test(test_input_finite_difference_split_2_moments, 1.e-3, 1.e-11) end - @testset_skip "grids need shift/scale for collisions" "FD split 3" begin + @testset "FD split 3" begin run_test(test_input_finite_difference_split_3_moments, 1.e-3, 1.e-11) end @@ -494,10 +510,10 @@ function runtests() @testset "Chebyshev split 1" begin run_test(test_input_chebyshev_split_1_moment, 1.e-3, 1.e-15) end - @testset_skip "grids need shift/scale for collisions" "Chebyshev split 2" begin + @testset "Chebyshev split 2" begin run_test(test_input_chebyshev_split_2_moments, 1.e-3, 1.e-15) end - @testset_skip "grids need shift/scale for collisions" "Chebyshev split 3" begin + @testset "Chebyshev split 3" begin run_test(test_input_chebyshev_split_3_moments, 1.e-3, 1.e-15) end end diff --git a/test/runtests.jl b/test/runtests.jl index c84944dd9..4bab7c641 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,8 +7,10 @@ function runtests() include(joinpath(@__DIR__, "calculus_tests.jl")) include(joinpath(@__DIR__, "interpolation_tests.jl")) include(joinpath(@__DIR__, "loop_setup_tests.jl")) + include(joinpath(@__DIR__, "velocity_integral_tests.jl")) include(joinpath(@__DIR__, "sound_wave_tests.jl")) include(joinpath(@__DIR__, "nonlinear_sound_wave_tests.jl")) + include(joinpath(@__DIR__, "Krook_collisions_tests.jl")) include(joinpath(@__DIR__, "harrisonthompson.jl")) include(joinpath(@__DIR__, "wall_bc_tests.jl")) end diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl new file mode 100644 index 000000000..07d0c5027 --- /dev/null +++ b/test/velocity_integral_tests.jl @@ -0,0 +1,144 @@ +module VelocityIntegralTests + +include("setup.jl") + +using moment_kinetics.input_structs: grid_input, advection_input +using moment_kinetics.coordinates: define_coordinate +using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure +using moment_kinetics.array_allocation: allocate_float + +using MPI + +function runtests() + @testset "velocity integral tests" verbose=use_verbose begin + println("velocity integral tests") + + # Tolerance for tests + atol = 1.0e-13 + + # define inputs needed for the test + ngrid = 17 #number of points per element + nelement_local = 20 # number of elements per rank + nelement_global = nelement_local # total number of elements + Lvpa = 18.0 #physical box size in reference units + Lvperp = 9.0 #physical box size in reference units + bc = "" #not required to take a particular value, not used + # fd_option and adv_input not actually used so given values unimportant + discretization = "chebyshev_pseudospectral" + fd_option = "fourth_order_centered" + adv_input = advection_input("default", 1.0, 0.0, 0.0) + nrank = 1 + irank = 0 + comm = MPI.COMM_NULL + # create the 'input' struct containing input info needed to create a + # coordinate + vr_input = grid_input("vperp", 1, 1, 1, nrank, irank, 1.0, discretization, + fd_option, bc, adv_input, comm, "uniform") + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, irank, + Lvpa, discretization, fd_option, bc, adv_input, comm, + "uniform") + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, + irank, Lvpa, discretization, fd_option, bc, adv_input, + comm, "uniform") + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, nrank, + irank, Lvperp, discretization, fd_option, bc, adv_input, + comm, "uniform") + # create the coordinate struct 'x' + vpa, vpa_spectral = define_coordinate(vpa_input) + vperp, vperp_spectral = define_coordinate(vperp_input) + vz, vz_spectral = define_coordinate(vz_input) + vr, vr_spectral = define_coordinate(vr_input) + + dfn = allocate_float(vpa.n,vperp.n) + dfn1D = allocate_float(vz.n, vr.n) + + function pressure(ppar,pperp) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres + end + + @testset "2D isotropic Maxwellian" begin + # assign a known isotropic Maxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + pperp = 2.0/3.0 + pres = get_pressure(ppar,pperp) + mass = 1.0 + vth = sqrt(2.0*pres/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/vth^3)*exp( - ((vpa_val-upar)^2 + vperp_val^2)/vth^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + pres_test = pressure(ppar_test,pperp_test) + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(pres_test, pres; atol=atol) + end + + @testset "1D Maxwellian" begin + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + mass = 1.0 + vth = sqrt(2.0*ppar/(dens*mass)) + for ivz in 1:vz.n + for ivr in 1:vr.n + vz_val = vz.grid[ivz] + dfn1D[ivz,ivr] = (dens/vth)*exp( - ((vz_val-upar)^2)/vth^2 ) + end + end + dens_test = get_density(dfn1D,vz,vr) + upar_test = get_upar(dfn1D,vz,vr,dens_test) + ppar_test = get_ppar(dfn1D,vz,vr,upar_test) + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + end + + @testset "biMaxwellian" begin + # assign a known biMaxwellian distribution in normalised units + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 4.0/5.0 + pperp = 1.0/4.0 + mass = 1.0 + vthpar = sqrt(2.0*ppar/(dens*mass)) + vthperp = sqrt(2.0*pperp/(dens*mass)) + for ivperp in 1:vperp.n + for ivpa in 1:vpa.n + vpa_val = vpa.grid[ivpa] + vperp_val = vperp.grid[ivperp] + dfn[ivpa,ivperp] = (dens/(vthpar*vthperp^2))*exp( - ((vpa_val-upar)^2)/vthpar^2 - (vperp_val^2)/vthperp^2 ) + end + end + + # now check that we can extract the correct moments from the distribution + + dens_test = get_density(dfn,vpa,vperp) + upar_test = get_upar(dfn,vpa,vperp,dens_test) + ppar_test = get_ppar(dfn,vpa,vperp,upar_test) + pperp_test = get_pperp(dfn,vpa,vperp) + + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + end + end +end + +end # VelocityIntegralTests + +using .VelocityIntegralTests + +VelocityIntegralTests.runtests() diff --git a/util/compare_collision_frequencies.jl b/util/compare_collision_frequencies.jl new file mode 100644 index 000000000..8396e5cb7 --- /dev/null +++ b/util/compare_collision_frequencies.jl @@ -0,0 +1,159 @@ +using moment_kinetics +using Plots +using Unitful + +function compare_collision_frequencies(input_file::String, + output_file::Union{String,Nothing}=nothing) + + input = moment_kinetics.moment_kinetics_input.read_input_file(input_file) + + io_input, evolve_moments, t_input, z_input, r_input, vpa_input, vperp_input, + gyrophase_input, vz_input, vr_input, vzeta_input, composition, species, collisions, + geometry, drive_input, num_diss_params, manufactured_solns_input = + moment_kinetics.moment_kinetics_input.mk_input(input) + + reference_params = moment_kinetics.reference_parameters.setup_reference_parameters(input) + + dimensional_parameters = moment_kinetics.utils.get_unnormalized_parameters( + input_file; Bnorm=reference_params.Bref, Lnorm=reference_params.Lref, + Nnorm=reference_params.Nref, Tnorm=reference_params.Tref) + + println("Omega_i0 ", dimensional_parameters["Omega_i0"]) + println("rho_i0 ", dimensional_parameters["rho_i0"]) + println("Omega_e0 ", dimensional_parameters["Omega_e0"]) + println("rho_e0 ", dimensional_parameters["rho_e0"]) + + # Effective collision frequency for dissipation? + # v_∥ dissipation term is D d^2f/dv_∥^2. Inserting factors of c_ref, this is a bit like + # pitch angle scattering D cref^2 d^2f/dv_∥^2 ~ D d^2f/dξ^2, so D is similar to a + # (normalised) collision frequency. + if num_diss_params.vpa_dissipation_coefficient < 0.0 + nu_vpa_diss = 0.0 + else + nu_vpa_diss = num_diss_params.vpa_dissipation_coefficient / + dimensional_parameters["timenorm"] + end + + println("ionization rate coefficient = ", + dimensional_parameters["ionization_rate_coefficient"]) + println("charge_exchange rate coefficient = ", + dimensional_parameters["CX_rate_coefficient"]) + + println("nu_ei0 ", dimensional_parameters["nu_ei0"]) + println("nu_ii0 ", dimensional_parameters["nu_ii0"]) + println("nu_ie0 ", dimensional_parameters["nu_ie0"]) + println("nu_vpa_diss ", nu_vpa_diss) + + # Neutral collison rates: + # The ionization term in the ion/neutral kinetic equations is ±R_ion*n_e*f_n. + # R_ion*n_e is an 'ionization rate' that just needs unnormalising - it gives the + # (inverse of the) characteristic time that it takes a neutral atom to be ionized. + nu_ionization0 = @. collisions.ionization * dimensional_parameters["Nnorm"] / + dimensional_parameters["timenorm"] + println("nu_ionization0 ", nu_ionization0) + # The charge-exchange term in the ion kinetic equation is -R_in*(n_n*f_i-n_i*f_n). + # So the rate at which ions experience CX reactions is R_in*n_n + nu_cx0 = @. collisions.charge_exchange * dimensional_parameters["Nnorm"] / + dimensional_parameters["timenorm"] + println("nu_cx0 ", nu_cx0) + + # Estimate classical particle and ion heat diffusion coefficients, for comparison to + # numerical dissipation + # Classical particle diffusivity estimate from Helander, D_⟂ on p.7. + # D_⟂ ∼ nu_ei * rho_e^2 / 2 + classical_particle_D0 = Unitful.upreferred(dimensional_parameters["nu_ei0"] * + dimensional_parameters["rho_e0"]^2 / 2.0) + println("classical_particle_D0 ", classical_particle_D0) + + # Classical thermal diffusivity estimate from Helander, eq. (1.8) + # chi_i = rho_i^2 / tau_ii / 2 = nu_ii * rho_i^2 / 2 + classical_heat_chi_i0 = Unitful.upreferred(dimensional_parameters["nu_ii0"] * + dimensional_parameters["rho_i0"]^2 / 2.0) + println("classical_heat_chi_i0 ", classical_heat_chi_i0) + + # rhostar is set as an input parameter. For the purposes of cross field transport, + # effective rho_i is rhostar*R rather than rho_i0. Might as well take R∼Lnorm for now, + # as difference will be O(1), if any. + rho_i_effective = Unitful.upreferred(geometry.rhostar * + dimensional_parameters["Lnorm"]) + rho_e_effective = + Unitful.upreferred(sqrt(dimensional_parameters["me"]/dimensional_parameters["mi"]) + * rho_i_effective) + effective_classical_particle_D0 = + Unitful.upreferred(dimensional_parameters["nu_ei0"] * rho_e_effective^2 / 2.0) + effective_classical_heat_chi_i0 = + Unitful.upreferred(dimensional_parameters["nu_ii0"] * rho_i_effective^2 / 2.0) + println("rho_i_effective ", rho_i_effective) + println("rho_e_effective ", rho_e_effective) + println("classical particle D0 with effective rho_e ", effective_classical_particle_D0) + println("classical heat chi_i0 with effective rho_i ", effective_classical_heat_chi_i0) + + # Get numerical diffusion parameters + if num_diss_params.r_dissipation_coefficient < 0.0 + D_r = 0.0 + else + D_r = Unitful.upreferred(num_diss_params.r_dissipation_coefficient * + dimensional_parameters["Lnorm"]^2 / + dimensional_parameters["timenorm"]) + end + if num_diss_params.z_dissipation_coefficient < 0.0 + D_z = 0.0 + else + D_z = Unitful.upreferred(num_diss_params.z_dissipation_coefficient * + dimensional_parameters["Lnorm"]^2 / + dimensional_parameters["timenorm"]) + end + println("numerical D_r ", D_r) + println("numerical D_z ", D_z) + + if output_file !== nothing + println("") + + temp, file_ext = splitext(output_file) + temp, ext = splitext(temp) + basename, iblock = splitext(temp) + iblock = parse(moment_kinetics.type_definitions.mk_int, iblock[2:end]) + fid = moment_kinetics.load_data.open_readonly_output_file(basename, ext[2:end]; + iblock=iblock) + + z = moment_kinetics.load_data.load_coordinate_data(fid, "z") + + density, parallel_flow, parallel_pressure, parallel_heat_flux, thermal_speed, + evolve_ppar = moment_kinetics.load_data.load_charged_particle_moments_data(fid) + + neutral_density, neutral_uz, neutral_pz, neutral_qz, neutral_thermal_speed = + moment_kinetics.load_data.load_neutral_particle_moments_data(fid) + + parallel_temperature = parallel_pressure ./ density + + # Ignoring variations in logLambda... + nu_ii = @. dimensional_parameters["nu_ii0"] * density / parallel_temperature^1.5 + println("nu_ii ", nu_ii[z.n_global÷2,1,1,end]) + + # Neutral collison rates: + # The ionization term in the ion/neutral kinetic equations is ±R_ion*n_e*f_n. + # R_ion*n_e is an 'ionization rate' that just needs unnormalising - it gives the + # (inverse of the) characteristic time that it takes a neutral atom to be ionized. + nu_ionization = @. collisions.ionization * density[:,:,1,:] / + dimensional_parameters["timenorm"] + println("nu_ionization ", nu_ionization[z.n_global÷2,1,end]) + # The charge-exchange term in the ion kinetic equation is -R_in*(n_n*f_i-n_i*f_n). + # So the rate at which ions experience CX reactions is R_in*n_n + nu_cx = @. collisions.charge_exchange * neutral_density[:,:,1,:] / + dimensional_parameters["timenorm"] + println("nu_cx ", nu_cx[z.n_global÷2,1,end]) + + # Make plot (using values from the final time point) + plot(legend=:outerright, xlabel="z", ylabel="frequency", ylims=(0.0, :auto)) + @views plot!(z.grid, nu_ii[:,1,1,end], label="nu_ii") + @views plot!(z.grid, nu_ionization[:,1,end], label="nu_ionization") + @views plot!(z.grid, nu_cx[:,1,end], label="nu_cx") + hline!([nu_vpa_diss], label="nu_vpa_diss") + ylabel!("frequency (s^-1)") + + savefig(joinpath(io_input.output_dir, + io_input.run_name * "_collision_frequencies.pdf")) + end + + return nothing +end