diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index ba32dd836d7c5b19c54abe5f20dd98f2686609af..832581f3afee4101640df133cfd494d176891d61 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -37,13 +37,11 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(
_rhs_id);
- auto& x_new =
- NumLib::GlobalVectorProvider::provider.getVector(
- _x_new_id);
+ auto& x_new = NumLib::GlobalVectorProvider::provider.getVector(_x_new_id);
bool error_norms_met = false;
- LinAlg::copy(x, x_new); // set initial guess, TODO save the copy
+ LinAlg::copy(x, x_new); // set initial guess
_convergence_criterion->preFirstIteration();
@@ -51,27 +49,34 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
for (; iteration <= _maxiter;
++iteration, _convergence_criterion->reset())
{
+ BaseLib::RunTime timer_dirichlet;
+ double time_dirichlet = 0.0;
+
BaseLib::RunTime time_iteration;
time_iteration.start();
- sys.preIteration(iteration, x);
+ timer_dirichlet.start();
+ sys.computeKnownSolutions(x_new);
+ sys.applyKnownSolutions(x_new);
+ time_dirichlet += timer_dirichlet.elapsed();
+
+ sys.preIteration(iteration, x_new);
BaseLib::RunTime time_assembly;
time_assembly.start();
- sys.assemble(x);
+ sys.assemble(x_new);
sys.getA(A);
sys.getRhs(rhs);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
- BaseLib::RunTime time_dirichlet;
- time_dirichlet.start();
- // Here _x_new has to be used and it has to be equal to x!
+ timer_dirichlet.start();
sys.applyKnownSolutionsPicard(A, rhs, x_new);
- INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
+ time_dirichlet += timer_dirichlet.elapsed();
+ INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
GlobalVector res;
- LinAlg::matMult(A, x_new, res); // res = A * x_new
+ LinAlg::matMult(A, x_new, res); // res = A * x_new
LinAlg::axpy(res, -1.0, rhs); // res -= rhs
_convergence_criterion->checkResidual(res);
}
@@ -90,7 +95,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
if (postIterationCallback)
postIterationCallback(iteration, x_new);
- switch(sys.postIteration(x_new))
+ switch (sys.postIteration(x_new))
{
case IterationResult::SUCCESS:
// Don't copy here. The old x might still be used further
@@ -109,7 +114,8 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
INFO(
"Picard: The postIteration() hook decided that this "
"iteration has to be repeated.");
- continue; // That throws the iteration result away.
+ LinAlg::copy(x, x_new); // throw the iteration result away
+ continue;
}
}
@@ -184,7 +190,7 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
bool error_norms_met = false;
// TODO be more efficient
- // init _minus_delta_x to the right size and 0.0
+ // init minus_delta_x to the right size
LinAlg::copy(x, minus_delta_x);
_convergence_criterion->preFirstIteration();
@@ -193,9 +199,17 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
for (; iteration <= _maxiter;
++iteration, _convergence_criterion->reset())
{
+ BaseLib::RunTime timer_dirichlet;
+ double time_dirichlet = 0.0;
+
BaseLib::RunTime time_iteration;
time_iteration.start();
+ timer_dirichlet.start();
+ sys.computeKnownSolutions(x);
+ sys.applyKnownSolutions(x);
+ time_dirichlet += timer_dirichlet.elapsed();
+
sys.preIteration(iteration, x);
BaseLib::RunTime time_assembly;
@@ -207,10 +221,10 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
minus_delta_x.setZero();
- BaseLib::RunTime time_dirichlet;
- time_dirichlet.start();
- sys.applyKnownSolutionsNewton(J, res, minus_delta_x, x);
- INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
+ timer_dirichlet.start();
+ sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
+ time_dirichlet += timer_dirichlet.elapsed();
+ INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
_convergence_criterion->checkResidual(res);
@@ -255,8 +269,6 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
continue; // That throws the iteration result away.
}
- // TODO could be done via swap. Note: that also requires swapping
- // the ids. Same for the Picard scheme.
LinAlg::copy(x_new, x); // copy new solution to x
NumLib::GlobalVectorProvider::provider.releaseVector(x_new);
}
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index d43582d38da3d619e30f1d2a47009a2f5408f708..df992ffe87dd1ef25cd967804c40a51628cc9a1e 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -179,7 +179,7 @@ private:
std::size_t _A_id = 0u; //!< ID of the \f$ A \f$ matrix.
std::size_t _rhs_id = 0u; //!< ID of the right-hand side vector.
std::size_t _x_new_id = 0u; //!< ID of the vector storing the solution of
- //!the linearized equation.
+ //! the linearized equation.
};
/*! Creates a new nonlinear solver from the given configuration.
diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
index b3bf35c59a6a99be5ab17ef1d0f69226b62f9efc..8f78d33151712272f2fe718186e5774627ebbbee 100644
--- a/NumLib/ODESolver/NonlinearSystem.h
+++ b/NumLib/ODESolver/NonlinearSystem.h
@@ -53,14 +53,19 @@ public:
*/
virtual void getJacobian(GlobalMatrix& Jac) const = 0;
+ //! Pre-compute known solutions and possibly store them internally.
+ virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+
//! Apply known solutions to the solution vector \c x.
+ //! \pre computeKnownSolutions() must have been called before.
virtual void applyKnownSolutions(GlobalVector& x) const = 0;
//! Apply known solutions to the linearized equation system
//! \f$ \mathit{Jac} \cdot (-\Delta x) = \mathit{res} \f$.
- virtual void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x,
- GlobalVector& x) = 0;
+ //! \pre computeKnownSolutions() must have been called before.
+ virtual void applyKnownSolutionsNewton(
+ GlobalMatrix& Jac, GlobalVector& res,
+ GlobalVector& minus_delta_x) const = 0;
};
/*! A System of nonlinear equations to be solved with the Picard fixpoint
@@ -86,13 +91,18 @@ public:
//! \pre assemble() must have been called before.
virtual void getRhs(GlobalVector& rhs) const = 0;
+ //! Pre-compute known solutions and possibly store them internally.
+ virtual void computeKnownSolutions(GlobalVector const& x) = 0;
+
//! Apply known solutions to the solution vector \c x.
+ //! \pre computeKnownSolutions() must have been called before.
virtual void applyKnownSolutions(GlobalVector& x) const = 0;
//! Apply known solutions to the linearized equation system
//! \f$ A \cdot x = \mathit{rhs} \f$.
+ //! \pre computeKnownSolutions() must have been called before.
virtual void applyKnownSolutionsPicard(GlobalMatrix& A, GlobalVector& rhs,
- GlobalVector& x) = 0;
+ GlobalVector& x) const = 0;
};
//! @}
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 617c722a65fe4f41124d7b7712fcf457b35bb1e7..f2392d04a1ac77e26ec483a54fb857ae824fb4bc 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -120,28 +120,31 @@ void TimeDiscretizedODESystem<
_mat_trans->computeJacobian(*_Jac, Jac);
}
+void TimeDiscretizedODESystem<
+ ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Newton>::computeKnownSolutions(GlobalVector const& x)
+{
+ _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+}
+
void TimeDiscretizedODESystem<
ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::applyKnownSolutions(GlobalVector& x) const
{
- ::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
+ ::detail::applyKnownSolutions(_known_solutions, x);
}
void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
NonlinearSolverTag::Newton>::
applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x, GlobalVector& x)
+ GlobalVector& minus_delta_x) const
{
- auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
-
- if (!known_solutions || known_solutions->empty())
+ if (!_known_solutions || _known_solutions->empty())
return;
using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
std::vector<IndexType> ids;
- for (auto const& bc : *known_solutions)
+ for (auto const& bc : *_known_solutions)
{
std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
}
@@ -149,8 +152,6 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
// For the Newton method the values must be zero
std::vector<double> values(ids.size(), 0);
MathLib::applyKnownSolution(Jac, res, minus_delta_x, ids, values);
-
- ::detail::applyKnownSolutions(known_solutions, x);
}
TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
@@ -201,34 +202,37 @@ void TimeDiscretizedODESystem<
void TimeDiscretizedODESystem<
ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
+ NonlinearSolverTag::Picard>::computeKnownSolutions(GlobalVector const& x)
{
- ::detail::applyKnownSolutions(
- _ode.getKnownSolutions(_time_disc.getCurrentTime(), x), x);
+ _known_solutions = _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
}
void TimeDiscretizedODESystem<
ODESystemTag::FirstOrderImplicitQuasilinear,
- NonlinearSolverTag::Picard>::applyKnownSolutionsPicard(GlobalMatrix& A,
- GlobalVector& rhs,
- GlobalVector& x)
+ NonlinearSolverTag::Picard>::applyKnownSolutions(GlobalVector& x) const
{
- auto const* known_solutions =
- _ode.getKnownSolutions(_time_disc.getCurrentTime(), x);
+ ::detail::applyKnownSolutions(_known_solutions, x);
+}
+
+void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
+ NonlinearSolverTag::Picard>::
+ applyKnownSolutionsPicard(GlobalMatrix& A,
+ GlobalVector& rhs,
+ GlobalVector& x) const
+{
+ if (!_known_solutions || _known_solutions->empty())
+ return;
- if (known_solutions)
+ using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+ std::vector<IndexType> ids;
+ std::vector<double> values;
+ for (auto const& bc : *_known_solutions)
{
- using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
- std::vector<IndexType> ids;
- std::vector<double> values;
- for (auto const& bc : *known_solutions)
- {
- std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
- std::copy(bc.values.cbegin(), bc.values.cend(),
- std::back_inserter(values));
+ std::copy(bc.ids.cbegin(), bc.ids.cend(), std::back_inserter(ids));
+ std::copy(bc.values.cbegin(), bc.values.cend(),
+ std::back_inserter(values));
}
MathLib::applyKnownSolution(A, rhs, x, ids, values);
- }
}
} // namespace NumLib
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 59ead507782483e9247ee37c7fd2d91f13e20287..498bb8ef808e90072ffbaee6455ef1053d00d60c 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -89,11 +89,12 @@ public:
void getJacobian(GlobalMatrix& Jac) const override;
+ void computeKnownSolutions(GlobalVector const& x) override;
+
void applyKnownSolutions(GlobalVector& x) const override;
void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
- GlobalVector& minus_delta_x,
- GlobalVector& x) override;
+ GlobalVector& minus_delta_x) const override;
bool isLinear() const override
{
@@ -129,6 +130,10 @@ private:
//! the object used to compute the matrix/vector for the nonlinear solver
std::unique_ptr<MatTrans> _mat_trans;
+ using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+ std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
+ 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$.
@@ -185,10 +190,12 @@ public:
_mat_trans->computeRhs(*_M, *_K, *_b, rhs);
}
+ void computeKnownSolutions(GlobalVector const& x) override;
+
void applyKnownSolutions(GlobalVector& x) const override;
void applyKnownSolutionsPicard(GlobalMatrix& A, GlobalVector& rhs,
- GlobalVector& x) override;
+ GlobalVector& x) const override;
bool isLinear() const override
{
@@ -224,6 +231,10 @@ private:
//! the object used to compute the matrix/vector for the nonlinear solver
std::unique_ptr<MatTrans> _mat_trans;
+ using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+ std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
+ nullptr; //!< stores precomputed values for known solutions
+
GlobalMatrix* _M; //!< Matrix \f$ M \f$.
GlobalMatrix* _K; //!< Matrix \f$ K \f$.
GlobalVector* _b; //!< Matrix \f$ b \f$.