Jacobian test for Exciter::SexsPti

GridKit/Model/PhasorDynamics/README.md

@@ -64,15 +64,16 @@ We recommend developers follow these steps when adding new component models:
 5. Once model is tested, add it to the system composer. This requires following steps:
     1. Add header file `MyModel.hpp` to `ComponentLibrary.hpp`, so that
        `MyModel` declaration is visible to the `SystemModel` class.
-    3. Modify `SystemModelJsonParser.hpp` so that `MyModel` is recognized by the
+    2. Modify `SystemModelJsonParser.hpp` so that `MyModel` is recognized by the
        parser.
-    4. Modify `SystemModelData.hpp` so that `MyModelData` is visible to the system
+    3. Modify `SystemModelData.hpp` so that `MyModelData` is visible to the system
        model.
-    5. Modify `SystemModel.hpp` so that `MyModel` components can be connected by the
+    4. Modify `SystemModel.hpp` so that `MyModel` components can be connected by the
        system composer.
-7. Recommended: Create an example in `examples/PhasorDynamics` using the new component.
-
-
+    5. Recommended: If applicable, add a smoke test to
+       `tests/UnitTesting/PhasorDynamics/SystemSingleComponentTests.hpp` for the
+       basic construction of `MyModel` through `SystemModel`.
+6. Recommended: Create an example in `examples/PhasorDynamics` using the new component.
 
 ## Input file parser
 

GridKit/Model/PowerElectronics/NodeBase.hpp

@@ -103,6 +103,7 @@ namespace GridKit
       int setBusID(IdxT bus_id)
       {
         bus_id_ = bus_id;
+        return 0;
       }
 
       IdxT busID() const

buildsystem/spack_repo/gridkit/packages/gridkit/package.py

@@ -20,7 +20,7 @@ class Gridkit(CMakePackage):
     variant("ipopt", default=False, description="Enable/Disable Ipopt")
     variant("klu", default=False, description="Enable/Disable KLU")
     variant("sundials", default=False, description="Enable/Disable SUNDIALS")
-    variant("ubsan", default=False, description="Enable/Disable undefined behavir sanitizer")
+    variant("ubsan", default=False, description="Enable/Disable undefined behavior sanitizer")
 
     conflicts("+klu", when="~sundials")
 

tests/UnitTests/PhasorDynamics/ExciterSexsPtiTests.hpp

@@ -13,6 +13,7 @@
 #include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
 #include <GridKit/Testing/Testing.hpp>
+#include <GridKit/Utilities/MapFromCOO.hpp>
 
 namespace GridKit
 {
@@ -289,6 +290,160 @@ namespace GridKit
         return success.report(__func__);
       }
 
+#ifdef GRIDKIT_ENABLE_ENZYME
+      /**
+       * @brief Checks Jacobian evaluation.
+       */
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        auto tol = 10 * std::numeric_limits<RealT>::epsilon();
+
+        // Jacobian via DependencyTracking
+        std::vector<DependencyTracking::Variable::DependencyMap> dependency_tracking_jacobian = DependencyTrackingJacobian();
+
+        // Jacobian via Enzyme
+        std::vector<DependencyTracking::Variable::DependencyMap> enzyme_jacobian = EnzymeJacobian();
+
+        /// Compare DependencyTracking dependencies to Enzyme's
+        for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i)
+        {
+          success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol));
+        }
+
+        return success.report(__func__);
+      }
+
+    private:
+      std::vector<DependencyTracking::Variable::DependencyMap> DependencyTrackingJacobian()
+      {
+        auto data = makeTestData();
+
+        DependencyTracking::Variable                                         Vr1{3.0};
+        DependencyTracking::Variable                                         Vi1{4.0};
+        PhasorDynamics::Bus<DependencyTracking::Variable, IdxT>              bus(Vr1, Vi1);
+        PhasorDynamics::Exciter::SexsPti<DependencyTracking::Variable, IdxT> exciter(&bus, data);
+
+        bus.allocate();
+        exciter.allocate();
+
+        // Get d/dy
+        bus.initialize();
+        exciter.initialize();
+
+        for (size_t i = 0; i < exciter.size(); ++i)
+        {
+          exciter.y()[i].setVariableNumber(i); ///< Exciter independent variables
+        }
+        for (size_t i = 0; i < bus.size(); ++i)
+        {
+          bus.y()[i].setVariableNumber(i + exciter.size()); // Bus independent variables
+        }
+
+        bus.evaluateResidual();
+        exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
+                                    ///< the dependencies
+        std::vector<DependencyTracking::Variable> residual_y = exciter.getResidual();
+
+        // Get d/dy'
+        bus.initialize();
+        exciter.initialize();
+
+        for (size_t i = 0; i < exciter.size(); ++i)
+        {
+          exciter.yp()[i].setVariableNumber(i); ///< Exciter independent variables
+        }
+
+        bus.evaluateResidual();
+        exciter.evaluateResidual(); ///< Computes the residual and the Jacobian values by tracking
+                                    ///< the dependencies
+        std::vector<DependencyTracking::Variable> residual_yp = exciter.getResidual();
+
+        // Print the dependencies
+        for (size_t i = 0; i < residual_y.size(); ++i)
+        {
+          std::cout << i << "th residual, y: ";
+          (residual_y[i]).print(std::cout);
+          std::cout << "\n";
+          std::cout << i << "th residual, yp: ";
+          (residual_yp[i]).print(std::cout);
+          std::cout << "\n";
+        }
+
+        // Extract the dependencies and add d/dy' to d/dy
+        std::vector<DependencyTracking::Variable::DependencyMap> dependencies(residual_y.size());
+        for (IdxT i = 0; i < residual_y.size(); ++i)
+        {
+          DependencyTracking::Variable::DependencyMap dependency_y  = (residual_y[i]).getDependencies();
+          DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies();
+
+          for (const auto& pair_y : dependency_y)
+          {
+            auto index_y = pair_y.first;
+            auto value_y = pair_y.second;
+            auto it_yp   = dependency_yp.find(index_y);
+            if (it_yp != dependency_yp.end())
+            {
+              auto value_yp = it_yp->second;
+              dependencies[i].insert(std::make_pair(index_y, value_y + value_yp));
+            }
+            else
+            {
+              dependencies[i].insert(std::make_pair(index_y, value_y));
+            }
+          }
+
+          // Insert yp dependencies that did not exist in the y dependencies
+          for (const auto& pair_yp : dependency_yp)
+          {
+            auto index_yp = pair_yp.first;
+            auto value_yp = pair_yp.second;
+            auto it_y     = dependency_y.find(index_yp);
+            if (it_y == dependency_y.end())
+            {
+              dependencies[i].insert(std::make_pair(index_yp, value_yp));
+            }
+          }
+        }
+
+        return dependencies;
+      }
+
+      std::vector<DependencyTracking::Variable::DependencyMap> EnzymeJacobian()
+      {
+        auto data = makeTestData();
+
+        PhasorDynamics::Bus<ScalarT, IdxT>              bus(3.0, 4.0);
+        PhasorDynamics::Exciter::SexsPti<ScalarT, IdxT> exciter(&bus, data);
+
+        bus.allocate();
+        exciter.allocate();
+
+        bus.initialize();
+        exciter.initialize();
+
+        exciter.updateTime(0.0, 1.0); // Set alpha to 1.0 to verify d/dy' term
+
+        for (size_t i = 0; i < bus.size(); ++i)
+        {
+          bus.setVariableIndex(i, i + exciter.size()); // Reset bus variable indices
+          bus.setResidualIndex(i, i + exciter.size()); // Reset bus residual indices
+        }
+
+        bus.evaluateResidual();
+        exciter.evaluateResidual();
+
+        bus.evaluateJacobian();
+        exciter.evaluateJacobian();
+        GridKit::LinearAlgebra::COO_Matrix<ScalarT, IdxT>& model_jacobian = exciter.getJacobian();
+        model_jacobian.deduplicate();
+        model_jacobian.printMatrix("Model Jacobian");
+
+        return GridKit::Testing::MapFromCOO(model_jacobian);
+      }
+#endif
+
     private:
       auto makeTestData() -> PhasorDynamics::Exciter::SexsPtiData<RealT, IdxT>
       {

tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp

@@ -142,6 +142,30 @@ namespace GridKit
         return success.report(__func__);
       }
 
+      TestOutcome loadZIP()
+      {
+        TestStatus success = true;
+
+        PhasorDynamics::SystemModel<ScalarT, IdxT>* system = new PhasorDynamics::SystemModel<ScalarT, IdxT>();
+
+        PhasorDynamics::BusInfinite<ScalarT, IdxT> bus;
+        system->addBus(&bus);
+
+        PhasorDynamics::LoadZIP<ScalarT, IdxT> load(&bus);
+        system->addComponent(&load);
+
+        success *= system->allocate() == 0;
+        success *= system->initialize() == 0;
+        success *= system->evaluateResidual() == 0;
+        success *= system->evaluateJacobian() == 0;
+        success *= system->size() == load.size();
+
+        delete system;
+        system = nullptr;
+
+        return success.report(__func__);
+      }
+
       TestOutcome genrou()
       {
         TestStatus success = true;

tests/UnitTests/PhasorDynamics/runExciterSexsPtiTests.cpp

@@ -13,6 +13,9 @@ int main()
   result += test.parameterValidation();
   result += test.systemAssembly();
   result += test.systemSignalFreshness();
+#ifdef GRIDKIT_ENABLE_ENZYME
+  result += test.jacobian();
+#endif
 
   return result.summary();
 }

tests/UnitTests/PhasorDynamics/runGenrouTests.cpp

@@ -10,9 +10,9 @@ int main()
   result += test.accessors();
   result += test.hard_coded_residual();
   result += test.residual();
-
 #ifdef GRIDKIT_ENABLE_ENZYME
   result += test.jacobian();
 #endif
+
   return result.summary();
 }

tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp

@@ -13,9 +13,14 @@ int main()
   result += test.busFault();
   result += test.ieeet1();
   result += test.load();
+  result += test.loadZIP();
   result += test.genrou();
   result += test.genClassical();
   result += test.tgov1();
 
+  // @todo The following components are not tested here because they require non-trivial constructors
+  // PhasorDynamics::Exciter::SexsPti
+  // PhasorDynamics::Stabilizer::Ieeest
+
   return result.summary();
 }