Allow empty partitions

ebos/eclbasevanguard.hh

@@ -97,6 +97,11 @@ struct OwnerCellsFirst {
     using type = UndefinedProperty;
 };
 
+template<class TypeTag, class MyTypeTag>
+struct AllowEmptyPartitions {
+    using type = UndefinedProperty;
+};
+
 template<class TypeTag, class MyTypeTag>
 struct SerialPartitioning {
     using type = UndefinedProperty;
@@ -158,6 +163,10 @@ struct OwnerCellsFirst<TypeTag, TTag::EclBaseVanguard> {
     static constexpr bool value = true;
 };
 template<class TypeTag>
+struct AllowEmptyPartitions<TypeTag, TTag::EclBaseVanguard> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
 struct SerialPartitioning<TypeTag, TTag::EclBaseVanguard> {
     static constexpr bool value = false;
 };
@@ -245,6 +254,8 @@ public:
 
         EWOMS_REGISTER_PARAM(TypeTag, bool, OwnerCellsFirst,
                              "Order cells owned by rank before ghost/overlap cells.");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, AllowEmptyPartitions,
+                             "Allow ranks with no grid cells.");
 #if HAVE_MPI
         EWOMS_REGISTER_PARAM(TypeTag, bool, SerialPartitioning,
                              "Perform partitioning for parallel runs on a single process.");
@@ -258,6 +269,7 @@ public:
                              "See https://sandialabs.github.io/Zoltan/ug_html/ug.html "
                              "for available Zoltan options.");
         EWOMS_HIDE_PARAM(TypeTag, ZoltanParams);
+        EWOMS_HIDE_PARAM(TypeTag, AllowEmptyPartitions);
 #endif
         EWOMS_REGISTER_PARAM(TypeTag, bool, AllowDistributedWells,
                              "Allow the perforations of a well to be distributed to interior of multiple processes");
@@ -283,6 +295,7 @@ public:
 #endif
 
         ownersFirst_ = EWOMS_GET_PARAM(TypeTag, bool, OwnerCellsFirst);
+        allowEmptyPartitions_ = EWOMS_GET_PARAM(TypeTag, bool, AllowEmptyPartitions);
 #if HAVE_MPI
         serialPartitioning_ = EWOMS_GET_PARAM(TypeTag, bool, SerialPartitioning);
         zoltanImbalanceTol_ = EWOMS_GET_PARAM(TypeTag, double, ZoltanImbalanceTol);
@@ -624,6 +637,9 @@ protected:
     /*! \brief Whether a cells is in the interior.
      */
     std::vector<int> is_interior_;
+
+    //! \brief Allow empty grid partitions.
+    bool allowEmptyPartitions_;
 };
 
 } // namespace Opm

ebos/eclcpgridvanguard.hh

@@ -101,6 +101,7 @@ public:
         : EclBaseVanguard<TypeTag>(simulator)
     {
         this->callImplementationInit();
+        this->grid().setAllowEmptyPartitions(this->allowEmptyPartitions_);
     }
 
     /*!

ebos/eclgenericoutputblackoilmodule.cc

@@ -20,6 +20,7 @@
 #include <config.h>
 
 #include <ebos/eclgenericoutputblackoilmodule.hh>
+#include <ebos/eclmpiserializer.hh>
 
 #include <opm/common/OpmLog/OpmLog.hpp>
 
@@ -642,7 +643,9 @@ addRftDataToWells(data::Wells& wellDatas, size_t reportStepNum)
 
 template<class FluidSystem, class Scalar>
 void EclGenericOutputBlackoilModule<FluidSystem,Scalar>::
-assignToSolution(data::Solution& sol)
+assignToSolution(data::Solution& sol,
+                 const Parallel::Communication& comm,
+                 const size_t num_cells)
 {
     using DataEntry = std::tuple<std::string,
                                  UnitSystem::measure,
@@ -709,49 +712,114 @@ assignToSolution(data::Solution& sol)
         {"WAT_VISC", UnitSystem::measure::viscosity, data::TargetType::RESTART_AUXILIARY,     viscosity_[waterPhaseIdx]},
     };
 
-    for (const auto& entry : data)
-        doInsert(entry);
+    // if we have processes with zero cells, we need to broadcast
+    // keys to include since local checks for zero vectors will cause
+    // differences across processes. we here assume process 0 owns cells
+    std::vector<std::size_t> process_cells(comm.size(), 1);
+    bool has_zeros = false;
+    int first_nonzero = 0;
+    if (comm.size() > 1) {
+        comm.allgather(&num_cells, 1, process_cells.data());
+        has_zeros = *std::min_element(process_cells.begin(), process_cells.end()) == 0;
+        first_nonzero = std::distance(process_cells.begin(),
+                            std::find_if(process_cells.begin(), process_cells.end(),
+                                         [](std::size_t r) { return r != 0; }));
+    }
+
+    using ExtraData = std::tuple<std::string,
+                                 UnitSystem::measure,
+                                 data::TargetType>;
+    std::vector<ExtraData> extras;
+    auto insertExtra = [has_zeros, first_nonzero,
+                        &extras,& sol, comm](const std::string& name,
+                                                        UnitSystem::measure measure,
+                                                        std::vector<Scalar>&& vec,
+                                                        data::TargetType type)
+    {
+        if (has_zeros && comm.rank() == first_nonzero)
+            extras.emplace_back(std::make_tuple(name, measure, type));
+        sol.insert(name, measure, vec, type);
 
-    if (!temperature_.empty()) {
-        if (enableEnergy_)
-            sol.insert("TEMP", UnitSystem::measure::temperature, std::move(temperature_), data::TargetType::RESTART_SOLUTION);
-        else {
-            // Flow allows for initializing of non-constant initial temperature.
-            // For output of this temperature for visualization and restart set --enable-opm-restart=true
-            assert(enableTemperature_);
-            sol.insert("TEMP", UnitSystem::measure::temperature, std::move(temperature_), data::TargetType::RESTART_AUXILIARY);
+    };
+    if (process_cells[comm.rank()] > 0) {
+        for (const auto& entry : data)
+            doInsert(entry);
+
+        if (!temperature_.empty()) {
+            if (enableEnergy_)
+                insertExtra("TEMP", UnitSystem::measure::temperature, std::move(temperature_), data::TargetType::RESTART_SOLUTION);
+            else {
+                // Flow allows for initializing of non-constant initial temperature.
+                // For output of this temperature for visualization and restart set --enable-opm-restart=true
+                assert(enableTemperature_);
+                insertExtra("TEMP", UnitSystem::measure::temperature, std::move(temperature_), data::TargetType::RESTART_AUXILIARY);
+            }
         }
-    }
 
-    if (FluidSystem::phaseIsActive(waterPhaseIdx) && !saturation_[waterPhaseIdx].empty()) {
-        sol.insert("SWAT", UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), data::TargetType::RESTART_SOLUTION);
-    }
-    if (FluidSystem::phaseIsActive(gasPhaseIdx) && !saturation_[gasPhaseIdx].empty()) {
-        sol.insert("SGAS", UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), data::TargetType::RESTART_SOLUTION);
-    }
+        if (FluidSystem::phaseIsActive(waterPhaseIdx) && !saturation_[waterPhaseIdx].empty()) {
+            insertExtra("SWAT", UnitSystem::measure::identity, std::move(saturation_[waterPhaseIdx]), data::TargetType::RESTART_SOLUTION);
+        }
+        if (FluidSystem::phaseIsActive(gasPhaseIdx) && !saturation_[gasPhaseIdx].empty()) {
+            insertExtra("SGAS", UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), data::TargetType::RESTART_SOLUTION);
+        }
 
-    // Fluid in place
-    for (const auto& phase : Inplace::phases()) {
-        if (outputFipRestart_ && !fip_[phase].empty()) {
-            sol.insert(EclString(phase),
-                       UnitSystem::measure::volume,
-                       fip_[phase],
-                       data::TargetType::SUMMARY);
+        // Fluid in place
+        for (const auto& phase : Inplace::phases()) {
+            if (outputFipRestart_ && !fip_[phase].empty()) {
+                ScalarBuffer f = fip_[phase];
+                insertExtra(EclString(phase),
+                            UnitSystem::measure::volume,
+                            std::move(f),
+                            data::TargetType::SUMMARY);
+            }
         }
-    }
 
-    // tracers
-    if (!tracerConcentrations_.empty()) {
-        const auto& tracers = eclState_.tracer();
-        for (std::size_t tracerIdx = 0; tracerIdx < tracers.size(); tracerIdx++) {
-            const auto& tracer = tracers[tracerIdx];
-            sol.insert(tracer.fname(), UnitSystem::measure::identity, std::move(tracerConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION);
+        // tracers
+        if (!tracerConcentrations_.empty()) {
+            const auto& tracers = eclState_.tracer();
+            for (std::size_t tracerIdx = 0; tracerIdx < tracers.size(); tracerIdx++) {
+                const auto& tracer = tracers[tracerIdx];
+                insertExtra(tracer.fname(), UnitSystem::measure::identity, std::move(tracerConcentrations_[tracerIdx]), data::TargetType::RESTART_TRACER_SOLUTION);
+            }
+            // We need put tracerConcentrations into a valid state.
+            // Otherwise next time we end up here outside of a restart write we will
+            // move invalidated data above (as it was moved away before and never
+            // reallocated)
+            tracerConcentrations_.resize(0);
+        }
+    }
+    std::vector<std::string> keys;
+    if (has_zeros) {
+        if (comm.rank() == 0) {
+            for (const auto& entry : sol) {
+                keys.push_back(entry.first);
+            }
+        }
+        EclMpiSerializer serializer(comm);
+        serializer.broadcast(first_nonzero, keys, extras);
+
+        if (process_cells[comm.rank()] == 0) {
+            for (const auto& key : keys) {
+                const auto entry =
+                    std::find_if(data.begin(), data.end(),
+                                [key](const DataEntry& e)
+                                {
+                                    return std::get<0>(e) == key;
+                                });
+                if (entry != data.end()) {
+                    sol.insert(key, std::get<1>(*entry),
+                               std::move(std::get<3>(*entry)), std::get<2>(*entry));
+                } else {
+                    const auto entry2  =
+                        std::find_if(extras.begin(), extras.end(),
+                                    [key](const ExtraData& e)
+                                    {
+                                        return std::get<0>(e) == key;
+                                    });
+                    sol.insert(key, std::get<1>(*entry2), {}, std::get<2>(*entry2));
+                }
+            }
         }
-        // We need put tracerConcentrations into a valid state.
-        // Otherwise next time we end up here outside of a restart write we will
-        // move invalidated data above (as it was moved away before and never
-        // reallocated)
-        tracerConcentrations_.resize(0);
     }
 }
 
@@ -839,11 +907,10 @@ regionSum(const ScalarBuffer& property,
 {
         ScalarBuffer totals(maxNumberOfRegions, 0.0);
 
-        if (property.empty())
+        if (property.empty() && comm.size() == 1)
             return totals;
 
-        assert(regionId.size() == property.size());
-        for (size_t j = 0; j < regionId.size(); ++j) {
+        for (size_t j = 0; j < regionId.size() && !property.empty(); ++j) {
             const int regionIdx = regionId[j] - 1;
             // the cell is not attributed to any region. ignore it!
             if (regionIdx < 0)

ebos/eclgenericoutputblackoilmodule.hh

@@ -92,7 +92,9 @@ public:
     /*!
      * \brief Move all buffers to data::Solution.
      */
-    void assignToSolution(data::Solution& sol);
+    void assignToSolution(data::Solution& sol,
+                          const Parallel::Communication& comm,
+                          const size_t num_cells);
 
     void setRestart(const data::Solution& sol,
                     unsigned elemIdx,

ebos/eclwriter.hh

@@ -334,8 +334,9 @@ public:
 
         data::Solution localCellData = {};
         if (! isSubStep) {
-            this->eclOutputModule_->assignToSolution(localCellData);
-
+            this->eclOutputModule_->assignToSolution(localCellData,
+                                                     simulator_.gridView().comm(),
+                                                     simulator_.gridView().size(/*codim=*/0));
             // add cell data to perforations for Rft output
             this->eclOutputModule_->addRftDataToWells(localWellData, reportStepNum);
         }

opm/simulators/linalg/PressureBhpTransferPolicy.hpp

@@ -50,7 +50,6 @@ namespace Opm
         }
         const int global_max = comm.communicator().max(glo_max);
         // used append the welldofs at the end
-        assert(loc_max + 1 == indset.size());
         size_t local_ind = loc_max + 1;
         for (int i = 0; i < nw; ++i) {
             // need to set unique global number
@@ -109,13 +108,20 @@ namespace Opm
             const auto& nw = fineOperator.getNumberOfExtraEquations();
             if (prm_.get<bool>("add_wells")) {
                 const size_t average_elements_per_row
-                    = static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() / fineLevelMatrix.N()));
+                    = static_cast<size_t>(std::ceil(fineLevelMatrix.nonzeroes() /
+                                                    std::max(1ul,fineLevelMatrix.N())));
                 const double overflow_fraction = 1.2;
-                coarseLevelMatrix_.reset(new CoarseMatrix(fineLevelMatrix.N() + nw,
-                                                          fineLevelMatrix.M() + nw,
-                                                          average_elements_per_row,
-                                                          overflow_fraction,
-                                                          CoarseMatrix::implicit));
+                if (fineLevelMatrix.N() + nw > 0) {
+                    coarseLevelMatrix_ = std::make_unique<CoarseMatrix>(fineLevelMatrix.N() + nw,
+                                                                        fineLevelMatrix.M() + nw,
+                                                                        average_elements_per_row,
+                                                                        overflow_fraction,
+                                                                        CoarseMatrix::implicit);
+                } else { // workaround bcrsmatrix bug with implicit build mode + empty matrix
+                    coarseLevelMatrix_ = std::make_unique<CoarseMatrix>(fineLevelMatrix.N() + nw,
+                                                                        fineLevelMatrix.M() + nw,
+                                                                        CoarseMatrix::random);
+                }
                 int rownum = 0;
                 for (const auto& row : fineLevelMatrix) {
                     for (auto col = row.begin(), cend = row.end(); col != cend; ++col) {
@@ -151,7 +157,12 @@ namespace Opm
 #endif
         if (prm_.get<bool>("add_wells")) {
             fineOperator.addWellPressureEquationsStruct(*coarseLevelMatrix_);
-            coarseLevelMatrix_->compress(); // all elemenst should be set
+            if (coarseLevelMatrix_->N() > 0) {
+                coarseLevelMatrix_->compress(); // all elements should be set
+            } else { // workaround for bug in compress()
+                coarseLevelMatrix_->endrowsizes();
+                coarseLevelMatrix_->endindices();
+            }
             if constexpr (!std::is_same_v<Communication, Dune::Amg::SequentialInformation>) {
                 extendCommunicatorWithWells(*communication_, coarseLevelCommunication_, nw);
             }
@@ -172,6 +183,9 @@ namespace Opm
 
     virtual void calculateCoarseEntries(const FineOperator& fineOperator) override
     {
+        if (coarseLevelMatrix_->N() == 0)
+            return;
+
         const auto& fineMatrix = fineOperator.getmat();
         *coarseLevelMatrix_ = 0;
         auto rowCoarse = coarseLevelMatrix_->begin();

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -441,8 +441,8 @@ checkGroupHigherConstraints(const Group& group,
     // What is the proper approach?
     const int fipnum = 0;
     int pvtreg = well_perf_data_.empty() || well_perf_data_[0].empty()
-        ? pvt_region_idx_[0]
-        : pvt_region_idx_[well_perf_data_[0][0].cell_index];
+        ? (pvt_region_idx_.empty() ? -1 : pvt_region_idx_[0])
+        : (pvt_region_idx_.empty() ? -1 : pvt_region_idx_[well_perf_data_[0][0].cell_index]);
 
     bool changed = false;
     if ( comm_.size() > 1)
@@ -452,8 +452,10 @@ checkGroupHigherConstraints(const Group& group,
         // is decided by the order in the Schedule using Well::seqIndex()
         int firstWellIndex = well_perf_data_.empty() ?
             std::numeric_limits<int>::max() : wells_ecl_[0].seqIndex();
-        auto regIndexPair = std::make_pair(pvtreg, firstWellIndex);
-        std::vector<decltype(regIndexPair)> pairs(comm_.size());
+        std::vector<std::pair<int,int>> pairs(comm_.size());
+        std::pair<int,int> regIndexPair{-1,std::numeric_limits<int>::max()};
+        if (pvtreg > -1)
+            regIndexPair = std::make_pair(pvtreg, firstWellIndex);
         comm_.allgather(&regIndexPair, 1, pairs.data());
         pvtreg = std::min_element(pairs.begin(), pairs.end(),
                                   [](const auto& p1, const auto& p2){ return p1.second < p2.second;})

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -214,8 +214,9 @@ namespace Opm {
             WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState());
 
             // Compute reservoir volumes for RESV controls.
-            rateConverter_.reset(new RateConverterType (phase_usage_,
-                                                        std::vector<int>(local_num_cells_, 0)));
+            rateConverter_ = std::make_unique<RateConverterType>(phase_usage_,
+                                                                 std::vector<int>(local_num_cells_, 0),
+                                                                 ebosSimulator_.gridView().comm());
             rateConverter_->template defineState<ElementContext>(ebosSimulator_);
 
             // Compute regional average pressures used by gpmaint