From b72a167c762a71d7cd9d49dc06462e9d490422c5 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 12 Apr 2017 13:12:36 +0200 Subject: [PATCH 1/4] Correctly compute maximum in a parallel flow_ebos run. --- opm/autodiff/BlackoilModelEbos.hpp | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index eb00c14d2..83fe624ed 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -791,6 +791,8 @@ namespace Opm { #if HAVE_MPI if ( comm.size() > 1 ) { + // Initialize maxCoeff with minimum. + std::fill(maxCoeff.begin(), maxCoeff.end(), std::numeric_limits::lowest()); // mask[c] is 1 if we need to compute something in parallel const auto & pinfo = boost::any_cast(istlSolver().parallelInformation()); @@ -803,7 +805,22 @@ namespace Opm { { B_avg[idx] = accumulateMaskedValues(B[ idx ], mask) / double(ncGlobal); R_sum[idx] = accumulateMaskedValues(R[ idx ], mask); - maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); + + if(comm.size()>1) + { + auto mi = mask->begin(); + for(auto elem = tempV[idx].begin(), end = tempV[idx].end(); elem != end; ++elem, ++mi) + { + if ( *mi ) + { + maxCoeff[idx] = std::max( maxCoeff[idx], *elem); + } + } + } + else + { + maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); + } assert(np >= np); if (idx < np) { From 0db663fe51555f95e3b0b7944a8a5e4ba92bd0fd Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 12 Apr 2017 13:27:24 +0200 Subject: [PATCH 2/4] Only compute convergence markers for interior elements. --- opm/autodiff/BlackoilModelEbos.hpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 83fe624ed..67f0cda41 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -901,18 +901,22 @@ namespace Opm { Vector maxCoeff(np); Vector maxNormWell(np); - std::vector< Vector > B( np, Vector( nc ) ); - std::vector< Vector > R( np, Vector( nc ) ); - std::vector< Vector > R2( np, Vector( nc ) ); - std::vector< Vector > tempV( np, Vector( nc ) ); + // As we will not initialize values in the following arrays + // for the non-interior elements, we have to make sure + // (at least for tempV) that the values there do not influence + // our reduction. + std::vector< Vector > B( np, Vector( nc, 0.0) ); + //std::vector< Vector > R( np, Vector( nc, 0.0) ) ); + std::vector< Vector > R2( np, Vector( nc, 0.0 ) ); + std::vector< Vector > tempV( np, Vector( nc, std::numeric_limits::lowest() ) ); const auto& ebosResid = ebosSimulator_.model().linearizer().residual(); ElementContext elemCtx(ebosSimulator_); const auto& gridView = ebosSimulator().gridView(); - const auto& elemEndIt = gridView.template end(); + const auto& elemEndIt = gridView.template end(); - for (auto elemIt = gridView.template begin(); + for (auto elemIt = gridView.template begin(); elemIt != elemEndIt; ++elemIt) { From b5612806ac467dda67d2eb27b2982b2e2de1aa40 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 12 Apr 2017 13:34:13 +0200 Subject: [PATCH 3/4] Revert to using standard algorithms instead of using masks. This is now possible as the values stored for ghost/overlap elements (minimum where we compute the maxiumum, zero where we sum up) will not influence the result of the computation any more. --- opm/autodiff/BlackoilModelEbos.hpp | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp index 67f0cda41..7e3c1ad7e 100644 --- a/opm/autodiff/BlackoilModelEbos.hpp +++ b/opm/autodiff/BlackoilModelEbos.hpp @@ -803,24 +803,12 @@ namespace Opm { // computation for ( int idx = 0; idx < np; ++idx ) { - B_avg[idx] = accumulateMaskedValues(B[ idx ], mask) / double(ncGlobal); - R_sum[idx] = accumulateMaskedValues(R[ idx ], mask); + B_avg[idx] = std::accumulate( B[ idx ].begin(), B[ idx ].end(), + 0.0 ) / double(ncGlobal); + R_sum[idx] = std::accumulate( R[ idx ].begin(), R[ idx ].end(), + 0.0 ); - if(comm.size()>1) - { - auto mi = mask->begin(); - for(auto elem = tempV[idx].begin(), end = tempV[idx].end(); elem != end; ++elem, ++mi) - { - if ( *mi ) - { - maxCoeff[idx] = std::max( maxCoeff[idx], *elem); - } - } - } - else - { - maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); - } + maxCoeff[idx] = *(std::max_element( tempV[ idx ].begin(), tempV[ idx ].end() )); assert(np >= np); if (idx < np) { From 9f5a904382b70ee771c400fcc87cb25950f781ec Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Wed, 12 Apr 2017 17:37:34 +0200 Subject: [PATCH 4/4] Removes deadlocks in the case where only few processes have wells. The problem was that updateWellControls was not called on all processes. But this is mandatory as the well switching output requires global communication. --- opm/autodiff/StandardWellsDense_impl.hpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index 7cc7fb04e..330ff6476 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -135,7 +135,7 @@ namespace Opm { } SimulatorReport report; - if ( ! localWellsActive() ) { + if ( ! wellsActive() ) { return report; } @@ -891,6 +891,12 @@ namespace Opm { invDuneD_.mv(resWell_, dx_well); updateWellState(dx_well, well_state); + } + // updateWellControls uses communication + // Therefore the following is executed if there + // are active wells anywhere in the global domain. + if( wellsActive() ) + { updateWellControls(well_state); setWellVariables(well_state); } @@ -1483,7 +1489,11 @@ namespace Opm { StandardWellsDense:: updateWellControls(WellState& xw) const { - if( !localWellsActive() ) return ; + // Even if there no wells active locally, we cannot + // return as the Destructor of the WellSwitchingLogger + // uses global communication. For no well active globally + // we simply return. + if( !wellsActive() ) return ; const int np = wells().number_of_phases; const int nw = wells().number_of_wells;