diff --git a/.github/workflows/spack_default_build.yaml b/.github/workflows/spack_default_build.yaml index fcdb28c9b..eb345aefb 100644 --- a/.github/workflows/spack_default_build.yaml +++ b/.github/workflows/spack_default_build.yaml @@ -101,6 +101,7 @@ jobs: spack_spec: - gridkit@develop +enzyme ^enzyme@0.0.173 + ^sundials@develop steps: - name: Add LLVM diff --git a/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp b/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp index bab5f20df..0565734ed 100644 --- a/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp +++ b/examples/Experimental/AdjointSensitivity/AdjointSensitivity.cpp @@ -138,5 +138,10 @@ int main() std::cout << "The two results differ beyond solver tolerance!\n"; } + delete idas; + delete bus; + delete gen; + delete model; + return retval; } diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index 841395c3a..1bea6597b 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; @@ -32,7 +32,7 @@ int main(int argc, char const* argv[]) double vinit = 1.0; // inductor - GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit); + GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit); // Form index to node uid realations // input induct->setExternalConnectionNodes(0, 1); @@ -41,22 +41,22 @@ int main(int argc, char const* argv[]) // internal induct->setExternalConnectionNodes(2, 2); // add component - sysmodel->addComponent(induct); + sysmodel.addComponent(induct); // resistor idoff++; - GridKit::Resistor* resis = new GridKit::Resistor(idoff, rinit); + GridKit::Resistor* resis = new GridKit::Resistor(idoff, rinit); // Form index to node uid realations // input resis->setExternalConnectionNodes(0, 0); // output resis->setExternalConnectionNodes(1, 1); // add - sysmodel->addComponent(resis); + sysmodel.addComponent(resis); // voltage source idoff++; - GridKit::VoltageSource* vsource = new GridKit::VoltageSource(idoff, vinit); + GridKit::VoltageSource* vsource = new GridKit::VoltageSource(idoff, vinit); // Form index to node uid realations // input vsource->setExternalConnectionNodes(0, -1); @@ -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 7cce626df..5bde5db44 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -25,7 +25,6 @@ namespace AnalysisManager retval = SUNContext_Create(SUN_COMM_NULL, &context_); checkOutput(retval, "SUNContext"); solver_ = IDACreate(context_); - tag_ = NULL; } /** @@ -39,22 +38,10 @@ namespace AnalysisManager template Ida::~Ida() { - N_VDestroy(yy_); - N_VDestroy(yp_); - N_VDestroy(yy0_); - N_VDestroy(yp0_); - if (model_->hasJacobian()) - { - SUNLinSolFree_KLU(linearSolver_); - SUNMatDestroy_Sparse(JacobianMat_); - } - else - { - SUNLinSolFree_Dense(linearSolver_); - SUNMatDestroy_Dense(JacobianMat_); - } - ///@todo this free is needed but on geninfbus this seg faults - // IDAFree(&solver_); + deleteQuadrature(); + deleteAdjoint(); + deleteSimulation(); + deleteBackwardSimulation(); SUNContext_Free(&context_); } @@ -79,7 +66,7 @@ namespace AnalysisManager checkAllocation((void*) yp0_, "N_VClone"); // Dummy initial time; will be overridden. - const sunrealtype t0 = SUN_RCONST(0.0); + const real_type t0 = 0.0; // Allocate and initialize IDA workspace retval = IDAInit(solver_, this->Residual, t0, yy_, yp_); @@ -90,8 +77,8 @@ namespace AnalysisManager checkOutput(retval, "IDASetUserData"); // Set tolerances - sunrealtype rel_tol; - sunrealtype abs_tol; + real_type rel_tol; + real_type abs_tol; model_->setTolerances(rel_tol, abs_tol); ///< \todo Function name should be "getTolerances"! retval = IDASStolerances(solver_, rel_tol, abs_tol); @@ -120,9 +107,7 @@ namespace AnalysisManager } // Set up linear solver - this->configureLinearSolver(); - - return retval; + return this->configureLinearSolver(); } template @@ -207,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) { @@ -216,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) { @@ -226,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()); @@ -242,7 +249,7 @@ namespace AnalysisManager } } - // Final copy out. No gaurentee last residual evaluation is final step. + // Final copy out. No guarantee last residual evaluation is final step. model_->updateTime(tf, 0.0); copyVec(yy_, model_->y()); copyVec(yp_, model_->yp()); @@ -254,10 +261,14 @@ namespace AnalysisManager template int Ida::deleteSimulation() { - IDAFree(&solver_); - SUNLinSolFree(linearSolver_); N_VDestroy(yy_); N_VDestroy(yp_); + N_VDestroy(tag_); + N_VDestroy(yy0_); + N_VDestroy(yp0_); + SUNLinSolFree(linearSolver_); + SUNMatDestroy(JacobianMat_); + IDAFree(&solver_); return 0; } @@ -295,7 +306,7 @@ namespace AnalysisManager int retval = 0; // Set all quadratures to zero - N_VConst(SUN_RCONST(0.0), q_); + N_VConst(0.0, q_); // Initialize quadratures retval = IDAQuadReInit(solver_, q_); @@ -380,9 +391,9 @@ namespace AnalysisManager template int Ida::initializeBackwardSimulation(real_type tf) { - int retval = 0; - sunrealtype rel_tol; - sunrealtype abs_tol; + int retval = 0; + real_type rel_tol; + real_type abs_tol; model_->initializeAdjoint(); @@ -522,7 +533,19 @@ namespace AnalysisManager } template - int Ida::Residual(sunrealtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) + 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) { GridKit::Model::Evaluator* model = static_cast*>(user_data); @@ -538,7 +561,7 @@ namespace AnalysisManager } template - int Ida::Jac(sunrealtype t, sunrealtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) + int Ida::Jac(real_type t, real_type cj, N_Vector yy, N_Vector yp, N_Vector, SUNMatrix J, void* user_data, N_Vector, N_Vector, N_Vector) { GridKit::Model::Evaluator* model = static_cast*>(user_data); @@ -561,25 +584,19 @@ namespace AnalysisManager // Set row pointers sunindextype* rowptrs = SUNSparseMatrix_IndexPointers(J); - for (unsigned int i = 0; i < csrrowdata.size(); i++) - { - rowptrs[i] = csrrowdata[i]; - } + std::copy(csrrowdata.cbegin(), csrrowdata.cend(), rowptrs); sunindextype* colvals = SUNSparseMatrix_IndexValues(J); - sunrealtype* data = SUNSparseMatrix_Data(J); + real_type* data = SUNSparseMatrix_Data(J); // Copy data from model jac to sundials - for (unsigned int i = 0; i < c.size(); i++) - { - colvals[i] = c[i]; - data[i] = val[i]; - } + std::copy(c.cbegin(), c.cend(), colvals); + std::copy(val.cbegin(), val.cend(), data); return 0; } template - int Ida::Integrand(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void* user_data) + int Ida::Integrand(real_type tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void* user_data) { GridKit::Model::Evaluator* model = static_cast*>(user_data); @@ -595,7 +612,7 @@ namespace AnalysisManager } template - int Ida::adjointResidual(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void* user_data) + int Ida::adjointResidual(real_type tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rrB, void* user_data) { GridKit::Model::Evaluator* model = static_cast*>(user_data); @@ -613,7 +630,7 @@ namespace AnalysisManager } template - int Ida::adjointIntegrand(sunrealtype tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void* user_data) + int Ida::adjointIntegrand(real_type tt, N_Vector yy, N_Vector yp, N_Vector yyB, N_Vector ypB, N_Vector rhsQB, void* user_data) { GridKit::Model::Evaluator* model = static_cast*>(user_data); @@ -633,41 +650,29 @@ namespace AnalysisManager template void Ida::copyVec(const N_Vector x, std::vector& y) { - const ScalarT* xdata = NV_DATA_S(x); - for (unsigned int i = 0; i < y.size(); ++i) - { - y[i] = xdata[i]; - } + const ScalarT* xdata = N_VGetArrayPointer(x); + std::copy_n(xdata, y.size(), y.begin()); } template void Ida::copyVec(const std::vector& x, N_Vector y) { - ScalarT* ydata = NV_DATA_S(y); - for (unsigned int i = 0; i < x.size(); ++i) - { - ydata[i] = x[i]; - } + ScalarT* ydata = N_VGetArrayPointer(y); + std::copy(x.cbegin(), x.cend(), ydata); } template void Ida::copyVec(const std::vector& x, N_Vector y) { - ScalarT* ydata = NV_DATA_S(y); - for (unsigned int i = 0; i < x.size(); ++i) - { - if (x[i]) - ydata[i] = 1.0; - else - ydata[i] = 0.0; - } + ScalarT* ydata = N_VGetArrayPointer(y); + std::copy(x.cbegin(), x.cend(), ydata); } template - void Ida::printOutput(sunrealtype t) + void Ida::printOutput(real_type t) { - sunrealtype* yval = N_VGetArrayPointer_Serial(yy_); - sunrealtype* ypval = N_VGetArrayPointer_Serial(yp_); + real_type* yval = N_VGetArrayPointer(yy_); + real_type* ypval = N_VGetArrayPointer(yp_); std::cout << std::setprecision(5) << std::setw(7) << t << " "; for (IdxT i = 0; i < model_->size(); ++i) @@ -682,10 +687,10 @@ namespace AnalysisManager } template - void Ida::printSpecial(sunrealtype t, N_Vector y) + void Ida::printSpecial(real_type t, N_Vector y) { - sunrealtype* yval = N_VGetArrayPointer_Serial(y); - IdxT N = static_cast(N_VGetLength_Serial(y)); + real_type* yval = N_VGetArrayPointer(y); + IdxT N = static_cast(N_VGetLength(y)); std::cout << "{"; std::cout << std::setprecision(5) << std::setw(7) << t; for (IdxT i = 0; i < N; ++i) @@ -698,35 +703,8 @@ namespace AnalysisManager template void Ida::printFinalStats() { - int retval = 0; - void* mem = solver_; - long int nst; - long int nre; - long int nje; - long int nni; - long int netf; - long int ncfn; - - retval = IDAGetNumSteps(mem, &nst); - checkOutput(retval, "IDAGetNumSteps"); - retval = IDAGetNumResEvals(mem, &nre); - checkOutput(retval, "IDAGetNumResEvals"); - retval = IDAGetNumJacEvals(mem, &nje); - checkOutput(retval, "IDAGetNumJacEvals"); - retval = IDAGetNumNonlinSolvIters(mem, &nni); - checkOutput(retval, "IDAGetNumNonlinSolvIters"); - retval = IDAGetNumErrTestFails(mem, &netf); - checkOutput(retval, "IDAGetNumErrTestFails"); - retval = IDAGetNumNonlinSolvConvFails(mem, &ncfn); - checkOutput(retval, "IDAGetNumNonlinSolvConvFails"); - - // std::cout << "\nFinal Run Statistics: \n\n"; - std::cout << "Number of steps = " << nst << "\n"; - std::cout << "Number of residual evaluations = " << nre << "\n"; - // std::cout << "Number of Jacobian evaluations = " << nje << "\n"; - std::cout << "Number of nonlinear iterations = " << nni << "\n"; - std::cout << "Number of error test failures = " << netf << "\n"; - std::cout << "Number of nonlinear conv. failures = " << ncfn << "\n"; + int retval = IDAPrintAllStats(solver_, stdout, SUN_OUTPUTFORMAT_TABLE); + checkOutput(retval, "IDAPrintAllStats"); } template diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index e6d07f579..cd42fe232 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() { @@ -95,93 +84,93 @@ namespace AnalysisManager const real_type* getIntegral() const { - return NV_DATA_S(q_); + return N_VGetArrayPointer(q_); } real_type* getIntegral() { - return NV_DATA_S(q_); + return N_VGetArrayPointer(q_); } const real_type* getAdjointIntegral() const { - return NV_DATA_S(qB_); + return N_VGetArrayPointer(qB_); } real_type* getAdjointIntegral() { - return NV_DATA_S(qB_); + return N_VGetArrayPointer(qB_); } - void printOutput(sunrealtype t); - void printSpecial(sunrealtype t, N_Vector x); + void printOutput(real_type t); + void printSpecial(real_type t, N_Vector x); void printFinalStats(); private: - static int Residual(sunrealtype t, - N_Vector yy, - N_Vector yp, - N_Vector rr, - void* user_data); - - static int Jac(sunrealtype t, - sunrealtype cj, - N_Vector yy, - N_Vector yp, - N_Vector resvec, - SUNMatrix J, - void* user_data, - N_Vector tmp1, - N_Vector tmp2, - N_Vector tmp3); - - static int Integrand(sunrealtype t, - N_Vector yy, - N_Vector yp, - N_Vector rhsQ, - void* user_data); - - static int adjointResidual(sunrealtype t, - N_Vector yy, - N_Vector yp, - N_Vector yyB, - N_Vector ypB, - N_Vector rrB, - void* user_data); - - static int adjointIntegrand(sunrealtype t, - N_Vector yy, - N_Vector yp, - N_Vector yyB, - N_Vector ypB, - N_Vector rhsQB, - void* user_data); + static int Residual(real_type t, + N_Vector yy, + N_Vector yp, + N_Vector rr, + void* user_data); + + static int Jac(real_type t, + real_type cj, + N_Vector yy, + N_Vector yp, + N_Vector resvec, + SUNMatrix J, + void* user_data, + N_Vector tmp1, + N_Vector tmp2, + N_Vector tmp3); + + static int Integrand(real_type t, + N_Vector yy, + N_Vector yp, + N_Vector rhsQ, + void* user_data); + + static int adjointResidual(real_type t, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rrB, + void* user_data); + + static int adjointIntegrand(real_type t, + N_Vector yy, + N_Vector yp, + N_Vector yyB, + N_Vector ypB, + N_Vector rhsQB, + void* user_data); private: - void* solver_; - SUNContext context_; - SUNMatrix JacobianMat_; - SUNMatrix JacobianMatB_; - SUNLinearSolver linearSolver_; - SUNLinearSolver linearSolverB_; - - real_type t_init_; - real_type t_final_; - int nout_; ///< Number of integration outputs - - N_Vector yy_; ///< Solution vector - N_Vector yp_; ///< Solution derivatives vector - N_Vector tag_; ///< Tags differential variables - N_Vector q_; ///< Integrand vector - - N_Vector yy0_; ///< Storage for initial values - N_Vector yp0_; ///< Storage for initial derivatives - - N_Vector yyB_; ///< Adjoint solution vector - N_Vector ypB_; ///< Adjoint solution derivatives vector - N_Vector qB_; ///< Backward integrand vector - - int backwardID_; + void* solver_{}; + SUNContext context_{}; + SUNMatrix JacobianMat_{}; + SUNMatrix JacobianMatB_{}; + SUNLinearSolver linearSolver_{}; + SUNLinearSolver linearSolverB_{}; + + real_type t_init_{}; + real_type t_final_{}; + int nout_{}; ///< Number of integration outputs + + N_Vector yy_{}; ///< Solution vector + N_Vector yp_{}; ///< Solution derivatives vector + N_Vector tag_{}; ///< Tags differential variables + N_Vector q_{}; ///< Integrand vector + + N_Vector yy0_{}; ///< Storage for initial values + N_Vector yp0_{}; ///< Storage for initial derivatives + + N_Vector yyB_{}; ///< Adjoint solution vector + N_Vector ypB_{}; ///< Adjoint solution derivatives vector + N_Vector qB_{}; ///< Backward integrand vector + + int backwardID_{}; private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); diff --git a/src/Solver/Optimization/DynamicConstraint.cpp b/src/Solver/Optimization/DynamicConstraint.cpp index 5d2e04582..784d4e3f6 100644 --- a/src/Solver/Optimization/DynamicConstraint.cpp +++ b/src/Solver/Optimization/DynamicConstraint.cpp @@ -28,7 +28,7 @@ namespace AnalysisManager bool DynamicConstraint::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag, IndexStyleEnum& index_style) { // This code handles one objective function - assert(model_->size_quad() == 1); + assert(model_->sizeQuadrature() == 1); // Number of parameters is size of the system plus 1 fictitious parameter // to store the objective value. diff --git a/src/Solver/Optimization/DynamicObjective.cpp b/src/Solver/Optimization/DynamicObjective.cpp index 0b76d82e0..81ad25cde 100644 --- a/src/Solver/Optimization/DynamicObjective.cpp +++ b/src/Solver/Optimization/DynamicObjective.cpp @@ -28,7 +28,7 @@ namespace AnalysisManager bool DynamicObjective::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag, IndexStyleEnum& index_style) { // This code handles one objective function - assert(model_->size_quad() == 1); + assert(model_->sizeQuadrature() == 1); // Number of optimization variables. n = model_->sizeParams(); diff --git a/src/Solver/SteadyState/Kinsol.cpp b/src/Solver/SteadyState/Kinsol.cpp index ef39bf04e..67f33982d 100644 --- a/src/Solver/SteadyState/Kinsol.cpp +++ b/src/Solver/SteadyState/Kinsol.cpp @@ -37,23 +37,14 @@ namespace AnalysisManager checkOutput(retval, "SUNContext"); solver_ = KINCreate(context_); - tag_ = NULL; } template Kinsol::~Kinsol() { + deleteSimulation(); SUNContext_Free(&context_); - KINFree(&solver_); - - N_VDestroy_Serial(this->yy_); - N_VDestroy_Serial(this->yy0_); - N_VDestroy_Serial(this->scale_); - - SUNMatDestroy(this->JacobianMat_); - SUNLinSolFree_Dense(this->linearSolver_); - - solver_ = 0; + solver_ = nullptr; } template @@ -93,16 +84,7 @@ namespace AnalysisManager checkOutput(retval, "KINSetScaledStepTol"); // Set up linear solver - JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_); - checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); - - linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); - checkAllocation((void*) linearSolver_, "SUNLinSol_Dense"); - - retval = KINSetLinearSolver(solver_, linearSolver_, JacobianMat_); - checkOutput(retval, "KINSetLinearSolver"); - - return retval; + return this->configureLinearSolver(); } template @@ -148,10 +130,12 @@ namespace AnalysisManager template int Kinsol::deleteSimulation() { - SUNLinSolFree(linearSolver_); KINFree(&solver_); - N_VDestroy(yy_); - N_VDestroy(scale_); + N_VDestroy(this->yy_); + N_VDestroy(this->yy0_); + N_VDestroy(this->scale_); + SUNMatDestroy(this->JacobianMat_); + SUNLinSolFree_Dense(this->linearSolver_); return 0; } @@ -173,40 +157,21 @@ namespace AnalysisManager template void Kinsol::copyVec(const N_Vector x, std::vector& y) { - const ScalarT* xdata = NV_DATA_S(x); - for (unsigned int i = 0; i < y.size(); ++i) - { - y[i] = xdata[i]; - } + const ScalarT* xdata = N_VGetArrayPointer(x); + std::copy_n(xdata, y.size(), y.begin()); } template void Kinsol::copyVec(const std::vector& x, N_Vector y) { - ScalarT* ydata = NV_DATA_S(y); - for (unsigned int i = 0; i < x.size(); ++i) - { - ydata[i] = x[i]; - } - } - - template - void Kinsol::copyVec(const std::vector& x, N_Vector y) - { - ScalarT* ydata = NV_DATA_S(y); - for (unsigned int i = 0; i < x.size(); ++i) - { - if (x[i]) - ydata[i] = 1.0; - else - ydata[i] = 0.0; - } + ScalarT* ydata = N_VGetArrayPointer(y); + std::copy_n(x.cbegin(), x.size(), ydata); } template void Kinsol::printOutput() { - sunrealtype* yval = N_VGetArrayPointer_Serial(yy_); + sunrealtype* yval = N_VGetArrayPointer(yy_); std::cout << std::setprecision(5) << std::setw(7); for (IdxT i = 0; i < model_->size(); ++i) @@ -220,7 +185,7 @@ namespace AnalysisManager void Kinsol::printSpecial(sunrealtype t, N_Vector y) { sunrealtype* yval = N_VGetArrayPointer_Serial(y); - IdxT N = N_VGetLength_Serial(y); + IdxT N = static_cast(N_VGetLength_Serial(y)); std::cout << "{"; std::cout << std::setprecision(5) << std::setw(7) << t; for (IdxT i = 0; i < N; ++i) @@ -233,29 +198,8 @@ namespace AnalysisManager template void Kinsol::printFinalStats() { - int retval = 0; - void* mem = solver_; - long int nni; - long int nfe; - long int nje; - long int nlfe; - - // retval = KINGetNumSteps(mem, &nst); - // checkOutput(retval, "KINGetNumSteps"); - retval = KINGetNumNonlinSolvIters(mem, &nni); - checkOutput(retval, "KINGetNumNonlinSolvIters"); - retval = KINGetNumFuncEvals(mem, &nfe); - checkOutput(retval, "KINGetNumFuncEvals"); - retval = KINGetNumJacEvals(mem, &nje); - checkOutput(retval, "KINGetNumJacEvals"); - retval = KINGetNumLinFuncEvals(mem, &nlfe); - checkOutput(retval, "KINGetNumLinFuncEvals"); - - // std::cout << "\nFinal Run Statistics: \n\n"; - std::cout << "Number of nonlinear iterations = " << nni << "\n"; - std::cout << "Number of function evaluations = " << nfe << "\n"; - std::cout << "Number of Jacobian evaluations = " << nje << "\n"; - std::cout << "Number of linear function evals. = " << nlfe << "\n"; + int retval = KINPrintAllStats(solver_, stdout, SUN_OUTPUTFORMAT_TABLE); + checkOutput(retval, "IDAPrintAllStats"); } template diff --git a/src/Solver/SteadyState/Kinsol.hpp b/src/Solver/SteadyState/Kinsol.hpp index b4cb190a2..c43feaee8 100644 --- a/src/Solver/SteadyState/Kinsol.hpp +++ b/src/Solver/SteadyState/Kinsol.hpp @@ -31,6 +31,8 @@ namespace AnalysisManager typedef typename GridKit::ScalarTraits::real_type real_type; + static_assert(std::is_same_v, "real_type must be the same type as sunrealtype"); + public: Kinsol(GridKit::Model::Evaluator* model); ~Kinsol(); @@ -85,22 +87,22 @@ namespace AnalysisManager // const real_type* getIntegral() const // { - // return NV_DATA_S(q_); + // return N_VGetArrayPointer(q_); // } // real_type* getIntegral() // { - // return NV_DATA_S(q_); + // return N_VGetArrayPointer(q_); // } // const real_type* getAdjointIntegral() const // { - // return NV_DATA_S(qB_); + // return N_VGetArrayPointer(qB_); // } // real_type* getAdjointIntegral() // { - // return NV_DATA_S(qB_); + // return N_VGetArrayPointer(qB_); // } void printOutput(); @@ -125,17 +127,16 @@ namespace AnalysisManager // N_Vector rhsQB, void *user_data); private: - void* solver_; - SUNContext context_; - SUNMatrix JacobianMat_; - SUNLinearSolver linearSolver_; + void* solver_{}; + SUNContext context_{}; + SUNMatrix JacobianMat_{}; + SUNLinearSolver linearSolver_{}; - N_Vector yy_; ///< Solution vector - N_Vector scale_; ///< Scaling vector - N_Vector tag_; ///< Tags differential variables - N_Vector q_; ///< Integrand vector + N_Vector yy_{}; ///< Solution vector + N_Vector scale_{}; ///< Scaling vector + N_Vector q_{}; ///< Integrand vector - N_Vector yy0_; ///< Storage for initial values + N_Vector yy0_{}; ///< Storage for initial values private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida);