LGR: rebuild Cartesian->compressed map after refinement, enable serial summary output, fix rank-asymmetric throw

opm/simulators/flow/CpGridVanguard.hpp

@@ -27,6 +27,7 @@
 #ifndef OPM_CPGRID_VANGUARD_HPP
 #define OPM_CPGRID_VANGUARD_HPP
 
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/models/common/multiphasebaseproperties.hh>
@@ -43,6 +44,8 @@
 #include <tuple>
 #include <vector>
 
+#include <fmt/format.h>
+
 namespace Opm {
 template <class TypeTag>
 class CpGridVanguard;
@@ -118,14 +121,57 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
 
     int compressedIndexForInteriorLGR(const std::string& lgr_tag, const Connection& conn) const override
     {
-        const std::array<int,3> lgr_ijk = {conn.getI(), conn.getJ(), conn.getK()};
-        const auto& lgr_level = this->grid().getLgrNameToLevel().at(lgr_tag);
+        // Every rank registers every requested LGR name, with an empty level
+        // grid on ranks that hold no cell of the box (interior or overlap) --
+        // the level structure is identical on all ranks.  A name that fails to
+        // resolve is therefore a programming error, not a distribution effect,
+        // and must be fatal rather than silently skipped.
+        const auto& nameToLevel = this->grid().getLgrNameToLevel();
+        const auto levelIt = nameToLevel.find(lgr_tag);
+        if (levelIt == nameToLevel.end()) {
+            OPM_THROW(std::logic_error,
+                      fmt::format("Internal error: LGR '{}' is not known to the grid. "
+                                  "The level structure must be identical on all ranks.",
+                                  lgr_tag));
+        }
+        const int lgr_level = levelIt->second;
+
         if (ParentType::lgrMappers_.has_value() == false) {
             ParentType::lgrMappers_.emplace(this->grid().mapLocalCartesianIndexSetsToLeafIndexSet());
         }
+
+        // An out-of-range Cartesian position within the level is likewise a
+        // bug (a COMPDATL record addressing outside its LGR box has already
+        // been validated at parse time), so it is fatal too.
         const auto& lgr_dim = this->grid().currentData()[lgr_level]->logicalCartesianSize();
+        const std::array<int,3> lgr_ijk = {conn.getI(), conn.getJ(), conn.getK()};
+        if (lgr_ijk[0] < 0 || lgr_ijk[0] >= lgr_dim[0] ||
+            lgr_ijk[1] < 0 || lgr_ijk[1] >= lgr_dim[1] ||
+            lgr_ijk[2] < 0 || lgr_ijk[2] >= lgr_dim[2])
+        {
+            OPM_THROW(std::logic_error,
+                      fmt::format("Internal error: connection ({},{},{}) is outside "
+                                  "LGR '{}' with dimensions {}x{}x{}.",
+                                  lgr_ijk[0], lgr_ijk[1], lgr_ijk[2], lgr_tag,
+                                  lgr_dim[0], lgr_dim[1], lgr_dim[2]));
+        }
         const auto lgr_cartesian_index = (lgr_ijk[2]*lgr_dim[0]*lgr_dim[1]) + (lgr_ijk[1]*lgr_dim[0]) + (lgr_ijk[0]);
-        return ParentType::lgrMappers_.value()[lgr_level].at(lgr_cartesian_index);
+
+        // A cell that is absent from this rank's level mapper is the one
+        // legitimate miss: the box lives elsewhere and this rank's level grid
+        // is empty (or holds another part of it).  Mirror
+        // compressedIndexForInterior and return -1; the global existence of
+        // every connection cell is checked collectively afterwards
+        // (checkAllConnectionsFound).  Using .at() here threw
+        // std::out_of_range on such ranks -- asymmetrically, which deadlocked
+        // runs at higher rank counts where more ranks hold no part of the box.
+        const auto& mapper = ParentType::lgrMappers_.value()[lgr_level];
+        const auto it = mapper.find(lgr_cartesian_index);
+        if (it == mapper.end()) {
+            return -1;
+        }
+
+        return static_cast<int>(it->second);
     }
     /*!
      * Checking consistency of simulator
@@ -253,9 +299,7 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
                              this->numJacobiBlocks(), this->enableEclOutput());
 #endif
 
-        this->updateGridView_();
-        this->updateCartesianToCompressedMapping_();
-        this->updateCellThickness_();
+        this->updateDerivedGridState_();
 
 #if HAVE_MPI
         this->distributeFieldProps_(this->eclState());
@@ -266,6 +310,22 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
         this->updateCellDepths_();
     }
 
+    /*!
+     * \brief Recompute everything the vanguard derives from the grid.
+     *
+     * Needed after every change to the leaf grid -- load balancing and local
+     * refinement both renumber the leaf cells.  Kept in one method so a
+     * future grid-changing step cannot miss one of the updates.  Cell depths
+     * are not included: they need the distributed field properties, so each
+     * caller updates them once those are available.
+     */
+    void updateDerivedGridState_()
+    {
+        this->updateGridView_();
+        this->updateCartesianToCompressedMapping_();
+        this->updateCellThickness_();
+    }
+
     /*!
      * \brief Add LGRs and update Leaf Grid View in the simulation grid.
      */
@@ -277,9 +337,13 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
             OpmLog::info("\nAdding LGRs to the grid and updating its leaf grid view");
             this->addLgrsUpdateLeafView(lgrs, lgrs.size(), *this->grid_);
 
-            this->updateGridView_();
+            // Refinement changed the leaf cell count and ordering, so the
+            // state derived at load-balance time -- in particular the
+            // (level-zero-only) Cartesian->compressed map used to resolve
+            // coarse well connections -- is stale and must be rebuilt before
+            // well connections are resolved.
+            this->updateDerivedGridState_();
             this->updateCellDepths_();
-            this->updateCellThickness_();
 
             if (this->grid_->comm().size()>1) {
                 // Add LGRs and update the leaf grid view in the global (undistributed) simulation grid.

opm/simulators/flow/FlowBaseVanguard.hpp

@@ -338,12 +338,28 @@ class FlowBaseVanguard : public BaseVanguard<TypeTag>,
         std::size_t num_cells = asImp_().grid().leafGridView().size(0);
         is_interior_.resize(num_cells);
 
+        // May run again after the grid changed (local refinement), so start
+        // from scratch rather than leaving entries for cells that no longer
+        // exist in the leaf.
+        cartesianToCompressed_.clear();
+
         ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
         for (const auto& element : elements(this->gridView()))
         {
             const auto elemIdx = elemMapper.index(element);
-            unsigned cartesianCellIdx = cartesianIndex(elemIdx);
-            cartesianToCompressed_[cartesianCellIdx] = elemIdx;
+            // On a refined grid a level-zero Cartesian index is not a unique
+            // key: every child of a refined cell reports its ancestor's
+            // index, so inserting them would make the winner arbitrary.
+            // Only unrefined cells enter the map, making it a mapping for
+            // existing cells on level zero only; cells inside a refinement
+            // are addressed through the LGR-aware lookup instead, and a
+            // level-zero index that has been refined away resolves to
+            // "not present" rather than to an arbitrary child.
+            if (!element.hasFather())
+            {
+                unsigned cartesianCellIdx = cartesianIndex(elemIdx);
+                cartesianToCompressed_[cartesianCellIdx] = elemIdx;
+            }
             if (element.partitionType() == Dune::InteriorEntity)
             {
                 is_interior_[elemIdx] = 1;

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -515,12 +515,7 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // also updated.
         this->eclWriter().mutableOutputModule().invalidateLocalData();
 
-        // For CpGrid with LGRs, ecl/vtk output is not supported yet.
-        const auto& grid = this->simulator().vanguard().gridView().grid();
-
-        using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
-        constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
-        if (!isCpGrid || (grid.maxLevel() == 0)) {
+        if (this->eclOutputEvalSupported_()) {
             this->eclWriter_->evalSummaryState(!this->episodeWillBeOver());
         }
 
@@ -632,16 +627,10 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         // the initial solution.
         this->thresholdPressures_.finishInit();
 
-        // For CpGrid with LGRs, ecl-output is not supported yet.
-        const auto& grid = this->simulator().vanguard().gridView().grid();
-
-        using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
-        constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
-        // Skip - for now -  calculate the initial fip values for CpGrid with LGRs.
-        if (!isCpGrid || (grid.maxLevel() == 0)) {
-            if (this->simulator().episodeIndex() == 0) {
-                eclWriter_->writeInitialFIPReport();
-            }
+        if (this->eclOutputEvalSupported_() &&
+            (this->simulator().episodeIndex() == 0))
+        {
+            eclWriter_->writeInitialFIPReport();
         }
     }
 
@@ -1767,6 +1756,26 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
     HybridNewton hybridNewton_;
 
 private:
+    /// Whether ECL summary/FIP evaluation is wired up for this grid
+    /// configuration.
+    ///
+    /// The only unsupported case is a *distributed* CpGrid with LGRs:
+    /// CollectDataOnIORank does not build its index maps for a distributed
+    /// refined grid yet.  In serial the index maps are the identity and the
+    /// leaf-grid well/cell data is written directly, so refined serial runs
+    /// are fine.  Kept in one place so no call site can miss a condition.
+    bool eclOutputEvalSupported_() const
+    {
+        const auto& grid = this->simulator().vanguard().gridView().grid();
+
+        using GridType = std::remove_cv_t<std::remove_reference_t<decltype(grid)>>;
+        constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
+
+        return !isCpGrid
+            || (grid.maxLevel() == 0)
+            || (grid.comm().size() == 1);
+    }
+
     /// Whether or not the current epsiode will end at the end of the
     /// current time step.
     ///