From b3fa2f07f0caeff8db1b692ad881b3190a09c8bc Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 27 Apr 2020 17:14:53 +0200
Subject: [PATCH] Store previous time step solution in TimeLoop.
This simplifies the TimeDiscretization interface:
pushState and popState are no longer needed. The corresponding
copies are done in the time loop.
Since the previous time step solution is no longer stored
in the time discretization, it is passed to all required
functions like assemble, getRHS and the like.
---
NumLib/ODESolver/MatrixTranslator.cpp | 5 +-
NumLib/ODESolver/MatrixTranslator.h | 8 +-
NumLib/ODESolver/NonlinearSolver.cpp | 28 ++++---
NumLib/ODESolver/NonlinearSolver.h | 15 +++-
NumLib/ODESolver/NonlinearSystem.h | 6 +-
NumLib/ODESolver/TimeDiscretization.cpp | 6 +-
NumLib/ODESolver/TimeDiscretization.h | 74 +++++--------------
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 11 ++-
NumLib/ODESolver/TimeDiscretizedODESystem.h | 7 +-
ProcessLib/TimeLoop.cpp | 66 ++++++++++-------
ProcessLib/TimeLoop.h | 1 +
Tests/NumLib/TimeLoopSingleODE.h | 11 ++-
12 files changed, 123 insertions(+), 115 deletions(-)
diff --git a/NumLib/ODESolver/MatrixTranslator.cpp b/NumLib/ODESolver/MatrixTranslator.cpp
index e6cf7ee6dda..ff7d31b25cb 100644
--- a/NumLib/ODESolver/MatrixTranslator.cpp
+++ b/NumLib/ODESolver/MatrixTranslator.cpp
@@ -29,12 +29,13 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
- const GlobalVector& b, GlobalVector& rhs) const
+ const GlobalVector& b, const GlobalVector& x_prev,
+ GlobalVector& rhs) const
{
namespace LinAlg = MathLib::LinAlg;
auto& tmp = NumLib::GlobalVectorProvider::provider.getVector(_tmp_id);
- _time_disc.getWeightedOldX(tmp);
+ _time_disc.getWeightedOldX(tmp, x_prev);
// rhs = M * weighted_old_x + b
LinAlg::matMultAdd(M, tmp, b, rhs);
diff --git a/NumLib/ODESolver/MatrixTranslator.h b/NumLib/ODESolver/MatrixTranslator.h
index 2008c173bd5..33a4ba00098 100644
--- a/NumLib/ODESolver/MatrixTranslator.h
+++ b/NumLib/ODESolver/MatrixTranslator.h
@@ -44,9 +44,10 @@ public:
virtual void computeA(GlobalMatrix const& M, GlobalMatrix const& K,
GlobalMatrix& A) const = 0;
- //! Computes \c rhs from \c M, \c K and \c b.
+ //! Computes \c rhs from \c M, \c K, \c b and \c x_prev.
virtual void computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
- const GlobalVector& b, GlobalVector& rhs) const = 0;
+ const GlobalVector& b, const GlobalVector& x_prev,
+ GlobalVector& rhs) const = 0;
/*! Computes \c res from \c M, \c K, \c b, \f$ \hat x \f$ and \f$ x_N \f$.
* You might also want read the remarks on
@@ -101,7 +102,8 @@ public:
//! Computes \f$ \mathtt{rhs} = M \cdot x_O + b \f$.
void computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
- const GlobalVector& b, GlobalVector& rhs) const override;
+ const GlobalVector& b, const GlobalVector& x_prev,
+ GlobalVector& rhs) const override;
//! Computes \f$ r = M \cdot \hat x + K \cdot x_C - b \f$.
void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 8d7127a480f..cd0195360de 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -23,8 +23,9 @@
namespace NumLib
{
void NonlinearSolver<NonlinearSolverTag::Picard>::
- calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
- int const process_id)
+ calculateNonEquilibriumInitialResiduum(
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id)
{
if (!_compensate_non_equilibrium_initial_residuum)
{
@@ -33,9 +34,9 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::
auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
- _equation_system->assemble(x, process_id);
+ _equation_system->assemble(x, x_prev, process_id);
_equation_system->getA(A);
- _equation_system->getRhs(rhs);
+ _equation_system->getRhs(*x_prev[process_id], rhs);
// r_neq = A * x - rhs
_r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
@@ -45,6 +46,7 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
std::function<void(int, std::vector<GlobalVector*> const&)> const&
postIterationCallback,
int const process_id)
@@ -85,9 +87,9 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
BaseLib::RunTime time_assembly;
time_assembly.start();
- sys.assemble(x_new, process_id);
+ sys.assemble(x_new, x_prev, process_id);
sys.getA(A);
- sys.getRhs(rhs);
+ sys.getRhs(*x_prev[process_id], rhs);
INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
// Subract non-equilibrium initial residuum if set
@@ -206,21 +208,23 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
}
void NonlinearSolver<NonlinearSolverTag::Newton>::
- calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
- int const process_id)
+ calculateNonEquilibriumInitialResiduum(
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id)
{
if (!_compensate_non_equilibrium_initial_residuum)
{
return;
}
- _equation_system->assemble(x, process_id);
+ _equation_system->assemble(x, x_prev, process_id);
_r_neq = &NumLib::GlobalVectorProvider::provider.getVector();
- _equation_system->getResidual(*x[process_id], *_r_neq);
+ _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
}
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
std::function<void(int, std::vector<GlobalVector*> const&)> const&
postIterationCallback,
int const process_id)
@@ -265,7 +269,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
time_assembly.start();
try
{
- sys.assemble(x, process_id);
+ sys.assemble(x, x_prev, process_id);
}
catch (AssemblyException const& e)
{
@@ -275,7 +279,7 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
iteration = _maxiter;
break;
}
- sys.getResidual(*x[process_id], res);
+ sys.getResidual(*x[process_id], *x_prev[process_id], res);
sys.getJacobian(J);
INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 6b250966857..26f2c94d228 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -36,11 +36,13 @@ class NonlinearSolverBase
{
public:
virtual void calculateNonEquilibriumInitialResiduum(
- std::vector<GlobalVector*> const& x, int const process_id) = 0;
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev, int const process_id) = 0;
/*! Assemble and solve the equation system.
*
* \param x in: the initial guess, out: the solution.
+ * \param x_prev previous time step solution.
* \param postIterationCallback called after each iteration if set.
* \param process_id usually used in staggered schemes.
*
@@ -49,6 +51,7 @@ public:
*/
virtual NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
std::function<void(int, std::vector<GlobalVector*> const&)> const&
postIterationCallback,
int const process_id) = 0;
@@ -100,10 +103,13 @@ public:
}
void calculateNonEquilibriumInitialResiduum(
- std::vector<GlobalVector*> const& x, int const process_id) override;
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
+ int const process_id) override;
NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
std::function<void(int, std::vector<GlobalVector*> const&)> const&
postIterationCallback,
int const process_id) override;
@@ -175,10 +181,13 @@ public:
}
void calculateNonEquilibriumInitialResiduum(
- std::vector<GlobalVector*> const& x, int const process_id) override;
+ std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
+ int const process_id) override;
NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
+ std::vector<GlobalVector*> const& x_prev,
std::function<void(int, std::vector<GlobalVector*> const&)> const&
postIterationCallback,
int const process_id) override;
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index 85c3787af2f..6094520a50c 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -38,6 +38,7 @@ public:
//! The linearized system is \f$A(x) \cdot x = b(x)\f$. Here the matrix
//! \f$A(x)\f$ and the vector \f$b(x)\f$ are assembled.
virtual void assemble(std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id) = 0;
/*! Writes the residual at point \c x to \c res.
@@ -47,6 +48,7 @@ public:
* \todo Remove argument \c x.
*/
virtual void getResidual(GlobalVector const& x,
+ GlobalVector const& x_prev,
GlobalVector& res) const = 0;
/*! Writes the Jacobian of the residual to \c Jac.
@@ -85,6 +87,7 @@ public:
//! The linearized system is \f$J(x) \cdot \Delta x = (x)\f$. Here the
//! residual vector \f$r(x)\f$ and its Jacobian \f$J(x)\f$ are assembled.
virtual void assemble(std::vector<GlobalVector*> const& x,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id) = 0;
//! Writes the linearized equation system matrix to \c A.
@@ -93,7 +96,8 @@ public:
//! Writes the linearized equation system right-hand side to \c rhs.
//! \pre assemble() must have been called before.
- virtual void getRhs(GlobalVector& rhs) const = 0;
+ virtual void getRhs(GlobalVector const& x_prev,
+ GlobalVector& rhs) const = 0;
//! Pre-compute known solutions and possibly store them internally.
virtual void computeKnownSolutions(GlobalVector const& x,
diff --git a/NumLib/ODESolver/TimeDiscretization.cpp b/NumLib/ODESolver/TimeDiscretization.cpp
index 08f9c684a6d..551384e4b4f 100644
--- a/NumLib/ODESolver/TimeDiscretization.cpp
+++ b/NumLib/ODESolver/TimeDiscretization.cpp
@@ -54,9 +54,11 @@ double TimeDiscretization::computeRelativeChangeFromPreviousTimestep(
}
double BackwardEuler::getRelativeChangeFromPreviousTimestep(
- GlobalVector const& x, MathLib::VecNormType norm_type)
+ GlobalVector const& x,
+ GlobalVector const& x_old,
+ MathLib::VecNormType norm_type)
{
- return computeRelativeChangeFromPreviousTimestep(x, _x_old, norm_type);
+ return computeRelativeChangeFromPreviousTimestep(x, x_old, norm_type);
}
} // end of namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index f53765d4f6c..f5dcbab099f 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -90,7 +90,7 @@ public:
TimeDiscretization() = default;
//! Sets the initial condition.
- virtual void setInitialState(const double t0, GlobalVector const& x0) = 0;
+ virtual void setInitialState(const double t0) = 0;
/*! Get the relative change of solutions between two successive time steps
* by \f$ e_n = \|u^{n+1}-u^{n}\|/\|u^{n+1}\| \f$.
@@ -99,30 +99,9 @@ public:
* \param norm_type The type of global vector norm.
*/
virtual double getRelativeChangeFromPreviousTimestep(
- GlobalVector const& x, MathLib::VecNormType norm_type) = 0;
-
- /*! Indicate that the current timestep is done and that you will proceed to
- * the next one.
- *
- * \warning Do not use this method for setting the initial condition,
- * rather use setInitialState()!
- *
- * \param t The current timestep.
- * \param x The solution at the current timestep.
- */
- virtual void pushState(const double t, GlobalVector const& x) = 0;
-
- /*!
- * Restores the given vector x to its old value.
- * Used only for repeating of the time step with new time step size when
- * the current time step is rejected. The restored x is only used as
- * an initial guess for linear solver in the first Picard nonlinear
- * iteration.
- *
- * \param x The solution at the current time step, which is going to be
- * reset to its previous value.
- */
- virtual void popState(GlobalVector& x) = 0;
+ GlobalVector const& x,
+ GlobalVector const& x_old,
+ MathLib::VecNormType norm_type) = 0;
/*! Indicate that the computation of a new timestep is being started now.
*
@@ -142,14 +121,16 @@ public:
//! Returns \f$ \hat x \f$, i.e. the discretized approximation of \f$ \dot x
//! \f$.
- void getXdot(GlobalVector const& x_at_new_timestep, GlobalVector& xdot) const
+ void getXdot(GlobalVector const& x_at_new_timestep,
+ GlobalVector const& x_old,
+ GlobalVector& xdot) const
{
namespace LinAlg = MathLib::LinAlg;
auto const dxdot_dx = getNewXWeight();
// xdot = dxdot_dx * x_at_new_timestep - x_old
- getWeightedOldX(xdot);
+ getWeightedOldX(xdot, x_old);
LinAlg::axpby(xdot, dxdot_dx, -1.0, x_at_new_timestep);
}
@@ -157,7 +138,8 @@ public:
virtual double getNewXWeight() const = 0;
//! Returns \f$ x_O \f$.
- virtual void getWeightedOldX(GlobalVector& y) const = 0; // = x_old
+ virtual void getWeightedOldX(
+ GlobalVector& y, GlobalVector const& x_old) const = 0; // = x_old
virtual ~TimeDiscretization() = default;
@@ -186,34 +168,12 @@ protected:
class BackwardEuler final : public TimeDiscretization
{
public:
- BackwardEuler()
- : _x_old(NumLib::GlobalVectorProvider::provider.getVector())
- {
- }
-
- ~BackwardEuler() override
- {
- NumLib::GlobalVectorProvider::provider.releaseVector(_x_old);
- }
-
- void setInitialState(const double t0, GlobalVector const& x0) override
- {
- _t = t0;
- MathLib::LinAlg::copy(x0, _x_old);
- }
+ void setInitialState(const double t0) override { _t = t0; }
double getRelativeChangeFromPreviousTimestep(
- GlobalVector const& x, MathLib::VecNormType norm_type) override;
-
- void pushState(const double /*t*/, GlobalVector const& x) override
- {
- MathLib::LinAlg::copy(x, _x_old);
- }
-
- void popState(GlobalVector& x) override
- {
- MathLib::LinAlg::copy(_x_old, x);
- }
+ GlobalVector const& x,
+ GlobalVector const& x_old,
+ MathLib::VecNormType norm_type) override;
void nextTimestep(const double t, const double delta_t) override
{
@@ -224,12 +184,13 @@ public:
double getCurrentTime() const override { return _t; }
double getCurrentTimeIncrement() const override { return _delta_t; }
double getNewXWeight() const override { return 1.0 / _delta_t; }
- void getWeightedOldX(GlobalVector& y) const override
+ void getWeightedOldX(GlobalVector& y,
+ GlobalVector const& x_old) const override
{
namespace LinAlg = MathLib::LinAlg;
// y = x_old / delta_t
- LinAlg::copy(_x_old, y);
+ LinAlg::copy(x_old, y);
LinAlg::scale(y, 1.0 / _delta_t);
}
@@ -237,7 +198,6 @@ private:
double _t = std::numeric_limits<double>::quiet_NaN(); //!< \f$ t_C \f$
double _delta_t =
std::numeric_limits<double>::quiet_NaN(); //!< the timestep size
- GlobalVector& _x_old; //!< the solution from the preceding timestep
};
//! @}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 85bd3bd4e08..21e60327490 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -73,6 +73,7 @@ TimeDiscretizedODESystem<
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id)
{
namespace LinAlg = MathLib::LinAlg;
@@ -86,7 +87,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
for (auto& v : xdot)
{
v = &NumLib::GlobalVectorProvider::provider.getVector();
- _time_disc.getXdot(*x_new_timestep[process_id], *v);
+ _time_disc.getXdot(*x_new_timestep[process_id], *x_prev[process_id],
+ *v);
}
_M->setZero();
@@ -124,13 +126,14 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
void TimeDiscretizedODESystem<
ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::getResidual(GlobalVector const& x_new_timestep,
+ GlobalVector const& x_prev,
GlobalVector& res) const
{
// TODO Maybe the duplicate calculation of xdot here and in assembleJacobian
// can be optimuized. However, that would make the interface a bit more
// fragile.
auto& xdot = NumLib::GlobalVectorProvider::provider.getVector(_xdot_id);
- _time_disc.getXdot(x_new_timestep, xdot);
+ _time_disc.getXdot(x_new_timestep, x_prev, xdot);
_mat_trans->computeResidual(*_M, *_K, *_b, x_new_timestep, xdot, res);
@@ -210,6 +213,7 @@ TimeDiscretizedODESystem<
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Picard>::
assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id)
{
namespace LinAlg = MathLib::LinAlg;
@@ -221,7 +225,8 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
for (auto& v : xdot)
{
v = &NumLib::GlobalVectorProvider::provider.getVector();
- _time_disc.getXdot(*x_new_timestep[process_id], *v);
+ _time_disc.getXdot(*x_new_timestep[process_id], *x_prev[process_id],
+ *v);
}
_M->setZero();
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index e28470e0f32..fbbe1a84a46 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -83,9 +83,11 @@ public:
~TimeDiscretizedODESystem() override;
void assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id) override;
void getResidual(GlobalVector const& x_new_timestep,
+ GlobalVector const& x_prev,
GlobalVector& res) const override;
void getJacobian(GlobalMatrix& Jac) const override;
@@ -173,6 +175,7 @@ public:
~TimeDiscretizedODESystem() override;
void assemble(std::vector<GlobalVector*> const& x_new_timestep,
+ std::vector<GlobalVector*> const& x_prev,
int const process_id) override;
void getA(GlobalMatrix& A) const override
@@ -180,9 +183,9 @@ public:
_mat_trans->computeA(*_M, *_K, A);
}
- void getRhs(GlobalVector& rhs) const override
+ void getRhs(GlobalVector const& x_prev, GlobalVector& rhs) const override
{
- _mat_trans->computeRhs(*_M, *_K, *_b, rhs);
+ _mat_trans->computeRhs(*_M, *_K, *_b, x_prev, rhs);
}
void computeKnownSolutions(GlobalVector const& x,
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index e7f43b569ed..53d3da26793 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -121,11 +121,13 @@ void setTimeDiscretizedODESystem(ProcessData& process_data)
setTimeDiscretizedODESystem(process_data, process_data.process);
}
-std::vector<GlobalVector*> setInitialConditions(
+std::pair<std::vector<GlobalVector*>, std::vector<GlobalVector*>>
+setInitialConditions(
double const t0,
std::vector<std::unique_ptr<ProcessData>> const& per_process_data)
{
std::vector<GlobalVector*> process_solutions;
+ std::vector<GlobalVector*> process_solutions_prev;
for (auto& process_data : per_process_data)
{
@@ -139,21 +141,27 @@ std::vector<GlobalVector*> setInitialConditions(
process_solutions.emplace_back(
&NumLib::GlobalVectorProvider::provider.getVector(
ode_sys.getMatrixSpecifications(process_id)));
+ process_solutions_prev.emplace_back(
+ &NumLib::GlobalVectorProvider::provider.getVector(
+ ode_sys.getMatrixSpecifications(process_id)));
- auto& x0 = *process_solutions[process_id];
- pcs.setInitialConditions(process_id, t0, x0);
- MathLib::LinAlg::finalizeAssembly(x0);
+ auto& x = *process_solutions[process_id];
+ auto& x_prev = *process_solutions_prev[process_id];
+ pcs.setInitialConditions(process_id, t0, x);
+ MathLib::LinAlg::finalizeAssembly(x);
- time_disc.setInitialState(t0, x0); // push IC
+ time_disc.setInitialState(t0); // push IC
+ MathLib::LinAlg::copy(x, x_prev); // pushState
}
- return process_solutions;
+ return {process_solutions, process_solutions_prev};
}
void calculateNonEquilibriumInitialResiduum(
std::vector<std::unique_ptr<ProcessData>> const& per_process_data,
std::vector<GlobalVector*>
- process_solutions)
+ process_solutions,
+ std::vector<GlobalVector*> const& process_solutions_prev)
{
INFO("Calculate non-equilibrium initial residuum.");
for (auto& process_data : per_process_data)
@@ -172,14 +180,15 @@ void calculateNonEquilibriumInitialResiduum(
double const dt = 1;
time_disc.nextTimestep(t, dt);
nonlinear_solver.calculateNonEquilibriumInitialResiduum(
- process_solutions, process_data->process_id);
+ process_solutions, process_solutions_prev,
+ process_data->process_id);
}
}
NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
- std::vector<GlobalVector*>& x, std::size_t const timestep, double const t,
- double const delta_t, ProcessData const& process_data,
- Output& output_control)
+ 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)
{
auto& process = process_data.process;
int const process_id = process_data.process_id;
@@ -205,7 +214,7 @@ NumLib::NonlinearSolverStatus solveOneTimeStepOneProcess(
};
auto const nonlinear_solver_status =
- nonlinear_solver.solve(x, post_iteration_callback, process_id);
+ nonlinear_solver.solve(x, x_prev, post_iteration_callback, process_id);
if (nonlinear_solver_status.error_norms_met)
{
@@ -262,6 +271,7 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
auto& time_disc = ppd.time_disc;
auto const& x = *_process_solutions[i];
+ auto const& x_prev = *_process_solutions_prev[i];
auto const& conv_crit = ppd.conv_crit;
const MathLib::VecNormType norm_type =
@@ -273,7 +283,7 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
? ((t == timestepper->begin())
? 0. // Always accepts the zeroth step
: time_disc->getRelativeChangeFromPreviousTimestep(
- x, norm_type))
+ x, x_prev, norm_type))
: 0.;
if (!ppd.nonlinear_solver_status.error_norms_met)
@@ -324,7 +334,7 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
bool is_initial_step = false;
// Reset the time step with the minimum step size, dt
- // Update the solution of the previous time step in time_disc.
+ // Update the solution of the previous time step.
for (std::size_t i = 0; i < _per_process_data.size(); i++)
{
const auto& ppd = *_per_process_data[i];
@@ -337,11 +347,11 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
continue;
}
- auto& time_disc = ppd.time_disc;
auto& x = *_process_solutions[i];
+ auto& x_prev = *_process_solutions_prev[i];
if (all_process_steps_accepted)
{
- time_disc->pushState(t, x);
+ MathLib::LinAlg::copy(x, x_prev); // pushState
}
else
{
@@ -352,7 +362,7 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
"Time step {:d} was rejected {:d} times "
"and it will be repeated with a reduced step size.",
accepted_steps + 1, _repeating_times_of_rejected_step);
- time_disc->popState(x);
+ MathLib::LinAlg::copy(x_prev, x); // popState
}
}
}
@@ -431,7 +441,8 @@ void TimeLoop::initialize()
}
// init solution storage
- _process_solutions = setInitialConditions(_start_time, _per_process_data);
+ std::tie(_process_solutions, _process_solutions_prev) =
+ setInitialConditions(_start_time, _per_process_data);
if (_chemical_system != nullptr)
{
@@ -504,8 +515,8 @@ bool TimeLoop::loop()
if (!non_equilibrium_initial_residuum_computed)
{
- calculateNonEquilibriumInitialResiduum(_per_process_data,
- _process_solutions);
+ calculateNonEquilibriumInitialResiduum(
+ _per_process_data, _process_solutions, _process_solutions_prev);
non_equilibrium_initial_residuum_computed = true;
}
@@ -581,13 +592,13 @@ void preTimestepForAllProcesses(
static NumLib::NonlinearSolverStatus solveMonolithicProcess(
const double t, const double dt, const std::size_t timestep_id,
ProcessData const& process_data, std::vector<GlobalVector*>& x,
- Output& output)
+ std::vector<GlobalVector*> const& x_prev, Output& output)
{
BaseLib::RunTime time_timestep_process;
time_timestep_process.start();
- auto const nonlinear_solver_status =
- solveOneTimeStepOneProcess(x, timestep_id, t, dt, process_data, output);
+ auto const nonlinear_solver_status = solveOneTimeStepOneProcess(
+ x, x_prev, timestep_id, t, dt, process_data, output);
INFO("[time] Solving process #{:d} took {:g} s in time step #{:d} ",
process_data.process_id, time_timestep_process.elapsed(), timestep_id);
@@ -634,7 +645,8 @@ NumLib::NonlinearSolverStatus TimeLoop::solveUncoupledEquationSystems(
{
auto const process_id = process_data->process_id;
nonlinear_solver_status = solveMonolithicProcess(
- t, dt, timestep_id, *process_data, _process_solutions, *_output);
+ t, dt, timestep_id, *process_data, _process_solutions,
+ _process_solutions_prev, *_output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
if (!nonlinear_solver_status.error_norms_met)
@@ -704,9 +716,9 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
process_data->process.setCoupledSolutionsForStaggeredScheme(
&coupled_solutions);
- nonlinear_solver_status =
- solveOneTimeStepOneProcess(_process_solutions, timestep_id, t,
- dt, *process_data, *_output);
+ nonlinear_solver_status = solveOneTimeStepOneProcess(
+ _process_solutions, _process_solutions_prev, timestep_id, t, dt,
+ *process_data, *_output);
process_data->nonlinear_solver_status = nonlinear_solver_status;
INFO(
diff --git a/ProcessLib/TimeLoop.h b/ProcessLib/TimeLoop.h
index c97f2297552..a555112d361 100644
--- a/ProcessLib/TimeLoop.h
+++ b/ProcessLib/TimeLoop.h
@@ -110,6 +110,7 @@ private:
private:
std::vector<GlobalVector*> _process_solutions;
+ std::vector<GlobalVector*> _process_solutions_prev;
std::unique_ptr<Output> _output;
std::vector<std::unique_ptr<ProcessData>> _per_process_data;
diff --git a/Tests/NumLib/TimeLoopSingleODE.h b/Tests/NumLib/TimeLoopSingleODE.h
index f0e7b0f46b3..b8986f99962 100644
--- a/Tests/NumLib/TimeLoopSingleODE.h
+++ b/Tests/NumLib/TimeLoopSingleODE.h
@@ -90,9 +90,14 @@ NumLib::NonlinearSolverStatus TimeLoopSingleODE<NLTag>::loop(
xs.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x0));
GlobalVector& x = *xs.back();
+ std::vector<GlobalVector*> xs_prev;
+ xs_prev.push_back(&NumLib::GlobalVectorProvider::provider.getVector(x0));
+ GlobalVector& x_prev = *xs_prev.back();
+
auto& time_disc = _ode_sys.getTimeDiscretization();
- time_disc.setInitialState(t0, x0); // push IC
+ time_disc.setInitialState(t0); // push IC
+ MathLib::LinAlg::copy(x, x_prev); // pushState
_nonlinear_solver->setEquationSystem(_ode_sys, *_convergence_criterion);
@@ -109,13 +114,13 @@ NumLib::NonlinearSolverStatus TimeLoopSingleODE<NLTag>::loop(
int const process_id = 0;
nonlinear_solver_status =
- _nonlinear_solver->solve(xs, nullptr, process_id);
+ _nonlinear_solver->solve(xs, xs_prev, nullptr, process_id);
if (!nonlinear_solver_status.error_norms_met)
{
break;
}
- time_disc.pushState(t, x);
+ MathLib::LinAlg::copy(x, x_prev); // pushState
auto const t_cb =
t; // make sure the callback cannot overwrite anything.
--
GitLab