Refactor SteadystateProblem, simplify solver creation (#2810)

Avoid unnecessary entanglement of forward/bwd solver setup. Less complex control flow.

Author: dweindl

include/amici/steadystateproblem.h

@@ -493,21 +493,6 @@ class SteadystateProblem {
      */
     void runSteadystateSimulationBwd(Solver const& solver, Model& model);
 
-    /**
-     * @brief Create and initialize a CVodeSolver instance for
-     * preequilibration simulation.
-     * @param solver Solver instance
-     * @param model Model instance.
-     * @param forwardSensis flag switching on integration with FSA
-     * @param backward flag switching on quadrature computation
-     * @param t0 Initial time for the steady state simulation.
-     * @return A unique pointer to the created Solver instance.
-     */
-    std::unique_ptr<Solver> createSteadystateSimSolver(
-        Solver const& solver, Model& model, bool forwardSensis, bool backward,
-        realtype t0
-    ) const;
-
     /**
      * @brief Initialize forward computation
      * @param it Index of the current output time point.

src/newton_solver.cpp

@@ -30,7 +30,8 @@ NewtonSolver::NewtonSolver(
             break;
         default:
             throw NewtonFailure(
-                AMICI_NOT_IMPLEMENTED, "Unknown linear solver type"
+                AMICI_NOT_IMPLEMENTED,
+                "Invalid solver for steady state simulation"
             );
         }
     } catch (NewtonFailure const&) {

src/steadystateproblem.cpp

@@ -348,18 +348,27 @@ SteadyStateStatus SteadystateProblem::findSteadyStateBySimulation(
 ) {
     try {
         if (it < 0) {
-            // Preequilibration? -> Create a new solver instance for simulation
-            bool integrateSensis = requires_state_sensitivities(
-                model, solver, it, SteadyStateContext::solverCreation
-            );
-            auto newtonSimSolver = createSteadystateSimSolver(
-                solver, model, integrateSensis, false, t0
-            );
-            runSteadystateSimulationFwd(*newtonSimSolver, model);
+            // Preequilibration -> Create a new solver instance for simulation
+            auto sim_solver = std::unique_ptr<Solver>(solver.clone());
+            sim_solver->logger = solver.logger;
+
+            // do we need sensitivities?
+            if (requires_state_sensitivities(
+                    model, solver, it, SteadyStateContext::solverCreation
+                )) {
+                // need forward to compute sx0
+                sim_solver->setSensitivityMethod(SensitivityMethod::forward);
+            } else {
+                sim_solver->setSensitivityMethod(SensitivityMethod::none);
+                sim_solver->setSensitivityOrder(SensitivityOrder::none);
+            }
+            sim_solver->setup(t0, &model, state_.x, state_.dx, state_.sx, sdx_);
+            runSteadystateSimulationFwd(*sim_solver, model);
         } else {
-            // Solver was already created, use this one
+            // Postequilibration -> Solver was already created, use that one
             runSteadystateSimulationFwd(solver, model);
         }
+
         return SteadyStateStatus::success;
     } catch (IntegrationFailure const& ex) {
         switch (ex.error_code) {
@@ -498,11 +507,18 @@ void SteadystateProblem::getQuadratureBySimulation(
     // xQ was written in getQuadratureByLinSolve() -> set to zero
     xQ_.zero();
 
-    auto simSolver = createSteadystateSimSolver(solver, model, false, true, t0);
+    auto sim_solver = std::unique_ptr<Solver>(solver.clone());
+    sim_solver->logger = solver.logger;
+    sim_solver->setSensitivityMethod(SensitivityMethod::none);
+    sim_solver->setSensitivityOrder(SensitivityOrder::none);
+    sim_solver->setup(t0, &model, xB_, xB_, state_.sx, sdx_);
+    sim_solver->setupSteadystate(
+        t0, &model, state_.x, state_.dx, xB_, xB_, xQ_
+    );
 
     // perform integration and quadrature
     try {
-        runSteadystateSimulationBwd(*simSolver, model);
+        runSteadystateSimulationBwd(*sim_solver, model);
         hasQuadrature_ = true;
     } catch (NewtonFailure const&) {
         hasQuadrature_ = false;
@@ -779,44 +795,6 @@ void SteadystateProblem::runSteadystateSimulationBwd(
     }
 }
 
-std::unique_ptr<Solver> SteadystateProblem::createSteadystateSimSolver(
-    Solver const& solver, Model& model, bool forwardSensis, bool backward,
-    realtype t0
-) const {
-    switch (solver.getLinearSolver()) {
-    case LinearSolver::dense:
-    case LinearSolver::KLU:
-        break;
-    default:
-        throw AmiException("Invalid solver for steady state simulation");
-    }
-
-    auto sim_solver = std::unique_ptr<Solver>(solver.clone());
-
-    sim_solver->logger = solver.logger;
-
-    // do we need sensitivities?
-    if (forwardSensis) {
-        // need forward to compute sx0
-        sim_solver->setSensitivityMethod(SensitivityMethod::forward);
-    } else {
-        sim_solver->setSensitivityMethod(SensitivityMethod::none);
-        sim_solver->setSensitivityOrder(SensitivityOrder::none);
-    }
-    // use x and sx as dummies for dx and sdx
-    // (they won't get touched in a CVodeSolver)
-    if (backward) {
-        sim_solver->setup(t0, &model, xB_, xB_, state_.sx, sdx_);
-        sim_solver->setupSteadystate(
-            t0, &model, state_.x, state_.dx, xB_, xB_, xQ_
-        );
-    } else {
-        sim_solver->setup(t0, &model, state_.x, state_.dx, state_.sx, sdx_);
-    }
-
-    return sim_solver;
-}
-
 void SteadystateProblem::getAdjointUpdates(
     Model& model, ExpData const& edata, std::vector<realtype>& dJydx
 ) {