Add IDA interface callbacks

examples/PhasorDynamics/Example1/example1.cpp

@@ -67,83 +67,91 @@ int main()
 
   double dt = 1.0 / 4.0 / 60.0;
 
-  std::stringstream buffer;
+  // A data structure to keep track of the data we want to
+  // compare to the reference solution. Rather than keeping
+  // the entire solution vector at every time step around,
+  // we instead narrow down exactly what we want to keep.
+  //
+  // Since this struct is "simple" enough (no constructors or
+  // assignment operators, and "simple" members), it is a POD
+  // (plain ol' data), which have some benefits in C++.
+  struct OutputData
+  {
+    double ti, Vr, Vi, dw;
+  };
+
+  // A list of output for each time step.
+  std::vector<OutputData> output;
+
+  // A callback which will be called by the integrator after
+  // each time step. It will be told the time of the current
+  // state, and it is allowed to access the up-to-date state
+  // of the components, which are captured by a closure
+  // due to the [&] notation (every variable that is referenced
+  // by the callback that is external to the callback itself -
+  // here output, bus1, and gen - will be considered a
+  // reference to that variable inside the callback). We select
+  // the subset of the output we're interested in recording and
+  // push it into output, which is updated outside the callback.
+  auto output_cb = [&](double t)
+  {
+    std::vector<double>& yval = sys.y();
 
-  /* Set up simulation */
+    output.push_back(OutputData{t, yval[0], yval[1], yval[3]});
+  };
+
+  // Set up simulation
   Ida<double, size_t> ida(&sys);
   ida.configureSimulation();
 
-  /* Run simulation */
+  // Run simulation - making sure to pass the callback to record output
   double start = static_cast<double>(clock());
-  // ida.printOutputF(0, 0, buffer);
+
+  // Run for 1s
   ida.initializeSimulation(0.0, false);
-  ida.runSimulationFixed(0.0, dt, 1.0, buffer);
+  int nout = static_cast<int>(std::round((1.0 - 0.0) / dt));
+  ida.runSimulation(1.0, nout, output_cb);
+
+  // Introduce fault and run for the next 0.1s
   fault.setStatus(1);
   ida.initializeSimulation(1.0, false);
-  ida.runSimulationFixed(1.0, dt, 1.1, buffer);
+  nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
+  ida.runSimulation(1.1, nout, output_cb);
+
+  // Clear the fault and run until t = 10s.
   fault.setStatus(0);
   ida.initializeSimulation(1.1, false);
-  ida.runSimulationFixed(1.1, dt, 10.0, buffer);
+  nout = static_cast<int>(std::round((10.0 - 1.1) / dt));
+  ida.runSimulation(10.0, nout, output_cb);
   double stop = static_cast<double>(clock());
 
-  // Go to the beginning of the data buffer
-  buffer.seekg(0, std::ios::beg);
-
-  double data;
-
-  size_t i       = 0;   // data row counter
-  size_t j       = 0;   // data column counter
-  double Vr      = 0.0; // Bus real voltage
-  double Vi      = 0.0; // Bus imaginary voltage
-  double dw      = 0.0; // Generator frequency deviation [rad/s]
-  double ti      = 0.0; // time
   double error_V = 0.0; // error in |V|
 
-  // Read through the simulation data storred in the buffer
-  while (buffer >> data)
+  // Read through the simulation data stored in the buffer.
+  // Since we captured by reference, output should be available
+  // for us to read here, outside the callback.
+  for (size_t i = 0; i < output.size(); i++)
   {
-    // At the end of each data line compare computed data to Powerworld results
-    // and reset column counter to zero.
-    if ((i % 48) == 0)
-    {
-      double err =
-          std::abs(std::sqrt(Vr * Vr + Vi * Vi) - reference_solution[i / 48][2])
-          / (1.0 + std::abs(reference_solution[i / 48][2]));
-      if (err > error_V)
-        error_V = err;
-      std::cout << "GridKit: t = " << ti
-                << ", |V| = " << std::sqrt(Vr * Vr + Vi * Vi)
-                << ", w = " << (1.0 + dw) << "\n";
-      std::cout << "Ref    : t = " << reference_solution[i / 48][0]
-                << ", |V| = " << reference_solution[i / 48][2]
-                << ", w = " << reference_solution[i / 48][1]
-                << "\n";
-      std::cout << "Error in |V| = "
-                << err
-                << "\n";
-      j  = 0;
-      Vr = 0.0;
-      Vi = 0.0;
-      std::cout << "\n";
-    }
-    if (j == 0)
-    {
-      ti = data;
-    }
-    if (j == 2)
-    {
-      Vr = data;
-    }
-    if (j == 3)
-    {
-      Vi = data;
-    }
-    if (j == 5)
-    {
-      dw = data;
-    }
-    ++j;
-    ++i;
+    OutputData           data    = output[i];
+    std::vector<double>& ref_sol = reference_solution[i + 1];
+
+    double err =
+        std::abs(std::sqrt(data.Vr * data.Vr + data.Vi * data.Vi) - ref_sol[2])
+        / (1.0 + std::abs(ref_sol[2]));
+    if (err > error_V)
+      error_V = err;
+
+    std::cout << "GridKit: t = " << data.ti
+              << ", |V| = " << std::sqrt(data.Vr * data.Vr + data.Vi * data.Vi)
+              << ", w = " << (1.0 + data.dw) << "\n";
+    std::cout << "Ref    : t = " << ref_sol[0]
+              << ", |V| = " << ref_sol[2]
+              << ", w = " << ref_sol[1]
+              << "\n";
+    std::cout << "Error in |V| = "
+              << err
+              << "\n";
+    std::cout << "\n";
   }
 
   int status = 0;

examples/PhasorDynamics/Example2/example2.cpp

@@ -25,13 +25,58 @@
 #include <Solver/Dynamic/Ida.hpp>
 #include <Utilities/Testing.hpp>
 
+using scalar_type = double;
+using real_type   = double;
+using index_type  = size_t;
+
+struct OutputData
+{
+  real_type   t;
+  scalar_type gen2speed;
+  scalar_type gen3speed;
+  scalar_type v2mag;
+  scalar_type v3mag;
+
+  OutputData& operator-=(const OutputData& other)
+  {
+    assert(t == other.t);
+    gen2speed -= other.gen2speed;
+    gen3speed -= other.gen3speed;
+    v2mag     -= other.v2mag;
+    v3mag     -= other.v3mag;
+    return *this;
+  }
+
+  double norm() const
+  {
+    return std::max({
+        std::abs(gen2speed),
+        std::abs(gen3speed),
+        std::abs(v2mag),
+        std::abs(v3mag),
+    });
+  }
+};
+
+const OutputData operator-(const OutputData& lhs, const OutputData& rhs)
+{
+  return OutputData(lhs) -= rhs;
+}
+
+std::ostream& operator<<(std::ostream& out, const OutputData& data)
+{
+  out << data.t << ","
+      << data.gen2speed << ","
+      << data.gen3speed << ","
+      << data.v2mag << ","
+      << data.v3mag;
+  return out;
+}
+
 int main()
 {
   using namespace GridKit::PhasorDynamics;
   using namespace AnalysisManager::Sundials;
-  using scalar_type = double;
-  using real_type   = double;
-  using index_type  = size_t;
 
   auto error_allowed = static_cast<real_type>(1e-4);
 
@@ -65,33 +110,48 @@ int main()
 
   real_type dt = 1.0 / 4.0 / 60.0;
 
-  std::stringstream buffer;
+  std::vector<OutputData> output;
+
+  auto output_cb = [&](real_type t)
+  {
+    std::vector<double>& yval = sys.y();
 
-  /* Set up simulation */
+    output.push_back(OutputData{t,
+                                1.0 + yval[5],
+                                1.0 + yval[26],
+                                std::sqrt(yval[0] * yval[0] + yval[1] * yval[1]),
+                                std::sqrt(yval[2] * yval[2] + yval[3] * yval[3])});
+  };
+
+  // Set up simulation
   Ida<scalar_type, index_type> ida(&sys);
   ida.configureSimulation();
 
-  /* Run simulation */
+  // Run simulation, output each `dt` interval
   scalar_type start = static_cast<scalar_type>(clock());
   ida.initializeSimulation(0.0, false);
-  ida.runSimulationFixed(0.0, dt, 1.0, buffer);
+
+  // Run for 1s
+  int nout = static_cast<int>(std::round((1.0 - 0.0) / dt));
+  ida.runSimulation(1.0, nout, output_cb);
+
+  // Introduce fault to ground and run for 0.1s
   fault.setStatus(1);
   ida.initializeSimulation(1.0, false);
-  ida.runSimulationFixed(1.0, dt, 1.1, buffer);
+  nout = static_cast<int>(std::round((1.1 - 1.0) / dt));
+  ida.runSimulation(1.1, nout, output_cb);
+
+  // Clear fault and run until t = 10s.
   fault.setStatus(0);
   ida.initializeSimulation(1.1, false);
-  ida.runSimulationFixed(1.1, dt, 10.0, buffer);
+  nout = static_cast<int>(std::round((10.0 - 1.1) / dt));
+  ida.runSimulation(10.0, nout, output_cb);
   double stop = static_cast<double>(clock());
 
   /* Check worst-case error */
   real_type worst_error      = 0;
   real_type worst_error_time = 0;
 
-  const index_type stride = 94;
-  const index_type nt     = 2401;
-  scalar_type      results[stride];
-  buffer.seekg(0, std::ios::beg);
-
   std::ostream  nullout(nullptr);
   std::ostream& out = nullout;
 
@@ -103,32 +163,22 @@ int main()
   out << "Time,gen2speed,gen3speed,v2mag,v3mag\n";
   out << 0. << "," << 1. << "," << 1. << "," << 1. << "," << 1. << "\n";
 
-  for (index_type i = 0; i < nt - 1; ++i)
+  for (index_type i = 0; i < output.size(); ++i)
   {
-    for (index_type j = 0; j < stride; ++j)
-    {
-      buffer >> results[j];
-    }
-    real_type   t             = results[0];
-    scalar_type gen2speed     = 1 + results[7];
-    scalar_type gen2speed_ref = reference_solution[i + 1][1];
-    scalar_type gen3speed     = 1 + results[28];
-    scalar_type gen3speed_ref = reference_solution[i + 1][2];
-    scalar_type v2mag         = sqrt(results[2] * results[2] + results[3] * results[3]);
-    scalar_type v2mag_ref     = reference_solution[i + 1][4];
-    scalar_type v3mag         = sqrt(results[4] * results[4] + results[5] * results[5]);
-    scalar_type v3mag_ref     = reference_solution[i + 1][5];
-
-    out << t << "," << gen2speed << "," << gen3speed << "," << v2mag << "," << v3mag << "\n";
-
-    real_type err = std::max({std::abs(gen2speed - gen2speed_ref),
-                              std::abs(gen3speed - gen3speed_ref),
-                              std::abs(v2mag - v2mag_ref),
-                              std::abs(v3mag - v3mag_ref)});
+    OutputData ref{reference_solution[i + 1][0],
+                   reference_solution[i + 1][1],
+                   reference_solution[i + 1][2],
+                   reference_solution[i + 1][4],
+                   reference_solution[i + 1][5]};
+    OutputData out_data = output[i];
+
+    out << out_data << '\n';
+
+    real_type err = (out_data - ref).norm();
     if (err > worst_error)
     {
       worst_error      = err;
-      worst_error_time = t;
+      worst_error_time = out_data.t;
     }
   }
   // fileout.close();