Initial Implementation of Subsystem Model for Partition Simulation

CHANGELOG.md

@@ -80,6 +80,7 @@
 - Added `HYGOV` governor model implementation for PhasorDynamics.
 - Added `REPCA` controller model implementation for PhasorDynamics.
 - Added `REECB` electrical-control model implementation for PhasorDynamics.
+- Added subsystem partitioning support for `PowerElectronics` models, including partition interfaces, and independent residual and Jacobian evaluation.
 
 ## v0.1
 

GridKit/Model/PowerElectronics/CMakeLists.txt

@@ -21,6 +21,7 @@ add_subdirectory(TransmissionLine)
 add_subdirectory(MicrogridLoad)
 add_subdirectory(MicrogridLine)
 add_subdirectory(MicrogridBusDQ)
+add_subdirectory(PartitionInterface)
 
 install(
   FILES CircuitComponent.hpp
@@ -29,4 +30,5 @@ install(
         SystemModelPowerElectronics.hpp
         NodeBase.hpp
         ExternalConnection.hpp
+        SubsystemModel.hpp
   DESTINATION include/GridKit/Model/PowerElectronics)

GridKit/Model/PowerElectronics/Capacitor/Capacitor.cpp

@@ -52,6 +52,25 @@ namespace GridKit
     return 0;
   }
 
+  /**
+   * @brief Compute the absolute tolerance for each variable in the model
+   *
+   * @param rel_tol The relative tolerance which can be used to pick the
+   *        absolute tolerance.
+   * @tparam ScalarT Scalar data type
+   * @tparam IdxT Index data type
+   * @return int 0 if successful, non-zero otherwise.
+   *
+   * This represents a "noise" level close to zero for which pure relative
+   * error cannot be used.
+   */
+  template <class ScalarT, typename IdxT>
+  int Capacitor<ScalarT, IdxT>::setAbsoluteTolerance(RealT rel_tol)
+  {
+    abs_tol_.setToConst(static_cast<ScalarT>(rel_tol));
+    return 0;
+  }
+
   /**
    * @brief Evaluate the resisdual of the Capcitor
    *
@@ -117,6 +136,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool Capacitor<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* Capacitor<ScalarT, IdxT>::clone() const
+  {
+    return new Capacitor<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class Capacitor<double, long int>;
   template class Capacitor<double, size_t>;

GridKit/Model/PowerElectronics/Capacitor/Capacitor.hpp

@@ -29,6 +29,7 @@ namespace GridKit
     using CircuitComponent<ScalarT, IdxT>::y_int_;
     using CircuitComponent<ScalarT, IdxT>::yp_ext_;
     using CircuitComponent<ScalarT, IdxT>::yp_int_;
+    using CircuitComponent<ScalarT, IdxT>::abs_tol_;
     using CircuitComponent<ScalarT, IdxT>::tag_;
     using CircuitComponent<ScalarT, IdxT>::f_ext_;
     using CircuitComponent<ScalarT, IdxT>::f_int_;
@@ -50,15 +51,19 @@ namespace GridKit
 
     int initialize();
     int tagDifferentiable();
+    int setAbsoluteTolerance(RealT);
     int evaluateInternalResidual() final;
     int evaluateExternalResidual() final;
     int evaluateJacobian();
     int evaluateIntegrand();
 
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT C_;

GridKit/Model/PowerElectronics/CircuitComponent.hpp

@@ -27,6 +27,138 @@ namespace GridKit
 
     CircuitComponent() = default;
 
+    CircuitComponent(const CircuitComponent& other)
+      : n_extern_(other.n_extern_),
+        n_intern_(other.n_intern_),
+        extern_indices_(other.extern_indices_),
+        size_(other.size_),
+        nnz_(other.nnz_),
+        size_quad_(other.size_quad_),
+        size_opt_(other.size_opt_),
+        current_jac_size_(other.current_jac_size_),
+
+        // These pointers refer to storage supplied by a parent system.
+        // The copied component must be connected to its own storage later.
+        y_int_(nullptr),
+        yp_int_(nullptr),
+        f_int_(nullptr),
+
+        tag_(other.tag_),
+        time_(other.time_),
+        alpha_(other.alpha_),
+        max_steps_(other.max_steps_),
+        idc_(other.idc_),
+        allocated_(other.allocated_)
+    {
+      /*
+       * VectorT disables its normal copy constructor and copy-assignment
+       * operator. Use its provided copyFromExternal() operation to perform
+       * an independent copy of the vector data.
+       */
+      auto copyVector = [](VectorT& destination, const VectorT& source)
+      {
+        const IdxT source_size = source.getSize();
+
+        if (source_size == 0)
+        {
+          return;
+        }
+
+        destination.resize(source_size);
+        destination.copyFromExternal(source);
+      };
+
+      /*
+       * Deep-copy the local-to-global connection mapping.
+       */
+      if (other.connection_nodes_)
+      {
+        connection_nodes_ = std::make_unique<IdxT[]>(static_cast<size_t>(size_));
+
+        for (size_t i = 0; i < static_cast<size_t>(size_); ++i)
+        {
+          connection_nodes_[i] = other.connection_nodes_[i];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian row indices.
+       */
+      if (other.jacobian_coo_rows_)
+      {
+        jacobian_coo_rows_ = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
+
+        for (size_t i = 0; i < static_cast<size_t>(nnz_); ++i)
+        {
+          jacobian_coo_rows_[i] = other.jacobian_coo_rows_[i];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian column indices.
+       */
+      if (other.jacobian_coo_cols_)
+      {
+        jacobian_coo_cols_ = std::make_unique<IdxT[]>(static_cast<size_t>(nnz_));
+
+        for (size_t i = 0; i < static_cast<size_t>(nnz_); ++i)
+        {
+          jacobian_coo_cols_[i] = other.jacobian_coo_cols_[i];
+        }
+      }
+
+      /*
+       * Deep-copy the COO Jacobian values.
+       */
+      if (other.jacobian_coo_values_)
+      {
+        jacobian_coo_values_ = std::make_unique<RealT[]>(static_cast<size_t>(nnz_));
+
+        for (size_t i = 0; i < static_cast<size_t>(nnz_); ++i)
+        {
+          jacobian_coo_values_[i] = other.jacobian_coo_values_[i];
+        }
+      }
+
+      if (size_ > 0)
+      {
+        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_));
+
+        for (size_t i = 0; i < static_cast<size_t>(size_); ++i)
+        {
+          y_ext_[i]  = nullptr;
+          yp_ext_[i] = nullptr;
+          f_ext_[i]  = nullptr;
+        }
+      }
+
+      // State, state derivative, residual, and absolute tolerance.
+      copyVector(y_, other.y_);
+      copyVector(yp_, other.yp_);
+      copyVector(f_, other.f_);
+      copyVector(abs_tol_, other.abs_tol_);
+      copyVector(g_, other.g_);
+      copyVector(yB_, other.yB_);
+      copyVector(ypB_, other.ypB_);
+      copyVector(fB_, other.fB_);
+      copyVector(gB_, other.gB_);
+      copyVector(param_, other.param_);
+      copyVector(param_up_, other.param_up_);
+      copyVector(param_lo_, other.param_lo_);
+    }
+
+    virtual CircuitComponent<ScalarT, IdxT>* clone() const
+    {
+      return nullptr;
+    }
+
+    virtual bool isCloneable() const
+    {
+      return false;
+    }
+
     /**
      * @note Cannot be marked final, since it is overriden to recurse in the system model.
      */
@@ -51,7 +183,7 @@ namespace GridKit
       return this->n_intern_;
     }
 
-    std::set<size_t> getExternIndices()
+    std::set<IdxT> getExternIndices()
     {
       return this->extern_indices_;
     }
@@ -69,7 +201,7 @@ namespace GridKit
     int setInternalConnectionNodes(size_t local_index, IdxT global_index)
     {
       assert(!extern_indices_.contains(static_cast<IdxT>(local_index)));
-      connection_nodes_[local_index] = global_index;
+      setConnectionNodes(local_index, global_index);
       return 0;
     }
 
@@ -88,10 +220,27 @@ namespace GridKit
     int setExternalConnectionNodes(size_t local_index, ExternalConnection<ScalarT, IdxT> connection)
     {
       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_;
+      y_ext_[local_index]  = connection.y_;
+      yp_ext_[local_index] = connection.yp_;
+      f_ext_[local_index]  = connection.f_;
+      setConnectionNodes(local_index, connection.idx_);
+      return 0;
+    }
+
+    /**
+     * @brief Update the connection index for a variable.
+     *
+     * Sets only the connection index without modifying the variable's
+     * internal/external classification or its associated data pointers.
+     *
+     * @param local_index Index of the local variable.
+     * @param connection_index New connection index for the variable.
+     *
+     * @return int 0 if successful.
+     */
+    int setConnectionNodes(size_t local_index, IdxT connection_index)
+    {
+      connection_nodes_[local_index] = connection_index;
       return 0;
     }
 
@@ -132,9 +281,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_));
 
-      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_));
+      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<IdxT[]>(static_cast<size_t>(size_));
 
       if (!allocated_)
@@ -439,6 +589,16 @@ namespace GridKit
       return idc_;
     }
 
+    /**
+     * @brief Check whether the component has already been allocated.
+     *
+     * @return true if allocate() has previously completed, false otherwise.
+     */
+    bool isAllocated() const
+    {
+      return allocated_;
+    }
+
   protected:
     /**
      * @brief Allocate state and residual storage owned by this component.

GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp

@@ -421,6 +421,18 @@ namespace GridKit
     return 0;
   }
 
+  template <class ScalarT, typename IdxT>
+  bool DistributedGenerator<ScalarT, IdxT>::isCloneable() const
+  {
+    return true;
+  }
+
+  template <class ScalarT, typename IdxT>
+  CircuitComponent<ScalarT, IdxT>* DistributedGenerator<ScalarT, IdxT>::clone() const
+  {
+    return new DistributedGenerator<ScalarT, IdxT>(*this);
+  }
+
   // Available template instantiations
   template class DistributedGenerator<double, long int>;
   template class DistributedGenerator<double, size_t>;

GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp

@@ -75,18 +75,21 @@ namespace GridKit
                          NodeT*                                      node_bus);
     virtual ~DistributedGenerator();
 
-    int initialize();
-    int allocate() final;
-    int tagDifferentiable();
-    int setAbsoluteTolerance(RealT);
-    int evaluateInternalResidual() final;
-    int evaluateExternalResidual() final;
-    int evaluateJacobian();
-    int evaluateIntegrand();
-    int initializeAdjoint();
-    int evaluateAdjointResidual();
+    int  initialize();
+    int  allocate() final;
+    int  tagDifferentiable();
+    int  setAbsoluteTolerance(RealT);
+    int  evaluateInternalResidual() final;
+    int  evaluateExternalResidual() final;
+    int  evaluateJacobian();
+    int  evaluateIntegrand();
+    int  initializeAdjoint();
+    int  evaluateAdjointResidual();
     // int evaluateAdjointJacobian();
-    int evaluateAdjointIntegrand();
+    int  evaluateAdjointIntegrand();
+    bool isCloneable() const;
+
+    CircuitComponent<ScalarT, IdxT>* clone() const;
 
   private:
     RealT wb_;