Add module with physical constants, mk_input return reference_parameters

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

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 = mk_input(input)
+        manufactured_solns_input, reference_parameters = 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

@@ -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")
@@ -352,7 +353,8 @@ 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 = input
+        drive_input, num_diss_params, manufactured_solns_input,
+        reference_parameters = input
 
     # Create loop range variables for shared-memory-parallel loops
     if debug_loop_type === nothing

src/moment_kinetics_input.jl

@@ -12,6 +12,7 @@ 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 ..finite_differences: fd_check_option
 using ..input_structs
 using ..numerical_dissipation: setup_numerical_dissipation
@@ -90,27 +91,26 @@ function mk_input(scan_input=Dict(); save_inputs_to_txt=false, ignore_MPI=true)
     composition.Er_constant = get(scan_input, "Er_constant", 0.0)
 
     # Get reference parameters for normalizations
-    reference_parameter_section = set_defaults_and_check_section!(
+    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)
 
-    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", reference_parameters.cref/reference_parameters.Lref/reference_parameters.Omegaref)
     #println("Info: Bzed is ",geometry.Bzed)
     #println("Info: Bmag is ",geometry.Bmag)
     #println("Info: rhostar is ",geometry.rhostar)
@@ -518,7 +518,8 @@ 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)
+                  num_diss_params, manufactured_solns_input,
+                  reference_parameters)
     println(io, "\nAll inputs returned from mk_input():")
     println(io, all_inputs)
     close(io)

src/plot_MMS_sequence.jl

@@ -66,7 +66,8 @@ 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 = mk_input(scan_input)
+        #    drive_input, num_diss_params, manufactured_solns_input,
+        #    reference_parameters = 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"