Double initialize

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);

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);

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);

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);

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<void>(fclose(m_sundialsInfoFile));
-        }
-        if (LS != nullptr) {
-            SUNLinSolFree(LS);
-        }
-        if (J != nullptr) {
-            SUNMatDestroy(J);
-        }
+    if (m_sundialsInfoFile != nullptr) {
+        static_cast<void>(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
@@ -386,6 +392,8 @@ int sundialsJac(sunrealtype time,
                 N_Vector /*tmp2*/)
 {
     auto sd = reinterpret_cast<SundialsInterface*>(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);
@@ -396,22 +404,12 @@ int sundialsJac(sunrealtype time,
         if (sd->flags[USE_MASK_FLAG]) {
             MatrixDataFilter<double> filterAd(*(a1));
             filterAd.addFilter(sd->maskElements);
-            sd->m_gds->jacobianFunction(time,
-                                        nvecdata(sd->use_omp, state),
-                                        nvecdata(sd->use_omp, dstateDt),
-                                        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,
-                                        nvecdata(sd->use_omp, state),
-                                        nvecdata(sd->use_omp, dstateDt),
-                                        *a1,
-                                        cj,
-                                        sd->mode);
+            sd->m_gds->jacobianFunction(time, stateData, dstateData, *a1, cj, sd->mode);
         }
 
         ++sd->jacCallCount;
@@ -437,22 +435,12 @@ int sundialsJac(sunrealtype time,
         if (sd->flags[USE_MASK_FLAG]) {
             MatrixDataFilter<double> filterAd(*a1);
             filterAd.addFilter(sd->maskElements);
-            sd->m_gds->jacobianFunction(time,
-                                        nvecdata(sd->use_omp, state),
-                                        nvecdata(sd->use_omp, dstateDt),
-                                        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,
-                                        nvecdata(sd->use_omp, state),
-                                        nvecdata(sd->use_omp, dstateDt),
-                                        *a1,
-                                        cj,
-                                        sd->mode);
+            sd->m_gds->jacobianFunction(time, stateData, dstateData, *a1, cj, sd->mode);
         }
 
         sd->jacCallCount++;

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,