Adding LoadZIP model component type

CHANGELOG.md

@@ -48,6 +48,7 @@
 - Refactored Jacobian assembly in `PhasorDyanmcics` module to reuse the CSR pattern.
 - Removed `COO_Matrix` class use in `PowerElectronics` module.
 - Added phasor dynamics application to generalize examples
+- Added LoadZIP model component type.
 
 ## v0.1
 

GridKit/AutomaticDifferentiation/Enzyme/DfDy.hpp

@@ -32,6 +32,79 @@ namespace GridKit
         using RealT   = typename GridKit::ScalarTraits<ScalarT>::RealT;
         using MatrixT = GridKit::LinearAlgebra::COO_Matrix<RealT, IdxT>;
 
+        /**
+         * @param[in] model - Pointer to the model to be differentiated
+         * @param[in] n_res - Number of residual functions
+         * @param[in] n_var - Number of independent variables
+         * @param[in] res_indices - Global residual indices
+         * @param[in] var_indices - Global variable indices
+         * @param[in] y - Internal variables
+         * @param[in] yp - Internal variable derivatives
+         * @param[in] wb - Bus variables
+         * @param[in,out] jac - Jacobian
+         */
+        static void eval(ModelT*     model,
+                         size_t      n_res,
+                         size_t      n_var,
+                         const IdxT* res_indices,
+                         const IdxT* var_indices,
+                         ScalarT*    y,
+                         ScalarT*    yp,
+                         ScalarT*    wb,
+                         IdxT*       rows,
+                         IdxT*       cols,
+                         RealT*      vals,
+                         MatrixT&    jac)
+        {
+          if (n_res > 0 && n_var > 0)
+          {
+            // df/dy
+            std::vector<ScalarT> elementary_v(n_var);
+            IdxT                 nnz = 0;
+            for (size_t var_i = 0; var_i < n_var; ++var_i)
+            {
+              // Sparse storage. @see LowerSparseStorage.hpp
+              ScalarT* output   = __enzyme_todense<ScalarT*>((void*) ident_load<ScalarT, IdxT>,
+                                                           (void*) ident_store<ScalarT, IdxT>,
+                                                           var_i);
+              ScalarT* d_output = __enzyme_todense<ScalarT*>((void*) sparse_load<ScalarT, IdxT>,
+                                                             (void*) sparse_store<ScalarT, IdxT>,
+                                                             var_i,
+                                                             res_indices,
+                                                             var_indices,
+                                                             rows,
+                                                             cols,
+                                                             vals,
+                                                             &nnz);
+
+              // Elementary vector for Jacobian-vector product
+              std::ranges::fill(elementary_v, 0.0);
+              elementary_v[var_i] = 1.0;
+
+              // Core automatic differentiaation intrinsic that will be replaced by a derivative
+              __enzyme_fwddiff<void>((void*) ModelWrapper<ModelT, function, ScalarT>::eval,
+                                     enzyme_const,
+                                     model,
+                                     enzyme_dup,
+                                     y,
+                                     output,
+                                     enzyme_const,
+                                     yp,
+                                     enzyme_const,
+                                     wb,
+                                     enzyme_dupnoneed,
+                                     elementary_v.data(),
+                                     d_output);
+            }
+
+            // Store result
+            jac.setValues(1.0, rows, cols, vals, nnz); //< @todo: Update once sparse storage format changes
+
+            // There is no df/dy' when alpha is not passed as an argument
+            // @todo: Implement a generic way to identify these cases at compile time
+          }
+        }
+
         /**
          * @param[in] model - Pointer to the model to be differentiated
          * @param[in] n_res - Number of residual functions

GridKit/Model/PhasorDynamics/CMakeLists.txt

@@ -36,6 +36,7 @@ add_subdirectory(BusFault)
 add_subdirectory(Exciter)
 add_subdirectory(Governor)
 add_subdirectory(Load)
+add_subdirectory(LoadZIP)
 add_subdirectory(SynchronousMachine)
 
 add_subdirectory(SignalNode)

GridKit/Model/PhasorDynamics/ComponentLibrary.hpp

@@ -6,6 +6,7 @@
 #include <GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp>
 #include <GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp>
 #include <GridKit/Model/PhasorDynamics/Load/Load.hpp>
+#include <GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
 #include <GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
 #include <GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp>

GridKit/Model/PhasorDynamics/LoadZIP/CMakeLists.txt

@@ -0,0 +1,43 @@
+
+set(_install_headers
+  LoadZIP.hpp
+  LoadZIPData.hpp)
+
+if(GRIDKIT_ENABLE_ENZYME)
+  gridkit_add_library(phasor_dynamics_loadzip
+    SOURCES
+      LoadZIPEnzyme.cpp
+    HEADERS
+      ${_install_headers}
+    INCLUDE_DIRECTORIES
+      PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+    LINK_LIBRARIES
+      PRIVATE ClangEnzymeFlags
+      PUBLIC GridKit::phasor_dynamics_core
+    COMPILE_OPTIONS
+      PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno)
+else()
+  gridkit_add_library(phasor_dynamics_loadzip
+    SOURCES
+      LoadZIP.cpp
+    INCLUDE_DIRECTORIES
+      PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+    LINK_LIBRARIES
+      PUBLIC GridKit::phasor_dynamics_core
+    HEADERS
+      ${_install_headers})
+endif()
+
+gridkit_add_library(phasor_dynamics_loadzip_dependency_tracking
+  SOURCES 
+    LoadZIPDependencyTracking.cpp
+  INCLUDE_DIRECTORIES 
+    PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+  LINK_LIBRARIES
+    PUBLIC GridKit::phasor_dynamics_core)
+
+# Link to interface target for all components
+target_link_libraries(phasor_dynamics_components
+  INTERFACE GridKit::phasor_dynamics_loadzip)
+target_link_libraries(phasor_dynamics_components_dependency_tracking
+  INTERFACE GridKit::phasor_dynamics_loadzip_dependency_tracking)

GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.cpp

@@ -0,0 +1,29 @@
+
+#include "LoadZIPImpl.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /**
+     * @brief Jacobian evaluation not implemented
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam IdxT    - matrix index data type
+     * @return int - error code, 0 = success
+     */
+    template <class ScalarT, typename IdxT>
+    int LoadZIP<ScalarT, IdxT>::evaluateJacobian()
+    {
+      Log::misc() << "Evaluate Jacobian for LoadZIP..." << std::endl;
+      Log::misc() << "Jacobian evaluation is not implemented!" << std::endl;
+
+      return 0;
+    }
+
+    // Available template instantiations
+    template class LoadZIP<double, long int>;
+    template class LoadZIP<double, size_t>;
+
+  } // namespace PhasorDynamics
+} // namespace GridKit

GridKit/Model/PhasorDynamics/LoadZIP/LoadZIP.hpp

@@ -0,0 +1,140 @@
+#pragma once
+
+#include <GridKit/Model/PhasorDynamics/Component.hpp>
+#include <GridKit/Model/PhasorDynamics/ComponentSignals.hpp>
+#include <GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPData.hpp>
+#include <GridKit/Model/VariableMonitor.hpp>
+
+// Forward declarations.
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    template <class ScalarT, typename IdxT>
+    class BusBase;
+
+    template <typename RealT, typename IdxT>
+    struct LoadZIPData;
+  } // namespace PhasorDynamics
+} // namespace GridKit
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /*!
+     * @brief Implementation of a ZIP load.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    class LoadZIP : public Component<ScalarT, IdxT>
+    {
+      using Component<ScalarT, IdxT>::gridkit_component_id_;
+      using Component<ScalarT, IdxT>::size_;
+      using Component<ScalarT, IdxT>::nnz_;
+      using Component<ScalarT, IdxT>::time_;
+      using Component<ScalarT, IdxT>::alpha_;
+      using Component<ScalarT, IdxT>::y_;
+      using Component<ScalarT, IdxT>::yp_;
+      using Component<ScalarT, IdxT>::tag_;
+      using Component<ScalarT, IdxT>::wb_;
+      using Component<ScalarT, IdxT>::h_;
+      using Component<ScalarT, IdxT>::f_;
+      using Component<ScalarT, IdxT>::J_;
+      using Component<ScalarT, IdxT>::J_rows_buffer_;
+      using Component<ScalarT, IdxT>::J_cols_buffer_;
+      using Component<ScalarT, IdxT>::J_vals_buffer_;
+      using Component<ScalarT, IdxT>::variable_indices_;
+      using Component<ScalarT, IdxT>::residual_indices_;
+
+    public:
+      using RealT           = typename Component<ScalarT, IdxT>::RealT;
+      using bus_type        = BusBase<ScalarT, IdxT>;
+      using model_data_type = LoadZIPData<RealT, IdxT>;
+      using MonitorT        = Model::VariableMonitor<LoadZIP, LoadZIPData>;
+
+      LoadZIP(bus_type* bus);
+      LoadZIP(bus_type* bus, RealT P0, RealT Q0, RealT V0, RealT alphaI, RealT alphaP);
+      LoadZIP(bus_type* bus, const model_data_type& data);
+      ~LoadZIP();
+
+      int setGridKitComponentID(IdxT) override final;
+      int allocate() override final;
+      int initialize() override final;
+      int tagDifferentiable() override final;
+      int evaluateResidual() override final;
+      int evaluateJacobian() override final;
+
+      int verify() const override final
+      {
+        return 0;
+      }
+
+    public:
+      void setP0(RealT P0)
+      {
+        P0_ = P0;
+      }
+
+      void setQ0(RealT Q0)
+      {
+        Q0_ = Q0;
+      }
+
+      void setV0(RealT V0)
+      {
+        V0_ = V0;
+      }
+
+      void setalphaI(RealT alphaI)
+      {
+        alphaI_ = alphaI;
+      }
+
+      void setalphaP(RealT alphaP)
+      {
+        alphaP_ = alphaP;
+      }
+
+    private:
+      void setDerivedParams();
+
+      ScalarT& Vr()
+      {
+        return bus_->Vr();
+      }
+
+      ScalarT& Vi()
+      {
+        return bus_->Vi();
+      }
+
+      ScalarT& Ir()
+      {
+        return bus_->Ir();
+      }
+
+      ScalarT& Ii()
+      {
+        return bus_->Ii();
+      }
+
+      const Model::VariableMonitorBase* getMonitor() const override;
+
+    public:
+      __attribute__((always_inline)) inline int evaluateBusResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
+      __attribute__((always_inline)) inline int evaluateInternalResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*);
+
+    private:
+      bus_type* bus_{nullptr};
+      RealT     P0_{0};
+      RealT     Q0_{0};
+      RealT     V0_{1.0};
+      RealT     alphaI_{0};
+      RealT     alphaP_{0};
+
+      std::unique_ptr<MonitorT> monitor_;
+    };
+
+  } // namespace PhasorDynamics
+} // namespace GridKit

GridKit/Model/PhasorDynamics/LoadZIP/LoadZIPData.hpp

@@ -0,0 +1,58 @@
+
+#pragma once
+
+#include <GridKit/Model/PhasorDynamics/ComponentData.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /// Initial parameters for a loadZIP
+    enum class LoadZIPParameters
+    {
+      P0,     ///< Load nominal real power
+      Q0,     ///< Load nominal reactive power
+      V0,     ///< Load nominal reactive power
+      alphaI, ///< Fraction of load to be represented as constant current
+      alphaP, ///< Fraction of load to be represented as constant power
+    };
+
+    /// Ports for a loadZIP
+    enum class LoadZIPPorts
+    {
+      bus, ///< Unique ID of the bus to which the loadZIP is connected
+    };
+
+    /// Variables able to be monitored for a loadZIP
+    enum class LoadZIPMonitorableVariables
+    {
+      Ir,
+      Ii,
+      Im,
+      P,
+      Q
+    };
+
+    /**
+     * @brief Contains modeling data for a load
+     *
+     * @tparam RealT Real parameter data type
+     * @tparam IdxT  Integer parameter data type
+     *
+     * Integer parameters are of the same type as matrix and vector indices.
+     */
+    template <typename RealT, typename IdxT>
+    struct LoadZIPData : public ComponentData<RealT,
+                                              IdxT,
+                                              LoadZIPParameters,
+                                              LoadZIPPorts,
+                                              LoadZIPMonitorableVariables>
+    {
+      LoadZIPData() = default;
+
+      using Parameters           = LoadZIPParameters;
+      using Ports                = LoadZIPPorts;
+      using MonitorableVariables = LoadZIPMonitorableVariables;
+    };
+  } // namespace PhasorDynamics
+} // namespace GridKit