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;
```
This commit is contained in:
Markus Blatt
2020-11-26 21:33:01 +01:00
parent 8ee58096ba
commit 3d92e41cad
2 changed files with 24 additions and 8 deletions

View File

@@ -104,23 +104,38 @@ namespace Opm {
} }
#endif #endif
B_->mv(x, y); B_->mv(x, y);
// The B matrix is basically a component-wise multiplication
// with a vector followed by a parallel reduction. We do that if (this->pinfo_->communication().size() > 1)
// reduction to all ranks computing for the well to save the {
// broadcast when applying C^T. // Only do communication if we must.
using YField = typename Y::block_type::value_type; // The B matrix is basically a component-wise multiplication
assert(y.size() == 1); // with a vector followed by a parallel reduction. We do that
this->pinfo_->communication().allreduce<std::plus<YField>>(y[0].container().data(), // reduction to all ranks computing for the well to save the
y[0].container().size()); // 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 //! y = A x
template<class X, class Y> template<class X, class Y>
void mmv (const X& x, Y& y) const 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); Y temp(y);
mv(x, temp); // includes parallel reduction mv(x, temp); // includes parallel reduction
y -= temp; y -= temp;
}
} }
private: private:
const Matrix* B_; const Matrix* B_;

View File

@@ -52,6 +52,7 @@ namespace Opm
, first_perf_(first_perf_index) , first_perf_(first_perf_index)
, perf_data_(&perf_data) , perf_data_(&perf_data)
{ {
assert(well.name()==pw_info.name());
assert(std::is_sorted(perf_data.begin(), perf_data.end(), assert(std::is_sorted(perf_data.begin(), perf_data.end(),
[](const auto& perf1, const auto& perf2){ [](const auto& perf1, const auto& perf2){
return perf1.ecl_index < perf2.ecl_index; return perf1.ecl_index < perf2.ecl_index;