From 0739bb2c16a5918ac1e878d765bbd0b0724fe3e4 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Thu, 8 Feb 2024 16:48:47 +0100
Subject: [PATCH] [NL] Change Dirichlet procedure for Newton
Instead of setting the solution in the Dirichlet nodes, the solution is
enforced through the r.h.s. This usually saves one global iteration and
more importantly avoids big jumps in the first iteration for mechanics
processes.
---
NumLib/ODESolver/NonlinearSolver.cpp | 1 -
NumLib/ODESolver/TimeDiscretizedODESystem.cpp | 22 ++++++++++++++++---
2 files changed, 19 insertions(+), 4 deletions(-)
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 59aac171411..324c11905c6 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -345,7 +345,6 @@ NonlinearSolverStatus NonlinearSolver<NonlinearSolverTag::Newton>::solve(
timer_dirichlet.start();
sys.computeKnownSolutions(*x[process_id], process_id);
- sys.applyKnownSolutions(*x[process_id]);
time_dirichlet += timer_dirichlet.elapsed();
sys.preIteration(iteration, *x[process_id]);
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
index 08d086c07e8..f45e91ffbef 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.cpp
@@ -10,6 +10,9 @@
#include "TimeDiscretizedODESystem.h"
+#include <range/v3/numeric/accumulate.hpp>
+#include <range/v3/view/transform.hpp>
+
#include "MathLib/LinAlg/ApplyKnownSolution.h"
#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
#include "NumLib/Exceptions.h"
@@ -142,14 +145,27 @@ void TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
}
using IndexType = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
+ std::size_t const size = ranges::accumulate(
+ *_known_solutions | ranges::views::transform([](auto const& bc)
+ { return bc.ids.size(); }),
+ 0);
std::vector<IndexType> ids;
+ ids.reserve(size);
+ std::vector<double> values;
+ values.reserve(size);
+
for (auto const& bc : *_known_solutions)
{
- ids.insert(end(ids), begin(bc.ids), end(bc.ids));
+ for (std::size_t i = 0; i < bc.ids.size(); ++i)
+ {
+ auto const id = bc.ids[i];
+ ids.push_back(id);
+ // minus_delta_x will be set to the difference between the current
+ // value and the Dirichlet BC value.
+ values.push_back(x[id] - bc.values[i]);
+ }
}
- // 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);
}
--
GitLab