Perform topological sort to order components

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)

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
+   source_prefactor = Pe:source_prefactor

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
 ------

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:

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:
 

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
 

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
 

include/component.hxx

@@ -1,17 +1,7 @@
 #pragma once
-
 #ifndef HERMES_COMPONENT_H
 #define HERMES_COMPONENT_H
 
-#include <bout/assert.hxx>
-#include <bout/bout_types.hxx>
-#include <bout/boutexception.hxx>
-#include <bout/field2d.hxx>
-#include <bout/field3d.hxx>
-#include <bout/generic_factory.hxx>
-#include <bout/options.hxx>
-#include <bout/unused.hxx>
-
 #include <cmath>
 #include <initializer_list>
 #include <map>
@@ -22,6 +12,15 @@
 #include <utility>
 #include <vector>
 
+#include <bout/assert.hxx>
+#include <bout/bout_types.hxx>
+#include <bout/boutexception.hxx>
+#include <bout/field2d.hxx>
+#include <bout/field3d.hxx>
+#include <bout/generic_factory.hxx>
+#include <bout/options.hxx>
+#include <bout/unused.hxx>
+
 #include "guarded_options.hxx"
 #include "permissions.hxx"
 #include "hermes_utils.hxx"
@@ -34,8 +33,10 @@ struct SpeciesInformation {
   SpeciesInformation(const std::vector<std::string>& electrons,
                      const std::vector<std::string>& neutrals,
                      const std::vector<std::string>& positive_ions,
-                     const std::vector<std::string> & negative_ions)
-    : electrons(electrons), neutrals(neutrals), positive_ions(positive_ions), negative_ions(negative_ions), ions(positive_ions) {
+                     const std::vector<std::string>& 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<Component> 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<Component> 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
@@ -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,
@@ -449,14 +451,15 @@ template<typename T>
 Options& add(Options& option, T value) {
   if (!option.isSet()) {
     return set(option, value);
-  } else {
-    try {
-      return set(option, value + bout::utils::variantStaticCastOrThrow<Options::ValueType, T>(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<Options::ValueType, T>(
+                               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());
   }
 }
 
@@ -475,14 +478,15 @@ template<typename T>
 Options& subtract(Options& option, T value) {
   if (!option.isSet()) {
     return set(option, -value);
-  } else {
-    try {
-      return set(option, bout::utils::variantStaticCastOrThrow<Options::ValueType, T>(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<Options::ValueType, T>(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());
   }
 }
 

include/component_scheduler.hxx

@@ -1,10 +1,11 @@
 #pragma once
-
 #ifndef COMPONENT_SCHEDULER_H
 #define COMPONENT_SCHEDULER_H
 
-#include <vector>
 #include <memory>
+#include <set>
+#include <string>
+#include <vector>
 
 #include <bout/bout_types.hxx>
 #include <bout/options.hxx>
@@ -57,6 +58,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<std::string> predeclared_variables;
+
 private:
   /// The components to be executed in order
   std::vector<std::unique_ptr<Component>> components;