Skip to content

Commit

Permalink
Merge branch 'UseXPrevInsteadOfXDot' into 'master'
Browse files Browse the repository at this point in the history 
Replace x_dot with x_prev

See merge request ogs/ogs!4699
  • Loading branch information
@endJunction
endJunction committed Aug 5, 2023
2 parents 87635c0 + 699ef8d commit 95ab096
Show file tree
Hide file tree
Showing 153 changed files with 995 additions and 1,164 deletions.
16 changes: 10 additions & 6 deletions NumLib/ODESolver/MatrixTranslator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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>::
Expand Down
8 changes: 4 additions & 4 deletions NumLib/ODESolver/MatrixTranslator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down
4 changes: 2 additions & 2 deletions NumLib/ODESolver/ODESystem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;

Expand Down Expand Up @@ -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;
Expand Down
13 changes: 0 additions & 13 deletions NumLib/ODESolver/TimeDiscretization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down
6 changes: 0 additions & 6 deletions NumLib/ODESolver/TimeDiscretization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
56 changes: 5 additions & 51 deletions NumLib/ODESolver/TimeDiscretizedODESystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<
Expand All @@ -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<
Expand Down Expand Up @@ -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<
Expand Down
5 changes: 0 additions & 5 deletions NumLib/ODESolver/TimeDiscretizedODESystem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
};

//! @}
Expand Down
4 changes: 2 additions & 2 deletions ProcessLib/AbstractJacobianAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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*/,
Expand Down
8 changes: 4 additions & 4 deletions ProcessLib/AnalyticalJacobianAssembler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}

Expand Down
4 changes: 2 additions & 2 deletions ProcessLib/AnalyticalJacobianAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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;
Expand Down
20 changes: 10 additions & 10 deletions ProcessLib/Assembly/ParallelVectorMatrixAssembler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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())
{
Expand Down Expand Up @@ -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 //////////////////////////////////////////////////////////////////
Expand All @@ -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 ///////////////////////////////////////////////////////////////

Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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);
}
Expand Down
2 changes: 1 addition & 1 deletion ProcessLib/Assembly/ParallelVectorMatrixAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Loading

0 comments on commit 95ab096

Please sign in to comment.