diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index 885134ec..1bea6597 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -14,7 +14,7 @@ #include #include -int main(int argc, char const* argv[]) +int main(int /* argc */, char const** /* argv */) { double abs_tol = 1.0e-8; double rel_tol = 1.0e-8; @@ -22,7 +22,7 @@ int main(int argc, char const* argv[]) // TODO:setup as named parameters // Create circuit model - auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac); + GridKit::PowerElectronicsModel sysmodel(rel_tol, abs_tol, use_jac); size_t idoff = 0; @@ -41,7 +41,7 @@ int main(int argc, char const* argv[]) // internal induct->setExternalConnectionNodes(2, 2); // add component - sysmodel->addComponent(induct); + sysmodel.addComponent(induct); // resistor idoff++; @@ -52,7 +52,7 @@ int main(int argc, char const* argv[]) // output resis->setExternalConnectionNodes(1, 1); // add - sysmodel->addComponent(resis); + sysmodel.addComponent(resis); // voltage source idoff++; @@ -65,54 +65,54 @@ int main(int argc, char const* argv[]) // internal vsource->setExternalConnectionNodes(2, 3); - sysmodel->addComponent(vsource); + sysmodel.addComponent(vsource); - sysmodel->allocate(4); + sysmodel.allocate(4); - std::cout << sysmodel->y().size() << std::endl; + std::cout << sysmodel.y().size() << std::endl; // Grounding for IDA. If no grounding then circuit is \mu > 1 // v_0 (grounded) // Create Intial points - sysmodel->y()[0] = vinit; // v_1 - sysmodel->y()[1] = vinit; // v_2 - sysmodel->y()[2] = 0.0; // i_L - sysmodel->y()[3] = 0.0; // i_s + sysmodel.y()[0] = vinit; // v_1 + sysmodel.y()[1] = vinit; // v_2 + sysmodel.y()[2] = 0.0; // i_L + sysmodel.y()[3] = 0.0; // i_s - sysmodel->yp()[0] = 0.0; // v'_1 - sysmodel->yp()[1] = 0.0; // v'_2 - sysmodel->yp()[2] = -vinit / linit; // i'_s - sysmodel->yp()[3] = -vinit / linit; // i'_L + sysmodel.yp()[0] = 0.0; // v'_1 + sysmodel.yp()[1] = 0.0; // v'_2 + sysmodel.yp()[2] = -vinit / linit; // i'_s + sysmodel.yp()[3] = -vinit / linit; // i'_L - sysmodel->initialize(); - sysmodel->evaluateResidual(); + sysmodel.initialize(); + sysmodel.evaluateResidual(); std::cout << "Verify Intial Resisdual is Zero: {"; - for (double i : sysmodel->getResidual()) + for (double i : sysmodel.getResidual()) { std::cout << i << ", "; } std::cout << "}\n"; - sysmodel->updateTime(0.0, 1.0); - sysmodel->evaluateJacobian(); + sysmodel.updateTime(0.0, 1.0); + sysmodel.evaluateJacobian(); std::cout << "Intial Jacobian with alpha = 1:\n"; - sysmodel->getJacobian().printMatrix(); + sysmodel.getJacobian().printMatrix(); // Create numerical integrator and configure it for the generator model - AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(sysmodel); + AnalysisManager::Sundials::Ida idas(&sysmodel); double t_init = 0.0; double t_final = 1.0; // setup simulation - idas->configureSimulation(); - idas->getDefaultInitialCondition(); - idas->initializeSimulation(t_init); + idas.configureSimulation(); + idas.getDefaultInitialCondition(); + idas.initializeSimulation(t_init); - idas->runSimulation(t_final); + idas.runSimulation(t_final); - std::vector& yfinial = sysmodel->y(); + std::vector& yfinial = sysmodel.y(); std::cout << "Final Vector y\n"; for (size_t i = 0; i < yfinial.size(); i++) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index d00df651..5bde5db4 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -40,12 +40,8 @@ namespace AnalysisManager { deleteQuadrature(); deleteAdjoint(); - N_VDestroy(yyB_); - N_VDestroy(ypB_); - N_VDestroy(qB_); - SUNLinSolFree(linearSolverB_); - SUNMatDestroy(JacobianMatB_); deleteSimulation(); + deleteBackwardSimulation(); SUNContext_Free(&context_); } @@ -196,6 +192,27 @@ namespace AnalysisManager return retval; } + /** + * @brief Run the IDA solver on the given model and produce a solution at + * the given final time. + * + * @tparam ScalarT Scalar data type + * @tparam IdxT Matrix and vector index data type + * @param tf The final simulation time. + * @param nout The number of integration segmentstimes. + * @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. + * @return int zero if successful, error code otherwise. + * + * @note 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`. + * + * @todo Consider adding initial time as the function argument, as well. + */ template int Ida::runSimulation(real_type tf, int nout, const std::optional> step_callback) { @@ -205,8 +222,8 @@ namespace AnalysisManager 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. */ + // In loop, call IDASolve, print results, and test for error. + // Break out of loop when NOUT preset output times have been reached. // printOutput(0.0); while (nout > iout) { @@ -215,8 +232,9 @@ namespace AnalysisManager 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) + // 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()); @@ -514,6 +532,18 @@ namespace AnalysisManager return 0; } + template + int Ida::deleteBackwardSimulation() + { + N_VDestroy(yyB_); + N_VDestroy(ypB_); + N_VDestroy(qB_); + SUNLinSolFree(linearSolverB_); + SUNMatDestroy(JacobianMatB_); + + return 0; + } + template int Ida::Residual(real_type tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) { diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index ac151963..cd42fe23 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -36,18 +36,6 @@ namespace AnalysisManager int setIntegrationTime(real_type t_init, real_type t_final, int nout); int initializeSimulation(real_type t0, bool findConsistent = false); - /// @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(); @@ -63,6 +51,7 @@ namespace AnalysisManager int runForwardSimulation(real_type tf, int nout = 1); int runBackwardSimulation(real_type t0); int deleteAdjoint(); + int deleteBackwardSimulation(); int saveInitialCondition() {