From 507197c75f31aa223e358ae531b66d7c39d4ab7f Mon Sep 17 00:00:00 2001 From: Philip Top Date: Thu, 30 Jul 2026 06:25:14 -0700 Subject: [PATCH 1/3] clean up some possible double initialize in the sundials interfaces --- src/griddyn/solvers/ArkodeInterface.cpp | 1 + src/griddyn/solvers/CvodeInterface.cpp | 1 + src/griddyn/solvers/IdaInterface.cpp | 2 + src/griddyn/solvers/KinsolInterface.cpp | 51 ++++++++++++----------- src/griddyn/solvers/SundialsInterface.cpp | 26 +++++++----- src/griddyn/solvers/SundialsInterface.h | 1 + 6 files changed, 48 insertions(+), 34 deletions(-) diff --git a/src/griddyn/solvers/ArkodeInterface.cpp b/src/griddyn/solvers/ArkodeInterface.cpp index 2dc9a1246..82c8002b8 100644 --- a/src/griddyn/solvers/ArkodeInterface.cpp +++ b/src/griddyn/solvers/ArkodeInterface.cpp @@ -325,6 +325,7 @@ void ArkodeInterface::initialize(CoreTime time0) retval = ARKodeSetMaxNumSteps(solverMem, max_iterations); checkFlag(&retval, "ARKodeSetMaxNumSteps", 1); + freeLinearSolver(); #ifdef ENABLE_KLU if (flags[DENSE_FLAG]) { J = SUNDenseMatrix(svsize, svsize); diff --git a/src/griddyn/solvers/CvodeInterface.cpp b/src/griddyn/solvers/CvodeInterface.cpp index 6b8241b51..0b0f39114 100644 --- a/src/griddyn/solvers/CvodeInterface.cpp +++ b/src/griddyn/solvers/CvodeInterface.cpp @@ -305,6 +305,7 @@ void CvodeInterface::initialize(CoreTime time0) retval = CVodeSetMaxNumSteps(solverMem, max_iterations); checkFlag(&retval, "CVodeSetMaxNumSteps", 1); + freeLinearSolver(); #ifdef GRIDDYN_ENABLE_KLU if (flags[DENSE_FLAG]) { J = SUNDenseMatrix(svsize, svsize, sunctx); diff --git a/src/griddyn/solvers/IdaInterface.cpp b/src/griddyn/solvers/IdaInterface.cpp index c70485dff..23f7256eb 100644 --- a/src/griddyn/solvers/IdaInterface.cpp +++ b/src/griddyn/solvers/IdaInterface.cpp @@ -286,6 +286,8 @@ void IdaInterface::initialize(CoreTime t0) retval = IDASetMaxNumSteps(solverMem, max_iterations); checkFlag(&retval, "IDASetMaxNumSteps", 1); + + freeLinearSolver(); #ifdef GRIDDYN_ENABLE_KLU if (flags[DENSE_FLAG]) { J = SUNDenseMatrix(svsize, svsize, sunctx); diff --git a/src/griddyn/solvers/KinsolInterface.cpp b/src/griddyn/solvers/KinsolInterface.cpp index 8597e2e97..d732bcaf6 100644 --- a/src/griddyn/solvers/KinsolInterface.cpp +++ b/src/griddyn/solvers/KinsolInterface.cpp @@ -96,6 +96,7 @@ void KinsolInterface::allocate(count_t stateCount, count_t /*numRoots*/) if (solverMem != nullptr) { KINFree(&(solverMem)); } + freeLinearSolver(); solverMem = KINCreate(sunctx); checkFlag(solverMem, "KINCreate", 0); @@ -172,40 +173,42 @@ void KinsolInterface::initialize(CoreTime /*t0*/) retval = KINSetNoInitSetup(solverMem, SUNTRUE); checkFlag(&retval, "KINSetNoInitSetup", 1); - retval = KINInit(solverMem, kinsolFunc, state); + if (!flags[INITIALIZED_FLAG]) { + retval = KINInit(solverMem, kinsolFunc, state); - checkFlag(&retval, "KINInit", 1); + checkFlag(&retval, "KINInit", 1); #ifdef GRIDDYN_ENABLE_KLU - if (flags[DENSE_FLAG]) { + if (flags[DENSE_FLAG]) { + J = SUNDenseMatrix(svsize, svsize, sunctx); + checkFlag(J, "SUNDenseMatrix", 0); + /* Create KLU solver object */ + LS = SUNLinSol_Dense(state, J, sunctx); + checkFlag(LS, "SUNLinSol_Dense", 0); + } else { + /* Create sparse SUNMatrix */ + J = SUNSparseMatrix(svsize, svsize, maxNNZ, CSR_MAT, sunctx); + checkFlag(J, "SUNSparseMatrix", 0); + + /* Create KLU solver object */ + LS = SUNLinSol_KLU(state, J, sunctx); + checkFlag(LS, "SUNLinSol_KLU", 0); + + retval = SUNLinSol_KLUSetOrdering(LS, 0); + checkFlag(&retval, "SUNLinSol_KLUSetOrdering", 1); + } +#else J = SUNDenseMatrix(svsize, svsize, sunctx); - checkFlag(J, "SUNDenseMatrix", 0); + checkFlag(J, "SUNSparseMatrix", 0); /* Create KLU solver object */ LS = SUNLinSol_Dense(state, J, sunctx); checkFlag(LS, "SUNLinSol_Dense", 0); - } else { - /* Create sparse SUNMatrix */ - J = SUNSparseMatrix(svsize, svsize, maxNNZ, CSR_MAT, sunctx); - checkFlag(J, "SUNSparseMatrix", 0); - - /* Create KLU solver object */ - LS = SUNLinSol_KLU(state, J, sunctx); - checkFlag(LS, "SUNLinSol_KLU", 0); - - retval = SUNLinSol_KLUSetOrdering(LS, 0); - checkFlag(&retval, "SUNLinSol_KLUSetOrdering", 1); - } -#else - J = SUNDenseMatrix(svsize, svsize, sunctx); - checkFlag(J, "SUNSparseMatrix", 0); - /* Create KLU solver object */ - LS = SUNLinSol_Dense(state, J, sunctx); - checkFlag(LS, "SUNLinSol_Dense", 0); #endif - retval = KINSetLinearSolver(solverMem, LS, J); + retval = KINSetLinearSolver(solverMem, LS, J); - checkFlag(&retval, "KINSetLinearSolver", 1); + checkFlag(&retval, "KINSetLinearSolver", 1); + } retval = KINSetJacFn(solverMem, kinsolJac); checkFlag(&retval, "KINSetJacFn", 1); diff --git a/src/griddyn/solvers/SundialsInterface.cpp b/src/griddyn/solvers/SundialsInterface.cpp index 7d1aa7371..a0eb7ec47 100644 --- a/src/griddyn/solvers/SundialsInterface.cpp +++ b/src/griddyn/solvers/SundialsInterface.cpp @@ -92,17 +92,10 @@ SundialsInterface::~SundialsInterface() if (types != nullptr) { NVECTOR_DESTROY(use_omp, types); } - if (flags[INITIALIZED_FLAG]) { - if (m_sundialsInfoFile != nullptr) { - static_cast(fclose(m_sundialsInfoFile)); - } - if (LS != nullptr) { - SUNLinSolFree(LS); - } - if (J != nullptr) { - SUNMatDestroy(J); - } + if (m_sundialsInfoFile != nullptr) { + static_cast(fclose(m_sundialsInfoFile)); } + freeLinearSolver(); if (sunctx != nullptr) { SUNContext_Free(&sunctx); } @@ -143,6 +136,7 @@ void SundialsInterface::allocate(count_t stateCount, count_t /*numRoots*/) [[maybe_unused]] bool prevOmp = use_omp; // looks unused if OPENMP is not available use_omp = flags[USE_OMP_FLAG]; flags.reset(INITIALIZED_FLAG); + freeLinearSolver(); if (state != nullptr) { NVECTOR_DESTROY(prevOmp, state); } @@ -245,6 +239,18 @@ void SundialsInterface::registerErrorHandler() checkFlag(&retval, "SUNContext_PushErrHandler", 1); } +void SundialsInterface::freeLinearSolver() +{ + if (LS != nullptr) { + SUNLinSolFree(LS); + LS = nullptr; + } + if (J != nullptr) { + SUNMatDestroy(J); + J = nullptr; + } +} + void SundialsInterface::kluReInit(SparseReinitMode sparseReInitModes) { #ifdef GRIDDYN_ENABLE_KLU diff --git a/src/griddyn/solvers/SundialsInterface.h b/src/griddyn/solvers/SundialsInterface.h index 19cd1ddaf..99224afea 100644 --- a/src/griddyn/solvers/SundialsInterface.h +++ b/src/griddyn/solvers/SundialsInterface.h @@ -128,6 +128,7 @@ class SundialsInterface: public SolverInterface { protected: void kluReInit(SparseReinitMode sparseReinitMode); void registerErrorHandler(); + void freeLinearSolver(); }; int sundialsJac(sunrealtype time, From 0c979bd341b83f8d31de01acf521fcf8752d151e Mon Sep 17 00:00:00 2001 From: Philip Top Date: Thu, 30 Jul 2026 06:48:46 -0700 Subject: [PATCH 2/3] clean up vector extraction in sundials --- src/griddyn/solvers/SundialsInterface.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/griddyn/solvers/SundialsInterface.cpp b/src/griddyn/solvers/SundialsInterface.cpp index a0eb7ec47..a92f3edcb 100644 --- a/src/griddyn/solvers/SundialsInterface.cpp +++ b/src/griddyn/solvers/SundialsInterface.cpp @@ -392,6 +392,8 @@ int sundialsJac(sunrealtype time, N_Vector /*tmp2*/) { auto sd = reinterpret_cast(userData); + auto* stateData = nvecdata(sd->use_omp, state); + auto* dstateData = nvecdata(sd->use_omp, dstateDt); if (matrixNeedsSetup(sd->jacCallCount, j)) { auto a1 = makeSparseMatrix(sd->svsize, sd->maxNNZ); @@ -403,8 +405,8 @@ int sundialsJac(sunrealtype time, MatrixDataFilter filterAd(*(a1)); filterAd.addFilter(sd->maskElements); sd->m_gds->jacobianFunction(time, - nvecdata(sd->use_omp, state), - nvecdata(sd->use_omp, dstateDt), + stateData, + dstateData, filterAd, cj, sd->mode); @@ -413,8 +415,8 @@ int sundialsJac(sunrealtype time, } } else { sd->m_gds->jacobianFunction(time, - nvecdata(sd->use_omp, state), - nvecdata(sd->use_omp, dstateDt), + stateData, + dstateData, *a1, cj, sd->mode); @@ -444,8 +446,8 @@ int sundialsJac(sunrealtype time, MatrixDataFilter filterAd(*a1); filterAd.addFilter(sd->maskElements); sd->m_gds->jacobianFunction(time, - nvecdata(sd->use_omp, state), - nvecdata(sd->use_omp, dstateDt), + stateData, + dstateData, filterAd, cj, sd->mode); @@ -454,8 +456,8 @@ int sundialsJac(sunrealtype time, } } else { sd->m_gds->jacobianFunction(time, - nvecdata(sd->use_omp, state), - nvecdata(sd->use_omp, dstateDt), + stateData, + dstateData, *a1, cj, sd->mode); From dca360e29fcedaedf1934487dc31d4c5b508dd66 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 30 Jul 2026 14:26:38 +0000 Subject: [PATCH 3/3] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/griddyn/solvers/SundialsInterface.cpp | 28 ++++------------------- 1 file changed, 4 insertions(+), 24 deletions(-) diff --git a/src/griddyn/solvers/SundialsInterface.cpp b/src/griddyn/solvers/SundialsInterface.cpp index a92f3edcb..0a05558c4 100644 --- a/src/griddyn/solvers/SundialsInterface.cpp +++ b/src/griddyn/solvers/SundialsInterface.cpp @@ -404,22 +404,12 @@ int sundialsJac(sunrealtype time, if (sd->flags[USE_MASK_FLAG]) { MatrixDataFilter filterAd(*(a1)); filterAd.addFilter(sd->maskElements); - sd->m_gds->jacobianFunction(time, - stateData, - dstateData, - filterAd, - cj, - sd->mode); + sd->m_gds->jacobianFunction(time, stateData, dstateData, filterAd, cj, sd->mode); for (auto& v : sd->maskElements) { a1->assign(v, v, 1.0); } } else { - sd->m_gds->jacobianFunction(time, - stateData, - dstateData, - *a1, - cj, - sd->mode); + sd->m_gds->jacobianFunction(time, stateData, dstateData, *a1, cj, sd->mode); } ++sd->jacCallCount; @@ -445,22 +435,12 @@ int sundialsJac(sunrealtype time, if (sd->flags[USE_MASK_FLAG]) { MatrixDataFilter filterAd(*a1); filterAd.addFilter(sd->maskElements); - sd->m_gds->jacobianFunction(time, - stateData, - dstateData, - filterAd, - cj, - sd->mode); + sd->m_gds->jacobianFunction(time, stateData, dstateData, filterAd, cj, sd->mode); for (auto& v : sd->maskElements) { a1->assign(v, v, 1.0); } } else { - sd->m_gds->jacobianFunction(time, - stateData, - dstateData, - *a1, - cj, - sd->mode); + sd->m_gds->jacobianFunction(time, stateData, dstateData, *a1, cj, sd->mode); } sd->jacCallCount++;