Add automatic CSR Jacobian construction with DependencyTracking in PhasorDyanmics

GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp

@@ -13,7 +13,6 @@
 #include <map>
 #include <string>
 #include <utility>
-#include <vector>
 
 #include <GridKit/Constants.hpp>
 #include <GridKit/ScalarTraits.hpp>
@@ -209,7 +208,8 @@ namespace GridKit
       */
       double der(size_t i) const
       {
-        return dependencies_[i];
+        auto it = dependencies_.find(i);
+        return it != dependencies_.end() ? it->second : 0.0;
       }
 
       /**
@@ -255,6 +255,10 @@ namespace GridKit
 
       /**
           @brief Turns variable into parameter, or vice versa.
+
+          @todo is_fixed_ is currently not contributing to the semantics of
+                the derivatives. Leaving as-is for now, as it is not used
+                for anything other than printed diagnostics.
        */
       void setFixed(bool b = false)
       {
@@ -265,10 +269,6 @@ namespace GridKit
       using DependencyMap = std::map<size_t, double>;
       inline const DependencyMap& getDependencies() const;
 
-      // set as the independent state variable and assign ID to it
-      inline void registerVariable(std::vector<Variable*>& x,
-                                   const size_t&           offset);
-
       // adds all dependencies of v to *this
       inline void addDependencies(const Variable& v);
 
@@ -305,8 +305,8 @@ namespace GridKit
       size_t variable_number_; ///< Independent variable ID
       bool   is_fixed_;        ///< Constant parameter flag.
 
-      mutable DependencyMap dependencies_;
-      static const size_t   INVALID_VAR_NUMBER = INVALID_INDEX<size_t>;
+      DependencyMap       dependencies_;
+      static const size_t INVALID_VAR_NUMBER = INVALID_INDEX<size_t>;
     };
 
     //------------------------------------

GridKit/AutomaticDifferentiation/DependencyTracking/VariableImplementation.hpp

@@ -16,19 +16,6 @@ namespace GridKit
       return dependencies_;
     }
 
-    /**
-        @brief Registers a variable as an unknown of the system and
-        adds a pointer to the global @a x vector.
-    */
-    void Variable::registerVariable(std::vector<Variable*>& x,
-                                    const size_t&           offset)
-    {
-      setVariableNumber(offset); // define global variable number
-      setFixed(false);           // not a constant
-
-      x[offset] = this;
-    }
-
     /**
         @brief Adds all dependencies of v to *this.
     */

GridKit/Model/PhasorDynamics/Branch/BranchDependencyTracking.cpp

@@ -11,15 +11,19 @@ namespace GridKit
   namespace PhasorDynamics
   {
     /**
-     * @brief Jacobian evaluation not implemented
+     * @brief Evaluate DependencyTracking::Variable Jacobian.
+     *
+     * @note Currently only used for testing.
+     *
+     * No-op, because branch currently does not own residual equations.
      *
      * @return int - error code, 0 = success
      */
     template <typename scalar_type, typename index_type>
     int Branch<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Branch..." << std::endl;
-      Log::misc() << "Jacobian evaluation is not implemented!" << std::endl;
+      Log::misc() << "Evaluate DependencyTracking Jacobian for Branch...\n";
+      Log::misc() << "Jacobian evaluation is experimental!\n";
 
       return 0;
     }

GridKit/Model/PhasorDynamics/Bus/Bus.hpp

@@ -124,6 +124,35 @@ namespace GridKit
         return 0;
       }
 
+      /**
+       * @brief Initialize DependencyTracking variable numbers.
+       *
+       * @note Assigns even indices to y and odd indices to yp.
+       *       Should be called in intialize(), after variables have been set (and updated as needed).
+       */
+      int initializeDependencyTrackingVariableNumbers()
+        requires std::is_same_v<ScalarT, DependencyTracking::Variable>
+      {
+        auto* y  = y_.getData();
+        auto* yp = yp_.getData();
+
+        for (IdxT j = 0; j < size_; ++j)
+        {
+          const IdxT var_idx = this->getVariableIndex(j);
+          if (var_idx != INVALID_INDEX<IdxT>)
+          {
+            // Even indices for y and odd indices for yp
+            y[j].setVariableNumber(static_cast<size_t>(2 * var_idx));
+            yp[j].setVariableNumber(static_cast<size_t>(2 * var_idx + 1));
+          }
+        }
+
+        y_.setDataUpdated();
+        yp_.setDataUpdated();
+
+        return 0;
+      }
+
       IdxT*  J_rows_buffer_{nullptr};
       IdxT*  J_cols_buffer_{nullptr};
       RealT* J_vals_buffer_{nullptr};

GridKit/Model/PhasorDynamics/Bus/BusDependencyTracking.cpp

@@ -1,3 +1,8 @@
+/**
+ * @file BusDependencyTracking.cpp
+ * @author Slaven Peles (peless@ornl.gov)
+ *
+ */
 
 #include "BusImpl.hpp"
 
@@ -6,15 +11,21 @@ namespace GridKit
   namespace PhasorDynamics
   {
     /**
-     * @brief Jacobian evaluation not implemented
+     * @brief Evaluate DependencyTracking::Variable Jacobian.
+     *
+     * @note Currently only used for testing.
+     *
+     * DependencyTracking::Variable stores the Jacobian as dependency maps,
+     * updated during calls to evaluateResidual().
+     * Not yet implemented for bus residuals.
      *
      * @return int - error code
      */
     template <typename scalar_type, typename index_type>
     int Bus<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for Bus..." << std::endl;
-      Log::misc() << "Jacobian evaluation is not implemented!" << std::endl;
+      Log::misc() << "Evaluate DependencyTracking Jacobian for Bus...\n";
+      Log::misc() << "Jacobian evaluation is not implemented!\n";
 
       return 0;
     }

GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp

@@ -185,6 +185,12 @@ namespace GridKit
       y_.setDataUpdated();
       yp_.setDataUpdated();
 
+      // For DependencyTracking::Variable, set variable numbers
+      if constexpr (std::is_same_v<ScalarT, DependencyTracking::Variable>)
+      {
+        this->initializeDependencyTrackingVariableNumbers();
+      }
+
       return 0;
     }
 

GridKit/Model/PhasorDynamics/Bus/BusInfiniteDependencyTracking.cpp

@@ -1,3 +1,8 @@
+/**
+ * @file BusInfiniteDependencyTracking.cpp
+ * @author Slaven Peles (peless@ornl.gov)
+ *
+ */
 
 #include "BusInfiniteImpl.hpp"
 

GridKit/Model/PhasorDynamics/BusBase.hpp

@@ -252,7 +252,6 @@ namespace GridKit
        */
       void allocateVectors(IdxT n)
       {
-
         y_.resize(n);
         yp_.resize(n);
         f_.resize(n);

GridKit/Model/PhasorDynamics/BusFault/BusFaultDependencyTracking.cpp

@@ -1,23 +1,37 @@
+/**
+ * @file BusFaultDependencyTracking.cpp
+ * @author Slaven Peles (peless@ornl.gov)
+ *
+ */
+
 #include "BusFaultImpl.hpp"
 
 namespace GridKit
 {
   namespace PhasorDynamics
   {
     /**
-     * @brief Jacobian evaluation not implemented yet
+     * @brief Evaluate DependencyTracking::Variable Jacobian.
+     *
+     * @note Currently only used for testing.
+     *
+     * DependencyTracking::Variable stores the Jacobian as dependency maps,
+     * updated during calls to evaluateResidual().
      *
      * @return int - error code, 0 = success
      */
     template <class scalar_type, typename index_type>
     int BusFault<scalar_type, index_type>::evaluateJacobian()
     {
-      Log::misc() << "Evaluate Jacobian for BusFault..." << std::endl;
-      Log::misc() << "Jacobian evaluation not implemented!" << std::endl;
+      Log::misc() << "Evaluate DependencyTracking Jacobian for BusFault...\n";
+      Log::misc() << "Jacobian evaluation is experimental!\n";
+
+      this->constructCsr();
+
       return 0;
     }
 
-    // Additional template instantiations
+    // Available template instantiations
     template class BusFault<DependencyTracking::Variable, long int>;
     template class BusFault<DependencyTracking::Variable, size_t>;
   } // namespace PhasorDynamics

GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp

@@ -170,6 +170,12 @@ namespace GridKit
       y_.setDataUpdated();
       yp_.setDataUpdated();
 
+      // For DependencyTracking::Variable, set variable numbers
+      if constexpr (std::is_same_v<scalar_type, DependencyTracking::Variable>)
+      {
+        this->initializeDependencyTrackingVariableNumbers();
+      }
+
       return 0;
     }
 

GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapterDependencyTracking.cpp

@@ -2,8 +2,6 @@
  * @file BusToSignalAdapterDependencyTracking.cpp
  * @author Philip Fackler (facklerpw@ornl.gov)
  *
- * @brief Deinition of BusToSignalAdapter connector interface.
- *
  */
 
 #include "BusToSignalAdapterImpl.hpp"