Krook 'collisions'

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"

docs/src/utils.md

@@ -0,0 +1,6 @@
+`utils`
+===============
+
+```@autodocs
+Modules = [moment_kinetics.utils]
+```

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/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
@@ -303,6 +306,8 @@ 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
 end
 
 """
@@ -502,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

src/krook_collisions.jl

@@ -0,0 +1,134 @@
+"""
+"""
+module 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_parameters)
+    Nref = reference_parameters.Nref
+    Tref = reference_parameters.Tref
+    mref = reference_parameters.mref
+    timeref = reference_parameters.timeref
+
+    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).
+    T0_Joules = Tref * proton_charge # Tref in Joules
+    mi_kg = mref # mi in kg
+    vth_i0 = sqrt(2.0 * T0_Joules / mi_kg) # vth_i at reference parameters in m.s^-1
+    nu_ii0_per_s = Nref * proton_charge^4 * logLambda_ii  /
+                   (4.0 * π * epsilon0^2 * mi_kg^2 * vth_i0^3) # s^-1
+    nu_ii0 = nu_ii0_per_s * timeref
+
+    return nu_ii0
+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
+        error("Krook collisions not implemented for 2V case 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]
+            T = (moments.charged.vth[iz,ir,is])^2
+            nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5)
+            @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]
+            T = vth^2
+            nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5)
+            @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]
+            T = vth^2
+            nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5)
+            @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]
+            T = vth^2
+            nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5)
+            @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]
+            T = vth^2
+            if vperp.n == 1
+                vth_prefactor = 1.0 / vth
+            else
+                vth_prefactor = 1.0 / vth^3
+            end
+            nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5)
+            @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

src/makie_post_processing.jl

@@ -18,11 +18,11 @@ 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 ..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,
@@ -778,10 +778,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, 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)