From 36d8f49b008f1440e01aa8e7fabf0b03e31365fe Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Mon, 2 Sep 2019 00:03:14 +0200
Subject: [PATCH] [NL] Out-of-balance forces. Picard.
---
NumLib/ODESolver/NonlinearSolver.cpp | 23 +++++++++++++++++++++++
NumLib/ODESolver/NonlinearSolver.h | 10 +++-------
2 files changed, 26 insertions(+), 7 deletions(-)
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index bf258e3d263..59f406cf3b0 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -27,6 +27,24 @@ void NonlinearSolver<NonlinearSolverTag::Picard>::assemble(
_equation_system->assemble(x, process_id);
}
+void NonlinearSolver<NonlinearSolverTag::Picard>::calculateOutOfBalanceForces(
+ std::vector<GlobalVector*> const& x, int const process_id)
+{
+ if (!_compensate_initial_out_of_balance_forces)
+ return;
+
+ auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
+ auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
+ _equation_system->assemble(x, process_id);
+ _equation_system->getA(A);
+ _equation_system->getRhs(rhs);
+
+ // F_oob = A * x - rhs
+ _f_oob = &NumLib::GlobalVectorProvider::provider.getVector();
+ MathLib::LinAlg::matMult(A, *x[process_id], *_f_oob);
+ MathLib::LinAlg::axpy(*_f_oob, -1.0, rhs); // res -= rhs
+}
+
NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
std::vector<GlobalVector*>& x,
std::function<void(int, GlobalVector const&)> const& postIterationCallback,
@@ -73,6 +91,11 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Picard>::solve(
sys.getRhs(rhs);
INFO("[time] Assembly took %g s.", time_assembly.elapsed());
+
+ // Apply out-of-balance forces if set.
+ if (_f_oob != nullptr)
+ LinAlg::axpy(rhs, -1, *_f_oob);
+
timer_dirichlet.start();
sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
time_dirichlet += timer_dirichlet.elapsed();
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 04c5815372a..dcb1b5c7543 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -193,13 +193,8 @@ public:
void assemble(std::vector<GlobalVector*> const& x,
int const process_id) const override;
- void calculateOutOfBalanceForces(GlobalVector const& /*x*/) override
- {
- if (!_compensate_initial_out_of_balance_forces)
- return;
-
- OGS_FATAL("Calculation of out-of-balance forces is not implemented.");
- }
+ void calculateOutOfBalanceForces(std::vector<GlobalVector*> const& x,
+ int const process_id) override;
NonlinearSolverStatus solve(
std::vector<GlobalVector*>& x,
@@ -220,6 +215,7 @@ private:
ConvergenceCriterion* _convergence_criterion = nullptr;
const int _maxiter; //!< maximum number of iterations
+ GlobalVector* _f_oob = nullptr; //!< Out-of-balance forces vector.
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
--
GitLab