diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 1671326c5..6267cf209 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -10,11 +10,11 @@ Calculate normalized collision frequency at reference parameters for Coulomb col Currently valid only for hydrogenic ions (Z=1) """ -function setup_krook_collisions(reference_parameters) - Nref = reference_parameters.Nref - Tref = reference_parameters.Tref - mref = reference_parameters.mref - timeref = reference_parameters.timeref +function setup_krook_collisions(reference_params) + Nref = reference_params.Nref + Tref = reference_params.Tref + mref = reference_params.mref + timeref = reference_params.timeref Nref_per_cm3 = Nref * 1.0e-6 diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index c2a9e2704..c235af19c 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -780,7 +780,7 @@ function get_run_info(run_dir, restart_index=nothing; itime_min=1, itime_max=-1, # and check input to catch errors io_input, evolve_moments, t_input, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, _, composition, species, collisions, geometry, drive_input, num_diss_params, - manufactured_solns_input, reference_parameters = mk_input(input) + 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) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index 03c0a557f..0b475f183 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -27,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") @@ -355,8 +356,7 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, 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, - reference_parameters = input + drive_input, num_diss_params, manufactured_solns_input = input # Create loop range variables for shared-memory-parallel loops if debug_loop_type === nothing diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 752fb68c5..b194b8be8 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -12,11 +12,11 @@ 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 ..constants 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 @@ -91,19 +91,9 @@ 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 = copy(set_defaults_and_check_section!( - scan_input, "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_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) ## set geometry_input geometry.Bzed = get(scan_input, "Bzed", 1.0) @@ -111,7 +101,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) 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", reference_parameters.cref/reference_parameters.Lref/reference_parameters.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) @@ -180,7 +170,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) collisions.constant_ionization_rate = get(scan_input, "constant_ionization_rate", false) include_krook_collisions = get(scan_input, "krook_collisions", false) if include_krook_collisions - collisions.krook_collision_frequency_prefactor = setup_krook_collisions(reference_parameters) + collisions.krook_collision_frequency_prefactor = setup_krook_collisions(reference_params) else collisions.krook_collision_frequency_prefactor = -1.0 end @@ -525,8 +515,7 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true) 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, - reference_parameters) + num_diss_params, manufactured_solns_input) println(io, "\nAll inputs returned from mk_input():") println(io, all_inputs) close(io) @@ -1124,4 +1113,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/plot_MMS_sequence.jl b/src/plot_MMS_sequence.jl index 4b242b137..4d48e98a8 100644 --- a/src/plot_MMS_sequence.jl +++ b/src/plot_MMS_sequence.jl @@ -66,8 +66,7 @@ function get_MMS_error_data(path_list,scan_type,scan_name) #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, - # reference_parameters = mk_input(scan_input) + # 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/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/utils.jl b/src/utils.jl index 5f9c50636..ad8e5f19e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -5,13 +5,15 @@ 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" 1.602176634e-19*Unitful.J true +Unitful.@unit eV "eV" "electron volt" proton_charge*Unitful.J true function __init__() Unitful.register(utils) @@ -29,18 +31,20 @@ 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, - reference_parameters = mk_input(input) + drive_input, external_source_settings, num_diss_params, manufactured_solns_input = + mk_input(input) - Nnorm = reference_parameters.Nref * Unitful.m^(-3) - Tnorm = reference_parameters.Tref * eV - Lnorm = reference_parameters.Lref * Unitful.m - Bnorm = reference_parameters.Bref * Unitful.T - cnorm = reference_parameters.cref * Unitful.m / Unitful.s - timenorm = reference_parameters.timeref * Unitful.s + 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_parameters.mnorm * Unitful.kg + mi = reference_params.mref * Unitful.kg parameters = OrderedDict{String,Any}() parameters["run_name"] = run_name diff --git a/util/compare_collision_frequencies.jl b/util/compare_collision_frequencies.jl index 806aac54b..8396e5cb7 100644 --- a/util/compare_collision_frequencies.jl +++ b/util/compare_collision_frequencies.jl @@ -9,12 +9,14 @@ function compare_collision_frequencies(input_file::String, 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, - reference_parameters = moment_kinetics.moment_kinetics_input.mk_input(input) + 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_parameters.Bref, Lnorm=reference_parameters.Lref, - Nnorm=reference_parameters.Nref, Tnorm=reference_parameters.Tref) + 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"]) @@ -149,7 +151,7 @@ function compare_collision_frequencies(input_file::String, hline!([nu_vpa_diss], label="nu_vpa_diss") ylabel!("frequency (s^-1)") - savefig(joinpath("runs", io_input.run_name, + savefig(joinpath(io_input.output_dir, io_input.run_name * "_collision_frequencies.pdf")) end