Minor suggestions to SUNDIALS interface

examples/PowerElectronics/RLCircuit/RLCircuit.cpp

@@ -14,15 +14,15 @@
 #include <Solver/Dynamic/DynamicSolver.hpp>
 #include <Solver/Dynamic/Ida.hpp>
 
-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;
   bool   use_jac = true;
 
   // TODO:setup as named parameters
   // Create circuit model
-  auto* sysmodel = new GridKit::PowerElectronicsModel<double, size_t>(rel_tol, abs_tol, use_jac);
+  GridKit::PowerElectronicsModel<double, size_t> 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<double, size_t>* idas = new AnalysisManager::Sundials::Ida<double, size_t>(sysmodel);
+  AnalysisManager::Sundials::Ida<double, size_t> 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<double>& yfinial = sysmodel->y();
+  std::vector<double>& yfinial = sysmodel.y();
 
   std::cout << "Final Vector y\n";
   for (size_t i = 0; i < yfinial.size(); i++)

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 <class ScalarT, typename IdxT>
     int Ida<ScalarT, IdxT>::runSimulation(real_type tf, int nout, const std::optional<std::function<void(real_type)>> step_callback)
     {
@@ -205,8 +222,8 @@ namespace AnalysisManager
       real_type dt   = (tf - t_init_) / static_cast<real_type>(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 <class ScalarT, typename IdxT>
+    int Ida<ScalarT, IdxT>::deleteBackwardSimulation()
+    {
+      N_VDestroy(yyB_);
+      N_VDestroy(ypB_);
+      N_VDestroy(qB_);
+      SUNLinSolFree(linearSolverB_);
+      SUNMatDestroy(JacobianMatB_);
+
+      return 0;
+    }
+
     template <class ScalarT, typename IdxT>
     int Ida<ScalarT, IdxT>::Residual(real_type tres, N_Vector yy, N_Vector yp, N_Vector rr, void* user_data)
     {

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<std::function<void(real_type)>> 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()
       {