From 0fe82174c8eb3f46f9b3ab668c8ef60264a249eb Mon Sep 17 00:00:00 2001
From: Thomas Fischer <thomas.fischer@ufz.de>
Date: Thu, 4 Mar 2021 14:26:03 +0100
Subject: [PATCH] [NL/ODESolver/Newton] Store ids to reuse the xdot vectors.

---
 NumLib/ODESolver/NonlinearSolver.cpp          | 4 ++--
 NumLib/ODESolver/NonlinearSolver.h            | 4 ++++
 NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 3 +--
 3 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index fc1caf567ac..6a8ffcd420a 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -41,7 +41,7 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::
     _equation_system->getRhs(*x_prev[process_id], rhs);
 
     // r_neq = A * x - rhs
-    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
     MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
     MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs);  // res -= rhs
 }
@@ -220,7 +220,7 @@ void NonlinearSolver<NonlinearSolverTag::Newton>::
     }
 
     _equation_system->assemble(x, x_prev, process_id);
-    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
+    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
     _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
 }
 
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 73ad95157e8..183368e2fab 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -141,6 +141,8 @@ private:
     std::size_t _minus_delta_x_id = 0u;  //!< ID of the \f$ -\Delta x\f$ vector.
     std::size_t _x_new_id =
         0u;  //!< ID of the vector storing \f$ x - (-\Delta x) \f$.
+    std::size_t _r_neq_id = 0u;  //!< ID of the non-equilibrium initial
+                                 //! residuum vector.
 
     /// Enables computation of the non-equilibrium initial residuum \f$ r_{\rm
     /// neq} \f$ before the first time step. The forces are zero if the external
@@ -214,6 +216,8 @@ private:
     std::size_t _rhs_id = 0u;    //!< ID of the right-hand side vector.
     std::size_t _x_new_id = 0u;  //!< ID of the vector storing the solution of
                                  //! the linearized equation.
+    std::size_t _r_neq_id = 0u;  //!< ID of the non-equilibrium initial
+                                 //! residuum vector.
 
     // clang-format off
     /// \copydoc NumLib::NonlinearSolver<NonlinearSolverTag::Newton>::_compensate_non_equilibrium_initial_residuum
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index c17197fbb51..d28ea348a1c 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -87,8 +87,7 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
     _xdot_ids.resize(x_new_timestep.size());
     for (std::size_t i = 0; i < xdot.size(); i++)
     {
-        xdot[i] =
-            &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
+        xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
         _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
     }
 
-- 
GitLab