From fead87e148b1922416cedfa2ce78574ec57a408f Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 31 Oct 2025 16:44:47 +0000 Subject: [PATCH 01/12] Switch transform method to using GuardedOptions Note that components still aren't setting up any permissions, so there would be runtime errors when executing the transform methods. Furthermore, there are some missing methods I realise I need for GuardedOptions. The unit tests are also failing to compile, although I don't understand what the problem is for them, yet. --- src/amjuel_reaction.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/amjuel_reaction.cxx b/src/amjuel_reaction.cxx index bbb9bff9e..f0fe74b98 100644 --- a/src/amjuel_reaction.cxx +++ b/src/amjuel_reaction.cxx @@ -120,7 +120,7 @@ void AmjuelReaction::transform_additional(GuardedOptions& state, Field3D& reacti Field3D T_e = get(electron["temperature"]); const int e_pop_change = this->parser->get_stoich().at("e"); if (e_pop_change != 0) { - if (electron.isSet("velocity")) { + if (IS_SET(electron["velocity"])) { // Transfer of electron kinetic to thermal energy due to density source // For ionisation: // Electrons with zero average velocity are created, diluting the kinetic energy. From 9e5775d765442ff5161b451ff9551ada707527ab Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Wed, 26 Nov 2025 13:44:02 +0000 Subject: [PATCH 02/12] Added unit test for topological sorting of components --- include/permissions.hxx | 2 + src/component_scheduler.cxx | 12 ++ src/permissions.cxx | 6 + tests/unit/test_component_scheduler.cxx | 188 ++++++++++++++++++++++++ 4 files changed, 208 insertions(+) diff --git a/include/permissions.hxx b/include/permissions.hxx index f648e85a9..33e4531b7 100644 --- a/include/permissions.hxx +++ b/include/permissions.hxx @@ -337,3 +337,5 @@ std::ostream& operator<<(std::ostream& os, const Permissions& permissions); /// behaviour if the input is corrupted; an exception may be thrown or /// the permissions that are read may be incomplete. std::istream& operator>>(std::istream& is, Permissions& permissions); + +std::string toString(const Permissions& value); diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index abaf9e6b5..4fb39e8f8 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -79,6 +79,18 @@ ComponentScheduler::ComponentScheduler(Options &scheduler_options, } } +// Use index to identify component +// std::map> dependencies; +// std::map var_written_by, var_final_written_by, read_dependencies; + +// Assemble the maps of which components write each variable +// Assemble read_dependencies using values from var_final_written_by, if present, else +// from var_written_by For each component +// For each write-final var, create dependency between component and all components that +// write it For each read-if-set var that has a write-final, create dependency between +// component and (final) writer(s) For each read-var create a dependency between +// component and (final) writers; if none then raise exception + std::unique_ptr ComponentScheduler::create(Options &scheduler_options, Options &component_options, Solver *solver) { diff --git a/src/permissions.cxx b/src/permissions.cxx index 04787b879..fa4adb845 100644 --- a/src/permissions.cxx +++ b/src/permissions.cxx @@ -242,3 +242,9 @@ std::istream& operator>>(std::istream& is, Permissions& permissions) { return is; } + +std::string toString(const Permissions& value) { + std::ostringstream ss; + ss << value; + return ss.str(); +} diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index c4e3c6f03..d8b106005 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -27,8 +27,26 @@ struct TestMultiply : public Component { } }; +struct OrderChecker : public Component { + OrderChecker(const std::string& name, Options& alloptions, Solver*) + : Component(getPermissions(name, alloptions)), name(name) {} + static Permissions getPermissions(const std::string& name, Options& alloptions) { + return alloptions[name]["permissions"].as(); + } + static void resetOrderInfo() { execution_order.clear(); } + + std::string name; + static std::vector execution_order; + +private: + void transform_impl(GuardedOptions&) override { execution_order.push_back(name); } +}; + +std::vector OrderChecker::execution_order; + RegisterComponent registertestcomponent("testcomponent"); RegisterComponent registertestcomponent2("multiply"); +RegisterComponent registercomponentorderchecker("orderchecker"); } // namespace TEST(SchedulerTest, OneComponent) { @@ -66,3 +84,173 @@ TEST(SchedulerTest, SubComponents) { ASSERT_TRUE(options["answer"] == 42 * 2); } +using Parameter = std::pair>; + +class ComponentOrderTest : public testing::TestWithParam { + void SetUp() override { OrderChecker::resetOrderInfo(); } +}; + +TEST_P(ComponentOrderTest, Sorted) { + Options options = GetParam().first.copy(); + auto scheduler = ComponentScheduler::create(options, options, nullptr); + scheduler->transform(options); + EXPECT_EQ(OrderChecker::execution_order, GetParam().second); +} + +INSTANTIATE_TEST_SUITE_P( + TopologicalSort, ComponentOrderTest, + testing::Values( + Parameter({{"components", ""}}, {}), + Parameter( + {{"components", "a"}, + {"a", {{"type", "orderchecker"}, {"permissions", toString(Permissions())}}}}, + {"a"}), + Parameter({{"components", "a"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1"), readWrite("2")}))}}}}, + {"a"}), + Parameter( + {{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1"), readWrite("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1"), readOnly("2")}))}}}}, + {"a", "b"}), + Parameter({{"components", "b,a"}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("1"), readOnly("2")}))}}}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1"), readWrite("2")}))}}}}, + {"a", "b"}), + Parameter( + {{"components", "b,a,c"}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1"), readOnly("2")}))}}}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1"), readWrite("2")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("2"), readOnly("1")}))}}}}, + {"a", "c", "b"}), + Parameter({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readIfSet("1")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1")}))}}}}, + {"b", "a"}), + Parameter({{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({writeFinal("1"), readIfSet("3")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1"), readOnly("2")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1"), readWrite("2"), + readIfSet("3")}))}}}}, + {"c", "b", "a"}), + Parameter({{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({writeBoundary("1")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1", Regions::Interior)}))}}}}, + {"c", "a", "b"}), + Parameter({{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readIfSet("1", Regions::Interior), + readWrite("2", Regions::All)}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({writeFinal("1", Regions::Boundaries), + readIfSet("2", Regions::Interior)}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("1", Regions::Boundaries), + readOnly("2", Regions::Boundaries)}))}}}}, + {"b", "a", "c"}))); + +class InvalidComponentOrderTest : public testing::TestWithParam { + void SetUp() override { OrderChecker::resetOrderInfo(); } +}; + +TEST_P(InvalidComponentOrderTest, BadDAG) { + Options options = GetParam().copy(); + EXPECT_THROW(ComponentScheduler::create(options, options, nullptr), BoutException); +} + +INSTANTIATE_TEST_SUITE_P( + InvalidTopologicalSort, InvalidComponentOrderTest, + testing::Values( + // Unsatisfiable dependency + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("1"), readWrite("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("3")}))}}}}), + // Circular dependency + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("1"), readWrite("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("2"), readOnly("1")}))}}}}), + // Circular dependency from readIfSet + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readIfSet("1"), readWrite("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("2"), readWrite("1")}))}}}}), + // Unsatisfiable dependency due to only setting one region + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1", Regions::Interior)}))}}}}), + // Circular dependency on only one region + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1", Regions::Interior), + readWrite("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1"), readOnly("2")}))}}}}))); From a02d845e864d8e1bab39f47d4f321200ac6e1641 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Wed, 26 Nov 2025 16:41:21 +0000 Subject: [PATCH 03/12] Add ability to order components automatically --- include/component.hxx | 2 + src/component_scheduler.cxx | 131 +++++++++++++++++++++--- src/permissions.cxx | 2 +- tests/unit/test_component_scheduler.cxx | 4 +- 4 files changed, 124 insertions(+), 15 deletions(-) diff --git a/include/component.hxx b/include/component.hxx index bd3dfe351..4440a6d19 100644 --- a/include/component.hxx +++ b/include/component.hxx @@ -142,6 +142,8 @@ struct Component { /// must be completed or else an exception will be thrown. void declareAllSpecies(const SpeciesInformation & info); + const Permissions& getPermissions() const { return state_variable_access; } + protected: /// Set the level of access needed by this component for a particular variable. void setPermissions(const std::string& variable, diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index 4fb39e8f8..99f91a02b 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -9,6 +9,123 @@ #include "../include/component.hxx" #include "../include/component_scheduler.hxx" +/// Perform a depth-first topological sort, starting from `item`. +void topological_sort(const std::vector>& dependencies, size_t item, + std::vector& sorted, std::vector& processing, + std::vector& processed) { + if (processed[item]) { + return; + } + if (processing[item]) { + throw BoutException("Circular dependency among components."); + } + processing[item] = true; + + for (const auto dep : dependencies[item]) { + topological_sort(dependencies, dep, sorted, processing, processed); + } + processed[item] = true; + sorted.push_back(item); +} + +/// Consumes a list of components and returns a new one that has been +/// topolgically sorted to ensure variables are written and read in +/// the right order. +std::vector> +sortComponents(std::vector>&& components) { + using Var = std::pair; + std::map> nonfinal_writes, final_writes; + // Build up information on which components write each variable + for (size_t i = 0; i < components.size(); i++) { + const Permissions& permissions = components[i]->getPermissions(); + for (const auto& [name, regions] : + permissions.getVariablesWithPermission(PermissionTypes::Write)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + nonfinal_writes[{name, region}].insert(i); + } + } + } + for (const auto& [name, regions] : + permissions.getVariablesWithPermission(PermissionTypes::Final)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + final_writes[{name, region}].insert(i); + } + } + } + } + + std::vector> component_dependencies(components.size()); + + // Components which do a final write on a variable depend on all + // components which do non-final writes on that variable + for (const auto& [var, comp_indices] : final_writes) { + for (size_t i : comp_indices) { + const auto item = nonfinal_writes.find(var); + if (item != nonfinal_writes.end()) { + component_dependencies[i].merge(item->second); + } + } + } + + // Work out which component(s) last write a variable before it may be read + std::map> variable_writers = std::move(final_writes); + variable_writers.merge(std::move(nonfinal_writes)); + + for (size_t i = 0; i < components.size(); i++) { + const Permissions& permissions = components[i]->getPermissions(); + // Create dependencies between components that read variables and those that write + // them + for (const auto& [name, regions] : + permissions.getVariablesWithPermission(PermissionTypes::Read)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + const auto item = variable_writers.find({name, region}); + if (item == variable_writers.end()) { + throw BoutException( + "Required variable {} (in region {}) is not written by any component.", + name, Permissions::fundamental_regions.at(region)); + } + component_dependencies[i].merge(item->second); + } + } + } + + // Create dependencies for ReadIfSet variables only if there exist another component + // which sets them + for (const auto& [name, regions] : + permissions.getVariablesWithPermission(PermissionTypes::ReadIfSet)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + const auto item = variable_writers.find({name, region}); + if (item != variable_writers.end()) { + component_dependencies[i].merge(item->second); + } + } + } + } + } + + std::vector processing(components.size(), false), + processed(components.size(), false); + std::vector order; + + for (size_t i = 0; i < components.size(); i++) { + if (!processed[i]) { + topological_sort(component_dependencies, i, order, processing, processed); + } + } + + // Create the result with components in the desired order + std::vector> result(components.size()); + for (size_t i = 0; i < components.size(); i++) { + std::swap(result[i], components[order[i]]); + } + + return result; +} + ComponentScheduler::ComponentScheduler(Options &scheduler_options, Options &component_options, Solver *solver) { @@ -77,19 +194,9 @@ ComponentScheduler::ComponentScheduler(Options &scheduler_options, for (auto& component : components) { component->declareAllSpecies(species); } -} -// Use index to identify component -// std::map> dependencies; -// std::map var_written_by, var_final_written_by, read_dependencies; - -// Assemble the maps of which components write each variable -// Assemble read_dependencies using values from var_final_written_by, if present, else -// from var_written_by For each component -// For each write-final var, create dependency between component and all components that -// write it For each read-if-set var that has a write-final, create dependency between -// component and (final) writer(s) For each read-var create a dependency between -// component and (final) writers; if none then raise exception + components = sortComponents(std::move(components)); +} std::unique_ptr ComponentScheduler::create(Options &scheduler_options, Options &component_options, diff --git a/src/permissions.cxx b/src/permissions.cxx index fa4adb845..3fdc08482 100644 --- a/src/permissions.cxx +++ b/src/permissions.cxx @@ -187,7 +187,7 @@ Permissions::getVariablesWithPermission(PermissionTypes permission, } std::string Permissions::regionNames(const Regions regions) { - std::vector regions_present(fundamental_regions.size()); + std::vector regions_present; for (auto & [region, name] : fundamental_regions) { if ((regions & region) == region) { regions_present.push_back(name); diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index d8b106005..07a437226 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -176,7 +176,7 @@ INSTANTIATE_TEST_SUITE_P( {"permissions", toString(Permissions({readWrite("1", Regions::Interior)}))}}}}, {"c", "a", "b"}), - Parameter({{"components", "a,b,c"}, + Parameter({{"components", "b,a,c"}, {"a", {{"type", "orderchecker"}, {"permissions", @@ -192,7 +192,7 @@ INSTANTIATE_TEST_SUITE_P( {"permissions", toString(Permissions({readOnly("1", Regions::Boundaries), readOnly("2", Regions::Boundaries)}))}}}}, - {"b", "a", "c"}))); + {"a", "b", "c"}))); class InvalidComponentOrderTest : public testing::TestWithParam { void SetUp() override { OrderChecker::resetOrderInfo(); } From c903e91bbed0d981598a0fecf8e6105f805ff59a Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Wed, 26 Nov 2025 18:33:34 +0000 Subject: [PATCH 04/12] Ensure ordering of components accounts for pre-set state variables --- include/component_scheduler.hxx | 5 +++++ src/component_scheduler.cxx | 9 +++++++++ tests/unit/test_component_scheduler.cxx | 11 +++++++++++ 3 files changed, 25 insertions(+) diff --git a/include/component_scheduler.hxx b/include/component_scheduler.hxx index f848614a8..a67c27b06 100644 --- a/include/component_scheduler.hxx +++ b/include/component_scheduler.hxx @@ -57,6 +57,11 @@ public: /// Preconditioning void precon(const Options &state, BoutReal gamma); + + /// All the variable names which are pre-set in the state, before + /// any components are applied. + static const std::set predeclared_variables; + private: /// The components to be executed in order std::vector> components; diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index 99f91a02b..c016f9ae5 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -9,6 +9,10 @@ #include "../include/component.hxx" #include "../include/component_scheduler.hxx" +const std::set ComponentScheduler::predeclared_variables = { + "time", "linear", "units:inv_meters_cubed", "units:eV", "units:Tesla", + "units:seconds", "units:meters"}; + /// Perform a depth-first topological sort, starting from `item`. void topological_sort(const std::vector>& dependencies, size_t item, std::vector& sorted, std::vector& processing, @@ -79,6 +83,9 @@ sortComponents(std::vector>&& components) { // them for (const auto& [name, regions] : permissions.getVariablesWithPermission(PermissionTypes::Read)) { + std::cout << name << "\n"; + if (ComponentScheduler::predeclared_variables.count(name) > 0) + continue; for (const auto& [region, _] : Permissions::fundamental_regions) { if ((regions & region) == region) { const auto item = variable_writers.find({name, region}); @@ -96,6 +103,8 @@ sortComponents(std::vector>&& components) { // which sets them for (const auto& [name, regions] : permissions.getVariablesWithPermission(PermissionTypes::ReadIfSet)) { + if (ComponentScheduler::predeclared_variables.count(name) > 0) + continue; for (const auto& [region, _] : Permissions::fundamental_regions) { if ((regions & region) == region) { const auto item = variable_writers.find({name, region}); diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index 07a437226..d1568f3ea 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -120,6 +120,17 @@ INSTANTIATE_TEST_SUITE_P( {{"type", "orderchecker"}, {"permissions", toString(Permissions({readOnly("1"), readOnly("2")}))}}}}, {"a", "b"}), + Parameter({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1"), readWrite("2"), + readOnly("time")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1"), readOnly("2"), + readIfSet("linear"), + readOnly("units:eV")}))}}}}, + {"a", "b"}), Parameter({{"components", "b,a"}, {"b", {{"type", "orderchecker"}, From d37d2ad40d9e947e1a90631355c4996a7093181e Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Thu, 27 Nov 2025 18:15:56 +0000 Subject: [PATCH 05/12] Added support for sorting components based on section permissions --- include/permissions.hxx | 2 +- src/component_scheduler.cxx | 234 ++++++++++++++++++------ tests/unit/test_component_scheduler.cxx | 52 +++++- 3 files changed, 233 insertions(+), 55 deletions(-) diff --git a/include/permissions.hxx b/include/permissions.hxx index 33e4531b7..aa9353c4b 100644 --- a/include/permissions.hxx +++ b/include/permissions.hxx @@ -230,7 +230,6 @@ public: friend std::ostream& operator<<(std::ostream& os, const Permissions& permissions); friend std::istream& operator>>(std::istream& is, Permissions& permissions); -private: /// Returns the access rights for the most specific entry in this /// object which matches the variable name. If there are no matching /// entries then the result will indicate no access rights. The @@ -239,6 +238,7 @@ private: /// entries. VarRights bestMatchRights(const std::string& variable) const; +private: std::map variable_permissions; static const std::regex LABEL_RE; diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index c016f9ae5..62ce153f2 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -1,3 +1,4 @@ +#include #include #include #include @@ -5,6 +6,7 @@ #include #include #include // for trim, strsplit +#include #include "../include/component.hxx" #include "../include/component_scheduler.hxx" @@ -32,33 +34,176 @@ void topological_sort(const std::vector>& dependencies, size_t sorted.push_back(item); } -/// Consumes a list of components and returns a new one that has been -/// topolgically sorted to ensure variables are written and read in -/// the right order. -std::vector> -sortComponents(std::vector>&& components) { - using Var = std::pair; - std::map> nonfinal_writes, final_writes; - // Build up information on which components write each variable +std::set getParents(const std::string& name) { + std::set result; + size_t start = 0, position = name.find(":", start); + while (position != std::string::npos) { + result.insert(name.substr(0, position)); + start = position + 1; + position = name.find(":", start); + } + return result; +} + +/// Produce a map between Option paths and all variable names held within +/// that path. If a path does not refer to a section then it just maps +/// to itself. +std::map> +getVariableHierarchy(const std::vector>& components) { + // Build up a set of all variable names which are read only if they + // are set by another component + std::set conditional_names; + for (const auto& component : components) { + const Permissions& permissions = component->getPermissions(); + for (const auto& [varname, _] : + permissions.getVariablesWithPermission(PermissionTypes::ReadIfSet, true)) { + conditional_names.insert(varname); + } + } + + // Build up a set of all variable names which are definitely + // read/written by components, and the sections which they imply + // exist + std::set unconditional_names, unconditional_sections; + for (const auto& component : components) { + const Permissions& permissions = component->getPermissions(); + for (const auto& [varname, _] : + permissions.getVariablesWithPermission(PermissionTypes::Read, false)) { + unconditional_names.insert(varname); + unconditional_sections.merge(getParents(varname)); + } + } + + // Split the set of all variable names used by components into + // those that are sections and those that are not. + std::set sections_present, non_sections; + std::set_intersection(unconditional_names.begin(), unconditional_names.end(), + unconditional_sections.begin(), unconditional_sections.end(), + std::inserter(sections_present, sections_present.begin())); + std::set_difference(unconditional_names.begin(), unconditional_names.end(), + sections_present.begin(), sections_present.end(), + std::inserter(non_sections, non_sections.begin())); + + std::map> result; + + // ReadIfSet variables will only actually be used if they are + // reference elsewhere. We create them with empty sets, which will + // get filled if they are present. + for (const auto& name : conditional_names) { + result.insert({name, {name}}); + } + + // Non-sections map to themselves + for (const auto& name : non_sections) { + result[name] = {name}; + } + // Sections map to those variables which they contain + for (const auto& section : sections_present) { + auto& children = result[section]; + const std::string sec_suffixed = section + ':'; + for (const auto& name : non_sections) { + if (name.rfind(sec_suffixed, 0) == 0) { + children.insert(name); + } + } + } + + return result; +} + +/// Get all variables to which the name could be referring (e.g., its +/// children if it is a section name). These will be filtered to +/// remove any variables for which more specific permissions are +/// given. +std::set +expandVariableName(const std::map>& hierarchy, + const Permissions& permissions, const std::string& name) { + const std::set& candidates = hierarchy.at(name); + std::set result; + // Only return the values that do not have a more specific permission + std::copy_if(candidates.begin(), candidates.end(), + std::inserter(result, result.begin()), + [&permissions, &name](const std::string& candidate) -> bool { + return permissions.bestMatchRights(candidate).name == name; + }); + return result; +} + +using Var = std::pair; + +std::map> +getPermissionComponentMap(const std::vector>& components, + const std::map>& hierarchy, + PermissionTypes permission) { + std::map> result; for (size_t i = 0; i < components.size(); i++) { const Permissions& permissions = components[i]->getPermissions(); for (const auto& [name, regions] : - permissions.getVariablesWithPermission(PermissionTypes::Write)) { - for (const auto& [region, _] : Permissions::fundamental_regions) { - if ((regions & region) == region) { - nonfinal_writes[{name, region}].insert(i); + permissions.getVariablesWithPermission(permission)) { + for (const auto& sub_name : expandVariableName(hierarchy, permissions, name)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + result[{sub_name, region}].insert(i); + } } } } + } + return result; +} + +/// Modifies component_dependencies to include information on which +/// variables each component reads. Returns a set of any read +/// variables which are not written by any component. +std::set +setReadDependencies(const std::vector>& components, + const std::map>& hierarchy, + const std::map>& writers, + PermissionTypes permission, + std::vector>& component_dependencies) { + std::set missing; + for (size_t i = 0; i < components.size(); i++) { + const Permissions& permissions = components[i]->getPermissions(); + // Create dependencies between components that read variables and those that write + // them for (const auto& [name, regions] : - permissions.getVariablesWithPermission(PermissionTypes::Final)) { - for (const auto& [region, _] : Permissions::fundamental_regions) { - if ((regions & region) == region) { - final_writes[{name, region}].insert(i); + permissions.getVariablesWithPermission(permission)) { + if (ComponentScheduler::predeclared_variables.count(name) > 0) + continue; + for (const auto& sub_name : expandVariableName(hierarchy, permissions, name)) { + for (const auto& [region, _] : Permissions::fundamental_regions) { + if ((regions & region) == region) { + const auto item = writers.find({sub_name, region}); + if (item == writers.end()) { + missing.insert( + fmt::format("{} ({})", sub_name, Permissions::regionNames(region))); + } else { + component_dependencies[i].insert(item->second.begin(), item->second.end()); + } + } } } } } + return missing; +} + +/// Consumes a list of components and returns a new one that has been +/// topolgically sorted to ensure variables are written and read in +/// the right order. +std::vector> +sortComponents(std::vector>&& components) { + // Map variables to the components that write them + std::map> variable_hierarchy = + getVariableHierarchy(components); + + // Get information on which components write each variable + std::map> nonfinal_writes = getPermissionComponentMap( + components, variable_hierarchy, + PermissionTypes::Write), + final_writes = getPermissionComponentMap( + components, variable_hierarchy, + PermissionTypes::Final); std::vector> component_dependencies(components.size()); @@ -68,6 +213,12 @@ sortComponents(std::vector>&& components) { for (size_t i : comp_indices) { const auto item = nonfinal_writes.find(var); if (item != nonfinal_writes.end()) { + // Note that calling merge actually removes the items from + // nonfinal_writes. This is fine because the only remaining + // thing for which we will use nonfinal_writes is setting up + // variable_writers. That doesn't use any information on + // nonfinal writes for variables which have a final write, so + // it won't do any harm to remove it. component_dependencies[i].merge(item->second); } } @@ -77,49 +228,26 @@ sortComponents(std::vector>&& components) { std::map> variable_writers = std::move(final_writes); variable_writers.merge(std::move(nonfinal_writes)); - for (size_t i = 0; i < components.size(); i++) { - const Permissions& permissions = components[i]->getPermissions(); - // Create dependencies between components that read variables and those that write - // them - for (const auto& [name, regions] : - permissions.getVariablesWithPermission(PermissionTypes::Read)) { - std::cout << name << "\n"; - if (ComponentScheduler::predeclared_variables.count(name) > 0) - continue; - for (const auto& [region, _] : Permissions::fundamental_regions) { - if ((regions & region) == region) { - const auto item = variable_writers.find({name, region}); - if (item == variable_writers.end()) { - throw BoutException( - "Required variable {} (in region {}) is not written by any component.", - name, Permissions::fundamental_regions.at(region)); - } - component_dependencies[i].merge(item->second); - } - } - } - - // Create dependencies for ReadIfSet variables only if there exist another component - // which sets them - for (const auto& [name, regions] : - permissions.getVariablesWithPermission(PermissionTypes::ReadIfSet)) { - if (ComponentScheduler::predeclared_variables.count(name) > 0) - continue; - for (const auto& [region, _] : Permissions::fundamental_regions) { - if ((regions & region) == region) { - const auto item = variable_writers.find({name, region}); - if (item != variable_writers.end()) { - component_dependencies[i].merge(item->second); - } - } - } - } + // Create dependency information for read variables + std::set missing = + setReadDependencies(components, variable_hierarchy, variable_writers, + PermissionTypes::Read, component_dependencies); + if (missing.size() > 0) { + throw BoutException( + "The following required variables are not written by any component:\n\t{}\n", + fmt::format("{}", fmt::join(missing, "\n\t"))); } + // If can not find a place where a ReadIfSet variable is written, it will just be + // skipped + setReadDependencies(components, variable_hierarchy, variable_writers, + PermissionTypes::ReadIfSet, component_dependencies); + // Create ancillary variables for sorting process std::vector processing(components.size(), false), processed(components.size(), false); std::vector order; + // Perform the sort for (size_t i = 0; i < components.size(); i++) { if (!processed[i]) { topological_sort(component_dependencies, i, order, processing, processed); diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index d1568f3ea..5882dbd64 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -175,6 +175,45 @@ INSTANTIATE_TEST_SUITE_P( {"permissions", toString(Permissions({readWrite("1"), readWrite("2"), readIfSet("3")}))}}}}, {"c", "b", "a"}), + Parameter({{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({writeFinal("1:1_1"), readWrite("1:1_2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1:1_1")}))}}}}, + {"b", "a", "c"}), + Parameter( + {{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({writeFinal("1"), readOnly("2")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("2:2_1"), writeFinal("2:2_2")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readWrite("1"), readOnly("2:2_1"), + readIfSet("3")}))}}}}, + {"b", "c", "a"}), + Parameter({{"components", "a,b,c"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readOnly("1"), readWrite("1:1_1")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", + toString(Permissions({readWrite("1:1_1"), readWrite("1:1_2")}))}}}, + {"c", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({readOnly("1:1_2")}))}}}}, + {"b", "a", "c"}), Parameter({{"components", "a,b,c"}, {"a", {{"type", "orderchecker"}, @@ -203,7 +242,18 @@ INSTANTIATE_TEST_SUITE_P( {"permissions", toString(Permissions({readOnly("1", Regions::Boundaries), readOnly("2", Regions::Boundaries)}))}}}}, - {"a", "b", "c"}))); + {"a", "b", "c"}), + Parameter({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({ + readOnly("1"), + }))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({writeFinal("1:1_1"), + readIfSet("1:1_2")}))}}}}, + {"b", "a"}))); class InvalidComponentOrderTest : public testing::TestWithParam { void SetUp() override { OrderChecker::resetOrderInfo(); } From a6c9454735ad05b6ccdad15c4c14bbbcdb5dd9de Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Tue, 2 Dec 2025 14:45:39 +0000 Subject: [PATCH 06/12] Mention new sorting and remove references to component ordering --- docs/sphinx/boundary_conditions.rst | 4 ++-- docs/sphinx/closure.rst | 7 +++--- docs/sphinx/developer.rst | 33 +++++++++++++++++++++++------ docs/sphinx/equations.rst | 33 ++++++++++------------------- docs/sphinx/examples.rst | 10 ++++----- docs/sphinx/inputs.rst | 5 +++-- 6 files changed, 50 insertions(+), 42 deletions(-) diff --git a/docs/sphinx/boundary_conditions.rst b/docs/sphinx/boundary_conditions.rst index 9cfc2ef33..8a5c4798f 100644 --- a/docs/sphinx/boundary_conditions.rst +++ b/docs/sphinx/boundary_conditions.rst @@ -250,7 +250,7 @@ The boundary fluxes might be set by sheath boundary conditions, which potentially depend on the density and temperature of all species. Recycling therefore can't be calculated until all species boundary conditions have been set. It is therefore expected that this component is a top-level -component (i.e. in the `Hermes` section) which comes after boundary conditions are set. +component (i.e. in the `Hermes` section). Recycling has been implemented at the target, the SOL edge and the PFR edge. Each is off by default and must be activated with a separate flag. Each can be @@ -620,4 +620,4 @@ Note that if you have the density controller enabled, it will work to counteract function = 0.01 source = Pe:source source_time_dependent = true - source_prefactor = Pe:source_prefactor \ No newline at end of file + source_prefactor = Pe:source_prefactor diff --git a/docs/sphinx/closure.rst b/docs/sphinx/closure.rst index b235ea71d..e9c0c4b6e 100644 --- a/docs/sphinx/closure.rst +++ b/docs/sphinx/closure.rst @@ -262,8 +262,7 @@ Input This top-level component calculates the frictional forces between each pair of species for which collisional frequencies have been calculated -(see `Braginskii Collisions component`_). As such, it must be run after -`braginskii_collisions`. If the option `frictional_heating` is +(see `Braginskii Collisions component`_). If the option `frictional_heating` is enabled then it will also calculate the energy source arising from friction. @@ -346,8 +345,8 @@ Braginskii Heat Exchange Input ----- This top-level component calculates the heat exchange between species -due to collisions (see `Braginskii Collisions component`_). As such, it must be run after -`braginskii_collisions`. There are no configurations for this component. +due to collisions (see `Braginskii Collisions component`_). There are no +configurations for this component. Theory ------ diff --git a/docs/sphinx/developer.rst b/docs/sphinx/developer.rst index f2d9a9a12..c89f046b2 100644 --- a/docs/sphinx/developer.rst +++ b/docs/sphinx/developer.rst @@ -577,6 +577,17 @@ The `name` is a string labelling the instance. The `alloptions` tree contains at * `alloptions[name]` options for this instance * `alloptions['units']` + +All component constructors must pass a `Permissions` object to the +constructor on the `Component::Component` base class. This specifies +which variables will be read/written by the `Component::transform` +method and will be used to construct a `GuardedOptions` object to be +passed into `Component::transform_impl`. The `Permissions` object +(`Component::state_variable_access`) can be further updated in the +body of the constructor of your component, based on what +configurations were specified in ``alloptions``. You should give read +and write permissions to the minimum number of variables necessary, to +avoid circular dependencies arising among components. Component Permissions @@ -620,9 +631,14 @@ and then in `Hermes::rhs` the components are run by a call:: scheduler->transform(state); The call to `ComponentScheduler::create` treats the "components" -option as a comma-separated list of names. The order of the components -is the order that they are run in. For each name in the list, the -scheduler looks up the options under the section of that name. +option as a comma-separated list of names. For each name in the list, +the scheduler looks up the options under the section of that name. The +``ComponentScheduler`` will use permission information stored by each +component in `Component::state_variable_access` to work out the order +to execute components. It will ensure that all writes to a variable +have occurred before the first time it is read. If there is a variable +needed by some component which is never set or if there is a circular +dependency then it will throw an exception. .. code-block:: ini @@ -639,8 +655,9 @@ scheduler looks up the options under the section of that name. This would create two `Component` objects, of type `component1` and `component2`. Each time `Hermes::rhs` is run, the `transform` -functions of `component1` and then `component2` will be called, -followed by their `finally` functions. +functions of `component1` and `component2` will be called, with the +order depending on what state variables each reads and writes. This is +followed by a call to their `finally` functions. It is often useful to group components together, for example to define the governing equations for different species. A `type` setting @@ -662,8 +679,10 @@ of components # options to control component3 This will create three components, which will be run in the order -`component1`, `component2`, `component3`: First all the components -in `group1`, and then `component3`. +`component1`, `component2`, `component3`: First all the components in +`group1`, and then `component3`. Grouped components will be sorted +individually when determining the run order; they may not be run +together. .. doxygenclass:: ComponentScheduler :members: diff --git a/docs/sphinx/equations.rst b/docs/sphinx/equations.rst index f04795ea2..eb59424de 100644 --- a/docs/sphinx/equations.rst +++ b/docs/sphinx/equations.rst @@ -120,9 +120,8 @@ quasineutral ~~~~~~~~~~~~ This component sets the density of one species, so that the overall -charge density is zero everywhere. This must therefore be done after -all other charged species densities have been calculated. It only -makes sense to use this component for species with a non-zero charge. +charge density is zero everywhere. It only makes sense to use this +component for species with a non-zero charge. .. doxygenstruct:: Quasineutral :members: @@ -229,8 +228,8 @@ energy, :math:`\mathcal{E}`: \mathcal{E} = \frac{1}{\gamma - 1} P + \frac{1}{2}m nv_{||}^2 Note that this component requires the parallel velocity :math:`v_{||}` -to calculate the pressure. It must therefore be listed after a component -that sets the velocity, such as `evolve_momentum`: +to calculate the pressure. It must therefore be listed alongside a +component that sets the velocity, such as `evolve_momentum`: .. code-block:: ini @@ -294,13 +293,6 @@ conduction for all species that use :ref:`evolve_pressure` or desired for a particular species then it can be turned off by setting `thermal_conduction = false` in the input options for that species. -This component requires a collision time to have been calculated -(i.e., with the :ref:`Braginskii Collisions` component). It is -recommended that this be one of the last component to run, to ensure density, -pressure, and temperature have their final values. However, it must be -run before :ref:`Recycling`, as that component will need to use the -`energy_flow_ylow` value, to which conduction contributes. - The choice of collision frequency used for conduction is set by the flag `conduction_collisions_mode`: `multispecies` uses all available collision frequencies involving the chosen species, while `braginskii` @@ -382,7 +374,7 @@ using flows already calculated for other species. It is used like `quasineutral` .. code-block:: ini [hermes] - components = h+, ..., e, ... # Note: e after all other species + components = h+, ..., e, ... [e] type = ..., zero_current,... # Set e:velocity @@ -436,7 +428,7 @@ The electron parallel viscosity is \eta_e = \frac{4}{3} 0.73 p_e \tau_e where :math:`\tau_e` is the electron collision time. The collisions between electrons -and all other species therefore need to be calculated before this component is run: +and all other species therefore needs to be calculated: .. code-block:: ini @@ -451,9 +443,8 @@ and all other species therefore need to be calculated before this component is r braginskii_ion_viscosity ~~~~~~~~~~~~~~~~~~~~~~~~ -Adds ion viscosity terms to all charged species that are not electrons. -The collision frequency is required so this is a top-level component that -must be calculated after collisions: +Adds ion viscosity terms to all charged species that are not +electrons, calculated using collision frequencies. .. code-block:: ini @@ -623,8 +614,7 @@ has cross-field transport. This discrepancy is due to historical reasons and wil 1D: neutral_parallel_diffusion ~~~~~~~~~~~~~~~~~~~~~~~~~~ -This adds diffusion to **all** neutral species (those with no or zero charge), -because it needs to be calculated after the collision frequencies are known. +This adds diffusion to **all** neutral species (those with no or zero charge). .. code-block:: ini @@ -982,9 +972,8 @@ electrostatic potential :math:`\phi` in the frame of the fluid, with an ion diamagnetic contribution. This is calculated by inverting a Laplacian equation similar to that solved in the vorticity equation. -This component needs to be run after all other currents have been -calculated. It marks currents as used, so out-of-order modifications -should raise errors. +This component will be run after all other currents have been +calculated. See the `examples/blob2d-vpol` example, which contains: diff --git a/docs/sphinx/examples.rst b/docs/sphinx/examples.rst index 92345d870..e9d3d588f 100644 --- a/docs/sphinx/examples.rst +++ b/docs/sphinx/examples.rst @@ -276,9 +276,10 @@ so that the equation solved is where :math:`T_e` is the fixed electron temperature (5eV). -The :ref:`vorticity` component uses the pressure to calculate the diamagnetic current, -so must come after the `e` component. This component then calculates the potential. -Options to control the vorticity component are set in the `[vorticity]` section. +The :ref:`vorticity` component uses the pressure to calculate the +diamagnetic current. This component then calculates the potential. +Options to control the vorticity component are set in the +`[vorticity]` section. .. math:: @@ -287,8 +288,7 @@ Options to control the vorticity component are set in the `[vorticity]` section. \nabla\cdot\left(\frac{1}{B^2}\nabla_\perp\phi\right) = \omega \end{aligned} -The `sheath_closure` component uses the potential, so must come after :ref:`vorticity`. -Options are also set as +The `sheath_closure` component uses the potential. Options are also set as .. code-block:: ini diff --git a/docs/sphinx/inputs.rst b/docs/sphinx/inputs.rst index ff7a936b7..fbaac1d7f 100644 --- a/docs/sphinx/inputs.rst +++ b/docs/sphinx/inputs.rst @@ -19,8 +19,9 @@ number of output timesteps ``nout`` and the output timestep size Note that the solver timestep is adaptive and not user-settable. This is followed by ``[mesh]``, ``[solver]`` and ``[hermes]`` headers, where the ``[hermes]`` -section defines the list of components used. The component order -matters, as the components are executed in order. +section defines the list of components used. The component order doesn't +matters, as they will be sorted to ensure state variables are set +before they are used. .. code-block:: ini From a16afd61e64f91d1e09c547954d79a06aa969a4a Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Tue, 2 Dec 2025 15:46:36 +0000 Subject: [PATCH 07/12] Check results are the same regardless of component order --- CMakeLists.txt | 1 + .../integrated/component-order/data/BOUT.inp | 164 ++++++++++++++++++ tests/integrated/component-order/runtest | 93 ++++++++++ 3 files changed, 258 insertions(+) create mode 100644 tests/integrated/component-order/data/BOUT.inp create mode 100755 tests/integrated/component-order/runtest diff --git a/CMakeLists.txt b/CMakeLists.txt index f2e1bfaec..91e0c4ac3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -333,6 +333,7 @@ if(HERMES_TESTS) hermes_add_integrated_test(alfven-wave) hermes_add_integrated_test(collfreq-braginskii-afn) hermes_add_integrated_test(collfreq-multispecies) + hermes_add_integrated_test(component-order) # Unit tests option(HERMES_UNIT_TESTS "Build the unit tests" ON) diff --git a/tests/integrated/component-order/data/BOUT.inp b/tests/integrated/component-order/data/BOUT.inp new file mode 100644 index 000000000..2614b19c8 --- /dev/null +++ b/tests/integrated/component-order/data/BOUT.inp @@ -0,0 +1,164 @@ +nout = 0 +timestep = 1 + + +[mesh] +# 1D simulation, use "y" as the dimension along the fieldline +nx = 10 +ny = 10 # Resolution along field-line +nz = 10 +ixseps1 = 0 +J = 1 + + +[hermes] +components = (d+, he+, he, d, e, + sheath_boundary_simple, braginskii_collisions, braginskii_friction, + braginskii_heat_exchange, reactions, electron_force_balance, + neutral_parallel_diffusion, braginskii_ion_viscosity, + braginskii_conduction, recycling) + + +[solver] +type = pvode + + +[sheath_boundary_simple] +lower_y = false +upper_y = true +diagnose = true + + +[neutral_parallel_diffusion] +dneut = 10 +diffusion_collisions_mode = multispecies + + +[braginskii_collisions] +electron_ion = true +electron_electron = true +electron_neutral = true +ion_ion = true +ion_neutral = true +neutral_neutral = true +diagnose = true + + +[braginskii_ion_viscosity] +viscosity_collisions_mode = multispecies + +[recycling] +species = d+, he+ + + +#################################### +[d+] +type = (evolve_density, evolve_pressure, evolve_momentum, + noflow_boundary) + +charge = 1 +AA = 2 +thermal_conduction = true +conduction_collisions_mode = multispecies +noflow_upper_y = false +recycle_as = d +target_recycle = true +diagnose = true + +[Nd+] +function = 1 + +[Pd+] +function = 1 + +[NVd+] +function = 1 + +#################################### +[he+] +type = (evolve_density, evolve_pressure, evolve_momentum, + noflow_boundary) + +charge = 2 +AA = 4 +thermal_conduction = true +conduction_collisions_mode = multispecies +noflow_upper_y = false +recycle_as = he +target_recycle = true +diagnose = true + +[Nhe+] +function = 0.01 + +[Phe+] +function = 0.001 + +[NVhe+] +function = 0.001 + +#################################### +[he] +type = (neutral_mixed, noflow_boundary) + +AA = 4 +diffusion_collisions_mode = multispecies +diagnose = true + +[Nhe] +function = 0.001 + +[Phe] +function = 0.0001 + +[NVhe] +function = 0.0001 + +#################################### +[d] +type = (neutral_mixed, noflow_boundary) + +AA = 2 +diffusion_collisions_mode = multispecies +diagnose = true + +[Nd] +function = 0.001 + +[Pd] +function = 0.0001 + +[NVd] +function = 0.0001 + +#################################### +[e] # Electrons +type = (quasineutral, evolve_pressure, zero_current, + noflow_boundary) + +noflow_upper_y = false + +charge = -1 +AA = 1/1836 +thermal_conduction = true # in evolve_pressure +conduction_collisions_mode = multispecies + +diagnose = true + +[Pe] + +function = `Pd+:function` + +#################################### + + +[reactions] +diagnose = true +type = ( + d + e -> d+ + 2e, # Deuterium ionisation + d+ + e -> d, # Deuterium recombination + d + d+ -> d+ + d, # Charge exchange + + he + e -> he+ + 2e, # Helium ionisation + he+ + e -> he, # Helium+ recombination + ) diff --git a/tests/integrated/component-order/runtest b/tests/integrated/component-order/runtest new file mode 100755 index 000000000..8b653027d --- /dev/null +++ b/tests/integrated/component-order/runtest @@ -0,0 +1,93 @@ +#!/usr/bin/env python3 + +# Python script to run and analyse two test cases that are identical +# except the order in which the components are specified. + +from __future__ import division +from __future__ import print_function + +try: + from builtins import str +except: + pass + +from boututils.run_wrapper import shell, launch_safe, getmpirun +from boutdata.collect import collect +import numpy +import xhermes + +from numpy import sqrt, max, abs, mean, array, log, concatenate +from pathlib import Path + + +def check_value(name, val, target, **kws): + # override default tolerances to get desired behaviour + kws["atol"] = kws.pop("atol", 0.0) + kws["rtol"] = kws.pop("rtol", 0.0) + + success = numpy.isclose(val, target, **kws) + if not success: + # On failure, report the effective abs tolerance used by isclose + atol_eff = kws["atol"] + kws["rtol"] * abs(target) + print( + f"Expected\n {target-atol_eff} < {name} < {target+atol_eff}\nBut actual value was\n {val}" + ) + return success + + +cases = [ + "d+, he+, he, d, e, sheath_boundary_simple, braginskii_collisions, braginskii_friction, braginskii_heat_exchange, reactions, electron_force_balance, neutral_parallel_diffusion, braginskii_ion_viscosity, braginskii_conduction, recycling", + "recycling, d+, braginskii_friction, braginskii_heat_exchange, he, reactions, electron_force_balance, neutral_parallel_diffusion, d, braginskii_ion_viscosity, sheath_boundary_simple, he+, braginskii_conduction, braginskii_collisions, e" +] + +# Remove old test results and symlink executable + +if not Path("hermes-3").is_file(): + shell("ln -s ../../../hermes-3 hermes-3") + +results = [] + +# Run and exit if unsuccessful +for i, c in enumerate(cases): + for file in ["BOUT.dmp.0.nc", "BOUT.log.0", "BOUT.restart.0.nc", "BOUT.settings", ".BOUT.pid.0"]: + path = Path("data") / file + if path.is_file(): + print(f"Removing old file data/{file}") + path.unlink() + + s, out = launch_safe(f"./hermes-3 hermes:components='{c}'", nproc=1, pipe=True) + + if s != 0: + print(" => Test failed: ") + print(out) + exit(1) + + # Save output to log file + with open(f"run.{i}.log", "w") as f: + f.write(out) + + results.append(xhermes.open("data", unnormalise=False).isel(t=-1)) + +success = True + +rtol = 1e-10 +atol = 1e-10 +diagnostics = ["Nd+", "Pd+", "Vd+", "NVd+", "Nd", "Pd", "Vd", "NVd", "Nhe+", "Phe+", "Vhe+", "NVhe+", "Nhe", "Phe", "Vhe", "NVhe", "Ne", "Pe", "Ve"] +idx = (3, 5, 7) + +failing = [] + +# Check output is the same between the simulations +for d in diagnostics: + if not check_value(d, results[0][d].values[idx], results[1][d].values[idx], rtol=rtol, atol=atol): + success = False + failing.append(d) + + +# Final output +if success: + print(" => Test passed") + exit(0) +else: + print(f" => Test failed: \nDisagreement in diagnostic(s) {', '.join(failing)}") + exit(1) From 7d1851d8f665817caa6a6e97939d208dd2fcac7b Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Tue, 2 Dec 2025 15:57:01 +0000 Subject: [PATCH 08/12] Check only one component can make final write --- src/component_scheduler.cxx | 4 ++++ tests/unit/test_component_scheduler.cxx | 8 ++++++++ 2 files changed, 12 insertions(+) diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index 62ce153f2..16a540521 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -210,6 +210,10 @@ sortComponents(std::vector>&& components) { // Components which do a final write on a variable depend on all // components which do non-final writes on that variable for (const auto& [var, comp_indices] : final_writes) { + if (comp_indices.size() > 1) { + throw BoutException( + "Multiple components have permission to make final write to variable {}", var); + } for (size_t i : comp_indices) { const auto item = nonfinal_writes.find(var); if (item != nonfinal_writes.end()) { diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index 5882dbd64..c7e40be13 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -276,6 +276,14 @@ INSTANTIATE_TEST_SUITE_P( {"b", {{"type", "orderchecker"}, {"permissions", toString(Permissions({readWrite("3")}))}}}}), + // Multiple final writes + Options({{"components", "a,b"}, + {"a", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({writeFinal("1")}))}}}, + {"b", + {{"type", "orderchecker"}, + {"permissions", toString(Permissions({writeFinal("1")}))}}}}), // Circular dependency Options({{"components", "a,b"}, {"a", From 7559d27c5d74f3dee388aa0dc9128107f3460975 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Wed, 3 Dec 2025 13:10:03 +0000 Subject: [PATCH 09/12] Fixed some incorrect access controls causing errors --- include/permissions.hxx | 27 +++++++++++++++++++++++-- src/braginskii_conduction.cxx | 4 ++-- src/neutral_boundary.cxx | 4 ++-- src/sheath_boundary.cxx | 11 +++++----- src/sheath_boundary_insulating.cxx | 13 ++++++------ src/sheath_boundary_simple.cxx | 13 ++++++------ tests/unit/test_component_scheduler.cxx | 2 +- tests/unit/test_permissions.cxx | 19 ++++++++++++++++- 8 files changed, 68 insertions(+), 25 deletions(-) diff --git a/include/permissions.hxx b/include/permissions.hxx index aa9353c4b..0d9b7f197 100644 --- a/include/permissions.hxx +++ b/include/permissions.hxx @@ -290,16 +290,25 @@ inline Permissions::VarRights writeFinal(std::string varname, } /// Convenience function to return an object expressing that the -/// variable should have Final write permissions on the boundaries. It +/// variable should have write permissions on the boundaries. It /// will have Read permissions in the interior, as this is normally /// required to set the boundaries correctly. inline Permissions::VarRights writeBoundary(std::string varname) { + return {std::move(varname), + {Regions::Nowhere, Regions::Interior, Regions::Boundaries, Regions::Nowhere}}; +} + +/// Convenience function to return an object expressing that the +/// variable should have Final write permissions on the boundaries. It +/// will have Read permissions in the interior, as this is normally +/// required to set the boundaries correctly. +inline Permissions::VarRights writeBoundaryFinal(std::string varname) { return {std::move(varname), {Regions::Nowhere, Regions::Interior, Regions::Nowhere, Regions::Boundaries}}; } /// Convenience function to return an object expressing that, if the -/// interior has been set, then the variable should have Final write +/// interior has been set, then the variable should have write /// permissions on the boundaries and read permissions for the /// interior. /// @@ -308,6 +317,20 @@ inline Permissions::VarRights writeBoundary(std::string varname) { /// have write permission regardless of whether or not the interior is /// set. inline Permissions::VarRights writeBoundaryIfSet(std::string varname) { + return {std::move(varname), + {Regions::Interior, Regions::Nowhere, Regions::Boundaries, Regions::Nowhere}}; +} + +/// Convenience function to return an object expressing that, if the +/// interior has been set, then the variable should have Final write +/// permissions on the boundaries and read permissions for the +/// interior. +/// +/// FIXME: Currently these permissiosn are not expressed properly, due +/// to limitations in how the permission system. The boundary will +/// have write permission regardless of whether or not the interior is +/// set. +inline Permissions::VarRights writeBoundaryFinalIfSet(std::string varname) { return {std::move(varname), {Regions::Interior, Regions::Nowhere, Regions::Nowhere, Regions::Boundaries}}; } diff --git a/src/braginskii_conduction.cxx b/src/braginskii_conduction.cxx index b9a5692c7..f38d08ed2 100644 --- a/src/braginskii_conduction.cxx +++ b/src/braginskii_conduction.cxx @@ -28,7 +28,6 @@ using bout::globals::mesh; BraginskiiConduction::BraginskiiConduction(const std::string&, Options& alloptions, Solver*) : Component({readOnly("species:{sp}:{input_vars}"), - writeBoundary("species:{sp}:pressure"), readWrite("species:{sp}:{output_vars}")}) { AUTO_TRACE(); @@ -91,7 +90,8 @@ BraginskiiConduction::BraginskiiConduction(const std::string&, Options& alloptio std::vector coll_types; - substitutePermissions("input_vars", {"AA", "density", "temperature"}); + substitutePermissions("input_vars", + {"AA", "density", "temperature", "pressure"}); substitutePermissions("output_vars", {"energy_source", "kappa_par", "energy_flow_ylow"}); std::vector species; diff --git a/src/neutral_boundary.cxx b/src/neutral_boundary.cxx index c4d68f68a..53f48d760 100644 --- a/src/neutral_boundary.cxx +++ b/src/neutral_boundary.cxx @@ -7,8 +7,8 @@ using bout::globals::mesh; NeutralBoundary::NeutralBoundary(std::string name, Options& alloptions, [[maybe_unused]] Solver* solver) - : Component({writeBoundary("species:{name}:{outputs}"), - writeBoundaryIfSet("species:{name}:{conditional_outputs}"), + : Component({writeBoundaryFinal("species:{name}:{outputs}"), + writeBoundaryFinalIfSet("species:{name}:{conditional_outputs}"), readWrite("species:{name}:energy_source")}), name(name) { AUTO_TRACE(); diff --git a/src/sheath_boundary.cxx b/src/sheath_boundary.cxx index 9ade3dc03..d0d07108a 100644 --- a/src/sheath_boundary.cxx +++ b/src/sheath_boundary.cxx @@ -49,16 +49,16 @@ SheathBoundary::SheathBoundary(std::string name, Options& alloptions, Solver*) : Component({ readIfSet("species:{all_species}:charge"), readIfSet("species:e:{e_whole_domain}"), - writeBoundary("species:e:{e_boundary}"), + writeBoundaryFinal("species:e:{e_boundary}"), readWrite("species:e:energy_source"), - writeBoundaryIfSet("species:e:{e_optional}"), + writeBoundaryFinalIfSet("species:e:{e_optional}"), writeBoundaryReadInteriorIfSet("species:e:pressure"), - readIfSet("species:{ions}:adiabatic"), + readIfSet("species:{ions}:{ion_whole_domain}"), readOnly("species:{ions}:AA"), readWrite("species:{ions}:energy_source"), - writeBoundary("species:{ions}:{ion_boundary}"), + writeBoundaryFinal("species:{ions}:{ion_boundary}"), writeBoundaryReadInteriorIfSet("species:{ions}:pressure"), - writeBoundaryIfSet("species:{ions}:{ion_optional}"), + writeBoundaryFinalIfSet("species:{ions}:{ion_optional}"), }) { AUTO_TRACE(); @@ -112,6 +112,7 @@ SheathBoundary::SheathBoundary(std::string name, Options& alloptions, Solver*) substitutePermissions("e_whole_domain", {"AA", "adiabatic"}); substitutePermissions("e_boundary", {"density", "temperature"}); substitutePermissions("e_optional", {"velocity", "momentum"}); + substitutePermissions("ion_whole_domain", {"charge", "adiabatic"}); substitutePermissions("ion_boundary", {"density", "temperature"}); substitutePermissions("ion_optional", {"velocity", "momentum"}); setPermissions(always_set_phi ? writeBoundaryReadInteriorIfSet("fields:phi") diff --git a/src/sheath_boundary_insulating.cxx b/src/sheath_boundary_insulating.cxx index f3a9047f3..17098e3da 100644 --- a/src/sheath_boundary_insulating.cxx +++ b/src/sheath_boundary_insulating.cxx @@ -53,16 +53,17 @@ SheathBoundaryInsulating::SheathBoundaryInsulating(std::string name, Options& al : Component({ readIfSet("species:{all_species}:charge"), readIfSet("species:e:{e_whole_domain}"), - writeBoundary("species:e:{e_boundary}"), + writeBoundaryFinal("species:e:{e_boundary}"), readWrite("species:e:energy_source"), - writeBoundaryIfSet("species:e:{e_optional}"), + writeBoundaryFinalIfSet("species:e:{e_optional}"), writeBoundaryReadInteriorIfSet("species:e:pressure"), - readIfSet("species:{ions}:adiabatic"), + readIfSet("species:{ions}:{ion_whole_domain}"), readOnly("species:{ions}:AA"), readWrite("species:{ions}:energy_source"), - writeBoundary("species:{ions}:{ion_boundary}"), - writeBoundaryIfSet("species:{ions}:{ion_optional}"), + writeBoundaryFinal("species:{ions}:{ion_boundary}"), + writeBoundaryFinalIfSet("species:{ions}:{ion_optional}"), writeBoundaryReadInteriorIfSet("species:{ions}:pressure"), + writeBoundaryFinalIfSet("fields:phi") }) { AUTO_TRACE(); @@ -95,9 +96,9 @@ SheathBoundaryInsulating::SheathBoundaryInsulating(std::string name, Options& al substitutePermissions("e_whole_domain", {"AA", "adiabatic"}); substitutePermissions("e_boundary", {"density", "temperature"}); substitutePermissions("e_optional", {"velocity", "momentum"}); + substitutePermissions("ion_whole_domain", {"charge", "adiabatic"}); substitutePermissions("ion_boundary", {"density", "temperature"}); substitutePermissions("ion_optional", {"velocity", "momentum"}); - setPermissions(writeBoundaryIfSet("fields:phi")); } void SheathBoundaryInsulating::transform_impl(GuardedOptions& state) { diff --git a/src/sheath_boundary_simple.cxx b/src/sheath_boundary_simple.cxx index e1fd58ce5..d4057986d 100644 --- a/src/sheath_boundary_simple.cxx +++ b/src/sheath_boundary_simple.cxx @@ -60,18 +60,18 @@ BoutReal limitFree(BoutReal fm, BoutReal fc, BoutReal mode) { SheathBoundarySimple::SheathBoundarySimple(std::string name, Options& alloptions, Solver*) : Component({ readIfSet("species:e:{e_whole_domain}"), - writeBoundary("species:e:{e_boundary}"), + writeBoundaryFinal("species:e:{e_boundary}"), readWrite("species:e:energy_source"), readWrite("species:e:energy_flow_ylow"), - writeBoundaryIfSet("species:e:{e_optional}"), + writeBoundaryFinalIfSet("species:e:{e_optional}"), writeBoundaryReadInteriorIfSet("species:e:pressure"), - readIfSet("species:{all_species}:charge"), + readIfSet("species:{ions}:{ion_whole_domain}"), readOnly("species:{ions}:AA"), readWrite("species:{ions}:energy_source"), readWrite("species:{ions}:energy_flow_ylow"), - writeBoundary("species:{ions}:{ion_boundary}"), + writeBoundaryFinal("species:{ions}:{ion_boundary}"), writeBoundaryReadInteriorIfSet("species:{ions}:pressure"), - writeBoundaryIfSet("species:{ions}:{ion_optional}"), + writeBoundaryFinalIfSet("species:{ions}:{ion_optional}"), }) { AUTO_TRACE(); @@ -145,9 +145,10 @@ SheathBoundarySimple::SheathBoundarySimple(std::string name, Options& alloptions .doc("Save additional output diagnostics") .withDefault(false); - substitutePermissions("e_whole_domain", {"AA", "charge"}); + substitutePermissions("e_whole_domain", {"AA", "adiabatic"}); substitutePermissions("e_boundary", {"density", "temperature"}); substitutePermissions("e_optional", {"velocity", "momentum"}); + substitutePermissions("ion_whole_domain", {"charge", "adiabatic"}); substitutePermissions("ion_boundary", {"density", "temperature"}); substitutePermissions("ion_optional", {"velocity", "momentum"}); setPermissions(always_set_phi ? writeBoundaryReadInteriorIfSet("fields:phi") diff --git a/tests/unit/test_component_scheduler.cxx b/tests/unit/test_component_scheduler.cxx index c7e40be13..9e2df5ac0 100644 --- a/tests/unit/test_component_scheduler.cxx +++ b/tests/unit/test_component_scheduler.cxx @@ -217,7 +217,7 @@ INSTANTIATE_TEST_SUITE_P( Parameter({{"components", "a,b,c"}, {"a", {{"type", "orderchecker"}, - {"permissions", toString(Permissions({writeBoundary("1")}))}}}, + {"permissions", toString(Permissions({writeBoundaryFinal("1")}))}}}, {"b", {{"type", "orderchecker"}, {"permissions", toString(Permissions({readOnly("1")}))}}}, diff --git a/tests/unit/test_permissions.cxx b/tests/unit/test_permissions.cxx index f5796bcb3..4988fe4ff 100644 --- a/tests/unit/test_permissions.cxx +++ b/tests/unit/test_permissions.cxx @@ -126,6 +126,8 @@ TEST(PermissionsTests, TestGetHighestPermission) { {Regions::Nowhere, Regions::Nowhere, Regions::Interior, Regions::Nowhere}}, {"species:d:collision_frequencies", {Regions::Nowhere, Regions::Nowhere, Regions::Boundaries, Regions::Nowhere}}, + writeBoundary("fields:phi"), + writeBoundaryIfSet("species:he+:temperature"), }); auto no_permission = make_permission(PermissionTypes::None, ""); @@ -148,6 +150,10 @@ TEST(PermissionsTests, TestGetHighestPermission) { EXPECT_EQ(example.getHighestPermission("species:d:collision_frequencies:d_d_coll"), make_permission(PermissionTypes::None, "species:d:collision_frequencies")); EXPECT_EQ(example.getHighestPermission("unset"), no_permission); + EXPECT_EQ(example.getHighestPermission("fields:phi"), + make_permission(PermissionTypes::Read, "fields:phi")); + EXPECT_EQ(example.getHighestPermission("species:he+:temperature"), + make_permission(PermissionTypes::ReadIfSet, "species:he+:temperature")); // Get the highest permission on the boundaries EXPECT_EQ(example.getHighestPermission("species:he:charge", Regions::Boundaries), @@ -169,6 +175,10 @@ TEST(PermissionsTests, TestGetHighestPermission) { Regions::Boundaries), make_permission(PermissionTypes::Write, "species:d:collision_frequencies")); EXPECT_EQ(example.getHighestPermission("unset", Regions::Boundaries), no_permission); + EXPECT_EQ(example.getHighestPermission("fields:phi", Regions::Boundaries), + make_permission(PermissionTypes::Write, "fields:phi")); + EXPECT_EQ(example.getHighestPermission("species:he+:temperature", Regions::Boundaries), + make_permission(PermissionTypes::Write, "species:he+:temperature")); // Get the highest permission on the interior EXPECT_EQ(example.getHighestPermission("species:he:charge", Regions::Interior), @@ -190,6 +200,10 @@ TEST(PermissionsTests, TestGetHighestPermission) { Regions::Interior), make_permission(PermissionTypes::None, "species:d:collision_frequencies")); EXPECT_EQ(example.getHighestPermission("unset", Regions::Interior), no_permission); + EXPECT_EQ(example.getHighestPermission("fields:phi", Regions::Interior), + make_permission(PermissionTypes::Read, "fields:phi")); + EXPECT_EQ(example.getHighestPermission("species:he+:temperature", Regions::Interior), + make_permission(PermissionTypes::ReadIfSet, "species:he+:temperature")); // Check the permission for the "Nowhere" region is always "None" EXPECT_EQ(example.getHighestPermission("species:he:charge", Regions::Nowhere), @@ -211,6 +225,9 @@ TEST(PermissionsTests, TestGetHighestPermission) { Regions::Nowhere), no_permission); EXPECT_EQ(example.getHighestPermission("unset", Regions::Nowhere), no_permission); + EXPECT_EQ(example.getHighestPermission("fields:phi", Regions::Nowhere), no_permission); + EXPECT_EQ(example.getHighestPermission("species:he+:temperature", Regions::Nowhere), + no_permission); // Check permissions for a species that might be mistaken for one of // the sections we've given permissions for @@ -358,7 +375,7 @@ TEST(PermissionsTests, TestIO) { const Permissions empty({}); const Permissions single({readOnly("test")}); const Permissions multiple( - {readIfSet("a", Regions::Interior), writeBoundary("b"), readWrite("c:d")}); + {readIfSet("a", Regions::Interior), writeBoundaryFinal("b"), readWrite("c:d")}); Permissions new_perm; std::stringstream ss1; From 904d98748c5ab85853a531775259271c74243730 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Mon, 15 Dec 2025 15:02:59 +0000 Subject: [PATCH 10/12] Add better comments explaining component sorting process --- src/component_scheduler.cxx | 80 +++++++++++++++++++++++++++---------- 1 file changed, 58 insertions(+), 22 deletions(-) diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index 16a540521..e97581a5b 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -15,7 +15,15 @@ const std::set ComponentScheduler::predeclared_variables = { "time", "linear", "units:inv_meters_cubed", "units:eV", "units:Tesla", "units:seconds", "units:meters"}; -/// Perform a depth-first topological sort, starting from `item`. +/// Perform a depth-first topological sort, starting from `item`. It +/// will finish once it reaches the end of `item`'s dependency chain, +/// so this needs to be called in a loop for all items. Information +/// about `item` is stored in the corresponding index of the vector +/// arguments. +/// +/// In pratice, `item` here represents the index of a particular +/// component. The indices of the dependencies of `item` are stored in +/// the corresponding element of `dependencies`. void topological_sort(const std::vector>& dependencies, size_t item, std::vector& sorted, std::vector& processing, std::vector& processed) { @@ -34,6 +42,8 @@ void topological_sort(const std::vector>& dependencies, size_t sorted.push_back(item); } +/// Get all the parent sections of a variable "path". Sections are +/// separated by colons in the path. std::set getParents(const std::string& name) { std::set result; size_t start = 0, position = name.find(":", start); @@ -45,9 +55,12 @@ std::set getParents(const std::string& name) { return result; } -/// Produce a map between Option paths and all variable names held within -/// that path. If a path does not refer to a section then it just maps -/// to itself. +/// Produce a map between Option paths and all variable names held +/// within that path. If the path refers to a section then it maps to +/// the set of all variables contained in that section and any +/// sub-sections. Otherwise the path corresponds to a variable and +/// just maps to itself. Only paths which are explicitly given a +/// permission by at least one component will be present. std::map> getVariableHierarchy(const std::vector>& components) { // Build up a set of all variable names which are read only if they @@ -61,7 +74,7 @@ getVariableHierarchy(const std::vector>& components) } } - // Build up a set of all variable names which are definitely + // Build up a set of all section/variable names which are definitely // read/written by components, and the sections which they imply // exist std::set unconditional_names, unconditional_sections; @@ -74,12 +87,15 @@ getVariableHierarchy(const std::vector>& components) } } - // Split the set of all variable names used by components into - // those that are sections and those that are not. - std::set sections_present, non_sections; + /// Assemble the set of all section names which are referred to + /// explicitly in the component permissions. + std::set sections_present; std::set_intersection(unconditional_names.begin(), unconditional_names.end(), unconditional_sections.begin(), unconditional_sections.end(), std::inserter(sections_present, sections_present.begin())); + /// Assemble the set of all variable names which are definitlye + /// read/written by components (i.e., not including sections) + std::set non_sections; std::set_difference(unconditional_names.begin(), unconditional_names.end(), sections_present.begin(), sections_present.end(), std::inserter(non_sections, non_sections.begin())); @@ -90,7 +106,9 @@ getVariableHierarchy(const std::vector>& components) // reference elsewhere. We create them with empty sets, which will // get filled if they are present. for (const auto& name : conditional_names) { - result.insert({name, {name}}); + // FIXME: this isn't an empty set + //result.insert({name, {name}}); + result.insert({name, {}}); } // Non-sections map to themselves @@ -131,6 +149,8 @@ expandVariableName(const std::map>& hierarchy using Var = std::pair; +/// Create a map between a variable and the set of components that +/// access it with the specified permission level. std::map> getPermissionComponentMap(const std::vector>& components, const std::map>& hierarchy, @@ -153,8 +173,11 @@ getPermissionComponentMap(const std::vector>& compone } /// Modifies component_dependencies to include information on which -/// variables each component reads. Returns a set of any read -/// variables which are not written by any component. +/// components depend on each other. It does this by making components +/// which read a variable depend on whichever component(s) write that +/// variable (information contained in the `writers` +/// argument). Returns a set of any read variables which are not +/// written by any component. std::set setReadDependencies(const std::vector>& components, const std::map>& hierarchy, @@ -193,7 +216,10 @@ setReadDependencies(const std::vector>& components, /// the right order. std::vector> sortComponents(std::vector>&& components) { - // Map variables to the components that write them + // Map between variable/section names specified by component + // permissions and the variables they contain. In the case of + // sections this is all variables within the section and any + // sub-sections. Non-section viarables map to themselves. std::map> variable_hierarchy = getVariableHierarchy(components); @@ -205,6 +231,10 @@ sortComponents(std::vector>&& components) { components, variable_hierarchy, PermissionTypes::Final); + // Object mapping between components (reprsented by the index of + // that component in the `components` argument) and the components + // each of these depends upon (represented by a set of the indices + // for those components). std::vector> component_dependencies(components.size()); // Components which do a final write on a variable depend on all @@ -217,22 +247,27 @@ sortComponents(std::vector>&& components) { for (size_t i : comp_indices) { const auto item = nonfinal_writes.find(var); if (item != nonfinal_writes.end()) { - // Note that calling merge actually removes the items from - // nonfinal_writes. This is fine because the only remaining - // thing for which we will use nonfinal_writes is setting up - // variable_writers. That doesn't use any information on - // nonfinal writes for variables which have a final write, so - // it won't do any harm to remove it. + // Note that calling merge actually removes the items from the + // sets stored in nonfinal_writes. This is fine because the + // only remaining thing for which we will use nonfinal_writes + // is setting up variable_writers and that doesn't use any + // information on variables which have a final + // write. Therefore it won't do any harm hear to remove + // information about variables which have a final write.. component_dependencies[i].merge(item->second); } } } - // Work out which component(s) last write a variable before it may be read + // Work out which component(s) last write a variable before it may + // be read. For variables with a final-write, it is whichever + // component performs that final write. std::map> variable_writers = std::move(final_writes); + // For other variables, it is the set of all components which have + // write permission. variable_writers.merge(std::move(nonfinal_writes)); - // Create dependency information for read variables + // Insert dependency information for components that (unconditionally) read variables std::set missing = setReadDependencies(components, variable_hierarchy, variable_writers, PermissionTypes::Read, component_dependencies); @@ -241,8 +276,9 @@ sortComponents(std::vector>&& components) { "The following required variables are not written by any component:\n\t{}\n", fmt::format("{}", fmt::join(missing, "\n\t"))); } - // If can not find a place where a ReadIfSet variable is written, it will just be - // skipped + // Insert dependency information for components that read variables + // if those variables have been set. If can not find a place where + // the variable is written, it will just be skipped. setReadDependencies(components, variable_hierarchy, variable_writers, PermissionTypes::ReadIfSet, component_dependencies); From 975a912dd640d384813694e433895c628b592e08 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Fri, 2 Jan 2026 13:38:02 +0000 Subject: [PATCH 11/12] Responded to some comments by @ZedThree The comments were in PR #445, but some of them actually relate to work done in #428. Only the latter are dealt with in this PR. --- include/component.hxx | 66 +++++++++++++++++---------------- include/component_scheduler.hxx | 5 ++- src/component_scheduler.cxx | 47 +++++++++++++---------- 3 files changed, 64 insertions(+), 54 deletions(-) diff --git a/include/component.hxx b/include/component.hxx index 4440a6d19..d0b18aaf0 100644 --- a/include/component.hxx +++ b/include/component.hxx @@ -1,17 +1,7 @@ #pragma once - #ifndef HERMES_COMPONENT_H #define HERMES_COMPONENT_H -#include -#include -#include -#include -#include -#include -#include -#include - #include #include #include @@ -22,6 +12,15 @@ #include #include +#include +#include +#include +#include +#include +#include +#include +#include + #include "guarded_options.hxx" #include "permissions.hxx" #include "hermes_utils.hxx" @@ -34,8 +33,10 @@ struct SpeciesInformation { SpeciesInformation(const std::vector& electrons, const std::vector& neutrals, const std::vector& positive_ions, - const std::vector & negative_ions) - : electrons(electrons), neutrals(neutrals), positive_ions(positive_ions), negative_ions(negative_ions), ions(positive_ions) { + const std::vector& negative_ions) + : electrons(std::move(electrons)), neutrals(std::move(neutrals)), + positive_ions(std::move(positive_ions)), negative_ions(std::move(negative_ions)), + ions(std::move(positive_ions)) { finish_construction(); } @@ -112,10 +113,9 @@ struct Component { /// @param name The species/name for this instance. /// @param options Component settings: options[name] are specific to this component /// @param solver Time-integration solver - static std::unique_ptr create(const std::string &type, // The type to create - const std::string &name, // The species/name for this instance - Options &options, // Component settings: options[name] are specific to this component - Solver *solver); // Time integration solver + static std::unique_ptr create(const std::string& type, + const std::string& name, Options& options, + Solver* solver); /// Tell the component the name of all species in the simulation, by type. It /// will use this information to substitute the following placeholders in @@ -451,14 +451,15 @@ template Options& add(Options& option, T value) { if (!option.isSet()) { return set(option, value); - } else { - try { - return set(option, value + bout::utils::variantStaticCastOrThrow(option.value)); - } catch (const std::bad_cast &e) { - // Convert to a more useful error message - throw BoutException("Could not convert {:s} to type {:s}", - option.str(), typeid(T).name()); - } + } + try { + return set(option, value + + bout::utils::variantStaticCastOrThrow( + option.value)); + } catch (const std::bad_cast& e) { + // Convert to a more useful error message + throw BoutException("Could not convert {:s} to type {:s}", option.str(), + typeid(T).name()); } } @@ -477,14 +478,15 @@ template Options& subtract(Options& option, T value) { if (!option.isSet()) { return set(option, -value); - } else { - try { - return set(option, bout::utils::variantStaticCastOrThrow(option.value) - value); - } catch (const std::bad_cast &e) { - // Convert to a more useful error message - throw BoutException("Could not convert {:s} to type {:s}", - option.str(), typeid(T).name()); - } + } + try { + return set(option, + bout::utils::variantStaticCastOrThrow(option.value) + - value); + } catch (const std::bad_cast& e) { + // Convert to a more useful error message + throw BoutException("Could not convert {:s} to type {:s}", option.str(), + typeid(T).name()); } } diff --git a/include/component_scheduler.hxx b/include/component_scheduler.hxx index a67c27b06..c02cbd2ad 100644 --- a/include/component_scheduler.hxx +++ b/include/component_scheduler.hxx @@ -1,10 +1,11 @@ #pragma once - #ifndef COMPONENT_SCHEDULER_H #define COMPONENT_SCHEDULER_H -#include #include +#include +#include +#include #include #include diff --git a/src/component_scheduler.cxx b/src/component_scheduler.cxx index e97581a5b..3eb83baf7 100644 --- a/src/component_scheduler.cxx +++ b/src/component_scheduler.cxx @@ -1,15 +1,23 @@ #include +#include +#include +#include #include +#include #include +#include #include #include +#include #include #include // for trim, strsplit +#include #include #include "../include/component.hxx" #include "../include/component_scheduler.hxx" +#include "../include/permissions.hxx" const std::set ComponentScheduler::predeclared_variables = { "time", "linear", "units:inv_meters_cubed", "units:eV", "units:Tesla", @@ -46,7 +54,8 @@ void topological_sort(const std::vector>& dependencies, size_t /// separated by colons in the path. std::set getParents(const std::string& name) { std::set result; - size_t start = 0, position = name.find(":", start); + size_t start = 0; + size_t position = name.find(":", start); while (position != std::string::npos) { result.insert(name.substr(0, position)); start = position + 1; @@ -77,7 +86,8 @@ getVariableHierarchy(const std::vector>& components) // Build up a set of all section/variable names which are definitely // read/written by components, and the sections which they imply // exist - std::set unconditional_names, unconditional_sections; + std::set unconditional_names; + std::set unconditional_sections; for (const auto& component : components) { const Permissions& permissions = component->getPermissions(); for (const auto& [varname, _] : @@ -191,8 +201,9 @@ setReadDependencies(const std::vector>& components, // them for (const auto& [name, regions] : permissions.getVariablesWithPermission(permission)) { - if (ComponentScheduler::predeclared_variables.count(name) > 0) + if (ComponentScheduler::predeclared_variables.count(name) > 0) { continue; + } for (const auto& sub_name : expandVariableName(hierarchy, permissions, name)) { for (const auto& [region, _] : Permissions::fundamental_regions) { if ((regions & region) == region) { @@ -211,25 +222,21 @@ setReadDependencies(const std::vector>& components, return missing; } -/// Consumes a list of components and returns a new one that has been -/// topolgically sorted to ensure variables are written and read in -/// the right order. -std::vector> -sortComponents(std::vector>&& components) { +/// Topologically sorts the list of components to ensure variables are +/// written and read in the right order. +void sortComponents(std::vector>& components) { // Map between variable/section names specified by component // permissions and the variables they contain. In the case of // sections this is all variables within the section and any // sub-sections. Non-section viarables map to themselves. - std::map> variable_hierarchy = + const std::map> variable_hierarchy = getVariableHierarchy(components); // Get information on which components write each variable - std::map> nonfinal_writes = getPermissionComponentMap( - components, variable_hierarchy, - PermissionTypes::Write), - final_writes = getPermissionComponentMap( - components, variable_hierarchy, - PermissionTypes::Final); + std::map> nonfinal_writes = + getPermissionComponentMap(components, variable_hierarchy, PermissionTypes::Write); + std::map> final_writes = + getPermissionComponentMap(components, variable_hierarchy, PermissionTypes::Final); // Object mapping between components (reprsented by the index of // that component in the `components` argument) and the components @@ -244,7 +251,7 @@ sortComponents(std::vector>&& components) { throw BoutException( "Multiple components have permission to make final write to variable {}", var); } - for (size_t i : comp_indices) { + for (const size_t i : comp_indices) { const auto item = nonfinal_writes.find(var); if (item != nonfinal_writes.end()) { // Note that calling merge actually removes the items from the @@ -283,8 +290,8 @@ sortComponents(std::vector>&& components) { PermissionTypes::ReadIfSet, component_dependencies); // Create ancillary variables for sorting process - std::vector processing(components.size(), false), - processed(components.size(), false); + std::vector processing(components.size(), false); + std::vector processed(components.size(), false); std::vector order; // Perform the sort @@ -300,7 +307,7 @@ sortComponents(std::vector>&& components) { std::swap(result[i], components[order[i]]); } - return result; + components = std::move(result); } ComponentScheduler::ComponentScheduler(Options &scheduler_options, @@ -372,7 +379,7 @@ ComponentScheduler::ComponentScheduler(Options &scheduler_options, component->declareAllSpecies(species); } - components = sortComponents(std::move(components)); + sortComponents(components); } std::unique_ptr ComponentScheduler::create(Options &scheduler_options, From 8fc17e0cad037999b97a31cdf53d6558cdd859d8 Mon Sep 17 00:00:00 2001 From: Chris MacMackin Date: Mon, 5 Jan 2026 16:50:21 +0000 Subject: [PATCH 12/12] Add missed variable permission --- src/braginskii_conduction.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/braginskii_conduction.cxx b/src/braginskii_conduction.cxx index f38d08ed2..ff4f9a73a 100644 --- a/src/braginskii_conduction.cxx +++ b/src/braginskii_conduction.cxx @@ -27,7 +27,7 @@ using bout::globals::mesh; BraginskiiConduction::BraginskiiConduction(const std::string&, Options& alloptions, Solver*) - : Component({readOnly("species:{sp}:{input_vars}"), + : Component({readIfSet("fields:Apar_flutter"), readOnly("species:{sp}:{input_vars}"), readWrite("species:{sp}:{output_vars}")}) { AUTO_TRACE();