diff --git a/NumLib/ODESolver/MatrixTranslator.cpp b/NumLib/ODESolver/MatrixTranslator.cpp
index e6cf7ee6dda85c3b10ded49bcc98fcba1dbffe65..ff7d31b25cb944d3ecc1e939ff86d2f4d6486676 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 2008c173bd5a76978dfe872b1c8455764e89f0a7..33a4ba00098e75762c14ff59aa170848fa4f5e4b 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 8d7127a480fb7a9847afb542544833fa72c79273..cd0195360dec1ca65a6dc8e481796e400ccc02fb 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 6b2509668574ab998f5045f8dd0c7e25b761d986..26f2c94d228fa91dda5ad1f7b0add34a5a8e7079 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 85c3787af2f60725b6681b62ea8ff05bbef629e8..6094520a50ca11f5d51b162b0f21abd1c0e6c0bb 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 08f9c684a6d16b8ae610e958d6ec893fb150adf5..551384e4b4f435449bff5ac92f587511b92a3a1b 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 f53765d4f6c2604b6eab80b82dd5d640b6fd6ee2..f5dcbab099f2a5cfabe56719e4056eb70695f6e9 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 85bd3bd4e080491b0f6377752c58d219dad53a69..21e603274906820cfc0842e54b68c743b7d1a8d2 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 e28470e0f3293a196f157cf412356c4d76b2942f..fbbe1a84a4692d8237dde58a09132cf8248de4ff 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 e7f43b569ed97d0b532430f941d596b7e16dbe8e..53d3da26793d0c817418f1ee3cc17cd9064626b4 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 c97f2297552fc162960943a2d922c508a1e0cd53..a555112d361b1fc7309620c529cbd4be6fab1db5 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 f0e7b0f46b3ad3cecd19c2a474d558ceeaa8c831..b8986f99962c40dd02f4d8279f47c98e60d954c7 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.