Remove BusFault + Replace with Bus::setFault

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 - Added `BusSignalVoltageOut` bus model with voltage signal outlets and current signal inlets.
 - Added `BusSignalVoltageIn` bus model with voltage signal inlets and current signal outlets.
+- Added `Bus::setFault` to apply or clear a fault to ground directly at a bus.
 
 ## v0.2
 

GridKit/Model/PhasorDynamics/Bus/Bus.hpp

@@ -47,6 +47,7 @@ namespace GridKit
       virtual ~Bus();
 
       virtual int setBusID(IdxT) override final;
+      virtual int setFault(bool status, RealT R, RealT X) override final;
       virtual int allocate() override final;
       virtual int tagDifferentiable() override final;
       virtual int setAbsoluteTolerance(RealT rel_tol) override final;
@@ -160,6 +161,10 @@ namespace GridKit
     private:
       ScalarT Vr0_{0.0};
       ScalarT Vi0_{0.0};
+
+      /* Fault admittance */
+      RealT fault_g_{0.0};
+      RealT fault_b_{0.0};
     };
 
   } // namespace PhasorDynamics

GridKit/Model/PhasorDynamics/Bus/BusEnzyme.cpp

@@ -13,7 +13,7 @@ namespace GridKit
     /**
      * @brief Jacobian evaluation experimental.
      *
-     * This sets values to 0, and these remain unchanged. It is needed to get
+     * This sets values to the fault admittance, zero when cleared. It is needed to get
      * the indices into the list of entries that will later be deduplicated.
      * Contributions to bus Jacobians from other components are stored in those components.
      *
@@ -36,14 +36,15 @@ namespace GridKit
         J_cols_buffer_[1] = variable_indices_.at(1);
         J_cols_buffer_[2] = variable_indices_.at(0);
         J_cols_buffer_[3] = variable_indices_.at(1);
-        J_vals_buffer_[0] = 0.0;
-        J_vals_buffer_[1] = 0.0;
-        J_vals_buffer_[2] = 0.0;
-        J_vals_buffer_[3] = 0.0;
 
         nnz_ = 4;
         this->constructCoo();
       }
+
+      J_vals_buffer_[0] = -fault_g_;
+      J_vals_buffer_[1] = fault_b_;
+      J_vals_buffer_[2] = -fault_b_;
+      J_vals_buffer_[3] = -fault_g_;
       return 0;
     }
 

GridKit/Model/PhasorDynamics/Bus/BusImpl.hpp

@@ -136,6 +136,32 @@ namespace GridKit
       return 0;
     }
 
+    /**
+     * @brief Apply or clear a fault to ground
+     *
+     * @param[in] status - true applies the fault, false clears it
+     * @param[in] R - fault resistance [p.u.]
+     * @param[in] X - fault reactance [p.u.]
+     */
+    template <typename scalar_type, typename index_type>
+    int Bus<scalar_type, index_type>::setFault(bool status, RealT R, RealT X)
+    {
+      if (status && !(ZERO<RealT> < R * R + X * X))
+      {
+        Log::error() << "Bus: fault impedance R + jX must be nonzero\n";
+        return 1;
+      }
+
+      fault_g_ = 0.0;
+      fault_b_ = 0.0;
+      if (status)
+      {
+        fault_g_ = R / (X * X + R * R);
+        fault_b_ = -X / (X * X + R * R);
+      }
+      return 0;
+    }
+
     /*!
      * @brief Bus variables are algebraic.
      */
@@ -193,7 +219,7 @@ namespace GridKit
     }
 
     /*!
-     * @brief PQ bus does not compute residuals, so here we just reset residual values.
+     * @brief Reset the current balance to the fault current.
      *
      * @warning This implementation assumes bus residuals are always evaluated
      * _before_ component model residuals.
@@ -204,8 +230,8 @@ namespace GridKit
     {
       auto* f = f_.getData();
 
-      f[0] = 0.0;
-      f[1] = 0.0;
+      f[0] = -(fault_g_ * Vr() - fault_b_ * Vi());
+      f[1] = -(fault_b_ * Vr() + fault_g_ * Vi());
       f_.setDataUpdated();
       return 0;
     }

GridKit/Model/PhasorDynamics/Bus/README.md

@@ -69,12 +69,12 @@ None.
 
 #### Algebraic
 
-Let $\mathcal{D}$ denote the set of components connected to the bus.
+Let $\mathcal{D}$ denote the set of components connected to the bus, and let $G_f$, $B_f$ be the conductance and susceptance of a fault applied to the Bus, respectively, with $G_f + jB_f = 1/(R_f + jX_f)$, both zero when no fault is applied.
 
 ```math
 \begin{aligned}
-0 &= \sum_{d \in \mathcal{D}} I_{r,d} \\
-0 &= \sum_{d \in \mathcal{D}} I_{i,d}
+0 &= -(G_f V_r - B_f V_i) + \sum_{d \in \mathcal{D}} I_{r,d} \\
+0 &= -(B_f V_r + G_f V_i) + \sum_{d \in \mathcal{D}} I_{i,d}
 \end{aligned}
 ```
 

GridKit/Model/PhasorDynamics/BusBase.hpp

@@ -239,6 +239,13 @@ namespace GridKit
 
       virtual int setBusID(IdxT) = 0;
 
+      /// Apply (status true) or clear a fault to ground with impedance R + jX.
+      virtual int setFault(bool /* status */, RealT /* R */, RealT /* X */)
+      {
+        Log::error() << "Faults are not supported by this bus type\n";
+        return 1;
+      }
+
       virtual const IdxT busID() const
       {
         return bus_id_;

tests/IntegrationTests/PhasorDynamics/TenGen/Genrou/TenGenGenrou.cpp

@@ -71,8 +71,6 @@ int main()
   Genrou<scalar_type, index_type> gen9(&bus9, 0.5, -0.09662372, 3., 0., 0., 7., .04, .05, .75, 2.1, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.);
   Genrou<scalar_type, index_type> gen10(&bus10, 0.5, -0.09932297, 3., 0., 0., 7., .04, .05, .75, 2.1, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.);
 
-  BusFault<scalar_type, index_type> fault(&bus10, 0, 1e-5, 0);
-
   /* Connect everything together */
   SystemModel<scalar_type, index_type> sys;
 
@@ -104,7 +102,6 @@ int main()
   sys.addComponent(&gen8);
   sys.addComponent(&gen9);
   sys.addComponent(&gen10);
-  sys.addComponent(&fault);
   sys.allocate();
 
   real_type dt = 1.0 / 4.0 / 60.0;
@@ -171,7 +168,7 @@ int main()
   }
 
   // Introduce fault to ground and run for 0.1s
-  fault.setStatus(1);
+  bus10.setFault(true, 0.0, 1e-5);
   ida.initializeSimulation(1.0);
   ida.runSimulation(1.1, dt, output_cb);
 
@@ -183,7 +180,7 @@ int main()
   success                     *= isEqual(gen10.y().getData()[omega_index], omega_ref, 5e-5);
 
   // Clear fault and run until t = 10s.
-  fault.setStatus(0);
+  bus10.setFault(false, 0.0, 1e-5);
   ida.initializeSimulation(1.1);
   ida.runSimulation(10.0, dt, output_cb);
   real_type stop = static_cast<real_type>(clock());

tests/UnitTests/PhasorDynamics/BusTests.hpp

@@ -1,6 +1,10 @@
+#include <complex>
 #include <iomanip>
 #include <iostream>
+#include <vector>
 
+#include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
+#include <GridKit/Definitions.hpp>
 #include <GridKit/Model/PhasorDynamics/Bus/Bus.hpp>
 #include <GridKit/Model/PhasorDynamics/Bus/BusInfinite.hpp>
 #include <GridKit/Testing/TestHelpers.hpp>
@@ -13,6 +17,9 @@ namespace GridKit
     template <class ScalarT, typename IdxT>
     class BusTests
     {
+    private:
+      using RealT = typename PhasorDynamics::Bus<ScalarT, IdxT>::RealT;
+
     public:
       BusTests()  = default;
       ~BusTests() = default;
@@ -96,6 +103,113 @@ namespace GridKit
 
         return success.report(__func__);
       }
+
+      /// Fault current applied and cleared through setFault
+      TestOutcome fault()
+      {
+        TestStatus success = true;
+
+        ScalarT Vr{1.0};
+        ScalarT Vi{2.0};
+        RealT   R{1.0};
+        RealT   X{2.0};
+
+        PhasorDynamics::Bus<ScalarT, IdxT> bus(Vr, Vi);
+        bus.allocate();
+        bus.initialize();
+
+        // Fault current I = -V / (R + jX)
+        const std::complex<RealT> current = -std::complex<RealT>(Vr, Vi) / std::complex<RealT>(R, X);
+
+        success *= bus.setFault(true, R, X) == 0;
+        bus.evaluateResidual();
+        success *= isEqual(bus.Ir(), current.real());
+        success *= isEqual(bus.Ii(), current.imag());
+
+        success *= bus.setFault(false, R, X) == 0;
+        bus.evaluateResidual();
+        success *= isEqual(bus.Ir(), 0.0);
+        success *= isEqual(bus.Ii(), 0.0);
+
+        // A zero fault impedance is rejected
+        success *= bus.setFault(true, 0.0, 0.0) != 0;
+
+        // An infinite bus cannot be faulted
+        PhasorDynamics::BusInfinite<ScalarT, IdxT> bus_inf;
+        success *= bus_inf.setFault(true, R, X) != 0;
+
+        return success.report(__func__);
+      }
+
+#ifdef GRIDKIT_ENABLE_ENZYME
+      /// Fault Jacobian matches DependencyTracking and keeps its pattern when cleared
+      TestOutcome jacobian()
+      {
+        TestStatus success = true;
+
+        RealT R{1.0};
+        RealT X{2.0};
+
+        // Jacobian via DependencyTracking, which numbers y entries 2 * index
+        DependencyTracking::Variable Vr{1.0};
+        DependencyTracking::Variable Vi{2.0};
+
+        PhasorDynamics::Bus<DependencyTracking::Variable, IdxT> dependency_tracking_bus(Vr, Vi);
+        dependency_tracking_bus.allocate();
+        dependency_tracking_bus.initialize();
+        dependency_tracking_bus.setFault(true, R, X);
+        dependency_tracking_bus.evaluateResidual();
+
+        const auto* f    = dependency_tracking_bus.getResidual().getData();
+        const auto  size = static_cast<size_t>(dependency_tracking_bus.size());
+
+        std::vector<DependencyTracking::Variable::DependencyMap> dependency_tracking_jacobian(size);
+        for (size_t row = 0; row < size; ++row)
+        {
+          for (const auto& [number, value] : f[row].getDependencies())
+          {
+            dependency_tracking_jacobian[row][number / 2] = value;
+          }
+        }
+
+        // Jacobian from the bus
+        PhasorDynamics::Bus<ScalarT, IdxT> bus(1.0, 2.0);
+        bus.allocate();
+        bus.initialize();
+        bus.setFault(true, R, X);
+        bus.evaluateResidual();
+        bus.evaluateJacobian();
+
+        auto*        jacobian = bus.getCooJacobian();
+        const IdxT*  rows     = jacobian->getRowData();
+        const IdxT*  cols     = jacobian->getColData();
+        const RealT* vals     = jacobian->getValues();
+
+        std::vector<DependencyTracking::Variable::DependencyMap> bus_jacobian(size);
+        for (IdxT i = 0; i < jacobian->getNnz(); ++i)
+        {
+          bus_jacobian[static_cast<size_t>(rows[i])][static_cast<size_t>(cols[i])] = vals[i];
+        }
+
+        for (size_t row = 0; row < size; ++row)
+        {
+          success *= isEqual(dependency_tracking_jacobian[row], bus_jacobian[row]);
+        }
+
+        // Clearing the fault keeps the same entries with zero values
+        const IdxT nnz = jacobian->getNnz();
+        bus.setFault(false, R, X);
+        bus.evaluateResidual();
+        bus.evaluateJacobian();
+        success *= bus.nnz() == nnz;
+        for (IdxT i = 0; i < nnz; ++i)
+        {
+          success *= isEqual(vals[i], 0.0);
+        }
+
+        return success.report(__func__);
+      }
+#endif
     };
 
   } // namespace Testing

tests/UnitTests/PhasorDynamics/CMakeLists.txt

@@ -5,7 +5,11 @@
 
 add_executable(test_phasor_bus runBusTests.cpp)
 target_link_libraries(
-  test_phasor_bus GridKit::phasor_dynamics_bus GridKit::testing)
+  test_phasor_bus
+  GridKit::definitions
+  GridKit::phasor_dynamics_bus
+  GridKit::phasor_dynamics_bus_dependency_tracking
+  GridKit::testing)
 
 add_executable(test_phasor_bus_fault runBusFaultTests.cpp)
 target_link_libraries(

tests/UnitTests/PhasorDynamics/runBusTests.cpp

@@ -10,6 +10,10 @@ int main()
 
   result += test.constructor();
   result += test.residual();
+  result += test.fault();
+#ifdef GRIDKIT_ENABLE_ENZYME
+  result += test.jacobian();
+#endif
 
   return result.summary();
 }