From af494328893af9eb56d75d1c85b84e0e09c85242 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Fri, 6 Sep 2024 14:02:17 +0100 Subject: [PATCH 01/23] Split get_settings function for neutrals and ions, in preparation for multiple source profiles --- moment_kinetics/src/external_sources.jl | 99 +++++++++++++++++-------- 1 file changed, 69 insertions(+), 30 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 98fe4169c..a264a1191 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -40,13 +40,13 @@ Returns a NamedTuple `(ion=ion_source_settings, neutral=neutral_source_settings) containing two NamedTuples of settings. """ function setup_external_sources!(input_dict, r, z, electron_physics) - function get_settings(neutrals) + function get_settings_ions() input = set_defaults_and_check_section!( - input_dict, neutrals ? "neutral_source" : "ion_source"; + input_dict, "ion_source"; active=false, source_strength=1.0, source_n=1.0, - source_T=neutrals ? get(input_dict, "T_wall", 1.0) : 1.0, + source_T=1.0, source_v0=0.0, # birth speed for "alphas" option source_vpa0=0.0, # birth vpa for "beam" option source_vperp0=0.0, # birth vperp for "beam" option @@ -137,7 +137,65 @@ function setup_external_sources!(input_dict, r, z, electron_physics) else PI_density_target_rank = nothing end - elseif neutrals && input["source_type"] == "recycling" + elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses") + PI_density_target = nothing + PI_controller_amplitude = nothing + controller_source_profile = nothing + PI_density_target_ir = nothing + PI_density_target_iz = nothing + PI_density_target_rank = nothing + else + error("Unrecognised ion source_type=$(input["source_type"])." + * "Possible values are: Maxwellian, density_profile_control, " + * "density_midpoint_control") + end + + return (; (Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude=r_amplitude, + z_amplitude=z_amplitude, PI_density_target=PI_density_target, + PI_controller_amplitude=PI_controller_amplitude, + controller_source_profile=controller_source_profile, + PI_density_target_ir=PI_density_target_ir, + PI_density_target_iz=PI_density_target_iz, + PI_density_target_rank=PI_density_target_rank) + end + + function get_settings_neutrals() + input = set_defaults_and_check_section!( + input_dict, "neutral_source"; + active=false, + source_strength=1.0, + source_n=1.0, + source_T=get(input_dict, "T_wall", 1.0), + source_v0=0.0, # birth speed for "alphas" option + source_vpa0=0.0, # birth vpa for "beam" option + source_vperp0=0.0, # birth vperp for "beam" option + sink_strength=1.0, # strength of sink in "alphas-with-losses" & "beam-with-losses" option + sink_vth=0.0, # thermal speed for sink in "alphas-with-losses" & "beam-with-losses" option + r_profile="constant", + r_width=1.0, + r_relative_minimum=0.0, + z_profile="constant", + z_width=1.0, + z_relative_minimum=0.0, + source_type="recycling", # "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" + PI_density_controller_P=0.0, + PI_density_controller_I=0.0, + PI_density_target_amplitude=1.0, + PI_density_target_r_profile="constant", + PI_density_target_r_width=1.0, + PI_density_target_r_relative_minimum=0.0, + PI_density_target_z_profile="constant", + PI_density_target_z_width=1.0, + PI_density_target_z_relative_minimum=0.0, + recycling_controller_fraction=0.5, + ) + + r_amplitude = get_source_profile(input["r_profile"], input["r_width"], + input["r_relative_minimum"], r) + z_amplitude = get_source_profile(input["z_profile"], input["z_width"], + input["z_relative_minimum"], z) + + if input["source_type"] == "recycling" recycling = input["recycling_controller_fraction"] if recycling ≤ 0.0 # Don't allow 0.0 as this is the default value, but makes no sense to have @@ -171,34 +229,15 @@ function setup_external_sources!(input_dict, r, z, electron_physics) end controller_source_profile ./= controller_source_integral end - - PI_density_target = nothing - PI_controller_amplitude = nothing - PI_density_target_ir = nothing - PI_density_target_iz = nothing - PI_density_target_rank = nothing - elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses") - PI_density_target = nothing - PI_controller_amplitude = nothing - controller_source_profile = nothing - PI_density_target_ir = nothing - PI_density_target_iz = nothing - PI_density_target_rank = nothing else - error("Unrecognised source_type=$(input["source_type"])." - * "Possible values are: Maxwellian, density_profile_control, " - * "density_midpoint_control, recycling (for neutrals only)") + error("Unrecognised neutral source_type=$(input["source_type"])." + * "Possible values are: recycling ") end - return (; (Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude=r_amplitude, - z_amplitude=z_amplitude, PI_density_target=PI_density_target, - PI_controller_amplitude=PI_controller_amplitude, - controller_source_profile=controller_source_profile, - PI_density_target_ir=PI_density_target_ir, - PI_density_target_iz=PI_density_target_iz, - PI_density_target_rank=PI_density_target_rank) + return (; (Symbol(k)=>v for (k,v) ∈ input)..., + r_amplitude=r_amplitude, z_amplitude=z_amplitude, + controller_source_profile=controller_source_profile) end - function get_electron_settings(ion_settings) # Note most settings for the electron source are copied from the ion source, # because we require that the particle sources are the same for ions and @@ -226,14 +265,14 @@ function setup_external_sources!(input_dict, r, z, electron_physics) source_type=ion_settings.source_type) end - ion_settings = get_settings(false) + ion_settings = get_settings_ions() if electron_physics ∈ (braginskii_fluid, kinetic_electrons, kinetic_electrons_with_temperature_equation) electron_settings = get_electron_settings(ion_settings) else electron_settings = (active=false,) end - neutral_settings = get_settings(true) + neutral_settings = get_settings_neutrals() return (ion=ion_settings, electron=electron_settings, neutral=neutral_settings) end From 0c67d262a5e7ff1d8ca908bb14c3da247661216b Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Fri, 6 Sep 2024 20:21:10 +0100 Subject: [PATCH 02/23] Convert form of ion_settings, electron_settings and neutral_settings from NamedTuple to structs --- moment_kinetics/src/external_sources.jl | 118 +++++++++++++++---- moment_kinetics/src/input_structs.jl | 143 ++++++++++++++++++++++++ 2 files changed, 239 insertions(+), 22 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index a264a1191..5295176f3 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -149,14 +149,11 @@ function setup_external_sources!(input_dict, r, z, electron_physics) * "Possible values are: Maxwellian, density_profile_control, " * "density_midpoint_control") end - - return (; (Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude=r_amplitude, + println(input) + return ion_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, z_amplitude=z_amplitude, PI_density_target=PI_density_target, - PI_controller_amplitude=PI_controller_amplitude, - controller_source_profile=controller_source_profile, - PI_density_target_ir=PI_density_target_ir, - PI_density_target_iz=PI_density_target_iz, - PI_density_target_rank=PI_density_target_rank) + PI_controller_amplitude, controller_source_profile, + PI_density_target_ir, PI_density_target_iz, PI_density_target_rank) end function get_settings_neutrals() @@ -177,7 +174,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics) z_profile="constant", z_width=1.0, z_relative_minimum=0.0, - source_type="recycling", # "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" + source_type="Maxwellian", # "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" PI_density_controller_P=0.0, PI_density_controller_I=0.0, PI_density_target_amplitude=1.0, @@ -187,15 +184,76 @@ function setup_external_sources!(input_dict, r, z, electron_physics) PI_density_target_z_profile="constant", PI_density_target_z_width=1.0, PI_density_target_z_relative_minimum=0.0, - recycling_controller_fraction=0.5, + recycling_controller_fraction=0.0, ) r_amplitude = get_source_profile(input["r_profile"], input["r_width"], input["r_relative_minimum"], r) z_amplitude = get_source_profile(input["z_profile"], input["z_width"], input["z_relative_minimum"], z) + if input["source_type"] == "density_profile_control" + PI_density_target_amplitude = input["PI_density_target_amplitude"] + PI_density_target_r_factor = + get_source_profile(input["PI_density_target_r_profile"], + input["PI_density_target_r_width"], + input["PI_density_target_r_relative_minimum"], r) + PI_density_target_z_factor = + get_source_profile(input["PI_density_target_z_profile"], + input["PI_density_target_z_width"], + input["PI_density_target_z_relative_minimum"], z) + PI_density_target = allocate_shared_float(z.n,r.n) + @serial_region begin + for ir ∈ 1:r.n, iz ∈ 1:z.n + PI_density_target[iz,ir] = + PI_density_target_amplitude * PI_density_target_r_factor[ir] * + PI_density_target_z_factor[iz] + end + end + PI_controller_amplitude = nothing + controller_source_profile = nothing + PI_density_target_ir = nothing + PI_density_target_iz = nothing + PI_density_target_rank = nothing + elseif input["source_type"] == "density_midpoint_control" + PI_density_target = input["PI_density_target_amplitude"] + + if comm_block[] != MPI.COMM_NULL + PI_controller_amplitude = allocate_shared_float(1) + controller_source_profile = allocate_shared_float(z.n, r.n) + else + PI_controller_amplitude = allocate_float(1) + controller_source_profile = allocate_float(z.n, r.n) + end + for ir ∈ 1:r.n, iz ∈ 1:z.n + controller_source_profile[iz,ir] = r_amplitude[ir] * z_amplitude[iz] + end - if input["source_type"] == "recycling" + # Find the indices, and process rank of the point at r=0, z=0. + # The result of findfirst() will be `nothing` if the point was not found. + PI_density_target_ir = findfirst(x->abs(x)<1.e-14, r.grid) + PI_density_target_iz = findfirst(x->abs(x)<1.e-14, z.grid) + if block_rank[] == 0 + # Only need to do communications from the root process of each + # shared-memory block + if PI_density_target_ir !== nothing && PI_density_target_iz !== nothing + PI_density_target_rank = iblock_index[] + else + PI_density_target_rank = 0 + end + if comm_inter_block[] != MPI.COMM_NULL + PI_density_target_rank = MPI.Allreduce(PI_density_target_rank, +, + comm_inter_block[]) + end + if PI_density_target_rank == 0 && iblock_index[] == 0 && + (PI_density_target_ir === nothing || + PI_density_target_iz === nothing) + error("No grid point with r=0 and z=0 was found for the " + * "'density_midpoint' controller.") + end + else + PI_density_target_rank = nothing + end + elseif input["source_type"] == "recycling" recycling = input["recycling_controller_fraction"] if recycling ≤ 0.0 # Don't allow 0.0 as this is the default value, but makes no sense to have @@ -229,16 +287,31 @@ function setup_external_sources!(input_dict, r, z, electron_physics) end controller_source_profile ./= controller_source_integral end + + PI_density_target = nothing + PI_controller_amplitude = nothing + PI_density_target_ir = nothing + PI_density_target_iz = nothing + PI_density_target_rank = nothing + elseif input["source_type"] ∈ ("Maxwellian", "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses") + PI_density_target = nothing + PI_controller_amplitude = nothing + controller_source_profile = nothing + PI_density_target_ir = nothing + PI_density_target_iz = nothing + PI_density_target_rank = nothing else error("Unrecognised neutral source_type=$(input["source_type"])." - * "Possible values are: recycling ") + * "Possible values are: Maxwellian, density_profile_control, " + * "density_midpoint_control, recycling (for neutrals only)") end - return (; (Symbol(k)=>v for (k,v) ∈ input)..., - r_amplitude=r_amplitude, z_amplitude=z_amplitude, - controller_source_profile=controller_source_profile) + return neutral_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, + z_amplitude=z_amplitude, PI_density_target=PI_density_target, + PI_controller_amplitude, controller_source_profile, + PI_density_target_ir, PI_density_target_iz, PI_density_target_rank) end - function get_electron_settings(ion_settings) + function get_settings_electrons(ion_settings) # Note most settings for the electron source are copied from the ion source, # because we require that the particle sources are the same for ions and # electrons. `source_T` can be set independently, and when using @@ -259,21 +332,22 @@ function setup_external_sources!(input_dict, r, z, electron_physics) end input["source_strength"] = ion_settings.source_strength end - return (; (Symbol(k)=>v for (k,v) ∈ input)..., active=ion_settings.active, - r_amplitude=ion_settings.r_amplitude, - z_amplitude=ion_settings.z_amplitude, - source_type=ion_settings.source_type) + return electron_source_profile(input["source_strength"], input["source_T"], + ion_settings.active, ion_settings.r_amplitude, + ion_settings.z_amplitude, ion_settings.source_type) end ion_settings = get_settings_ions() if electron_physics ∈ (braginskii_fluid, kinetic_electrons, kinetic_electrons_with_temperature_equation) - electron_settings = get_electron_settings(ion_settings) + electron_settings = get_settings_electrons(ion_settings) else - electron_settings = (active=false,) + electron_settings = electron_source_profile(ion_settings.source_strength, + ion_settings.source_T, false, ion_settings.r_amplitude, + ion_settings.z_amplitude, ion_settings.source_type) end neutral_settings = get_settings_neutrals() - + return (ion=ion_settings, electron=electron_settings, neutral=neutral_settings) end diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index 30bc39cb2..ebf735495 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -12,6 +12,7 @@ export mk_to_toml export species_parameters, species_parameters_mutable export species_composition export drive_input, drive_input_mutable +export ion_source_profile, electron_source_profile, neutral_source_profile export collisions_input, krook_collisions_input, fkpl_collisions_input export io_input export pp_input @@ -383,6 +384,148 @@ struct drive_input force_Er_zero_at_wall::Bool end +""" +Source profile structs for ions and electrons which allows them to have any number +of different sources (from wall perhaps, superposition of core sources, etc.). These +sources are then contained in a vector of structs. + +Since the ion source must be the same as the electron source in all respects (apart +from possibly a different electron temperature or source strength), the electron +vector of source profile structs will be a kind of mirror of the ion vector of structs. +""" +Base.@kwdef struct ion_source_profile + # struct containing source profile data for ions + # is the source active or not + active::Bool + # An overall multiplier for the strength (i.e. 0.0 would turn off the source) + source_strength::mk_float + # For use with "energy" option, Krook source (...) + source_n::mk_float + # Temperature of source (variation along z can be introduced later) + source_T::mk_float + # birth speed for "alphas" option + source_v0::mk_float + # birth vpa for "beam" option + source_vpa0::mk_float + # birth vperp for "beam" option + source_vperp0::mk_float + # strength of sink in "alphas-with-losses" & "beam-with-losses" option + sink_strength::mk_float + # thermal speed for sink in "alphas-with-losses" & "beam-with-losses" option + sink_vth::mk_float + # profile for source in r ('constant' or 'gaussian' or 'parabolic' etc.) + r_profile::String + # width of source in r (doesn't apply to constant profile) + r_width::mk_float + # relative minimum of source in r, acts as the baseline + r_relative_minimum::mk_float + # profile for source in z ('constant' or 'gaussian' or 'parabolic' etc.) + z_profile::String + # width of source in z (doesn't apply to constant profile) + z_width::mk_float + # relative minimum of source in z, acts as the baseline (so you can elevate + # your gaussian source above 0, for example) + z_relative_minimum::mk_float + # velocity profile for the source, Maxwellian would be default, can have beams too + source_type::String #"Maxwellian" "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" + # proportional integral controllers - for when you want to find a source profile that matches + # a target density + PI_density_controller_P::mk_float + PI_density_controller_I::mk_float + PI_density_target_amplitude::mk_float + PI_density_target_r_profile::String + PI_density_target_r_width::mk_float + PI_density_target_r_relative_minimum::mk_float + PI_density_target_z_profile::String + PI_density_target_z_width::mk_float + PI_density_target_z_relative_minimum::mk_float + recycling_controller_fraction::mk_float + # r_amplitude through the r coordinate (in 1D this can just be set to 1.0) + r_amplitude::Vector{mk_float} + # z_amplitude through the z coordinate, which will have your gaussian profile, + # constant profile, parabolic, etc.. + z_amplitude::Vector{mk_float} + PI_density_target::Union{mk_float, Nothing, MPISharedArray{mk_float,2}} + PI_controller_amplitude::Union{Nothing, MPISharedArray{mk_float,1}} + controller_source_profile::Union{Nothing, MPISharedArray{mk_float,2}, Array{mk_float, 2}} + PI_density_target_ir::Union{mk_int, Nothing} + PI_density_target_iz::Union{mk_int, Nothing} + PI_density_target_rank::Union{mk_int, Nothing} #possibly this should have Int64 as well, + # in the event that the code is running with mk_int = Int32 but the rank is set to 0::Int64 +end + +Base.@kwdef struct electron_source_profile + # most of the electron parameters must be the same as for ions, so only + # source strength (in the case of an ion energy source) and source Temperature + # can be different. The other four are set by the ion source profile. + source_strength::mk_float + source_T::mk_float + active::Bool + r_amplitude::Vector{mk_float} + z_amplitude::Vector{mk_float} + source_type::String +end + +Base.@kwdef struct neutral_source_profile + # struct containing source profile data for neutrals + # is the source active or not + active::Bool + # An overall multiplier for the strength (i.e. 0.0 would turn off the source) + source_strength::mk_float + # For use with "energy" option, Krook source (...) + source_n::mk_float + # Temperature of source (variation along z can be introduced later) + source_T::mk_float + # birth speed for "alphas" option + source_v0::mk_float + # birth vpa for "beam" option + source_vpa0::mk_float + # birth vperp for "beam" option + source_vperp0::mk_float + # strength of sink in "alphas-with-losses" & "beam-with-losses" option + sink_strength::mk_float + # thermal speed for sink in "alphas-with-losses" & "beam-with-losses" option + sink_vth::mk_float + # profile for source in r ('constant' or 'gaussian' or 'parabolic' etc.) + r_profile::String + # width of source in r (doesn't apply to constant profile) + r_width::mk_float + # relative minimum of source in r, acts as the baseline + r_relative_minimum::mk_float + # profile for source in z ('constant' or 'gaussian' or 'parabolic' etc.) + z_profile::String + # width of source in z (doesn't apply to constant profile) + z_width::mk_float + # relative minimum of source in z, acts as the baseline (so you can elevate + # your gaussian source above 0, for example) + z_relative_minimum::mk_float + # velocity profile for the source, Maxwellian would be default, can have beams too + source_type::String #"Maxwellian" "energy", "alphas", "alphas-with-losses", "beam", "beam-with-losses" + # proportional integral controllers - for when you want to find a source profile that matches + # a target density + PI_density_controller_P::mk_float + PI_density_controller_I::mk_float + PI_density_target_amplitude::mk_float + PI_density_target_r_profile::String + PI_density_target_r_width::mk_float + PI_density_target_r_relative_minimum::mk_float + PI_density_target_z_profile::String + PI_density_target_z_width::mk_float + PI_density_target_z_relative_minimum::mk_float + recycling_controller_fraction::mk_float + # r_amplitude through the r coordinate (in 1D this can just be set to 1.0) + r_amplitude::Vector{mk_float} + # z_amplitude through the z coordinate, which will have your gaussian profile, + # constant profile, parabolic, etc.. + z_amplitude::Vector{mk_float} + PI_density_target::Union{mk_float, Nothing, MPISharedArray{mk_float,2}} + PI_controller_amplitude::Union{Nothing, MPISharedArray{mk_float,1}} + controller_source_profile::Union{Nothing, MPISharedArray{mk_float,2}, Array{mk_float, 2}} + PI_density_target_ir::Union{mk_int, Nothing} + PI_density_target_iz::Union{mk_int, Nothing} + PI_density_target_rank::Union{mk_int, Nothing} #possibly this should have Int64 as well, + # in the event that the code is running with mk_int = Int32 but the rank is set to 0::Int64 +end """ Structs set up for the collision operators so far in use. These will each be contained in the main collisions_input struct below, as substructs. From 372988aa7a76b065b121584eb3f57657e14ca2ab Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Fri, 6 Sep 2024 20:28:15 +0100 Subject: [PATCH 03/23] Add options to ion and neutral source profile type error --- moment_kinetics/src/external_sources.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 5295176f3..dcaca6905 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -147,7 +147,8 @@ function setup_external_sources!(input_dict, r, z, electron_physics) else error("Unrecognised ion source_type=$(input["source_type"])." * "Possible values are: Maxwellian, density_profile_control, " - * "density_midpoint_control") + * "density_midpoint_control, energy, alphas, alphas-with-losses, " + * "beam, beam-with-losses") end println(input) return ion_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, @@ -303,7 +304,8 @@ function setup_external_sources!(input_dict, r, z, electron_physics) else error("Unrecognised neutral source_type=$(input["source_type"])." * "Possible values are: Maxwellian, density_profile_control, " - * "density_midpoint_control, recycling (for neutrals only)") + * "density_midpoint_control, energy, alphas, alphas-with-losses, " + * "beam, beam-with-losses, recycling (for neutrals only)") end return neutral_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, From 4a74181e55d116dfed18daf1f08f08edbd6128dc Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Fri, 6 Sep 2024 21:44:46 +0100 Subject: [PATCH 04/23] Change names of structs from 'profile' to 'data' --- moment_kinetics/src/input_structs.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/moment_kinetics/src/input_structs.jl b/moment_kinetics/src/input_structs.jl index ebf735495..284cf147a 100644 --- a/moment_kinetics/src/input_structs.jl +++ b/moment_kinetics/src/input_structs.jl @@ -12,7 +12,7 @@ export mk_to_toml export species_parameters, species_parameters_mutable export species_composition export drive_input, drive_input_mutable -export ion_source_profile, electron_source_profile, neutral_source_profile +export ion_source_data, electron_source_data, neutral_source_data export collisions_input, krook_collisions_input, fkpl_collisions_input export io_input export pp_input @@ -393,8 +393,8 @@ Since the ion source must be the same as the electron source in all respects (ap from possibly a different electron temperature or source strength), the electron vector of source profile structs will be a kind of mirror of the ion vector of structs. """ -Base.@kwdef struct ion_source_profile - # struct containing source profile data for ions +Base.@kwdef struct ion_source_data + # struct containing source data for ions # is the source active or not active::Bool # An overall multiplier for the strength (i.e. 0.0 would turn off the source) @@ -454,10 +454,10 @@ Base.@kwdef struct ion_source_profile # in the event that the code is running with mk_int = Int32 but the rank is set to 0::Int64 end -Base.@kwdef struct electron_source_profile +Base.@kwdef struct electron_source_data # most of the electron parameters must be the same as for ions, so only # source strength (in the case of an ion energy source) and source Temperature - # can be different. The other four are set by the ion source profile. + # can be different. The other four are set by the ion source data. source_strength::mk_float source_T::mk_float active::Bool @@ -466,8 +466,8 @@ Base.@kwdef struct electron_source_profile source_type::String end -Base.@kwdef struct neutral_source_profile - # struct containing source profile data for neutrals +Base.@kwdef struct neutral_source_data + # struct containing source data for neutrals # is the source active or not active::Bool # An overall multiplier for the strength (i.e. 0.0 would turn off the source) From a3812405e996e51e0a3c20b0b7609c7d527ed0f9 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sat, 7 Sep 2024 10:37:05 +0100 Subject: [PATCH 05/23] Change create_moments_ion and create_moments_neutral functions to be compatible with source struct vectors --- moment_kinetics/src/velocity_moments.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl index 29684c26a..0b48035d7 100644 --- a/moment_kinetics/src/velocity_moments.jl +++ b/moment_kinetics/src/velocity_moments.jl @@ -138,7 +138,7 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar, entropy_production = allocate_shared_float(nz, nr, n_species) - if ion_source_settings.active + if any(x -> x.active, ion_source_settings) external_source_amplitude = allocate_shared_float(nz, nr) if evolve_density external_source_density_amplitude = allocate_shared_float(nz, nr) @@ -363,7 +363,7 @@ function create_moments_neutral(nz, nr, n_species, evolve_density, evolve_upar, dvth_dz = nothing end - if neutral_source_settings.active + if any(x -> x.active, neutral_source_settings) external_source_amplitude = allocate_shared_float(nz, nr) if evolve_density external_source_density_amplitude = allocate_shared_float(nz, nr) From fb0b082c44e057ddf26854eb93d8738b96d41801 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 09:19:27 +0100 Subject: [PATCH 06/23] Restructure main initialisation files and update loops to view each source for each species as a vector of structs, rather than a NamedTuple for each. This means each amplitude array for the moments structs has one extra dimension, one for each source, so that each source can be applied to each species pdf in turn. --- .../src/electron_kinetic_equation.jl | 11 +- moment_kinetics/src/external_sources.jl | 633 +++++++++++------- moment_kinetics/src/initial_conditions.jl | 4 +- .../src/moment_kinetics_structs.jl | 30 +- moment_kinetics/src/time_advance.jl | 14 +- moment_kinetics/src/velocity_moments.jl | 74 +- 6 files changed, 450 insertions(+), 316 deletions(-) diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl index fe1aeeb1e..e88476c29 100644 --- a/moment_kinetics/src/electron_kinetic_equation.jl +++ b/moment_kinetics/src/electron_kinetic_equation.jl @@ -1819,12 +1819,11 @@ function electron_kinetic_equation_euler_update!(fvec_out, fvec_in, moments, z, vperp, vpa, dt) end - if external_source_settings.electron.active - external_electron_source!(fvec_out.pdf_electron, fvec_in.pdf_electron, - moments.electron.dens, moments.electron.upar, moments, - composition, external_source_settings.electron, vperp, - vpa, dt) - end + total_external_electron_sources!(fvec_out.pdf_electron, fvec_in.pdf_electron, + moments.electron.dens, moments.electron.upar, moments, + composition, external_source_settings.electron, vperp, + vpa, dt) + if evolve_ppar electron_energy_equation!(fvec_out.electron_ppar, fvec_in.electron_ppar, diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index dcaca6905..2fc218aec 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -14,7 +14,9 @@ module external_sources export setup_external_sources!, external_ion_source!, external_neutral_source!, external_ion_source_controller!, external_neutral_source_controller!, initialize_external_source_amplitude!, - initialize_external_source_controller_integral! + initialize_external_source_controller_integral!, + total_external_ion_sources!, total_external_neutral_sources!, + total_external_ion_source_controllers!, total_external_neutral_source_controllers! using ..array_allocation: allocate_float, allocate_shared_float using ..calculus @@ -40,10 +42,10 @@ Returns a NamedTuple `(ion=ion_source_settings, neutral=neutral_source_settings) containing two NamedTuples of settings. """ function setup_external_sources!(input_dict, r, z, electron_physics) - function get_settings_ions() + function get_settings_ions(source_index, active_flag) input = set_defaults_and_check_section!( - input_dict, "ion_source"; - active=false, + input_dict, "ion_source_$source_index"; + active=active_flag, source_strength=1.0, source_n=1.0, source_T=1.0, @@ -150,17 +152,16 @@ function setup_external_sources!(input_dict, r, z, electron_physics) * "density_midpoint_control, energy, alphas, alphas-with-losses, " * "beam, beam-with-losses") end - println(input) - return ion_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, + return ion_source_data(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, z_amplitude=z_amplitude, PI_density_target=PI_density_target, PI_controller_amplitude, controller_source_profile, PI_density_target_ir, PI_density_target_iz, PI_density_target_rank) end - function get_settings_neutrals() + function get_settings_neutrals(source_index, active_flag) input = set_defaults_and_check_section!( - input_dict, "neutral_source"; - active=false, + input_dict, "neutral_source_$source_index"; + active=active_flag, source_strength=1.0, source_n=1.0, source_T=get(input_dict, "T_wall", 1.0), @@ -308,7 +309,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics) * "beam, beam-with-losses, recycling (for neutrals only)") end - return neutral_source_profile(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, + return neutral_source_data(; Dict(Symbol(k)=>v for (k,v) ∈ input)..., r_amplitude, z_amplitude=z_amplitude, PI_density_target=PI_density_target, PI_controller_amplitude, controller_source_profile, PI_density_target_ir, PI_density_target_iz, PI_density_target_rank) @@ -334,23 +335,51 @@ function setup_external_sources!(input_dict, r, z, electron_physics) end input["source_strength"] = ion_settings.source_strength end - return electron_source_profile(input["source_strength"], input["source_T"], + return electron_source_data(input["source_strength"], input["source_T"], ion_settings.active, ion_settings.r_amplitude, ion_settings.z_amplitude, ion_settings.source_type) end - ion_settings = get_settings_ions() + # put all ion sources into ion_source_data struct vector + ion_sources = ion_source_data[] + counter = 1 + while "ion_source_$counter" ∈ keys(input_dict) + push!(ion_sources, get_settings_ions(counter, true)) + counter += 1 + end + + # If there are no ion sources, add an inactive ion source to the vector + if counter == 1 + inactive_ion_source = get_settings_ions(1, false) + push!(ion_sources, inactive_ion_source) + counter += 1 + end + + # put all electron sources into electron_source_data struct vector, where + # each entry is a mirror of the ion source vector. + electron_sources = electron_source_data[] if electron_physics ∈ (braginskii_fluid, kinetic_electrons, kinetic_electrons_with_temperature_equation) - electron_settings = get_settings_electrons(ion_settings) + electron_sources .= get_settings_electrons.(ion_sources) else - electron_settings = electron_source_profile(ion_settings.source_strength, - ion_settings.source_T, false, ion_settings.r_amplitude, - ion_settings.z_amplitude, ion_settings.source_type) + # this is not very efficient because we're copying the same electron_settings + # for each ion source + push!(electron_sources, get_settings_electrons(inactive_ion_source)) end - neutral_settings = get_settings_neutrals() - - return (ion=ion_settings, electron=electron_settings, neutral=neutral_settings) + + # put all neutral sources into neutral_source_data struct vector + neutral_sources = neutral_source_data[] + counter = 1 + while "neutral_source_$counter" ∈ keys(input_dict) + push!(neutral_sources, get_settings_neutrals(counter, true)) + counter += 1 + end + # If there are no neutral sources, add an inactive neutral source to the vector + if counter == 1 + inactive_neutral_source = get_settings_neutrals(1, false) + push!(neutral_sources, inactive_neutral_source) + end + return (ion=ion_sources, electron=electron_sources, neutral=neutral_sources) end """ @@ -407,195 +436,204 @@ function initialize_external_source_amplitude!(moments, external_source_settings begin_r_z_region() ion_source_settings = external_source_settings.ion - if ion_source_settings.active - if ion_source_settings.source_type == "energy" - @loop_r_z ir iz begin - moments.ion.external_source_amplitude[iz,ir] = - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] - end - if moments.evolve_density + # The electron loop must be in the same as the ion loop so that each electron source + # can be matched to the corresponding ion source. + electron_source_settings = external_source_settings.electron + + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active + if ion_source_settings[index].source_type == "energy" @loop_r_z ir iz begin - moments.ion.external_source_density_amplitude[iz,ir] = 0.0 + moments.ion.external_source_amplitude[iz,ir,index] = + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] end - end - if moments.evolve_upar - @loop_r_z ir iz begin - moments.ion.external_source_momentum_amplitude[iz,ir] = - - moments.ion.dens[iz,ir] * moments.ion.upar[iz,ir] * - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] + if moments.evolve_density + @loop_r_z ir iz begin + moments.ion.external_source_density_amplitude[iz,ir,index] = 0.0 + end end - end - if moments.evolve_ppar + if moments.evolve_upar + @loop_r_z ir iz begin + moments.ion.external_source_momentum_amplitude[iz,ir,index]= + - moments.ion.dens[iz,ir] * moments.ion.upar[iz,ir] * + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end + end + if moments.evolve_ppar + @loop_r_z ir iz begin + moments.ion.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * ion_source.source_T + + moments.ion.upar[iz,ir]^2 - moments.ion.ppar[iz,ir]) * + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end + end + else @loop_r_z ir iz begin - moments.ion.external_source_pressure_amplitude[iz,ir] = - (0.5 * ion_source_settings.source_T + - moments.ion.upar[iz,ir]^2 - moments.ion.ppar[iz,ir]) * - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] + moments.ion.external_source_amplitude[iz,ir,index] = + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end + if moments.evolve_density + @loop_r_z ir iz begin + moments.ion.external_source_density_amplitude[iz,ir,index] = + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end + end + if moments.evolve_upar + @loop_r_z ir iz begin + moments.ion.external_source_momentum_amplitude[iz,ir,index] = 0.0 + end + end + if moments.evolve_ppar + @loop_r_z ir iz begin + moments.ion.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * ion_source_settings[index].source_T + + moments.ion.upar[iz,ir]^2) * + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end end end - else - @loop_r_z ir iz begin - moments.ion.external_source_amplitude[iz,ir] = - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] - end - if moments.evolve_density + end + + # now do same for electron sources, which (if present) are mostly mirrors of ion sources + println(electron_source_settings) + if electron_source_settings[index].active + if electron_source_settings[index].source_type == "energy" @loop_r_z ir iz begin - moments.ion.external_source_density_amplitude[iz,ir] = - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] + moments.electron.external_source_amplitude[iz,ir,index] = + electron_source_settings[index].source_strength * + electron_source_settings[index].r_amplitude[ir] * + electron_source_settings[index].z_amplitude[iz] end - end - if moments.evolve_upar @loop_r_z ir iz begin - moments.ion.external_source_momentum_amplitude[iz,ir] = 0.0 + moments.electron.external_source_density_amplitude[iz,ir,index] = 0.0 end - end - if moments.evolve_ppar @loop_r_z ir iz begin - moments.ion.external_source_pressure_amplitude[iz,ir] = - (0.5 * ion_source_settings.source_T + - moments.ion.upar[iz,ir]^2) * - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] + moments.electron.external_source_momentum_amplitude[iz,ir,index] = + - moments.electron.dens[iz,ir] * moments.electron.upar[iz,ir] * + electron_source_settings[index].source_strength * + electron_source_settings[index].r_amplitude[ir] * + electron_source_settings[index].z_amplitude[iz] end - end - end - end - - electron_source_settings = external_source_settings.electron - if electron_source_settings.active - if electron_source_settings.source_type == "energy" - @loop_r_z ir iz begin - moments.electron.external_source_amplitude[iz,ir] = - electron_source_settings.source_strength * - electron_source_settings.r_amplitude[ir] * - electron_source_settings.z_amplitude[iz] - end - @loop_r_z ir iz begin - moments.electron.external_source_density_amplitude[iz,ir] = 0.0 - end - @loop_r_z ir iz begin - moments.electron.external_source_momentum_amplitude[iz,ir] = - - moments.electron.dens[iz,ir] * moments.electron.upar[iz,ir] * - electron_source_settings.source_strength * - electron_source_settings.r_amplitude[ir] * - electron_source_settings.z_amplitude[iz] - end - @loop_r_z ir iz begin - moments.electron.external_source_pressure_amplitude[iz,ir] = - (0.5 * electron_source_settings.source_T + - moments.electron.upar[iz,ir]^2 - moments.electron.ppar[iz,ir]) * - electron_source_settings.source_strength * - electron_source_settings.r_amplitude[ir] * - electron_source_settings.z_amplitude[iz] - end - else - @loop_r_z ir iz begin - moments.electron.external_source_amplitude[iz,ir] = - moments.ion.external_source_amplitude[iz,ir] - end - if moments.evolve_density @loop_r_z ir iz begin - moments.electron.external_source_density_amplitude[iz,ir] = - moments.ion.external_source_density_amplitude[iz,ir] + moments.electron.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * electron_source_settings[index].source_T + + moments.electron.upar[iz,ir]^2 - moments.electron.ppar[iz,ir]) * + electron_source_settings[index].source_strength * + electron_source_settings[index].r_amplitude[ir] * + electron_source_settings[index].z_amplitude[iz] end else @loop_r_z ir iz begin - # Note set this using *ion* settings to force electron density source - # to always be equal to ion density source (even when - # evolve_density=false) to ensure the source does not inject charge - # into the simulation. - moments.electron.external_source_density_amplitude[iz,ir] = - ion_source_settings.source_strength * - ion_source_settings.r_amplitude[ir] * - ion_source_settings.z_amplitude[iz] + moments.electron.external_source_amplitude[iz,ir,index] = + moments.ion.external_source_amplitude[iz,ir,index] + end + if moments.evolve_density + @loop_r_z ir iz begin + moments.electron.external_source_density_amplitude[iz,ir,index] = + moments.ion.external_source_density_amplitude[iz,ir,index] + end + else + @loop_r_z ir iz begin + # Note set this using *ion* settings to force electron density source + # to always be equal to ion density source (even when + # evolve_density=false) to ensure the source does not inject charge + # into the simulation. + moments.electron.external_source_density_amplitude[iz,ir,index] = + ion_source_settings[index].source_strength * + ion_source_settings[index].r_amplitude[ir] * + ion_source_settings[index].z_amplitude[iz] + end + end + @loop_r_z ir iz begin + moments.electron.external_source_momentum_amplitude[iz,ir,index] = 0.0 + end + @loop_r_z ir iz begin + moments.electron.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * electron_source_settings[index].source_T + + moments.electron.upar[iz,ir]^2) * + moments.electron.external_source_amplitude[iz,ir,index] end - end - @loop_r_z ir iz begin - moments.electron.external_source_momentum_amplitude[iz,ir] = 0.0 - end - @loop_r_z ir iz begin - moments.electron.external_source_pressure_amplitude[iz,ir] = - (0.5 * electron_source_settings.source_T + - moments.electron.upar[iz,ir]^2) * - moments.electron.external_source_amplitude[iz,ir] end end end if n_neutral_species > 0 neutral_source_settings = external_source_settings.neutral - if neutral_source_settings.active - if neutral_source_settings.source_type == "energy" - @loop_r_z ir iz begin - moments.neutral.external_source_amplitude[iz,ir] = - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] - end - if moments.evolve_density + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + if neutral_source.source_type == "energy" @loop_r_z ir iz begin - moments.neutral.external_source_density_amplitude[iz,ir] = 0.0 + moments.neutral.external_source_amplitude[iz,ir,index] = + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] end - end - if moments.evolve_upar - @loop_r_z ir iz begin - moments.neutral.external_source_momentum_amplitude[iz,ir] = - - moments.neutral.dens[iz,ir] * moments.neutral.upar[iz,ir] * - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] + if moments.evolve_density + @loop_r_z ir iz begin + moments.neutral.external_source_density_amplitude[iz,ir,index] = 0.0 + end end - end - if moments.evolve_ppar - @loop_r_z ir iz begin - moments.neutral.external_source_pressure_amplitude[iz,ir] = - (0.5 * neutral_source_settings.source_T + - moments.neutral.upar[iz,ir]^2 - - moments.neutral.ppar[iz,ir]) * - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] + if moments.evolve_upar + @loop_r_z ir iz begin + moments.neutral.external_source_momentum_amplitude[iz,ir,index] = + - moments.neutral.dens[iz,ir] * moments.neutral.upar[iz,ir] * + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] + end end - end - else - @loop_r_z ir iz begin - moments.neutral.external_source_amplitude[iz,ir] = - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] + if moments.evolve_ppar + @loop_r_z ir iz begin + moments.neutral.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * neutral_source_settings[index].source_T + + moments.neutral.upar[iz,ir]^2 - + moments.neutral.ppar[iz,ir]) * + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] + end end - if moments.evolve_density + else @loop_r_z ir iz begin - moments.neutral.external_source_density_amplitude[iz,ir] = - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] + moments.neutral.external_source_amplitude[iz,ir,index] = + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] end - end - if moments.evolve_upar - @loop_r_z ir iz begin - moments.neutral.external_source_momentum_amplitude[iz,ir] = 0.0 + if moments.evolve_density + @loop_r_z ir iz begin + moments.neutral.external_source_density_amplitude[iz,ir,index] = + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] + end end - end - if moments.evolve_ppar - @loop_r_z ir iz begin - moments.neutral.external_source_pressure_amplitude[iz,ir] = - (0.5 * neutral_source_settings.source_T + - moments.neutral.uz[iz,ir]^2) * - neutral_source_settings.source_strength * - neutral_source_settings.r_amplitude[ir] * - neutral_source_settings.z_amplitude[iz] + if moments.evolve_upar + @loop_r_z ir iz begin + moments.neutral.external_source_momentum_amplitude[iz,ir,index] = 0.0 + end + end + if moments.evolve_ppar + @loop_r_z ir iz begin + moments.neutral.external_source_pressure_amplitude[iz,ir,index] = + (0.5 * neutral_source_settings[index].source_T + + moments.neutral.uz[iz,ir]^2) * + neutral_source_settings[index].source_strength * + neutral_source_settings[index].r_amplitude[ir] * + neutral_source_settings[index].z_amplitude[iz] + end end end end @@ -618,21 +656,21 @@ function initialize_external_source_controller_integral!( begin_serial_region() @serial_region begin ion_source_settings = external_source_settings.ion - if ion_source_settings.active - if ion_source_settings.PI_density_controller_I != 0.0 && - ion_source_settings.source_type ∈ ("density_profile_control", - "density_midpoint_control") - moments.ion.external_source_controller_integral .= 0.0 + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active && + ion_source_settings[index].PI_density_controller_I != 0.0 && + ion_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control") + moments.ion.external_source_controller_integral[:, :, index] .= 0.0 end end if n_neutral_species > 0 neutral_source_settings = external_source_settings.neutral - if neutral_source_settings.active - if neutral_source_settings.PI_density_controller_I != 0.0 && - neutral_source_settings.source_type ∈ ("density_profile_control", - "density_midpoint_control") - moments.neutral.external_source_controller_integral .= 0.0 + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active && + neutral_source_settings[index].PI_density_controller_I != 0.0 && + neutral_source_settings[index].source_type ∈ ("density_profile_control", "density_midpoint_control") + moments.neutral.external_source_controller_integral[:, :, index] .= 0.0 end end end @@ -641,17 +679,33 @@ function initialize_external_source_controller_integral!( return nothing end +""" + total_external_ion_sources!(pdf, fvec, moments, ion_sources, vperp, vpa, dt, scratch_dummy) + +Contribute all of the ion sources to the ion pdf, one by one. +""" +function total_external_ion_sources!(pdf, fvec, moments, ion_sources, vperp, + vpa, dt, scratch_dummy) + for index ∈ eachindex(ion_sources) + if ion_sources[index].active + external_ion_source!(pdf, fvec, moments, ion_sources[index], + index, vperp, vpa, dt, scratch_dummy) + end + end + return nothing +end + """ external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vpa, dt) Add external source term to the ion kinetic equation. """ -function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vpa, dt, scratch_dummy) +function external_ion_source!(pdf, fvec, moments, ion_source, index, vperp, vpa, dt, scratch_dummy) - source_type = ion_source_settings.source_type - source_amplitude = moments.ion.external_source_amplitude - source_T = ion_source_settings.source_T - source_n = ion_source_settings.source_n + source_type = ion_source.source_type + @views source_amplitude = moments.ion.external_source_amplitude[:,:,index] + source_T = ion_source.source_T + source_n = ion_source.source_n if vperp.n == 1 vth_factor = 1.0 / sqrt(source_T) else @@ -862,6 +916,24 @@ function external_ion_source!(pdf, fvec, moments, ion_source_settings, vperp, vp return nothing end +""" + total_external_electron_sources!(pdf, fvec, moments, electron_sources, vperp, vpa, dt, scratch_dummy) + +Contribute all of the electron sources to the electron pdf, one by one. +""" +function total_external_electron_sources!(pdf_out, pdf_in, electron_density, electron_upar, + moments, composition, electron_sources, vperp, + vpa, dt) + for index ∈ eachindex(electron_sources) + if electron_sources[index].active + external_electron_source!(pdf_out, pdf_in, electron_density, electron_upar, + moments, composition, electron_sources[index], index, + vperp, vpa, dt) + end + end + return nothing +end + """ external_electron_source!(pdf, fvec, moments, electron_source_settings, vperp, vpa, dt) @@ -869,14 +941,14 @@ end Add external source term to the electron kinetic equation. """ function external_electron_source!(pdf_out, pdf_in, electron_density, electron_upar, - moments, composition, electron_source_settings, vperp, + moments, composition, electron_source, index, vperp, vpa, dt) begin_r_z_vperp_region() me_over_mi = composition.me_over_mi - source_amplitude = moments.electron.external_source_amplitude - source_T = electron_source_settings.source_T + @views source_amplitude = moments.electron.external_source_amplitude[:, :, index] + source_T = electron_source.source_T if vperp.n == 1 vth_factor = 1.0 / sqrt(source_T / me_over_mi) else @@ -913,18 +985,36 @@ function external_electron_source!(pdf_out, pdf_in, electron_density, electron_u return nothing end + +""" + total_external_neutral_sources!(pdf, fvec, moments, neutral_sources, vperp, vpa, dt, scratch_dummy) + +Contribute all of the neutral sources to the neutral pdf, one by one. +""" +function total_external_neutral_sources!(pdf, fvec, moments, neutral_sources, + vzeta, vr, vz, dt) + for index ∈ eachindex(neutral_sources) + if neutral_sources[index].active + external_neutral_source!(pdf, fvec, moments, neutral_sources[index], + index, vzeta, vr, vz, dt) + end + end + return nothing +end + + """ external_neutral_source!(pdf, fvec, moments, neutral_source_settings, vzeta, vr, vz, dt) Add external source term to the neutral kinetic equation. """ -function external_neutral_source!(pdf, fvec, moments, neutral_source_settings, vzeta, vr, - vz, dt) +function external_neutral_source!(pdf, fvec, moments, neutral_source, index, vzeta, vr, + vz, dt) begin_sn_r_z_vzeta_vr_region() - source_amplitude = moments.neutral.external_source_amplitude - source_T = neutral_source_settings.source_T + @views source_amplitude = moments.neutral.external_source_amplitude[:, :, index] + source_T = neutral_source.source_T if vzeta.n == 1 && vr.n == 1 vth_factor = 1.0 / sqrt(source_T) else @@ -1009,13 +1099,27 @@ function external_neutral_source!(pdf, fvec, moments, neutral_source_settings, v return nothing end +""" + total_external_ion_source_controllers!(fvec_in, moments, ion_sources, dt) + +Contribute all of the ion source controllers to fvec_in, one by one. +""" +function total_external_ion_source_controllers!(fvec_in, moments, ion_sources, dt) + for index ∈ eachindex(ion_sources) + if ion_sources[index].active + external_ion_source_controller!(fvec_in, moments, ion_sources[index], index, dt) + end + end + return nothing +end + """ external_ion_source_controller!(fvec_in, moments, ion_source_settings, dt) Calculate the amplitude when using a PI controller for the density to set the external source amplitude. """ -function external_ion_source_controller!(fvec_in, moments, ion_source_settings, dt) +function external_ion_source_controller!(fvec_in, moments, ion_source_settings, index, dt) begin_r_z_region() is = 1 @@ -1027,15 +1131,15 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, if ion_source_settings.source_type == "Maxwellian" if moments.evolve_ppar @loop_r_z ir iz begin - ion_moments.external_source_pressure_amplitude[iz,ir] = + ion_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * ion_source_settings.source_T + upar[iz,ir,is]^2) * - ion_moments.external_source_amplitude[iz,ir] + ion_moments.external_source_amplitude[iz,ir,index] end end elseif ion_source_settings.source_type == "energy" if moments.evolve_upar @loop_r_z ir iz begin - ion_moments.external_source_momentum_amplitude[iz,ir] = + ion_moments.external_source_momentum_amplitude[iz,ir,index] = - density[iz,ir] * upar[iz,ir] * ion_source_settings.source_strength * ion_source_settings.r_amplitude[ir] * @@ -1044,7 +1148,7 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, end if moments.evolve_ppar @loop_r_z ir iz begin - ion_moments.external_source_pressure_amplitude[iz,ir] = + ion_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * ion_source_settings.source_T + upar[iz,ir]^2 - ppar[iz,ir]) * ion_source_settings.source_strength * ion_source_settings.r_amplitude[ir] * @@ -1065,13 +1169,13 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, ion_source_settings.PI_density_target_ir, is] n_error = ion_source_settings.PI_density_target - n_mid - ion_moments.external_source_controller_integral[1,1] += + ion_moments.external_source_controller_integral[1,1,index] += dt * ion_source_settings.PI_density_controller_I * n_error # Only want a source, so never allow amplitude to be negative amplitude = max( ion_source_settings.PI_density_controller_P * n_error + - ion_moments.external_source_controller_integral[1,1], + ion_moments.external_source_controller_integral[1,1,index], 0) else amplitude = nothing @@ -1085,18 +1189,18 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, amplitude = controller_amplitude[1] @loop_r_z ir iz begin - ion_moments.external_source_amplitude[iz,ir] = + ion_moments.external_source_amplitude[iz,ir,index] = amplitude * ion_source_settings.controller_source_profile[iz,ir] end if moments.evolve_density @loop_r_z ir iz begin - ion_moments.external_source_density_amplitude[iz,ir] = + ion_moments.external_source_density_amplitude[iz,ir,index] = amplitude * ion_source_settings.controller_source_profile[iz,ir] end end if moments.evolve_ppar @loop_r_z ir iz begin - ion_moments.external_source_pressure_amplitude[iz,ir] = + ion_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * ion_source_settings.source_T + upar[iz,ir,is]^2) * amplitude * ion_source_settings.controller_source_profile[iz,ir] end @@ -1111,20 +1215,20 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, amplitude = ion_moments.external_source_amplitude @loop_r_z ir iz begin n_error = target[iz,ir] - density[iz,ir,is] - integral[iz,ir] += dt * I * n_error + integral[iz,ir,index] += dt * I * n_error # Only want a source, so never allow amplitude to be negative - amplitude[iz,ir] = max(P * n_error + integral[iz,ir], 0) + amplitude[iz,ir,index] = max(P * n_error + integral[iz,ir,index], 0) end if moments.evolve_density @loop_r_z ir iz begin - ion_moments.external_source_density_amplitude[iz,ir] = amplitude[iz,ir] + ion_moments.external_source_density_amplitude[iz,ir,index] = amplitude[iz,ir,index] end end if moments.evolve_ppar @loop_r_z ir iz begin - ion_moments.external_source_pressure_amplitude[iz,ir] = + ion_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * ion_source_settings.source_T + upar[iz,ir,is]^2) * - amplitude[iz,ir] + amplitude[iz,ir,index] end end elseif ion_source_settings.source_type == "alphas" @@ -1142,6 +1246,20 @@ function external_ion_source_controller!(fvec_in, moments, ion_source_settings, return nothing end +""" + total_external_electron_source_controllers!(fvec_in, moments, electron_sources, dt) + +Contribute all of the electron source controllers to fvec_in, one by one. +""" +function total_external_electron_source_controllers!(fvec_in, moments, electron_sources, dt) + for index ∈ eachindex(electron_sources) + if electron_sources[index].active + external_electron_source_controller!(fvec_in, moments, electron_sources[index], index, dt) + end + end + return nothing +end + """ external_electron_source_controller!(fvec_in, moments, electron_source_settings, dt) @@ -1155,30 +1273,30 @@ called before this function is called so that `moments.ion.external_source_ampli up to date. """ function external_electron_source_controller!(fvec_in, moments, electron_source_settings, - dt) + index, dt) begin_r_z_region() is = 1 electron_moments = moments.electron - ion_source_amplitude = moments.ion.external_source_amplitude + @views ion_source_amplitude = moments.ion.external_source_amplitude[:, :, index] if electron_source_settings.source_type == "Maxwellian" @loop_r_z ir iz begin - electron_moments.external_source_pressure_amplitude[iz,ir] = + electron_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * electron_source_settings.source_T + fvec_in.electron_upar[iz,ir,is]^2) * - electron_moments.external_source_amplitude[iz,ir] + electron_moments.external_source_amplitude[iz,ir,index] end elseif electron_source_settings.source_type == "energy" @loop_r_z ir iz begin - electron_moments.external_source_momentum_amplitude[iz,ir] = + electron_moments.external_source_momentum_amplitude[iz,ir,index] = - electron_moments.density[iz,ir] * electron_moments.upar[iz,ir] * electron_source_settings.source_strength * electron_source_settings.r_amplitude[ir] * electron_source_settings.z_amplitude[iz] end @loop_r_z ir iz begin - electron_moments.external_source_pressure_amplitude[iz,ir] = + electron_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * electron_source_settings.source_T + electron_moments.upar[iz,ir]^2 - electron_moments.ppar[iz,ir]) * electron_source_settings.source_strength * @@ -1187,30 +1305,45 @@ function external_electron_source_controller!(fvec_in, moments, electron_source_ end else @loop_r_z ir iz begin - electron_moments.external_source_amplitude[iz,ir] = ion_source_amplitude[iz,ir] + electron_moments.external_source_amplitude[iz,ir,index] = ion_source_amplitude[iz,ir,index] end @loop_r_z ir iz begin - electron_moments.external_source_momentum_amplitude[iz,ir] = + electron_moments.external_source_momentum_amplitude[iz,ir,index] = - electron_moments.density[iz,ir] * electron_moments.upar[iz,ir] * - electron_moments.external_source_amplitude[iz,ir] + electron_moments.external_source_amplitude[iz,ir,index] end @loop_r_z ir iz begin - electron_moments.external_source_pressure_amplitude[iz,ir] = + electron_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * electron_source_settings.source_T + electron_moments.upar[iz,ir]^2 - electron_moments.ppar[iz,ir]) * - electron_moments.external_source_amplitude[iz,ir] + electron_moments.external_source_amplitude[iz,ir,index] end end # Density source is always the same as the ion one @loop_r_z ir iz begin - electron_moments.external_source_density_amplitude[iz,ir] = - moments.ion.external_source_density_amplitude[iz,ir] + electron_moments.external_source_density_amplitude[iz,ir,index] = + moments.ion.external_source_density_amplitude[iz,ir,index] end return nothing end +""" + total_external_neutral_source_controllers!(fvec_in, moments, neutral_sources, dt) + +Contribute all of the neutral source controllers to fvec_in, one by one. +""" +function total_external_neutral_source_controllers!(fvec_in, moments, neutral_sources, r, z, dt) + for index ∈ eachindex(neutral_sources) + if neutral_sources[index].active + external_neutral_source_controller!(fvec_in, moments, neutral_sources[index], + index, r, z, dt) + end + end + return nothing +end + """ external_neutral_source_controller!(fvec_in, moments, neutral_source_settings, r, z, dt) @@ -1218,8 +1351,8 @@ end Calculate the amplitude when using a PI controller for the density to set the external source amplitude. """ -function external_neutral_source_controller!(fvec_in, moments, neutral_source_settings, r, - z, dt) +function external_neutral_source_controller!(fvec_in, moments, neutral_source_settings, index, + r, z, dt) begin_r_z_region() is = 1 @@ -1228,15 +1361,15 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se if neutral_source_settings.source_type == "Maxwellian" if moments.evolve_ppar @loop_r_z ir iz begin - neutral_moments.external_source_pressure_amplitude[iz,ir] = + neutral_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * neutral_source_settings.source_T + fvec_in.upar[iz,ir,is]^2) * - neutral_moments.external_source_amplitude[iz,ir] + neutral_moments.external_source_amplitude[iz,ir,index] end end elseif neutral_source_settings.source_type == "energy" if moments.evolve_upar @loop_r_z ir iz begin - neutral_moments.external_source_momentum_amplitude[iz,ir] = + neutral_moments.external_source_momentum_amplitude[iz,ir,index] = - neutral_moments.density[iz,ir] * neutral_moments.uz[iz,ir] * neutral_source_settings.source_strength * neutral_source_settings.r_amplitude[ir] * @@ -1245,7 +1378,7 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se end if moments.evolve_ppar @loop_r_z ir iz begin - neutral_moments.external_source_pressure_amplitude[iz,ir] = + neutral_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * neutral_source_settings.source_T + neutral_moments.uz[iz,ir]^2 - neutral_moments.pz[iz,ir]) * neutral_source_settings.source_strength * @@ -1268,13 +1401,13 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se is] n_error = neutral_source_settings.PI_density_target - n_mid - neutral_moments.external_source_controller_integral[1,1] += + neutral_moments.external_source_controller_integral[1,1,index] += dt * neutral_source_settings.PI_density_controller_I * n_error # Only want a source, so never allow amplitude to be negative amplitude = max( neutral_source_settings.PI_density_controller_P * n_error + - neutral_moments.external_source_controller_integral[1,1], + neutral_moments.external_source_controller_integral[1,1,index], 0) else amplitude = nothing @@ -1288,20 +1421,20 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se amplitude = controller_amplitude[1] @loop_r_z ir iz begin - neutral_moments.external_source_amplitude[iz,ir] = - amplitude * neutral_source_settings.controller_source_profile[iz,ir] + neutral_moments.external_source_amplitude[iz,ir,index] = + amplitude * neutral_source_settings.controller_source_profile[iz,ir,index] end if moments.evolve_density @loop_r_z ir iz begin - neutral_moments.external_source_density_amplitude[iz,ir] = - amplitude * neutral_source_settings.controller_source_profile[iz,ir] + neutral_moments.external_source_density_amplitude[iz,ir,index] = + amplitude * neutral_source_settings.controller_source_profile[iz,ir,index] end end if moments.evolve_ppar @loop_r_z ir iz begin - neutral_moments.external_source_pressure_amplitude[iz,ir] = + neutral_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * neutral_source_settings.source_T + fvec_in.upar[iz,ir,is]^2) * - amplitude * neutral_source_settings.controller_source_profile[iz,ir] + amplitude * neutral_source_settings.controller_source_profile[iz,ir,index] end end elseif neutral_source_settings.source_type == "density_profile_control" @@ -1315,19 +1448,19 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se amplitude = neutral_moments.external_source_amplitude @loop_r_z ir iz begin n_error = target[iz,ir] - density[iz,ir,is] - PI_integral[iz,ir] += dt * I * n_error - amplitude[iz,ir] = P * n_error + PI_integral[iz,ir] + PI_integral[iz,ir,index] += dt * I * n_error + amplitude[iz,ir,index] = P * n_error + PI_integral[iz,ir,index] end if moments.evolve_density @loop_r_z ir iz begin - neutral_moments.external_source_density_amplitude[iz,ir] = amplitude[iz,ir] + neutral_moments.external_source_density_amplitude[iz,ir,index] = amplitude[iz,ir,index] end end if moments.evolve_ppar @loop_r_z ir iz begin - neutral_moments.external_source_pressure_amplitude[iz,ir] = + neutral_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * neutral_source_settings.source_T + fvec_in.upar[iz,ir,is]^2) * - amplitude[iz,ir] + amplitude[iz,ir,index] end end elseif neutral_source_settings.source_type == "recycling" @@ -1365,18 +1498,18 @@ function external_neutral_source_controller!(fvec_in, moments, neutral_source_se profile = neutral_source_settings.controller_source_profile prefactor = target_flux * neutral_source_settings.recycling_controller_fraction @loop_r_z ir iz begin - amplitude[iz,ir] = prefactor * profile[iz,ir] + amplitude[iz,ir,index] = prefactor * profile[iz,ir] end if moments.evolve_density @loop_r_z ir iz begin - neutral_moments.external_source_density_amplitude[iz,ir] = amplitude[iz,ir] + neutral_moments.external_source_density_amplitude[iz,ir,index] = amplitude[iz,ir,index] end end if moments.evolve_ppar @loop_r_z ir iz begin - neutral_moments.external_source_pressure_amplitude[iz,ir] = + neutral_moments.external_source_pressure_amplitude[iz,ir,index] = (0.5 * neutral_source_settings.source_T + fvec_in.upar[iz,ir,is]^2) * - amplitude[iz,ir] + amplitude[iz,ir,index] end end else diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 8d7db7598..2549e07f4 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -67,8 +67,8 @@ function allocate_pdf_and_moments(composition, r, z, vperp, vpa, vzeta, vr, vz, evolve_moments.density, evolve_moments.parallel_flow, evolve_moments.parallel_pressure, external_source_settings.ion, num_diss_params) - electron = create_moments_electron(z.n, r.n, - composition.electron_physics, num_diss_params) + electron = create_moments_electron(z.n, r.n, composition.electron_physics, + num_diss_params, length(external_source_settings.electron)) neutral = create_moments_neutral(z.n, r.n, composition.n_neutral_species, evolve_moments.density, evolve_moments.parallel_flow, evolve_moments.parallel_pressure, external_source_settings.neutral, diff --git a/moment_kinetics/src/moment_kinetics_structs.jl b/moment_kinetics/src/moment_kinetics_structs.jl index acb1a378e..d16138fd1 100644 --- a/moment_kinetics/src/moment_kinetics_structs.jl +++ b/moment_kinetics/src/moment_kinetics_structs.jl @@ -132,18 +132,18 @@ struct moments_ion_substruct dvth_dz::Union{MPISharedArray{mk_float,3},Nothing} # this is the entropy production dS/dt = - int (ln f sum_s' C_ss' [f_s,f_s']) d^3 v dSdt::MPISharedArray{mk_float,3} - # Spatially varying amplitude of the external source term - external_source_amplitude::MPISharedArray{mk_float,2} + # Spatially varying amplitude of the external source term (third index is for different sources) + external_source_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the density moment of the external source term - external_source_density_amplitude::MPISharedArray{mk_float,2} + external_source_density_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel momentum moment of the external source # term - external_source_momentum_amplitude::MPISharedArray{mk_float,2} + external_source_momentum_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel pressure moment of the external source # term - external_source_pressure_amplitude::MPISharedArray{mk_float,2} + external_source_pressure_amplitude::MPISharedArray{mk_float,3} # Integral term for the PID controller of the external source term - external_source_controller_integral::MPISharedArray{mk_float,2} + external_source_controller_integral::MPISharedArray{mk_float,3} # Store coefficient 'A' from applying moment constraints so we can write it out as a # diagnostic constraints_A_coefficient::Union{MPISharedArray{mk_float,3},Nothing} @@ -184,15 +184,15 @@ struct moments_electron_substruct # this is the parallel friction force between ions and electrons parallel_friction::MPISharedArray{mk_float,2} # Spatially varying amplitude of the external source term - external_source_amplitude::MPISharedArray{mk_float,2} + external_source_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the density moment of the external source term - external_source_density_amplitude::MPISharedArray{mk_float,2} + external_source_density_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel momentum moment of the external source # term - external_source_momentum_amplitude::MPISharedArray{mk_float,2} + external_source_momentum_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel pressure moment of the external source # term - external_source_pressure_amplitude::MPISharedArray{mk_float,2} + external_source_pressure_amplitude::MPISharedArray{mk_float,3} # if evolve_ppar = true, then the velocity variable is (vpa - upa)/vth, which introduces # a factor of vth for each power of wpa in velocity space integrals. # v_norm_fac accounts for this: it is vth if using the above definition for the parallel velocity, @@ -297,17 +297,17 @@ struct moments_neutral_substruct # this is the z-derivative of the heat flux along z dqz_dz::Union{MPISharedArray{mk_float,3},Nothing} # Spatially varying amplitude of the external source term - external_source_amplitude::MPISharedArray{mk_float,2} + external_source_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the density moment of the external source term - external_source_density_amplitude::MPISharedArray{mk_float,2} + external_source_density_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel momentum moment of the external source # term - external_source_momentum_amplitude::MPISharedArray{mk_float,2} + external_source_momentum_amplitude::MPISharedArray{mk_float,3} # Spatially varying amplitude of the parallel pressure moment of the external source # term - external_source_pressure_amplitude::MPISharedArray{mk_float,2} + external_source_pressure_amplitude::MPISharedArray{mk_float,3} # Integral term for the PID controller of the external source term - external_source_controller_integral::MPISharedArray{mk_float,2} + external_source_controller_integral::MPISharedArray{mk_float,3} # Store coefficient 'A' from applying moment constraints so we can write it out as a # diagnostic constraints_A_coefficient::Union{MPISharedArray{mk_float,3},Nothing} diff --git a/moment_kinetics/src/time_advance.jl b/moment_kinetics/src/time_advance.jl index 774028179..40dcca66d 100644 --- a/moment_kinetics/src/time_advance.jl +++ b/moment_kinetics/src/time_advance.jl @@ -1197,8 +1197,8 @@ function setup_advance_flags(moments, composition, t_params, collisions, if collisions.mxwl_diff.D_nn > 0.0 advance_maxwell_diffusion_nn = true end - advance_external_source = external_source_settings.ion.active && !t_params.implicit_ion_advance - advance_neutral_external_source = external_source_settings.neutral.active + advance_external_source = any(x -> x.active, external_source_settings.ion) && !t_params.implicit_ion_advance + advance_neutral_external_source = any(x -> x.active, external_source_settings.neutral) advance_ion_numerical_dissipation = !(t_params.implicit_ion_advance || t_params.implicit_vpa_advection) advance_neutral_numerical_dissipation = true # if evolving the density, must advance the continuity equation, @@ -1369,7 +1369,7 @@ function setup_implicit_advance_flags(moments, composition, t_params, collisions end end advance_krook_collisions_ii = collisions.krook.nuii0 > 0.0 - advance_external_source = external_source_settings.ion.active + advance_external_source = any(x -> x.active, external_source_settings.ion) advance_ion_numerical_dissipation = true advance_sources = moments.evolve_density || moments.evolve_upar || moments.evolve_ppar explicit_weakform_fp_collisions = collisions.fkpl.nuii > 0.0 && vperp.n > 1 @@ -3222,11 +3222,11 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, neutral_z_advect, neutral_r_advect, neutral_vz_advect = advect_objects.neutral_z_advect, advect_objects.neutral_r_advect, advect_objects.neutral_vz_advect if advance.external_source - external_ion_source_controller!(fvec_in, moments, external_source_settings.ion, + total_external_ion_source_controllers!(fvec_in, moments, external_source_settings.ion, dt) end if advance.neutral_external_source - external_neutral_source_controller!(fvec_in, moments, + total_external_neutral_source_controllers!(fvec_in, moments, external_source_settings.neutral, r, z, dt) end @@ -3350,11 +3350,11 @@ function euler_time_advance!(fvec_out, fvec_in, pdf, fields, moments, end if advance.external_source - external_ion_source!(fvec_out.pdf, fvec_in, moments, external_source_settings.ion, + total_external_ion_sources!(fvec_out.pdf, fvec_in, moments, external_source_settings.ion, vperp, vpa, dt, scratch_dummy) end if advance.neutral_external_source - external_neutral_source!(fvec_out.pdf_neutral, fvec_in, moments, + total_external_neutral_sources!(fvec_out.pdf_neutral, fvec_in, moments, external_source_settings.neutral, vzeta, vr, vz, dt) end diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl index 0b48035d7..b362b721b 100644 --- a/moment_kinetics/src/velocity_moments.jl +++ b/moment_kinetics/src/velocity_moments.jl @@ -138,39 +138,40 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar, entropy_production = allocate_shared_float(nz, nr, n_species) + n_sources = length(ion_source_settings) if any(x -> x.active, ion_source_settings) - external_source_amplitude = allocate_shared_float(nz, nr) + external_source_amplitude = allocate_shared_float(nz, nr, n_sources) if evolve_density - external_source_density_amplitude = allocate_shared_float(nz, nr) + external_source_density_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_density_amplitude = allocate_shared_float(1, 1) + external_source_density_amplitude = allocate_shared_float(1, 1, n_sources) end if evolve_upar - external_source_momentum_amplitude = allocate_shared_float(nz, nr) + external_source_momentum_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_momentum_amplitude = allocate_shared_float(1, 1) + external_source_momentum_amplitude = allocate_shared_float(1, 1, n_sources) end if evolve_ppar - external_source_pressure_amplitude = allocate_shared_float(nz, nr) + external_source_pressure_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_pressure_amplitude = allocate_shared_float(1, 1) + external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) end if ion_source_settings.PI_density_controller_I != 0.0 && ion_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") if ion_source_settings.source_type == "density_profile_control" - external_source_controller_integral = allocate_shared_float(nz, nr) + external_source_controller_integral = allocate_shared_float(nz, nr, n_sources) else - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end else - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end else - external_source_amplitude = allocate_shared_float(1, 1) - external_source_density_amplitude = allocate_shared_float(1, 1) - external_source_momentum_amplitude = allocate_shared_float(1, 1) - external_source_pressure_amplitude = allocate_shared_float(1, 1) - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_density_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_momentum_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end if evolve_density || evolve_upar || evolve_ppar @@ -199,7 +200,7 @@ end """ create a moment struct containing information about the electron moments """ -function create_moments_electron(nz, nr, electron_model, num_diss_params) +function create_moments_electron(nz, nr, electron_model, num_diss_params, n_sources) # allocate array used for the particle density density = allocate_shared_float(nz, nr) # initialise Bool variable that indicates if the density is updated for each species @@ -222,11 +223,11 @@ function create_moments_electron(nz, nr, electron_model, num_diss_params) parallel_heat_flux_updated = Ref(false) # allocate array used for the election-ion parallel friction force parallel_friction_force = allocate_shared_float(nz, nr) - # allocate arrays used for external sources - external_source_amplitude = allocate_shared_float(nz, nr) - external_source_density_amplitude = allocate_shared_float(nz, nr) - external_source_momentum_amplitude = allocate_shared_float(nz, nr) - external_source_pressure_amplitude = allocate_shared_float(nz, nr) + # allocate arrays used for external sources (third index is for the different sources) + external_source_amplitude = allocate_shared_float(nz, nr, n_sources) + external_source_density_amplitude = allocate_shared_float(nz, nr, n_sources) + external_source_momentum_amplitude = allocate_shared_float(nz, nr, n_sources) + external_source_pressure_amplitude = allocate_shared_float(nz, nr, n_sources) # allocate array used for the thermal speed thermal_speed = allocate_shared_float(nz, nr) # if evolving the electron pdf, it will be a function of the vth-normalised peculiar velocity @@ -363,39 +364,40 @@ function create_moments_neutral(nz, nr, n_species, evolve_density, evolve_upar, dvth_dz = nothing end + n_sources = length(neutral_source_settings) if any(x -> x.active, neutral_source_settings) - external_source_amplitude = allocate_shared_float(nz, nr) + external_source_amplitude = allocate_shared_float(nz, nr, n_sources) if evolve_density - external_source_density_amplitude = allocate_shared_float(nz, nr) + external_source_density_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_density_amplitude = allocate_shared_float(1, 1) + external_source_density_amplitude = allocate_shared_float(1, 1, n_sources) end if evolve_upar - external_source_momentum_amplitude = allocate_shared_float(nz, nr) + external_source_momentum_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_momentum_amplitude = allocate_shared_float(1, 1) + external_source_momentum_amplitude = allocate_shared_float(1, 1, n_sources) end if evolve_ppar - external_source_pressure_amplitude = allocate_shared_float(nz, nr) + external_source_pressure_amplitude = allocate_shared_float(nz, nr, n_sources) else - external_source_pressure_amplitude = allocate_shared_float(1, 1) + external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) end if neutral_source_settings.PI_density_controller_I != 0.0 && neutral_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") if neutral_source_settings.source_type == "density_profile_control" - external_source_controller_integral = allocate_shared_float(nz, nr) + external_source_controller_integral = allocate_shared_float(nz, nr, n_sources) else - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end else - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end else - external_source_amplitude = allocate_shared_float(1, 1) - external_source_density_amplitude = allocate_shared_float(1, 1) - external_source_momentum_amplitude = allocate_shared_float(1, 1) - external_source_pressure_amplitude = allocate_shared_float(1, 1) - external_source_controller_integral = allocate_shared_float(1, 1) + external_source_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_density_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_momentum_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) + external_source_controller_integral = allocate_shared_float(1, 1, n_sources) end if evolve_density || evolve_upar || evolve_ppar From 2e5527ab3ffee1ecdbfd14456635edbe0784ce31 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 10:16:49 +0100 Subject: [PATCH 07/23] Update file_io.jl to write out external_source_amplitudes for each moment and species as an r, z, n array instead of just r, z. --- moment_kinetics/src/file_io.jl | 37 ++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl index 2515be182..0951973a6 100644 --- a/moment_kinetics/src/file_io.jl +++ b/moment_kinetics/src/file_io.jl @@ -1332,14 +1332,15 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, units="") ion_source_settings = external_source_settings.ion - if ion_source_settings.active + if any(x -> x.active, ion_source_settings) + n = length(ion_source_settings) external_source_amplitude = create_dynamic_variable!( - dynamic, "external_source_amplitude", mk_float, z, r; + dynamic, "external_source_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external source for ions", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_density_amplitude", mk_float, z, r; + dynamic, "external_source_density_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external density source for ions", units="n_ref*c_ref/L_ref") else @@ -1347,7 +1348,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_upar external_source_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_momentum_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external momentum source for ions", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1355,7 +1356,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_ppar external_source_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_pressure_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external pressure source for ions", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1365,7 +1366,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, ion_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") if ion_source_settings.source_type == "density_profile_control" external_source_controller_integral = create_dynamic_variable!( - dynamic, "external_source_controller_integral", mk_float, z, r; + dynamic, "external_source_controller_integral", mk_float, z, r, n; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for ions") else @@ -1544,21 +1545,22 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi units="c_ref") electron_source_settings = external_source_settings.electron - if electron_source_settings.active + if any(x -> x.active, electron_source_settings) + n = length(electron_source_settings) external_source_electron_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_amplitude", mk_float, z, r; + dynamic, "external_source_electron_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external source for electrons", units="n_ref/c_ref^3*c_ref/L_ref") external_source_electron_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_density_amplitude", mk_float, z, r; + dynamic, "external_source_electron_density_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external density source for electrons", units="n_ref*c_ref/L_ref") external_source_electron_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external momentum source for electrons", units="m_ref*n_ref*c_ref*c_ref/L_ref") external_source_electron_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external pressure source for electrons", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1752,14 +1754,15 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo units="c_ref") neutral_source_settings = external_source_settings.neutral - if n_neutral_species > 0 && neutral_source_settings.active + if n_neutral_species > 0 && any(x -> x.active, neutral_source_settings) + n = length(neutral_source_settings) external_source_neutral_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external source for neutrals", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_neutral_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_density_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_density_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external density source for neutrals", units="n_ref*c_ref/L_ref") else @@ -1767,7 +1770,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_upar external_source_neutral_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external momentum source for neutrals", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1775,7 +1778,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_ppar external_source_neutral_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r, n; parallel_io=parallel_io, description="Amplitude of the external pressure source for neutrals", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1785,7 +1788,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo neutral_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") if neutral_source_settings.source_type == "density_profile_control" external_source_neutral_controller_integral = create_dynamic_variable!( - dynamic, "external_source_neutral_controller_integral", mk_float, z, r; + dynamic, "external_source_neutral_controller_integral", mk_float, z, r, n; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for neutrals") else From 0cf1c36af6158d73df48d278db7eb5d52dec1c71 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 10:37:42 +0100 Subject: [PATCH 08/23] Update load_data.jl for multiple sources --- moment_kinetics/src/load_data.jl | 72 +++++++++++++++++--------------- 1 file changed, 39 insertions(+), 33 deletions(-) diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl index 936301dc1..c2e6a53e7 100644 --- a/moment_kinetics/src/load_data.jl +++ b/moment_kinetics/src/load_data.jl @@ -4336,28 +4336,30 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dppar_dz = get_z_derivative(run_info, "parallel_pressure") dvth_dz = get_z_derivative(run_info, "thermal_speed") dqpar_dz = get_z_derivative(run_info, "parallel_heat_flux") - if run_info.external_source_settings.ion.active + if any(x -> x.active, run_info.external_source_settings.ion) + n_sources = length(run_info.external_source_settings.ion) external_source_amplitude = get_variable(run_info, "external_source_amplitude") if run_info.evolve_density external_source_density_amplitude = get_variable(run_info, "external_source_density_amplitude") else - external_source_density_amplitude = zeros(0,0,run_info.nt) + external_source_density_amplitude = zeros(0,0,n_sources,run_info.nt) end if run_info.evolve_upar external_source_momentum_amplitude = get_variable(run_info, "external_source_momentum_amplitude") else - external_source_momentum_amplitude = zeros(0,0,run_info.nt) + external_source_momentum_amplitude = zeros(0,0,n_sources,run_info.nt) end if run_info.evolve_ppar external_source_pressure_amplitude = get_variable(run_info, "external_source_pressure_amplitude") else - external_source_pressure_amplitude = zeros(0,0,run_info.nt) + external_source_pressure_amplitude = zeros(0,0,n_sources,run_info.nt) end else - external_source_amplitude = zeros(0,0,run_info.nt) - external_source_density_amplitude = zeros(0,0,run_info.nt) - external_source_momentum_amplitude = zeros(0,0,run_info.nt) - external_source_pressure_amplitude = zeros(0,0,run_info.nt) + n_sources = 0 + external_source_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_density_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_momentum_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_pressure_amplitude = zeros(0,0,n_sources,run_info.nt) end nz, nr, nspecies, nt = size(vth) @@ -4389,10 +4391,10 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dvth_dz=dvth_dz[:,:,:,it], dqpar_dz=dqpar_dz[:,:,:,it], vth=vth[:,:,:,it], - external_source_amplitude=external_source_amplitude[:,:,it], - external_source_density_amplitude=external_source_density_amplitude[:,:,it], - external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,it], - external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,it]), + external_source_amplitude=external_source_amplitude[:,:,n_sources,it], + external_source_density_amplitude=external_source_density_amplitude[:,:,n_sources,it], + external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,n_sources,it], + external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,n_sources,it]), evolve_density=run_info.evolve_density, evolve_upar=run_info.evolve_upar, evolve_ppar=run_info.evolve_ppar) @@ -4475,16 +4477,18 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dppar_dz = get_z_derivative(run_info, "electron_parallel_pressure") dvth_dz = get_z_derivative(run_info, "electron_thermal_speed") dqpar_dz = get_z_derivative(run_info, "electron_parallel_heat_flux") - if run_info.external_source_settings.electron.active + if any(x -> x.active, run_info.external_source_settings.electron) + n_sources = length(run_info.external_source_settings.electron) external_source_amplitude = get_variable(run_info, "external_source_electron_amplitude") external_source_density_amplitude = get_variable(run_info, "external_source_electron_density_amplitude") external_source_momentum_amplitude = get_variable(run_info, "external_source_electron_momentum_amplitude") external_source_pressure_amplitude = get_variable(run_info, "external_source_electron_pressure_amplitude") else - external_source_amplitude = zeros(0,0,run_info.nt) - external_source_density_amplitude = zeros(0,0,run_info.nt) - external_source_momentum_amplitude = zeros(0,0,run_info.nt) - external_source_pressure_amplitude = zeros(0,0,run_info.nt) + n_sources = 0 + external_source_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_density_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_momentum_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_pressure_amplitude = zeros(0,0,n_sources,run_info.nt) end nz, nr, nt = size(vth) @@ -4508,10 +4512,10 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dppar_dz=dppar_dz[:,:,it], dqpar_dz=dqpar_dz[:,:,it], dvth_dz=dvth_dz[:,:,it], - external_source_amplitude=external_source_amplitude[:,:,it], - external_source_density_amplitude=external_source_density_amplitude[:,:,it], - external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,it], - external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,it]),) + external_source_amplitude=external_source_amplitude[:,:,n_sources,it], + external_source_density_amplitude=external_source_density_amplitude[:,:,n_sources,it], + external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,n_sources,it], + external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,n_sources,it]),) @views update_electron_speed_vpa!(advect, density[:,:,it], upar[:,:,it], ppar[:,:,it], moments, run_info.vpa.grid, run_info.external_source_settings.electron) @@ -4593,28 +4597,30 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dpz_dz = get_z_derivative(run_info, "pz_neutral") dvth_dz = get_z_derivative(run_info, "thermal_speed_neutral") dqz_dz = get_z_derivative(run_info, "qz_neutral") - if run_info.external_source_settings.neutral.active + if any(x -> x.active, run_info.external_source_settings.neutral) + n_sources = length(run_info.external_source_settings.neutral) external_source_amplitude = get_variable(run_info, "external_source_neutral_amplitude") if run_info.evolve_density external_source_density_amplitude = get_variable(run_info, "external_source_neutral_density_amplitude") else - external_source_density_amplitude = zeros(0,0,run_info.nt) + external_source_density_amplitude = zeros(0,0,n_sources,run_info.nt) end if run_info.evolve_upar external_source_momentum_amplitude = get_variable(run_info, "external_source_neutral_momentum_amplitude") else - external_source_momentum_amplitude = zeros(0,0,run_info.nt) + external_source_momentum_amplitude = zeros(0,0,n_sources,run_info.nt) end if run_info.evolve_ppar external_source_pressure_amplitude = get_variable(run_info, "external_source_neutral_pressure_amplitude") else - external_source_pressure_amplitude = zeros(0,0,run_info.nt) + external_source_pressure_amplitude = zeros(0,0,n_sources,run_info.nt) end else - external_source_amplitude = zeros(0,0,run_info.nt) - external_source_density_amplitude = zeros(0,0,run_info.nt) - external_source_momentum_amplitude = zeros(0,0,run_info.nt) - external_source_pressure_amplitude = zeros(0,0,run_info.nt) + n_sources = 0 + external_source_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_density_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_momentum_amplitude = zeros(0,0,n_sources,run_info.nt) + external_source_pressure_amplitude = zeros(0,0,n_sources,run_info.nt) end nz, nr, nspecies, nt = size(vth) @@ -4644,10 +4650,10 @@ function get_variable(run_info, variable_name; normalize_advection_speed_shape=t dvth_dz=dvth_dz[:,:,:,it], dqz_dz=dqz_dz[:,:,:,it], vth=vth[:,:,:,it], - external_source_amplitude=external_source_amplitude[:,:,it], - external_source_density_amplitude=external_source_density_amplitude[:,:,it], - external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,it], - external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,it]), + external_source_amplitude=external_source_amplitude[:,:,n_sources,it], + external_source_density_amplitude=external_source_density_amplitude[:,:,n_sources,it], + external_source_momentum_amplitude=external_source_momentum_amplitude[:,:,n_sources,it], + external_source_pressure_amplitude=external_source_pressure_amplitude[:,:,n_sources,it]), evolve_density=run_info.evolve_density, evolve_upar=run_info.evolve_upar, evolve_ppar=run_info.evolve_ppar) From 0362f2155fc308c5aaaabe8a5581af386ba8158a Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 11:44:18 +0100 Subject: [PATCH 09/23] Update source functionality for all moment kinetic advances of density, flow and pressure for each species --- moment_kinetics/src/continuity.jl | 24 +++--- .../src/electron_fluid_equations.jl | 26 ++++--- .../src/electron_kinetic_equation.jl | 22 +++--- moment_kinetics/src/electron_vpa_advection.jl | 29 ++++---- moment_kinetics/src/energy_equation.jl | 20 +++-- moment_kinetics/src/force_balance.jl | 24 +++--- moment_kinetics/src/neutral_vz_advection.jl | 63 ++++++++-------- moment_kinetics/src/source_terms.jl | 74 ++++++++++--------- moment_kinetics/src/vpa_advection.jl | 66 +++++++++-------- 9 files changed, 193 insertions(+), 155 deletions(-) diff --git a/moment_kinetics/src/continuity.jl b/moment_kinetics/src/continuity.jl index 26e9e1314..6a7c328e8 100644 --- a/moment_kinetics/src/continuity.jl +++ b/moment_kinetics/src/continuity.jl @@ -30,11 +30,13 @@ function continuity_equation!(dens_out, fvec_in, moments, composition, dt, spect end end - if ion_source_settings.active - source_amplitude = moments.ion.external_source_density_amplitude - @loop_s_r_z is ir iz begin - dens_out[iz,ir,is] += - dt * source_amplitude[iz,ir] + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active + @views source_amplitude = moments.ion.external_source_density_amplitude[:, :, index] + @loop_s_r_z is ir iz begin + dens_out[iz,ir,is] += + dt * source_amplitude[iz,ir] + end end end @@ -71,11 +73,13 @@ function neutral_continuity_equation!(dens_out, fvec_in, moments, composition, d end end - if neutral_source_settings.active - source_amplitude = moments.neutral.external_source_density_amplitude - @loop_s_r_z is ir iz begin - dens_out[iz,ir,is] += - dt * source_amplitude[iz,ir] + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_amplitude = moments.neutral.external_source_density_amplitude[:, :, index] + @loop_s_r_z is ir iz begin + dens_out[iz,ir,is] += + dt * source_amplitude[iz,ir] + end end end diff --git a/moment_kinetics/src/electron_fluid_equations.jl b/moment_kinetics/src/electron_fluid_equations.jl index 395a94f4c..1ddc44046 100644 --- a/moment_kinetics/src/electron_fluid_equations.jl +++ b/moment_kinetics/src/electron_fluid_equations.jl @@ -225,13 +225,15 @@ function electron_energy_equation!(ppar_out, ppar_in, electron_density, electron end end - if electron_source_settings.active - pressure_source_amplitude = moments.external_source_pressure_amplitude - density_source_amplitude = moments.external_source_density_amplitude - @loop_r_z ir iz begin - ppar_out[iz,ir] += dt * (2.0 * pressure_source_amplitude[iz,ir] - - T_in[iz,ir] * density_source_amplitude[iz,ir]) / - electron_density[iz,ir] + for index ∈ eachindex(electron_source_settings) + if electron_source_settings[index].active + @views pressure_source_amplitude = moments.external_source_pressure_amplitude[:, :, index] + @views density_source_amplitude = moments.external_source_density_amplitude[:, :, index] + @loop_r_z ir iz begin + ppar_out[iz,ir] += dt * (2.0 * pressure_source_amplitude[iz,ir] + - T_in[iz,ir] * density_source_amplitude[iz,ir]) / + electron_density[iz,ir] + end end end @@ -306,10 +308,12 @@ function electron_energy_equation!(ppar_out, ppar_in, electron_density, electron end end - if electron_source_settings.active - source_amplitude = moments.external_source_pressure_amplitude - @loop_r_z ir iz begin - ppar_out[iz,ir] += dt * source_amplitude[iz,ir] + for index ∈ eachindex(electron_source_settings) + if electron_source_settings[index].active + @views source_amplitude = moments.external_source_pressure_amplitude[:, :, index] + @loop_r_z ir iz begin + ppar_out[iz,ir] += dt * source_amplitude[iz,ir] + end end end end diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl index e88476c29..ce89fdc81 100644 --- a/moment_kinetics/src/electron_kinetic_equation.jl +++ b/moment_kinetics/src/electron_kinetic_equation.jl @@ -2273,16 +2273,18 @@ function add_contribution_from_pdf_term!(pdf_out, pdf_in, ppar, dens, upar, mome end end - if electron_source_settings.active - source_density_amplitude = moments.electron.external_source_density_amplitude - source_momentum_amplitude = moments.electron.external_source_momentum_amplitude - source_pressure_amplitude = moments.electron.external_source_pressure_amplitude - @loop_r_z ir iz begin - term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - - (0.5 * source_pressure_amplitude[iz,ir] + - source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) - @loop_vperp_vpa ivperp ivpa begin - pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + for index ∈ eachindex(electron_source_settings) + if electron_source_settings[index].active + @views source_density_amplitude = moments.electron.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.electron.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.electron.external_source_pressure_amplitude[:, :, index] + @loop_r_z ir iz begin + term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - + (0.5 * source_pressure_amplitude[iz,ir] + + source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + end end end end diff --git a/moment_kinetics/src/electron_vpa_advection.jl b/moment_kinetics/src/electron_vpa_advection.jl index 2d7968286..5f8abe6b5 100644 --- a/moment_kinetics/src/electron_vpa_advection.jl +++ b/moment_kinetics/src/electron_vpa_advection.jl @@ -66,19 +66,22 @@ function update_electron_speed_vpa!(advect, density, upar, ppar, moments, vpa, advect.speed[ivpa,ivperp,iz,ir] = ((vth[iz,ir] * dppar_dz[iz,ir] + vpa[ivpa] * dqpar_dz[iz,ir]) / (2 * ppar[iz,ir]) - vpa[ivpa]^2 * dvth_dz[iz,ir]) end - if electron_source_settings.active - source_density_amplitude = moments.electron.external_source_density_amplitude - source_momentum_amplitude = moments.electron.external_source_momentum_amplitude - source_pressure_amplitude = moments.electron.external_source_pressure_amplitude - @loop_r_z ir iz begin - term1 = source_density_amplitude[iz,ir] * upar[iz,ir]/(density[iz,ir]*vth[iz,ir]) - term2_over_vpa = - -0.5 * (source_pressure_amplitude[iz,ir] + - 2.0 * upar[iz,ir] * source_momentum_amplitude[iz,ir]) / - ppar[iz,ir] + - 0.5 * source_density_amplitude[iz,ir] / density[iz,ir] - @loop_vperp_vpa ivperp ivpa begin - advect.speed[ivpa,ivperp,iz,ir] += term1 + vpa[ivpa] * term2_over_vpa + + for index ∈ eachindex(electron_source_settings) + if electron_source_settings[index].active + @views source_density_amplitude = moments.electron.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.electron.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.electron.external_source_pressure_amplitude[:, :, index] + @loop_r_z ir iz begin + term1 = source_density_amplitude[iz,ir] * upar[iz,ir]/(density[iz,ir]*vth[iz,ir]) + term2_over_vpa = + -0.5 * (source_pressure_amplitude[iz,ir] + + 2.0 * upar[iz,ir] * source_momentum_amplitude[iz,ir]) / + ppar[iz,ir] + + 0.5 * source_density_amplitude[iz,ir] / density[iz,ir] + @loop_vperp_vpa ivperp ivpa begin + advect.speed[ivpa,ivperp,iz,ir] += term1 + vpa[ivpa] * term2_over_vpa + end end end end diff --git a/moment_kinetics/src/energy_equation.jl b/moment_kinetics/src/energy_equation.jl index 8890be3b6..7e708a61e 100644 --- a/moment_kinetics/src/energy_equation.jl +++ b/moment_kinetics/src/energy_equation.jl @@ -23,10 +23,12 @@ function energy_equation!(ppar, fvec, moments, collisions, dt, spectral, composi end - if ion_source_settings.active - source_amplitude = moments.ion.external_source_pressure_amplitude - @loop_s_r_z is ir iz begin - ppar[iz,ir,is] += dt * source_amplitude[iz,ir] + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active + @views source_amplitude = moments.ion.external_source_pressure_amplitude[:, :, index] + @loop_s_r_z is ir iz begin + ppar[iz,ir,is] += dt * source_amplitude[iz,ir] + end end end @@ -75,10 +77,12 @@ function neutral_energy_equation!(pz, fvec, moments, collisions, dt, spectral, - 3.0*fvec.pz_neutral[iz,ir,isn]*moments.neutral.duz_dz[iz,ir,isn]) end - if neutral_source_settings.active - source_amplitude = moments.neutral.external_source_pressure_amplitude - @loop_s_r_z isn ir iz begin - pz[iz,ir,isn] += dt * source_amplitude[iz,ir] + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_amplitude = moments.neutral.external_source_pressure_amplitude[:, :, index] + @loop_s_r_z isn ir iz begin + pz[iz,ir,isn] += dt * source_amplitude[iz,ir] + end end end diff --git a/moment_kinetics/src/force_balance.jl b/moment_kinetics/src/force_balance.jl index 902c3901e..d96301711 100644 --- a/moment_kinetics/src/force_balance.jl +++ b/moment_kinetics/src/force_balance.jl @@ -29,11 +29,13 @@ function force_balance!(pflx, density_out, fvec, moments, fields, collisions, dt 0.5*geometry.bzed[iz,ir]*fields.Ez[iz,ir]*density[iz,ir,is]) end - if ion_source_settings.active && false - source_amplitude = moments.ion.external_source_momentum_amplitude - @loop_s_r_z is ir iz begin - pflx[iz,ir,is] += - dt * source_amplitude[iz,ir] + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active && false + @views source_amplitude = moments.ion.external_source_momentum_amplitude[:, :, index] + @loop_s_r_z is ir iz begin + pflx[iz,ir,is] += + dt * source_amplitude[iz,ir] + end end end @@ -83,11 +85,13 @@ function neutral_force_balance!(pflx, density_out, fvec, moments, fields, collis 2.0*density[iz,ir,isn]*uz[iz,ir,isn]*moments.neutral.duz_dz_upwind[iz,ir,isn]) end - if neutral_source_settings.active && false - source_amplitude = moments.neutral.external_source_momentum_amplitude - @loop_sn_r_z isn ir iz begin - pflx[iz,ir,isn] += - dt * source_amplitude[iz,ir] + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active && false + @views source_amplitude = moments.neutral.external_source_momentum_amplitude[:, :, index] + @loop_sn_r_z isn ir iz begin + pflx[iz,ir,isn] += + dt * source_amplitude[iz,ir] + end end end diff --git a/moment_kinetics/src/neutral_vz_advection.jl b/moment_kinetics/src/neutral_vz_advection.jl index 96b61cd39..a0f9dbdc5 100644 --- a/moment_kinetics/src/neutral_vz_advection.jl +++ b/moment_kinetics/src/neutral_vz_advection.jl @@ -126,25 +126,28 @@ function update_speed_n_u_p_evolution_neutral!(advect, fvec, moments, vz, z, r, end end end - if neutral_source_settings.active - source_density_amplitude = moments.neutral.external_source_density_amplitude - source_momentum_amplitude = moments.neutral.external_source_momentum_amplitude - source_pressure_amplitude = moments.neutral.external_source_pressure_amplitude - density = fvec.density_neutral - uz = fvec.uz_neutral - pz = fvec.pz_neutral - vth = moments.neutral.vth - vz_grid = vz.grid - @loop_s_r_z is ir iz begin - term1 = source_density_amplitude[iz,ir] * uz[iz,ir,is]/(density[iz,ir,is]*vth[iz,ir,is]) - term2_over_vpa = - -0.5 * (source_pressure_amplitude[iz,ir] + - 2.0 * uz[iz,ir,is] * source_momentum_amplitude[iz,ir]) / - pz[iz,ir,is] + - 0.5 * source_density_amplitude[iz,ir] / density[iz,ir,is] - @loop_vzeta_vr_vz ivzeta ivr ivz begin - advect[is].speed[ivz,ivr,ivzeta,iz,ir] += term1 + - vz_grid[ivz] * term2_over_vpa + + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_density_amplitude = moments.neutral.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.neutral.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.neutral.external_source_pressure_amplitude[:, :, index] + density = fvec.density_neutral + uz = fvec.uz_neutral + pz = fvec.pz_neutral + vth = moments.neutral.vth + vz_grid = vz.grid + @loop_s_r_z is ir iz begin + term1 = source_density_amplitude[iz,ir] * uz[iz,ir,is]/(density[iz,ir,is]*vth[iz,ir,is]) + term2_over_vpa = + -0.5 * (source_pressure_amplitude[iz,ir] + + 2.0 * uz[iz,ir,is] * source_momentum_amplitude[iz,ir]) / + pz[iz,ir,is] + + 0.5 * source_density_amplitude[iz,ir] / density[iz,ir,is] + @loop_vzeta_vr_vz ivzeta ivr ivz begin + advect[is].speed[ivz,ivr,ivzeta,iz,ir] += term1 + + vz_grid[ivz] * term2_over_vpa + end end end end @@ -184,7 +187,7 @@ function update_speed_n_p_evolution_neutral!(advect, fields, fvec, moments, vz, end end end - if ion_source_settings.active + if any(x -> x.active, neutral_source_settings) error("External source not implemented for evolving n and ppar case") end end @@ -219,15 +222,17 @@ function update_speed_n_u_evolution_neutral!(advect, fvec, moments, vz, z, r, co end end end - if neutral_source_settings.active - source_density_amplitude = moments.neutral.external_source_density_amplitude - density = fvec.density_neutral - uz = fvec.uz_neutral - vth = moments.neutral.vth - @loop_sn_r_z isn ir iz begin - term = source_density_amplitude[iz,ir] * uz[iz,ir,isn] / density[iz,ir,isn] - @loop_vzeta_vr_vz ivzeta ivr ivz begin - advect[isn].speed[ivz,ivr,ivzeta,iz,ir] += term + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_density_amplitude = moments.neutral.external_source_density_amplitude[:, :, index] + density = fvec.density_neutral + uz = fvec.uz_neutral + vth = moments.neutral.vth + @loop_sn_r_z isn ir iz begin + term = source_density_amplitude[iz,ir] * uz[iz,ir,isn] / density[iz,ir,isn] + @loop_vzeta_vr_vz ivzeta ivr ivz begin + advect[isn].speed[ivz,ivr,ivzeta,iz,ir] += term + end end end end diff --git a/moment_kinetics/src/source_terms.jl b/moment_kinetics/src/source_terms.jl index 9850d6f23..8a0a1dbe8 100644 --- a/moment_kinetics/src/source_terms.jl +++ b/moment_kinetics/src/source_terms.jl @@ -64,12 +64,14 @@ function source_terms_evolve_density!(pdf_out, pdf_in, dens, upar, ddens_dz, dup end end - if ion_source_settings.active - source_density_amplitude = moments.ion.external_source_density_amplitude - @loop_r_z ir iz begin - term = dt * source_density_amplitude[iz,ir] / dens[iz,ir] - @loop_vperp_vpa ivperp ivpa begin - pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active + @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] + @loop_r_z ir iz begin + term = dt * source_density_amplitude[iz,ir] / dens[iz,ir] + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + end end end end @@ -96,16 +98,18 @@ function source_terms_evolve_ppar_no_collisions!(pdf_out, pdf_in, dens, upar, pp end end - if ion_source_settings.active - source_density_amplitude = moments.ion.external_source_density_amplitude - source_momentum_amplitude = moments.ion.external_source_momentum_amplitude - source_pressure_amplitude = moments.ion.external_source_pressure_amplitude - @loop_r_z ir iz begin - term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - - (0.5 * source_pressure_amplitude[iz,ir] + - source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) - @loop_vperp_vpa ivperp ivpa begin - pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + for index ∈ eachindex(ion_source_settings) + if ion_source_settings[index].active + @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.ion.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.ion.external_source_pressure_amplitude[:, :, index] + @loop_r_z ir iz begin + term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - + (0.5 * source_pressure_amplitude[iz,ir] + + source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) + @loop_vperp_vpa ivperp ivpa begin + pdf_out[ivpa,ivperp,iz,ir] -= term * pdf_in[ivpa,ivperp,iz,ir] + end end end end @@ -192,12 +196,14 @@ function source_terms_evolve_density_neutral!(pdf_out, pdf_in, dens, upar, ddens end end - if neutral_source_settings.active - source_density_amplitude = moments.neutral.external_source_density_amplitude - @loop_r_z ir iz begin - term = dt * source_density_amplitude[iz,ir] / dens[iz,ir] - @loop_vzeta_vr_vz ivzeta ivr ivz begin - pdf_out[ivz,ivr,ivzeta,iz,ir] -= term * pdf_in[ivz,ivr,ivzeta,iz,ir] + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_density_amplitude = moments.neutral.external_source_density_amplitude[:, :, index] + @loop_r_z ir iz begin + term = dt * source_density_amplitude[iz,ir] / dens[iz,ir] + @loop_vzeta_vr_vz ivzeta ivr ivz begin + pdf_out[ivz,ivr,ivzeta,iz,ir] -= term * pdf_in[ivz,ivr,ivzeta,iz,ir] + end end end end @@ -223,20 +229,22 @@ function source_terms_evolve_ppar_no_collisions_neutral!(pdf_out, pdf_in, dens, end end - if neutral_source_settings.active - source_density_amplitude = moments.neutral.external_source_density_amplitude - source_momentum_amplitude = moments.neutral.external_source_momentum_amplitude - source_pressure_amplitude = moments.neutral.external_source_pressure_amplitude - @loop_r_z ir iz begin - term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - - (0.5 * source_pressure_amplitude[iz,ir] + - source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) - @loop_vzeta_vr_vz ivzeta ivr ivz begin - pdf_out[ivz,ivr,ivzeta,iz,ir] -= term * pdf_in[ivz,ivr,ivzeta,iz,ir] + for index ∈ eachindex(neutral_source_settings) + if neutral_source_settings[index].active + @views source_density_amplitude = moments.neutral.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.neutral.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.neutral.external_source_pressure_amplitude[:, :, index] + @loop_r_z ir iz begin + term = dt * (1.5 * source_density_amplitude[iz,ir] / dens[iz,ir] - + (0.5 * source_pressure_amplitude[iz,ir] + + source_momentum_amplitude[iz,ir]) / ppar[iz,ir]) + @loop_vzeta_vr_vz ivzeta ivr ivz begin + pdf_out[ivz,ivr,ivzeta,iz,ir] -= term * pdf_in[ivz,ivr,ivzeta,iz,ir] + end end end end - + return nothing end diff --git a/moment_kinetics/src/vpa_advection.jl b/moment_kinetics/src/vpa_advection.jl index bd70b3503..328395e9a 100644 --- a/moment_kinetics/src/vpa_advection.jl +++ b/moment_kinetics/src/vpa_advection.jl @@ -433,24 +433,26 @@ function update_speed_n_u_p_evolution!(advect, fvec, moments, vpa, z, r, composi end end end - if ion_source_settings.active - source_density_amplitude = moments.ion.external_source_density_amplitude - source_momentum_amplitude = moments.ion.external_source_momentum_amplitude - source_pressure_amplitude = moments.ion.external_source_pressure_amplitude - density = fvec.density - upar = fvec.upar - ppar = fvec.ppar - vth = moments.ion.vth - vpa_grid = vpa.grid - @loop_s_r_z is ir iz begin - term1 = source_density_amplitude[iz,ir] * upar[iz,ir,is]/(density[iz,ir,is]*vth[iz,ir,is]) - term2_over_vpa = - -0.5 * (source_pressure_amplitude[iz,ir] + - 2.0 * upar[iz,ir,is] * source_momentum_amplitude[iz,ir]) / - ppar[iz,ir,is] + - 0.5 * source_density_amplitude[iz,ir] / density[iz,ir,is] - @loop_vperp_vpa ivperp ivpa begin - advect[is].speed[ivpa,ivperp,iz,ir] += term1 + vpa_grid[ivpa] * term2_over_vpa + for index ∈ eachindex(ion_source_settings.ion) + if ion_source_settings[index].active + @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] + @views source_momentum_amplitude = moments.ion.external_source_momentum_amplitude[:, :, index] + @views source_pressure_amplitude = moments.ion.external_source_pressure_amplitude[:, :, index] + density = fvec.density + upar = fvec.upar + ppar = fvec.ppar + vth = moments.ion.vth + vpa_grid = vpa.grid + @loop_s_r_z is ir iz begin + term1 = source_density_amplitude[iz,ir] * upar[iz,ir,is]/(density[iz,ir,is]*vth[iz,ir,is]) + term2_over_vpa = + -0.5 * (source_pressure_amplitude[iz,ir] + + 2.0 * upar[iz,ir,is] * source_momentum_amplitude[iz,ir]) / + ppar[iz,ir,is] + + 0.5 * source_density_amplitude[iz,ir] / density[iz,ir,is] + @loop_vperp_vpa ivperp ivpa begin + advect[is].speed[ivpa,ivperp,iz,ir] += term1 + vpa_grid[ivpa] * term2_over_vpa + end end end end @@ -494,7 +496,7 @@ function update_speed_n_p_evolution!(advect, fields, fvec, moments, vpa, z, r, end end end - if ion_source_settings.active + if any(x -> x.active, ion_source_settings) error("External source not implemented for evolving n and ppar case") end end @@ -538,18 +540,20 @@ function update_speed_n_u_evolution!(advect, fvec, moments, vpa, z, r, compositi end end end - if ion_source_settings.active - source_density_amplitude = moments.ion.external_source_density_amplitude - source_strength = ion_source_settings.source_strength - r_amplitude = ion_source_settings.r_amplitude - z_amplitude = ion_source_settings.z_amplitude - density = fvec.density - upar = fvec.upar - vth = moments.ion.vth - @loop_s_r_z is ir iz begin - term = source_density_amplitude[iz,ir] * upar[iz,ir,is] / density[iz,ir,is] - @loop_vperp_vpa ivperp ivpa begin - advect[is].speed[ivpa,ivperp,iz,ir] += term + for index ∈ eachindex(ion_source_settings.ion) + if ion_source_settings[index].active + @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] + source_strength = ion_source_settings[index].source_strength + r_amplitude = ion_source_settings[index].r_amplitude + z_amplitude = ion_source_settings[index].z_amplitude + density = fvec.density + upar = fvec.upar + vth = moments.ion.vth + @loop_s_r_z is ir iz begin + term = source_density_amplitude[iz,ir] * upar[iz,ir,is] / density[iz,ir,is] + @loop_vperp_vpa ivperp ivpa begin + advect[is].speed[ivpa,ivperp,iz,ir] += term + end end end end From 41759d9a3f9ef75b7241837c9e3ceebd5f6ac478 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 11:48:58 +0100 Subject: [PATCH 10/23] clean up external_sources.jl --- moment_kinetics/src/external_sources.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 2fc218aec..c3a7d8b7e 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -507,7 +507,6 @@ function initialize_external_source_amplitude!(moments, external_source_settings end # now do same for electron sources, which (if present) are mostly mirrors of ion sources - println(electron_source_settings) if electron_source_settings[index].active if electron_source_settings[index].source_type == "energy" @loop_r_z ir iz begin From 2622a24ec073a2808d77bad7cf4db29b753cd5fd Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 16:08:26 +0100 Subject: [PATCH 11/23] Clean up small typos --- moment_kinetics/src/external_sources.jl | 2 +- moment_kinetics/src/vpa_advection.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index c3a7d8b7e..48546a02f 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -364,7 +364,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics) else # this is not very efficient because we're copying the same electron_settings # for each ion source - push!(electron_sources, get_settings_electrons(inactive_ion_source)) + push!(electron_sources, get_settings_electrons(get_settings_ions(1, false))) end # put all neutral sources into neutral_source_data struct vector diff --git a/moment_kinetics/src/vpa_advection.jl b/moment_kinetics/src/vpa_advection.jl index 328395e9a..32f43007a 100644 --- a/moment_kinetics/src/vpa_advection.jl +++ b/moment_kinetics/src/vpa_advection.jl @@ -433,7 +433,7 @@ function update_speed_n_u_p_evolution!(advect, fvec, moments, vpa, z, r, composi end end end - for index ∈ eachindex(ion_source_settings.ion) + for index ∈ eachindex(ion_source_settings) if ion_source_settings[index].active @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] @views source_momentum_amplitude = moments.ion.external_source_momentum_amplitude[:, :, index] @@ -540,7 +540,7 @@ function update_speed_n_u_evolution!(advect, fvec, moments, vpa, z, r, compositi end end end - for index ∈ eachindex(ion_source_settings.ion) + for index ∈ eachindex(ion_source_settings) if ion_source_settings[index].active @views source_density_amplitude = moments.ion.external_source_density_amplitude[:, :, index] source_strength = ion_source_settings[index].source_strength From 6123063d5ed901c9ce516d746f3980b748a5b663 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Sun, 8 Sep 2024 16:09:17 +0100 Subject: [PATCH 12/23] Update all example input files and test files --- .../fokker-planck-relaxation-slowing-down-alphas.toml | 2 +- ...ck-source-sink-slowing-down-alphas-no-self-collisions.toml | 2 +- .../fokker-planck-source-sink-slowing-down-alphas.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_boltzmann-vpadiss0.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_boltzmann.toml | 2 +- ...bc_recyclefraction0.5_split3_braginskii-vpadiss0-IMEX.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_kinetic-vpadiss0.toml | 2 +- .../recycling-fraction/wall-bc_recyclefraction0.5-init.toml | 2 +- examples/recycling-fraction/wall-bc_recyclefraction0.5.toml | 2 +- .../recycling-fraction/wall-bc_recyclefraction0.5_split1.toml | 2 +- .../recycling-fraction/wall-bc_recyclefraction0.5_split2.toml | 2 +- .../wall-bc_recyclefraction0.5_split3-init.toml | 2 +- .../recycling-fraction/wall-bc_recyclefraction0.5_split3.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_SSPRK4.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_fekete104.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_fekete42.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_fekete64.toml | 2 +- .../wall-bc_recyclefraction0.5_split3_rkf54.toml | 2 +- examples/wall-bc/wall-bc_no-neutrals.toml | 2 +- examples/wall-bc/wall-bc_no-neutrals_split1.toml | 2 +- examples/wall-bc/wall-bc_no-neutrals_split2.toml | 2 +- examples/wall-bc/wall-bc_no-neutrals_split3.toml | 2 +- examples/wall-bc/wall-bc_volumerecycle.toml | 4 ++-- examples/wall-bc/wall-bc_volumerecycle_split1.toml | 4 ++-- examples/wall-bc/wall-bc_volumerecycle_split2.toml | 4 ++-- examples/wall-bc/wall-bc_volumerecycle_split3.toml | 4 ++-- moment_kinetics/debug_test/kinetic_electron_inputs.jl | 2 +- moment_kinetics/debug_test/recycling_fraction_inputs.jl | 2 +- moment_kinetics/test/harrisonthompson.jl | 2 +- moment_kinetics/test/recycling_fraction_tests.jl | 2 +- .../wall-bc/wall-bc_recyclefraction0.5.toml | 2 +- .../wall-bc/wall-bc_recyclefraction0.5_split1.toml | 2 +- .../wall-bc/wall-bc_recyclefraction0.5_split2.toml | 2 +- .../wall-bc/wall-bc_recyclefraction0.5_split3.toml | 2 +- 34 files changed, 38 insertions(+), 38 deletions(-) diff --git a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml index e6b75e29c..5218fbf4f 100644 --- a/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml +++ b/examples/fokker-planck/fokker-planck-relaxation-slowing-down-alphas.toml @@ -65,7 +65,7 @@ sd_mi = 0.5 # mD/malpha sd_me = 0.000013616 # 0.25/1836.0 me/malpha sd_q = 1.0 -[ion_source] +[ion_source_1] active=false source_strength=0.0 source_T=0.005 diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml index 1a95fc210..bb69178fb 100644 --- a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas-no-self-collisions.toml @@ -64,7 +64,7 @@ sd_mi = 0.5 # mD/malpha sd_me = 0.000013616 # 0.25/1836.0 me/malpha sd_q = 1.0 -[ion_source] +[ion_source_1] active=true source_strength=1.0 source_T=0.005 diff --git a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml index cd49c82e7..fc53ee470 100644 --- a/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml +++ b/examples/fokker-planck/fokker-planck-source-sink-slowing-down-alphas.toml @@ -64,7 +64,7 @@ sd_mi = 0.5 # mD/malpha sd_me = 0.000013616 # 0.25/1836.0 me/malpha sd_q = 1.0 -[ion_source] +[ion_source_1] active=true source_strength=1.0 source_T=0.005 diff --git a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann-vpadiss0.toml b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann-vpadiss0.toml index 58998f2f0..b09ee80c6 100644 --- a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann-vpadiss0.toml +++ b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann-vpadiss0.toml @@ -74,7 +74,7 @@ nwrite_dfns = 100000 steady_state_residual = true converged_residual_value = 1.0e-3 -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann.toml b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann.toml index 4b36701ae..f993e7561 100644 --- a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann.toml +++ b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_boltzmann.toml @@ -74,7 +74,7 @@ nwrite_dfns = 100000 steady_state_residual = true converged_residual_value = 1.0e-3 -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_braginskii-vpadiss0-IMEX.toml b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_braginskii-vpadiss0-IMEX.toml index c9230a33a..c6f91785d 100644 --- a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_braginskii-vpadiss0-IMEX.toml +++ b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_braginskii-vpadiss0-IMEX.toml @@ -84,7 +84,7 @@ steady_state_residual = true converged_residual_value = 1.0e-3 #write_after_fixed_step_count = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_kinetic-vpadiss0.toml b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_kinetic-vpadiss0.toml index 0a66ed75e..9af371c6a 100644 --- a/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_kinetic-vpadiss0.toml +++ b/examples/kinetic-electrons/wall-bc_recyclefraction0.5_split3_kinetic-vpadiss0.toml @@ -90,7 +90,7 @@ minimum_dt = 1.0e-9 initialization_residual_value = 2.5 converged_residual_value = 0.1 #1.0e-3 -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5-init.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5-init.toml index 9442b1dde..6c73e7b25 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5-init.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5-init.toml @@ -75,7 +75,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5.toml index 70e81e234..85ab587b9 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5.toml @@ -80,7 +80,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split1.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split1.toml index 634959e88..60a6b8671 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split1.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split1.toml @@ -80,7 +80,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split2.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split2.toml index f1ba6bc68..770c49197 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split2.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split2.toml @@ -80,7 +80,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml index 7c6e17f90..4ea030e52 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3-init.toml @@ -75,7 +75,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml index 1352d2fdd..96dfe0d41 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3.toml @@ -80,7 +80,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_SSPRK4.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_SSPRK4.toml index 73be81551..30f0ab14d 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_SSPRK4.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_SSPRK4.toml @@ -74,7 +74,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml index 26a9fda8b..a489cfbd7 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete104.toml @@ -82,7 +82,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete42.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete42.toml index 7478a386e..f8f00d7e3 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete42.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete42.toml @@ -82,7 +82,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete64.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete64.toml index 3c86487ac..5814d393e 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete64.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_fekete64.toml @@ -82,7 +82,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_rkf54.toml b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_rkf54.toml index 8372db4cc..5534ffc3d 100644 --- a/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_rkf54.toml +++ b/examples/recycling-fraction/wall-bc_recyclefraction0.5_split3_rkf54.toml @@ -82,7 +82,7 @@ converged_residual_value = 1.0e-3 [krook_collisions] use_krook = true -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/wall-bc/wall-bc_no-neutrals.toml b/examples/wall-bc/wall-bc_no-neutrals.toml index ce5da5ba8..ab5b26595 100644 --- a/examples/wall-bc/wall-bc_no-neutrals.toml +++ b/examples/wall-bc/wall-bc_no-neutrals.toml @@ -70,7 +70,7 @@ nwrite = 5000 nwrite_dfns = 5000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/wall-bc/wall-bc_no-neutrals_split1.toml b/examples/wall-bc/wall-bc_no-neutrals_split1.toml index cae94a045..0e6086ed0 100644 --- a/examples/wall-bc/wall-bc_no-neutrals_split1.toml +++ b/examples/wall-bc/wall-bc_no-neutrals_split1.toml @@ -70,7 +70,7 @@ nwrite = 5000 nwrite_dfns = 5000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/wall-bc/wall-bc_no-neutrals_split2.toml b/examples/wall-bc/wall-bc_no-neutrals_split2.toml index ce74ec368..2159851d8 100644 --- a/examples/wall-bc/wall-bc_no-neutrals_split2.toml +++ b/examples/wall-bc/wall-bc_no-neutrals_split2.toml @@ -70,7 +70,7 @@ nwrite = 5000 nwrite_dfns = 5000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/wall-bc/wall-bc_no-neutrals_split3.toml b/examples/wall-bc/wall-bc_no-neutrals_split3.toml index 23316e611..25664c688 100644 --- a/examples/wall-bc/wall-bc_no-neutrals_split3.toml +++ b/examples/wall-bc/wall-bc_no-neutrals_split3.toml @@ -70,7 +70,7 @@ nwrite = 5000 nwrite_dfns = 5000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/examples/wall-bc/wall-bc_volumerecycle.toml b/examples/wall-bc/wall-bc_volumerecycle.toml index 1cbf7f1b2..ee6a36dd6 100644 --- a/examples/wall-bc/wall-bc_volumerecycle.toml +++ b/examples/wall-bc/wall-bc_volumerecycle.toml @@ -72,14 +72,14 @@ nwrite = 10000 nwrite_dfns = 10000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 source_strength = 2.0 source_T = 2.0 -[neutral_source] +[neutral_source_1] active = true source_type = "recycling" recycling_controller_fraction = 0.25 diff --git a/examples/wall-bc/wall-bc_volumerecycle_split1.toml b/examples/wall-bc/wall-bc_volumerecycle_split1.toml index 442ce033d..a2fba2249 100644 --- a/examples/wall-bc/wall-bc_volumerecycle_split1.toml +++ b/examples/wall-bc/wall-bc_volumerecycle_split1.toml @@ -72,14 +72,14 @@ nwrite = 10000 nwrite_dfns = 10000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 source_strength = 2.0 source_T = 2.0 -[neutral_source] +[neutral_source_1] active = true source_type = "recycling" recycling_controller_fraction = 0.25 diff --git a/examples/wall-bc/wall-bc_volumerecycle_split2.toml b/examples/wall-bc/wall-bc_volumerecycle_split2.toml index 08eb2ed6c..1a5c2e968 100644 --- a/examples/wall-bc/wall-bc_volumerecycle_split2.toml +++ b/examples/wall-bc/wall-bc_volumerecycle_split2.toml @@ -72,14 +72,14 @@ nwrite = 10000 nwrite_dfns = 10000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 source_strength = 2.0 source_T = 2.0 -[neutral_source] +[neutral_source_1] active = true source_type = "recycling" recycling_controller_fraction = 0.25 diff --git a/examples/wall-bc/wall-bc_volumerecycle_split3.toml b/examples/wall-bc/wall-bc_volumerecycle_split3.toml index 0c8dad98d..6b3e03cb9 100644 --- a/examples/wall-bc/wall-bc_volumerecycle_split3.toml +++ b/examples/wall-bc/wall-bc_volumerecycle_split3.toml @@ -72,14 +72,14 @@ nwrite = 10000 nwrite_dfns = 10000 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 source_strength = 2.0 source_T = 2.0 -[neutral_source] +[neutral_source_1] active = true source_type = "recycling" recycling_controller_fraction = 0.25 diff --git a/moment_kinetics/debug_test/kinetic_electron_inputs.jl b/moment_kinetics/debug_test/kinetic_electron_inputs.jl index 684b114ff..5b8d4f7b0 100644 --- a/moment_kinetics/debug_test/kinetic_electron_inputs.jl +++ b/moment_kinetics/debug_test/kinetic_electron_inputs.jl @@ -81,7 +81,7 @@ test_input = Dict("n_ion_species" => 1, "vz_L" => 6.0, "vz_bc" => "zero", "vz_discretization" => "chebyshev_pseudospectral", - "ion_source" => Dict{String,Any}("active" => true, + "ion_source_1" => Dict{String,Any}("active" => true, "z_profile" => "gaussian", "z_width" => 0.125, "source_strength" => 2.0, diff --git a/moment_kinetics/debug_test/recycling_fraction_inputs.jl b/moment_kinetics/debug_test/recycling_fraction_inputs.jl index 7f4b24266..a49015970 100644 --- a/moment_kinetics/debug_test/recycling_fraction_inputs.jl +++ b/moment_kinetics/debug_test/recycling_fraction_inputs.jl @@ -73,7 +73,7 @@ test_input = Dict("n_ion_species" => 1, "vz_L" => 6.0, "vz_bc" => "zero", "vz_discretization" => "chebyshev_pseudospectral", - "ion_source" => Dict("active" => true, + "ion_source_1" => Dict("active" => true, "z_profile" => "gaussian", "z_width" => 0.125, "source_strength" => 2.0, diff --git a/moment_kinetics/test/harrisonthompson.jl b/moment_kinetics/test/harrisonthompson.jl index 06206a263..4830cb216 100644 --- a/moment_kinetics/test/harrisonthompson.jl +++ b/moment_kinetics/test/harrisonthompson.jl @@ -112,7 +112,7 @@ test_input_finite_difference = Dict("n_ion_species" => 1, "vz_L" => 8.0, "vz_bc" => "zero", "vz_discretization" => "finite_difference", - "ion_source" => Dict("active" => true, + "ion_source_1" => Dict("active" => true, "source_strength" => ionization_frequency, "source_T" => 0.25, "z_profile" => "constant", diff --git a/moment_kinetics/test/recycling_fraction_tests.jl b/moment_kinetics/test/recycling_fraction_tests.jl index fed259261..8add56a1f 100644 --- a/moment_kinetics/test/recycling_fraction_tests.jl +++ b/moment_kinetics/test/recycling_fraction_tests.jl @@ -83,7 +83,7 @@ test_input = Dict("n_ion_species" => 1, "vz_L" => 18.0, "vz_bc" => "zero", "vz_discretization" => "chebyshev_pseudospectral", - "ion_source" => Dict("active" => true, + "ion_source_1" => Dict("active" => true, "z_profile" => "gaussian", "z_width" => 0.125, "source_strength" => 2.0, diff --git a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5.toml b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5.toml index 8943ac402..a25476a08 100644 --- a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5.toml +++ b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5.toml @@ -75,7 +75,7 @@ nwrite_dfns = 10000 n_rk_stages = 4 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split1.toml b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split1.toml index f7a1484d7..cf192f884 100644 --- a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split1.toml +++ b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split1.toml @@ -75,7 +75,7 @@ nwrite_dfns = 10000 n_rk_stages = 4 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split2.toml b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split2.toml index ef1a221cb..d6bba39f0 100644 --- a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split2.toml +++ b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split2.toml @@ -75,7 +75,7 @@ nwrite_dfns = 10000 n_rk_stages = 4 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 diff --git a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split3.toml b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split3.toml index 79c01a3be..0c0c2711d 100644 --- a/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split3.toml +++ b/publication_inputs/2023_EFTC_jto-poster/wall-bc/wall-bc_recyclefraction0.5_split3.toml @@ -75,7 +75,7 @@ nwrite_dfns = 10000 n_rk_stages = 4 split_operators = false -[ion_source] +[ion_source_1] active = true z_profile = "gaussian" z_width = 0.125 From dddad290b8620b37ee44301ac1f26d75dd1cdc9b Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Mon, 9 Sep 2024 09:52:53 +0100 Subject: [PATCH 13/23] Fix bug in velocity_moments.jl PI_density_controller loop --- moment_kinetics/src/velocity_moments.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/moment_kinetics/src/velocity_moments.jl b/moment_kinetics/src/velocity_moments.jl index b362b721b..096dae7b2 100644 --- a/moment_kinetics/src/velocity_moments.jl +++ b/moment_kinetics/src/velocity_moments.jl @@ -156,9 +156,9 @@ function create_moments_ion(nz, nr, n_species, evolve_density, evolve_upar, else external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) end - if ion_source_settings.PI_density_controller_I != 0.0 && - ion_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") - if ion_source_settings.source_type == "density_profile_control" + if any(x -> x.PI_density_controller_I != 0.0 && x.source_type ∈ + ("density_profile_control", "density_midpoint_control"), ion_source_settings) + if any(x -> x.source_type == "density_profile_control", ion_source_settings) external_source_controller_integral = allocate_shared_float(nz, nr, n_sources) else external_source_controller_integral = allocate_shared_float(1, 1, n_sources) @@ -382,9 +382,9 @@ function create_moments_neutral(nz, nr, n_species, evolve_density, evolve_upar, else external_source_pressure_amplitude = allocate_shared_float(1, 1, n_sources) end - if neutral_source_settings.PI_density_controller_I != 0.0 && - neutral_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") - if neutral_source_settings.source_type == "density_profile_control" + if any(x -> x.PI_density_controller_I != 0.0 && x.source_type ∈ + ("density_profile_control", "density_midpoint_control"), neutral_source_settings) + if any(x -> x.source_type == "density_profile_control", neutral_source_settings) external_source_controller_integral = allocate_shared_float(nz, nr, n_sources) else external_source_controller_integral = allocate_shared_float(1, 1, n_sources) From da0f366ee874bc1dcde966f86ca792fa19fd3acf Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Mon, 9 Sep 2024 11:08:13 +0100 Subject: [PATCH 14/23] Add exponential wall decay profile option to simulate neutral ionisation for ions when n_neutrals = 0, and fix some more small errors introduced during restructuring of source system --- moment_kinetics/src/external_sources.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 48546a02f..b53e8b71e 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -350,9 +350,7 @@ function setup_external_sources!(input_dict, r, z, electron_physics) # If there are no ion sources, add an inactive ion source to the vector if counter == 1 - inactive_ion_source = get_settings_ions(1, false) - push!(ion_sources, inactive_ion_source) - counter += 1 + push!(ion_sources, get_settings_ions(1, false)) end # put all electron sources into electron_source_data struct vector, where @@ -360,11 +358,9 @@ function setup_external_sources!(input_dict, r, z, electron_physics) electron_sources = electron_source_data[] if electron_physics ∈ (braginskii_fluid, kinetic_electrons, kinetic_electrons_with_temperature_equation) - electron_sources .= get_settings_electrons.(ion_sources) + electron_sources = [get_settings_electrons(ion_sources[i]) for i ∈ 1:counter] else - # this is not very efficient because we're copying the same electron_settings - # for each ion source - push!(electron_sources, get_settings_electrons(get_settings_ions(1, false))) + electron_sources = [get_settings_electrons(get_settings_ions(1, false)) for i ∈ 1:counter] end # put all neutral sources into neutral_source_data struct vector @@ -408,6 +404,10 @@ function get_source_profile(profile_type, width, relative_minimum, coord) end end return profile + elseif profile_type == "wall_exp_decay" + x = coord.grid + return @. (1.0 - relative_minimum) * exp(-(x-x[1]) / width) + relative_minimum + + (1.0 - relative_minimum) * exp(-(x[end]-x) / width) + relative_minimum else error("Unrecognised source profile type '$profile_type'.") end @@ -572,7 +572,7 @@ function initialize_external_source_amplitude!(moments, external_source_settings neutral_source_settings = external_source_settings.neutral for index ∈ eachindex(neutral_source_settings) if neutral_source_settings[index].active - if neutral_source.source_type == "energy" + if neutral_source_settings[index].source_type == "energy" @loop_r_z ir iz begin moments.neutral.external_source_amplitude[iz,ir,index] = neutral_source_settings[index].source_strength * @@ -1087,7 +1087,7 @@ function external_neutral_source!(pdf, fvec, moments, neutral_source, index, vze end - if neutral_source_settings.source_type == "energy" + if neutral_source.source_type == "energy" # Take particles out of pdf so source does not change density @loop_sn_r_z_vzeta_vr_vz isn ir iz ivzeta ivr ivz begin pdf[iveta,ivr,ivz,iz,ir,isn] -= dt * source_amplitude[iz,ir] * From 9bfb6310a9c69685a363fe1d3997e7277ec1ae00 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Mon, 9 Sep 2024 11:09:32 +0100 Subject: [PATCH 15/23] Temporary fix for file_io.jl, where only the first source profile of each species has its information written out to the HDF5 file. Need to update this when possible, so that each source_amplitude can be written out for each source instead. --- moment_kinetics/src/file_io.jl | 67 ++++++++++++++++------------------ 1 file changed, 32 insertions(+), 35 deletions(-) diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl index 0951973a6..48b60dee2 100644 --- a/moment_kinetics/src/file_io.jl +++ b/moment_kinetics/src/file_io.jl @@ -1333,14 +1333,13 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, ion_source_settings = external_source_settings.ion if any(x -> x.active, ion_source_settings) - n = length(ion_source_settings) external_source_amplitude = create_dynamic_variable!( - dynamic, "external_source_amplitude", mk_float, z, r, n; + dynamic, "external_source_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external source for ions", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_density_amplitude", mk_float, z, r, n; + dynamic, "external_source_density_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external density source for ions", units="n_ref*c_ref/L_ref") else @@ -1348,7 +1347,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_upar external_source_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_momentum_amplitude", mk_float, z, r, n; + dynamic, "external_source_momentum_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external momentum source for ions", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1356,17 +1355,17 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_ppar external_source_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_pressure_amplitude", mk_float, z, r, n; + dynamic, "external_source_pressure_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external pressure source for ions", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else external_source_pressure_amplitude = nothing end - if ion_source_settings.PI_density_controller_I != 0.0 && - ion_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") - if ion_source_settings.source_type == "density_profile_control" + if any(x -> x.PI_density_controller_I != 0.0 && x.source_type ∈ + ("density_profile_control", "density_midpoint_control"), ion_source_settings) + if any(x -> x.source_type == "density_profile_control", ion_source_settings) external_source_controller_integral = create_dynamic_variable!( - dynamic, "external_source_controller_integral", mk_float, z, r, n; + dynamic, "external_source_controller_integral", mk_float, z, r; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for ions") else @@ -1546,21 +1545,20 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi electron_source_settings = external_source_settings.electron if any(x -> x.active, electron_source_settings) - n = length(electron_source_settings) external_source_electron_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_amplitude", mk_float, z, r, n; + dynamic, "external_source_electron_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external source for electrons", units="n_ref/c_ref^3*c_ref/L_ref") external_source_electron_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_density_amplitude", mk_float, z, r, n; + dynamic, "external_source_electron_density_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external density source for electrons", units="n_ref*c_ref/L_ref") external_source_electron_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r, n; + dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external momentum source for electrons", units="m_ref*n_ref*c_ref*c_ref/L_ref") external_source_electron_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r, n; + dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external pressure source for electrons", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1755,14 +1753,13 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo neutral_source_settings = external_source_settings.neutral if n_neutral_species > 0 && any(x -> x.active, neutral_source_settings) - n = length(neutral_source_settings) external_source_neutral_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_amplitude", mk_float, z, r, n; + dynamic, "external_source_neutral_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external source for neutrals", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_neutral_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_density_amplitude", mk_float, z, r, n; + dynamic, "external_source_neutral_density_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external density source for neutrals", units="n_ref*c_ref/L_ref") else @@ -1770,7 +1767,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_upar external_source_neutral_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r, n; + dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external momentum source for neutrals", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1778,17 +1775,17 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_ppar external_source_neutral_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r, n; + dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r; parallel_io=parallel_io, description="Amplitude of the external pressure source for neutrals", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else external_source_neutral_pressure_amplitude = nothing end - if neutral_source_settings.PI_density_controller_I != 0.0 && - neutral_source_settings.source_type ∈ ("density_profile_control", "density_midpoint_control") - if neutral_source_settings.source_type == "density_profile_control" + if any(x -> x.PI_density_controller_I != 0.0 && x.source_type ∈ + ("density_profile_control", "density_midpoint_control"), neutral_source_settings) + if any(x -> x.source_type == "density_profile_control", neutral_source_settings) external_source_neutral_controller_integral = create_dynamic_variable!( - dynamic, "external_source_neutral_controller_integral", mk_float, z, r, n; + dynamic, "external_source_neutral_controller_integral", mk_float, z, r; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for neutrals") else @@ -2515,21 +2512,21 @@ function write_ion_moments_data_to_binary(scratch, moments, n_ion_species, t_par end if io_moments.external_source_amplitude !== nothing append_to_dynamic_var(io_moments.external_source_amplitude, - moments.ion.external_source_amplitude, t_idx, + moments.ion.external_source_amplitude[:, :, 1], t_idx, parallel_io, z, r) if moments.evolve_density append_to_dynamic_var(io_moments.external_source_density_amplitude, - moments.ion.external_source_density_amplitude, + moments.ion.external_source_density_amplitude[:, :, 1], t_idx, parallel_io, z, r) end if moments.evolve_upar append_to_dynamic_var(io_moments.external_source_momentum_amplitude, - moments.ion.external_source_momentum_amplitude, + moments.ion.external_source_momentum_amplitude[:, :, 1], t_idx, parallel_io, z, r) end if moments.evolve_ppar append_to_dynamic_var(io_moments.external_source_pressure_amplitude, - moments.ion.external_source_pressure_amplitude, + moments.ion.external_source_pressure_amplitude[:, :, 1], t_idx, parallel_io, z, r) end end @@ -2613,16 +2610,16 @@ function write_electron_moments_data_to_binary(scratch, moments, t_params, elect t_idx, parallel_io, z, r) if io_moments.external_source_electron_amplitude !== nothing append_to_dynamic_var(io_moments.external_source_electron_amplitude, - moments.electron.external_source_amplitude, t_idx, + moments.electron.external_source_amplitude[:, :, 1], t_idx, parallel_io, z, r) append_to_dynamic_var(io_moments.external_source_electron_density_amplitude, - moments.electron.external_source_density_amplitude, + moments.electron.external_source_density_amplitude[:, :, 1], t_idx, parallel_io, z, r) append_to_dynamic_var(io_moments.external_source_electron_momentum_amplitude, - moments.electron.external_source_momentum_amplitude, + moments.electron.external_source_momentum_amplitude[:, :, 1], t_idx, parallel_io, z, r) append_to_dynamic_var(io_moments.external_source_electron_pressure_amplitude, - moments.electron.external_source_pressure_amplitude, + moments.electron.external_source_pressure_amplitude[:, :, 1], t_idx, parallel_io, z, r) end append_to_dynamic_var(io_moments.electron_constraints_A_coefficient, @@ -2725,21 +2722,21 @@ function write_neutral_moments_data_to_binary(scratch, moments, n_neutral_specie t_idx, parallel_io, z, r, n_neutral_species) if io_moments.external_source_neutral_amplitude !== nothing append_to_dynamic_var(io_moments.external_source_neutral_amplitude, - moments.neutral.external_source_amplitude, t_idx, + moments.neutral.external_source_amplitude[:, :, 1], t_idx, parallel_io, z, r) if moments.evolve_density append_to_dynamic_var(io_moments.external_source_neutral_density_amplitude, - moments.neutral.external_source_density_amplitude, + moments.neutral.external_source_density_amplitude[:, :, 1], t_idx, parallel_io, z, r) end if moments.evolve_upar append_to_dynamic_var(io_moments.external_source_neutral_momentum_amplitude, - moments.neutral.external_source_momentum_amplitude, + moments.neutral.external_source_momentum_amplitude[:, :, 1], t_idx, parallel_io, z, r) end if moments.evolve_ppar append_to_dynamic_var(io_moments.external_source_neutral_pressure_amplitude, - moments.neutral.external_source_pressure_amplitude, + moments.neutral.external_source_pressure_amplitude[:, :, 1], t_idx, parallel_io, z, r) end end From 25bdbfe4a59490929a8b62af458bf36d2adbba55 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Mon, 9 Sep 2024 14:47:20 +0100 Subject: [PATCH 16/23] Add collisionality_scan option to krook operator for ions --- moment_kinetics/src/krook_collisions.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/moment_kinetics/src/krook_collisions.jl b/moment_kinetics/src/krook_collisions.jl index 1880610e9..c113f0bca 100644 --- a/moment_kinetics/src/krook_collisions.jl +++ b/moment_kinetics/src/krook_collisions.jl @@ -41,6 +41,10 @@ function setup_krook_collisions_input(toml_input::Dict, reference_params) input_section["nuii0"] = nuii_krook_default input_section["nuee0"] = nuee_krook_default input_section["nuei0"] = nuei_krook_default + elseif frequency_option == "collisionality_scan" + input_section["nuii0"] *= nuii_krook_default + input_section["nuee0"] *= nuee_krook_default + input_section["nuei0"] *= nuei_krook_default elseif frequency_option == "manual" # use the frequency from the input file # do nothing @@ -74,7 +78,7 @@ function get_collision_frequency_ii(collisions, n, vth) colk = collisions.krook nuii0 = colk.nuii0 frequency_option = colk.frequency_option - if frequency_option == "reference_parameters" + if frequency_option ∈ ("reference_parameters", "collisionality_scan") return @. nuii0 * n * vth^(-3) elseif frequency_option == "manual" # Include 0.0*n so that the result gets promoted to an array if n is an array, From 8dd90254c183b7642e9c1b6b1792f8155471b0d8 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Tue, 10 Sep 2024 14:51:11 +0100 Subject: [PATCH 17/23] Allow NamedTuple 'coords' in create_dynamic_variable!() This change makes it simple to add an extra dimension for which there is not a `coordinate` struct, replacing the special-purpose arguments for n_ion_species, n_neutral_species and diagnostic_var_size. --- moment_kinetics/ext/file_io_netcdf.jl | 103 ++------- moment_kinetics/src/file_io.jl | 290 ++++++++++++-------------- moment_kinetics/src/file_io_hdf5.jl | 94 ++------- 3 files changed, 176 insertions(+), 311 deletions(-) diff --git a/moment_kinetics/ext/file_io_netcdf.jl b/moment_kinetics/ext/file_io_netcdf.jl index bb7a86f2b..fa22f99c3 100644 --- a/moment_kinetics/ext/file_io_netcdf.jl +++ b/moment_kinetics/ext/file_io_netcdf.jl @@ -97,9 +97,17 @@ end function write_single_value!(file_or_group::NCDataset, name, value::Union{Number, AbstractString, AbstractArray{T,N}}, - coords...; parallel_io, n_ion_species=nothing, - n_neutral_species=nothing, description=nothing, - units=nothing) where {T,N} + coords::Union{coordinate,NamedTuple}...; parallel_io, + description=nothing, units=nothing) where {T,N} + + if any(c.n < 0 for c ∈ coords) + error("Got a negative `n` in $coords") + end + if any(c.n == 0 for c ∈ coords) + # No data to write + return nothing + end + if description !== nothing || units !== nothing attributes = Dict{String, Any}() if description !== nothing @@ -112,11 +120,6 @@ function write_single_value!(file_or_group::NCDataset, name, attributes = () end - if n_ion_species !== nothing && n_neutral_species != nothing - error("Cannot have both ion-species and neutral species dimensions." * - "Got n_ion_species=$n_ion_species, n_neutral_species=$n_neutral_species") - end - if isa(value, Number) || isa(value, AbstractString) coords !== () && error("cannot pass coordinates with a scalar") @@ -134,26 +137,6 @@ function write_single_value!(file_or_group::NCDataset, name, maybe_create_netcdf_dim(file_or_group, c) end dims = Tuple(c.name for c in coords) - - if n_ion_species !== nothing - if n_ion_species < 0 - error("n_ion_species must be non-negative, got $n_ion_species") - elseif n_ion_species == 0 - # No data to write - return nothing - end - maybe_create_netcdf_dim(file_or_group, "ion_species", n_ion_species) - dims = tuple(dims..., "ion_species") - elseif n_neutral_species !== nothing - if n_neutral_species < 0 - error("n_neutral_species must be non-negative, got $n_neutral_species") - elseif n_neutral_species == 0 - # No data to write - return nothing - end - maybe_create_netcdf_dim(file_or_group, "neutral_species", n_neutral_species) - dims = tuple(dims..., "neutral_species") - end end if isa(value, Bool) # As a hack, write bools to NetCDF as Char, as NetCDF does not support bools (?), @@ -169,50 +152,20 @@ function write_single_value!(file_or_group::NCDataset, name, end function create_dynamic_variable!(file_or_group::NCDataset, name, type, - coords::coordinate...; parallel_io, - n_ion_species=nothing, n_neutral_species=nothing, - diagnostic_var_size=nothing, description=nothing, - units=nothing) - - if n_ion_species !== nothing && n_neutral_species !== nothing - error("Variable should not contain both ion and neutral species dimensions. " - * "Got n_ion_species=$n_ion_species and " - * "n_neutral_species=$n_neutral_species") - end - if diagnostic_var_size !== nothing && n_ion_species !== nothing - error("Diagnostic variable should not contain both ion species dimension. Got " - * "diagnostic_var_size=$diagnostic_var_size and " - * "n_ion_species=$n_ion_species") + coords::Union{coordinate,NamedTuple}...; parallel_io, + description=nothing, units=nothing) + + if any(c.n < 0 for c ∈ coords) + error("Got a negative `n` in $coords") end - if diagnostic_var_size !== nothing && n_neutral_species !== nothing - error("Diagnostic variable should not contain both neutral species dimension. " - * "Got diagnostic_var_size=$diagnostic_var_size and " - * "n_neutral_species=$n_neutral_species") + if any(c.n == 0 for c ∈ coords) + # No data to write + return nothing end # Create time dimension if necessary maybe_create_netcdf_dim(file_or_group, "time", Inf) - # Create species dimension if necessary - if n_ion_species !== nothing - if n_ion_species < 0 - error("n_ion_species must be non-negative, got $n_ion_species") - elseif n_ion_species == 0 - # No data to write - return nothing - end - maybe_create_netcdf_dim(file_or_group, "ion_species", n_ion_species) - end - if n_neutral_species !== nothing - if n_neutral_species < 0 - error("n_neutral_species must be non-negative, got $n_neutral_species") - elseif n_neutral_species == 0 - # No data to write - return nothing - end - maybe_create_netcdf_dim(file_or_group, "neutral_species", n_neutral_species) - end - # Create other dimensions if necessary for c ∈ coords maybe_create_netcdf_dim(file_or_group, c) @@ -221,23 +174,7 @@ function create_dynamic_variable!(file_or_group::NCDataset, name, type, # create the variable so it can be expanded indefinitely (up to the largest unsigned # integer in size) in the time dimension coord_dims = Tuple(c.name for c ∈ coords) - if diagnostic_var_size !== nothing - if isa(diagnostic_var_size, Number) - # Make diagnostic_var_size a Tuple - diagnostic_var_size = (diagnostic_var_size,) - end - for (i,dim_size) ∈ enumerate(diagnostic_var_size) - maybe_create_netcdf_dim(file_or_group, "$name$i", dim_size) - end - fixed_dims = Tuple("$name$i" for i ∈ 1:length(diagnostic_var_size)) - elseif n_ion_species !== nothing - fixed_dims = tuple(coord_dims..., "ion_species") - elseif n_neutral_species !== nothing - fixed_dims = tuple(coord_dims..., "neutral_species") - else - fixed_dims = coord_dims - end - dims = tuple(fixed_dims..., "time") + dims = tuple(coord_dims..., "time") # create the variable so it can be expanded indefinitely (up to the largest unsigned # integer in size) in the time dimension diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl index 5f6adab0d..0b058278f 100644 --- a/moment_kinetics/src/file_io.jl +++ b/moment_kinetics/src/file_io.jl @@ -459,8 +459,7 @@ function setup_electron_io(io_input, vpa, vperp, z, r, composition, collisions, if io_input.write_electron_error_diagnostics io_f_electron_loworder = create_dynamic_variable!(dynamic, "f_electron_loworder", mk_float, - vpa, vperp, z, r; - n_ion_species=composition.n_ion_species, + vpa, vperp, z, r, parallel_io=parallel_io, description="low-order approximation to electron distribution function, used to diagnose timestepping error") else @@ -469,8 +468,7 @@ function setup_electron_io(io_input, vpa, vperp, z, r, composition, collisions, if io_input.write_electron_steady_state_diagnostics io_f_electron_start_last_timestep = create_dynamic_variable!(dynamic, "f_electron_start_last_timestep", - mk_float, vpa, vperp, z, r; - n_ion_species=composition.n_ion_species, + mk_float, vpa, vperp, z, r, parallel_io=parallel_io, description="electron distribution function at the start of the last electron pseudo-timestep before output, used to measure steady state residual") else @@ -599,7 +597,10 @@ Get names of all variables function get_variable_keys end """ - write_single_value!(file_or_group, name, value; description=nothing, units=nothing) + write_single_value!(file_or_group, name, + data::Union{Number, AbstractString, AbstractArray{T,N}}, + coords::Union{coordinate,mk_int,NamedTuple}...; parallel_io, + description=nothing, units=nothing) where {T,N} Write a single variable to a file or group. If a description or units are passed, add as attributes of the variable. @@ -773,21 +774,24 @@ function write_boundary_distributions!(fid, boundary_distributions, parallel_io, @serial_region begin boundary_distributions_io = create_io_group(fid, "boundary_distributions") + n_ion_species_coord = (name="n_ion_species", n=composition.n_ion_species) + n_neutral_species_coord = (name="n_neutral_species", + n=composition.n_neutral_species) write_single_value!(boundary_distributions_io, "pdf_rboundary_ion_left", boundary_distributions.pdf_rboundary_ion[:,:,:,1,:], vpa, vperp, z, - parallel_io=parallel_io, n_ion_species=composition.n_ion_species, + n_ion_species_coord; parallel_io=parallel_io, description="Initial ion-particle pdf at left radial boundary") write_single_value!(boundary_distributions_io, "pdf_rboundary_ion_right", boundary_distributions.pdf_rboundary_ion[:,:,:,2,:], vpa, vperp, z, - parallel_io=parallel_io, n_ion_species=composition.n_ion_species, + n_ion_species_coord; parallel_io=parallel_io, description="Initial ion-particle pdf at right radial boundary") write_single_value!(boundary_distributions_io, "pdf_rboundary_neutral_left", boundary_distributions.pdf_rboundary_neutral[:,:,:,:,1,:], vz, vr, vzeta, z, - parallel_io=parallel_io, n_neutral_species=composition.n_neutral_species, + n_neutral_species_coord; parallel_io=parallel_io, description="Initial neutral-particle pdf at left radial boundary") write_single_value!(boundary_distributions_io, "pdf_rboundary_neutral_right", boundary_distributions.pdf_rboundary_neutral[:,:,:,:,2,:], vz, vr, vzeta, z, - parallel_io=parallel_io, n_neutral_species=composition.n_neutral_species, + n_neutral_species_coord; parallel_io=parallel_io, description="Initial neutral-particle pdf at right radial boundary") end return nothing @@ -967,21 +971,16 @@ function define_io_coordinate!(parent, coord, coord_name, description, parallel_ end """ - create_dynamic_variable!(file_or_group, name, type, coords::coordinate...; - n_ion_species=1, n_neutral_species=1, - diagnostic_var_size=nothing, description=nothing, - units=nothing) + create_dynamic_variable!(file_or_group, name, type, + coords::{coordinate,NamedTuple}...; + description=nothing, units=nothing) Create a time-evolving variable in `file_or_group` named `name` of type `type`. `coords` are the coordinates corresponding to the dimensions of the array, in the order of the -array dimensions. The species dimension does not have a `coordinate`, so the number of -species is passed as `nspecies`. A description and/or units can be added with the keyword -arguments. - -If a Tuple giving an array size is passed to `diagnostic_var_size`, a 'diagnostic' -variable is created - i.e. one that does not depend on the coordinates, so is assumed to -be the same on all processes and only needs to be written from the root process (for each -output file). +array dimensions - they may be either `coordinate` structs or `NamedTuple`s that contain +at least the fields `name`, `n`. The species dimension does not have a `coordinate`, so +the number of species is passed as `nspecies`. A description and/or units can be added +with the keyword arguments. """ function create_dynamic_variable! end @@ -1078,13 +1077,13 @@ function define_dynamic_moment_variables!(fid, n_ion_species, n_neutral_species, n_failure_vars = length(t_params.failure_caused_by) io_failure_caused_by = create_dynamic_variable!( - dynamic, "failure_caused_by", mk_int; diagnostic_var_size=n_failure_vars, + dynamic, "failure_caused_by", mk_int, (name="n_failure_vars", n=n_failure_vars); parallel_io=parallel_io, description="cumulative count of how many times each variable caused a " * "timestep failure for the run") n_limit_vars = length(t_params.limit_caused_by) io_limit_caused_by = create_dynamic_variable!( - dynamic, "limit_caused_by", mk_int; diagnostic_var_size=n_limit_vars, + dynamic, "limit_caused_by", mk_int, (name="n_limit_vars", n=n_limit_vars); parallel_io=parallel_io, description="cumulative count of how many times each factor limited the " * "timestep for the run") @@ -1203,17 +1202,17 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, evolve_ppar, write_error_diagnostics, write_steady_state_diagnostics) dynamic = get_group(fid, "dynamic_data") + n_ion_species_coord = (name="n_ion_species", n=n_ion_species) # io_density is the handle for the ion particle density - io_density = create_dynamic_variable!(dynamic, "density", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_density = create_dynamic_variable!(dynamic, "density", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species density", units="n_ref") if write_error_diagnostics io_density_loworder = - create_dynamic_variable!(dynamic, "density_loworder", mk_float, z, r; - n_ion_species=n_ion_species, parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "density_loworder", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="low-order approximation to ion species density, used to diagnose timestepping error", units="n_ref") else @@ -1221,8 +1220,8 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if write_steady_state_diagnostics io_density_start_last_timestep = - create_dynamic_variable!(dynamic, "density_start_last_timestep", mk_float, z, r; - n_ion_species=n_ion_species, parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "density_start_last_timestep", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species density at the start of the last timestep before output, used to measure steady state residual", units="n_ref") else @@ -1230,15 +1229,14 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end # io_upar is the handle for the ion parallel flow density - io_upar = create_dynamic_variable!(dynamic, "parallel_flow", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_upar = create_dynamic_variable!(dynamic, "parallel_flow", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species parallel flow", units="c_ref = sqrt(2*T_ref/mi)") if write_error_diagnostics io_upar_loworder = - create_dynamic_variable!(dynamic, "parallel_flow_loworder", mk_float, z, r; - n_ion_species=n_ion_species, parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "parallel_flow_loworder", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="low-order approximation to ion species parallel flow, used to diagnose timestepping error", units="c_ref = sqrt(2*T_ref/mi)") else @@ -1247,7 +1245,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, if write_steady_state_diagnostics io_upar_start_last_timestep = create_dynamic_variable!(dynamic, "parallel_flow_start_last_timestep", - mk_float, z, r; n_ion_species=n_ion_species, + mk_float, z, r, n_ion_species_coord; parallel_io=parallel_io, description="ion species parallel flow at the start of the last timestep before output, used to measure steady state residual", units="c_ref = sqrt(2*T_ref/mi)") @@ -1256,15 +1254,14 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end # io_ppar is the handle for the ion parallel pressure - io_ppar = create_dynamic_variable!(dynamic, "parallel_pressure", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_ppar = create_dynamic_variable!(dynamic, "parallel_pressure", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species parallel pressure", units="n_ref*T_ref") if write_error_diagnostics io_ppar_loworder = - create_dynamic_variable!(dynamic, "parallel_pressure_loworder", mk_float, z, r; - n_ion_species=n_ion_species, parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "parallel_pressure_loworder", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="low-order approximation to ion species parallel pressure, used to diagnose timestepping error", units="n_ref*T_ref") else @@ -1273,7 +1270,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, if write_steady_state_diagnostics io_ppar_start_last_timestep = create_dynamic_variable!(dynamic, "parallel_pressure_start_last_timestep", - mk_float, z, r; n_ion_species=n_ion_species, + mk_float, z, r, n_ion_species_coord; parallel_io=parallel_io, description="ion species parallel pressure at the start of the last timestep before output, used to measure steady state residual", units="n_ref*T_ref") @@ -1282,15 +1279,14 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end # 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, + io_pperp = create_dynamic_variable!(dynamic, "perpendicular_pressure", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species perpendicular pressure", units="n_ref*T_ref") if write_error_diagnostics io_pperp_loworder = - create_dynamic_variable!(dynamic, "perpendicular_pressure_loworder", mk_float, z, - r; n_ion_species=n_ion_species, parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "perpendicular_pressure_loworder", mk_float, + z, r, n_ion_species_coord; parallel_io=parallel_io, description="low-order approximation to ion species perpendicular pressure, used to diagnose timestepping error", units="n_ref*T_ref") else @@ -1299,7 +1295,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, if write_steady_state_diagnostics io_pperp_start_last_timestep = create_dynamic_variable!(dynamic, "perpendicular_pressure_start_last_timestep", - mk_float, z, r; n_ion_species=n_ion_species, + mk_float, z, r, n_ion_species_coord; parallel_io=parallel_io, description="ion species perpendicular pressure at the start of the last timestep before output, used to measure steady state residual", units="n_ref*T_ref") @@ -1308,23 +1304,20 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end # io_qpar is the handle for the ion parallel heat flux - io_qpar = create_dynamic_variable!(dynamic, "parallel_heat_flux", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_qpar = create_dynamic_variable!(dynamic, "parallel_heat_flux", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species parallel heat flux", units="n_ref*T_ref*c_ref") # io_vth is the handle for the ion thermal speed - io_vth = create_dynamic_variable!(dynamic, "thermal_speed", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_vth = create_dynamic_variable!(dynamic, "thermal_speed", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species thermal speed", units="c_ref") # io_dSdt is the handle for the entropy production (due to collisions) - io_dSdt = create_dynamic_variable!(dynamic, "entropy_production", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, + io_dSdt = create_dynamic_variable!(dynamic, "entropy_production", mk_float, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species entropy production", units="") @@ -1384,41 +1377,38 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, if parallel_io || z.irank == 0 # io_chodura_lower is the handle for the ion thermal speed - io_chodura_lower = create_dynamic_variable!(dynamic, "chodura_integral_lower", mk_float, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, - description="Generalised Chodura integral lower sheath entrance", - units="c_ref") + io_chodura_lower = create_dynamic_variable!(dynamic, "chodura_integral_lower", mk_float, r, + n_ion_species_coord; + parallel_io=parallel_io, + description="Generalised Chodura integral lower sheath entrance", + units="c_ref") else io_chodura_lower = nothing end if parallel_io || z.irank == z.nrank - 1 # io_chodura_upper is the handle for the ion thermal speed - io_chodura_upper = create_dynamic_variable!(dynamic, "chodura_integral_upper", mk_float, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, - description="Generalised Chodura integral upper sheath entrance", - units="c_ref") + io_chodura_upper = create_dynamic_variable!(dynamic, "chodura_integral_upper", mk_float, r, + n_ion_species_coord; + parallel_io=parallel_io, + description="Generalised Chodura integral upper sheath entrance", + units="c_ref") else io_chodura_upper = nothing end if evolve_density || evolve_upar || evolve_ppar ion_constraints_A_coefficient = - create_dynamic_variable!(dynamic, "ion_constraints_A_coefficient", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, - description="'A' coefficient enforcing density constraint for ions") + create_dynamic_variable!(dynamic, "ion_constraints_A_coefficient", mk_float, + z, r, n_ion_species_coord; parallel_io=parallel_io, + description="'A' coefficient enforcing density constraint for ions") ion_constraints_B_coefficient = - create_dynamic_variable!(dynamic, "ion_constraints_B_coefficient", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, - description="'B' coefficient enforcing flow constraint for ions") + create_dynamic_variable!(dynamic, "ion_constraints_B_coefficient", mk_float, + z, r, n_ion_species_coord; parallel_io=parallel_io, + description="'B' coefficient enforcing flow constraint for ions") ion_constraints_C_coefficient = - create_dynamic_variable!(dynamic, "ion_constraints_C_coefficient", mk_float, z, r; - n_ion_species=n_ion_species, - parallel_io=parallel_io, - description="'C' coefficient enforcing pressure constraint for ions") + create_dynamic_variable!(dynamic, "ion_constraints_C_coefficient", mk_float, + z, r, n_ion_species_coord; parallel_io=parallel_io, + description="'C' coefficient enforcing pressure constraint for ions") else ion_constraints_A_coefficient = nothing ion_constraints_B_coefficient = nothing @@ -1449,9 +1439,9 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi if !electron_only_io # io_density is the handle for the ion particle density io_electron_density = create_dynamic_variable!(dynamic, "electron_density", mk_float, z, r; - parallel_io=parallel_io, - description="electron species density", - units="n_ref") + parallel_io=parallel_io, + description="electron species density", + units="n_ref") if write_error_diagnostics io_electron_density_loworder = create_dynamic_variable!(dynamic, "electron_density_loworder", mk_float, z, r; @@ -1473,9 +1463,9 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi # io_electron_upar is the handle for the electron parallel flow density io_electron_upar = create_dynamic_variable!(dynamic, "electron_parallel_flow", mk_float, z, r; - parallel_io=parallel_io, - description="electron species parallel flow", - units="c_ref = sqrt(2*T_ref/mi)") + parallel_io=parallel_io, + description="electron species parallel flow", + units="c_ref = sqrt(2*T_ref/mi)") if write_error_diagnostics io_electron_upar_loworder = create_dynamic_variable!(dynamic, "electron_parallel_flow_loworder", mk_float, z, @@ -1530,15 +1520,15 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi # io_electron_qpar is the handle for the electron parallel heat flux io_electron_qpar = create_dynamic_variable!(dynamic, "electron_parallel_heat_flux", mk_float, z, r; - parallel_io=parallel_io, - description="electron species parallel heat flux", - units="n_ref*T_ref*c_ref") + parallel_io=parallel_io, + description="electron species parallel heat flux", + units="n_ref*T_ref*c_ref") # io_electron_vth is the handle for the electron thermal speed io_electron_vth = create_dynamic_variable!(dynamic, "electron_thermal_speed", mk_float, z, r; - parallel_io=parallel_io, - description="electron species thermal speed", - units="c_ref") + parallel_io=parallel_io, + description="electron species thermal speed", + units="c_ref") electron_source_settings = external_source_settings.electron if electron_source_settings.active @@ -1567,16 +1557,16 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi electron_constraints_A_coefficient = create_dynamic_variable!(dynamic, "electron_constraints_A_coefficient", mk_float, z, r; - parallel_io=parallel_io, - description="'A' coefficient enforcing density constraint for electrons") + parallel_io=parallel_io, + description="'A' coefficient enforcing density constraint for electrons") electron_constraints_B_coefficient = create_dynamic_variable!(dynamic, "electron_constraints_B_coefficient", mk_float, z, r; - parallel_io=parallel_io, - description="'B' coefficient enforcing flow constraint for electrons") + parallel_io=parallel_io, + description="'B' coefficient enforcing flow constraint for electrons") electron_constraints_C_coefficient = create_dynamic_variable!(dynamic, "electron_constraints_C_coefficient", mk_float, z, r; - parallel_io=parallel_io, - description="'C' coefficient enforcing pressure constraint for electrons") + parallel_io=parallel_io, + description="'C' coefficient enforcing pressure constraint for electrons") if electron_physics ∈ (kinetic_electrons, kinetic_electrons_with_temperature_equation) io_electron_step_counter = create_dynamic_variable!( @@ -1601,15 +1591,15 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi n_failure_vars = length(t_params.failure_caused_by) io_electron_failure_caused_by = create_dynamic_variable!( - dynamic, "electron_failure_caused_by", mk_int; - diagnostic_var_size=n_failure_vars, parallel_io=parallel_io, + dynamic, "electron_failure_caused_by", mk_int, + (name="n_failure_vars", n=n_failure_vars); parallel_io=parallel_io, description="cumulative count of how many times each variable caused an " * "electron pseudo-timestep failure for the run") n_limit_vars = length(t_params.limit_caused_by) io_electron_limit_caused_by = create_dynamic_variable!( - dynamic, "electron_limit_caused_by", mk_int; diagnostic_var_size=n_limit_vars, - parallel_io=parallel_io, + dynamic, "electron_limit_caused_by", mk_int, + (name="n_limit_vars", n=n_limit_vars); parallel_io=parallel_io, description="cumulative count of how many times each factor limited the " * "electron pseudo-timestep for the run") @@ -1652,18 +1642,18 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo evolve_ppar, write_error_diagnostics, write_steady_state_diagnostics) dynamic = get_group(fid, "dynamic_data") + n_neutral_species_coord = (name="n_neutral_species", n=n_neutral_species) # io_density_neutral is the handle for the neutral particle density - io_density_neutral = create_dynamic_variable!(dynamic, "density_neutral", mk_float, z, r; - n_neutral_species=n_neutral_species, + io_density_neutral = create_dynamic_variable!(dynamic, "density_neutral", mk_float, z, + r, n_neutral_species_coord; parallel_io=parallel_io, description="neutral species density", units="n_ref") if write_error_diagnostics io_density_neutral_loworder = - create_dynamic_variable!(dynamic, "density_neutral_loworder", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "density_neutral_loworder", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="low-order approximation to neutral species density, used to diagnose timestepping error", units="n_ref") else @@ -1671,9 +1661,8 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if write_steady_state_diagnostics io_density_neutral_start_last_timestep = - create_dynamic_variable!(dynamic, "density_neutral_start_last_timestep", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "density_neutral_start_last_timestep", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species density at the start of the last timestep before output, used to measure steady state residual", units="n_ref") else @@ -1681,16 +1670,15 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end # io_uz_neutral is the handle for the neutral z momentum density - io_uz_neutral = create_dynamic_variable!(dynamic, "uz_neutral", mk_float, z, r; - n_neutral_species=n_neutral_species, + io_uz_neutral = create_dynamic_variable!(dynamic, "uz_neutral", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species mean z velocity", units="c_ref = sqrt(2*T_ref/mi)") if write_error_diagnostics io_uz_neutral_loworder = - create_dynamic_variable!(dynamic, "uz_neutral_loworder", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "uz_neutral_loworder", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="low-order approximation to neutral species mean z velocity, used to diagnose timestepping error", units="c_ref = sqrt(2*T_ref/mi)") else @@ -1698,8 +1686,8 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if write_steady_state_diagnostics io_uz_neutral_start_last_timestep = - create_dynamic_variable!(dynamic, "uz_neutral_start_last_timestep", mk_float, z, r; - n_neutral_species=n_neutral_species, + create_dynamic_variable!(dynamic, "uz_neutral_start_last_timestep", mk_float, + z, r, n_neutral_species_coord; parallel_io=parallel_io, description="neutral species mean z velocity at the start of the last timestep before output, used to measure steady state residual", units="c_ref = sqrt(2*T_ref/mi)") @@ -1708,16 +1696,15 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end # io_pz_neutral is the handle for the neutral species zz pressure - io_pz_neutral = create_dynamic_variable!(dynamic, "pz_neutral", mk_float, z, r; - n_neutral_species=n_neutral_species, + io_pz_neutral = create_dynamic_variable!(dynamic, "pz_neutral", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species mean zz pressure", units="n_ref*T_ref") if write_error_diagnostics io_pz_neutral_loworder = - create_dynamic_variable!(dynamic, "pz_neutral_loworder", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, + create_dynamic_variable!(dynamic, "pz_neutral_loworder", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="low-order approximation to neutral species mean zz pressure, used to diagnose timestepping error", units="n_ref*T_ref") else @@ -1725,8 +1712,8 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if write_steady_state_diagnostics io_pz_neutral_start_last_timestep = - create_dynamic_variable!(dynamic, "pz_neutral_start_last_timestep", mk_float, z, r; - n_neutral_species=n_neutral_species, + create_dynamic_variable!(dynamic, "pz_neutral_start_last_timestep", mk_float, + z, r, n_neutral_species_coord; parallel_io=parallel_io, description="neutral species mean zz pressure at the start of the last timestep before output, used to measure steady state residual", units="n_ref*T_ref") @@ -1735,16 +1722,15 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end # io_qz_neutral is the handle for the neutral z heat flux - io_qz_neutral = create_dynamic_variable!(dynamic, "qz_neutral", mk_float, z, r; - n_neutral_species=n_neutral_species, + io_qz_neutral = create_dynamic_variable!(dynamic, "qz_neutral", mk_float, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species z heat flux", units="n_ref*T_ref*c_ref") # io_thermal_speed_neutral is the handle for the neutral thermal speed io_thermal_speed_neutral = create_dynamic_variable!( - dynamic, "thermal_speed_neutral", mk_float, z, r; - n_neutral_species=n_neutral_species, + dynamic, "thermal_speed_neutral", mk_float, z, r, n_neutral_species_coord; parallel_io=parallel_io, description="neutral species thermal speed", units="c_ref") @@ -1804,20 +1790,20 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo if evolve_density || evolve_upar || evolve_ppar neutral_constraints_A_coefficient = - create_dynamic_variable!(dynamic, "neutral_constraints_A_coefficient", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, - description="'A' coefficient enforcing density constraint for neutrals") + create_dynamic_variable!(dynamic, "neutral_constraints_A_coefficient", + mk_float, z, r, n_neutral_species_coord; + parallel_io=parallel_io, + description="'A' coefficient enforcing density constraint for neutrals") neutral_constraints_B_coefficient = - create_dynamic_variable!(dynamic, "neutral_constraints_B_coefficient", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, - description="'B' coefficient enforcing flow constraint for neutrals") + create_dynamic_variable!(dynamic, "neutral_constraints_B_coefficient", + mk_float, z, r, n_neutral_species_coord; + parallel_io=parallel_io, + description="'B' coefficient enforcing flow constraint for neutrals") neutral_constraints_C_coefficient = - create_dynamic_variable!(dynamic, "neutral_constraints_C_coefficient", mk_float, z, r; - n_neutral_species=n_neutral_species, - parallel_io=parallel_io, - description="'C' coefficient enforcing pressure constraint for neutrals") + create_dynamic_variable!(dynamic, "neutral_constraints_C_coefficient", + mk_float, z, r, n_neutral_species_coord; + parallel_io=parallel_io, + description="'C' coefficient enforcing pressure constraint for neutrals") else neutral_constraints_A_coefficient = nothing neutral_constraints_B_coefficient = nothing @@ -1856,16 +1842,15 @@ function define_dynamic_dfn_variables!(fid, r, z, vperp, vpa, vzeta, vr, vz, com nl_solver_params) dynamic = get_group(fid, "dynamic_data") + n_ion_species_coord = (name="n_ion_species", n=composition.n_ion_species) # io_f is the handle for the ion pdf - io_f = create_dynamic_variable!(dynamic, "f", mk_float, vpa, vperp, z, r; - n_ion_species=composition.n_ion_species, - parallel_io=parallel_io, + io_f = create_dynamic_variable!(dynamic, "f", mk_float, vpa, vperp, z, r, + n_ion_species_coord; parallel_io=parallel_io, description="ion species distribution function") if io_input.write_error_diagnostics io_f_loworder = create_dynamic_variable!(dynamic, "f_loworder", mk_float, vpa, - vperp, z, r; - n_ion_species=composition.n_ion_species, + vperp, z, r, n_ion_species_coord; parallel_io=parallel_io, description="low-order approximation to ion species distribution function, used to diagnose timestepping error") else @@ -1874,8 +1859,7 @@ function define_dynamic_dfn_variables!(fid, r, z, vperp, vpa, vzeta, vr, vz, com if io_input.write_steady_state_diagnostics io_f_start_last_timestep = create_dynamic_variable!(dynamic, "f_start_last_timestep", mk_float, vpa, - vperp, z, r; - n_ion_species=composition.n_ion_species, + vperp, z, r, n_ion_species_coord; parallel_io=parallel_io, description="ion species distribution function at the start of the last timestep before output, used to measure steady state residual") else @@ -1913,16 +1897,17 @@ function define_dynamic_dfn_variables!(fid, r, z, vperp, vpa, vzeta, vr, vz, com io_f_electron_start_last_timestep = nothing end + n_neutral_species_coord = (name="n_neutral_species", n=composition.n_neutral_species) + # io_f_neutral is the handle for the neutral pdf - io_f_neutral = create_dynamic_variable!(dynamic, "f_neutral", mk_float, vz, vr, vzeta, z, r; - n_neutral_species=composition.n_neutral_species, + io_f_neutral = create_dynamic_variable!(dynamic, "f_neutral", mk_float, vz, vr, vzeta, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species distribution function") if io_input.write_error_diagnostics io_f_neutral_loworder = create_dynamic_variable!(dynamic, "f_neutral_loworder", mk_float, vz, vr, - vzeta, z, r; - n_ion_species=composition.n_ion_species, + vzeta, z, r, n_neutral_species_coord; parallel_io=parallel_io, description="low-order approximation to neutral species distribution function, used to diagnose timestepping error") else @@ -1931,9 +1916,8 @@ function define_dynamic_dfn_variables!(fid, r, z, vperp, vpa, vzeta, vr, vz, com if io_input.write_steady_state_diagnostics io_f_neutral_start_last_timestep = create_dynamic_variable!(dynamic, "f_neutral_start_last_timestep", - mk_float, vz, vr, vzeta, z, r; - n_ion_species=composition.n_ion_species, - parallel_io=parallel_io, + mk_float, vz, vr, vzeta, z, r, + n_neutral_species_coord; parallel_io=parallel_io, description="neutral species distribution function at the start of the last timestep before output, used to measure steady state residual") else io_f_neutral_start_last_timestep = nothing diff --git a/moment_kinetics/src/file_io_hdf5.jl b/moment_kinetics/src/file_io_hdf5.jl index 271d79fee..53675051c 100644 --- a/moment_kinetics/src/file_io_hdf5.jl +++ b/moment_kinetics/src/file_io_hdf5.jl @@ -85,8 +85,7 @@ end # HDF5.H5DataStore is the supertype for HDF5.File and HDF5.Group function write_single_value!(file_or_group::HDF5.H5DataStore, name, data::Union{Number, AbstractString, AbstractArray{T,N}}, - coords::Union{coordinate,mk_int}...; parallel_io, - n_ion_species=nothing, n_neutral_species=nothing, + coords::Union{coordinate,mk_int,NamedTuple}...; parallel_io, description=nothing, units=nothing) where {T,N} if isa(data, Union{Number, AbstractString}) file_or_group[name] = data @@ -99,32 +98,18 @@ function write_single_value!(file_or_group::HDF5.H5DataStore, name, return nothing end - if n_ion_species !== nothing && n_neutral_species != nothing - error("Cannot have both ion-species and neutral species dimensions." * - "Got n_ion_species=$n_ion_species, n_neutral_species=$n_neutral_species") + if any(isa(c, mk_int) ? c < 0 : c.n < 0 for c ∈ coords) + error("Got a negative `n` in $coords") end - - if n_ion_species !== nothing - if n_ion_species < 0 - error("n_ion_species must be non-negative, got $n_ion_species") - elseif n_ion_species == 0 - # No data to write - return nothing - end - coords = tuple(coords..., n_ion_species) - elseif n_neutral_species !== nothing - if n_neutral_species < 0 - error("n_neutral_species must be non-negative, got $n_neutral_species") - elseif n_neutral_species == 0 - # No data to write - return nothing - end - coords = tuple(coords..., n_neutral_species) + if any(isa(c, mk_int) ? c == 0 : c.n == 0 for c ∈ coords) + # No data to write + return nothing end + dim_sizes, chunk_sizes = hdf5_get_fixed_dim_sizes(coords, parallel_io) io_var = create_dataset(file_or_group, name, T, dim_sizes, chunk=chunk_sizes) - local_ranges = Tuple(isa(c, mk_int) ? (1:c) : c.local_io_range for c ∈ coords) - global_ranges = Tuple(isa(c, mk_int) ? (1:c) : c.global_io_range for c ∈ coords) + local_ranges = Tuple(isa(c, mk_int) ? (1:c) : isa(c, coordinate) ? c.local_io_range : c.n for c ∈ coords) + global_ranges = Tuple(isa(c, mk_int) ? (1:c) : isa(c, coordinate) ? c.global_io_range : c.n for c ∈ coords) if N == 1 io_var[global_ranges[1]] = @view data[local_ranges[1]] @@ -176,7 +161,7 @@ of species). """ function hdf5_get_fixed_dim_sizes(coords, parallel_io) if parallel_io - dim_sizes = Tuple(isa(c, mk_int) ? c : c.n_global for c in coords) + dim_sizes = Tuple(isa(c, mk_int) ? c : (isa(c, coordinate) ? c.n_global : c.n) for c in coords) else dim_sizes = Tuple(isa(c, mk_int) ? c : c.n for c in coords) end @@ -209,65 +194,24 @@ function hdf5_get_dynamic_dim_sizes(fixed_coords, parallel_io) end function create_dynamic_variable!(file_or_group::HDF5.H5DataStore, name, type, - coords::coordinate...; parallel_io, - n_ion_species=nothing, n_neutral_species=nothing, - diagnostic_var_size=nothing, description=nothing, - units=nothing) - - if n_ion_species !== nothing && n_neutral_species !== nothing - error("Variable should not contain both ion and neutral species dimensions. " - * "Got n_ion_species=$n_ion_species and " - * "n_neutral_species=$n_neutral_species") - end - if diagnostic_var_size !== nothing && n_ion_species !== nothing - error("Diagnostic variable should not contain both ion species dimension. Got " - * "diagnostic_var_size=$diagnostic_var_size and " - * "n_ion_species=$n_ion_species") + coords::Union{coordinate,NamedTuple}...; parallel_io, + description=nothing, units=nothing) + + if any(isa(c, mk_int) ? c < 0 : c.n < 0 for c ∈ coords) + error("Got a negative `n` in $coords") end - if diagnostic_var_size !== nothing && n_neutral_species !== nothing - error("Diagnostic variable should not contain both neutral species dimension. " - * "Got diagnostic_var_size=$diagnostic_var_size and " - * "n_neutral_species=$n_neutral_species") + if any(isa(c, mk_int) ? c == 0 : c.n == 0 for c ∈ coords) + # No data to write + return nothing end - # Add the number of species to the spatial/velocity-space coordinates - if diagnostic_var_size !== nothing - if isa(diagnostic_var_size, Number) - # Make diagnostic_var_size a Tuple - diagnostic_var_size = (diagnostic_var_size,) - end - fixed_coords = diagnostic_var_size - elseif n_ion_species !== nothing - if n_ion_species < 0 - error("n_ion_species must be non-negative, got $n_ion_species") - elseif n_ion_species == 0 - # No data to write - return nothing - end - fixed_coords = tuple(coords..., n_ion_species) - elseif n_neutral_species !== nothing - if n_neutral_species < 0 - error("n_neutral_species must be non-negative, got $n_neutral_species") - elseif n_neutral_species == 0 - # No data to write - return nothing - end - fixed_coords = tuple(coords..., n_neutral_species) - else - fixed_coords = coords - end initial_dim_sizes, max_dim_sizes, chunk_size = - hdf5_get_dynamic_dim_sizes(fixed_coords, parallel_io) + hdf5_get_dynamic_dim_sizes(coords, parallel_io) var = create_dataset(file_or_group, name, type, (initial_dim_sizes, max_dim_sizes), chunk=chunk_size) # Add attribute listing the dimensions belonging to this variable dim_names = Tuple(c.name for c ∈ coords) - if n_ion_species !== nothing - dim_names = tuple(dim_names..., "ion_species") - elseif n_neutral_species !== nothing - dim_names = tuple(dim_names..., "neutral_species") - end add_attribute!(var, "dims", join(dim_names, ",")) if description !== nothing From 8d61d1e895210725015835a0030bf4647f991728 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Wed, 11 Sep 2024 15:21:08 +0100 Subject: [PATCH 18/23] Allow output of all sources as one variable (after John's update to allow create_dynamic_variable() to take in not just coords but NamedTuples to help define dimensions of an array) --- moment_kinetics/src/file_io.jl | 81 ++++++++++++++++++---------------- 1 file changed, 43 insertions(+), 38 deletions(-) diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl index b81243b9c..614caede7 100644 --- a/moment_kinetics/src/file_io.jl +++ b/moment_kinetics/src/file_io.jl @@ -1323,13 +1323,14 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, ion_source_settings = external_source_settings.ion if any(x -> x.active, ion_source_settings) + n_sources = (name="n_ion_sources", n=length(ion_source_settings)) external_source_amplitude = create_dynamic_variable!( - dynamic, "external_source_amplitude", mk_float, z, r; + dynamic, "external_source_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external source for ions", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_density_amplitude", mk_float, z, r; + dynamic, "external_source_density_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external density source for ions", units="n_ref*c_ref/L_ref") else @@ -1337,7 +1338,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_upar external_source_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_momentum_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external momentum source for ions", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1345,7 +1346,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, end if evolve_ppar external_source_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_pressure_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external pressure source for ions", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1355,7 +1356,7 @@ function define_dynamic_ion_moment_variables!(fid, n_ion_species, r::coordinate, ("density_profile_control", "density_midpoint_control"), ion_source_settings) if any(x -> x.source_type == "density_profile_control", ion_source_settings) external_source_controller_integral = create_dynamic_variable!( - dynamic, "external_source_controller_integral", mk_float, z, r; + dynamic, "external_source_controller_integral", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for ions") else @@ -1532,20 +1533,21 @@ function define_dynamic_electron_moment_variables!(fid, r::coordinate, z::coordi electron_source_settings = external_source_settings.electron if any(x -> x.active, electron_source_settings) + n_sources = (name="n_electron_sources", n=length(electron_source_settings)) external_source_electron_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_amplitude", mk_float, z, r; + dynamic, "external_source_electron_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external source for electrons", units="n_ref/c_ref^3*c_ref/L_ref") external_source_electron_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_density_amplitude", mk_float, z, r; + dynamic, "external_source_electron_density_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external density source for electrons", units="n_ref*c_ref/L_ref") external_source_electron_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_electron_momentum_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external momentum source for electrons", units="m_ref*n_ref*c_ref*c_ref/L_ref") external_source_electron_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_electron_pressure_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external pressure source for electrons", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1736,13 +1738,14 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo neutral_source_settings = external_source_settings.neutral if n_neutral_species > 0 && any(x -> x.active, neutral_source_settings) + n_sources = (name="n_neutral_sources", n=length(neutral_source_settings)) external_source_neutral_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external source for neutrals", units="n_ref/c_ref^3*c_ref/L_ref") if evolve_density external_source_neutral_density_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_density_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_density_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external density source for neutrals", units="n_ref*c_ref/L_ref") else @@ -1750,7 +1753,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_upar external_source_neutral_momentum_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_momentum_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external momentum source for neutrals", units="m_ref*n_ref*c_ref*c_ref/L_ref") else @@ -1758,7 +1761,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo end if evolve_ppar external_source_neutral_pressure_amplitude = create_dynamic_variable!( - dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r; + dynamic, "external_source_neutral_pressure_amplitude", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Amplitude of the external pressure source for neutrals", units="m_ref*n_ref*c_ref^2*c_ref/L_ref") else @@ -1768,7 +1771,7 @@ function define_dynamic_neutral_moment_variables!(fid, n_neutral_species, r::coo ("density_profile_control", "density_midpoint_control"), neutral_source_settings) if any(x -> x.source_type == "density_profile_control", neutral_source_settings) external_source_neutral_controller_integral = create_dynamic_variable!( - dynamic, "external_source_neutral_controller_integral", mk_float, z, r; + dynamic, "external_source_neutral_controller_integral", mk_float, z, r, n_sources; parallel_io=parallel_io, description="Integral term for the PID controller of the external source for neutrals") else @@ -2492,34 +2495,35 @@ function write_ion_moments_data_to_binary(scratch, moments, n_ion_species, t_par parallel_io, 0, n_ion_species) end if io_moments.external_source_amplitude !== nothing + n_sources = size(moments.ion.external_source_amplitude)[3] append_to_dynamic_var(io_moments.external_source_amplitude, - moments.ion.external_source_amplitude[:, :, 1], t_idx, - parallel_io, z, r) + moments.ion.external_source_amplitude, t_idx, + parallel_io, z, r, n_sources) if moments.evolve_density append_to_dynamic_var(io_moments.external_source_density_amplitude, - moments.ion.external_source_density_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.ion.external_source_density_amplitude, + t_idx, parallel_io, z, r, n_sources) end if moments.evolve_upar append_to_dynamic_var(io_moments.external_source_momentum_amplitude, - moments.ion.external_source_momentum_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.ion.external_source_momentum_amplitude, + t_idx, parallel_io, z, r, n_sources) end if moments.evolve_ppar append_to_dynamic_var(io_moments.external_source_pressure_amplitude, - moments.ion.external_source_pressure_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.ion.external_source_pressure_amplitude, + t_idx, parallel_io, z, r, n_sources) end end if io_moments.external_source_controller_integral !== nothing - if size(moments.ion.external_source_controller_integral) == (1,1) + if size(moments.ion.external_source_controller_integral) == (1,1, n_sources) append_to_dynamic_var(io_moments.external_source_controller_integral, - moments.ion.external_source_controller_integral[1,1], + moments.ion.external_source_controller_integral, t_idx, parallel_io) else append_to_dynamic_var(io_moments.external_source_controller_integral, moments.ion.external_source_controller_integral, - t_idx, parallel_io, z, r) + t_idx, parallel_io, z, r, n_sources) end end if moments.evolve_density || moments.evolve_upar || moments.evolve_ppar @@ -2590,18 +2594,19 @@ function write_electron_moments_data_to_binary(scratch, moments, t_params, elect append_to_dynamic_var(io_moments.electron_thermal_speed, moments.electron.vth, t_idx, parallel_io, z, r) if io_moments.external_source_electron_amplitude !== nothing + n_sources = size(moments.electron.external_source_amplitude)[3] append_to_dynamic_var(io_moments.external_source_electron_amplitude, - moments.electron.external_source_amplitude[:, :, 1], t_idx, - parallel_io, z, r) + moments.electron.external_source_amplitude, t_idx, + parallel_io, z, r, n_sources) append_to_dynamic_var(io_moments.external_source_electron_density_amplitude, - moments.electron.external_source_density_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.electron.external_source_density_amplitude, + t_idx, parallel_io, z, r, n_sources) append_to_dynamic_var(io_moments.external_source_electron_momentum_amplitude, - moments.electron.external_source_momentum_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.electron.external_source_momentum_amplitude, + t_idx, parallel_io, z, r, n_sources) append_to_dynamic_var(io_moments.external_source_electron_pressure_amplitude, - moments.electron.external_source_pressure_amplitude[:, :, 1], - t_idx, parallel_io, z, r) + moments.electron.external_source_pressure_amplitude, + t_idx, parallel_io, z, r, n_sources) end append_to_dynamic_var(io_moments.electron_constraints_A_coefficient, moments.electron.constraints_A_coefficient, t_idx, @@ -2703,28 +2708,28 @@ function write_neutral_moments_data_to_binary(scratch, moments, n_neutral_specie t_idx, parallel_io, z, r, n_neutral_species) if io_moments.external_source_neutral_amplitude !== nothing append_to_dynamic_var(io_moments.external_source_neutral_amplitude, - moments.neutral.external_source_amplitude[:, :, 1], t_idx, + moments.neutral.external_source_amplitude, t_idx, parallel_io, z, r) if moments.evolve_density append_to_dynamic_var(io_moments.external_source_neutral_density_amplitude, - moments.neutral.external_source_density_amplitude[:, :, 1], + moments.neutral.external_source_density_amplitude, t_idx, parallel_io, z, r) end if moments.evolve_upar append_to_dynamic_var(io_moments.external_source_neutral_momentum_amplitude, - moments.neutral.external_source_momentum_amplitude[:, :, 1], + moments.neutral.external_source_momentum_amplitude, t_idx, parallel_io, z, r) end if moments.evolve_ppar append_to_dynamic_var(io_moments.external_source_neutral_pressure_amplitude, - moments.neutral.external_source_pressure_amplitude[:, :, 1], + moments.neutral.external_source_pressure_amplitude, t_idx, parallel_io, z, r) end end if io_moments.external_source_neutral_controller_integral !== nothing if size(moments.neutral.external_source_neutral_controller_integral) == (1,1) append_to_dynamic_var(io_moments.external_source_neutral_controller_integral, - moments.neutral.external_source_controller_integral[1,1], + moments.neutral.external_source_controller_integral, t_idx, parallel_io) else append_to_dynamic_var(io_moments.external_source_neutral_controller_integral, From d4dd56525edb57f18e7e7389979bd9361eab01b5 Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Wed, 11 Sep 2024 16:23:13 +0100 Subject: [PATCH 19/23] bug fixes for plotting sources in makie_post_processing --- .../src/makie_post_processing.jl | 27 ++++++++++--------- moment_kinetics/src/load_data.jl | 14 ++++++++-- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl index fb5aeefa0..65a9ab301 100644 --- a/makie_post_processing/makie_post_processing/src/makie_post_processing.jl +++ b/makie_post_processing/makie_post_processing/src/makie_post_processing.jl @@ -40,7 +40,8 @@ using moment_kinetics.load_data: close_run_info, get_run_info_no_setup, get_vari neutral_moment_variables, all_moment_variables, ion_dfn_variables, electron_dfn_variables, neutral_dfn_variables, all_dfn_variables, ion_variables, - neutral_variables, all_variables + neutral_variables, all_variables, ion_source_variables, + neutral_source_variables, electron_source_variables using moment_kinetics.initial_conditions: vpagrid_to_dzdt using .shared_utils: calculate_and_write_frequencies using moment_kinetics.type_definitions: mk_float, mk_int @@ -233,14 +234,6 @@ function makie_post_process(run_dir::Union{String,Tuple}, # Plots from moment variables ############################# - moment_variable_list = tuple(em_variables..., ion_moment_variables...) - if has_electrons - moment_variable_list = tuple(moment_variable_list..., electron_moment_variables...) - end - if has_neutrals - moment_variable_list = tuple(moment_variable_list..., neutral_moment_variables...) - end - if any(ri !== nothing for ri ∈ run_info_moments) has_moments = true @@ -276,7 +269,7 @@ function makie_post_process(run_dir::Union{String,Tuple}, end do_steady_state_residuals = any(input_dict[v]["steady_state_residual"] - for v ∈ moment_variable_list) + for v ∈ all_moment_variables) if do_steady_state_residuals textoutput_files = Tuple(ri.run_prefix * "_residuals.txt" for ri in run_info if ri !== nothing) @@ -299,7 +292,7 @@ function makie_post_process(run_dir::Union{String,Tuple}, steady_state_residual_fig_axes = nothing end - for variable_name ∈ moment_variable_list + for variable_name ∈ all_moment_variables plots_for_variable(run_info, variable_name; plot_prefix=plot_prefix, has_rdim=has_rdim, has_zdim=has_zdim, is_1V=is_1V, steady_state_residual_fig_axes=steady_state_residual_fig_axes) @@ -1030,8 +1023,18 @@ function plots_for_variable(run_info, variable_name; plot_prefix, has_rdim=true, elseif variable_name ∈ neutral_moment_variables || variable_name ∈ neutral_dfn_variables species_indices = 1:maximum(ri.n_neutral_species for ri ∈ run_info) - else + elseif variable_name ∈ ion_moment_variables || + variable_name ∈ ion_dfn_variables species_indices = 1:maximum(ri.n_ion_species for ri ∈ run_info) + elseif variable_name in ion_source_variables + species_indices = 1:maximum(length(ri.external_source_settings.ion) for ri ∈ run_info) + elseif variable_name in electron_source_variables + species_indices = 1:maximum(length(ri.external_source_settings.electron) for ri ∈ run_info) + elseif variable_name in neutral_source_variables + species_indices = 1:maximum(length(ri.external_source_settings.neutral) for ri ∈ run_info) + else + species_indices = 1:1 + #error("variable_name=$variable_name not found in any defined group") end for is ∈ species_indices if is !== nothing diff --git a/moment_kinetics/src/load_data.jl b/moment_kinetics/src/load_data.jl index 99c260372..cf9f56883 100644 --- a/moment_kinetics/src/load_data.jl +++ b/moment_kinetics/src/load_data.jl @@ -66,10 +66,20 @@ const neutral_moment_variables = ("density_neutral", "uz_neutral", "pz_neutral", "thermal_speed_neutral", "temperature_neutral", "qz_neutral", "total_energy_neutral", "total_energy_flux_neutral") +const ion_source_variables = ("external_source_amplitude", "external_source_density_amplitude", + "external_source_momentum_amplitude", "external_source_pressure_amplitude", + "external_source_controller_integral") +const neutral_source_variables = ("external_source_neutral_amplitude", "external_source_neutral_density_amplitude", + "external_source_neutral_momentum_amplitude", "external_source_neutral_pressure_amplitude", + "external_source_neutral_controller_integral") +const electron_source_variables = ("external_source_electron_amplitude", "external_source_electron_density_amplitude", + "external_source_electron_momentum_amplitude", "external_source_electron_pressure_amplitude") const all_moment_variables = tuple(em_variables..., ion_moment_variables..., electron_moment_variables..., - neutral_moment_variables...) - + neutral_moment_variables..., + ion_source_variables..., + electron_source_variables..., + neutral_source_variables...) const ion_dfn_variables = ("f",) const electron_dfn_variables = ("f_electron",) const neutral_dfn_variables = ("f_neutral",) From e8f29265a200c8cea6efa13783bac2f7387e5011 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 12 Sep 2024 20:33:15 +0100 Subject: [PATCH 20/23] Use ion_sources directly instead of `counter` to set up electron sources `counter` is not the same as `length(ion_sources)` when one or more ion source sections are found. It is simpler to directly iterate over `ion_sources`. --- moment_kinetics/src/external_sources.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index b53e8b71e..0aaa408e9 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -358,9 +358,9 @@ function setup_external_sources!(input_dict, r, z, electron_physics) electron_sources = electron_source_data[] if electron_physics ∈ (braginskii_fluid, kinetic_electrons, kinetic_electrons_with_temperature_equation) - electron_sources = [get_settings_electrons(ion_sources[i]) for i ∈ 1:counter] + electron_sources = [get_settings_electrons(this_source) for this_source ∈ ion_sources] else - electron_sources = [get_settings_electrons(get_settings_ions(1, false)) for i ∈ 1:counter] + electron_sources = [get_settings_electrons(get_settings_ions(1, false))] end # put all neutral sources into neutral_source_data struct vector From 150d05a039761acb422b127eee8462271d2efe25 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 12 Sep 2024 21:53:08 +0100 Subject: [PATCH 21/23] Fix import of external source function in electron_kinetic_equation --- moment_kinetics/src/electron_kinetic_equation.jl | 2 +- moment_kinetics/src/external_sources.jl | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/electron_kinetic_equation.jl b/moment_kinetics/src/electron_kinetic_equation.jl index ce89fdc81..78772beed 100644 --- a/moment_kinetics/src/electron_kinetic_equation.jl +++ b/moment_kinetics/src/electron_kinetic_equation.jl @@ -23,7 +23,7 @@ using ..electron_fluid_equations: electron_energy_equation!, electron_energy_res using ..electron_z_advection: electron_z_advection!, update_electron_speed_z! using ..electron_vpa_advection: electron_vpa_advection!, update_electron_speed_vpa! using ..em_fields: update_phi! -using ..external_sources: external_electron_source! +using ..external_sources: total_external_electron_sources! using ..file_io: get_electron_io_info, write_electron_state, finish_electron_io using ..krook_collisions: electron_krook_collisions! using ..moment_constraints: hard_force_moment_constraints!, diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 0aaa408e9..07474ca76 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -16,7 +16,8 @@ export setup_external_sources!, external_ion_source!, external_neutral_source!, initialize_external_source_amplitude!, initialize_external_source_controller_integral!, total_external_ion_sources!, total_external_neutral_sources!, - total_external_ion_source_controllers!, total_external_neutral_source_controllers! + total_external_ion_source_controllers!, total_external_neutral_source_controllers!, + external_electron_source!, total_external_electron_sources! using ..array_allocation: allocate_float, allocate_shared_float using ..calculus From fe922fa9df9290cd6dff7529883ed78f582e5d8e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 12 Sep 2024 22:10:03 +0100 Subject: [PATCH 22/23] Fix typo in external_electron_source!() --- moment_kinetics/src/external_sources.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 07474ca76..444833d11 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -935,7 +935,8 @@ function total_external_electron_sources!(pdf_out, pdf_in, electron_density, ele end """ - external_electron_source!(pdf, fvec, moments, electron_source_settings, vperp, + external_electron_source!(pdf_out, pdf_in, electron_density, electron_upar, + moments, composition, electron_source, index, vperp, vpa, dt) Add external source term to the electron kinetic equation. @@ -974,7 +975,7 @@ function external_electron_source!(pdf_out, pdf_in, electron_density, electron_u end end - if electron_source_settings.source_type == "energy" + if electron_source.source_type == "energy" # Take particles out of pdf so source does not change density @loop_r_z_vperp_vpa ir iz ivperp ivpa begin pdf_out[ivpa,ivperp,iz,ir] -= dt * source_amplitude[iz,ir] * From 5ecabac4e89173f3fa76e28a9250f76e487080c8 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Thu, 12 Sep 2024 22:24:18 +0100 Subject: [PATCH 23/23] Fix output of neutral source terms --- moment_kinetics/src/file_io.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/moment_kinetics/src/file_io.jl b/moment_kinetics/src/file_io.jl index 614caede7..3b02b6b7e 100644 --- a/moment_kinetics/src/file_io.jl +++ b/moment_kinetics/src/file_io.jl @@ -2707,23 +2707,24 @@ function write_neutral_moments_data_to_binary(scratch, moments, n_neutral_specie append_to_dynamic_var(io_moments.thermal_speed_neutral, moments.neutral.vth, t_idx, parallel_io, z, r, n_neutral_species) if io_moments.external_source_neutral_amplitude !== nothing + n_sources = size(moments.neutral.external_source_amplitude)[3] append_to_dynamic_var(io_moments.external_source_neutral_amplitude, moments.neutral.external_source_amplitude, t_idx, - parallel_io, z, r) + parallel_io, z, r, n_sources) if moments.evolve_density append_to_dynamic_var(io_moments.external_source_neutral_density_amplitude, moments.neutral.external_source_density_amplitude, - t_idx, parallel_io, z, r) + t_idx, parallel_io, z, r, n_sources) end if moments.evolve_upar append_to_dynamic_var(io_moments.external_source_neutral_momentum_amplitude, moments.neutral.external_source_momentum_amplitude, - t_idx, parallel_io, z, r) + t_idx, parallel_io, z, r, n_sources) end if moments.evolve_ppar append_to_dynamic_var(io_moments.external_source_neutral_pressure_amplitude, moments.neutral.external_source_pressure_amplitude, - t_idx, parallel_io, z, r) + t_idx, parallel_io, z, r, n_sources) end end if io_moments.external_source_neutral_controller_integral !== nothing