Allow components to call for creation of sub components

CMakeLists.txt

@@ -125,6 +125,7 @@ set(HERMES_SOURCES
     include/classical_diffusion.hxx
     include/cx_reaction.hxx
     include/binormal_stpm.hxx
+    include/braginskii_closure.hxx
     include/braginskii_collisions.hxx
     include/braginskii_conduction.hxx
     include/braginskii_electron_viscosity.hxx

docs/sphinx/closure.rst

@@ -40,6 +40,11 @@ is captured in the :ref:`sec-neutral_parallel_diffusion` top-level component, wh
 both parallel Braginskii transport and perpendicular pressure-diffusion for 2D/3D
 are captured in the :ref:`sec-neutral_mixed` species-level component.
 
+A user can automatically activate all of these components at once
+using the `BraginskiiClosure` component.
+
+.. doxygenclass:: BraginskiiClosure
+   :members:
 
 
 Collision frequency selection

docs/sphinx/developer.rst

@@ -617,6 +617,17 @@ species and density of ions::
 See the documentation for `Component::declareAllSpecies` for a list of
 all substitutions that will be performed.
 
+Components may request the creation of additional components, upon
+which they depend. This is done by overriding the
+`Component::additionalComponents` method, which returns a list of
+`ComponentInformation` structs that specify the names and types of
+components required. These extra components will use the default
+settings for components of this type, unless other values are
+specified in a section of the input file with the component name.
+
+.. doxygenstruct:: ComponentInformation
+   :members:
+
 
 Component scheduler
 ~~~~~~~~~~~~~~~~~~~
@@ -631,7 +642,10 @@ and then in `Hermes::rhs` the components are run by a call::
 
 The call to `ComponentScheduler::create` treats the "components"
 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
+the scheduler looks up the options under the section of that name. If
+any of the listed components request further components
+(via the `Component::AdditionalComponents` method) then these will be
+created too. 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

hermes-3.cxx

@@ -29,6 +29,7 @@
 #include "include/amjuel_data.hxx"
 #include "include/anomalous_diffusion.hxx"
 #include "include/binormal_stpm.hxx"
+#include "include/braginskii_closure.hxx"
 #include "include/braginskii_collisions.hxx"
 #include "include/braginskii_conduction.hxx"
 #include "include/braginskii_electron_viscosity.hxx"

include/braginskii_closure.hxx

@@ -0,0 +1,70 @@
+#pragma once
+#ifndef BRAGINSKII_H
+#define BRAGINSKII_H
+
+#include <bout/field3d.hxx>
+
+#include "component.hxx"
+
+/// Meta-component to set up all components necessary for the
+/// Braginskii closure: `braginskii_collisions`,
+/// `braginskii_friction`, `braginskii_heat_exchange`,
+/// `braginskii_conduction`, `braginskii_electron_viscosity`,
+/// `braginskii_ion_viscosity`, and `braginskii_thermal_force`. Each
+/// of these components will have the same name as its type.
+class BraginskiiClosure : public NamedComponent<BraginskiiClosure> {
+public:
+  /// @param alloptions Settings, which may include
+  ///    - <name>
+  ///      - electron_viscosity  : bool    Include electron viscosity? (default: true)
+  ///      - ion_viscosity  : bool    Include ion viscosity? (default: true)
+  ///      - thermal_force  : bool    Inlucde thermal force between species? (default: true)
+  BraginskiiClosure(std::string name, Options& alloptions, Solver*)
+      : NamedComponent(name, {}) {
+    Options& options = alloptions[name];
+    electron_viscosity = options["electron_viscosity"]
+                             .doc("Include electron viscosity terms?")
+                             .withDefault<bool>(true);
+    ion_viscosity = options["ion_viscosity"]
+                        .doc("Include ion viscosity terms?")
+                        .withDefault<bool>(true);
+    thermal_force = options["thermal_force"]
+                        .doc("Include thermal force terms?")
+                        .withDefault<bool>(true);
+  }
+
+  virtual std::vector<ComponentInformation> additionalComponents() override {
+    std::vector<ComponentInformation> result = {
+        {"braginskii_collisions", "braginskii_collisions"},
+        {"braginskii_friction", "braginskii_friction"},
+        {"braginskii_heat_exchange", "braginskii_heat_exchange"},
+        {"braginskii_conduction", "braginskii_conduction"}};
+    if (electron_viscosity) {
+      result.emplace_back("braginskii_electron_viscosity",
+                          "braginskii_electron_viscosity");
+    }
+    if (ion_viscosity) {
+      result.emplace_back("braginskii_ion_viscosity", "braginskii_ion_viscosity");
+    }
+    if (thermal_force) {
+      result.emplace_back("braginskii_thermal_force", "braginskii_thermal_force");
+    }
+    return result;
+  }
+
+  static constexpr auto type = "braginskii_closure";
+
+private:
+  bool electron_viscosity; /// Whether to include electron viscosity terms
+  bool ion_viscosity;      /// Whether to include ion viscosity terms
+  bool thermal_force;      /// Whether to include thermal force terms
+
+  /// Empty transform; all the work actually happens in the subcomponents.
+  void transform_impl(GuardedOptions&) override {}
+};
+
+namespace {
+RegisterComponent<BraginskiiClosure> registercomponentbraginskiiclosure;
+}
+
+#endif // BRAGINSKII_H

include/component.hxx

@@ -76,6 +76,40 @@ private:
   }
 };
 
+/// Lightweight type to store information used to build a component.
+struct ComponentInformation {
+  std::string name;
+  std::string type;
+
+  ComponentInformation() {};
+
+  ComponentInformation(const std::string& name_, const std::string& type_)
+      : name(name_), type(type_) {}
+
+  ComponentInformation(std::string&& name_, std::string&& type_)
+      : name(std::move(name_)), type(std::move(type_)) {}
+
+  bool operator<(const ComponentInformation& other) const {
+    return std::pair(name, type) < std::pair(other.name, other.type);
+  }
+
+  bool operator==(const ComponentInformation& other) const {
+    return std::pair(name, type) == std::pair(other.name, other.type);
+  }
+};
+
+/// Format `ComponentInformation` to string. Format string specification is the
+/// same as when formatting a string.
+/// See https://fmt.dev/12.0/syntax/#format-specification-mini-language.
+///
+/// TODO: provide custom formatting to configure exactly how the
+/// component name and type are displayed.
+template <>
+struct fmt::formatter<ComponentInformation> : formatter<std::string> {
+  auto format(const ComponentInformation& ci, format_context& ctx) const
+      -> format_context::iterator;
+};
+
 /// Interface for a component of a simulation model
 ///
 /// The constructor of derived types should have signature
@@ -92,6 +126,11 @@ struct Component {
 
   virtual ~Component() {}
 
+  /// Return a list of names/types of other components needed by this
+  /// component. All configurations for these components will take the
+  /// default value, unless set in the input file.
+  virtual std::vector<ComponentInformation> additionalComponents() { return {}; }
+
   /// Modify the given simulation state. This method will wrap the
   /// state in a GuardedOptions object and pass that to the private
   /// implementation of transform provided by each component.

include/electron_force_balance.hxx

@@ -20,16 +20,14 @@
 ///
 struct ElectronForceBalance : public NamedComponent<ElectronForceBalance> {
   ElectronForceBalance(std::string name, Options& alloptions, Solver*)
-      : NamedComponent(name,
-                       {readOnly("species:e:pressure"),
-                        readOnly("species:e:density", Regions::Interior),
-                        readOnly("species:e:charge"),
-                        // FIXME: Only writes if already exists
-                        readWrite("species:e:momentum_source"),
-                        readIfSet("species:{non_electrons}:density", Regions::Interior),
-                        readIfSet("species:{non_electrons}:charge"),
-                        // FIXME: Only written if density and charge have been set.
-                        readWrite("species:{non_electrons}:momentum_source")}) {
+      : NamedComponent(
+            name, {readOnly("species:e:pressure"),
+                   readOnly("species:e:density", Regions::Interior),
+                   readOnly("species:e:charge"), readIfSet("species:e:momentum_source"),
+                   readIfSet("species:{non_electrons}:density", Regions::Interior),
+                   readIfSet("species:{non_electrons}:charge"),
+                   // FIXME: Only written if density and charge have been set.
+                   readWrite("species:{non_electrons}:momentum_source")}) {
     auto& options = alloptions[name];
     diagnose = options["diagnose"]
                    .doc("Save additional output diagnostics")

src/braginskii_collisions.cxx

@@ -25,7 +25,8 @@
 BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& alloptions,
                                            Solver*)
     : NamedComponent(name,
-                     {readOnly("species:{non_electrons}:density", Regions::Interior),
+                     {readOnly("species:{all_species}:density", Regions::Interior),
+                      readOnly("species:{electrons}:temperature", Regions::Interior),
                       readIfSet("species:{non_electrons}:charge"),
                       readIfSet("species:{negative_ions}:temperature", Regions::Interior),
                       readOnly("species:{all_species}:AA")}) {
@@ -65,8 +66,6 @@ BraginskiiCollisions::BraginskiiCollisions(const std::string& name, Options& all
   diagnose =
       options["diagnose"].doc("Output additional diagnostics?").withDefault<bool>(false);
 
-  setPermissions(readOnly("species:{electrons}:temperature", Regions::Interior));
-  setPermissions(readOnly("species:{electrons}:density", Regions::Interior));
   if (electron_electron) {
     setPermissions(readWrite(
         "species:{electrons}:collision_frequencies:{electrons}_{electrons2}_coll"));

src/braginskii_ion_viscosity.cxx

@@ -28,6 +28,8 @@
 #include "../include/div_ops.hxx"
 #include "../include/hermes_utils.hxx"
 
+// FIXME: Why does it matter if upstream_density_feedback executed after this? Is there some wierd interaction of the readIfSet variables?
+
 using bout::globals::mesh;
 
 BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
@@ -36,10 +38,11 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
           name,
           {
               readIfSet("species:{non_electrons}:pressure"),
-              readIfSet("species:{non_electrons}:temperature"),
-              readIfSet("species:{non_electrons}:density"),
-              readIfSet("species:{non_electrons}:velocity"),
               readIfSet("species:{non_electrons}:charge"),
+              // FIXME: The rest of these (except DivJextra) only
+              // apply for species in which the above are set
+              readOnly("species:{non_electrons}:temperature"),
+              readOnly("species:{non_electrons}:density"),
               readIfSet("species:{non_electrons}:collision_frequencies:{coll_type}"),
               readWrite("species:{non_electrons}:momentum_source"),
               readWrite("species:{non_electrons}:energy_source"),
@@ -121,11 +124,20 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
     Curlb_B.z *= SQ(Lnorm);
 
     Curlb_B *= 2. / coord->Bxy;
+
+    setPermissions(readOnly("fields:phi"));
+  }
+
+  if (parallel) {
+    setPermissions(readOnly("species:{non_electrons}:velocity"));
   }
+
   if (bounce_frequency) {
     const Options& units = alloptions["units"];
     const BoutReal Lnorm = units["meters"];
     bounce_frequency_R /= Lnorm;
+
+    setPermissions(readOnly("species:{non_electrons}:pressure"));
   }
 
   std::vector<std::string> coll_types;
@@ -135,9 +147,6 @@ BraginskiiIonViscosity::BraginskiiIonViscosity(const std::string& name,
     coll_types.push_back("{non_electrons}_{all_species}_coll");
     coll_types.push_back("{non_electrons}_{all_species}_cx");
   }
-  if (perpendicular) {
-    setPermissions(readOnly("fields:phi"));
-  }
   substitutePermissions("coll_type", coll_types);
 }
 

src/component.cxx

@@ -9,6 +9,12 @@
 #include "../include/guarded_options.hxx"
 #include "../include/permissions.hxx"
 
+auto fmt::formatter<ComponentInformation>::format(const ComponentInformation& ci,
+                                                  format_context& ctx) const
+    -> format_context::iterator {
+  return formatter<std::string>::format(fmt::format("{} ({})", ci.name, ci.type), ctx);
+}
+
 std::unique_ptr<Component> Component::create(const std::string& type,
                                              const std::string& name, Options& alloptions,
                                              Solver* solver) {

src/component_scheduler.cxx

@@ -11,6 +11,7 @@
 #include <bout/bout_types.hxx>
 #include <bout/boutexception.hxx>
 #include <bout/options.hxx>
+#include <bout/output.hxx>
 #include <bout/utils.hxx> // for trim, strsplit
 #include <fmt/format.h>
 #include <fmt/ranges.h>
@@ -308,6 +309,15 @@ setReadDependencies(const std::vector<std::unique_ptr<Component>>& components,
   return missing;
 }
 
+void printComponents(const std::vector<std::unique_ptr<Component>>& components) {
+  if (!components.empty()) {
+    output_info << "Components will be executed in the following order:\n";
+  }
+  for (const auto& comp : components) {
+    output_info << fmt::format("\t{}\n", *comp);
+  }
+}
+
 /// Topologically sorts the list of components to ensure variables are
 /// written and read in the right order.
 ///
@@ -425,6 +435,7 @@ void sortComponents(std::vector<std::unique_ptr<Component>>& components) {
   }
 
   components = std::move(result);
+  printComponents(components);
 }
 } // namespace
 
@@ -434,6 +445,12 @@ ComponentScheduler::ComponentScheduler(Options& scheduler_options,
   const std::string component_names = scheduler_options["components"]
                                           .doc("Components in order of execution")
                                           .as<std::string>();
+  const bool autosort = scheduler_options["autosort"]
+                            .doc("Perform a topological sort to ensure components "
+                                 "executed in the right order?")
+                            .withDefault<bool>(true);
+
+  std::list<ComponentInformation> required_components;
 
   std::vector<std::string> electrons;
   std::vector<std::string> neutrals;
@@ -483,18 +500,52 @@ ComponentScheduler::ComponentScheduler(Options& scheduler_options,
         continue;
       }
 
-      components.push_back(
-          Component::create(type_trimmed, name_trimmed, component_options, solver));
+      required_components.emplace_back(name_trimmed, type_trimmed);
     }
   }
 
+  // Use sets for efficient lookup of whether a component is already
+  // in the queue to be created or already has been created.
+  //
+  // We use a vector to decide the order in which to create
+  // components, to keep this close to what is in the file.
+  std::set<ComponentInformation> unbuilt_components(required_components.begin(),
+                                                    required_components.end());
+  std::set<ComponentInformation> created_components;
+  // We use a vector to keep track of the actual order in which
+  // components are created
+  std::vector<ComponentInformation> component_order;
+
+  while (required_components.size() > 0) {
+    const ComponentInformation component = required_components.front();
+    required_components.pop_front();
+    auto comp =
+        Component::create(component.type, component.name, component_options, solver);
+    std::vector<ComponentInformation> sub_components = comp->additionalComponents();
+    for (auto sub_comp = sub_components.rbegin(); sub_comp != sub_components.rend();
+         ++sub_comp) {
+      if (unbuilt_components.count(*sub_comp) == 0
+          and created_components.count(*sub_comp) == 0) {
+        required_components.push_front(*sub_comp);
+        unbuilt_components.insert(*sub_comp);
+      }
+    }
+    component_order.push_back(component);
+    components.push_back(std::move(comp));
+    created_components.insert(unbuilt_components.extract(component));
+  }
+
   const SpeciesInformation species(electrons, neutrals, positive_ions, negative_ions);
 
   for (auto& component : components) {
     component->declareAllSpecies(species);
   }
 
-  ::sortComponents(components);
+  if (autosort) {
+    ::sortComponents(components);
+  } else {
+    printComponents(components);
+  }
 }
 
 std::unique_ptr<ComponentScheduler> ComponentScheduler::create(Options& scheduler_options,