From bffaea3f83d56afddbdeb24636f23fd9c6994d89 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 26 Aug 2023 12:53:18 +0100 Subject: [PATCH 01/49] Move generic settings definition functions to input_structs.jl This allows them to be called in other modules in a setup function, which can then be called in moment_kinetics_input.jl, avoiding an unresolvable cyclic import. --- src/input_structs.jl | 110 +++++++++++++++++++++++++++++++++++ src/makie_post_processing.jl | 6 +- src/moment_kinetics_input.jl | 108 ---------------------------------- 3 files changed, 113 insertions(+), 111 deletions(-) diff --git a/src/input_structs.jl b/src/input_structs.jl index b013f703b..b3cb44238 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 @@ -502,4 +504,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/makie_post_processing.jl b/src/makie_post_processing.jl index 1c4fd8600..3ce545ef2 100644 --- a/src/makie_post_processing.jl +++ b/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, diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 199941d6b..cce39cfb2 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -35,114 +35,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 From 6293496dc50fd38df4230f1b6248a4ea2557049b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 26 Aug 2023 17:34:03 +0100 Subject: [PATCH 02/49] Move define_coordinate() calls into mk_input() This is a litle bit cleaner, and allows setup of settings from input to use the coordinate objects. --- src/makie_post_processing.jl | 7 +++---- src/moment_kinetics.jl | 28 +++++----------------------- src/moment_kinetics_input.jl | 25 ++++++++++++++++++++++--- src/plot_MMS_sequence.jl | 9 ++++----- src/post_processing.jl | 9 +++------ 5 files changed, 37 insertions(+), 41 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 3ce545ef2..c235af19c 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -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 = 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 c6eb39282..d0100d16f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -71,7 +71,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!, @@ -350,28 +349,11 @@ function setup_moment_kinetics(input_dict::Dict; restart_prefix_iblock=nothing, 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 cce39cfb2..759841276 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -10,6 +10,7 @@ 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 ..finite_differences: fd_check_option using ..input_structs @@ -480,6 +481,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 @@ -496,9 +514,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) 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 9cee0bf5a..3e1a7f1c8 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1047,12 +1047,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 From fdac38e055a150766485e02b8acaaceaa5c326c0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 27 Aug 2023 15:47:29 +0100 Subject: [PATCH 03/49] Add module with physical constants, mk_input return reference_parameters --- src/constants.jl | 33 +++++++++++++++++++++++++++++++++ src/makie_post_processing.jl | 2 +- src/moment_kinetics.jl | 4 +++- src/moment_kinetics_input.jl | 19 ++++++++++--------- src/plot_MMS_sequence.jl | 3 ++- 5 files changed, 49 insertions(+), 12 deletions(-) create mode 100644 src/constants.jl 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/makie_post_processing.jl b/src/makie_post_processing.jl index c235af19c..c2a9e2704 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 = 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) diff --git a/src/moment_kinetics.jl b/src/moment_kinetics.jl index d0100d16f..bb0536db2 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") @@ -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 diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 759841276..31f0d61e4 100644 --- a/src/moment_kinetics_input.jl +++ b/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) diff --git a/src/plot_MMS_sequence.jl b/src/plot_MMS_sequence.jl index 4d48e98a8..4b242b137 100644 --- a/src/plot_MMS_sequence.jl +++ b/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" From df1bd5b2008ae133af6ebcfbf416eccccfde22d0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 30 Sep 2022 14:40:43 +0100 Subject: [PATCH 04/49] Function to print some dimensional parameters of a simulation --- Project.toml | 1 + docs/src/utils.md | 6 +++ src/moment_kinetics.jl | 1 + src/utils.jl | 93 ++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 101 insertions(+) create mode 100644 docs/src/utils.md create mode 100644 src/utils.jl 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/src/moment_kinetics.jl b/src/moment_kinetics.jl index bb0536db2..73bfd15b8 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -62,6 +62,7 @@ include("plot_sequence.jl") include("makie_post_processing.jl") include("time_advance.jl") +include("utils.jl") using TimerOutputs using Dates using Glob diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 000000000..fcf549b69 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,93 @@ +""" +Utility functions +""" +module utils + +export get_unnormalized_parameters, print_unnormalized_parameters + +using ..moment_kinetics_input: mk_input + +using OrderedCollections +using TOML +using Unitful + +Unitful.@unit eV "eV" "electron volt" 1.602176634e-19*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, + reference_parameters = 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 + + # Assume single ion species so normalised ion mass is always 1 + mi = reference_parameters.mnorm * 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 + + 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 From 33584ea23e181f9a6a0c4fd43b9ef8e467945dec Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 20 Mar 2023 20:39:59 +0000 Subject: [PATCH 05/49] Script to compare collision frequencies to numerical dissipation --- util/compare_collision_frequencies.jl | 157 ++++++++++++++++++++++++++ 1 file changed, 157 insertions(+) create mode 100644 util/compare_collision_frequencies.jl diff --git a/util/compare_collision_frequencies.jl b/util/compare_collision_frequencies.jl new file mode 100644 index 000000000..806aac54b --- /dev/null +++ b/util/compare_collision_frequencies.jl @@ -0,0 +1,157 @@ +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, + reference_parameters = moment_kinetics.moment_kinetics_input.mk_input(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) + + 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("runs", io_input.run_name, + io_input.run_name * "_collision_frequencies.pdf")) + end + + return nothing +end From 8d8bcd4995b1fd739d8b44f234e8a4e18587590c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 4 Oct 2022 15:27:48 +0100 Subject: [PATCH 06/49] Krook collision operator for ion-ion Coulomb collisions --- src/input_structs.jl | 3 + src/krook_collisions.jl | 129 +++++++++++++++++++++++++++++++++++ src/moment_kinetics.jl | 1 + src/moment_kinetics_input.jl | 11 ++- src/time_advance.jl | 31 +++++++-- src/utils.jl | 2 + 6 files changed, 170 insertions(+), 7 deletions(-) create mode 100644 src/krook_collisions.jl diff --git a/src/input_structs.jl b/src/input_structs.jl index b3cb44238..ed3d20914 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -61,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 @@ -305,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 """ diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl new file mode 100644 index 000000000..a4c40b5ee --- /dev/null +++ b/src/krook_collisions.jl @@ -0,0 +1,129 @@ +""" +""" +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 + 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 + * 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/moment_kinetics.jl b/src/moment_kinetics.jl index 73bfd15b8..03c0a557f 100644 --- a/src/moment_kinetics.jl +++ b/src/moment_kinetics.jl @@ -46,6 +46,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") diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 31f0d61e4..752fb68c5 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -13,6 +13,7 @@ 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 @@ -177,6 +178,12 @@ 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) + include_krook_collisions = get(scan_input, "krook_collisions", false) + if include_krook_collisions + collisions.krook_collision_frequency_prefactor = setup_krook_collisions(reference_parameters) + else + collisions.krook_collision_frequency_prefactor = -1.0 + end # parameters related to the time stepping nstep = get(scan_input, "nstep", 5) @@ -948,7 +955,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) Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength diff --git a/src/time_advance.jl b/src/time_advance.jl index 9d2a12d99..4892dd87e 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -40,6 +40,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 +405,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 +475,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 +526,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) @@ -1039,6 +1044,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 @@ -1652,6 +1665,12 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, collisions, dt) end + # Add a for Krook collision operoator 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 index fcf549b69..5f9c50636 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -63,6 +63,8 @@ function get_unnormalized_parameters(input::Dict) 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 From 94e14083b42ac79283218d94d35e8094aa6e69c1 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 21 Sep 2023 18:29:29 +0100 Subject: [PATCH 07/49] Implement "split 2" and "split 3" test for NonlinearSoundWaveTests These require interpolation in the velocity grids to compare to the expected distribution function values. --- test/nonlinear_sound_wave_tests.jl | 229 ++++++++++++++++------------- 1 file changed, 123 insertions(+), 106 deletions(-) diff --git a/test/nonlinear_sound_wave_tests.jl b/test/nonlinear_sound_wave_tests.jl index 2742760cd..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,26 +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 - 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) - 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) - 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 @@ -436,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 @@ -452,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 @@ -478,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 @@ -493,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 From fcdd763e57d72e51794225a1ff3779077275e57e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 22 Sep 2023 13:28:42 +0100 Subject: [PATCH 08/49] Add regression test for Krook collision operator Test is based on the NonlinearSoundWave test, but with Krook collisions added. Only test Chebyshev derivatives (not finite difference) to save time, as derivatives are not used in the Krook collision operator. --- test/Krook_collisions_tests.jl | 443 +++++++++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 444 insertions(+) create mode 100644 test/Krook_collisions_tests.jl diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl new file mode 100644 index 000000000..1ebae0805 --- /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" => true, + "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/runtests.jl b/test/runtests.jl index c84944dd9..17ba5030d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ function runtests() include(joinpath(@__DIR__, "loop_setup_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 From 9b21e092998320607989bbe228e90bb13df70512 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 22 Sep 2023 13:42:06 +0100 Subject: [PATCH 09/49] Handle 2V case in Krook collision operator If there is a vperp dimension (i.e. ion distribution function is not marginalised over vperp) then need a prefactor of 1/vth^3 rather than 1/vth (as for the 1V case) in front of the Maxwellian. --- src/krook_collisions.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index a4c40b5ee..751d6caa1 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -112,11 +112,16 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v 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 + - n * vth_prefactor * exp(-((vpa.grid[ivpa] - fvec_in.upar[iz,ir,is])/vth)^2 - (vperp.grid[ivperp]/vth)^2)) end From cf8a956293a23708d00e4589aafde9a55dabf844 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 25 Sep 2023 12:51:31 +0100 Subject: [PATCH 10/49] Restrict check for 1V in Krook collisions to moment-kinetic cases 2V Krook collisions are now implemented for full-f simulations. --- src/krook_collisions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 751d6caa1..1671326c5 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -42,8 +42,8 @@ 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") + 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 From 49f41f8536ce552fc05b8ab2d3a3fa6a031c215b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 25 Sep 2023 14:57:15 +0100 Subject: [PATCH 11/49] 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 ++++---- src/makie_post_processing.jl | 2 +- src/moment_kinetics.jl | 4 +-- src/moment_kinetics_input.jl | 34 ++++++++++++------------- src/plot_MMS_sequence.jl | 3 +-- src/reference_parameters.jl | 36 +++++++++++++++++++++++++++ src/utils.jl | 24 ++++++++++-------- util/compare_collision_frequencies.jl | 12 +++++---- 8 files changed, 82 insertions(+), 43 deletions(-) create mode 100644 src/reference_parameters.jl 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 From 030fbe780dd93d3d0893b2ed8dae54314297440a Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 17 May 2023 11:42:45 +0100 Subject: [PATCH 12/49] Ported capability to compare upar and ppar to expectations from MMS test -- currently the tests do not pass. --- src/manufactured_solns.jl | 61 +++++++++++++++++++++++++++++++++++++-- src/post_processing.jl | 13 +++++++++ 2 files changed, 71 insertions(+), 3 deletions(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 8f93faf2b..c432dd8c9 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,7 +201,48 @@ using IfElse * "manufactured_solns:type=$(manufactured_solns_input.type)") end end - + + # 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 = composition.epsilon_offset + alpha = composition.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) + 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 = composition.epsilon_offset + alpha = composition.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 + end + + # ion perpendicular pressure symbolic function + function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + pperpi = densi # simple vperp^2 dependence of dfni + return pperpi + end function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition, manufactured_solns_input) if z_bc == "periodic" @@ -318,6 +361,12 @@ 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) dfni = dfni_sym(Lr, Lz, r_bc, z_bc, composition, geometry, nr, manufactured_solns_input, charged_species) @@ -329,6 +378,9 @@ 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}) 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 +391,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) return manufactured_solns_list end diff --git a/src/post_processing.jl b/src/post_processing.jl index 3e1a7f1c8..2837bad8b 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1099,6 +1099,9 @@ function analyze_and_plot_data(prefix...; run_index=nothing) composition, species, r_global.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 dfnn_func = manufactured_solns_list.dfnn_func densn_func = manufactured_solns_list.densn_func manufactured_E_fields = @@ -1135,16 +1138,26 @@ 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[:,:,:,:]) 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]) 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_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") From 6a9664d2822fe9d847ea8d1cacdc4a70448a3f4d Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 16:17:00 +0100 Subject: [PATCH 13/49] Clean up manufactured solutions input Options for manufactured solutions were moved from 'composition' to a dedicated NamedTuple. However, the fields in the 'composition' struct were not removed - remove them, and fix bugs caused by those fields still having been used, although they were not being set correctly from the input options. --- src/input_structs.jl | 8 +------- src/manufactured_solns.jl | 13 ++++++------- src/moment_kinetics_input.jl | 9 +-------- 3 files changed, 8 insertions(+), 22 deletions(-) diff --git a/src/input_structs.jl b/src/input_structs.jl index ed3d20914..285e42bcd 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -259,13 +259,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 diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index c432dd8c9..000c342d7 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -211,8 +211,8 @@ using IfElse Er, Ez, phi = electric_fields(Lr,Lz,r_bc,z_bc,composition,nr,manufactured_solns_input,species) rhostar = geometry.rhostar bzed = geometry.bzed - epsilon = composition.epsilon_offset - alpha = composition.alpha_switch + 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 ) @@ -226,8 +226,8 @@ using IfElse ppari = 0.0 # not supported elseif z_bc == "wall" densi = densi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) - epsilon = composition.epsilon_offset - alpha = composition.alpha_switch + 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) @@ -243,8 +243,7 @@ using IfElse pperpi = densi # simple vperp^2 dependence of dfni return pperpi end - function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,composition, - manufactured_solns_input) + function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,manufactured_solns_input) if z_bc == "periodic" jpari_into_LHS_wall_sym = 0.0 elseif z_bc == "wall" @@ -316,7 +315,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 diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 752fb68c5..a15e1c733 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -844,20 +844,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) From 02ded9ac90a3b42f445d6bc863c86198da978224 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:07:56 +0100 Subject: [PATCH 14/49] Put 'show_element_boundaries' in postproc input for manufactured_solns --- src/makie_post_processing.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index c2a9e2704..1b0cebefe 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -673,6 +673,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"]) From bf27259ded892efa81c1ad29325ec9743bd7b01e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:08:50 +0100 Subject: [PATCH 15/49] Make MMS plots for ion parallel flow and parallel pressure --- src/makie_post_processing.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 1b0cebefe..3926cbec4 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -5438,6 +5438,7 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; 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,7 +5468,8 @@ 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, @@ -6353,6 +6355,8 @@ 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 From a1a17ec80fb942f2b9d7f6d831cc63e8bcfdcd7d Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:09:59 +0100 Subject: [PATCH 16/49] Skip MMS plots of Er for 1D cases These plots are uninteresting, and cause errors because the y-axis limits are both 0. --- src/makie_post_processing.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 3926cbec4..7b7dd62d2 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -6362,6 +6362,9 @@ function manufactured_solutions_analysis(run_info; plot_prefix) 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) From b539670d1c750253cd2885357f85568f04a92fd0 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Wed, 27 Sep 2023 17:11:03 +0100 Subject: [PATCH 17/49] Fix slicing of distribution functions in MMS plots --- src/makie_post_processing.jl | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index 7b7dd62d2..a98855b2a 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -5782,13 +5782,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 @@ -6079,7 +6085,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 @@ -6286,7 +6292,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 From f181300f27a94116cd8e94a014bf4ff5870e2513 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 28 Sep 2023 14:58:08 +0100 Subject: [PATCH 18/49] Fix post processing script consistent with updated composition struct. --- src/post_processing.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 2837bad8b..160a2a84d 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -311,8 +311,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 From d91c61866ca225543c009eb17af9d15e0a4a42f6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 29 Sep 2023 14:09:33 +0100 Subject: [PATCH 19/49] Addition of factor to compensate for the pressure normalisation convention in master, which seems to normalise pressure to 2 nref Tref rather than nref Tref = (1/2)mref nref cref^2. --- src/manufactured_solns.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 000c342d7..4cf0adcea 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -222,6 +222,8 @@ using IfElse # 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" @@ -234,7 +236,7 @@ using IfElse - (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 + return ppari*norm_fac end # ion perpendicular pressure symbolic function From 1932bf5dfe1c26bee57f5bb2fde67f76ce55adae Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Fri, 12 May 2023 11:24:07 +0100 Subject: [PATCH 20/49] Porting of commit cd07a26. Refactoring of the update_derived_moments function to use get_density, get_upar and get_ppar functions. In the split moment operation, the moments passed in to these functions along with the distribution are set to the enforced values of (density,upar) = (1.0,0.0) appropriately. This branch now supports the test_velocity_integrals.jl script. --- src/velocity_moments.jl | 76 ++++++++++++++++--------- test_velocity_integrals.jl | 112 +++++++++++++++++++++++++++++++++++++ 2 files changed, 161 insertions(+), 27 deletions(-) create mode 100644 test_velocity_integrals.jl diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 8db08f9c3..45e5aa9ae 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -22,6 +22,12 @@ 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 + using ..type_definitions: mk_float using ..array_allocation: allocate_shared_float, allocate_bool, allocate_float using ..calculus: integral @@ -431,13 +437,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 +484,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 +495,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 +505,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 +556,46 @@ 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 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 + """ NB: the incoming pdf is the normalized pdf """ diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl new file mode 100644 index 000000000..bb39123b1 --- /dev/null +++ b/test_velocity_integrals.jl @@ -0,0 +1,112 @@ +using Printf +using Plots +using LaTeXStrings +using MPI +using Measures + +if abspath(PROGRAM_FILE) == @__FILE__ + using Pkg + Pkg.activate(".") + + import moment_kinetics + 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 + using moment_kinetics.array_allocation: allocate_float + + # 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 = 12.0 #physical box size in reference units + Lvperp = 6.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 + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, + nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + # create the coordinate struct 'x' + println("made inputs") + println("vpa: ngrid: ",ngrid," nelement: ",nelement_local) + println("vperp: ngrid: ",ngrid," nelement: ",nelement_local) + vpa, vpa_spectral = define_coordinate(vpa_input) + vperp, vperp_spectral = define_coordinate(vperp_input) + + dfn = allocate_float(vpa.n,vperp.n) + + function pressure(ppar,pperp) + pres = (1.0/3.0)*(ppar + 2.0*pperp) + return pres + end + # 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 = 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) + # output test results + println("") + println("Isotropic Maxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) + println("") + + # 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) + # output test results + + println("") + println("biMaxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) + println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) + println("") +end From fa279702c201addd7840bea1218f6e0d04145497 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 17 May 2023 14:27:26 +0100 Subject: [PATCH 21/49] Cherry-pick of commit 8632cf7 to add perpendicular pressure to the body of the code --- src/file_io.jl | 28 ++++++++++++++++++++++++--- src/initial_conditions.jl | 6 ++++-- src/moment_kinetics_structs.jl | 1 + src/post_processing.jl | 11 +++++++++-- src/time_advance.jl | 12 +++++++----- src/velocity_moments.jl | 35 +++++++++++++++++++++++++++++++++- 6 files changed, 80 insertions(+), 13 deletions(-) diff --git a/src/file_io.jl b/src/file_io.jl index bb53a90f2..94f3e4123 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 @@ -587,6 +589,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; @@ -638,7 +647,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 @@ -777,6 +786,7 @@ function reopen_moments_io(file_info) 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("perpendicular_pressure"), getvar("thermal_speed"), getvar("density_neutral"), getvar("uz_neutral"), getvar("pz_neutral"), getvar("qz_neutral"), getvar("thermal_speed_neutral"), @@ -861,6 +871,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"), @@ -916,6 +927,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, @@ -1004,6 +1017,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, @@ -1267,7 +1282,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, @@ -1366,6 +1381,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 @@ -1442,6 +1462,7 @@ function debug_dump(fvec::Union{scratch_pdf,Nothing}, density = nothing upar = nothing ppar = nothing + pperp = nothing pdf_neutral = nothing density_neutral = nothing else @@ -1449,6 +1470,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 @@ -1462,7 +1484,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..92cb07ed8 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!, 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) @@ -924,6 +925,7 @@ 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, 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/post_processing.jl b/src/post_processing.jl index 160a2a84d..69287f396 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, @@ -488,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), @@ -531,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), @@ -1157,6 +1162,8 @@ function analyze_and_plot_data(prefix...; run_index=nothing) 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,perpendicular_pressure,pperp_sym,"ion",z,r,time,nz_global,nr_global,ntime, + L"\widetilde{p}_{\perp i}",L"\widetilde{p}_{\perp i}^{sym}",L"\varepsilon(\widetilde{p}_{\perp i})","pperp") 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/time_advance.jl b/src/time_advance.jl index 4892dd87e..46687ca34 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! 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! @@ -656,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...) @@ -665,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 @@ -673,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 @@ -1234,7 +1235,6 @@ 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 @@ -1406,7 +1406,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 """ @@ -1501,6 +1502,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 diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 45e5aa9ae..45e71cf72 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -10,6 +10,7 @@ export update_moments! export update_density! export update_upar! export update_ppar! +export update_pperp! export update_qpar! export reset_moments_status! export moments_chrg_substruct, moments_ntrl_substruct @@ -64,6 +65,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 @@ -188,6 +191,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 @@ -253,7 +258,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) @@ -590,6 +595,34 @@ function get_ppar(ff, vpa, vperp, upar) 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 From 82de4c4f911378a05d3f7807b0e2c1d47e1eccd0 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Thu, 18 May 2023 10:27:57 +0100 Subject: [PATCH 22/49] Attempt to cherry-pick the features from 6bd406e. Some bugs remain where pperp (and thus vthi) are computed incorrectly according to the manufactured solutions test in 1V and 2V. --- src/initial_conditions.jl | 11 +++------ src/manufactured_solns.jl | 40 +++++++++++++++++++++++++------ src/post_processing.jl | 10 ++++++-- src/time_advance.jl | 11 +++++---- src/velocity_moments.jl | 30 ++++++++++++++++++++++++ test_velocity_integrals.jl | 48 +++++++++++++++++++++++++++++++++----- 6 files changed, 122 insertions(+), 28 deletions(-) diff --git a/src/initial_conditions.jl b/src/initial_conditions.jl index 92cb07ed8..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!, update_pperp!, reset_moments_status! +using ..velocity_moments: update_ppar!, update_upar!, update_density!, update_pperp!, update_vth!, reset_moments_status! using ..manufactured_solns: manufactured_solutions @@ -899,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 @@ -931,12 +931,7 @@ function init_pdf_moments_manufactured_solns!(pdf, moments, vz, vr, vzeta, vpa, 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/manufactured_solns.jl b/src/manufactured_solns.jl index 4cf0adcea..ad00ec145 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -240,12 +240,35 @@ using IfElse end # ion perpendicular pressure symbolic function - function pperpi_sym(Lr,Lz,r_bc,z_bc,composition,manufactured_solns_input,species) + 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) - pperpi = densi # simple vperp^2 dependence of dfni - return pperpi + 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 - function jpari_into_LHS_wall_sym(Lr,Lz,r_bc,z_bc,manufactured_solns_input) + + # 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" @@ -352,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] @@ -367,7 +390,9 @@ using IfElse 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) + 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) @@ -382,6 +407,7 @@ using IfElse 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}) @@ -395,7 +421,7 @@ using IfElse 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) + pperpi_func = pperpi_func, vthi_func = vthi_func) return manufactured_solns_list end diff --git a/src/post_processing.jl b/src/post_processing.jl index 69287f396..1b701486d 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1017,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] @@ -1100,12 +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 = @@ -1145,6 +1147,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) 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 @@ -1153,6 +1156,7 @@ function analyze_and_plot_data(prefix...; run_index=nothing) 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 @@ -1162,8 +1166,10 @@ function analyze_and_plot_data(prefix...; run_index=nothing) 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,perpendicular_pressure,pperp_sym,"ion",z,r,time,nz_global,nr_global,ntime, + compare_moments_symbolic_test(run_name,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,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/time_advance.jl b/src/time_advance.jl index 46687ca34..e357b4025 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -15,7 +15,7 @@ using ..file_io: write_data_to_ascii, write_moments_data_to_binary, write_dfns_d using ..looping 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_pperp!, 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! @@ -1204,7 +1204,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 @@ -1240,9 +1242,7 @@ function rk_update!(scratch, pdf, moments, fields, boundary_distributions, vz, v # 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") @@ -1452,6 +1452,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 diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index 45e71cf72..da41a1fea 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -12,6 +12,7 @@ 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! @@ -28,6 +29,7 @@ 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 @@ -629,6 +631,34 @@ function get_pperp(ff, vpa, vperp) 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_velocity_integrals.jl b/test_velocity_integrals.jl index bb39123b1..24f674242 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -11,7 +11,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ import moment_kinetics 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 + using moment_kinetics.velocity_moments: get_density, get_upar, get_ppar, get_pperp, get_pressure using moment_kinetics.array_allocation: allocate_float # define inputs needed for the test @@ -30,29 +30,37 @@ if abspath(PROGRAM_FILE) == @__FILE__ comm = MPI.COMM_NULL # create the 'input' struct containing input info needed to create a # coordinate - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + vr_input = grid_input("vperp", 1, 1, 1, + nrank, irank, 1.0, discretization, fd_option, bc, adv_input,comm) + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) # create the coordinate struct 'x' println("made inputs") - println("vpa: ngrid: ",ngrid," nelement: ",nelement_local) - println("vperp: ngrid: ",ngrid," nelement: ",nelement_local) + println("vpa: ngrid: ",ngrid," nelement: ",nelement_local, " Lvpa: ",Lvpa) + println("vperp: ngrid: ",ngrid," nelement: ",nelement_local, " Lvperp: ",Lvperp) 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 + # 2D isotropic Maxwellian test # 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 = pressure(ppar,pperp) + pres = get_pressure(ppar,pperp) mass = 1.0 vth = sqrt(2.0*pres/(dens*mass)) for ivperp in 1:vperp.n @@ -72,12 +80,40 @@ if abspath(PROGRAM_FILE) == @__FILE__ pres_test = pressure(ppar_test,pperp_test) # output test results println("") - println("Isotropic Maxwellian") + println("Isotropic 2D Maxwellian") println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) println("") + ################### + # 1D Maxwellian test + + dens = 3.0/4.0 + upar = 2.0/3.0 + ppar = 2.0/3.0 + mass = 1.0 + vth = sqrt(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) + # output test results + println("") + println("1D Maxwellian") + println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) + println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) + println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) + println("") + + ################### + # biMaxwellian test + # assign a known biMaxwellian distribution in normalised units dens = 3.0/4.0 upar = 2.0/3.0 From ed979adc05a343dc134bd8ae92f79ebd63070d11 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Oct 2023 10:57:40 +0100 Subject: [PATCH 23/49] Include factor of 2 for expected vth in test_velocity_integrals.jl Needs the factor of 2 to be consistent with normalisations currently used. --- test_velocity_integrals.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl index 24f674242..9bcee9325 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -93,7 +93,7 @@ if abspath(PROGRAM_FILE) == @__FILE__ upar = 2.0/3.0 ppar = 2.0/3.0 mass = 1.0 - vth = sqrt(ppar/(dens*mass)) + vth = sqrt(2.0*ppar/(dens*mass)) for ivz in 1:vz.n for ivr in 1:vr.n vz_val = vz.grid[ivz] From 5bd78ecac315ba9c066697d9ccf55e7a11155daf Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 24 Oct 2023 10:58:57 +0100 Subject: [PATCH 24/49] Increase Lvpa and Lvperp in test_velocity_integrals.jl Need larger boxes to make moments accurate to roughly machine precision. --- test_velocity_integrals.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl index 9bcee9325..25a023539 100644 --- a/test_velocity_integrals.jl +++ b/test_velocity_integrals.jl @@ -18,8 +18,8 @@ if abspath(PROGRAM_FILE) == @__FILE__ ngrid = 17 #number of points per element nelement_local = 20 # number of elements per rank nelement_global = nelement_local # total number of elements - Lvpa = 12.0 #physical box size in reference units - Lvperp = 6.0 #physical box size in reference units + 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" From 4e7dccbca233d8dbc5a9bbf5cd14195b963e5b0b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 24 Oct 2023 16:13:01 +0100 Subject: [PATCH 25/49] Correct file_io so that perpendicular pressure is correctly saved to the output binary. --- src/file_io.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/file_io.jl b/src/file_io.jl index 94f3e4123..e6e132918 100644 --- a/src/file_io.jl +++ b/src/file_io.jl @@ -785,8 +785,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("perpendicular_pressure"), + 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"), From b7a76209f48959c5c8ba10fb5d241fe0e969a049 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Tue, 24 Oct 2023 19:08:32 +0100 Subject: [PATCH 26/49] Fix file name for output plots --- src/post_processing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/post_processing.jl b/src/post_processing.jl index 1b701486d..a625dec01 100644 --- a/src/post_processing.jl +++ b/src/post_processing.jl @@ -1166,9 +1166,9 @@ function analyze_and_plot_data(prefix...; run_index=nothing) 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,perpendicular_pressure,pperp_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + 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,thermal_speed,vthi_sym,"ion",z_global.grid,r_global.grid,time,z_global.n,r_global.n,ntime, + 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", From 88ad6d0610379268bd1ca534d5c4ea3ec2bcefe7 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 25 Oct 2023 13:08:22 +0100 Subject: [PATCH 27/49] Attempt to port manufactured solutions test for Krook operator to new implementation. Suited only for 1D1V with standard drift kinetics for now. Errors suggests an order unity factor difference somewhere in the two implementations. --- src/manufactured_solns.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index ad00ec145..7373bc9d0 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -456,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 @@ -519,6 +522,19 @@ using IfElse 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 + tempi = vthi^2 + nu_ii = nu_krook * densi * (tempi^(-3.0/2.0)) + if vperp_coord.n > 1 + pvth = 3 + else + pvth = 1 + end + FMaxwellian = (densi/vthi^pvth)*exp( -( ( vpa-upari)^2 + vperp^2 )/vthi^2) + Si += -nu_krook*(FMaxwellian - dfni) + end + Source_i = expand_derivatives(Si) From ba4d733dfbd8af9f7929ad6095884d391cfafbb3 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Wed, 25 Oct 2023 13:09:01 +0100 Subject: [PATCH 28/49] Put rfac preceding the numerical dissipation term in r to prevent usage when r.n = 1. --- src/manufactured_solns.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 7373bc9d0..7a88d640a 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -517,7 +517,7 @@ 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)) From e37a0a57024cdf6d97d748c0fcf104f887552df6 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 09:49:03 +0000 Subject: [PATCH 29/49] Clean up comment --- src/time_advance.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/time_advance.jl b/src/time_advance.jl index e357b4025..6802646cf 100644 --- a/src/time_advance.jl +++ b/src/time_advance.jl @@ -1668,7 +1668,7 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, collisions, dt) end - # Add a for Krook collision operoator for ions + # Add Krook collision operator for ions if advance.krook_collisions krook_collisions!(fvec_out.pdf, fvec_in, moments, composition, collisions, vperp, vpa, dt) From 68235462ea2ee98a0bca2c450d8e4ad293ecca8b Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 09:53:01 +0000 Subject: [PATCH 30/49] Modifications to port Krook operator manufactured solutions test to the present branch for the 1D1V case. It is noted that the manufactured solutions test passes when the Krook collision frequency is a constant independent of position. With this in mind, I have modified how the Krook operator is specified in the input file. Now one should specify "krook_collisions_option" ("default" -- use the frequency from the reference parameters and use the spatially dependent factor n/vth^3 -- or "manual" -- use a frequency set in the input file). The frequency can be set in the input file with the "nuii_krook" parameter, if the krook_collisions_option = "manual". The tests and previous functionality are preserved. --- src/input_structs.jl | 2 ++ src/krook_collisions.jl | 16 ++++++++++------ src/manufactured_solns.jl | 8 ++++++-- src/moment_kinetics_input.jl | 11 +++++++---- test/Krook_collisions_tests.jl | 2 +- 5 files changed, 26 insertions(+), 13 deletions(-) diff --git a/src/input_structs.jl b/src/input_structs.jl index 285e42bcd..75b216ca3 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -302,6 +302,8 @@ mutable struct collisions_input constant_ionization_rate::Bool # Coulomb collision rate at the reference density and temperature krook_collision_frequency_prefactor::mk_float + # Coulomb collision rate at the reference density and temperature + krook_collisions_option::String end """ diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 751d6caa1..db96a570c 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -48,14 +48,18 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # 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 collisions.krook_collisions_option == "manual" + cfac = 0.0 + else # default option, collisions.krook_collisions_option == "default" + cfac = 1.0 + end 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) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -69,7 +73,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v 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) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -83,7 +87,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v 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) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -98,7 +102,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v 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) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -117,7 +121,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v else vth_prefactor = 1.0 / vth^3 end - nu_ii = collisions.krook_collision_frequency_prefactor * n * T^(-1.5) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 7a88d640a..2c2f7d05a 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -525,14 +525,18 @@ using IfElse nu_krook = collisions.krook_collision_frequency_prefactor if nu_krook > 0.0 tempi = vthi^2 - nu_ii = nu_krook * densi * (tempi^(-3.0/2.0)) + if collisions.krook_collisions_option == "manual" + nuii_krook = nu_krook + else # default option + nuii_krook = nu_krook * densi * tempi^(-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 += -nu_krook*(FMaxwellian - dfni) + Si += -nuii_krook*(FMaxwellian - dfni) end diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index a15e1c733..f6bd3a93b 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -178,9 +178,12 @@ 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) - 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_collisions_option = get(scan_input, "krook_collisions_option", "reference_parameters") + nuii_krook_default = setup_krook_collisions(reference_parameters) + if collisions.krook_collisions_option == "default" + 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) else collisions.krook_collision_frequency_prefactor = -1.0 end @@ -950,7 +953,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) constant_ionization_rate = false krook_collision_frequency_prefactor = -1.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, - krook_collision_frequency_prefactor) + krook_collision_frequency_prefactor,"default") Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl index 1ebae0805..03c747a20 100644 --- a/test/Krook_collisions_tests.jl +++ b/test/Krook_collisions_tests.jl @@ -99,7 +99,7 @@ test_input_full_f = Dict("n_ion_species" => 1, "evolve_moments_parallel_flow" => false, "evolve_moments_parallel_pressure" => false, "evolve_moments_conservation" => true, - "krook_collisions" => true, + "krook_collisions_option" => "default", "T_e" => 1.0, "initial_density1" => 0.5, "initial_temperature1" => 1.0, From 09cb8fbb2386a8bd58614fd768e8fc2f09bde6f4 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 10:39:11 +0000 Subject: [PATCH 31/49] Allow 2V Krook operator providing that split moments are not used. --- src/krook_collisions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index db96a570c..7c68a7bbe 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -42,7 +42,7 @@ 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 + if vperp.n > 1 && (moments.evolve_density || moments.evolve_upar || moments.evolve_ppar) error("Krook collisions not implemented for 2V case yet") end From 58e7c8ee25c6e9a2d1947eae514af68460e38cb7 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:06:11 +0000 Subject: [PATCH 32/49] Correct update_moments! to use the update_vth! function so that the thermal_speed written to the HDF5/netCDF file is correct. --- src/velocity_moments.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/velocity_moments.jl b/src/velocity_moments.jl index da41a1fea..60ef582e4 100644 --- a/src/velocity_moments.jl +++ b/src/velocity_moments.jl @@ -394,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], @@ -408,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 From 6132fc6d6f566dbba421c4d9a781c8957734bd02 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:08:22 +0000 Subject: [PATCH 33/49] Example input files for MMS test with Krook operator (spatially constant frequency coefficient). --- ...new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 96 +++++++++++++++++++ ..._new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 96 +++++++++++++++++++ 2 files changed, 192 insertions(+) create mode 100644 runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml create mode 100644 runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml 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..5a329be7a --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml @@ -0,0 +1,96 @@ +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 = false +krook_collisions_option = "manual" +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 = "finite_difference" +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..5173fa230 --- /dev/null +++ b/runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml @@ -0,0 +1,96 @@ +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 = false +krook_collisions_option = "manual" +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 = "finite_difference" +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 From 1885c7998f84ad02773c6ac19cb198bd5f12f1f5 Mon Sep 17 00:00:00 2001 From: Michael Hardman Date: Mon, 30 Oct 2023 11:42:23 +0000 Subject: [PATCH 34/49] Change input options for Krook so that default is to not use the Krook operator. Check-in tests still fail for undiagnosed reasons. --- src/moment_kinetics_input.jl | 2 +- test/Krook_collisions_tests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index f6bd3a93b..3ae2eb309 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -180,7 +180,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) collisions.krook_collisions_option = get(scan_input, "krook_collisions_option", "reference_parameters") nuii_krook_default = setup_krook_collisions(reference_parameters) - if collisions.krook_collisions_option == "default" + 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) diff --git a/test/Krook_collisions_tests.jl b/test/Krook_collisions_tests.jl index 03c747a20..62e59c892 100644 --- a/test/Krook_collisions_tests.jl +++ b/test/Krook_collisions_tests.jl @@ -99,7 +99,7 @@ test_input_full_f = Dict("n_ion_species" => 1, "evolve_moments_parallel_flow" => false, "evolve_moments_parallel_pressure" => false, "evolve_moments_conservation" => true, - "krook_collisions_option" => "default", + "krook_collisions_option" => "reference_parameters", "T_e" => 1.0, "initial_density1" => 0.5, "initial_temperature1" => 1.0, From 5647180065ade485d9baf7cbf40e7de6c5a236d4 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 12:39:43 +0000 Subject: [PATCH 35/49] Fix figure saving for 1D plot in makie_post_processing --- src/makie_post_processing.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index a98855b2a..dbcd30da9 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -1647,10 +1647,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 +1686,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] From 7ed8a1c803653c2c3140eba2cba770458d6dc028 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 12:40:35 +0000 Subject: [PATCH 36/49] Pass nvperp to manufactured solutions funcs in makie_post_processing --- src/makie_post_processing.jl | 72 ++++++++++++++++++++++-------------- 1 file changed, 45 insertions(+), 27 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index dbcd30da9..6cf302019 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -313,7 +313,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 @@ -5433,7 +5438,7 @@ 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) @@ -5473,7 +5478,8 @@ function manufactured_solutions_get_field_and_field_sym(run_info, variable_name; 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]] @@ -5560,7 +5566,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) @@ -5570,7 +5577,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 @@ -5808,7 +5815,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) @@ -5833,7 +5840,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) @@ -5857,7 +5864,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) @@ -5885,7 +5892,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) @@ -5992,7 +5999,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 @@ -6012,8 +6020,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) @@ -6032,7 +6040,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) @@ -6053,7 +6062,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) @@ -6174,7 +6184,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 @@ -6194,8 +6205,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) @@ -6214,8 +6225,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) @@ -6236,8 +6247,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) @@ -6259,8 +6270,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) @@ -6319,7 +6330,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 @@ -6338,18 +6349,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 @@ -6373,7 +6390,8 @@ function manufactured_solutions_analysis(run_info; plot_prefix) 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 From a6771e528bc9bd1914b4ae26836395961784092a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 13:07:16 +0000 Subject: [PATCH 37/49] Convert test_velocity_integrals.jl to a CI test --- test/velocity_integral_tests.jl | 151 ++++++++++++++++++++++++++++++++ test_velocity_integrals.jl | 148 ------------------------------- 2 files changed, 151 insertions(+), 148 deletions(-) create mode 100644 test/velocity_integral_tests.jl delete mode 100644 test_velocity_integrals.jl diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl new file mode 100644 index 000000000..bdf6765d9 --- /dev/null +++ b/test/velocity_integral_tests.jl @@ -0,0 +1,151 @@ +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) + vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, + nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) + vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, + nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + # 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) + # output test results + @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) + # output test results + println("") + println("1D Maxwellian") + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + println("") + 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) + # output test results + + println("") + println("biMaxwellian") + @test isapprox(dens_test, dens; atol=atol) + @test isapprox(upar_test, upar; atol=atol) + @test isapprox(ppar_test, ppar; atol=atol) + println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) + println("") + end + end +end + +end # VelocityIntegralTests + +using .VelocityIntegralTests + +VelocityIntegralTests.runtests() diff --git a/test_velocity_integrals.jl b/test_velocity_integrals.jl deleted file mode 100644 index 25a023539..000000000 --- a/test_velocity_integrals.jl +++ /dev/null @@ -1,148 +0,0 @@ -using Printf -using Plots -using LaTeXStrings -using MPI -using Measures - -if abspath(PROGRAM_FILE) == @__FILE__ - using Pkg - Pkg.activate(".") - - import moment_kinetics - 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 - - # 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) - vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, - nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) - # create the coordinate struct 'x' - println("made inputs") - println("vpa: ngrid: ",ngrid," nelement: ",nelement_local, " Lvpa: ",Lvpa) - println("vperp: ngrid: ",ngrid," nelement: ",nelement_local, " Lvperp: ",Lvperp) - 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 - # 2D isotropic Maxwellian test - # 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) - # output test results - println("") - println("Isotropic 2D Maxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("pres_test: ", pres_test, " pres: ", pres, " error: ", abs(pres_test-pres)) - println("") - - ################### - # 1D Maxwellian test - - 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) - # output test results - println("") - println("1D Maxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) - println("") - - ################### - # biMaxwellian test - - # 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) - # output test results - - println("") - println("biMaxwellian") - println("dens_test: ", dens_test, " dens: ", dens, " error: ", abs(dens_test-dens)) - println("upar_test: ", upar_test, " upar: ", upar, " error: ", abs(upar_test-upar)) - println("ppar_test: ", ppar_test, " ppar: ", ppar, " error: ", abs(ppar_test-ppar)) - println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) - println("") -end From 6884d86ebca9f58479b4f3648c43c509f255b54e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 14:42:58 +0000 Subject: [PATCH 38/49] Deactivate Krook collisions by default --- src/moment_kinetics_input.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 4e75fc7b7..1132b961f 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -168,14 +168,18 @@ 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", "reference_parameters") + 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) - else + 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 @@ -957,7 +961,7 @@ function load_defaults(n_ion_species, n_neutral_species, electron_physics) constant_ionization_rate = false krook_collision_frequency_prefactor = -1.0 collisions = collisions_input(charge_exchange, ionization, constant_ionization_rate, - krook_collision_frequency_prefactor,"default") + krook_collision_frequency_prefactor,"none") Bzed = 1.0 # magnetic field component along z Bmag = 1.0 # magnetic field strength From 92a4e29ad9fa130101738c3033eb5625958e0d26 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 14:44:07 +0000 Subject: [PATCH 39/49] Pass 'element_spacing_option' to coords in velocity_integral_tests.jl This is needed since 'sqrt spacing' option was added. Also include velocity_integral_tests.jl in runtests.jl so that it is included in the CI tests. Also remove some println() statements that had been left in by mistake. --- test/runtests.jl | 1 + test/velocity_integral_tests.jl | 26 +++++++++++--------------- 2 files changed, 12 insertions(+), 15 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 17ba5030d..4bab7c641 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,6 +7,7 @@ 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")) diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl index bdf6765d9..dc84800f0 100644 --- a/test/velocity_integral_tests.jl +++ b/test/velocity_integral_tests.jl @@ -32,14 +32,17 @@ function runtests() 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) - vz_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vpa_input = grid_input("vpa", ngrid, nelement_global, nelement_local, - nrank, irank, Lvpa, discretization, fd_option, bc, adv_input,comm) - vperp_input = grid_input("vperp", ngrid, nelement_global, nelement_local, - nrank, irank, Lvperp, discretization, fd_option, bc, adv_input,comm) + 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) @@ -100,12 +103,9 @@ function runtests() upar_test = get_upar(dfn1D,vz,vr,dens_test) ppar_test = get_ppar(dfn1D,vz,vr,upar_test) # output test results - println("") - println("1D Maxwellian") @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) - println("") end @testset "biMaxwellian" begin @@ -133,13 +133,9 @@ function runtests() pperp_test = get_pperp(dfn,vpa,vperp) # output test results - println("") - println("biMaxwellian") @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) - println("pperp_test: ", pperp_test, " pperp: ", pperp, " error: ", abs(pperp_test-pperp)) - println("") end end end From afe37095ef4e1550b2d578ffb6e16d1b6c53441e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 29 Sep 2023 10:16:32 +0100 Subject: [PATCH 40/49] Better values for MPI variables when ignore_MPI=true in mk_input() Set irank/nrank as if running in serial, to avoid errors from set_element_scale_and_shift()). --- src/moment_kinetics_input.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/moment_kinetics_input.jl b/src/moment_kinetics_input.jl index 1132b961f..72e23d6f5 100644 --- a/src/moment_kinetics_input.jl +++ b/src/moment_kinetics_input.jl @@ -372,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) From 5d5ddafaa43c4bcf6819cfca7a3ac1e4b7368c95 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:17:46 +0000 Subject: [PATCH 41/49] Tidy up comments --- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 2 -- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 2 -- src/input_structs.jl | 2 +- src/krook_collisions.jl | 2 +- test/velocity_integral_tests.jl | 3 --- 5 files changed, 2 insertions(+), 9 deletions(-) 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 index 5a329be7a..7265dd542 100644 --- 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 @@ -31,7 +31,6 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -#krook_collisions = false krook_collisions_option = "manual" nuii_krook = 1.0 nstep = 2000 @@ -60,7 +59,6 @@ vperp_ngrid = 1 vperp_nelement = 1 vperp_L = 6.0 vperp_bc = "periodic" -#vperp_discretization = "finite_difference" vperp_discretization = "chebyshev_pseudospectral" vz_ngrid = 17 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 index 5173fa230..4d35a8a77 100644 --- 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 @@ -31,7 +31,6 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -#krook_collisions = false krook_collisions_option = "manual" nuii_krook = 1.0 nstep = 2000 @@ -60,7 +59,6 @@ vperp_ngrid = 17 vperp_nelement = 8 vperp_L = 6.0 vperp_bc = "periodic" -#vperp_discretization = "finite_difference" vperp_discretization = "chebyshev_pseudospectral" vz_ngrid = 17 diff --git a/src/input_structs.jl b/src/input_structs.jl index f7734a982..be8a19969 100644 --- a/src/input_structs.jl +++ b/src/input_structs.jl @@ -306,7 +306,7 @@ mutable struct collisions_input constant_ionization_rate::Bool # Coulomb collision rate at the reference density and temperature krook_collision_frequency_prefactor::mk_float - # Coulomb collision rate at the reference density and temperature + # Setting to switch between different options for Krook collision operator krook_collisions_option::String end diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index b77db5bf1..f6ceec9c4 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -50,7 +50,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # normalized by sqrt(pi) (see velocity_moments.integrate_over_vspace). if collisions.krook_collisions_option == "manual" cfac = 0.0 - else # default option, collisions.krook_collisions_option == "default" + else # default option, collisions.krook_collisions_option == "reference_parameters" cfac = 1.0 end if moments.evolve_ppar && moments.evolve_upar diff --git a/test/velocity_integral_tests.jl b/test/velocity_integral_tests.jl index dc84800f0..07d0c5027 100644 --- a/test/velocity_integral_tests.jl +++ b/test/velocity_integral_tests.jl @@ -81,7 +81,6 @@ function runtests() ppar_test = get_ppar(dfn,vpa,vperp,upar_test) pperp_test = get_pperp(dfn,vpa,vperp) pres_test = pressure(ppar_test,pperp_test) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(pres_test, pres; atol=atol) @@ -102,7 +101,6 @@ function runtests() 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) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) @test isapprox(ppar_test, ppar; atol=atol) @@ -131,7 +129,6 @@ function runtests() 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) - # output test results @test isapprox(dens_test, dens; atol=atol) @test isapprox(upar_test, upar; atol=atol) From 8a5e2986ab21e5159ebf25fb0590cf2b5b422d61 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:58:01 +0000 Subject: [PATCH 42/49] Make clearer that Ti_over_Tref is not normalised temperature --- src/manufactured_solns.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/manufactured_solns.jl b/src/manufactured_solns.jl index 2c2f7d05a..df267400e 100644 --- a/src/manufactured_solns.jl +++ b/src/manufactured_solns.jl @@ -524,11 +524,11 @@ using IfElse end nu_krook = collisions.krook_collision_frequency_prefactor if nu_krook > 0.0 - tempi = vthi^2 + Ti_over_Tref = vthi^2 if collisions.krook_collisions_option == "manual" nuii_krook = nu_krook else # default option - nuii_krook = nu_krook * densi * tempi^(-1.5) + nuii_krook = nu_krook * densi * Ti_over_Tref^(-1.5) end if vperp_coord.n > 1 pvth = 3 From 4c9e8109841f8f6c0dd27190f38352443b73c07c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 19:58:58 +0000 Subject: [PATCH 43/49] Krook MMS test inputs use spatially-varying collision frequency This option now passes the tests, so might as well make it the suggested option. --- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_16_vperp_1_krook.toml | 2 +- runs/1D-wall_MMS_new_nel_r_1_z_16_vpa_8_vperp_8_krook.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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 index 7265dd542..e9d44cc5b 100644 --- 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 @@ -31,7 +31,7 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -krook_collisions_option = "manual" +krook_collisions_option = "reference_parameters" nuii_krook = 1.0 nstep = 2000 dt = 0.0005 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 index 4d35a8a77..ae02bfd9c 100644 --- 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 @@ -31,7 +31,7 @@ z_IC_temperature_amplitude2 = 0.0 z_IC_temperature_phase2 = 0.0 charge_exchange_frequency = 0.0 ionization_frequency = 0.0 -krook_collisions_option = "manual" +krook_collisions_option = "reference_parameters" nuii_krook = 1.0 nstep = 2000 dt = 0.0005 From 1654f7ce23b75b3f11df5c78816f9749ece30d94 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 30 Oct 2023 20:05:55 +0000 Subject: [PATCH 44/49] Make clearer that T_over_Tref is not normalised temperature ...with the current normalisations used by the code, which normalise pressure by `mref*Nref*cref^2 = 2*Nref*Tref`. --- src/krook_collisions.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index f6ceec9c4..91a8e09a0 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -58,8 +58,8 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # 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) * cfac + (1.0 - cfac)) + T_over_Tref = (moments.charged.vth[iz,ir,is])^2 + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -72,8 +72,8 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @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) * cfac + (1.0 - cfac)) + T_over_Tref = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -86,8 +86,8 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @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) * cfac + (1.0 - cfac)) + T_over_Tref = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -101,8 +101,8 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @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) * cfac + (1.0 - cfac)) + T_over_Tref = vth^2 + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -115,13 +115,13 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @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 + T_over_Tref = 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) * cfac + (1.0 - cfac)) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] From 45fd04f57c220c62bbca5d6046e20fbb0cddff6a Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 31 Oct 2023 09:35:15 +0000 Subject: [PATCH 45/49] Define Krook collision frequency in terms of cref rather than Tref This means that to calculate the normalised, spatially-varying collision frequency we multiply by (the normalised) n/vth^3, and do not need to worry about potentially changing normalisation of T (to mref*cref^2 or to Tref). --- src/krook_collisions.jl | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 91a8e09a0..0d65b5c06 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -15,6 +15,7 @@ function setup_krook_collisions(reference_params) Tref = reference_params.Tref mref = reference_params.mref timeref = reference_params.timeref + cref = reference_params.cref Nref_per_cm3 = Nref * 1.0e-6 @@ -24,11 +25,8 @@ function setup_krook_collisions(reference_params) # 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 + (4.0 * π * epsilon0^2 * mref^2 * cref^3) # s^-1 nu_ii0 = nu_ii0_per_s * timeref return nu_ii0 @@ -58,8 +56,8 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # 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_over_Tref = (moments.charged.vth[iz,ir,is])^2 - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) + vth = moments.charged.vth[iz,ir,is] + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -72,8 +70,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - T_over_Tref = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -86,8 +83,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - T_over_Tref = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -101,8 +97,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - T_over_Tref = vth^2 - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * T_over_Tref^(-1.5) * cfac + (1.0 - cfac)) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] @@ -115,13 +110,12 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - T_over_Tref = 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_over_Tref^(-1.5) * cfac + (1.0 - cfac)) + nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) @loop_vperp_vpa ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir,is] -= dt * nu_ii * (fvec_in.pdf[ivpa,ivperp,iz,ir,is] From 37856fa45ab9a3674ef95a04050b2ed491861d51 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 31 Oct 2023 09:56:04 +0000 Subject: [PATCH 46/49] Move calculation of Krook collision frequency to a function Makes it easier to reuse in future (e.g. for post processing). --- src/krook_collisions.jl | 49 ++++++++++++++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 10 deletions(-) diff --git a/src/krook_collisions.jl b/src/krook_collisions.jl index 0d65b5c06..fe6477098 100644 --- a/src/krook_collisions.jl +++ b/src/krook_collisions.jl @@ -2,6 +2,8 @@ """ module krook_collisions +export setup_krook_collisions, get_collision_frequency, krook_collisions! + using ..constants: epsilon0, proton_charge using ..looping @@ -32,6 +34,38 @@ function setup_krook_collisions(reference_params) 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 @@ -46,18 +80,13 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v # 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 collisions.krook_collisions_option == "manual" - cfac = 0.0 - else # default option, collisions.krook_collisions_option == "reference_parameters" - cfac = 1.0 - end 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 = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) + 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] @@ -70,7 +99,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) + 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] @@ -83,7 +112,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) + 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] @@ -97,7 +126,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v @loop_s_r_z is ir iz begin n = fvec_in.density[iz,ir,is] vth = moments.charged.vth[iz,ir,is] - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) + 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] @@ -115,7 +144,7 @@ function krook_collisions!(pdf_out, fvec_in, moments, composition, collisions, v else vth_prefactor = 1.0 / vth^3 end - nu_ii = collisions.krook_collision_frequency_prefactor * ( n * vth^(-3) * cfac + (1.0 - cfac)) + 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] From eff26477311deb7ba20ea66dcdbd268d8890a79c Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 31 Oct 2023 10:03:28 +0000 Subject: [PATCH 47/49] Make plots of collision frequency if Krook collisions are used --- src/makie_post_processing.jl | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index e6835b1d3..fb6c7c5e8 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -20,6 +20,7 @@ using ..array_allocation: allocate_float using ..coordinates: define_coordinate 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 @@ -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") @@ -1300,6 +1302,16 @@ function plots_for_variable(run_info, variable_name; plot_prefix, is_1D=false, vth = Tuple(postproc_load_variable(ri, "thermal_speed") for ri ∈ run_info) variable = Tuple(v.^2 for v ∈ vth) + elseif variable_name == "collision_frequency" + if all(ri.collisions.krook_collisions_option == "none" for ri ∈ run_info) + return nothing + end + n = Tuple(postproc_load_variable(ri, "density") + for ri ∈ run_info) + vth = Tuple(postproc_load_variable(ri, "thermal_speed") + for ri ∈ run_info) + variable = Tuple(get_collision_frequency(ri.collisions, this_n, this_vth) + for (ri, this_n, this_vth) ∈ zip(run_info, n, vth)) elseif variable_name == "temperature_neutral" vth = Tuple(postproc_load_variable(ri, "thermal_speed_neutral") for ri ∈ run_info) From 326efbb2cbc28ebe41a306548575570b0ea3e3c5 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 31 Oct 2023 10:11:46 +0000 Subject: [PATCH 48/49] Fix coordinate creation in makie_post_processing Need to add an `element_spacing_option` argument for the dummy coordinate input. --- src/makie_post_processing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index fb6c7c5e8..b6841b596 100644 --- a/src/makie_post_processing.jl +++ b/src/makie_post_processing.jl @@ -817,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) From c7735e160c3082053fc1d2ba3b14ea5e201d6ec3 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 31 Oct 2023 10:22:25 +0000 Subject: [PATCH 49/49] Handle derived variables in get_variable() in makie_post_processing ...rather than handling ad-hoc in `plots_for_variable()`. This is cleaner, and makes it possible when using makie_post_processing interactively to easily get "temperature", "collision_frequency", etc. --- src/makie_post_processing.jl | 67 +++++++++++++++++++++++------------- 1 file changed, 44 insertions(+), 23 deletions(-) diff --git a/src/makie_post_processing.jl b/src/makie_post_processing.jl index b6841b596..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 @@ -1184,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} @@ -1293,33 +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 == "collision_frequency" - if all(ri.collisions.krook_collisions_option == "none" for ri ∈ run_info) - return nothing - end - n = Tuple(postproc_load_variable(ri, "density") - for ri ∈ run_info) - vth = Tuple(postproc_load_variable(ri, "thermal_speed") - for ri ∈ run_info) - variable = Tuple(get_collision_frequency(ri.collisions, this_n, this_vth) - for (ri, this_n, this_vth) ∈ zip(run_info, n, 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 ||