Unit test and example for Genrou class.

examples/DistributedGeneratorTest/DGTest.cpp

@@ -1,5 +1,5 @@
 
-
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/Grid3Bus/Grid3BusSys.cpp

@@ -9,6 +9,7 @@
  *
  */
 
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/Microgrid/Microgrid.cpp

@@ -1,5 +1,6 @@
 
 
+#define _USE_MATH_DEFINES
 #include <cmath>
 #include <filesystem>
 #include <fstream>

examples/PhasorDynamics/CMakeLists.txt

@@ -1 +1,2 @@
 add_subdirectory(Example1)
+add_subdirectory(Example2)

examples/PhasorDynamics/Example1/CMakeLists.txt

@@ -4,11 +4,7 @@ target_link_libraries(phasordynamics_example1
 	GRIDKIT::phasor_dynamics_bus_fault 
 	GRIDKIT::phasor_dynamics_branch
 	GRIDKIT::phasor_dynamics_genrou
-	SUNDIALS::sunlinsolklu
-	SUNDIALS::core
-	SUNDIALS::ida
-	SUNDIALS::idas
-	SUNDIALS::sunmatrixdense)
+	GRIDKIT::solvers_dyn)
 install(TARGETS phasordynamics_example1 RUNTIME DESTINATION bin)
 
-add_test(NAME GenrouTest COMMAND $<TARGET_FILE:phasordynamics_example1>)
+add_test(NAME GenrouTest1 COMMAND $<TARGET_FILE:phasordynamics_example1>)

examples/PhasorDynamics/Example1/example1.cpp

@@ -1,29 +1,18 @@
-#include <stdio.h>
-#define _USE_MATH_DEFINES
-#include <math.h>
+#include "example1.hpp"
+
+#include <cmath>
 #include <time.h>
 
-// #include <sundials_core.h>
-#include <idas/idas.h>
-#include <nvector/nvector_serial.h>
-#include <sunlinsol/sunlinsol_dense.h>
-#include <sunlinsol/sunlinsol_klu.h>
-#include <sunmatrix/sunmatrix_sparse.h>
-
-#include "Example1_Powerworld_Reference.hpp"
-#include "Model/PhasorDynamics/Branch/Branch.cpp"
-#include "Model/PhasorDynamics/Branch/Branch.hpp"
-#include "Model/PhasorDynamics/Bus/Bus.cpp"
-#include "Model/PhasorDynamics/Bus/Bus.hpp"
-#include "Model/PhasorDynamics/Bus/BusInfinite.cpp"
-#include "Model/PhasorDynamics/Bus/BusInfinite.hpp"
-#include "Model/PhasorDynamics/BusFault/BusFault.hpp"
-#include "Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp"
-#include "Model/PhasorDynamics/SystemModel.hpp"
-#include "Solver/Dynamic/Ida.cpp"
-#include "Solver/Dynamic/Ida.hpp"
+#include <Model/PhasorDynamics/Branch/Branch.hpp>
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFault.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
+#include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
 #include <Utilities/Testing.hpp>
 
+#define _USE_MATH_DEFINES
 #define _CRT_SECURE_NO_WARNINGS
 
 int main()

examples/PhasorDynamics/Example1/Example1_Powerworld_Reference.hpp → examples/PhasorDynamics/Example1/example1.hpp

@@ -1,3 +1,10 @@
+/**
+ * @file example1.hpp
+ * @author Adam Birchfield ([email protected])
+ * @author Slaven Peles ([email protected])
+ * @brief Reference solution for the 2-bus system obtained with Powerworld
+ *
+ */
 #include <vector>
 
 // Columns:

examples/PhasorDynamics/Example2/CMakeLists.txt

@@ -0,0 +1,11 @@
+add_executable(phasordynamics_example2 example2.cpp)
+target_link_libraries(phasordynamics_example2 
+	GRIDKIT::phasor_dynamics_bus 
+	GRIDKIT::phasor_dynamics_bus_fault 
+	GRIDKIT::phasor_dynamics_branch
+	GRIDKIT::phasor_dynamics_genrou
+	GRIDKIT::phasor_dynamics_load
+	GRIDKIT::solvers_dyn)
+install(TARGETS phasordynamics_example2 RUNTIME DESTINATION bin)
+
+add_test(NAME GenrouTest2 COMMAND $<TARGET_FILE:phasordynamics_example2>)

examples/PhasorDynamics/Example2/example2.cpp

@@ -0,0 +1,125 @@
+#include "example2.hpp"
+
+#include <cmath>
+#include <cstdio>
+#include <time.h>
+#include <vector>
+
+#include <Model/PhasorDynamics/Branch/Branch.hpp>
+#include <Model/PhasorDynamics/Bus/Bus.hpp>
+#include <Model/PhasorDynamics/Bus/BusInfinite.hpp>
+#include <Model/PhasorDynamics/BusFault/BusFault.hpp>
+#include <Model/PhasorDynamics/Load/Load.hpp>
+#include <Model/PhasorDynamics/SynchronousMachine/GENROUwS/Genrou.hpp>
+#include <Model/PhasorDynamics/SystemModel.hpp>
+#include <Solver/Dynamic/Ida.hpp>
+#include <Utilities/Testing.hpp>
+
+#define _CRT_SECURE_NO_WARNINGS
+#define _USE_MATH_DEFINES
+
+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);
+
+  std::cout << "Example 2 version 1\n";
+
+  /* Create model parts */
+  SystemModel<scalar_type, index_type> sys;
+  BusInfinite<scalar_type, index_type> bus1(1.06, 0.0);
+  Bus<scalar_type, index_type>         bus2(1.0599558398065716, -0.009675621941024773);
+  Bus<scalar_type, index_type>         bus3(0.9610827543495831, -0.13122476630506485);
+  Branch<scalar_type, index_type>      branch12(&bus1, &bus2, 0.05, 0.21, 0, 0.1);
+  Branch<scalar_type, index_type>      branch13(&bus1, &bus3, 0.06, 0.15, 0, 0.12);
+  Branch<scalar_type, index_type>      branch23(&bus2, &bus3, 0.08, 0.27, 0, 0.45);
+  Genrou<scalar_type, index_type>      gen2(&bus2, 1, 0.5, -0.07588, 2.7, 0., 0., 7., .04, .05, .75, 1.9, 0.17, 0.15, 0.4, 0.35, 0.15, 0.14999, 0., 0.);
+  Genrou<scalar_type, index_type>      gen3(&bus3, 1, 0.25, 0.26587, 1.6, 0., 0., 7.5, .04, .05, .75, 2.3, 0.2, 0.18, 0.5, 0.5, 0.18, 0.15, 0., 0.);
+  Load<scalar_type, index_type>        load3(&bus3, 0.4447197839297772, 0.20330047265361242);
+  BusFault<scalar_type, index_type>    fault(&bus3, 0, 1e-5, 0);
+
+  /* Connect everything together */
+  sys.addBus(&bus1);
+  sys.addBus(&bus2);
+  sys.addBus(&bus3);
+  sys.addComponent(&branch12);
+  sys.addComponent(&branch13);
+  sys.addComponent(&branch23);
+  sys.addComponent(&fault);
+  sys.addComponent(&load3);
+  sys.addComponent(&gen2);
+  sys.addComponent(&gen3);
+  sys.allocate();
+
+  real_type dt = 1.0 / 4.0 / 60.0;
+
+  std::stringstream buffer;
+
+  /* Set up simulation */
+  Ida<scalar_type, index_type> ida(&sys);
+  ida.configureSimulation();
+
+  /* Run simulation */
+  scalar_type start = static_cast<scalar_type>(clock());
+  // ida.printOutputF(0, 0, buffer);
+  ida.initializeSimulation(0.0, false);
+  ida.runSimulationFixed(0.0, dt, 1.0, buffer);
+  fault.setStatus(1);
+  ida.initializeSimulation(1.0, false);
+  ida.runSimulationFixed(1.0, dt, 1.1, buffer);
+  fault.setStatus(0);
+  ida.initializeSimulation(1.1, false);
+  ida.runSimulationFixed(1.1, dt, 10.0, buffer);
+  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);
+  FILE* f = fopen("example2_results.csv", "w");
+  fprintf(f, "Time,gen2speed,gen3speed,v2mag,v3mag\n");
+  fprintf(f, "0,1,1,1,1\n");
+  for (index_type i = 0; i < nt - 1; ++i)
+  {
+    for (index_type j = 0; j < stride; ++j)
+      buffer >> results[j];
+    // for (j = 0; j < stride; ++j) printf("%d %d %g\n", i, j, results[j]);
+    real_type   t             = results[0];
+    scalar_type tref          = reference_solution[i + 1][0];
+    // printf("Time GridKit %g PowerWorld %g\n", t, tref);
+    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];
+    fprintf(f, "%g,%g,%g,%g,%g\n", t, gen2speed, gen3speed, v2mag, v3mag);
+    real_type err = std::max(std::max(std::abs(gen2speed - gen2speed_ref),
+                                      std::abs(gen3speed - gen3speed_ref)),
+                             std::max(std::abs(v2mag - v2mag_ref),
+                                      std::abs(v3mag - v3mag_ref)));
+    if (err > worst_error)
+    {
+      worst_error      = err;
+      worst_error_time = t;
+    }
+  }
+  fclose(f);
+
+  std::cout << "Worst error " << worst_error
+            << " at time t = " << worst_error_time << "\n";
+  std::cout << "\n\nComplete in " << (stop - start) / CLOCKS_PER_SEC << " seconds\n";
+
+  return worst_error < error_allowed ? 0 : 1;
+}