diff --git a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp
index d53a653e6..6245ad76d 100644
--- a/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp
+++ b/GridKit/AutomaticDifferentiation/Enzyme/ModelWrappers.hpp
@@ -63,7 +63,7 @@ namespace GridKit
                          const ScalarT* wb,
                          ScalarT*       f)
         {
-          model->evaluateInternalResidual(y, yp, wb, f);
+          model->evaluateInternalResidualKernel(y, yp, wb, f);
         }
       };
 
@@ -91,7 +91,7 @@ namespace GridKit
                          const ScalarT* ws,
                          ScalarT*       f)
         {
-          model->evaluateInternalResidual(y, yp, wb, ws, f);
+          model->evaluateInternalResidualKernel(y, yp, wb, ws, f);
         }
       };
 
diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp
index 2a4bf6365..a2253b0f8 100644
--- a/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp
+++ b/GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp
@@ -118,7 +118,7 @@ namespace GridKit
     public:
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
     private:
diff --git a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp
index 2a2953360..3e84577d2 100644
--- a/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/BusFault/BusFaultImpl.hpp
@@ -228,7 +228,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) int BusFault<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) int BusFault<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT*                  y,
         [[maybe_unused]] const ScalarT* yp,
         const ScalarT*                  wb,
@@ -259,7 +259,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), f);
         evaluateBusResidual(y, yp, wb_.data(), h_.data());
         Ir() += h_[0];
         Ii() += h_[1];
@@ -275,7 +275,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), f);
       }
 
       f_.setDataUpdated();
diff --git a/GridKit/Model/PhasorDynamics/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp
index a930c8ffb..b27876bcc 100644
--- a/GridKit/Model/PhasorDynamics/Component.hpp
+++ b/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_;
diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp
index 5cb3d23a9..28f4293c5 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp
@@ -130,7 +130,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        [[gnu::always_inline]] inline int evaluateInternalResidual(
+        [[gnu::always_inline]] inline int evaluateInternalResidualKernel(
             const ScalarT* y,
             const ScalarT* yp,
             const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp
index 289e4465a..1cc425396 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp
@@ -849,7 +849,7 @@ namespace GridKit
         wb_[0] = Vr();
         wb_[1] = Vi();
 
-        evaluateInternalResidual(y_.getData(), yp_.getData(), wb_.data(), ws_.data(), f_.getData());
+        evaluateInternalResidualKernel(y_.getData(), yp_.getData(), wb_.data(), ws_.data(), f_.getData());
         f_.setDataUpdated();
         return 0;
       }
@@ -895,7 +895,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       [[gnu::always_inline]] inline int
-      Reecb<scalar_type, index_type>::evaluateInternalResidual(
+      Reecb<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp
index 73213d15d..b1ed62a47 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp
@@ -133,7 +133,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        [[gnu::always_inline]] inline int evaluateInternalResidual(
+        [[gnu::always_inline]] inline int evaluateInternalResidualKernel(
             const ScalarT* y,
             const ScalarT* yp,
             const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp
index 74a939ee1..fc1b4e623 100644
--- a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp
@@ -674,7 +674,7 @@ namespace GridKit
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
 
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         f_.setDataUpdated();
         return 0;
       }
@@ -723,7 +723,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       [[gnu::always_inline]] inline int
-      Repca<scalar_type, index_type>::evaluateInternalResidual(
+      Repca<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp
index 2bfb937ca..9306abf15 100644
--- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp
+++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp
@@ -121,7 +121,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT* y, const ScalarT* yp, const ScalarT* wb, const ScalarT* ws, ScalarT* f);
 
         __attribute__((always_inline)) inline int evaluateBusResidual(
diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp
index 456e5ef75..22f399771 100644
--- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp
@@ -611,7 +611,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       __attribute__((always_inline)) inline int
-      Regca<scalar_type, index_type>::evaluateInternalResidual(
+      Regca<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
@@ -754,7 +754,7 @@ namespace GridKit
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
 
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         evaluateBusResidual(y, yp, wb_.data(), h_.data());
         f_.setDataUpdated();
 
diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp
index 5d9e8f36d..7da10b527 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp
@@ -115,7 +115,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT* y,
             const ScalarT* yp,
             const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp
index 4789ad481..784993149 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp
@@ -491,7 +491,7 @@ namespace GridKit
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
 
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         f_.setDataUpdated();
         return 0;
       }
@@ -527,7 +527,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       __attribute__((always_inline)) inline int
-      Esdc1a<scalar_type, index_type>::evaluateInternalResidual(
+      Esdc1a<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp
index a86518cca..991db2071 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp
@@ -120,7 +120,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
       private:
diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
index 31f2760d5..f1d67176f 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp
@@ -333,7 +333,7 @@ namespace GridKit
        *
        */
       template <typename scalar_type, typename index_type>
-      __attribute__((always_inline)) inline int Ieeet1<scalar_type, index_type>::evaluateInternalResidual(
+      __attribute__((always_inline)) inline int Ieeet1<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
@@ -416,7 +416,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
 
         f_.setDataUpdated();
 
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp
index b1268515e..39d36b802 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp
@@ -105,7 +105,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
       private:
diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
index 8b8e53a11..4c2dc262b 100644
--- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp
@@ -211,7 +211,7 @@ namespace GridKit
       }
 
       template <typename scalar_type, typename index_type>
-      __attribute__((always_inline)) inline int SexsPti<scalar_type, index_type>::evaluateInternalResidual(
+      __attribute__((always_inline)) inline int SexsPti<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
@@ -253,7 +253,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
 
         f_.setDataUpdated();
 
diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp
index bc824576d..d4d60a569 100644
--- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp
@@ -100,7 +100,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        [[gnu::always_inline]] inline int evaluateInternalResidual(
+        [[gnu::always_inline]] inline int evaluateInternalResidualKernel(
             const ScalarT* y,
             const ScalarT* yp,
             const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp
index 8a58a875a..ff7491c6a 100644
--- a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp
@@ -480,7 +480,7 @@ namespace GridKit
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
 
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         f_.setDataUpdated();
         return 0;
       }
@@ -533,7 +533,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       [[gnu::always_inline]] inline int
-      GastPti<scalar_type, index_type>::evaluateInternalResidual(
+      GastPti<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT* y,
           const ScalarT* yp,
           const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp
index 9844dffd3..da7c8d146 100644
--- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp
@@ -115,7 +115,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT* y,
             const ScalarT* yp,
             const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp
index 9761ca9c6..3043e6335 100644
--- a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp
@@ -532,7 +532,7 @@ namespace GridKit
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
 
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         f_.setDataUpdated();
         return 0;
       }
@@ -567,7 +567,7 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       __attribute__((always_inline)) inline int
-      Hygov<scalar_type, index_type>::evaluateInternalResidual(
+      Hygov<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT*                  y,
           const ScalarT*                  yp,
           [[maybe_unused]] const ScalarT* wb,
diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp
index 89ce12994..ac1ade1cc 100644
--- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp
@@ -112,7 +112,7 @@ namespace GridKit
         }
 
       public:
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
       private:
diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
index e7210b462..93f44f54e 100644
--- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
+++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp
@@ -350,7 +350,7 @@ namespace GridKit
        *
        */
       template <typename scalar_type, typename index_type>
-      __attribute__((always_inline)) inline int Tgov1<scalar_type, index_type>::evaluateInternalResidual(
+      __attribute__((always_inline)) inline int Tgov1<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT*                  y,
           const ScalarT*                  yp,
           [[maybe_unused]] const ScalarT* wb,
@@ -412,7 +412,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
 
         f_.setDataUpdated();
 
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp
index e31259e19..c5bcbadab 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZ.hpp
@@ -116,7 +116,7 @@ namespace GridKit
     public:
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
     private:
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp
index f3708e491..9eafa8782 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZ/LoadZImpl.hpp
@@ -190,7 +190,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) int LoadZ<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) int LoadZ<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT*                  y,
         [[maybe_unused]] const ScalarT* yp,
         const ScalarT*                  wb,
@@ -218,7 +218,7 @@ namespace GridKit
       const auto* y  = y_.getData();
       const auto* yp = yp_.getData();
       auto*       f  = f_.getData();
-      evaluateInternalResidual(y, yp, wb_.data(), f);
+      evaluateInternalResidualKernel(y, yp, wb_.data(), f);
       evaluateBusResidual(y, yp, wb_.data(), h_.data());
       Ir() += h_[0];
       Ii() += h_[1];
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp
index 5387e6f27..13a9b90ef 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIP.hpp
@@ -129,7 +129,7 @@ namespace GridKit
     public:
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
     private:
diff --git a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp
index 58058b07c..73f1a7418 100644
--- a/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Load/LoadZIP/LoadZIPImpl.hpp
@@ -218,7 +218,7 @@ namespace GridKit
       const auto* y  = y_.getData();
       const auto* yp = yp_.getData();
       auto*       f  = f_.getData();
-      evaluateInternalResidual(y, yp, wb_.data(), f);
+      evaluateInternalResidualKernel(y, yp, wb_.data(), f);
       evaluateBusResidual(y, yp, wb_.data(), h_.data());
       Ir() += h_[0];
       Ii() += h_[1];
@@ -236,7 +236,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) int LoadZIP<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) int LoadZIP<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT*                  y,
         [[maybe_unused]] const ScalarT* yp,
         const ScalarT*                  wb,
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
index 0ed5893e6..aa7884cd3 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
@@ -112,7 +112,7 @@ namespace GridKit
 
         const Model::VariableMonitorBase* getMonitor() const override;
 
-        __attribute__((always_inline)) inline int evaluateInternalResidual(
+        __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
             const ScalarT*,
             const ScalarT*,
             const ScalarT*,
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
index 0b4454775..b8173a34f 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
@@ -298,7 +298,7 @@ namespace GridKit
       }
 
       template <typename scalar_type, typename index_type>
-      __attribute__((always_inline)) inline int Ieeest<scalar_type, index_type>::evaluateInternalResidual(
+      __attribute__((always_inline)) inline int Ieeest<scalar_type, index_type>::evaluateInternalResidualKernel(
           const ScalarT*                  y,
           const ScalarT*                  yp,
           [[maybe_unused]] const ScalarT* wb,
@@ -360,7 +360,7 @@ namespace GridKit
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
         auto*       f  = f_.getData();
-        evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+        evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
         f_.setDataUpdated();
 
         return 0;
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp
index 1dac5ac96..16377f820 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp
@@ -201,7 +201,7 @@ namespace GridKit
       }
 
     public:
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp
index 2169ea443..4fa3c3552 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp
@@ -518,7 +518,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) inline int Genrou<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) inline int Genrou<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT* y,
         const ScalarT* yp,
         const ScalarT* wb,
@@ -644,7 +644,7 @@ namespace GridKit
       const auto* y  = y_.getData();
       const auto* yp = yp_.getData();
       auto*       f  = f_.getData();
-      evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+      evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
       evaluateBusResidual(y, yp, wb_.data(), h_.data());
 
       // Genrou contribution to bus algebraic equations
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp
index 74655e7f1..72dd41800 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp
@@ -166,7 +166,7 @@ namespace GridKit
       }
 
     public:
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp
index 385eda4d6..c17dbd3a5 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp
@@ -380,7 +380,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) inline int Gensal<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) inline int Gensal<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT* y,
         const ScalarT* yp,
         const ScalarT* wb,
@@ -499,7 +499,7 @@ namespace GridKit
       const auto* y  = y_.getData();
       const auto* yp = yp_.getData();
       auto*       f  = f_.getData();
-      evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f);
+      evaluateInternalResidualKernel(y, yp, wb_.data(), ws_.data(), f);
       evaluateBusResidual(y, yp, wb_.data(), h_.data());
 
       // Gensal contribution to bus algebraic equations
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp
index 14ba8ce36..88ae60d4c 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp
@@ -145,7 +145,7 @@ namespace GridKit
       }
 
     public:
-      __attribute__((always_inline)) inline int evaluateInternalResidual(
+      __attribute__((always_inline)) inline int evaluateInternalResidualKernel(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
       __attribute__((always_inline)) inline int evaluateBusResidual(
           const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
index 5bd9ab0eb..dbc52660b 100644
--- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassicalImpl.hpp
@@ -280,7 +280,7 @@ namespace GridKit
      *
      */
     template <typename scalar_type, typename index_type>
-    __attribute__((always_inline)) int GenClassical<scalar_type, index_type>::evaluateInternalResidual(
+    __attribute__((always_inline)) int GenClassical<scalar_type, index_type>::evaluateInternalResidualKernel(
         const ScalarT* y,
         const ScalarT* yp,
         const ScalarT* wb,
@@ -350,7 +350,7 @@ namespace GridKit
       const auto* y  = y_.getData();
       const auto* yp = yp_.getData();
       auto*       f  = f_.getData();
-      evaluateInternalResidual(y, yp, wb_.data(), f);
+      evaluateInternalResidualKernel(y, yp, wb_.data(), f);
       evaluateBusResidual(y, yp, wb_.data(), h_.data());
 
       Ir() += h_[0];
diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp
index cfe031300..fd857ad82 100644
--- a/GridKit/Model/PhasorDynamics/SystemModel.hpp
+++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp
@@ -90,6 +90,8 @@ namespace GridKit
 
       int tagDifferentiable() override;
       int setAbsoluteTolerance(RealT rel_tol) override;
+      int evaluateInternalResidual() override;
+      int evaluateExternalResidual() override;
       int evaluateResidual() override;
       int evaluateJacobian() override;
 
diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp
index 8e640eac1..be45116c8 100644
--- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp
@@ -1062,20 +1062,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 adition. 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_)
       {
@@ -1084,9 +1077,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;
diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
index 02eba339d..added05fa 100644
--- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp
@@ -99,6 +99,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;
@@ -207,12 +228,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);