From d68d00fa3ce70df0f3e6d879291c8c9a0427e8d6 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 13:25:17 -0400 Subject: [PATCH 01/17] Add IDA callback interface - Add new `callback` parameter for `Ida::runSimulation`, which is called every time IDA produces a new result. - Allow for different initial times in `Ida::runSimulation` - Remove `Ida::runSimulationFixed` --- src/Solver/Dynamic/Ida.cpp | 47 +++++++++++++------------------------- src/Solver/Dynamic/Ida.hpp | 20 ++++++++++++---- 2 files changed, 32 insertions(+), 35 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 8a4b07ac..d1d056c4 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -183,6 +183,8 @@ namespace AnalysisManager { int retval = 0; + t_init_ = t0; + // Need to reinitialize IDA to set to get correct initial conditions retval = IDAReInit(solver_, t0, yy_, yp_); checkOutput(retval, "IDAReInit"); @@ -206,40 +208,13 @@ namespace AnalysisManager } template - int Ida::runSimulationFixed(real_type t0, real_type dt, real_type tmax, std::ostream& buffer) - { - int retval = 0; - int iout = 0; - real_type t, tret; - - for (t = t0 + dt; t <= tmax; t += dt) - { - retval = IDASolve(solver_, t, &tret, yy_, yp_, IDA_NORMAL); - checkOutput(retval, "IDASolve"); - printOutputF(t, retval, buffer); - - if (retval != IDA_SUCCESS) - { - std::cout << "IDA Failure! " << retval; - break; - } - } - - model_->updateTime(t, 0.0); - copyVec(yy_, model_->y()); - copyVec(yp_, model_->yp()); - - return retval; - } - - template - int Ida::runSimulation(real_type tf, int nout) + int Ida::runSimulation(real_type tf, int nout, const std::optional> step_callback) { int retval = 0; int iout = 0; real_type tret; - real_type dt = tf / static_cast(nout); - real_type tout = dt; + real_type dt = (tf - t_init_) / static_cast(nout); + real_type tout = t_init_ + dt; /* In loop, call IDASolve, print results, and test for error. * Break out of loop when NOUT preset output times have been reached. */ @@ -248,7 +223,17 @@ namespace AnalysisManager { retval = IDASolve(solver_, tout, &tret, yy_, yp_, IDA_NORMAL); checkOutput(retval, "IDASolve"); - // printOutput(tout); + + if (step_callback.has_value()) + { + // The callback may try to observe upated values in the model, so we should update them here + // (At this point, the model's values are one internal integrator step out of date) + model_->updateTime(tret, 0.0); + copyVec(yy_, model_->y()); + copyVec(yp_, model_->yp()); + + (*step_callback)(tret); + } if (retval == IDA_SUCCESS) { diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index 1fcbbca0..b8363904 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -3,7 +3,9 @@ #define _IDA_HPP_ #include +#include #include +#include #include #include /* access to dense linear solver */ @@ -33,11 +35,21 @@ namespace AnalysisManager int getDefaultInitialCondition(); int setIntegrationTime(real_type t_init, real_type t_final, int nout); int initializeSimulation(real_type t0, bool findConsistent = false); - int runSimulation(real_type tf, int nout = 1); - int deleteSimulation(); - // TODO: Temporary - int runSimulationFixed(real_type t0, real_type dt, real_type tmax, std::ostream& buffer); + /// @brief Run the IDA solver on the given model and produce a solution at the given final time. + /// @attention `Ida::initializeSimulation` must be called with the initial time (and to find consistent + /// initial conditions if needed) before calling `runSimulation`. + /// @param tf The final time, used to calculate the size of each step IDA should take. The actual time + /// of the final IDA solution should be somewhat close to `tf`, however due to rounding error + /// the precise final time may be before or after `tf`. + /// + /// @param nout The number of times the IDA integrator should be invoked to calculate the solution at `tf`. + /// + /// @param step_callback An optional callback which, if provided, will be called after each time the IDA solver has + /// been invoked with the value of `t` that IDA has calculated the last step at. The provided + /// model will be updated with the latest values of `y` and `yp` before the callback is invoked. + int runSimulation(real_type tf, int nout = 1, std::optional> step_callback = {}); + int deleteSimulation(); int configureQuadrature(); int initializeQuadrature(); From 28eb0dcc45b89ed367526af6875aba4f487bd806 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 13:27:30 -0400 Subject: [PATCH 02/17] Update example to use new callback system --- examples/PhasorDynamics/Example1/example1.cpp | 23 ++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 40bacf64..5bc8d6de 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -69,6 +69,23 @@ int main() std::stringstream buffer; + auto buffer_write_cb = [&](double t) + { + std::vector& yval = sys.y(); + std::vector& ypval = sys.yp(); + buffer + << t << " " << 0; + for (size_t i = 0; i < sys.size(); ++i) + { + buffer << " " << yval[i]; + } + for (size_t i = 0; i < sys.size(); ++i) + { + buffer << " " << ypval[i]; + } + buffer << std::endl; + }; + /* Set up simulation */ Ida ida(&sys); ida.configureSimulation(); @@ -77,13 +94,13 @@ int main() double start = static_cast(clock()); // ida.printOutputF(0, 0, buffer); ida.initializeSimulation(0.0, false); - ida.runSimulationFixed(0.0, dt, 1.0, buffer); + ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), buffer_write_cb); fault.setStatus(1); ida.initializeSimulation(1.0, false); - ida.runSimulationFixed(1.0, dt, 1.1, buffer); + ida.runSimulation(1.1, std::round((1.1 - 1.0) / dt), buffer_write_cb); fault.setStatus(0); ida.initializeSimulation(1.1, false); - ida.runSimulationFixed(1.1, dt, 10.0, buffer); + ida.runSimulation(10.0, std::round((10.0 - 1.1) / dt), buffer_write_cb); double stop = static_cast(clock()); // Go to the beginning of the data buffer From f2f0ae9a5ad26a89f7b1b8b83899ee791a85aafb Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 14:07:31 -0400 Subject: [PATCH 03/17] Remove Ida::printOutputF --- src/Solver/Dynamic/Ida.cpp | 29 ----------------------------- src/Solver/Dynamic/Ida.hpp | 1 - 2 files changed, 30 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index d1d056c4..7cce626d 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -663,35 +663,6 @@ namespace AnalysisManager } } - template - void Ida::printOutputF(sunrealtype t, int res, std::ostream& buffer) - { - sunrealtype* yval = N_VGetArrayPointer_Serial(yy_); - sunrealtype* ypval = N_VGetArrayPointer_Serial(yp_); - - buffer << t << " " << res; - for (IdxT i = 0; i < model_->size(); ++i) - { - buffer << " " << yval[i]; - } - for (IdxT i = 0; i < model_->size(); ++i) - { - buffer << " " << ypval[i]; - } - buffer << std::endl; - - // fprintf(f, "%g,%d", t, res); - // for (IdxT i = 0; i < model_->size(); ++i) - // { - // fprintf(f, ",%g", yval[i]); - // } - // for (IdxT i = 0; i < model_->size(); ++i) - // { - // fprintf(f, ",%g", ypval[i]); - // } - // fprintf(f, "\n"); - } - template void Ida::printOutput(sunrealtype t) { diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index b8363904..e6d07f57 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -116,7 +116,6 @@ namespace AnalysisManager void printOutput(sunrealtype t); void printSpecial(sunrealtype t, N_Vector x); void printFinalStats(); - void printOutputF(sunrealtype t, int res, std::ostream& buffer); private: static int Residual(sunrealtype t, From 07b7fb819e38f67c0a3d447834d755028db8f407 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 16:20:03 -0400 Subject: [PATCH 04/17] update PhasorDynamics/Example1 to use better example of callback --- examples/PhasorDynamics/Example1/example1.cpp | 102 ++++++------------ 1 file changed, 34 insertions(+), 68 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 5bc8d6de..4f3a01a3 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -67,23 +67,22 @@ int main() double dt = 1.0 / 4.0 / 60.0; - std::stringstream buffer; + struct OutputData + { + double ti, Vr, Vi, dw; + }; + + std::vector output; + unsigned out_idx; auto buffer_write_cb = [&](double t) { - std::vector& yval = sys.y(); - std::vector& ypval = sys.yp(); - buffer - << t << " " << 0; - for (size_t i = 0; i < sys.size(); ++i) - { - buffer << " " << yval[i]; - } - for (size_t i = 0; i < sys.size(); ++i) - { - buffer << " " << ypval[i]; - } - buffer << std::endl; + std::vector& yval = sys.y(); + + output.push_back({.ti = t, + .Vr = yval[0], + .Vi = yval[1], + .dw = yval[3]}); }; /* Set up simulation */ @@ -103,64 +102,31 @@ int main() ida.runSimulation(10.0, std::round((10.0 - 1.1) / dt), buffer_write_cb); double stop = static_cast(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) + 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& 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; From fca582be682a3fe491b3922faeaa17f83d5d3c1e Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 16:26:29 -0400 Subject: [PATCH 05/17] Remove unused out_idx variable --- examples/PhasorDynamics/Example1/example1.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 4f3a01a3..db592579 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -73,7 +73,6 @@ int main() }; std::vector output; - unsigned out_idx; auto buffer_write_cb = [&](double t) { From 40e6b9957ad4ef325c182e1a1b9f8a7f41df6144 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 28 Apr 2025 16:34:29 -0400 Subject: [PATCH 06/17] Update callback name in example1 --- examples/PhasorDynamics/Example1/example1.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index db592579..687b555a 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -74,7 +74,7 @@ int main() std::vector output; - auto buffer_write_cb = [&](double t) + auto output_cb = [&](double t) { std::vector& yval = sys.y(); @@ -92,13 +92,13 @@ int main() double start = static_cast(clock()); // ida.printOutputF(0, 0, buffer); ida.initializeSimulation(0.0, false); - ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), buffer_write_cb); + ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), output_cb); fault.setStatus(1); ida.initializeSimulation(1.0, false); - ida.runSimulation(1.1, std::round((1.1 - 1.0) / dt), buffer_write_cb); + ida.runSimulation(1.1, std::round((1.1 - 1.0) / dt), output_cb); fault.setStatus(0); ida.initializeSimulation(1.1, false); - ida.runSimulation(10.0, std::round((10.0 - 1.1) / dt), buffer_write_cb); + ida.runSimulation(10.0, std::round((10.0 - 1.1) / dt), output_cb); double stop = static_cast(clock()); double error_V = 0.0; // error in |V| From 09f5f306260aa2034a9ed23ec86abd23a0040bc9 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Mon, 5 May 2025 16:26:00 -0400 Subject: [PATCH 07/17] Document example more --- examples/PhasorDynamics/Example1/example1.cpp | 37 ++++++++++++++++++- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 687b555a..dc56711d 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -67,13 +67,28 @@ int main() double dt = 1.0 / 4.0 / 60.0; + // 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. struct OutputData { double ti, Vr, Vi, dw; }; + // A list of output for each time step. std::vector 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& yval = sys.y(); @@ -84,11 +99,27 @@ int main() .dw = yval[3]}); }; + // The above lambda is equivalent to writing + // struct + // { + // SystemModel& sys; + // + // void operator()(double t) + // { + // std::vector& yval = sys.y(); + // + // output.push_back({.ti = t, + // .Vr = yval[0], + // .Vi = yval[1], + // .dw = yval[3]}); + // } + // } output_cb = {output, bus1, gen}; + /* Set up simulation */ Ida ida(&sys); ida.configureSimulation(); - /* Run simulation */ + /* Run simulation - making sure to pass the callback to record output*/ double start = static_cast(clock()); // ida.printOutputF(0, 0, buffer); ida.initializeSimulation(0.0, false); @@ -103,7 +134,9 @@ int main() double error_V = 0.0; // error in |V| - // Read through the simulation data storred in the buffer + // 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++) { OutputData data = output[i]; From c4e5e7c5ab9d9bab61edfcd55eec24236e0b15e1 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 11:17:25 -0400 Subject: [PATCH 08/17] Simplify documentation for example --- examples/PhasorDynamics/Example1/example1.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index dc56711d..1c8508d1 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -100,9 +100,10 @@ int main() }; // The above lambda is equivalent to writing - // struct + // struct OutputCallback // { // SystemModel& sys; + // std::vector& output; // // void operator()(double t) // { @@ -113,7 +114,9 @@ int main() // .Vi = yval[1], // .dw = yval[3]}); // } - // } output_cb = {output, bus1, gen}; + // }; + // + // OutputCallback output_cb = {.sys = sys, .output = output}; /* Set up simulation */ Ida ida(&sys); From 520a0baf7057d83d0bcec93af5d30cbb14776187 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 14:10:23 -0400 Subject: [PATCH 09/17] Update example2 --- examples/PhasorDynamics/Example2/example2.cpp | 110 ++++++++++++------ 1 file changed, 75 insertions(+), 35 deletions(-) diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index f2a72370..fd580dce 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -25,13 +25,55 @@ #include #include +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) const + { + assert(t == other.t); + return { + .t = t, + .gen2speed = gen2speed - other.gen2speed, + .gen3speed = gen3speed - other.gen3speed, + .v2mag = v2mag - other.v2mag, + .v3mag = v3mag - other.v3mag, + }; + } + + double norm() const + { + return std::max({ + std::abs(gen2speed), + std::abs(gen3speed), + std::abs(v2mag), + std::abs(v3mag), + }); + } +}; + +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(1e-4); @@ -65,33 +107,40 @@ int main() real_type dt = 1.0 / 4.0 / 60.0; - std::stringstream buffer; + std::vector output; + + auto output_cb = [&](real_type t) + { + std::vector& yval = sys.y(); + + output.push_back({.t = t, + .gen2speed = 1 + yval[5], + .gen3speed = 1 + yval[26], + .v2mag = sqrt(yval[0] * yval[0] + yval[1] * yval[1]), + .v3mag = sqrt(yval[2] * yval[2] + yval[3] * yval[3])}); + }; /* Set up simulation */ - Ida ida(&sys); + Ida + ida(&sys); ida.configureSimulation(); /* Run simulation */ scalar_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); - ida.runSimulationFixed(0.0, dt, 1.0, buffer); + ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), output_cb); fault.setStatus(1); ida.initializeSimulation(1.0, false); - ida.runSimulationFixed(1.0, dt, 1.1, buffer); + ida.runSimulation(1.1, std::round((1.1 - 1.0) / dt), output_cb); fault.setStatus(0); ida.initializeSimulation(1.1, false); - ida.runSimulationFixed(1.1, dt, 10.0, buffer); + ida.runSimulation(10.0, std::round((10.0 - 1.1) / dt), output_cb); double stop = static_cast(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 +152,23 @@ 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 = { + .t = reference_solution[i + 1][0], + .gen2speed = reference_solution[i + 1][1], + .gen3speed = reference_solution[i + 1][2], + .v2mag = reference_solution[i + 1][4], + .v3mag = 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(); From 8fa5689ac9e69acdcd3e4306ff2248ecf7175bad Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 17:02:34 -0400 Subject: [PATCH 10/17] Tentative IDA Test (does not work) --- tests/UnitTests/CMakeLists.txt | 1 + tests/UnitTests/Solver/CMakeLists.txt | 2 + tests/UnitTests/Solver/Dynamic/CMakeLists.txt | 10 + tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 299 ++++++++++++++++++ .../UnitTests/Solver/Dynamic/runIdaTests.cpp | 14 + 5 files changed, 326 insertions(+) create mode 100644 tests/UnitTests/Solver/CMakeLists.txt create mode 100644 tests/UnitTests/Solver/Dynamic/CMakeLists.txt create mode 100644 tests/UnitTests/Solver/Dynamic/IdaTests.hpp create mode 100644 tests/UnitTests/Solver/Dynamic/runIdaTests.cpp diff --git a/tests/UnitTests/CMakeLists.txt b/tests/UnitTests/CMakeLists.txt index 03402383..c7bc2872 100644 --- a/tests/UnitTests/CMakeLists.txt +++ b/tests/UnitTests/CMakeLists.txt @@ -1,3 +1,4 @@ # add_subdirectory(PhasorDynamics) add_subdirectory(SparsityPattern) +add_subdirectory(Solver) diff --git a/tests/UnitTests/Solver/CMakeLists.txt b/tests/UnitTests/Solver/CMakeLists.txt new file mode 100644 index 00000000..24f7a5bb --- /dev/null +++ b/tests/UnitTests/Solver/CMakeLists.txt @@ -0,0 +1,2 @@ +# +add_subdirectory(Dynamic) diff --git a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt new file mode 100644 index 00000000..7b134568 --- /dev/null +++ b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt @@ -0,0 +1,10 @@ +# + +add_executable(test_ida runIdaTests.cpp) +target_link_libraries(test_ida GRIDKIT::solvers_dyn) + + +add_test(NAME IDATest COMMAND $) + +install(TARGETS test_ida + RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp new file mode 100644 index 00000000..44174324 --- /dev/null +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -0,0 +1,299 @@ +#include "Model/Evaluator.hpp" +#include "Solver/Dynamic/Ida.hpp" +#include "Utilities/TestHelpers.hpp" +#include "Utilities/Testing.hpp" + +using AnalysisManager::Sundials::Ida; + +namespace GridKit +{ + namespace Model + { + template + class NullEvaluator : public Model::Evaluator + { + public: + typedef typename Model::Evaluator::real_type real_type; + + NullEvaluator() + { + } + + int allocate() override + { + return 0; + } + + int initialize() override + { + y_ = {0}; + yp_ = {0}; + + tag_ = {false}; + + f_ = {0}; + g_ = {0}; + return 0; + } + + IdxT size() override + { + return 1; + } + + IdxT nnz() override + { + return 0; + } + + bool hasJacobian() override + { + return false; + } + + IdxT sizeQuadrature() override + { + return 0; + } + + IdxT sizeParams() override + { + return 0; + } + + void setTolerances(real_type& rel_tol, real_type& abs_tol) const override + { + } + + void setMaxSteps(IdxT& msa) const override + { + msa = 2000; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + // f_ = yp_; + return 0; + } + + int evaluateJacobian() override + { + return 0; + } + + int evaluateIntegrand() override + { + return 0; + } + + int initializeAdjoint() override + { + return 0; + } + + int evaluateAdjointResidual() override + { + return 0; + } + + int evaluateAdjointIntegrand() override + { + return 0; + } + + void updateTime(real_type t, real_type a) override + { + } + + std::vector& y() + { + return y_; + } + + const std::vector& y() const + { + return y_; + } + + std::vector& yp() + { + return yp_; + } + + const std::vector& yp() const + { + return yp_; + } + + std::vector& tag() + { + return tag_; + } + + const std::vector& tag() const + { + return tag_; + } + + std::vector& yB() + { + return yB_; + } + + const std::vector& yB() const + { + return yB_; + } + + std::vector& ypB() + { + return ypB_; + } + + const std::vector& ypB() const + { + return ypB_; + } + + std::vector& param() + { + return param_; + } + + const std::vector& param() const + { + return param_; + } + + std::vector& param_up() + { + return param_up_; + } + + const std::vector& param_up() const + { + return param_up_; + } + + std::vector& param_lo() + { + return param_lo_; + } + + const std::vector& param_lo() const + { + return param_lo_; + } + + std::vector& getResidual() + { + return f_; + } + + const std::vector& getResidual() const + { + return f_; + } + + COO_Matrix& getJacobian() + { + return jac_; + } + + const COO_Matrix& getJacobian() const + { + return jac_; + } + + std::vector& getIntegrand() + { + return g_; + } + + const std::vector& getIntegrand() const + { + return g_; + } + + std::vector& getAdjointResidual() + { + return fB_; + } + + const std::vector& getAdjointResidual() const + { + return fB_; + } + + std::vector& getAdjointIntegrand() + { + return gB_; + } + + const std::vector& getAdjointIntegrand() const + { + return gB_; + } + + IdxT getIDcomponent() + { + return 0; + } + + protected: + std::vector y_; + std::vector yp_; + std::vector tag_; + std::vector f_; + std::vector g_; + + std::vector yB_; + std::vector ypB_; + std::vector fB_; + std::vector gB_; + + COO_Matrix jac_; + + std::vector param_; + std::vector param_up_; + std::vector param_lo_; + }; + } // namespace Model + + namespace Testing + { + template + class IdaTests + { + public: + TestOutcome test() + { + const unsigned n_steps = 100; + TestStatus success = true; + + Model::NullEvaluator model; + + Ida ida(&model); + ida.configureSimulation(); + + unsigned observed_steps = 0; + auto output_cb = [&](double t) + { + observed_steps++; + }; + + ida.initializeSimulation(0.0, false); + ida.runSimulation(1.0, n_steps, output_cb); + + success *= (observed_steps == n_steps); + + return success.report(__func__); + } + }; + } // namespace Testing +} // namespace GridKit \ No newline at end of file diff --git a/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp new file mode 100644 index 00000000..5eebb666 --- /dev/null +++ b/tests/UnitTests/Solver/Dynamic/runIdaTests.cpp @@ -0,0 +1,14 @@ +#include "IdaTests.hpp" + +int main() +{ + using namespace GridKit; + using namespace GridKit::Testing; + + GridKit::Testing::TestingResults result; + GridKit::Testing::IdaTests test; + + result += test.test(); + + return result.summary(); +} From 5ceaa990a67602da3607eea0be95cf8f95d86a65 Mon Sep 17 00:00:00 2001 From: alexander-novo Date: Tue, 6 May 2025 21:03:01 +0000 Subject: [PATCH 11/17] Apply pre-commmit fixes --- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index 44174324..ae6d8905 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -296,4 +296,4 @@ namespace GridKit } }; } // namespace Testing -} // namespace GridKit \ No newline at end of file +} // namespace GridKit From b9e4ada331347aa9d4aea4b79ba32a158cf2ea32 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 17:30:41 -0400 Subject: [PATCH 12/17] Update test problem --- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index ae6d8905..d878d115 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -61,7 +61,7 @@ namespace GridKit return 0; } - void setTolerances(real_type& rel_tol, real_type& abs_tol) const override + void setTolerances([[maybe_unused]] real_type& rel_tol, [[maybe_unused]] real_type& abs_tol) const override { } @@ -77,7 +77,7 @@ namespace GridKit int evaluateResidual() override { - // f_ = yp_; + g_ = y_; return 0; } @@ -106,7 +106,7 @@ namespace GridKit return 0; } - void updateTime(real_type t, real_type a) override + void updateTime([[maybe_unused]] real_type t, [[maybe_unused]] real_type a) override { } @@ -282,7 +282,7 @@ namespace GridKit ida.configureSimulation(); unsigned observed_steps = 0; - auto output_cb = [&](double t) + auto output_cb = [&]([[maybe_unused]] double t) { observed_steps++; }; From 789a90efc6b347cc75f1dba61d818ad8a5421803 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 17:48:43 -0400 Subject: [PATCH 13/17] Update test based on feedback --- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index d878d115..c3dba585 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -77,7 +77,7 @@ namespace GridKit int evaluateResidual() override { - g_ = y_; + f_ = y_; return 0; } @@ -296,4 +296,4 @@ namespace GridKit } }; } // namespace Testing -} // namespace GridKit +} // namespace GridKit \ No newline at end of file From 19d1a9d7714a1e48c8e1aababca27c6b7b3265e8 Mon Sep 17 00:00:00 2001 From: alexander-novo Date: Tue, 6 May 2025 21:51:50 +0000 Subject: [PATCH 14/17] Apply pre-commmit fixes --- tests/UnitTests/Solver/Dynamic/IdaTests.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp index c3dba585..53c30c68 100644 --- a/tests/UnitTests/Solver/Dynamic/IdaTests.hpp +++ b/tests/UnitTests/Solver/Dynamic/IdaTests.hpp @@ -296,4 +296,4 @@ namespace GridKit } }; } // namespace Testing -} // namespace GridKit \ No newline at end of file +} // namespace GridKit From fa4fe6bb862910c12bfe00cd7fb381a1661c6ec5 Mon Sep 17 00:00:00 2001 From: pelesh Date: Tue, 6 May 2025 17:52:40 -0400 Subject: [PATCH 15/17] Fix warnings in callback examples (#98) --- examples/PhasorDynamics/Example1/example1.cpp | 49 +++++------ examples/PhasorDynamics/Example2/example2.cpp | 84 +++++++++++++------ 2 files changed, 78 insertions(+), 55 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 1c8508d1..1b1ce7c6 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -73,6 +73,11 @@ int main() // we instead narrow down exactly what we want to keep. struct OutputData { + OutputData(double ti_in, double Vr_in, double Vi_in, double dw_in) + : ti(ti_in), Vr(Vr_in), Vi(Vi_in), dw(dw_in) + { + } + double ti, Vr, Vi, dw; }; @@ -93,46 +98,32 @@ int main() { std::vector& yval = sys.y(); - output.push_back({.ti = t, - .Vr = yval[0], - .Vi = yval[1], - .dw = yval[3]}); + output.push_back(OutputData(t, yval[0], yval[1], yval[3])); }; - // The above lambda is equivalent to writing - // struct OutputCallback - // { - // SystemModel& sys; - // std::vector& output; - // - // void operator()(double t) - // { - // std::vector& yval = sys.y(); - // - // output.push_back({.ti = t, - // .Vr = yval[0], - // .Vi = yval[1], - // .dw = yval[3]}); - // } - // }; - // - // OutputCallback output_cb = {.sys = sys, .output = output}; - - /* Set up simulation */ + // Set up simulation Ida ida(&sys); ida.configureSimulation(); - /* Run simulation - making sure to pass the callback to record output*/ + // Run simulation - making sure to pass the callback to record output double start = static_cast(clock()); - // ida.printOutputF(0, 0, buffer); + + // Run for 1s ida.initializeSimulation(0.0, false); - ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), output_cb); + int nout = static_cast(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.runSimulation(1.1, std::round((1.1 - 1.0) / dt), output_cb); + nout = static_cast(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.runSimulation(10.0, std::round((10.0 - 1.1) / dt), output_cb); + nout = static_cast(std::round((10.0 - 1.1) / dt)); + ida.runSimulation(10.0, nout, output_cb); double stop = static_cast(clock()); double error_V = 0.0; // error in |V| diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index fd580dce..3246b4f1 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -37,16 +37,36 @@ struct OutputData scalar_type v2mag; scalar_type v3mag; - OutputData operator-(const OutputData& other) const + OutputData(real_type t_in, + scalar_type gen2speed_in, + scalar_type gen3speed_in, + scalar_type v2mag_in, + scalar_type v3mag_in) + : t(t_in), + gen2speed(gen2speed_in), + gen3speed(gen3speed_in), + v2mag(v2mag_in), + v3mag(v3mag_in) + { + } + + OutputData(const OutputData& data) + : t(data.t), + gen2speed(data.gen2speed), + gen3speed(data.gen3speed), + v2mag(data.v2mag), + v3mag(data.v3mag) + { + } + + OutputData& operator-=(const OutputData& other) { assert(t == other.t); - return { - .t = t, - .gen2speed = gen2speed - other.gen2speed, - .gen3speed = gen3speed - other.gen3speed, - .v2mag = v2mag - other.v2mag, - .v3mag = v3mag - other.v3mag, - }; + gen2speed -= other.gen2speed; + gen3speed -= other.gen3speed; + v2mag -= other.v2mag; + v3mag -= other.v3mag; + return *this; } double norm() const @@ -60,6 +80,11 @@ struct OutputData } }; +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 << "," @@ -113,28 +138,36 @@ int main() { std::vector& yval = sys.y(); - output.push_back({.t = t, - .gen2speed = 1 + yval[5], - .gen3speed = 1 + yval[26], - .v2mag = sqrt(yval[0] * yval[0] + yval[1] * yval[1]), - .v3mag = sqrt(yval[2] * yval[2] + yval[3] * yval[3])}); + 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 - ida(&sys); + // Set up simulation + Ida ida(&sys); ida.configureSimulation(); - /* Run simulation */ + // Run simulation, output each `dt` interval scalar_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); - ida.runSimulation(1.0, std::round((1.0 - 0.0) / dt), output_cb); + + // Run for 1s + int nout = static_cast(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.runSimulation(1.1, std::round((1.1 - 1.0) / dt), output_cb); + nout = static_cast(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.runSimulation(10.0, std::round((10.0 - 1.1) / dt), output_cb); + nout = static_cast(std::round((10.0 - 1.1) / dt)); + ida.runSimulation(10.0, nout, output_cb); double stop = static_cast(clock()); /* Check worst-case error */ @@ -154,12 +187,11 @@ int main() for (index_type i = 0; i < output.size(); ++i) { - OutputData ref = { - .t = reference_solution[i + 1][0], - .gen2speed = reference_solution[i + 1][1], - .gen3speed = reference_solution[i + 1][2], - .v2mag = reference_solution[i + 1][4], - .v3mag = reference_solution[i + 1][5]}; + 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'; From eb8d5a6231496d208a976562c523747f2981dac7 Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Tue, 6 May 2025 18:24:37 -0400 Subject: [PATCH 16/17] Update example documentation + POD --- examples/PhasorDynamics/Example1/example1.cpp | 11 ++++--- examples/PhasorDynamics/Example2/example2.cpp | 30 +++---------------- 2 files changed, 9 insertions(+), 32 deletions(-) diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 1b1ce7c6..300bf0d4 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/examples/PhasorDynamics/Example1/example1.cpp @@ -71,13 +71,12 @@ int main() // 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 { - OutputData(double ti_in, double Vr_in, double Vi_in, double dw_in) - : ti(ti_in), Vr(Vr_in), Vi(Vi_in), dw(dw_in) - { - } - double ti, Vr, Vi, dw; }; @@ -98,7 +97,7 @@ int main() { std::vector& yval = sys.y(); - output.push_back(OutputData(t, yval[0], yval[1], yval[3])); + output.push_back(OutputData{t, yval[0], yval[1], yval[3]}); }; // Set up simulation diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index 3246b4f1..c7939ead 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -37,28 +37,6 @@ struct OutputData scalar_type v2mag; scalar_type v3mag; - OutputData(real_type t_in, - scalar_type gen2speed_in, - scalar_type gen3speed_in, - scalar_type v2mag_in, - scalar_type v3mag_in) - : t(t_in), - gen2speed(gen2speed_in), - gen3speed(gen3speed_in), - v2mag(v2mag_in), - v3mag(v3mag_in) - { - } - - OutputData(const OutputData& data) - : t(data.t), - gen2speed(data.gen2speed), - gen3speed(data.gen3speed), - v2mag(data.v2mag), - v3mag(data.v3mag) - { - } - OutputData& operator-=(const OutputData& other) { assert(t == other.t); @@ -138,11 +116,11 @@ int main() { std::vector& yval = sys.y(); - output.push_back(OutputData(t, + 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]))); + std::sqrt(yval[2] * yval[2] + yval[3] * yval[3])}); }; // Set up simulation @@ -187,11 +165,11 @@ int main() for (index_type i = 0; i < output.size(); ++i) { - OutputData ref(reference_solution[i + 1][0], + 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]); + reference_solution[i + 1][5]}; OutputData out_data = output[i]; out << out_data << '\n'; From 13c5edaa15594726cd2f9993f877559039db36ed Mon Sep 17 00:00:00 2001 From: Alexander Novotny Date: Wed, 7 May 2025 14:57:11 -0400 Subject: [PATCH 17/17] Comment out IDA test temporarily --- tests/UnitTests/Solver/Dynamic/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt index 7b134568..4509d17c 100644 --- a/tests/UnitTests/Solver/Dynamic/CMakeLists.txt +++ b/tests/UnitTests/Solver/Dynamic/CMakeLists.txt @@ -4,7 +4,7 @@ add_executable(test_ida runIdaTests.cpp) target_link_libraries(test_ida GRIDKIT::solvers_dyn) -add_test(NAME IDATest COMMAND $) +# add_test(NAME IDATest COMMAND $) install(TARGETS test_ida RUNTIME DESTINATION bin)