Skip to content
Snippets Groups Projects
Commit 4466324b authored by Dmitri Naumov's avatar Dmitri Naumov
Browse files

Remove M and K from Newton's assembleWithJacobian

The matrices are not in use by any process.
parent 3b1362c1
No related branches found
No related tags found
No related merge requests found
Showing
with 57 additions and 148 deletions
......@@ -103,55 +103,18 @@ public:
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;
};
......
......@@ -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);
}
......@@ -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);
}
......@@ -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<
......
......@@ -149,13 +149,9 @@ private:
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.
};
......
......@@ -29,8 +29,6 @@ public:
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;
......@@ -40,8 +38,6 @@ public:
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*/)
{
......
......@@ -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()
......
......@@ -33,8 +33,6 @@ public:
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;
......@@ -42,7 +40,6 @@ public:
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;
......
......@@ -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;
......@@ -64,8 +60,6 @@ struct MultiStats
void print() const
{
M.print("M");
K.print("K");
b.print("b");
Jac.print("J");
}
......
......@@ -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_;
};
......
......@@ -175,8 +175,8 @@ GlobalMatrixOutput::GlobalMatrixOutput()
}
void GlobalMatrixOutput::operator()(double const t, int const process_id,
GlobalMatrix const& M,
GlobalMatrix const& K,
GlobalMatrix const* M,
GlobalMatrix const* K,
GlobalVector const& b,
GlobalMatrix const* const Jac)
{
......@@ -188,20 +188,22 @@ void GlobalMatrixOutput::operator()(double const t, int const process_id,
#ifndef USE_PETSC
++counter_;
if (M)
{
auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
process_id, "M", "mat");
fh << "M ";
outputGlobalMatrix(M, fh);
outputGlobalMatrix(*M, fh);
}
if (K)
{
auto fh = openGlobalMatrixOutputFile(filenamePrefix_, counter_, t,
process_id, "K", "mat");
fh << "K ";
outputGlobalMatrix(K, fh);
outputGlobalMatrix(*K, fh);
}
{
......@@ -286,8 +288,8 @@ LocalMatrixOutput::LocalMatrixOutput()
void LocalMatrixOutput::operator()(
double const t, int const process_id, std::size_t const element_id,
std::vector<double> const& local_M_data,
std::vector<double> const& local_K_data,
std::vector<double> const* local_M_data,
std::vector<double> const* local_K_data,
std::vector<double> const& local_b_data,
std::vector<double> const* const local_Jac_data)
{
......@@ -305,16 +307,16 @@ void LocalMatrixOutput::operator()(
fmt::print(fh, "## t = {:.15g}, process id = {}, element id = {}\n\n", t,
process_id, element_id);
if (!local_M_data.empty())
if (local_M_data && !local_M_data->empty())
{
DBUG("... M");
fmt::print(fh, "# M\n{}\n\n", toSquareMatrixRowMajor(local_M_data));
fmt::print(fh, "# M\n{}\n\n", toSquareMatrixRowMajor(*local_M_data));
}
if (!local_K_data.empty())
if (local_K_data && !local_K_data->empty())
{
DBUG("... K");
fmt::print(fh, "# K\n{}\n\n", toSquareMatrixRowMajor(local_K_data));
fmt::print(fh, "# K\n{}\n\n", toSquareMatrixRowMajor(*local_K_data));
}
if (!local_b_data.empty())
......
......@@ -24,8 +24,8 @@ struct GlobalMatrixOutput
{
GlobalMatrixOutput();
void operator()(double const t, int const process_id, GlobalMatrix const& M,
GlobalMatrix const& K, GlobalVector const& b,
void operator()(double const t, int const process_id, GlobalMatrix const* M,
GlobalMatrix const* K, GlobalVector const& b,
GlobalMatrix const* const Jac = nullptr);
private:
......@@ -45,8 +45,8 @@ struct LocalMatrixOutput
void operator()(double const t, int const process_id,
std::size_t const element_id,
std::vector<double> const& local_M_data,
std::vector<double> const& local_K_data,
std::vector<double> const* local_M_data,
std::vector<double> const* local_K_data,
std::vector<double> const& local_b_data,
std::vector<double> const* const local_Jac_data = nullptr);
......
......@@ -27,7 +27,6 @@ void assembleWithJacobianOneElement(
ProcessLib::LocalAssemblerInterface& local_assembler,
const NumLib::LocalToGlobalIndexMap& dof_table, const double t,
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,
ProcessLib::AbstractJacobianAssembler& jacobian_assembler,
......@@ -35,16 +34,14 @@ void assembleWithJacobianOneElement(
{
indices = NumLib::getIndices(mesh_item_id, dof_table);
local_M_data.clear();
local_K_data.clear();
local_b_data.clear();
local_Jac_data.clear();
auto const local_x = x.get(indices);
auto const local_x_prev = x_prev.get(indices);
jacobian_assembler.assembleWithJacobian(
local_assembler, t, dt, local_x, local_x_prev, local_M_data,
local_K_data, local_b_data, local_Jac_data);
jacobian_assembler.assembleWithJacobian(local_assembler, t, dt, local_x,
local_x_prev, local_b_data,
local_Jac_data);
if (local_Jac_data.empty())
{
......@@ -54,8 +51,7 @@ void assembleWithJacobianOneElement(
"current process.");
}
cache.add(local_M_data, local_K_data, local_b_data, local_Jac_data,
indices);
cache.add(local_b_data, local_Jac_data, indices);
}
int getNumberOfThreads()
......@@ -121,7 +117,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
const double t, double const dt, std::vector<GlobalVector*> const& xs,
std::vector<GlobalVector*> const& x_prevs, int const process_id,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac)
GlobalVector& b, GlobalMatrix& Jac)
{
// checks //////////////////////////////////////////////////////////////////
if (process_id != 0)
......@@ -150,8 +146,6 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
// algorithm ///////////////////////////////////////////////////////////////
auto stats = CumulativeStats<MultiStats>::create();
ConcurrentMatrixView M_view(M);
ConcurrentMatrixView K_view(K);
ConcurrentMatrixView b_view(b);
ConcurrentMatrixView Jac_view(Jac);
......@@ -168,8 +162,6 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
// temporary data only stored here in order to avoid frequent memory
// reallocations.
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;
......@@ -178,7 +170,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
auto const jac_asm = jacobian_assembler_.copy();
auto stats_this_thread = stats->clone();
MultiMatrixElementCache cache{M_view, K_view, b_view, Jac_view,
MultiMatrixElementCache cache{b_view, Jac_view,
stats_this_thread->data};
// TODO corner case: what if all elements on a submesh are deactivated?
......@@ -204,8 +196,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
{
assembleWithJacobianOneElement(
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);
local_b_data, local_Jac_data, indices, *jac_asm, cache);
}
catch (...)
{
......@@ -214,9 +205,8 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
continue;
}
local_matrix_output_(t, process_id, element_id, local_M_data,
local_K_data, local_b_data,
&local_Jac_data);
local_matrix_output_(t, process_id, element_id, nullptr,
nullptr, local_b_data, &local_Jac_data);
}
}
else
......@@ -242,8 +232,7 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
{
assembleWithJacobianOneElement(
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);
local_b_data, local_Jac_data, indices, *jac_asm, cache);
}
catch (...)
{
......@@ -252,16 +241,15 @@ void ParallelVectorMatrixAssembler::assembleWithJacobian(
continue;
}
local_matrix_output_(t, process_id, element_id, local_M_data,
local_K_data, local_b_data,
&local_Jac_data);
local_matrix_output_(t, process_id, element_id, nullptr,
nullptr, local_b_data, &local_Jac_data);
}
}
} // OpenMP parallel section
stats->print();
global_matrix_output_(t, process_id, M, K, b, &Jac);
global_matrix_output_(t, process_id, nullptr, nullptr, b, &Jac);
exception.rethrow();
}
} // namespace ProcessLib::Assembly
......@@ -30,7 +30,7 @@ public:
std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_tables,
const double t, double const dt, std::vector<GlobalVector*> const& xs,
std::vector<GlobalVector*> const& x_prevs, int const process_id,
GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac);
GlobalVector& b, GlobalMatrix& Jac);
private:
AbstractJacobianAssembler& jacobian_assembler_;
......
......@@ -155,8 +155,7 @@ public:
void assembleWithJacobian(double const 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)
{
DBUG("AssemblyMixin assembleWithJacobian(t={}, dt={}, process_id={}).",
......@@ -179,7 +178,7 @@ public:
pvma_.assembleWithJacobian(loc_asms, sad.active_element_ids,
dof_tables, t, dt, x, x_prev,
process_id, M, K, b_submesh, Jac);
process_id, b_submesh, Jac);
MathLib::LinAlg::axpy(b, 1.0, b_submesh);
......@@ -198,7 +197,7 @@ public:
pvma_.assembleWithJacobian(loc_asms, pv.getActiveElementIDs(),
dof_tables, t, dt, x, x_prev, process_id,
M, K, b, Jac);
b, Jac);
}
AssemblyMixinBase::copyResiduumVectorsToBulkMesh(
......
......@@ -39,7 +39,7 @@ public:
//! \c K and the vector \c b.
virtual void applyNaturalBC(const double /*t*/,
std::vector<GlobalVector*> const& /*x*/,
int const /*process_id*/, GlobalMatrix& /*K*/,
int const /*process_id*/, GlobalMatrix* /*K*/,
GlobalVector& /*b*/, GlobalMatrix* /*Jac*/)
{
// By default it is assumed that the BC is not a natural BC. Therefore
......
......@@ -14,7 +14,7 @@ namespace ProcessLib
{
void BoundaryConditionCollection::applyNaturalBC(
const double t, std::vector<GlobalVector*> const& x, int const process_id,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac) const
GlobalMatrix* K, GlobalVector& b, GlobalMatrix* Jac) const
{
for (auto const& bc : _boundary_conditions)
{
......
......@@ -28,7 +28,7 @@ public:
}
void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
int const process_id, GlobalMatrix& K, GlobalVector& b,
int const process_id, GlobalMatrix* K, GlobalVector& b,
GlobalMatrix* Jac) const;
std::vector<NumLib::IndexValueVector<GlobalIndexType>> const*
......
......@@ -82,7 +82,7 @@ void GenericNaturalBoundaryCondition<BoundaryConditionData,
applyNaturalBC(const double t,
std::vector<GlobalVector*> const& x,
int const process_id,
GlobalMatrix& K,
GlobalMatrix* K,
GlobalVector& b,
GlobalMatrix* Jac)
{
......
......@@ -36,7 +36,7 @@ public:
/// Calls local assemblers which calculate their contributions to the global
/// matrix and the right-hand-side.
void applyNaturalBC(const double t, std::vector<GlobalVector*> const& x,
int const process_id, GlobalMatrix& K, GlobalVector& b,
int const process_id, GlobalMatrix* K, GlobalVector& b,
GlobalMatrix* Jac) override;
private:
......
......@@ -27,7 +27,7 @@ public:
std::size_t const id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary, double const t,
std::vector<GlobalVector*> const& x, int const process_id,
GlobalMatrix& K, GlobalVector& b, GlobalMatrix* Jac) = 0;
GlobalMatrix* K, GlobalVector& b, GlobalMatrix* Jac) = 0;
};
template <typename ShapeFunction, int GlobalDim>
......
......@@ -57,7 +57,7 @@ public:
void assemble(std::size_t const mesh_item_id,
NumLib::LocalToGlobalIndexMap const& dof_table_boundary,
double const t, std::vector<GlobalVector*> const& x,
int const process_id, GlobalMatrix& /*K*/, GlobalVector& b,
int const process_id, GlobalMatrix* /*K*/, GlobalVector& b,
GlobalMatrix* /*Jac*/) override
{
NodalVectorType _local_rhs = NodalVectorType::Zero(_local_matrix_size);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment