Merge pull request #4933 from aritorto/lookupElemIndices

LookUp origin cell indices in eclgenericproblem

ebos/eclgenericproblem.cc

@@ -29,6 +29,7 @@
 #include <dune/grid/common/gridview.hh>
 
 #include <opm/grid/CpGrid.hpp>
+#include <opm/grid/LookUpData.hh>
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 

ebos/eclgenericproblem.hh

@@ -50,6 +50,7 @@ namespace Opm {
 class Deck;
 class EclipseState;
 class Schedule;
+template<typename Grid, typename GridView> class LookUpData;
 
 int eclPositionalParameter(Dune::ParameterTree& tree,
                            std::set<std::string>& seenParams,
@@ -371,6 +372,17 @@ protected:
     // equilibration parameters
     int numPressurePointsEquil_;
 
+    // To lookup origin cell indices
+    using Grid = std::remove_cv_t< typename std::remove_reference<decltype(gridView_.grid())>::type>;
+    using LookUpData = Opm::LookUpData<Grid,GridView>;
+    const LookUpData lookUpData_;
+
+    auto getLookUpData(unsigned elemIdx) const
+    {
+        using GridType = std::remove_cv_t< typename std::remove_reference<decltype(gridView_.grid())>::type>;
+        return lookUpData_.template getOriginIndex<GridType>(elemIdx);
+    }
+
 private:
     template<class T>
     void updateNum(const std::string& name, std::vector<T>& numbers, std::size_t num_regions);

ebos/eclgenericproblem_impl.hh

@@ -89,6 +89,7 @@ EclGenericProblem(const EclipseState& eclState,
     , schedule_(schedule)
     , gridView_(gridView)
     , mixControls_(schedule)
+    , lookUpData_(gridView)
 {
 }
 
@@ -168,14 +169,16 @@ readRockParameters_(const std::vector<Scalar>& cellCenterDepths,
         rockTableIdx_.resize(numElem);
         const auto& num = eclState_.fieldProps().get_int(rock_config.rocknum_property());
         for (std::size_t elemIdx = 0; elemIdx < numElem; ++ elemIdx) {
-            rockTableIdx_[elemIdx] = num[elemIdx] - 1;
+            auto coarseElemIdx = this->getLookUpData(elemIdx);
+            rockTableIdx_[elemIdx] = num[coarseElemIdx] - 1;
             auto fmtError =
-                [&num,elemIdx,&ijkIndex,&rock_config](const char* type, std::size_t size)
+                [&num,coarseElemIdx,&ijkIndex,&rock_config](const char* type, std::size_t size)
                 {
                     return fmt::format("{} table index {} for elem {} read from {}"
                                       " is is out of bounds for number of tables {}",
-                                      type, num[elemIdx], ijkIndex(elemIdx),
-                                      rock_config.rocknum_property(), size);
+                                       type, num[coarseElemIdx],
+                                       ijkIndex(coarseElemIdx),
+                                       rock_config.rocknum_property(), size);
                 };
             if (!rockCompPoroMult_.empty() &&
                 rockTableIdx_[elemIdx] >= rockCompPoroMult_.size()) {
@@ -210,7 +213,8 @@ readRockParameters_(const std::vector<Scalar>& cellCenterDepths,
             if (!rockTableIdx_.empty()) {
                 tableIdx = rockTableIdx_[elemIdx];
             }
-            overburdenPressure_[elemIdx] = overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true);
+            overburdenPressure_[elemIdx] =
+                overburdenTables[tableIdx].eval(cellCenterDepths[elemIdx], /*extrapolation=*/true);
         }
     }
 }
@@ -353,7 +357,7 @@ rockFraction(unsigned elementIdx, unsigned timeIdx) const
     // geometric volume of the element. Note that it can
     // be larger than 1.0 if porevolume multipliers are used
     // to for instance implement larger boundary cells
-    Scalar porosity = poro[elementIdx];
+    Scalar porosity = poro[this->getLookUpData(elementIdx)];
     return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
 }
 
@@ -370,10 +374,10 @@ updateNum(const std::string& name, std::vector<T>& numbers, std::size_t num_regi
     unsigned numElems = gridView_.size(/*codim=*/0);
     numbers.resize(numElems);
     for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx) {
-        if (numData[elemIdx] > (int)num_regions) {
+        if (numData[this->getLookUpData(elemIdx)] > (int)num_regions) {
             throw std::runtime_error("Values larger than maximum number of regions " + std::to_string(num_regions) + " provided in " + name);
-        } else if (numData[elemIdx] > 0) {
-            numbers[elemIdx] = static_cast<T>(numData[elemIdx]) - 1;
+        } else if (numData[this->getLookUpData(elemIdx)] > 0) {
+            numbers[elemIdx] = static_cast<T>(numData[this->getLookUpData(elemIdx)]) - 1;
         } else {
             throw std::runtime_error("zero or negative values provided for region array: " + name);
         }

ebos/eclproblem.hh

@@ -2026,6 +2026,11 @@ protected:
 
     void readEclRestartSolution_()
     {
+        // Throw an exception if the grid has LGRs. Refined grid are not supported for restart.
+        if(this->simulator().vanguard().grid().maxLevel() > 0) {
+            throw std::invalid_argument("Refined grids are not yet supported for restart ");
+        }
+
         // Set the start time of the simulation
         auto& simulator = this->simulator();
         const auto& schedule = simulator.vanguard().schedule();