Move handling of reference parameters into separate module

Emphasises that the main code does not depend (directly) on the
reference parameters, rather they are needed in a few places during
setup to set defaults to the physically-consistent values of normalised
parameters (rhostar, collision frequency).

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
 

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)

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

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,27 +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 = 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)
     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", 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

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"

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

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

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