diff --git a/NumLib/DOF/MatrixProviderUser.h b/NumLib/DOF/MatrixProviderUser.h
index cf3a53a2ef417cb938139c851872467f7d9e2243..440e3692facc1a7240fa4681951c8a48767c5f24 100644
--- a/NumLib/DOF/MatrixProviderUser.h
+++ b/NumLib/DOF/MatrixProviderUser.h
@@ -54,9 +54,6 @@ public:
class VectorProvider
{
public:
- //! Get an uninitialized vector.
- virtual GlobalVector& getVector() = 0;
-
//! Get an uninitialized vector with the given \c id.
virtual GlobalVector& getVector(std::size_t& id) = 0;
diff --git a/NumLib/DOF/SimpleMatrixVectorProvider.cpp b/NumLib/DOF/SimpleMatrixVectorProvider.cpp
index c1d27d5c1d0d49f96340a4a6eed6782af76372a8..67fe93e71545deeaa79bc6d978f9f34e984d8e5b 100644
--- a/NumLib/DOF/SimpleMatrixVectorProvider.cpp
+++ b/NumLib/DOF/SimpleMatrixVectorProvider.cpp
@@ -176,15 +176,6 @@ getVector_(std::size_t& id, Args&&... args)
return get_<do_search>(id, _unused_vectors, _used_vectors, std::forward<Args>(args)...);
}
-
-GlobalVector&
-SimpleMatrixVectorProvider::
-getVector()
-{
- std::size_t id = 0u;
- return *getVector_<false>(id).first;
-}
-
GlobalVector&
SimpleMatrixVectorProvider::
getVector(std::size_t& id)
diff --git a/NumLib/DOF/SimpleMatrixVectorProvider.h b/NumLib/DOF/SimpleMatrixVectorProvider.h
index 3498a52d08a4f1e7bd484c8c31803c1345aa5032..328333bae61a055656aafc0b5002cb81ba1a5e1d 100644
--- a/NumLib/DOF/SimpleMatrixVectorProvider.h
+++ b/NumLib/DOF/SimpleMatrixVectorProvider.h
@@ -36,7 +36,6 @@ public:
SimpleMatrixVectorProvider(SimpleMatrixVectorProvider const&) = delete;
SimpleMatrixVectorProvider& operator=(SimpleMatrixVectorProvider const&) = delete;
- GlobalVector& getVector() override;
GlobalVector& getVector(std::size_t& id) override;
GlobalVector& getVector(GlobalVector const& x) override;
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index fc1caf567ac7c4315145b182c734235be2e913c8..6a8ffcd420afcdd6813fc01fbe47eb3e62127571 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 73ad95157e83990ba41d48454524f99a973c09d4..183368e2fab9625256ab12f3d5802f4179a1c9a9 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 52c135787d63180d51221a240619721a65239edf..d28ea348a1c45f9a2e17e8b815ddb70edc619b3d 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -84,9 +84,10 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
auto const dxdot_dx = _time_disc.getNewXWeight();
std::vector<GlobalVector*> xdot(x_new_timestep.size());
+ _xdot_ids.resize(x_new_timestep.size());
for (std::size_t i = 0; i < xdot.size(); i++)
{
- xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector();
+ xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
_time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
}
@@ -220,9 +221,11 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
auto const dt = _time_disc.getCurrentTimeIncrement();
auto const& x_curr = *x_new_timestep[process_id];
std::vector<GlobalVector*> xdot(x_new_timestep.size());
+ _xdot_ids.resize(x_new_timestep.size());
+
for (std::size_t i = 0; i < xdot.size(); i++)
{
- xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector();
+ xdot[i] = &NumLib::GlobalVectorProvider::provider.getVector(_xdot_ids[i]);
_time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 0264082d0e1be4f441a711b0512e5004a6c0938a..a2b9ffe1f7293208cd04fab3df25180a8ff873bd 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -148,6 +148,7 @@ private:
//! ID of the vector storing xdot in intermediate computations.
mutable std::size_t _xdot_id = 0u;
+ mutable std::vector<std::size_t> _xdot_ids;
};
/*! Time discretized first order implicit quasi-linear ODE;
@@ -239,6 +240,7 @@ private:
std::size_t _M_id = 0u; //!< ID of the \c _M matrix.
std::size_t _K_id = 0u; //!< ID of the \c _K matrix.
std::size_t _b_id = 0u; //!< ID of the \c _b vector.
+ mutable std::vector<std::size_t> _xdot_ids;
};
//! @}
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index f501ea6e38f1ed95e47f992f12931bb3301f3c2e..eab386ff5381c4864a7ef5bc993bd04cf3388e5d 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -104,10 +104,14 @@ void postTimestepForAllProcesses(
double const t, double const dt,
std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
std::vector<GlobalVector*> const& process_solutions,
- std::vector<GlobalVector*> const& process_solutions_prev)
+ std::vector<GlobalVector*> const& process_solutions_prev,
+ std::vector<std::size_t>& xdot_vector_ids)
{
std::vector<GlobalVector*> x_dots;
x_dots.reserve(per_process_data.size());
+ xdot_vector_ids.resize(per_process_data.size());
+
+ std::size_t cnt = 0;
for (auto& process_data : per_process_data)
{
auto const process_id = process_data->process_id;
@@ -115,7 +119,8 @@ void postTimestepForAllProcesses(
auto const& time_discretization = *process_data->time_disc;
x_dots.emplace_back(&NumLib::GlobalVectorProvider::provider.getVector(
- ode_sys.getMatrixSpecifications(process_id)));
+ ode_sys.getMatrixSpecifications(process_id), xdot_vector_ids[cnt]));
+ cnt++;
time_discretization.getXdot(*process_solutions[process_id],
*process_solutions_prev[process_id],
@@ -274,7 +279,8 @@ void calculateNonEquilibriumInitialResiduum(
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
std::vector<GlobalVector*>& x, std::vector<GlobalVector*> const& x_prev,
std::size_t const timestep, double const t, double const delta_t,
- ProcessData const& process_data, Output& output_control)
+ ProcessData const& process_data, Output& output_control,
+ std::size_t& xdot_id)
{
auto& process = process_data.process;
int const process_id = process_data.process_id;
@@ -308,7 +314,7 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
}
GlobalVector& x_dot = NumLib::GlobalVectorProvider::provider.getVector(
- ode_sys.getMatrixSpecifications(process_id));
+ ode_sys.getMatrixSpecifications(process_id), xdot_id);
time_disc.getXdot(*x[process_id], *x_prev[process_id], x_dot);
@@ -622,9 +628,9 @@ bool TimeLoop::loop()
// iteration, an exception thrown in assembly, for example.
if (nonlinear_solver_status.error_norms_met)
{
- postTimestepForAllProcesses(t, dt, _per_process_data,
- _process_solutions,
- _process_solutions_prev);
+ postTimestepForAllProcesses(
+ t, dt, _per_process_data, _process_solutions,
+ _process_solutions_prev, _xdot_vector_ids);
}
INFO("[time] Time step #{:d} took {:g} s.", timesteps,
@@ -677,13 +683,14 @@ bool TimeLoop::loop()
static NumLib::NonlinearSolverStatus solveMonolithicProcess(
const double t, const double dt, const std::size_t timestep_id,
ProcessData const& process_data, std::vector<GlobalVector*>& x,
- std::vector<GlobalVector*> const& x_prev, Output& output)
+ std::vector<GlobalVector*> const& x_prev, Output& output,
+ std::size_t& xdot_id)
{
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
- x, x_prev, timestep_id, t, dt, process_data, output);
+ x, x_prev, timestep_id, t, dt, process_data, output, xdot_id);
INFO("[time] Solving process #{:d} took {:g} s in time step #{:d} ",
process_data.process_id, time_timestep_process.elapsed(), timestep_id);
@@ -703,7 +710,7 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
auto const process_id = process_data->process_id;
nonlinear_solver_status = solveMonolithicProcess(
t, dt, timestep_id, *process_data, _process_solutions,
- _process_solutions_prev, *_output);
+ _process_solutions_prev, *_output, _xdot_id);
process_data->nonlinear_solver_status = nonlinear_solver_status;
if (!nonlinear_solver_status.error_norms_met)
@@ -776,7 +783,7 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
nonlinear_solver_status = solveOneTimeStepOneProcess(
_process_solutions, _process_solutions_prev, timestep_id, t, dt,
- *process_data, *_output);
+ *process_data, *_output, _xdot_id);
process_data->nonlinear_solver_status = nonlinear_solver_status;
INFO(
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index bc52a8c1a34b3d56730210cda6934421356e1589..0992c4c1a8bd750b534143616957e49f599246aa 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -126,5 +126,14 @@ private:
/// Solutions of the previous coupling iteration for the convergence
/// criteria of the coupling iteration.
std::vector<GlobalVector*> _solutions_of_last_cpl_iteration;
+
+ /// store the ids of the global xdot vectors in the global vector
+ /// provider for reuse, needed in postTimestepForAllProcesses
+ /// the length of the vector is the size of _per_process_data
+ std::vector<std::size_t> _xdot_vector_ids;
+
+ // store the id of the global xdot vector in the global vector provider for
+ // reuse; needed in solveOneTimeStepOneProcess
+ std::size_t _xdot_id = 0;
};
} // namespace ProcessLib