Skip to content

Commit

Permalink
Merge branch 'PassGlobalMatrixPointersToAssemblers' into 'master'
Browse files Browse the repository at this point in the history 
Remove M and K from Newton's assembleWithJacobian

See merge request ogs/ogs!5050
  • Loading branch information
@endJunction
endJunction committed Jul 18, 2024
2 parents 481f189 + 3e9a70c commit dfb90fa
Show file tree
Hide file tree
Showing 129 changed files with 432 additions and 697 deletions.
41 changes: 2 additions & 39 deletions NumLib/ODESolver/ODESystem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -103,55 +103,18 @@ class ODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
void preAssemble(const double t, double const dt,
GlobalVector const& x) override = 0;

/*! Assemble \c M, \c K, \c b and the Jacobian
/*! Assemble \c b and the Jacobian
* \f$ \mathtt{Jac} := \partial r/\partial x_N \f$
* at the provided state (\c t, \c x).
*
* For the meaning of the other parameters refer to the the introductory
* remarks on
* \ref concept_time_discretization "time discretization".
*
* \remark
* \parblock
* The Jacobian will be generally of the following form:
* \f[ \mathtt{Jac} := \frac{\partial r(x_C, t_C)}{\partial x_N} =
* M \cdot \frac{\partial \hat x}{\partial x_N}
* + \frac{\partial M}{\partial x_N} \cdot \hat x
* + K \cdot \frac{\partial x_C}{\partial x_N}
* + \frac{\partial K}{\partial x_N} \cdot x_N
* + \frac{\partial b}{\partial x_N},
* \f]
* where \f$ M \f$, \f$ K \f$ and \f$ b \f$ are matrix-valued
* (vector-valued, respectively)
* functions that depend on \f$ x_C \f$ and \f$ t_C \f$.
*
* Due to the arguments provided to this method its implementation only has
* to
* compute the derivatives
* \f$ \frac{\partial M}{\partial x_N} \cdot \hat x \f$,
* \f$ \frac{\partial K}{\partial x_N} \cdot x_N \f$ and
* \f$ \frac{\partial b}{\partial x_N} \f$.
* The other terms can be readily taken from the method parameters.
*
* In particular for the ForwardEuler time discretization scheme the
* equation will
* collapse to
* \f$ \mathtt{Jac} =
* M \cdot \frac{\partial \hat x}{\partial x_N}
* \f$
* since in that scheme \f$ x_N \neq x_C \f$.
*
* Of course, the implementation of this method is allowed to compute the
* Jacobian in a
* different way, as long as that is consistent with the definition of \f$
* \mathtt{Jac} \f$.
* \endparblock
*/
virtual void assembleWithJacobian(const double t, double const dt,
std::vector<GlobalVector*> const& x,
std::vector<GlobalVector*> const& x_prev,
int const process_id, GlobalMatrix& M,
GlobalMatrix& K, GlobalVector& b,
int const process_id, GlobalVector& b,
GlobalMatrix& Jac) = 0;
};

Expand Down
32 changes: 10 additions & 22 deletions NumLib/ODESolver/TimeDiscretizedODESystem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,10 +55,6 @@ TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
{
_Jac = &NumLib::GlobalMatrixProvider::provider.getMatrix(
_ode.getMatrixSpecifications(process_id), _Jac_id);
_M = &NumLib::GlobalMatrixProvider::provider.getMatrix(
_ode.getMatrixSpecifications(process_id), _M_id);
_K = &NumLib::GlobalMatrixProvider::provider.getMatrix(
_ode.getMatrixSpecifications(process_id), _K_id);
_b = &NumLib::GlobalVectorProvider::provider.getVector(
_ode.getMatrixSpecifications(process_id), _b_id);
}
Expand All @@ -68,8 +64,6 @@ TimeDiscretizedODESystem<
NonlinearSolverTag::Newton>::~TimeDiscretizedODESystem()
{
NumLib::GlobalMatrixProvider::provider.releaseMatrix(*_Jac);
NumLib::GlobalMatrixProvider::provider.releaseMatrix(*_M);
NumLib::GlobalMatrixProvider::provider.releaseMatrix(*_K);
NumLib::GlobalVectorProvider::provider.releaseVector(*_b);
}

Expand All @@ -79,35 +73,29 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
std::vector<GlobalVector*> const& x_prev,
int const process_id)
{
namespace LinAlg = MathLib::LinAlg;

auto const t = _time_disc.getCurrentTime();
auto const dt = _time_disc.getCurrentTimeIncrement();
auto const& x_curr = *x_new_timestep[process_id];

_M->setZero();
_K->setZero();
_b->setZero();
_Jac->setZero();

_ode.preAssemble(t, dt, x_curr);
_ode.assembleWithJacobian(t, dt, x_new_timestep, x_prev, process_id, *_M,
*_K, *_b, *_Jac);
_ode.assembleWithJacobian(t, dt, x_new_timestep, x_prev, process_id, *_b,
*_Jac);

LinAlg::finalizeAssembly(*_M);
LinAlg::finalizeAssembly(*_K);
LinAlg::finalizeAssembly(*_b);
MathLib::LinAlg::finalizeAssembly(*_b);
MathLib::LinAlg::finalizeAssembly(*_Jac);
}

void TimeDiscretizedODESystem<
ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::getResidual(GlobalVector const& x_new_timestep,
GlobalVector const& x_prev,
GlobalVector& res) const
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
getResidual(GlobalVector const& /*x_new_timestep*/,
GlobalVector const& /*x_prev*/,
GlobalVector& res) const
{
double const dt = _time_disc.getCurrentTimeIncrement();
_mat_trans->computeResidual(*_M, *_K, *_b, dt, x_new_timestep, x_prev, res);
MathLib::LinAlg::copy(*_b, res); // res = b
MathLib::LinAlg::scale(res, -1.); // res = -b
}

void TimeDiscretizedODESystem<
Expand Down
4 changes: 0 additions & 4 deletions NumLib/ODESolver/TimeDiscretizedODESystem.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,13 +149,9 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
nullptr; //!< stores precomputed values for known solutions

GlobalMatrix* _Jac; //!< the Jacobian of the residual
GlobalMatrix* _M; //!< Matrix \f$ M \f$.
GlobalMatrix* _K; //!< Matrix \f$ K \f$.
GlobalVector* _b; //!< Matrix \f$ b \f$.

std::size_t _Jac_id = 0u; //!< ID of the \c _Jac matrix.
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.
};

Expand Down
4 changes: 0 additions & 4 deletions ProcessLib/AbstractJacobianAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,6 @@ class AbstractJacobianAssembler
double const t, double const dt,
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) = 0;

Expand All @@ -40,8 +38,6 @@ class AbstractJacobianAssembler
LocalAssemblerInterface& /*local_assembler*/, double const /*t*/,
double const /*dt*/, 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*/)
{
Expand Down
8 changes: 2 additions & 6 deletions ProcessLib/AnalyticalJacobianAssembler.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,24 +17,20 @@ 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_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_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_x_prev,
int const process_id, std::vector<double>& local_M_data,
std::vector<double>& local_K_data, std::vector<double>& local_b_data,
int const process_id, std::vector<double>& local_b_data,
std::vector<double>& local_Jac_data)
{
local_assembler.assembleWithJacobianForStaggeredScheme(
t, dt, local_x, local_x_prev, process_id, local_M_data, local_K_data,
local_b_data, local_Jac_data);
t, dt, local_x, local_x_prev, process_id, local_b_data, local_Jac_data);
}

std::unique_ptr<AbstractJacobianAssembler> AnalyticalJacobianAssembler::copy()
Expand Down
3 changes: 0 additions & 3 deletions ProcessLib/AnalyticalJacobianAssembler.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,16 +33,13 @@ class AnalyticalJacobianAssembler final : public AbstractJacobianAssembler
double const t, double const dt,
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) override;

void assembleWithJacobianForStaggeredScheme(
LocalAssemblerInterface& local_assembler, double const t,
double const dt, 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) override;

Expand Down
22 changes: 11 additions & 11 deletions ProcessLib/Assembly/AssembledMatrixCache.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ struct AssembledMatrixCache final
void assemble(
const double t, double const dt, std::vector<GlobalVector*> const& x,
std::vector<GlobalVector*> const& x_prev, int const process_id,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
GlobalMatrix* const M, GlobalMatrix* const K, GlobalVector* const b,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
VectorMatrixAssembler& global_assembler,
VectorOfLocalAssemblers const& local_assemblers,
Expand All @@ -53,7 +53,7 @@ template <typename VectorOfLocalAssemblers>
void AssembledMatrixCache::assemble(
const double t, double const dt, std::vector<GlobalVector*> const& x,
std::vector<GlobalVector*> const& x_prev, int const process_id,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b,
GlobalMatrix* const M, GlobalMatrix* const K, GlobalVector* const b,
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
VectorMatrixAssembler& global_assembler,
VectorOfLocalAssemblers const& local_assemblers,
Expand All @@ -70,9 +70,9 @@ void AssembledMatrixCache::assemble(
local_assemblers, active_element_ids, dof_tables, t, dt, x, x_prev,
process_id, M, K, b);

MathLib::finalizeMatrixAssembly(M);
MathLib::finalizeMatrixAssembly(K);
MathLib::finalizeVectorAssembly(b);
MathLib::finalizeMatrixAssembly(*M);
MathLib::finalizeMatrixAssembly(*K);
MathLib::finalizeVectorAssembly(*b);

INFO("[time] Calling local assemblers took {:g} s", time_asm.elapsed());

Expand All @@ -83,9 +83,9 @@ void AssembledMatrixCache::assemble(
BaseLib::RunTime time_save;
time_save.start();

K_ = MathLib::MatrixVectorTraits<GlobalMatrix>::newInstance(K);
M_ = MathLib::MatrixVectorTraits<GlobalMatrix>::newInstance(M);
b_ = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(b);
K_ = MathLib::MatrixVectorTraits<GlobalMatrix>::newInstance(*K);
M_ = MathLib::MatrixVectorTraits<GlobalMatrix>::newInstance(*M);
b_ = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(*b);

INFO("[time] Saving global K, M, b took {:g} s",
time_save.elapsed());
Expand All @@ -98,9 +98,9 @@ void AssembledMatrixCache::assemble(
BaseLib::RunTime time_restore;
time_restore.start();

MathLib::LinAlg::copy(*K_, K);
MathLib::LinAlg::copy(*M_, M);
MathLib::LinAlg::copy(*b_, b);
MathLib::LinAlg::copy(*K_, *K);
MathLib::LinAlg::copy(*M_, *M);
MathLib::LinAlg::copy(*b_, *b);

INFO("[time] Restoring global K, M, b took {:g} s",
time_restore.elapsed());
Expand Down
6 changes: 0 additions & 6 deletions ProcessLib/Assembly/MatrixAssemblyStats.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,11 @@ struct Stats

struct MultiStats
{
Stats M;
Stats K;
Stats b;
Stats Jac;

MultiStats& operator+=(MultiStats const& other)
{
M += other.M;
K += other.K;
b += other.b;
Jac += other.Jac;

Expand All @@ -64,8 +60,6 @@ struct MultiStats

void print() const
{
M.print("M");
K.print("K");
b.print("b");
Jac.print("J");
}
Expand Down
16 changes: 3 additions & 13 deletions ProcessLib/Assembly/MatrixElementCache.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -257,31 +257,21 @@ class MultiMatrixElementCache final
using GlobalVectorView = ConcurrentMatrixView<1>;

public:
MultiMatrixElementCache(GlobalMatrixView& M, GlobalMatrixView& K,
GlobalVectorView& b, GlobalMatrixView& Jac,
MultiMatrixElementCache(GlobalVectorView& b, GlobalMatrixView& Jac,
MultiStats& stats)
: cache_M_(M, stats.M),
cache_K_(K, stats.K),
cache_b_(b, stats.b),
cache_Jac_(Jac, stats.Jac)
: cache_b_(b, stats.b), cache_Jac_(Jac, stats.Jac)
{
}

void add(std::vector<double> const& local_M_data,
std::vector<double> const& local_K_data,
std::vector<double> const& local_b_data,
void add(std::vector<double> const& local_b_data,
std::vector<double> const& local_Jac_data,
std::vector<GlobalIndexType> const& indices)
{
cache_M_.add(local_M_data, indices);
cache_K_.add(local_K_data, indices);
cache_b_.add(local_b_data, indices);
cache_Jac_.add(local_Jac_data, indices);
}

private:
MatrixElementCache<2> cache_M_;
MatrixElementCache<2> cache_K_;
MatrixElementCache<1> cache_b_;
MatrixElementCache<2> cache_Jac_;
};
Expand Down
Loading

0 comments on commit dfb90fa

Please sign in to comment.