From 00ea0a1b3cffc6c91fb74526852d455e5007b4ff Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 9 May 2025 13:58:05 -0700 Subject: [PATCH 01/14] Improvements for IDA and KINSOL --- src/Solver/Dynamic/Ida.cpp | 124 ++++++++-------------------- src/Solver/Dynamic/Ida.hpp | 131 ++++++++++++++++-------------- src/Solver/SteadyState/Kinsol.cpp | 62 ++++---------- src/Solver/SteadyState/Kinsol.hpp | 28 ++++--- 4 files changed, 132 insertions(+), 213 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 7cce626df..b92ec815d 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -9,6 +9,7 @@ #include "Model/Evaluator.hpp" +//TODO: move implementations to header since they are templated? namespace AnalysisManager { @@ -25,7 +26,6 @@ namespace AnalysisManager retval = SUNContext_Create(SUN_COMM_NULL, &context_); checkOutput(retval, "SUNContext"); solver_ = IDACreate(context_); - tag_ = NULL; } /** @@ -43,18 +43,9 @@ namespace AnalysisManager 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_); + SUNLinSolFree(linearSolver_); + SUNMatDestroy(JacobianMat_); + IDAFree(&solver_); SUNContext_Free(&context_); } @@ -79,7 +70,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 +81,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); @@ -172,6 +163,7 @@ namespace AnalysisManager template int Ida::setIntegrationTime(real_type t_init, real_type t_final, int nout) { + //TODO: I don't see these variables used anywhere. Remove? t_init_ = t_init; t_final_ = t_final; nout_ = nout; @@ -183,8 +175,6 @@ 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"); @@ -197,6 +187,7 @@ namespace AnalysisManager if (tag_) initType = IDA_YA_YDP_INIT; + // TODO: this doesn't doing anything without calling IDAGetConsistentIC retval = IDACalcIC(solver_, initType, 0.1); checkOutput(retval, "IDACalcIC"); @@ -242,7 +233,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()); @@ -295,7 +286,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_); @@ -381,8 +372,8 @@ namespace AnalysisManager int Ida::initializeBackwardSimulation(real_type tf) { int retval = 0; - sunrealtype rel_tol; - sunrealtype abs_tol; + real_type rel_tol; + real_type abs_tol; model_->initializeAdjoint(); @@ -522,7 +513,7 @@ namespace AnalysisManager } template - int Ida::Residual(sunrealtype tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data) + 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 +529,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 +552,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 +580,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 +598,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 +618,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_Serial(yy_); + real_type* ypval = N_VGetArrayPointer_Serial(yp_); std::cout << std::setprecision(5) << std::setw(7) << t << " "; for (IdxT i = 0; i < model_->size(); ++i) @@ -682,9 +655,9 @@ 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); + real_type* yval = N_VGetArrayPointer_Serial(y); IdxT N = static_cast(N_VGetLength_Serial(y)); std::cout << "{"; std::cout << std::setprecision(5) << std::setw(7) << t; @@ -698,35 +671,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..677479e01 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -26,6 +26,10 @@ 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"); + // TODO: can't include this assert because cpp file instantiates a few values of IdxT + // static_assert(std::is_same_v, "IdxT must be the same type as sunindextype"); + public: Ida(GridKit::Model::Evaluator* model); ~Ida(); @@ -95,96 +99,97 @@ 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_; + 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 + 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 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 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 + N_Vector yyB_{}; ///< Adjoint solution vector + N_Vector ypB_{}; ///< Adjoint solution derivatives vector + N_Vector qB_{}; ///< Backward integrand vector - int backwardID_; + int backwardID_{}; private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); + //TODO: should template type be real_type rather than ScalarT? static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); static void copyVec(const std::vector& x, N_Vector y); diff --git a/src/Solver/SteadyState/Kinsol.cpp b/src/Solver/SteadyState/Kinsol.cpp index ef39bf04e..b7e09a0eb 100644 --- a/src/Solver/SteadyState/Kinsol.cpp +++ b/src/Solver/SteadyState/Kinsol.cpp @@ -37,7 +37,6 @@ namespace AnalysisManager checkOutput(retval, "SUNContext"); solver_ = KINCreate(context_); - tag_ = NULL; } template @@ -46,14 +45,14 @@ namespace AnalysisManager SUNContext_Free(&context_); KINFree(&solver_); - N_VDestroy_Serial(this->yy_); - N_VDestroy_Serial(this->yy0_); - N_VDestroy_Serial(this->scale_); + N_VDestroy(this->yy_); + N_VDestroy(this->yy0_); + N_VDestroy(this->scale_); SUNMatDestroy(this->JacobianMat_); SUNLinSolFree_Dense(this->linearSolver_); - solver_ = 0; + solver_ = nullptr; } template @@ -173,40 +172,28 @@ 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]; - } + ScalarT* ydata = N_VGetArrayPointer(y); + std::copy_n(x.cbegin(), x.size(), ydata); } 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(x.cbegin(), x.cend(), 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 +207,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 +220,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..028d8e05e 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,17 @@ 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 tag_{}; ///< Tags differential variables + 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); From 5ba5452a7f34ab74b28b983b789eaef6c55eb7ff Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 9 May 2025 14:32:36 -0700 Subject: [PATCH 02/14] Remove N_VGetArrayPointer_Serial calls --- src/Solver/Dynamic/Ida.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index b92ec815d..effc59c18 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -639,8 +639,8 @@ namespace AnalysisManager template void Ida::printOutput(real_type t) { - real_type* yval = N_VGetArrayPointer_Serial(yy_); - real_type* 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) @@ -657,8 +657,8 @@ namespace AnalysisManager template void Ida::printSpecial(real_type t, N_Vector y) { - real_type* 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) From 4ade8cb71690eb230434d11b943d73c3b2245941 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 9 May 2025 14:48:03 -0700 Subject: [PATCH 03/14] Resolve t_init_ issue --- src/Solver/Dynamic/Ida.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index effc59c18..d613e713c 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -163,7 +163,6 @@ namespace AnalysisManager template int Ida::setIntegrationTime(real_type t_init, real_type t_final, int nout) { - //TODO: I don't see these variables used anywhere. Remove? t_init_ = t_init; t_final_ = t_final; nout_ = nout; @@ -175,6 +174,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"); From cd467dc1be0417cb137ad65869a58dbe5397ede6 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Fri, 9 May 2025 15:01:00 -0700 Subject: [PATCH 04/14] Add missing N_VDestroy --- src/Solver/Dynamic/Ida.cpp | 3 ++- src/Solver/Dynamic/Ida.hpp | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index d613e713c..9565db78c 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -99,12 +99,13 @@ namespace AnalysisManager std::vector& tag = model_->tag(); if (static_cast(tag.size()) == model_->size()) { - tag_ = N_VClone(yy_); + N_Vector tag_ = N_VClone(yy_); checkAllocation((void*) tag_, "N_VClone"); model_->tagDifferentiable(); copyVec(tag, tag_); retval = IDASetId(solver_, tag_); + N_VDestroy(tag_); checkOutput(retval, "IDASetId"); retval = IDASetSuppressAlg(solver_, SUNTRUE); checkOutput(retval, "IDASetSuppressAlg"); diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index 677479e01..7d5ef6b45 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -175,7 +175,6 @@ namespace AnalysisManager 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 From c3560db8f28eddcfc1a9dcb752b2d1e87f3d5dcd Mon Sep 17 00:00:00 2001 From: Steven-Roberts Date: Fri, 9 May 2025 22:04:09 +0000 Subject: [PATCH 05/14] Apply pre-commmit fixes --- src/Solver/Dynamic/Ida.cpp | 6 +++--- src/Solver/Dynamic/Ida.hpp | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 9565db78c..fea9f3e73 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -9,7 +9,7 @@ #include "Model/Evaluator.hpp" -//TODO: move implementations to header since they are templated? +// TODO: move implementations to header since they are templated? namespace AnalysisManager { @@ -373,7 +373,7 @@ namespace AnalysisManager template int Ida::initializeBackwardSimulation(real_type tf) { - int retval = 0; + int retval = 0; real_type rel_tol; real_type abs_tol; @@ -660,7 +660,7 @@ namespace AnalysisManager void Ida::printSpecial(real_type t, N_Vector y) { real_type* yval = N_VGetArrayPointer(y); - IdxT N = static_cast(N_VGetLength(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) diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index 7d5ef6b45..b05b29bc0 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -173,9 +173,9 @@ namespace AnalysisManager real_type t_final_{}; int nout_{}; ///< Number of integration outputs - N_Vector yy_{}; ///< Solution vector - N_Vector yp_{}; ///< Solution derivatives vector - N_Vector q_{}; ///< Integrand vector + N_Vector yy_{}; ///< Solution vector + N_Vector yp_{}; ///< Solution derivatives vector + N_Vector q_{}; ///< Integrand vector N_Vector yy0_{}; ///< Storage for initial values N_Vector yp0_{}; ///< Storage for initial derivatives @@ -188,7 +188,7 @@ namespace AnalysisManager private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); - //TODO: should template type be real_type rather than ScalarT? + // TODO: should template type be real_type rather than ScalarT? static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); static void copyVec(const std::vector& x, N_Vector y); From 106f89ed3c1e221eb8c356dd65e7ce2894200e48 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 14 May 2025 17:40:14 -0700 Subject: [PATCH 06/14] Fix some adjoint memory leaks --- .../AdjointSensitivity/AdjointSensitivity.cpp | 5 ++++ .../PowerElectronics/RLCircuit/RLCircuit.cpp | 6 ++--- src/Solver/Dynamic/Ida.cpp | 27 ++++++++++--------- src/Solver/Dynamic/Ida.hpp | 1 + src/Solver/Optimization/DynamicConstraint.cpp | 2 +- src/Solver/Optimization/DynamicObjective.cpp | 2 +- 6 files changed, 26 insertions(+), 17 deletions(-) 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..885134ecd 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -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); @@ -45,7 +45,7 @@ int main(int argc, char const* argv[]) // 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); @@ -56,7 +56,7 @@ int main(int argc, char const* argv[]) // 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); diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 9565db78c..c86858fdb 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -39,13 +39,9 @@ namespace AnalysisManager template Ida::~Ida() { - N_VDestroy(yy_); - N_VDestroy(yp_); - N_VDestroy(yy0_); - N_VDestroy(yp0_); - SUNLinSolFree(linearSolver_); - SUNMatDestroy(JacobianMat_); - IDAFree(&solver_); + deleteQuadrature(); + deleteAdjoint(); + deleteSimulation(); SUNContext_Free(&context_); } @@ -99,13 +95,12 @@ namespace AnalysisManager std::vector& tag = model_->tag(); if (static_cast(tag.size()) == model_->size()) { - N_Vector tag_ = N_VClone(yy_); + tag_ = N_VClone(yy_); checkAllocation((void*) tag_, "N_VClone"); model_->tagDifferentiable(); copyVec(tag, tag_); retval = IDASetId(solver_, tag_); - N_VDestroy(tag_); checkOutput(retval, "IDASetId"); retval = IDASetSuppressAlg(solver_, SUNTRUE); checkOutput(retval, "IDASetSuppressAlg"); @@ -189,7 +184,6 @@ namespace AnalysisManager if (tag_) initType = IDA_YA_YDP_INIT; - // TODO: this doesn't doing anything without calling IDAGetConsistentIC retval = IDACalcIC(solver_, initType, 0.1); checkOutput(retval, "IDACalcIC"); @@ -247,10 +241,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; } @@ -510,6 +508,11 @@ namespace AnalysisManager template int Ida::deleteAdjoint() { + N_VDestroy(yyB_); + N_VDestroy(ypB_); + N_VDestroy(qB_); + SUNLinSolFree(linearSolverB_); + SUNMatDestroy(JacobianMatB_); IDAAdjFree(solver_); return 0; } diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index 7d5ef6b45..f24b11133 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -175,6 +175,7 @@ namespace AnalysisManager 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 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(); From 1baa96ba95e69870f03b23d3a3058bc9a3c0a774 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 14 May 2025 17:40:59 -0700 Subject: [PATCH 07/14] Remove unused function --- src/Solver/SteadyState/Kinsol.cpp | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/Solver/SteadyState/Kinsol.cpp b/src/Solver/SteadyState/Kinsol.cpp index b7e09a0eb..6a469d86e 100644 --- a/src/Solver/SteadyState/Kinsol.cpp +++ b/src/Solver/SteadyState/Kinsol.cpp @@ -183,13 +183,6 @@ namespace AnalysisManager std::copy_n(x.cbegin(), x.size(), ydata); } - template - void Kinsol::copyVec(const std::vector& x, N_Vector y) - { - ScalarT* ydata = N_VGetArrayPointer(y); - std::copy(x.cbegin(), x.cend(), ydata); - } - template void Kinsol::printOutput() { From 5024cff8a18e40cfe2b051faa4f6ca6e9a5b0e62 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 14 May 2025 17:43:01 -0700 Subject: [PATCH 08/14] merge --- src/Solver/Dynamic/Ida.cpp | 6 +++--- src/Solver/Dynamic/Ida.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index c86858fdb..8be632631 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -9,7 +9,7 @@ #include "Model/Evaluator.hpp" -//TODO: move implementations to header since they are templated? +// TODO: move implementations to header since they are templated? namespace AnalysisManager { @@ -371,7 +371,7 @@ namespace AnalysisManager template int Ida::initializeBackwardSimulation(real_type tf) { - int retval = 0; + int retval = 0; real_type rel_tol; real_type abs_tol; @@ -663,7 +663,7 @@ namespace AnalysisManager void Ida::printSpecial(real_type t, N_Vector y) { real_type* yval = N_VGetArrayPointer(y); - IdxT N = static_cast(N_VGetLength(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) diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index f24b11133..f3c511033 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -175,7 +175,7 @@ namespace AnalysisManager N_Vector yy_{}; ///< Solution vector N_Vector yp_{}; ///< Solution derivatives vector - N_Vector tag_{}; ///< Tags differential variables + N_Vector tag_{}; ///< Tags differential variables N_Vector q_{}; ///< Integrand vector N_Vector yy0_{}; ///< Storage for initial values @@ -189,7 +189,7 @@ namespace AnalysisManager private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); - //TODO: should template type be real_type rather than ScalarT? + // TODO: should template type be real_type rather than ScalarT? static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); static void copyVec(const std::vector& x, N_Vector y); From 6049604b0a479e0050a12d31500b1ad0c757ad8e Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 14 May 2025 20:50:52 -0400 Subject: [PATCH 09/14] Update src/Solver/SteadyState/Kinsol.hpp Co-authored-by: pelesh --- src/Solver/SteadyState/Kinsol.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Solver/SteadyState/Kinsol.hpp b/src/Solver/SteadyState/Kinsol.hpp index 028d8e05e..c43feaee8 100644 --- a/src/Solver/SteadyState/Kinsol.hpp +++ b/src/Solver/SteadyState/Kinsol.hpp @@ -134,7 +134,6 @@ namespace AnalysisManager N_Vector yy_{}; ///< Solution vector N_Vector scale_{}; ///< Scaling vector - N_Vector tag_{}; ///< Tags differential variables N_Vector q_{}; ///< Integrand vector N_Vector yy0_{}; ///< Storage for initial values From 2f42618eb022360feb38fda55e2cf286f58901e9 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 14 May 2025 18:16:58 -0700 Subject: [PATCH 10/14] Use helper functions --- src/Solver/Dynamic/Ida.cpp | 4 +--- src/Solver/SteadyState/Kinsol.cpp | 29 +++++++---------------------- 2 files changed, 8 insertions(+), 25 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 8be632631..10ff770b7 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -107,9 +107,7 @@ namespace AnalysisManager } // Set up linear solver - this->configureLinearSolver(); - - return retval; + return this->configureLinearSolver(); } template diff --git a/src/Solver/SteadyState/Kinsol.cpp b/src/Solver/SteadyState/Kinsol.cpp index 6a469d86e..67f33982d 100644 --- a/src/Solver/SteadyState/Kinsol.cpp +++ b/src/Solver/SteadyState/Kinsol.cpp @@ -42,16 +42,8 @@ namespace AnalysisManager template Kinsol::~Kinsol() { + deleteSimulation(); SUNContext_Free(&context_); - KINFree(&solver_); - - N_VDestroy(this->yy_); - N_VDestroy(this->yy0_); - N_VDestroy(this->scale_); - - SUNMatDestroy(this->JacobianMat_); - SUNLinSolFree_Dense(this->linearSolver_); - solver_ = nullptr; } @@ -92,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 @@ -147,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; } From 0ff5cb50ffd0bf07e75a8057476653f0ee544060 Mon Sep 17 00:00:00 2001 From: Nicholson Koukpaizan Date: Wed, 21 May 2025 13:01:05 -0400 Subject: [PATCH 11/14] sundials@develop in CI. --- .github/workflows/spack_default_build.yaml | 1 + 1 file changed, 1 insertion(+) 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 From b51d3bc79561fa581fce47169e813a27eaa08bc6 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 21 May 2025 17:45:25 -0700 Subject: [PATCH 12/14] Remove todos --- src/Solver/Dynamic/Ida.cpp | 1 - src/Solver/Dynamic/Ida.hpp | 5 ----- 2 files changed, 6 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 10ff770b7..3b1dd2b7c 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -9,7 +9,6 @@ #include "Model/Evaluator.hpp" -// TODO: move implementations to header since they are templated? namespace AnalysisManager { diff --git a/src/Solver/Dynamic/Ida.hpp b/src/Solver/Dynamic/Ida.hpp index f3c511033..ac1519639 100644 --- a/src/Solver/Dynamic/Ida.hpp +++ b/src/Solver/Dynamic/Ida.hpp @@ -26,10 +26,6 @@ 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"); - // TODO: can't include this assert because cpp file instantiates a few values of IdxT - // static_assert(std::is_same_v, "IdxT must be the same type as sunindextype"); - public: Ida(GridKit::Model::Evaluator* model); ~Ida(); @@ -189,7 +185,6 @@ namespace AnalysisManager private: // static void copyMat(Model::Evaluator::Mat& J, SlsMat Jida); - // TODO: should template type be real_type rather than ScalarT? static void copyVec(const N_Vector x, std::vector& y); static void copyVec(const std::vector& x, N_Vector y); static void copyVec(const std::vector& x, N_Vector y); From 5f65419a642db8f0f5e9efc184e3f7834a4fd0ed Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Wed, 21 May 2025 18:17:59 -0700 Subject: [PATCH 13/14] Fix segfault --- src/Solver/Dynamic/Ida.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Solver/Dynamic/Ida.cpp b/src/Solver/Dynamic/Ida.cpp index 3b1dd2b7c..d00df6517 100644 --- a/src/Solver/Dynamic/Ida.cpp +++ b/src/Solver/Dynamic/Ida.cpp @@ -40,6 +40,11 @@ namespace AnalysisManager { deleteQuadrature(); deleteAdjoint(); + N_VDestroy(yyB_); + N_VDestroy(ypB_); + N_VDestroy(qB_); + SUNLinSolFree(linearSolverB_); + SUNMatDestroy(JacobianMatB_); deleteSimulation(); SUNContext_Free(&context_); } @@ -505,11 +510,6 @@ namespace AnalysisManager template int Ida::deleteAdjoint() { - N_VDestroy(yyB_); - N_VDestroy(ypB_); - N_VDestroy(qB_); - SUNLinSolFree(linearSolverB_); - SUNMatDestroy(JacobianMatB_); IDAAdjFree(solver_); return 0; } From 46774cfa5c541dc808cc82c8fcf0f7bd90ddc62e Mon Sep 17 00:00:00 2001 From: pelesh Date: Thu, 22 May 2025 12:31:55 -0400 Subject: [PATCH 14/14] Minor suggestions to SUNDIALS interface (#111) * Minor suggestions to SUNDIALS interface - Create "symmetric" functions for creating and deleting SUNDIALS objects. - Modify RL example to test Ida destructor. --------- Co-authored-by: pelesh --- .../PowerElectronics/RLCircuit/RLCircuit.cpp | 54 +++++++++---------- src/Solver/Dynamic/Ida.cpp | 48 +++++++++++++---- src/Solver/Dynamic/Ida.hpp | 13 +---- 3 files changed, 67 insertions(+), 48 deletions(-) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index 885134ecd..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; @@ -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 d00df6517..5bde5db44 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 ac1519639..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() {