Map region arrays onto the leaf grid for LGR (EQLNUM, PVTNUM, SWATINIT, FIP)

opm/simulators/flow/FlowGenericProblem_impl.hpp

@@ -544,6 +544,22 @@ readBlackoilExtentionsInitialConditions_(std::size_t numDof,
                                          bool enableBioeffects,
                                          bool enableMICP)
 {
+    // LGR (local grid refinement / CARFIN) is supported for black-oil only. The
+    // solvent/polymer/biofilm/MICP initial conditions below are read straight
+    // from the input-grid field properties and indexed by leaf cell; they are
+    // NOT mapped onto refined cells (see the black-oil EQLNUM/FIPNUM handling via
+    // LookUpData), so refined cells would get wrong/out-of-range values. Fail
+    // early with a clear message rather than producing silently wrong results.
+    if ((enableSolvent || enablePolymer || enablePolymerMolarWeight ||
+         enableBioeffects || enableMICP) &&
+        (eclState_.getLgrs().size() > 0))
+    {
+        throw std::runtime_error(
+            "Local grid refinement (LGR/CARFIN) is only supported for black-oil "
+            "runs. It is not supported together with the solvent, polymer, "
+            "biofilm or MICP extensions.");
+    }
+
     auto getArray = [](const std::vector<double>& input)
     {
         if constexpr (std::is_same_v<Scalar,double>) {

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -30,6 +30,7 @@
 #include <dune/common/fvector.hh>
 
 #include <opm/grid/CpGrid.hpp>
+#include <opm/grid/LookUpData.hh>
 
 #include <opm/simulators/utils/moduleVersion.hpp>
 
@@ -174,6 +175,12 @@ class OutputBlackOilModule : public GenericOutputModule<GetPropType<TypeTag, Pro
         , simulator_(simulator)
         , collectOnIORank_(collectOnIORank)
     {
+        // The region arrays arrive on the (unrefined) input grid, but everything
+        // downstream indexes them by leaf cell.  With an LGR the leaf has more
+        // cells, so map each onto it here - a refined cell inherits its parent's
+        // region - before anything looks at them.  Identity without LGRs.
+        this->mapRegionsOntoLeaf_();
+
         for (auto& region_pair : this->regions_) {
             this->createLocalRegion_(region_pair.second);
         }
@@ -226,15 +233,27 @@ class OutputBlackOilModule : public GenericOutputModule<GetPropType<TypeTag, Pro
             auto rset = this->eclState_.fieldProps().fip_regions();
             rset.push_back("PVTNUM");
 
-            // Note: We explicitly use decltype(auto) here because the
-            // default scheme (-> auto) will deduce an undesirable type.  We
-            // need the "reference to vector" semantics in this instance.
+            // RegionPhasePoreVolAverage indexes by leaf cell.  The FIP entries of
+            // regions_ are already on the leaf (mapRegionsOntoLeaf_ does that), so
+            // serve those directly.  PVTNUM is not an FIP region and is deliberately
+            // kept out of regions_, so it needs its own leaf-mapped copy - with LGRs
+            // a refined cell inherits its parent's PVTNUM; identity without.
+            // decltype(auto) keeps the required "reference to vector" semantics.
+            const LookUpData<Grid, GridView> lookUpData(this->simulator_.gridView());
+            auto pvtnum = lookUpData.template assignFieldPropsIntOnLeaf<int>(
+                this->eclState_.fieldProps(), "PVTNUM", /*needsTranslation=*/false);
+
             this->regionAvgDensity_
                 .emplace(this->simulator_.gridView().comm(),
                          FluidSystem::numPhases, rset,
-                         [fp = std::cref(this->eclState_.fieldProps())]
+                         [&regions = std::as_const(this->regions_),
+                          pvtnum = std::move(pvtnum)]
                          (const std::string& rsetName) -> decltype(auto)
-                         { return fp.get().get_int(rsetName); });
+                         {
+                             return (rsetName == "PVTNUM")
+                                 ? static_cast<const std::vector<int>&>(pvtnum)
+                                 : static_cast<const std::vector<int>&>(regions.at(rsetName));
+                         });
         }
     }
 
@@ -894,13 +913,28 @@ class OutputBlackOilModule : public GenericOutputModule<GetPropType<TypeTag, Pro
         }
     }
 
+    /// \brief Put every region array on the leaf grid.
+    ///
+    /// Called once, before the arrays are used.  A refined cell inherits its
+    /// parent's region; without LGRs this is the identity.
+    void mapRegionsOntoLeaf_()
+    {
+        const LookUpData<Grid, GridView> lookUpData(simulator_.gridView());
+        const auto numLeaf = simulator_.gridView().size(0);
+
+        std::vector<int> onLeaf(numLeaf, 0);
+        for (auto& [name, region] : this->regions_) {
+            std::size_t elemIdx = 0;
+            for (const auto& elem : elements(simulator_.gridView())) {
+                onLeaf[elemIdx++] = lookUpData(elem, region);
+            }
+
+            region = onLeaf;
+        }
+    }
+
     void createLocalRegion_(std::vector<int>& region)
     {
-        // For CpGrid with LGRs, where level zero grid has been distributed,
-        // resize region is needed, since in this case the total amount of
-        // element - per process - in level zero grid and leaf grid do not
-        // coincide, in general.
-        region.resize(simulator_.gridView().size(0));
         std::size_t elemIdx = 0;
         for (const auto& elem : elements(simulator_.gridView())) {
             if (elem.partitionType() != Dune::InteriorEntity) {

opm/simulators/flow/equil/InitStateEquil.hpp

@@ -745,7 +745,7 @@ class InitialStateComputer
                                  const bool co2store_or_h2store);
 
     template<class RMap>
-    void setRegionPvtIdx(const EclipseState& eclState, const RMap& reg);
+    void setRegionPvtIdx(const EclipseState& eclState, const GridView& gridView, const RMap& reg);
 
     template <class RMap, class MaterialLawManager, class Comm>
     void calcPressSatRsRv(const RMap& reg,

opm/simulators/flow/equil/InitStateEquil_impl.hpp

@@ -28,6 +28,7 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/grid/utility/RegionMapping.hpp>
+#include <opm/grid/LookUpData.hh>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Tables/PbvdTable.hpp>
@@ -1449,8 +1450,18 @@ equilnum(const EclipseState& eclipseState,
     std::vector<int> eqlnum(gridview.size(0), 0);
 
     if (eclipseState.fieldProps().has_int("EQLNUM")) {
-        const auto& e = eclipseState.fieldProps().get_int("EQLNUM");
-        std::ranges::transform(e, eqlnum.begin(), [](int n) { return n - 1; });
+        // EQLNUM is given on the (unrefined) input grid, but the equilibration
+        // works on the leaf grid. With LGRs the leaf has more cells than the
+        // input grid and a different ordering, so copying the input array
+        // directly into a leaf-sized vector misaligns it (and leaves refined
+        // cells at region 1). LookUpData maps each leaf cell to its input-grid
+        // origin - a refined cell inherits its parent cell's EQLNUM - which for
+        // an unrefined grid reduces to the identity, so non-LGR cases are
+        // unchanged. needsTranslation == true applies the 1-based -> 0-based
+        // shift previously done by the transform.
+        const LookUpData<typename GridView::Grid, GridView> lookUpData(gridview);
+        eqlnum = lookUpData.template assignFieldPropsIntOnLeaf<int>(
+            eclipseState.fieldProps(), "EQLNUM", /*needsTranslation=*/true);
     }
     OPM_BEGIN_PARALLEL_TRY_CATCH();
     const int num_regions = eclipseState.getTableManager().getEqldims().getNumEquilRegions();
@@ -1501,12 +1512,17 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
     //Check for presence of kw SWATINIT
     if (applySwatInit) {
         if (eclipseState.fieldProps().has_double("SWATINIT")) {
-            if constexpr (std::is_same_v<Scalar,double>) {
-                swatInit_ = eclipseState.fieldProps().get_double("SWATINIT");
+            // SWATINIT is given on the (unrefined) input grid but is consumed per
+            // leaf cell; with LGRs the leaf is larger and reordered, so map it
+            // onto the leaf via LookUpData (a refined cell inherits its parent
+            // cell's value; identity without LGRs).
+            const LookUpData<Grid, GridView> lookUpData(gridView);
+            auto input =
+                lookUpData.assignFieldPropsDoubleOnLeaf(eclipseState.fieldProps(), "SWATINIT");
+            if constexpr (std::is_same_v<Scalar, double>) {
+                swatInit_ = std::move(input);
             } else {
-                const auto& input = eclipseState.fieldProps().get_double("SWATINIT");
-                swatInit_.resize(input.size());
-                std::ranges::copy(input, swatInit_.begin());
+                swatInit_.assign(input.begin(), input.end());
             }
         }
     }
@@ -1520,10 +1536,10 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
     const std::vector<EquilRecord> rec = getEquil(eclipseState);
     const auto& tables = eclipseState.getTableManager();
     // Create (inverse) region mapping.
-    const RegionMapping<> eqlmap(equilnum(eclipseState, grid));
+    const RegionMapping<> eqlmap(equilnum(eclipseState, gridView));
     const int invalidRegion = -1;
     regionPvtIdx_.resize(rec.size(), invalidRegion);
-    setRegionPvtIdx(eclipseState, eqlmap);
+    setRegionPvtIdx(eclipseState, gridView, eqlmap);
 
     // Create Rs functions.
     rsFunc_.reserve(rec.size());
@@ -1979,13 +1995,22 @@ void InitialStateComputer<FluidSystem,
                           GridView,
                           ElementMapper,
                           CartesianIndexMapper>::
-setRegionPvtIdx(const EclipseState& eclState, const RMap& reg)
-{
-    const auto& pvtnumData = eclState.fieldProps().get_int("PVTNUM");
+setRegionPvtIdx(const EclipseState& eclState, const GridView& gridView, const RMap& reg)
+{
+    // PVTNUM is given on the (unrefined) input grid, but reg.cells(r) are leaf
+    // cell indices. With LGRs the leaf has more cells (and a different ordering)
+    // than the input grid, so indexing the input PVTNUM array by a leaf index is
+    // wrong (and out of bounds for refined cells). Map PVTNUM onto the leaf via
+    // LookUpData - a refined cell inherits its parent cell's PVTNUM - which is
+    // the identity without LGRs. needsTranslation == true applies the 1-based ->
+    // 0-based shift previously done explicitly.
+    const LookUpData<typename GridView::Grid, GridView> lookUpData(gridView);
+    const auto pvtnumData = lookUpData.template assignFieldPropsIntOnLeaf<int>(
+        eclState.fieldProps(), "PVTNUM", /*needsTranslation=*/true);
 
     for (const auto& r : reg.activeRegions()) {
         const auto& cells = reg.cells(r);
-        regionPvtIdx_[r] = pvtnumData[*cells.begin()] - 1;
+        regionPvtIdx_[r] = pvtnumData[*cells.begin()];
     }
 }