From 3d92e41cadfa64966529757c006521098f1867c5 Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Thu, 26 Nov 2020 21:33:01 +0100
Subject: [PATCH] Recover prev. iteration count and curves for undistributed
 wells

Rounding errors for `B.mmv(x,y)` are slightly different from
```
Y temp(y);
B-mv(x, temp);
y -= temp;
```
---
 opm/simulators/wells/WellHelpers.hpp        | 31 +++++++++++++++------
 opm/simulators/wells/WellInterface_impl.hpp |  1 +
 2 files changed, 24 insertions(+), 8 deletions(-)

diff --git a/opm/simulators/wells/WellHelpers.hpp b/opm/simulators/wells/WellHelpers.hpp
index c0ae7edb3..53d15793f 100644
--- a/opm/simulators/wells/WellHelpers.hpp
+++ b/opm/simulators/wells/WellHelpers.hpp
@@ -104,23 +104,38 @@ namespace Opm {
                 }
 #endif
                 B_->mv(x, y);
-                // The B matrix is basically a component-wise multiplication
-                // with a vector followed by a parallel reduction. We do that
-                // reduction to all ranks computing for the well to save the
-                // broadcast when applying C^T.
-                using YField = typename Y::block_type::value_type;
-                assert(y.size() == 1);
-                this->pinfo_->communication().allreduce<std::plus<YField>>(y[0].container().data(),
-                                                                           y[0].container().size());
+
+                if (this->pinfo_->communication().size() > 1)
+                {
+                    // Only do communication if we must.
+                    // The B matrix is basically a component-wise multiplication
+                    // with a vector followed by a parallel reduction. We do that
+                    // reduction to all ranks computing for the well to save the
+                    // broadcast when applying C^T.
+                    using YField = typename Y::block_type::value_type;
+                    assert(y.size() == 1);
+                    this->pinfo_->communication().allreduce<std::plus<YField>>(y[0].container().data(),
+                                                                               y[0].container().size());
+                }
             }
 
             //! y = A x
             template<class X, class Y>
             void mmv (const X& x, Y& y) const
             {
+                if (this->pinfo_->communication().size() == 1)
+                {
+                    // Do the same thing as before. The else branch
+                    // produces different rounding errors and results
+                    // slightly different iteration counts / well curves
+                    B_->mmv(x, y);
+                }
+                else
+                {
                 Y temp(y);
                 mv(x, temp); // includes parallel reduction
                 y -= temp;
+                }
             }
         private:
             const Matrix* B_;
diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp
index adabca396..ce4343d5e 100644
--- a/opm/simulators/wells/WellInterface_impl.hpp
+++ b/opm/simulators/wells/WellInterface_impl.hpp
@@ -52,6 +52,7 @@ namespace Opm
       , first_perf_(first_perf_index)
       , perf_data_(&perf_data)
     {
+        assert(well.name()==pw_info.name());
         assert(std::is_sorted(perf_data.begin(), perf_data.end(),
                               [](const auto& perf1, const auto& perf2){
                                   return perf1.ecl_index < perf2.ecl_index;