Remove external variable copy for PowerElectronics

CHANGELOG.md

@@ -72,6 +72,7 @@
 - Added new `Rosenbrock` integrator.
 - Clarified naming conventions for macros.
 - Added `dt_fixed`, `rel_tol`, and `abs_tol` options to phasor dynamics solver JSON files and renamed the `dt` option to `dt_monitor`.
+- Remove unnecessary data copying while evaluating `PowerElectronics` models, speeding up large simulations by up to 3x
 
 ## v0.1
 

GridKit/Model/PowerElectronics/Bus/GroundedBus.hpp

@@ -9,38 +9,29 @@ namespace GridKit
     template <typename ScalarT, typename IdxT>
     class GroundedBus : public NodeBase<ScalarT, IdxT>
     {
-      using NodeBase<ScalarT, IdxT>::y;
+      using NodeBase<ScalarT, IdxT>::y_ext_;
+      using NodeBase<ScalarT, IdxT>::yp_ext_;
+      using NodeBase<ScalarT, IdxT>::f_ext_;
 
     public:
       GroundedBus(ScalarT voltage)
         : NodeBase<ScalarT, IdxT>(0, 1), voltage_(voltage)
       {
       }
 
-      int initialize() final
-      {
-        if (int err_code = NodeBase<ScalarT, IdxT>::initialize())
-          return err_code;
-
-        auto* y_data = y().getData();
-        y_data[0]    = voltage_;
-        y().setDataUpdated();
-
-        return 0;
-      }
-
       int allocate() final
       {
         if (int err_code = NodeBase<ScalarT, IdxT>::allocate())
           return err_code;
 
-        this->setExternalConnectionNodes(0, INVALID_INDEX<IdxT>);
+        this->setExternalConnectionNodes(0, ExternalConnection<ScalarT, IdxT>{.y_ = &voltage_, .yp_ = &dummy_, .f_ = &dummy_, .idx_ = INVALID_INDEX<IdxT>});
 
         return 0;
       }
 
     private:
       ScalarT voltage_;
+      ScalarT dummy_;
     };
   } // namespace PowerElectronics
 } // namespace GridKit

GridKit/Model/PowerElectronics/CMakeLists.txt

@@ -28,4 +28,5 @@ install(
         CircuitGraph.hpp
         SystemModelPowerElectronics.hpp
         NodeBase.hpp
+        ExternalConnection.hpp
   DESTINATION include/GridKit/Model/PowerElectronics)

GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp

@@ -59,22 +59,17 @@ namespace GridKit
   template <class ScalarT, typename IdxT>
   int Capacitor<ScalarT, IdxT>::evaluateInternalResidual()
   {
-    const auto* y = y_.getData();
-
-    f_int_[0] = -C_ * yp_int_[0] + y[0] - y[1] - y_int_[0];
+    f_int_[0] = -C_ * yp_int_[0] + *y_ext_[0] - *y_ext_[1] - y_int_[0];
     return 0;
   }
 
   template <class ScalarT, typename IdxT>
   int Capacitor<ScalarT, IdxT>::evaluateExternalResidual()
   {
-    auto* f = f_.getData();
-
     // input
-    f[0] = C_ * yp_int_[0];
+    *f_ext_[0] += C_ * yp_int_[0];
     // output
-    f[1] = -C_ * yp_int_[0];
-    f_.setDataUpdated();
+    *f_ext_[1] += -C_ * yp_int_[0];
     return 0;
   }
 

GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp

@@ -25,12 +25,12 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::nnz_;
     using CircuitComponent<ScalarT, IdxT>::time_;
     using CircuitComponent<ScalarT, IdxT>::alpha_;
-    using CircuitComponent<ScalarT, IdxT>::y_;
+    using CircuitComponent<ScalarT, IdxT>::y_ext_;
     using CircuitComponent<ScalarT, IdxT>::y_int_;
-    using CircuitComponent<ScalarT, IdxT>::yp_;
+    using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
-    using CircuitComponent<ScalarT, IdxT>::f_;
+    using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
     using CircuitComponent<ScalarT, IdxT>::g_;
     using CircuitComponent<ScalarT, IdxT>::yB_;

GridKit/Model/PowerElectronics/CircuitComponent.hpp

@@ -9,6 +9,7 @@
 
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
 #include <GridKit/Model/Evaluator.hpp>
+#include <GridKit/Model/PowerElectronics/ExternalConnection.hpp>
 
 namespace GridKit
 {
@@ -56,15 +57,41 @@ namespace GridKit
     }
 
     /**
-     * @brief Create the mappings from local to global indices
+     * @brief Create the mappings from local to global indices for an internal variable.
+     * Used for constructing system Jacobians \see connection_nodes_.
      *
-     * @param local_index
-     * @param global_index
-     * @return int
+     * @param local_index The index of the local variable
+     * @param global_index The index of the corresponding system variable.
+     *
+     * @pre `local_index` *must* be the index of an internal variable. Using this method for
+     * an external variable will not properly setup the data pointers for that variable.
+     */
+    int setInternalConnectionNodes(size_t local_index, size_t global_index)
+    {
+      assert(!extern_indices_.contains(local_index));
+      connection_nodes_[local_index] = global_index;
+      return 0;
+    }
+
+    /**
+     * @brief Create the mappings from local to global indices for an external variable.
+     * External variables need extra information than internal variables - their data
+     * pointers \ref y_ext_, \ref yp_ext_, and \ref f_ext_.
+     *
+     * @param local_index The index of the local variable
+     * @param connection The necessary connection information for the variable
+     *
+     * @pre `local_index` *must* be the index of an external variable. As of now, using this method
+     * to set information for a local variable will silently discard the unnecessary information, but
+     * this may change in the future.
      */
-    int setExternalConnectionNodes(IdxT local_index, IdxT global_index)
+    int setExternalConnectionNodes(size_t local_index, ExternalConnection<ScalarT, IdxT> connection)
     {
-      connection_nodes_[static_cast<size_t>(local_index)] = global_index;
+      assert(extern_indices_.contains(local_index));
+      y_ext_[local_index]            = connection.y_;
+      yp_ext_[local_index]           = connection.yp_;
+      f_ext_[local_index]            = connection.f_;
+      connection_nodes_[local_index] = connection.idx_;
       return 0;
     }
 
@@ -74,13 +101,21 @@ namespace GridKit
      * f(local_index) = global_index
      *
      * @param local_index index of local value in vector
-     * @return IdxT Index of the same value in the global vector
+     * @return size_t Index of the same value in the global vector
      */
-    IdxT getNodeConnection(IdxT local_index) const
+    IdxT getNodeConnection(size_t local_index) const
     {
       return connection_nodes_[local_index];
     }
 
+    int initialize() override
+    {
+      y_.setDataUpdated();
+      yp_.setDataUpdated();
+
+      return 0;
+    }
+
     /**
      * @brief Allocates all of the internal buffers for the component.
      * If a components needs a more specialized allocation (such as by having additional internal buffers),
@@ -97,7 +132,10 @@ namespace GridKit
       jacobian_coo_cols_   = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
       jacobian_coo_values_ = std::make_unique<RealT[]>(static_cast<size_t>(nnz_));
 
-      connection_nodes_ = std::make_unique<IdxT[]>(static_cast<size_t>(size_));
+      y_ext_            = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+      yp_ext_           = std::make_unique<const ScalarT*[]>(static_cast<size_t>(size_));
+      f_ext_            = std::make_unique<ScalarT*[]>(static_cast<size_t>(size_));
+      connection_nodes_ = std::make_unique<size_t[]>(static_cast<size_t>(size_));
 
       if (!allocated_)
       {
@@ -151,6 +189,8 @@ namespace GridKit
       if (int err_code = evaluateInternalResidual())
         return err_code;
 
+      f_.setDataUpdated();
+
       return evaluateExternalResidual();
     }
 
@@ -402,23 +442,57 @@ namespace GridKit
   protected:
     /**
      * @brief Allocate state and residual storage owned by this component.
+     *
+     * Most components do not need state and residual storages. The most notable exception
+     * is currently the system, so a separate flag is provided for the system.
+     * Systems still can't directly access \ref y_, \ref yp_, and \ref f_, so they need
+     * their corresponding \ref y_int_, \ref yp_int_, and \ref f_int_ set, since there isn't
+     * another system above them to set it.
+     *
+     * @todo This is a weird exception specifically for systems - and in a hierarchical setting
+     * will only be needed by the *topmost* system - subsystems shouldn't allocate and should have their
+     * internal pointers set by the system above them. Ideally we can remove this exception by having
+     * the integrator allocate these buffers instead of the system and set the internal pointers for the
+     * topmost system.
      */
-    void allocateVectors(IdxT n)
+    void allocateVectors(IdxT n, bool system = false)
     {
-      y_.resize(n);
-      yp_.resize(n);
-      f_.resize(n);
       abs_tol_.resize(n);
+
+      if (system)
+      {
+        y_.resize(n);
+        yp_.resize(n);
+        f_.resize(n);
+
+        y_int_  = y_.getData();
+        yp_int_ = yp_.getData();
+        f_int_  = f_.getData();
+      }
     }
 
-    size_t                  n_extern_;
-    size_t                  n_intern_;
-    std::set<IdxT>          extern_indices_;
-    ///@todo may want to replace the mapping of connection_nodes to Node objects instead of IdxT. Allows for container free setup
-    std::unique_ptr<IdxT[]> connection_nodes_;
+    /// Number of external variables in this component - ones which are referenced but not owned by this component.
+    size_t                    n_extern_;
+    /// Number of internal variables in this component - ones which are only referenced by this component.
+    size_t                    n_intern_;
+    /**
+     * @brief A set of variable indices which correspond to the external variables. Variables indices not in this set are internal.
+     *
+     * @invariant Must have a size of n_extern_. Each element must be in the range [0, `size_` - 1]. Not currently verified anywhere.
+     */
+    std::set<IdxT>            extern_indices_;
+    /**
+     * @brief A map from local variable indices to system (global) variable indices. Used for Jacobian construction in
+     * \ref PowerElectronicsModel::evaluateJacobian().
+     * @note If a variable does not map to a corresponding variable in the system (such as with reference nodes), a special
+     * sentinel value of \ref INVALID_INDEX is used. During Jacobian construction, such rows and columns will be pruned.
+     */
+    std::unique_ptr<size_t[]> connection_nodes_;
 
   protected:
+    /// The number of variables in this component. Should be equal to \ref n_extern_ plus \ref n_intern_ \see getSize()
     IdxT size_{0};
+    /// The number of nonzero elements in this component's Jacobian. \see nnz()
     IdxT nnz_{0};
     IdxT size_quad_{0};
     IdxT size_opt_{0};
@@ -438,11 +512,33 @@ namespace GridKit
     /// @brief A pointer to the internal residuals of this component
     ScalarT*       f_int_;
 
-    VectorT           y_;
-    VectorT           yp_;
+    /**
+     * An array of (input) pointers to state values for external variables.
+     * \note The size of this array is equal to \ref size_, allowing you to index it with the index
+     * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing
+     * and dereferencing the pointer in an internal variable index is undefined behavior.
+     * \see setExternalConnectionNodes()
+     */
+    std::unique_ptr<const ScalarT*[]> y_ext_;
+    /**
+     * An array of (input) pointers to derivative values for external variables.
+     * \note The size of this array is equal to \ref size_, allowing you to index it with the index
+     * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing
+     * and dereferencing the pointer in an internal variable index is undefined behavior.
+     * \see setExternalConnectionNodes()
+     */
+    std::unique_ptr<const ScalarT*[]> yp_ext_;
+    /**
+     * An array of (output) pointers to residuals for external variables.
+     * \note The size of this array is equal to \ref size_, allowing you to index it with the index
+     * of the variable in question (i.e. consisten with \ref extern_indices_). Therefore, accessing
+     * and dereferencing the pointer in an internal variable index is undefined behavior.
+     * \see setExternalConnectionNodes()
+     */
+    std::unique_ptr<ScalarT*[]>       f_ext_;
+
     std::vector<bool> tag_;
     VectorT           abs_tol_;
-    VectorT           f_;
 
     VectorT g_;
 
@@ -463,6 +559,29 @@ namespace GridKit
     IdxT idc_;
 
     bool allocated_{false};
+
+  private:
+    /**
+     * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed.
+     * Instead use \ref y_int_ for an internal variable or \ref y_ext_ for an external variable, respectively.
+     * For components which want an internal buffer (such as a system), make sure that \ref y_int_ points here.
+     * \see allocateVectors()
+     */
+    VectorT y_;
+    /**
+     * The internal buffer for derivatives for the component. For most components, it will be empty and shouldn't be accessed.
+     * Instead use \ref yp_int_ for an internal variable or \ref yp_ext_ for an external variable, respectively.
+     * For components which want an internal buffer (such as a system), make sure that \ref yp_int_ points here.
+     * \see allocateVectors()
+     */
+    VectorT yp_;
+    /**
+     * The internal buffer for state for the component. For most components, it will be empty and shouldn't be accessed.
+     * Instead use \ref f_int_ for an internal variable or \ref f_ext_ for an external variable, respectively.
+     * For components which want an internal buffer (such as a system), make sure that \ref f_int_ points here.
+     * \see allocateVectors()
+     */
+    VectorT f_;
   };
 
 } // namespace GridKit