From e60c61db3867474cd0eaa71fbb441a03fbfa193b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 15 Sep 2024 20:55:11 +0100 Subject: [PATCH 1/7] Set electron temperature to constant when using Boltzmann response Previously the electron temperature (and corresponding thermal speed) variables were unused, but set equal to the ion temperature, which has the potential to be confusing. --- moment_kinetics/src/initial_conditions.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/moment_kinetics/src/initial_conditions.jl b/moment_kinetics/src/initial_conditions.jl index 2549e07f4..ea9fc5f92 100644 --- a/moment_kinetics/src/initial_conditions.jl +++ b/moment_kinetics/src/initial_conditions.jl @@ -281,7 +281,7 @@ function initialize_electrons!(pdf, moments, fields, geometry, composition, r, z # Not restarting, so create initial profiles # initialise the electron thermal speed profile - init_electron_vth!(moments.electron.vth, moments.ion.vth, composition.T_e, composition.me_over_mi, z.grid) + init_electron_vth!(moments.electron.vth, moments.ion.vth, composition, z.grid) begin_r_z_region() # calculate the electron temperature from the thermal speed @loop_r_z ir iz begin @@ -1063,14 +1063,18 @@ initialise the electron thermal speed profile. for now the only initialisation option for the temperature is constant in z. returns vth0 = sqrt(2*Ts/Te) """ -function init_electron_vth!(vth_e, vth_i, T_e, me_over_mi, z) +function init_electron_vth!(vth_e, vth_i, composition, z) begin_r_z_region() - # @loop_r_z ir iz begin - # vth_e[iz,ir] = sqrt(T_e) - # end - @loop_r_z ir iz begin - vth_e[iz,ir] = vth_i[iz,ir,1] / sqrt(me_over_mi) - #vth_e[iz,ir] = exp(-5*(z[iz]/z[end])^2)/sqrt(me_over_mi) + if composition.electron_physics ∈ (boltzmann_electron_response, + boltzmann_electron_response_with_simple_sheath) + @loop_r_z ir iz begin + vth_e[iz,ir] = sqrt(composition.T_e / composition.me_over_mi) + end + else + @loop_r_z ir iz begin + vth_e[iz,ir] = vth_i[iz,ir,1] / sqrt(composition.me_over_mi) + #vth_e[iz,ir] = exp(-5*(z[iz]/z[end])^2)/sqrt(composition.me_over_mi) + end end end From bfe6211e313d29cbc13911850c4185d6b41e9b26 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sun, 15 Sep 2024 20:56:20 +0100 Subject: [PATCH 2/7] Set electron parallel flow equal to the ion one for Boltzmann response The electron parallel flow is not used when Boltzmann response is used, but most sensible thing is probably to set zero current. This version assumes a single ion species, and sets the electron parallel flow equal to the parallel flow of ion species 1. --- moment_kinetics/src/electron_fluid_equations.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/moment_kinetics/src/electron_fluid_equations.jl b/moment_kinetics/src/electron_fluid_equations.jl index 1ddc44046..1b92038c5 100644 --- a/moment_kinetics/src/electron_fluid_equations.jl +++ b/moment_kinetics/src/electron_fluid_equations.jl @@ -111,6 +111,11 @@ function calculate_electron_upar_from_charge_conservation!(upar_e, updated, dens # convert from parallel particle flux to parallel particle density upar_e[iz,ir] /= dens_e[iz,ir] end + else + begin_r_z_region() + @loop_r_z ir iz begin + upar_e[iz,ir] = upar_i[iz,ir,1] + end end updated[] = true end From afb6ba95826f6152ea5ce5a1ce479d93c49ca79e Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 16 Sep 2024 14:39:44 +0100 Subject: [PATCH 3/7] Skip electron moment plots for Boltzmann electron runs The electron moment plots are not helpful in Boltzmann electron runs, as they just replicate ion moments or a constant T_e, so skip those plots. --- .../makie_post_processing/src/makie_post_processing.jl | 4 ++++ 1 file changed, 4 insertions(+) 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 65a9ab301..00ccbe3fd 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 @@ -1004,6 +1004,10 @@ function plots_for_variable(run_info, variable_name; plot_prefix, has_rdim=true, all(ri.collisions.krook_collisions_option == "none" for ri ∈ run_info) # No Krook collisions active, so do not make plots. return nothing + elseif variable_name ∈ union(electron_moment_variables, electron_source_variables, electron_dfn_variables) && + all(ri.composition.electron_physics ∈ (boltzmann_electron_response, boltzmann_electron_response_with_simple_sheath) + for ri ∈ run_info) + return nothing end println("Making plots for $variable_name") From 9318084686d94dc7d058774d51487a4a5185b5bd Mon Sep 17 00:00:00 2001 From: Lucas McConnell Date: Mon, 16 Sep 2024 17:33:20 +0100 Subject: [PATCH 4/7] Add individual loop for electrons in external_sources initialisation, so that in the boltzmann case they can have a different number of sources than the ions. --- moment_kinetics/src/external_sources.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/moment_kinetics/src/external_sources.jl b/moment_kinetics/src/external_sources.jl index 444833d11..0c49c83b4 100644 --- a/moment_kinetics/src/external_sources.jl +++ b/moment_kinetics/src/external_sources.jl @@ -506,8 +506,10 @@ function initialize_external_source_amplitude!(moments, external_source_settings end end end + end - # now do same for electron sources, which (if present) are mostly mirrors of ion sources + # now do same for electron sources, which (if present) are mostly mirrors of ion sources + for index ∈ eachindex(electron_source_settings) if electron_source_settings[index].active if electron_source_settings[index].source_type == "energy" @loop_r_z ir iz begin From 29e3a9411ebd38972091aff4f842a48bc5e489c3 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Fri, 6 Sep 2024 20:19:20 +0100 Subject: [PATCH 5/7] Add some explanation of MPISharedArray type to the docs --- docs/src/developing.md | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/docs/src/developing.md b/docs/src/developing.md index 3ec2536de..ebe7b5f3c 100644 --- a/docs/src/developing.md +++ b/docs/src/developing.md @@ -28,6 +28,29 @@ julia> run_moment_kinetics("input.toml") It might be convenient to add `using Revise` to your `startup.jl` file (`~/julia/config/startup.jl`) so it's always loaded. +## Array types + +Most arrays in `moment_kinetics` are declared using a custom array type +`MPISharedArray`. Most of the time this type is just an alias for `Array`, and +so it needs the same template parameters (see [Julia's Array +documentation](https://docs.julialang.org/en/v1/manual/arrays/)) - the data +type and the number of dimensions, e.g. `MPISharedArray{mk_float,3}`. Although +these arrays use shared memory, Julia does not know about this. We use +`MPI.Win_allocate_shared()` to allocate the shared memory, then wrap it in an +`Array` in [`moment_kinetics.communication.allocate_shared`](@ref). + +The reason for using the alias, is that when the shared-memory debugging mode +is activated, we instead create arrays using a type `DebugMPISharedArray`, +which allows us to track some debugging information along with the array, see +[Shared memory debugging](@ref), and make `MPISharedArray` an alias for +`DebugMPISharedArray` instead. The reason for the alias is that if we declared +our structs with just `Array` type, then when debugging is activated we would +not be able to store `DebugMPISharedArray` instances in those structs, and if +we declared the structs with `AbstractArray`, they would not be concretely +typed, which could impact performance by creating code that is not 'type +stable' (i.e. all concrete types are known at compile time). + + ## Parallelization The code is parallelized at the moment using MPI and shared-memory arrays. Arrays representing the pdf, moments, etc. are shared between all processes. Using shared memory means, for example, we can take derivatives along one dimension while parallelising the other for any dimension without having to communicate to re-distribute the arrays. Using shared memory instead of (in future as well as) distributed memory parallelism has the advantage that it is easier to split up the points within each element between processors, giving a finer-grained parallelism which should let the code use larger numbers of processors efficiently. From f5949f8cfc3f3eb90d93c3c38faeb9f1b7992090 Mon Sep 17 00:00:00 2001 From: John Omotani Date: Sat, 7 Sep 2024 12:03:37 +0100 Subject: [PATCH 6/7] Docstring for MPISharedArray --- docs/src/developing.md | 5 +++-- moment_kinetics/src/communication.jl | 3 +++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/src/developing.md b/docs/src/developing.md index ebe7b5f3c..1ce3064f0 100644 --- a/docs/src/developing.md +++ b/docs/src/developing.md @@ -31,8 +31,9 @@ It might be convenient to add `using Revise` to your `startup.jl` file (`~/julia ## Array types Most arrays in `moment_kinetics` are declared using a custom array type -`MPISharedArray`. Most of the time this type is just an alias for `Array`, and -so it needs the same template parameters (see [Julia's Array +[`moment_kinetics.communication.MPISharedArray`](@ref). Most of the time this +type is just an alias for `Array`, and so it needs the same template parameters +(see [Julia's Array documentation](https://docs.julialang.org/en/v1/manual/arrays/)) - the data type and the number of dimensions, e.g. `MPISharedArray{mk_float,3}`. Although these arrays use shared memory, Julia does not know about this. We use diff --git a/moment_kinetics/src/communication.jl b/moment_kinetics/src/communication.jl index f8fa46fff..394446a4a 100644 --- a/moment_kinetics/src/communication.jl +++ b/moment_kinetics/src/communication.jl @@ -544,6 +544,9 @@ end end """ +Type used to declare a shared-memory array. When debugging is not active `MPISharedArray` +is just an alias for `Array`, but when `@debug_shared_array` is activated, it is instead +defined as an alias for `DebugMPISharedArray`. """ const MPISharedArray = @debug_shared_array_ifelse(DebugMPISharedArray, Array) From 5b036b82ff3babfa79639fd066dd4d90a0b6ff2b Mon Sep 17 00:00:00 2001 From: John Omotani Date: Mon, 16 Sep 2024 21:52:31 +0100 Subject: [PATCH 7/7] When we catch an exception, print it to stderr before calling MPI.Abort Previously, we were printing to stdout, but this makes it more likely for the error message and stack trace to get lost, e.g. when `quietoutput()` is used in the tests. --- moment_kinetics/src/moment_kinetics.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/moment_kinetics/src/moment_kinetics.jl b/moment_kinetics/src/moment_kinetics.jl index 4c57e6ed8..2fe755af2 100644 --- a/moment_kinetics/src/moment_kinetics.jl +++ b/moment_kinetics/src/moment_kinetics.jl @@ -148,8 +148,8 @@ function run_moment_kinetics(to::Union{TimerOutput,Nothing}, input_dict=Dict(); # Stop code from hanging when running on multiple processes if only one of them # throws an error if global_size[] > 1 - println("$(typeof(e)) on process $(global_rank[]):") - showerror(stdout, e, catch_backtrace()) + println(stderr, "$(typeof(e)) on process $(global_rank[]):") + showerror(stderr, e, catch_backtrace()) flush(stdout) flush(stderr) MPI.Abort(comm_world, 1)