integrate co-simulation into phasor dynamics applications

- add client/server modes to `DynamicSimulation` and `ContingencyAnalysis`
- introduce simulation runners to coordinate events, interface solves, and
  monitor output
- remove standalone co-simulation examples
- add a partitioned ieee39 example
- document co-simulation configuration and execution

Author: superwhiskers

GridKit/Utilities/Logger/Logger.cpp

@@ -9,6 +9,9 @@
 #include <GridKit/Definitions.hpp>
 #include <GridKit/Utilities/Colors.hpp>
 
+// TODO: allow instantiating a scoped logger that prefixes with some
+//       class/logical distinction name like `[CliArgs]` or `[GENROU]`.
+
 namespace GridKit
 {
   namespace Utilities
@@ -129,14 +132,16 @@ namespace GridKit
     std::ostream& Logger::warning()
     {
       using namespace Colors;
-      *(output_streams_[WARNINGS]) << "[" << YELLOW << "WARNING" << CLEAR << "] ";
+      *(output_streams_[WARNINGS])
+          << "[" << YELLOW << "WARNING" << CLEAR << "] ";
       return *(output_streams_[WARNINGS]);
     }
 
     /**
      * @brief Returns reference to analysis summary messages output stream.
      *
-     * @return Reference to analysis summary messages stream in `output_streams_`.
+     * @return Reference to analysis summary messages stream in
+     * `output_streams_`.
      *
      * @pre `output_streams_` vector is allocated and correctly initialized.
      */
@@ -177,7 +182,8 @@ namespace GridKit
     }
 
     /**
-     * @brief Set outputs of active streams to user provided `std::ostream` object.
+     * @brief Set outputs of active streams to user provided `std::ostream`
+     *        object.
      *
      * All active outputs are redirected to `out` stream. All inactive ones are
      * directed to null device.

GridKit/Utilities/Logger/Logger.hpp

@@ -16,7 +16,6 @@ namespace GridKit
      * @brief Manages and logs outputs from GridKit code.
      *
      * All methods and data in this class are static.
-     *
      */
     class Logger
     {

application/PhasorDynamics/AnalysisUtilities.hpp

@@ -1,8 +1,8 @@
 #pragma once
 
+#include <cmath>
 #include <cstddef>
 #include <filesystem>
-#include <format>
 #include <fstream>
 #include <iostream>
 #include <stdexcept>
@@ -30,7 +30,8 @@ namespace GridKit
      * @brief Describes an event that is used to modify the simulation at the
      * given time point
      *
-     * ContingencyAnalysis repeats fault events at each bus with the same R and X.
+     * ContingencyAnalysis repeats fault events at each bus with the same R and
+     * X.
      */
     struct SystemEvent
     {
@@ -72,6 +73,12 @@ namespace GridKit
       double                                         rel_tol;
       /// absolute tolerance for the solver
       double                                         abs_tol;
+      /// relative residual tolerance for the co-simulation interface
+      double                                         cosim_rel_tol;
+      /// scaled step tolerance for the co-simulation interface solver
+      double                                         cosim_step_tol;
+      /// estimated relative error of co-simulation residual evaluations
+      double                                         cosim_function_error;
       /// fixed solver time step size, or 0 for adaptive stepping
       double                                         dt_fixed;
       /// maximum number of solver time steps, or 0 for the IDA default
@@ -101,9 +108,12 @@ namespace GridKit
     using json = ::nlohmann::json;
     using Log  = ::GridKit::Utilities::Logger;
 
-    inline constexpr double DEFAULT_SOLVER_REL_TOL   = 1.0e-7;
-    inline constexpr double DEFAULT_SOLVER_ABS_TOL   = 1.0e-9;
-    inline constexpr double DEFAULT_VERIFICATION_TOL = 1.0e-4;
+    inline constexpr double DEFAULT_SOLVER_REL_TOL       = 1.0e-7;
+    inline constexpr double DEFAULT_SOLVER_ABS_TOL       = 1.0e-9;
+    inline constexpr double DEFAULT_VERIFICATION_TOL     = 1.0e-4;
+    inline constexpr double DEFAULT_COSIM_REL_TOL        = 1.0e-8;
+    inline constexpr double DEFAULT_COSIM_STEP_TOL       = 1.0e-12;
+    inline constexpr double DEFAULT_COSIM_FUNCTION_ERROR = 1.0e-10;
 
     /**
      * @brief JSON parser implementation for `StudyData`
@@ -117,20 +127,38 @@ namespace GridKit
       j.at("tmax").get_to(c.tmax);
       c.rel_tol            = j.value("rel_tol", DEFAULT_SOLVER_REL_TOL);
       c.abs_tol            = j.value("abs_tol", DEFAULT_SOLVER_ABS_TOL);
-      c.dt_fixed           = j.value("dt_fixed", 0.0);
-      c.max_steps          = j.value("max_steps", std::size_t{0});
-      c.max_order          = j.value("max_order", 5);
-      c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;
+      const auto tolerance = [&](const char* name, double fallback)
+      {
+        const double value = j.value(name, fallback);
+        if (!std::isfinite(value) || value <= 0.0)
+        {
+          throw std::invalid_argument(
+              std::string(name) + " must be positive and finite");
+        }
+        return value;
+      };
+      c.cosim_rel_tol  = tolerance("cosim_rel_tol", DEFAULT_COSIM_REL_TOL);
+      c.cosim_step_tol = tolerance("cosim_step_tol", DEFAULT_COSIM_STEP_TOL);
+      c.cosim_function_error =
+          tolerance("cosim_function_error", DEFAULT_COSIM_FUNCTION_ERROR);
+      c.dt_fixed  = j.value("dt_fixed", 0.0);
+      c.max_steps = j.value("max_steps", std::size_t{0});
+      c.max_order = j.value("max_order", 5);
+      c.consistent_ic_type =
+          AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;
       if (j.contains("consistent_ic_type"))
       {
-        const auto consistent_ic_type_str = j.at("consistent_ic_type").get<std::string>();
+        const auto consistent_ic_type_str =
+            j.at("consistent_ic_type").get<std::string>();
         if (consistent_ic_type_str == "y")
         {
-          c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::Y;
+          c.consistent_ic_type =
+              AnalysisManager::Sundials::IdaConsistentICType::Y;
         }
         else if (consistent_ic_type_str == "ya_ydp")
         {
-          c.consistent_ic_type = AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;
+          c.consistent_ic_type =
+              AnalysisManager::Sundials::IdaConsistentICType::YA_YDP;
         }
         else
         {
@@ -223,7 +251,8 @@ namespace GridKit
         c.error_type = ErrorType::RELATIVE;
       }
 
-      c.abs_err_threshold = j.value("abs_err_threshold", Testing::DEFAULT_ABS_ERROR_THRESHOLD);
+      c.abs_err_threshold =
+          j.value("abs_err_threshold", Testing::DEFAULT_ABS_ERROR_THRESHOLD);
     }
 
     /**
@@ -302,7 +331,8 @@ namespace GridKit
           // 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))
+            if ((!is_symlink(data.output_file))
+                || (read_symlink(data.output_file) != model_output_file))
             {
               Log::error() << "Study output file not usable" << std::endl;
             }
@@ -317,39 +347,23 @@ namespace GridKit
       return data;
     }
 
-    void checkCommandLine(int argc, const std::string& appName)
-    {
-      if (argc < 2)
-      {
-        Log::error() << "No input file provided\n"
-                     << std::format(
-                            "\n"
-                            "Usage:\n"
-                            "       {} <json-input-file>\n"
-                            "\n"
-                            "Please provide a json input file for the study to run.\n"
-                            "\n",
-                            appName);
-        exit(1);
-      }
-    }
-
     Testing::TestStatus checkErrors(
-        const StudyData& study_data,
-        bool             print_results = true)
+        const StudyData& study_data, bool print_results = true)
     {
-      // Generate aggregate errors comparing variable output to reference solution
+      // Generate aggregate errors comparing variable output to reference
+      // solution
       auto func   = std::string{"monitor file vs reference file"};
       auto status = Testing::TestStatus{func.c_str()};
 
       const auto& out_file = study_data.output_file;
       const auto& ref_file = study_data.reference_file;
       if (!out_file.empty() && !ref_file.empty())
       {
-        auto errorSet = Testing::compareCSV(out_file,
-                                            ref_file,
-                                            study_data.error_type,
-                                            study_data.abs_err_threshold);
+        auto errorSet = Testing::compareCSV(
+            out_file,
+            ref_file,
+            study_data.error_type,
+            study_data.abs_err_threshold);
 
         // Print the errors
         if (print_results)
@@ -363,7 +377,8 @@ namespace GridKit
           status *= study_data.error_tol.size() == errorSet->var_errors.size();
           for (std::size_t i = 0; i < study_data.error_tol.size(); ++i)
           {
-            status *= errorSet->var_errors[i].max_value < study_data.error_tol[i];
+            status *=
+                errorSet->var_errors[i].max_value < study_data.error_tol[i];
           }
         }
         else

application/PhasorDynamics/CMakeLists.txt

@@ -20,6 +20,11 @@ target_include_directories(
   ContingencyAnalysis
   PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/nlohmann-json/include)
 
+if(GRIDKIT_ENABLE_ZMQ)
+  target_link_libraries(DynamicSimulation PRIVATE GridKit::co_sim GridKit::solvers_steady)
+  target_link_libraries(ContingencyAnalysis PRIVATE GridKit::co_sim GridKit::solvers_steady)
+endif()
+
 if(GRIDKIT_ENABLE_OPENMP)
   target_link_libraries(ContingencyAnalysis PUBLIC OpenMP::OpenMP_CXX)
 endif()

application/PhasorDynamics/CoSimulation.hpp

@@ -0,0 +1,123 @@
+#pragma once
+
+#include "SimulationOptions.hpp"
+#include "SimulationRunner.hpp"
+
+#ifdef GRIDKIT_ENABLE_ZMQ
+
+#include <optional>
+#include <stdexcept>
+#include <utility>
+
+#include <GridKit/CoSim/CoSim.hpp>
+
+#include "CoSimulationTypes.hpp"
+#include "CoupledClient.hpp"
+#include "CoupledServer.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    inline CoSim::RunResult runCoupledClientStudies(
+        const SimulationOptions& options,
+        NextCoupledStudy         next,
+        CompleteCoupledStudy     complete,
+        bool                     shutdown_peer,
+        bool                     print_results = false)
+    {
+      CoupledClient client(
+          options,
+          std::move(next),
+          std::move(complete),
+          shutdown_peer,
+          print_results);
+
+      CoSim::CoSim<double, CoSim::Disposition::CLIENT> adapter(
+          options.endpoint,
+          CoSim::MessageStepFunction<double>(
+              [&](
+                  const CoSimulationMessage&                    reply,
+                  const CoSim::MessageExchangeFunction<double>& exchange)
+              { return client.handleReply(reply, exchange); }),
+          CoSim::MessageExtractFunction<double>([&](const auto&)
+                                                { return client.extract(); }),
+          options.timeout_ms);
+
+      const auto result = adapter.run();
+
+      if (result == CoSim::RunResult::COMPLETED && !client.finished())
+      {
+        throw std::runtime_error(
+            "co-simulation peer ended before the study completed");
+      }
+
+      return result;
+    }
+
+    inline Testing::TestStatus runCoupledServer(
+        const StudyData& base, const SimulationOptions& options)
+    {
+      CoupledServer server(base, options);
+
+      CoSim::CoSim<double, CoSim::Disposition::SERVER> adapter(
+          options.endpoint,
+          CoSim::MessageStepFunction<double>(
+              [&](const CoSimulationMessage& message, const auto&)
+              { return server.handleRequest(message); }),
+          CoSim::MessageExtractFunction<double>([&](const auto&)
+                                                { return server.extract(); }),
+          options.timeout_ms);
+
+      if (adapter.run() == CoSim::RunResult::STOPPED)
+      {
+        return false;
+      }
+
+      if (!server.finished())
+      {
+        throw std::runtime_error(
+            "co-simulation peer closed before the study completed");
+      }
+
+      return server.status();
+    }
+
+    inline Testing::TestStatus runCoupledStudy(
+        const StudyData& data, const SimulationOptions& options)
+    {
+      if (options.role == "server")
+      {
+        return runCoupledServer(data, options);
+      }
+
+      bool                first  = true;
+      Testing::TestStatus status = true;
+      const auto          next   = [&]() -> std::optional<StudyData>
+      {
+        if (!first)
+        {
+          return std::nullopt;
+        }
+
+        first = false;
+        return data;
+      };
+
+      if (runCoupledClientStudies(
+              options,
+              next,
+              [&](Testing::TestStatus result)
+              { status *= result; },
+              false,
+              true)
+          == CoSim::RunResult::STOPPED)
+      {
+        status = false;
+      }
+
+      return status;
+    }
+  } // namespace PhasorDynamics
+} // namespace GridKit
+#endif

application/PhasorDynamics/CoSimulationTypes.hpp

@@ -0,0 +1,25 @@
+#pragma once
+
+#include <GridKit/Definitions.hpp>
+
+#ifdef GRIDKIT_ENABLE_ZMQ
+
+#include <functional>
+#include <optional>
+
+#include <GridKit/CoSim/Message.hpp>
+
+#include "AnalysisUtilities.hpp"
+
+namespace GridKit
+{
+  namespace PhasorDynamics
+  {
+    /// Alias for the co-simulation message type used in the applications.
+    using CoSimulationMessage = CoSim::MessageData<double>;
+
+    using NextCoupledStudy     = std::function<std::optional<StudyData>()>;
+    using CompleteCoupledStudy = std::function<void(Testing::TestStatus)>;
+  } // namespace PhasorDynamics
+} // namespace GridKit
+#endif