From a221cec7ff4c6071fd8ecbd44fb09d048416b15f Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Wed, 11 Nov 2020 07:43:50 +0100
Subject: [PATCH] [nonlinear] fix norm calculation in the Newton solver
Use the norm provided by the linear solver also for the initial
residual. In subsequent steps, the linear solver calculated the
norm of the solution instead of the residual. Fix this.
---
dumux/nonlinear/newtonsolver.hh | 21 +++++++++++++++------
1 file changed, 15 insertions(+), 6 deletions(-)
diff --git a/dumux/nonlinear/newtonsolver.hh b/dumux/nonlinear/newtonsolver.hh
index b33f113a46..647a0c6ba2 100644
--- a/dumux/nonlinear/newtonsolver.hh
+++ b/dumux/nonlinear/newtonsolver.hh
@@ -444,12 +444,18 @@ public:
{
if (numSteps_ == 0)
{
- Scalar norm2 = b.two_norm2();
- if (comm_.size() > 1)
- norm2 = comm_.sum(norm2);
+ if constexpr (Detail::hasNorm<LinearSolver, SolutionVector>())
+ initialResidual_ = this->linearSolver().norm(b);
- using std::sqrt;
- initialResidual_ = sqrt(norm2);
+ else
+ {
+ Scalar norm2 = b.two_norm2();
+ if (comm_.size() > 1)
+ norm2 = comm_.sum(norm2);
+
+ using std::sqrt;
+ initialResidual_ = sqrt(norm2);
+ }
}
// solve by calling the appropriate implementation depending on whether the linear solver
@@ -789,7 +795,10 @@ protected:
void computeResidualReduction_(const SolutionVector &uCurrentIter)
{
if constexpr (Detail::hasNorm<LinearSolver, SolutionVector>())
- residualNorm_ = this->linearSolver().norm(uCurrentIter);
+ {
+ this->assembler().assembleResidual(uCurrentIter);
+ residualNorm_ = this->linearSolver().norm(this->assembler().residual());
+ }
else
residualNorm_ = this->assembler().residualNorm(uCurrentIter);
--
GitLab