Initial implementation of generic phasor dynamics app

CHANGELOG.md

@@ -47,6 +47,7 @@
 - Refactored Jacobian assembly in `PowerElectronics` module to reuse the CSR pattern.
 - Refactored Jacobian assembly in `PhasorDyanmcics` module to reuse the CSR pattern.
 - Removed `COO_Matrix` class use in `PowerElectronics` module.
+- Added phasor dynamics application to generalize examples
 
 ## v0.1
 

CMakeLists.txt

@@ -16,9 +16,9 @@ set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
 if (MSVC)
-  set(CMAKE_CXX_FLAGS "/Wall")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /Wall")
 else()
-  set(CMAKE_CXX_FLAGS "-Wall -Wextra -Wconversion -Wpedantic")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wconversion -Wpedantic")
 endif()
 
 set(GRIDKIT_THIRD_PARTY_DIR ${PROJECT_SOURCE_DIR}/third-party)
@@ -111,6 +111,9 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
 add_subdirectory(GridKit)
 
+# Create applications
+add_subdirectory(application)
+
 # Create examples and tests
 include(CTest)
 add_subdirectory(examples)

GridKit/Model/VariableMonitor.hpp

@@ -79,7 +79,7 @@ namespace GridKit
         /// Output format
         Format      format;
         /// Delimiter (used only with CSV format currently)
-        std::string delim;
+        std::string delim{","};
       };
 
       virtual ~VariableMonitorBase()

GridKit/Testing/Testing.hpp

@@ -55,6 +55,16 @@ namespace GridKit
         return *this;
       }
 
+      operator bool() const
+      {
+        return outcome_ == TestOutcome::PASS;
+      }
+
+      int get() const
+      {
+        return outcome_;
+      }
+
       void skipTest()
       {
         outcome_ = TestOutcome::SKIP;

application/CMakeLists.txt

@@ -0,0 +1,3 @@
+if(TARGET SUNDIALS::idas)
+  add_subdirectory(PhasorDynamics)
+endif()

application/PhasorDynamics/CMakeLists.txt

@@ -0,0 +1,12 @@
+add_executable(PDSim PDSim.cpp)
+target_link_libraries(PDSim
+    PUBLIC
+    GridKit::phasor_dynamics_components
+    GridKit::solvers_dyn
+    GridKit::Utilities
+    GridKit::testing)
+target_include_directories(PDSim PRIVATE
+    ${GRIDKIT_THIRD_PARTY_DIR}/nlohmann-json/include
+    ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include)
+
+install(TARGETS PDSim EXPORT gridkit-targets RUNTIME)

application/PhasorDynamics/PDSim.cpp

@@ -0,0 +1,106 @@
+#include "PDSim.hpp"
+
+#include <filesystem>
+#include <fstream>
+
+#include <GridKit/Model/PhasorDynamics/SystemModel.hpp>
+#include <GridKit/Solver/Dynamic/Ida.hpp>
+#include <GridKit/Testing/TestHelpers.hpp>
+#include <GridKit/Testing/Testing.hpp>
+
+using Log = GridKit::Utilities::Logger;
+
+using namespace GridKit::PhasorDynamics;
+using namespace GridKit::Testing;
+using namespace AnalysisManager::Sundials;
+
+using scalar_type = double;
+using real_type   = double;
+using index_type  = size_t;
+
+int main(int argc, const char* argv[])
+{
+  // Study file
+  if (argc < 2)
+  {
+    Log::error() << "No input file provided" << std::endl;
+    std::cout << "\n"
+                 "Usage:\n"
+                 "       pdsim <json-input-file>\n"
+                 "\n"
+                 "Please provide a json input file for the study to run.\n"
+                 "\n";
+    exit(1);
+  }
+
+  auto study = parseStudyData(argv[1]);
+
+  // Instantiate system
+  SystemModel<scalar_type, index_type> sys(study.model_data);
+  sys.allocate();
+
+  real_type dt = study.dt;
+
+  // Set up simulation
+  Ida<scalar_type, index_type> ida(&sys);
+  ida.configureSimulation();
+
+  // Start timer
+  real_type start = static_cast<real_type>(clock());
+
+  using EventType = SystemEvent::Type;
+
+  // Initilize simultation for first run
+  ida.initializeSimulation(0.0, false);
+  real_type curr_time = 0.0;
+  for (const auto& event : study.events)
+  {
+    // Run to event time
+    int nout = static_cast<int>(std::round((event.time - curr_time) / dt));
+    ida.runSimulation(event.time, nout);
+
+    // Set up run for event (to start at event time)
+    if (event.type == EventType::FAULT_ON)
+    {
+      sys.getBusFault(event.element_id)->setStatus(true);
+    }
+    else if (event.type == EventType::FAULT_OFF)
+    {
+      sys.getBusFault(event.element_id)->setStatus(false);
+    }
+
+    // Re-initialize simulation at event time
+    ida.initializeSimulation(event.time, false);
+    curr_time = event.time;
+  }
+
+  // Run to final time
+  int nout = static_cast<int>(std::round((study.tmax - curr_time) / dt));
+  ida.runSimulation(study.tmax, nout);
+
+  real_type stop = static_cast<real_type>(clock());
+
+  // Stop the variable monitor
+  sys.stopMonitor();
+
+  // Generate aggregate errors comparing variable output to reference solution
+  std::string func{"monitor file vs reference file"};
+  TestStatus  status{func.c_str()};
+  if (!study.output_file.empty() && !study.reference_file.empty())
+  {
+    auto errorSet = compareCSV(study.output_file, study.reference_file);
+
+    // Print the errors
+    errorSet.display();
+
+    // Check against specified tolerance
+    status *= errorSet.total.max_error < study.error_tol;
+
+    status.report();
+  }
+
+  // Report run time
+  std::cout << "\n\nComplete in " << (stop - start) / CLOCKS_PER_SEC << " seconds\n";
+
+  return status.get();
+}

application/PhasorDynamics/PDSim.hpp

@@ -0,0 +1,188 @@
+#pragma once
+
+#include <filesystem>
+#include <fstream>
+#include <vector>
+
+#include <magic_enum/magic_enum.hpp>
+#include <nlohmann/json.hpp>
+
+#include <GridKit/Model/PhasorDynamics/SystemModelData.hpp>
+#include <GridKit/Utilities/Logger/Logger.hpp>
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    namespace fs = ::std::filesystem;
+
+    /**
+     * @brief Describes an event that is used to modify the simulation at the
+     * given time point
+     */
+    struct SystemEvent
+    {
+      /// Type of event determines action performed
+      enum class Type
+      {
+        FAULT_ON,
+        FAULT_OFF
+      };
+
+      /// Time event takes place
+      double      time;
+      /// Event type
+      Type        type;
+      /// ID of element used in event (e.g., bus fault id)
+      std::size_t element_id;
+    };
+
+    /**
+     * @brief Data defined in JSON file for parameterized study
+     */
+    struct StudyData
+    {
+      /// path to system model JSON file
+      fs::path                 system_model_file;
+      /// time step size
+      double                   dt;
+      /// max time
+      double                   tmax;
+      /// set of system events
+      std::vector<SystemEvent> events;
+      /// path to output file
+      fs::path                 output_file;
+      /// path to reference file for validation
+      fs::path                 reference_file;
+      /// Error tolerance (between output file and reference file)
+      double                   error_tol;
+      /// Instance of model data
+      SystemModelData<>        model_data;
+    };
+
+    using json = ::nlohmann::json;
+    using Log  = ::GridKit::Utilities::Logger;
+
+    /**
+     * @brief JSON parser implemntation for `StudyData`
+     */
+    void from_json(const json& j, StudyData& c)
+    {
+      using namespace magic_enum;
+
+      j.at("system_model_file").get_to(c.system_model_file);
+      j.at("dt").get_to(c.dt);
+      j.at("tmax").get_to(c.tmax);
+
+      for (auto& raw_event : j.at("events"))
+      {
+        auto& event = c.events.emplace_back();
+        raw_event.at("time").get_to(event.time);
+        raw_event.at("element_id").get_to(event.element_id);
+
+        auto type_str   = raw_event.at("type").get<std::string>();
+        using EventType = SystemEvent::Type;
+        auto type_wrap  = enum_cast<EventType>(type_str, case_insensitive);
+        if (!type_wrap.has_value())
+        {
+          Log::error() << "Unable to parse event type \"" << type_str << "\"\n";
+        }
+        event.type = type_wrap.value();
+      }
+
+      if (j.contains("output_file"))
+      {
+        j.at("output_file").get_to(c.output_file);
+      }
+
+      if (j.contains("reference_file"))
+      {
+        j.at("reference_file").get_to(c.reference_file);
+      }
+
+      c.error_tol = j.value("error_tolerance", 1.0e-4);
+    }
+
+    /**
+     * @brief Check for existence and successful input file open
+     */
+    std::ifstream openFile(const fs::path& file_path)
+    {
+      if (!exists(file_path))
+      {
+        Log::error() << "File not found: " << file_path << std::endl;
+      }
+      auto fs = std::ifstream(file_path);
+      if (!fs)
+      {
+        Log::error() << "Failed to open file: " << file_path << std::endl;
+      }
+      return fs;
+    }
+
+    /**
+     * @brief Wrapper function to parse `StudyData` from JSON and perform
+     * follow-up configuration
+     */
+    StudyData parseStudyData(const fs::path& file_path)
+    {
+      auto data = StudyData(json::parse(openFile(file_path)));
+
+      auto loc = file_path.parent_path();
+      if (!data.system_model_file.is_absolute())
+      {
+        data.system_model_file = loc / data.system_model_file;
+      }
+      if (!data.reference_file.empty())
+      {
+        if (!data.reference_file.is_absolute())
+        {
+          data.reference_file = loc / data.reference_file;
+        }
+      }
+
+      auto csv        = ::GridKit::Model::VariableMonitorFormat::CSV;
+      data.model_data = parseSystemModelData(data.system_model_file);
+      std::string model_output_file;
+      // Find output file (CSV) specified in model input file
+      for (const auto& sink : data.model_data.monitor_sink)
+      {
+        if (sink.format == csv && sink.delim == ",")
+        {
+          model_output_file = sink.file_name;
+        }
+      }
+
+      if (model_output_file.empty())
+      {
+        // Add study output file to model if one did not already exist
+        data.model_data.monitor_sink.emplace_back(data.output_file, csv);
+      }
+      else
+      {
+        if (data.output_file.empty())
+        {
+          data.output_file = model_output_file;
+        }
+        else
+        {
+          // If model file already specifies a CSV output file, then the study
+          // output file must be a symlink to the model output file
+          if (exists(data.output_file))
+          {
+            if ((!is_symlink(data.output_file)) || (read_symlink(data.output_file) != model_output_file))
+            {
+              Log::error() << "Study output file not usable" << std::endl;
+            }
+          }
+          else
+          {
+            fs::create_symlink(model_output_file, data.output_file);
+          }
+        }
+      }
+
+      return data;
+    }
+  } // namespace PhasorDynamics
+} // namespace GridKit

application/PhasorDynamics/README.md

@@ -0,0 +1,25 @@
+# Input file for GridKit phasor dynamics application
+
+## Root elements
+
+   Name               | Value
+ ---------------------|-------------------------------------------------------
+  `system_model_file` | Path to the system model file[^1]
+  `dt`                | A floating point value for time step size
+  `tmax`              | A floating point value for max time
+  `events`            | An array of event groups (see [Events](#events) below)
+  `output_file`       | Path to output (CSV) file
+  `reference_file`    | A string containing the name of the case
+  `error_tolerance`   | A string containing the name of the case
+
+[^1]: See system model [case format](../../Model/PhasorDynamics/INPUT_FORMAT.md)
+
+## Events
+
+Each event group describes a system event that occurs at a given time point
+
+   Name              | Value
+ --------------------|-------------------------------------------------------
+  `time`             | A floating point value for time event occurs
+  `type`             | Event type (one of { "fault_on", "fault_off" })
+  `element_id`       | An integer value referencing the element associated with the event (e.g., bus fault id)