diff --git a/examples/PhasorDynamics/Example1/example1.cpp b/examples/PhasorDynamics/Example1/example1.cpp index 40bacf64..300bf0d4 100644 --- a/examples/PhasorDynamics/Example1/example1.cpp +++ b/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 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(); - /* Set up simulation */ + output.push_back(OutputData{t, yval[0], yval[1], yval[3]}); + }; + + // 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); + + // Run for 1s ida.initializeSimulation(0.0, false); - ida.runSimulationFixed(0.0, dt, 1.0, buffer); + 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.runSimulationFixed(1.0, dt, 1.1, buffer); + 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.runSimulationFixed(1.1, dt, 10.0, buffer); + nout = static_cast(std::round((10.0 - 1.1) / dt)); + ida.runSimulation(10.0, nout, output_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) + // 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& 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; diff --git a/examples/PhasorDynamics/Example2/example2.cpp b/examples/PhasorDynamics/Example2/example2.cpp index f2a72370..c7939ead 100644 --- a/examples/PhasorDynamics/Example2/example2.cpp +++ b/examples/PhasorDynamics/Example2/example2.cpp @@ -25,13 +25,58 @@ #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) + { + 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(1e-4); @@ -65,33 +110,48 @@ 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(); - /* 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 ida(&sys); ida.configureSimulation(); - /* Run simulation */ + // Run simulation, output each `dt` interval scalar_type start = static_cast(clock()); ida.initializeSimulation(0.0, false); - ida.runSimulationFixed(0.0, dt, 1.0, buffer); + + // 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.runSimulationFixed(1.0, dt, 1.1, buffer); + 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.runSimulationFixed(1.1, dt, 10.0, buffer); + 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 */ 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(); diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 8a4b07ac..7cce626d 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) { @@ -678,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 1fcbbca0..e6d07f57 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(); @@ -104,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, 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..4509d17c --- /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..53c30c68 --- /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([[maybe_unused]] real_type& rel_tol, [[maybe_unused]] real_type& abs_tol) const override + { + } + + void setMaxSteps(IdxT& msa) const override + { + msa = 2000; + } + + int tagDifferentiable() override + { + return 0; + } + + int evaluateResidual() override + { + f_ = y_; + 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([[maybe_unused]] real_type t, [[maybe_unused]] 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 = [&]([[maybe_unused]] 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 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(); +}