Merge branch 'UseXPrevInsteadOfXDot' into 'master'

Replace x_dot with x_prev

See merge request ogs/ogs!4699

NumLib/ODESolver/MatrixTranslator.cpp

@@ -45,16 +45,20 @@ void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
 
 void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::
  computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
- GlobalVector const& b, GlobalVector const& x_curr,
- GlobalVector const& xdot, GlobalVector& res) const
+ GlobalVector const& b, double const dt,
+ GlobalVector const& x_curr, GlobalVector const& x_prev,
+ GlobalVector& res) const
 {
  namespace LinAlg = MathLib::LinAlg;
 
  // res = M * x_dot + K * x_curr - b
- LinAlg::matMult(M, xdot, res); // the local vector x_dot seems to be
- // necessary because of this multiplication
- LinAlg::matMultAdd(K, x_curr, res, res);
- LinAlg::axpy(res, -1.0, b);
+ GlobalVector x_dot;
+ LinAlg::copy(x_curr, x_dot); // x_dot = x
+ LinAlg::axpy(x_dot, -1., x_prev); // x_dot = x - x_prev
+ LinAlg::scale(x_dot, 1. / dt); // x_dot = (x - x_prev)/dt
+ LinAlg::matMult(M, x_dot, res); // res = M*x_dot
+ LinAlg::matMultAdd(K, x_curr, res, res); // res = M*x_dot + K*x
+ LinAlg::axpy(res, -1., b); // res = M*x_dot + X*x - b
 }
 
 void MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>::

NumLib/ODESolver/MatrixTranslator.h

@@ -54,9 +54,9 @@ class MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
  * \ref concept_time_discretization "time discretization".
  */
  virtual void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
- GlobalVector const& b,
+ GlobalVector const& b, double dt,
  GlobalVector const& x_curr,
- GlobalVector const& xdot,
+ GlobalVector const& x_prev,
  GlobalVector& res) const = 0;
 
  //! Computes the Jacobian of the residual and writes it to \c Jac_out.
@@ -107,8 +107,8 @@ class MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>
 
  //! Computes \f$ r = M \cdot \hat x + K \cdot x_C - b \f$.
  void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
- GlobalVector const& b, GlobalVector const& x_curr,
- GlobalVector const& xdot,
+ GlobalVector const& b, double dt,
+ GlobalVector const& x_curr, GlobalVector const& x_prev,
  GlobalVector& res) const override;
 
  //! Writes \c Jac_in to \c Jac_out.

NumLib/ODESolver/ODESystem.h

@@ -59,7 +59,7 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
  //! Assemble \c M, \c K and \c b at the provided state (\c t, \c x).
  virtual void assemble(const double t, double const dt,
  std::vector<GlobalVector*> const& x,
- std::vector<GlobalVector*> const& xdot,
+ std::vector<GlobalVector*> const& x_prev,
  int const process_id, GlobalMatrix& M,
  GlobalMatrix& K, GlobalVector& b) = 0;
 
@@ -141,7 +141,7 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
  */
  virtual void assembleWithJacobian(const double t, double const dt,
  std::vector<GlobalVector*> const& x,
- std::vector<GlobalVector*> const& xdot,
+ std::vector<GlobalVector*> const& x_prev,
  int const process_id, GlobalMatrix& M,
  GlobalMatrix& K, GlobalVector& b,
  GlobalMatrix& Jac) = 0;

NumLib/ODESolver/TimeDiscretization.cpp

@@ -48,19 +48,6 @@ double computeRelativeChangeFromPreviousTimestep(GlobalVector const& x,
  return norm_dx / std::numeric_limits<double>::epsilon();
 }
 
-void TimeDiscretization::getXdot(GlobalVector const& x_at_new_timestep,
- GlobalVector const& x_old,
- GlobalVector& xdot) const
-{
- namespace LinAlg = MathLib::LinAlg;
-
- double const dt = getCurrentTimeIncrement();
-
- // xdot = 1/dt * x_at_new_timestep - x_old
- getWeightedOldX(xdot, x_old);
- LinAlg::axpby(xdot, 1. / dt, -1.0, x_at_new_timestep);
-}
-
 void BackwardEuler::getWeightedOldX(GlobalVector& y,
  GlobalVector const& x_old) const
 {

NumLib/ODESolver/TimeDiscretization.h

@@ -128,12 +128,6 @@ class TimeDiscretization
  //! assembled.
  virtual double getCurrentTimeIncrement() const = 0;
 
- //! Returns \f$ \hat x \f$, i.e. the discretized approximation of \f$ \dot x
- //! \f$.
- void getXdot(GlobalVector const& x_at_new_timestep,
- GlobalVector const& x_old,
- GlobalVector& xdot) const;
-
  //! Returns \f$ x_O \f$.
  virtual void getWeightedOldX(
  GlobalVector& y, GlobalVector const& x_old) const = 0; // = x_old

NumLib/ODESolver/TimeDiscretizedODESystem.cpp

@@ -82,44 +82,19 @@ 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_ids[i]);
- _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
- }
-
  _M->setZero();
  _K->setZero();
  _b->setZero();
  _Jac->setZero();
 
  _ode.preAssemble(t, dt, x_curr);
- try
- {
- _ode.assembleWithJacobian(t, dt, x_new_timestep, xdot, process_id, *_M,
- *_K, *_b, *_Jac);
- }
- catch (AssemblyException const&)
- {
- for (auto& v : xdot)
- {
- NumLib::GlobalVectorProvider::provider.releaseVector(*v);
- }
- throw;
- }
+ _ode.assembleWithJacobian(t, dt, x_new_timestep, x_prev, process_id, *_M,
+ *_K, *_b, *_Jac);
 
  LinAlg::finalizeAssembly(*_M);
  LinAlg::finalizeAssembly(*_K);
  LinAlg::finalizeAssembly(*_b);
  MathLib::LinAlg::finalizeAssembly(*_Jac);
-
- for (auto& v : xdot)
- {
- NumLib::GlobalVectorProvider::provider.releaseVector(*v);
- }
 }
 
 void TimeDiscretizedODESystem<
@@ -128,15 +103,8 @@ void TimeDiscretizedODESystem<
  GlobalVector const& x_prev,
  GlobalVector& res) const
 {
- // TODO Maybe the duplicate calculation of xdot here and in assembleJacobian
- // can be optimized. However, that would make the interface a bit more
- // fragile.
- auto& xdot = NumLib::GlobalVectorProvider::provider.getVector(_xdot_id);
- _time_disc.getXdot(x_new_timestep, x_prev, xdot);
-
- _mat_trans->computeResidual(*_M, *_K, *_b, x_new_timestep, xdot, res);
-
- NumLib::GlobalVectorProvider::provider.releaseVector(xdot);
+ double const dt = _time_disc.getCurrentTimeIncrement();
+ _mat_trans->computeResidual(*_M, *_K, *_b, dt, x_new_timestep, x_prev, res);
 }
 
 void TimeDiscretizedODESystem<
@@ -220,31 +188,17 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
  auto const t = _time_disc.getCurrentTime();
  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_ids[i]);
- _time_disc.getXdot(*x_new_timestep[i], *x_prev[i], *xdot[i]);
- }
 
  _M->setZero();
  _K->setZero();
  _b->setZero();
 
  _ode.preAssemble(t, dt, x_curr);
- _ode.assemble(t, dt, x_new_timestep, xdot, process_id, *_M, *_K, *_b);
+ _ode.assemble(t, dt, x_new_timestep, x_prev, process_id, *_M, *_K, *_b);
 
  LinAlg::finalizeAssembly(*_M);
  LinAlg::finalizeAssembly(*_K);
  LinAlg::finalizeAssembly(*_b);
-
- for (auto& v : xdot)
- {
- NumLib::GlobalVectorProvider::provider.releaseVector(*v);
- }
 }
 
 void TimeDiscretizedODESystem<

NumLib/ODESolver/TimeDiscretizedODESystem.h

@@ -153,10 +153,6 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
  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.
-
- //! 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;
@@ -256,7 +252,6 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
  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;
 };
 
 //! @}

ProcessLib/AbstractJacobianAssembler.h

@@ -29,7 +29,7 @@ class AbstractJacobianAssembler
  virtual void assembleWithJacobian(LocalAssemblerInterface& local_assembler,
  double const t, double const dt,
  std::vector<double> const& local_x,
- std::vector<double> const& local_xdot,
+ std::vector<double> const& local_x_prev,
  std::vector<double>& local_M_data,
  std::vector<double>& local_K_data,
  std::vector<double>& local_b_data,
@@ -40,7 +40,7 @@ class AbstractJacobianAssembler
  virtual void assembleWithJacobianForStaggeredScheme(
  LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
  double const /*dt*/, Eigen::VectorXd const& /*local_x*/,
- Eigen::VectorXd const& /*local_xdot*/, int const /*process_id*/,
+ Eigen::VectorXd const& /*local_x_prev*/, int const /*process_id*/,
  std::vector<double>& /*local_M_data*/,
  std::vector<double>& /*local_K_data*/,
  std::vector<double>& /*local_b_data*/,

ProcessLib/AnalyticalJacobianAssembler.cpp

@@ -17,24 +17,24 @@ namespace ProcessLib
 {
 void AnalyticalJacobianAssembler::assembleWithJacobian(
  LocalAssemblerInterface& local_assembler, double const t, double const dt,
- std::vector<double> const& local_x, std::vector<double> const& local_xdot,
+ std::vector<double> const& local_x, std::vector<double> const& local_x_prev,
  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data)
 {
- local_assembler.assembleWithJacobian(t, dt, local_x, local_xdot,
+ local_assembler.assembleWithJacobian(t, dt, local_x, local_x_prev,
  local_M_data, local_K_data,
  local_b_data, local_Jac_data);
 }
 
 void AnalyticalJacobianAssembler::assembleWithJacobianForStaggeredScheme(
  LocalAssemblerInterface& local_assembler, double const t, double const dt,
- Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_xdot,
+ Eigen::VectorXd const& local_x, Eigen::VectorXd const& local_x_prev,
  int const process_id, std::vector<double>& local_M_data,
  std::vector<double>& local_K_data, std::vector<double>& local_b_data,
  std::vector<double>& local_Jac_data)
 {
  local_assembler.assembleWithJacobianForStaggeredScheme(
- t, dt, local_x, local_xdot, process_id, local_M_data, local_K_data,
+ t, dt, local_x, local_x_prev, process_id, local_M_data, local_K_data,
  local_b_data, local_Jac_data);
 }
 

ProcessLib/AnalyticalJacobianAssembler.h

@@ -34,7 +34,7 @@ class AnalyticalJacobianAssembler final : public AbstractJacobianAssembler
  void assembleWithJacobian(LocalAssemblerInterface& local_assembler,
  double const t, double const dt,
  std::vector<double> const& local_x,
- std::vector<double> const& local_xdot,
+ std::vector<double> const& local_x_prev,
  std::vector<double>& local_M_data,
  std::vector<double>& local_K_data,
  std::vector<double>& local_b_data,
@@ -43,7 +43,7 @@ class AnalyticalJacobianAssembler final : public AbstractJacobianAssembler
  void assembleWithJacobianForStaggeredScheme(
  LocalAssemblerInterface& local_assembler, double const t,
  double const dt, Eigen::VectorXd const& local_x,
- Eigen::VectorXd const& local_xdot, int const process_id,
+ Eigen::VectorXd const& local_x_prev, int const process_id,
  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
  std::vector<double>& local_b_data,
  std::vector<double>& local_Jac_data) override;

ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp

@@ -26,7 +26,7 @@ void assembleWithJacobianOneElement(
  const std::size_t mesh_item_id,
  ProcessLib::LocalAssemblerInterface& local_assembler,
  const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
- const double dt, const GlobalVector& x, const GlobalVector& xdot,
+ const double dt, const GlobalVector& x, const GlobalVector& x_prev,
  std::vector<double>& local_M_data, std::vector<double>& local_K_data,
  std::vector<double>& local_b_data, std::vector<double>& local_Jac_data,
  std::vector<GlobalIndexType>& indices,
@@ -41,10 +41,10 @@ void assembleWithJacobianOneElement(
  local_Jac_data.clear();
 
  auto const local_x = x.get(indices);
- auto const local_xdot = xdot.get(indices);
+ auto const local_x_prev = x_prev.get(indices);
  jacobian_assembler.assembleWithJacobian(
- local_assembler, t, dt, local_x, local_xdot, local_M_data, local_K_data,
- local_b_data, local_Jac_data);
+ local_assembler, t, dt, local_x, local_x_prev, local_M_data,
+ local_K_data, local_b_data, local_Jac_data);
 
  if (local_Jac_data.empty())
  {
@@ -121,7 +121,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
  std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
  dof_tables,
  const double t, double const dt, std::vector<GlobalVector*> const& xs,
- std::vector<GlobalVector*> const& xdots, int const process_id,
+ std::vector<GlobalVector*> const& x_prevs, int const process_id,
  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac)
 {
  // checks //////////////////////////////////////////////////////////////////
@@ -142,11 +142,11 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
  }
  auto const& x = *xs.front();
 
- if (xdots.size() != 1)
+ if (x_prevs.size() != 1)
  {
- OGS_FATAL("More than 1 xdot vector");
+ OGS_FATAL("More than 1 x_prev vector");
  }
- auto const& xdot = *xdots.front();
+ auto const& x_prev = *x_prevs.front();
 
  // algorithm ///////////////////////////////////////////////////////////////
 
@@ -204,7 +204,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
  try
  {
  assembleWithJacobianOneElement(
- element_id, loc_asm, dof_table, t, dt, x, xdot,
+ element_id, loc_asm, dof_table, t, dt, x, x_prev,
  local_M_data, local_K_data, local_b_data,
  local_Jac_data, indices, *jac_asm, cache);
  }
@@ -242,7 +242,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
  try
  {
  assembleWithJacobianOneElement(
- element_id, loc_asm, dof_table, t, dt, x, xdot,
+ element_id, loc_asm, dof_table, t, dt, x, x_prev,
  local_M_data, local_K_data, local_b_data,
  local_Jac_data, indices, *jac_asm, cache);
  }

ProcessLib/Assembly/ParallelVectorMatrixAssembler.h

@@ -31,7 +31,7 @@ class ParallelVectorMatrixAssembler
  std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
  dof_tables,
  const double t, double const dt, std::vector<GlobalVector*> const& xs,
- std::vector<GlobalVector*> const& xdots, int const process_id,
+ std::vector<GlobalVector*> const& x_prevs, int const process_id,
  GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac);
 
 private: