Add hierarchical PhasorDynamics system model with a unit test

GridKit/Model/PhasorDynamics/INPUT_FORMAT.md

@@ -35,9 +35,9 @@ also be encoded as [MessagePack](https://msgpack.org).
 ## Format
 
 The root object contains `header`, `params`, `buses`, and `devices`, and may
-also contain `signals` and `monitors`. `header` contains information about
-the case, `params` contains system-wide parameters, `buses` is an array of
-electrical buses, and `devices` is an array of system devices.
+also contain `signals`, `monitors` and `subsystems`. `header` contains
+information about the case, `params` contains system-wide parameters, `buses`
+is an array of electrical buses, and `devices` is an array of system devices.
 
 ### Header
 
@@ -161,6 +161,37 @@ are specified:
   [Ieeest](Stabilizer/IEEEST/README.md) | the IEEEST stabilizer model
   [ConstantSignalSource](SignalSource/README.md) | Constant complex signal source
 
+### Subsystems
+
+Contained in the optional `subsystems` key is an array of objects, each of
+which describes a hierarchical subsystem with its own `buses` and `devices`,
+and optional `signals` and `subsystems`. Subsystems use the root's `params`.
+Each subsystem also requires:
+
+Key | Description
+-- | --
+`id` | Subsystem ID (string), unique across the whole hierarchy
+`level` | Hierarchy level (integer): the root system is level 0 (implied, not given in the file), its subsystems level 1, their subsystems level 2, and so on
+
+The `level` must match the depth of the subsystem. It is redundant in this
+hierarchical layout, but makes the hierarchy explicit, so that a flat layout
+can be supported later.
+That flat layout will also need a `subsystem` tag (the ID of the containing
+subsystem) on buses, devices and signals; this tag is not supported yet.
+
+The input must already be partitioned:
+
+-   A system with subsystems contains only buses and subsystems: its
+    `devices` array must be empty.
+-   Devices may use buses of their own system and of its ancestors, but not
+    of sibling subsystems.
+-   Devices may only use signals of their own system.
+-   Bus numbers and signal IDs must be unique across the whole hierarchy.
+
+Not all of these are checked yet: duplicate subsystem IDs, duplicate signal
+IDs across systems, and bus numbers repeated in sibling subsystems are not
+detected.
+
 ## Example File for a 2-Bus System
 
 ```json

GridKit/Model/PhasorDynamics/SystemModel.hpp

@@ -3,6 +3,7 @@
 #include <map>
 #include <memory>
 #include <string>
+#include <vector>
 
 #include <GridKit/Model/PhasorDynamics/Component.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeSet.hpp>
@@ -36,6 +37,14 @@ namespace GridKit
      * Model::Evaluator for the system model. This is still work in
      * progress and code is not optimized.
      *
+     * Subsystems own their buses and devices and may connect to the parent's
+     * buses (external buses). A system with subsystems contains only buses
+     * and subsystems. Only the top-level system allocates storage.
+     *
+     * evaluateJacobian() on the top-level system assembles the system Jacobian.
+     * On a subsystem, it assembles the contributions of the subsystem's buses
+     * and devices in global indices, of size offset() + size().
+     *
      * @todo Address thread safety for the system model methods.
      *
      */
@@ -74,6 +83,9 @@ namespace GridKit
 
       SystemModel(const SystemModelData<RealT, IdxT>& data);
 
+      SystemModel(const SystemModelData<RealT, IdxT>& data,
+                  const std::vector<BusT*>&           external_buses);
+
       virtual ~SystemModel();
 
       int setGridKitComponentID(IdxT component_id) override;
@@ -100,6 +112,7 @@ namespace GridKit
 
       void addBus(BusT* bus);
       void addComponent(ComponentT* component);
+      void addSubsystem(SystemModel* subsystem);
 
       void setSystemBase(RealT freq_system_base, RealT va_system_base);
 
@@ -108,17 +121,34 @@ namespace GridKit
       ComponentT*            getComponent(IdxT gridkit_component_id);
       Branch<ScalarT, IdxT>* getBranch(const std::string& id);
 
+      IdxT         numSubsystems() const;
+      SystemModel* getSubsystem(IdxT index);
+      IdxT         offset() const;
+
+      const std::vector<BusT*>& getExternalBuses() const;
+
     private:
-      std::vector<BusT*>       buses_;
-      SignalNodeSetT           signal_nodes_;
-      std::vector<ComponentT*> components_;
+      BusT* findBus(IdxT bus_id);
+      IdxT  computeSize();
+      void  bindChildren(IdxT global_offset);
+      void  addMonitorsTo(MonitorT& monitor) const;
+      void  evaluateCooJacobians();
 
-      std::map<IdxT, IdxT> gridkit_bus_indices_; ///< Map between gridkit_bus_id and bus_id
+      std::vector<BusT*>        buses_;
+      SignalNodeSetT            signal_nodes_;
+      std::vector<ComponentT*>  components_;
+      std::vector<SystemModel*> subsystems_;
 
       /// Map between branch id and gridkit_component_id
       /// @todo remove this map and use branch IDs to access individual branches.
       std::map<std::string, IdxT> gridkit_branch_indices_;
 
+      std::map<IdxT, IdxT>  gridkit_bus_indices_;    ///< Map between gridkit_bus_id and bus_id
+      std::map<IdxT, BusT*> external_buses_;         ///< Parent's buses by bus_id, not owned
+      std::vector<BusT*>    ordered_external_buses_; ///< External buses in construction order, not owned
+
+      IdxT offset_{0};           ///< Global index of this system's first variable
+      bool is_subsystem_{false}; ///< Allocated and bound by a parent system
       bool owns_components_{false};
 
       /// Variable monitor

GridKit/Model/PhasorDynamics/SystemModelData.hpp

@@ -124,6 +124,16 @@ namespace GridKit
 
       /// Monitor sink specs
       std::vector<MonitorSinkSpec> monitor_sink;
+
+      /// Subsystem ID, unique in the hierarchy; empty for the top system
+      std::string id;
+
+      /// Hierarchy level: 0 for the top system, parent level + 1 for a subsystem
+      IdxT level{0};
+
+      /// Hierarchical subsystems; devices may use buses of this system and
+      /// its ancestors. System bases must match this system's.
+      std::vector<SystemModelData> subsystems;
     };
 
     /**

GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp

@@ -2,6 +2,7 @@
 
 #include <sstream>
 #include <stdexcept>
+#include <utility>
 
 #include <nlohmann/json.hpp>
 
@@ -19,6 +20,9 @@ namespace GridKit
     using Log           = ::GridKit::Utilities::Logger;
     using MonitorFormat = ::GridKit::Model::VariableMonitorFormat;
 
+    template <typename RealT, typename IdxT>
+    void parseSystemContents(const json& j, SystemModelData<RealT, IdxT>& sm);
+
     /// JSON parser function implementation for the `SystemModelData` type
     ///
     /// See the `README.md` in `GridKit/Model/PhasorDynamics` for more information
@@ -77,6 +81,13 @@ namespace GridKit
         }
       }
 
+      parseSystemContents(j, sm);
+    }
+
+    /// Parse buses, signals, devices and subsystems of a system or subsystem
+    template <typename RealT, typename IdxT>
+    void parseSystemContents(const json& j, SystemModelData<RealT, IdxT>& sm)
+    {
       /// Gets all electrical buses
       j.at("buses").get_to(sm.bus);
 
@@ -206,6 +217,45 @@ namespace GridKit
           throw std::runtime_error("JSON parser failed");
         }
       }
+
+      /// Gets hierarchical subsystems, which use this system's bases
+      if (j.contains("subsystems"))
+      {
+        if (!j.at("subsystems").empty() && !j.at("devices").empty())
+        {
+          Log::error() << "\n\tA system with subsystems may only contain buses and subsystems; "
+                       << "move its devices into a subsystem.\n";
+          throw std::runtime_error("JSON parser failed");
+        }
+
+        for (auto& raw_subsystem : j.at("subsystems"))
+        {
+          SystemModelData<RealT, IdxT> subsystem;
+          subsystem.freq_base = sm.freq_base;
+          subsystem.va_base   = sm.va_base;
+          subsystem.level     = sm.level + 1;
+
+          if (!raw_subsystem.contains("id") || !raw_subsystem.at("id").is_string())
+          {
+            Log::error() << "\n\tEach subsystem requires a string \"id\".\n";
+            throw std::runtime_error("JSON parser failed");
+          }
+          raw_subsystem.at("id").get_to(subsystem.id);
+
+          // The level is redundant here; it allows a flat layout later.
+          if (!raw_subsystem.contains("level")
+              || !raw_subsystem.at("level").is_number_integer()
+              || raw_subsystem.at("level").get<long long>() != static_cast<long long>(subsystem.level))
+          {
+            Log::error() << "\n\tSubsystem \"" << subsystem.id << "\" requires \"level\": "
+                         << subsystem.level << " (parent level + 1).\n";
+            throw std::runtime_error("JSON parser failed");
+          }
+
+          parseSystemContents(raw_subsystem, subsystem);
+          sm.subsystems.push_back(std::move(subsystem));
+        }
+      }
     }
   } // namespace PhasorDynamics
 } // namespace GridKit

GridKit/Model/PhasorDynamics/SystemModelDependencyTracking.cpp

@@ -25,6 +25,14 @@ namespace GridKit
     template <typename scalar_type, typename index_type>
     int SystemModel<scalar_type, index_type>::evaluateJacobian()
     {
+      // The top-level system builds the Jacobian from its whole residual.
+      // @todo Subsystem Jacobian from the subsystem's own residual.
+      if (is_subsystem_)
+      {
+        Log::error() << "Subsystem Jacobians are only available with Enzyme\n";
+        return -1;
+      }
+
       this->constructCsr();
 
       // Log::misc() << "System DependencyTracking Jacobian\n";