System external/residual calls

Author: lukelowry

GridKit/Model/PhasorDynamics/Component.hpp

@@ -63,6 +63,26 @@ namespace GridKit
 
       virtual int verify() const = 0;
 
+      /**
+       * @brief Evaluate this component's internal residual.
+       *
+       * The default preserves the existing single-phase residual evaluation
+       * for components that have not yet separated internal and external
+       * residuals.
+       */
+      virtual int evaluateInternalResidual()
+      {
+        return this->evaluateResidual();
+      }
+
+      /**
+       * @brief Evaluate contributions to residuals owned elsewhere.
+       */
+      virtual int evaluateExternalResidual()
+      {
+        return 0;
+      }
+
       IdxT size() override final
       {
         return size_;

GridKit/Model/PhasorDynamics/SystemModel.hpp

@@ -88,6 +88,8 @@ namespace GridKit
 
       int tagDifferentiable() override;
       int setAbsoluteTolerance(RealT rel_tol) override;
+      int evaluateInternalResidual() override;
+      int evaluateExternalResidual() override;
       int evaluateResidual() override;
       int evaluateJacobian() override;
 

GridKit/Model/PhasorDynamics/SystemModelImpl.hpp

@@ -663,20 +663,13 @@ namespace GridKit
     }
 
     /**
-     * @brief Compute system residual vector
+     * @brief Evaluate bus residuals and component internal residuals.
      *
-     * Buses and components read and write their bound system-vector slices
-     * directly.
-     *
-     * @warning Residuals must be computed for buses, before component
-     * residuals are computed. Buses own residuals for currents
-     * Ir and Ii, but the contributions to these residuals come
-     * from components. Buses assign their residual values, while components
-     * add to those values by in-place addition. This is why (for now) bus
-     * residuals need to be computed first.
+     * Buses assign their residuals before components evaluate so that legacy
+     * component contributions continue to accumulate correctly.
      */
     template <typename scalar_type, typename index_type>
-    int SystemModel<scalar_type, index_type>::evaluateResidual()
+    int SystemModel<scalar_type, index_type>::evaluateInternalResidual()
     {
       for (const auto& bus : buses_)
       {
@@ -685,9 +678,35 @@ namespace GridKit
 
       for (const auto& component : components_)
       {
-        component->evaluateResidual();
+        component->evaluateInternalResidual();
+      }
+
+      return 0;
+    }
+
+    /**
+     * @brief Evaluate contributions to residuals owned elsewhere.
+     */
+    template <typename scalar_type, typename index_type>
+    int SystemModel<scalar_type, index_type>::evaluateExternalResidual()
+    {
+      for (const auto& component : components_)
+      {
+        component->evaluateExternalResidual();
       }
 
+      return 0;
+    }
+
+    /**
+     * @brief Compute and assemble the system residual vector.
+     */
+    template <typename scalar_type, typename index_type>
+    int SystemModel<scalar_type, index_type>::evaluateResidual()
+    {
+      evaluateInternalResidual();
+      evaluateExternalResidual();
+
       f_.setDataUpdated();
 
       return 0;

tests/UnitTests/PhasorDynamics/SystemTests.hpp

@@ -98,6 +98,27 @@ namespace GridKit
         }
       };
 
+      class ExternalResidualBranch final : public PhasorDynamics::Branch<ScalarT, IdxT>
+      {
+      public:
+        using BranchT = PhasorDynamics::Branch<ScalarT, IdxT>;
+        using BranchT::BranchT;
+
+        int evaluateExternalResidual() override final
+        {
+          ++external_residual_evaluations_;
+          return 0;
+        }
+
+        std::size_t externalResidualEvaluations() const
+        {
+          return external_residual_evaluations_;
+        }
+
+      private:
+        std::size_t external_residual_evaluations_{0};
+      };
+
     public:
       SystemTests()  = default;
       ~SystemTests() = default;
@@ -206,12 +227,24 @@ namespace GridKit
         PhasorDynamics::Bus<ScalarT, IdxT> bus2(Vr2, Vi2);
         system.addBus(&bus2);
 
-        PhasorDynamics::Branch<ScalarT, IdxT> branch(&bus1, &bus2, R, X, G, B);
+        ExternalResidualBranch branch(&bus1, &bus2, R, X, G, B);
         system.addComponent(&branch);
 
-        system.allocate();
-        system.initialize();
-        system.evaluateResidual();
+        success *= system.allocate() == 0;
+        success *= system.initialize() == 0;
+
+        const auto external_residual_evaluations = branch.externalResidualEvaluations();
+
+        success *= system.evaluateResidual() == 0;
+        success *= branch.externalResidualEvaluations() == external_residual_evaluations + 1;
+
+        success *= isEqual(bus1.Ir(), Ir1);
+        success *= isEqual(bus1.Ii(), Ii1);
+        success *= isEqual(bus2.Ir(), Ir2);
+        success *= isEqual(bus2.Ii(), Ii2);
+
+        success *= system.evaluateResidual() == 0;
+        success *= branch.externalResidualEvaluations() == external_residual_evaluations + 2;
 
         success *= isEqual(bus1.Ir(), Ir1);
         success *= isEqual(bus1.Ii(), Ii1);