Eliminate useless state in Genrou and Gensal

GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp

@@ -48,7 +48,7 @@ namespace GridKit
         S12_(0.),
         mva_base_(100.)
     {
-      size_ = 19;
+      size_ = 17;
       setDerivedParams();
     }
 
@@ -97,7 +97,7 @@ namespace GridKit
         S12_(S12),
         mva_base_(100.)
     {
-      size_ = 19;
+      size_ = 17;
       setDerivedParams();
     }
 
@@ -112,7 +112,7 @@ namespace GridKit
       initializeParameters(data);
       initializeMonitor();
 
-      size_ = 19;
+      size_ = 17;
       setDerivedParams();
     }
 
@@ -129,7 +129,7 @@ namespace GridKit
       initializeParameters(data);
       initializeMonitor();
 
-      size_ = 19;
+      size_ = 17;
       setDerivedParams();
     }
 
@@ -147,7 +147,7 @@ namespace GridKit
       initializeParameters(data);
       initializeMonitor();
 
-      size_ = 19;
+      size_ = 17;
       setDerivedParams();
     }
 
@@ -451,10 +451,6 @@ namespace GridKit
       y[14]       = iq;
       y[15]       = ir;
       y[16]       = ii;
-      y[17]       = G_ * (vd * std::sin(delta) + vq * std::cos(delta))
-              - B_ * (vd * -std::cos(delta) + vq * std::sin(delta));
-      y[18] = B_ * (vd * std::sin(delta) + vq * std::cos(delta))
-              + G_ * (vd * -std::cos(delta) + vq * std::sin(delta));
 
       ScalarT Te = y[12];
       // Convert Te to system base for governor PM signal.
@@ -543,8 +539,6 @@ namespace GridKit
       ScalarT iq     = y[14];
       ScalarT ir     = y[15];
       ScalarT ii     = y[16];
-      ScalarT inr    = y[17];
-      ScalarT ini    = y[18];
 
       /* Read derivatives */
       ScalarT delta_dot = yp[0];
@@ -582,12 +576,8 @@ namespace GridKit
       f[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id);
       f[13] = id - (ir * std::sin(delta) - ii * std::cos(delta));
       f[14] = iq - (ir * std::cos(delta) + ii * std::sin(delta));
-      f[15] = ir + G_ * vr - B_ * vi - inr;
-      f[16] = ii + B_ * vr + G_ * vi - ini;
-
-      /* 2 Genrou current source definitions */
-      f[17] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq));
-      f[18] = ini - (B_ * (std::sin(delta) * vd + std::cos(delta) * vq) + G_ * (-std::cos(delta) * vd + std::sin(delta) * vq));
+      f[15] = ir - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq - vr) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq - vi));
+      f[16] = ii - (B_ * (std::sin(delta) * vd + std::cos(delta) * vq - vr) + G_ * (-std::cos(delta) * vd + std::sin(delta) * vq - vi));
 
       return 0;
     }

GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp

@@ -46,8 +46,6 @@ namespace GridKit
       IQ,     ///< \f$I_q\f$ q-axis current
       IR,     ///< \f$I_r\f$ network real current
       II,     ///< \f$I_i\f$ network imaginary current
-      INR,    ///< \f$I_{\mathrm{N},r}\f$ Norton source real current
-      INI,    ///< \f$I_{\mathrm{N},i}\f$ Norton source imaginary current
       MAXIMUM,
     };
 

GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp

@@ -25,7 +25,7 @@ namespace GridKit
       initializeParameters(data);
       initializeMonitor();
 
-      size_ = 16;
+      size_ = 14;
       setDerivedParams();
     }
 
@@ -314,10 +314,6 @@ namespace GridKit
       y[11] = iq;
       y[12] = ir;
       y[13] = ii;
-      y[14] = G_ * (vd * std::sin(delta) + vq * std::cos(delta))
-              - B_ * (vd * -std::cos(delta) + vq * std::sin(delta));
-      y[15] = B_ * (vd * std::sin(delta) + vq * std::cos(delta))
-              + G_ * (vd * -std::cos(delta) + vq * std::sin(delta));
 
       // Convert Te to system base for governor PM signal.
       pmech_set_ = toSystemBase(Te);
@@ -402,8 +398,6 @@ namespace GridKit
       ScalarT iq     = y[11];
       ScalarT ir     = y[12];
       ScalarT ii     = y[13];
-      ScalarT inr    = y[14];
-      ScalarT ini    = y[15];
 
       /* Read derivatives */
       ScalarT delta_dot  = yp[0];
@@ -437,12 +431,8 @@ namespace GridKit
       f[9]  = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id);
       f[10] = id - (ir * std::sin(delta) - ii * std::cos(delta));
       f[11] = iq - (ir * std::cos(delta) + ii * std::sin(delta));
-      f[12] = ir + G_ * vr - B_ * vi - inr;
-      f[13] = ii + B_ * vr + G_ * vi - ini;
-
-      /* 2 Gensal current source definitions */
-      f[14] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq));
-      f[15] = ini - (B_ * (std::sin(delta) * vd + std::cos(delta) * vq) + G_ * (-std::cos(delta) * vd + std::sin(delta) * vq));
+      f[12] = ir - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq - vr) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq - vi));
+      f[13] = ii - (B_ * (std::sin(delta) * vd + std::cos(delta) * vq - vr) + G_ * (-std::cos(delta) * vd + std::sin(delta) * vq - vi));
 
       return 0;
     }

examples/PhasorDynamics/Tiny/TwoBus/Basic/TwoBusBasicJson.cpp

@@ -78,32 +78,30 @@ int main(int argc, const char* argv[])
   // Get access to the fault
   auto* fault = sys.getBusFault(0);
 
-  // Set time step to 1/4 of a 60Hz cycle
-  real_type dt = 1.0 / 4.0 / 60.0;
+  // Monitor every quarter cycle at 60 Hz
+  real_type dt_monitor = 1.0 / 4.0 / 60.0;
 
   // Set up simulation
   Ida<scalar_type, size_t> ida(&sys);
+  ida.setMaxSteps(10000);
   ida.configureSimulation();
 
-  // Run simulation - making sure to pass the callback to record output
+  // Run simulation
   real_type start = static_cast<real_type>(clock());
 
   // Run for 1s
   ida.initializeSimulation(0.0, false);
-  int nout = static_cast<int>(std::round((1.0 - 0.0) / dt));
-  ida.runSimulation(1.0, nout);
+  ida.runSimulation(1.0, dt_monitor);
 
   // Introduce fault and run for the next 0.1s
   fault->setStatus(true);
   ida.initializeSimulation(1.0);
-  nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
-  ida.runSimulation(1.1, nout);
+  ida.runSimulation(1.1, dt_monitor);
 
   // Clear the fault and run until t = 10s.
   fault->setStatus(false);
   ida.initializeSimulation(1.1);
-  nout = static_cast<int>(std::round((10.0 - 1.1) / dt));
-  ida.runSimulation(10.0, nout);
+  ida.runSimulation(10.0, dt_monitor);
   real_type stop = static_cast<real_type>(clock());
 
   std::cout << "\n\nComplete in " << (stop - start) / CLOCKS_PER_SEC << " seconds\n";

tests/IntegrationTests/PhasorDynamics/PDIntegrationTests.hpp

@@ -243,6 +243,8 @@ namespace GridKit
 
         // Set up simulation
         AnalysisManager::Sundials::Ida<RealT, size_t> ida(&sys);
+        ida.setTolerance(static_cast<RealT>(1.0e-7),
+                         static_cast<RealT>(1.0e-9));
         ida.configureSimulation();
 
         // Run for 1s

tests/UnitTests/PhasorDynamics/GenrouTests.hpp

@@ -236,10 +236,8 @@ namespace GridKit
             2.211,
             0.85,
             1.2,
-            static_cast<ScalarT>(64.) / static_cast<ScalarT>(65.),
-            -static_cast<ScalarT>(237.) / static_cast<ScalarT>(130.),
-            -static_cast<ScalarT>(141.) / static_cast<ScalarT>(130.),
-            -static_cast<ScalarT>(241.) / static_cast<ScalarT>(260.)};
+            -static_cast<ScalarT>(13.) / static_cast<ScalarT>(130.),
+            -static_cast<ScalarT>(143.) / static_cast<ScalarT>(52.)};
 
         bus.allocate();
         bus.initialize();
@@ -273,8 +271,6 @@ namespace GridKit
         y[14] = .3;  // iq
         y[15] = .9;  // ir
         y[16] = .25; // ii
-        y[17] = .3;  // inr
-        y[18] = .15; // ini
 
         // Set derivative values matching the answer key
         yp[0] = 2.0 * pi * 60.0; // delta_dot

tests/UnitTests/PhasorDynamics/GensalTests.hpp

@@ -277,10 +277,8 @@ namespace GridKit
             1.8125,
             0.5,
             0.25,
-            2.65,
-            -0.05,
-            0.3,
-            -1.2};
+            2.95,
+            -1.25};
 
         bus.allocate();
         bus.initialize();
@@ -306,8 +304,6 @@ namespace GridKit
         y[11] = -0.5;  // iq
         y[12] = 0.75;  // ir
         y[13] = -0.25; // ii
-        y[14] = 0.1;   // inr
-        y[15] = -0.2;  // ini
 
         yp[0] = 2.0 * pi * 60.0; // delta_dot
         yp[1] = -1.0;            // omega_dot