From c3edb0f7edf9f2e94ef933011675f1215df6743f Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Mon, 6 Apr 2020 17:52:36 +0200
Subject: [PATCH] [multidomain][box][bugfix] Fix residual norm calculation for
parallel+nonoverlapping+box
---
dumux/multidomain/fvassembler.hh | 42 +++++++++++++++++++++++++++++---
1 file changed, 38 insertions(+), 4 deletions(-)
diff --git a/dumux/multidomain/fvassembler.hh b/dumux/multidomain/fvassembler.hh
index 401ad03d74..f2141beed7 100644
--- a/dumux/multidomain/fvassembler.hh
+++ b/dumux/multidomain/fvassembler.hh
@@ -220,17 +220,51 @@ public:
}
//! compute the residual and return it's vector norm
- //! TODO: this needs to be adapted in parallel
Scalar residualNorm(const SolutionVector& curSol)
{
ResidualType residual;
setResidualSize(residual);
assembleResidual(residual, curSol);
- // calculate the square norm of the residual
- const Scalar result2 = residual.two_norm2();
+ // for box communicate the residual with the neighboring processes
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(residual)), [&](const auto domainId)
+ {
+ const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
+ if (GridGeometry<domainId>::discMethod == DiscretizationMethod::box
+ && gridGeometry->gridView().comm().size() > 1)
+ {
+ using VertexMapper = typename GridGeometry<domainId>::VertexMapper;
+ using SolutionVector = std::remove_reference_t<decltype(residual[domainId])>;
+ using DataHandle = VectorCommDataHandleSum<VertexMapper, SolutionVector,
+ GridGeometry<domainId>::GridView::dimension>;
+ DataHandle sumResidualHandle(gridGeometry->vertexMapper(), residual[domainId]);
+ gridGeometry->gridView().communicate(sumResidualHandle,
+ Dune::InteriorBorder_InteriorBorder_Interface,
+ Dune::ForwardCommunication);
+ }
+ });
+
+ // calculate the squared norm of the residual
+ Scalar resultSquared = 0.0;
+
+ // add up the result across processes
+ using namespace Dune::Hybrid;
+ forEach(integralRange(Dune::Hybrid::size(residual)), [&](const auto domainId)
+ {
+ Scalar localNormSquared = residual[domainId].two_norm2();
+
+ const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
+ if (gridGeometry->gridView().comm().size() > 1)
+ {
+ localNormSquared = gridGeometry->gridView().comm().sum(localNormSquared);
+ }
+
+ resultSquared += localNormSquared;
+ });
+
using std::sqrt;
- return sqrt(result2);
+ return sqrt(resultSquared);
}
/*!
--
GitLab