Add support for loading unstructured grids from external files

CMakeLists.txt

@@ -685,6 +685,7 @@ macro(opm-simulators_targets_hook)
 
   target_sources(test_outputdir PRIVATE $<TARGET_OBJECTS:moduleVersion>)
   target_sources(test_equil PRIVATE $<TARGET_OBJECTS:moduleVersion>)
+  target_sources(test_grid_from_file PRIVATE $<TARGET_OBJECTS:moduleVersion>)
   target_sources(test_group_higher_constraints PRIVATE $<TARGET_OBJECTS:moduleVersion>)
   target_sources(test_injection_topup_phase_validation PRIVATE $<TARGET_OBJECTS:moduleVersion>)
   target_sources(test_RestartSerialization PRIVATE $<TARGET_OBJECTS:moduleVersion>)

CMakeLists_files.cmake

@@ -507,6 +507,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_OilSatfuncConsistencyChecks.cpp
   tests/test_outputdir.cpp
   tests/test_parametersystem.cpp
+  tests/test_grid_from_file.cpp
   tests/test_parallel_wbp_sourcevalues.cpp
   tests/test_parallelwellinfo.cpp
   tests/test_partitionCells.cpp

opm/simulators/flow/CpGridVanguard.hpp

@@ -308,7 +308,25 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
     std::function<std::array<double,dimensionworld>(int)>
     cellCentroids() const
     {
-        return this->cellCentroids_(this->cartesianIndexMapper(), true);
+        if (!cellCentroidsFromGrid_) {
+            return this->cellCentroids_(this->cartesianIndexMapper(), true);
+        } else {
+            // Use centroids from the grid file if available
+            auto centroidIter = this->grid().beginCellCentroids();
+            int num_cells = this->gridView().size(0);
+            return [centroidIter, num_cells](int i) -> std::array<double, dimensionworld> {
+                if (i < 0 || i >= num_cells) return std::array<double, dimensionworld>{};
+                auto it = centroidIter;
+                std::advance(it, i);
+                const auto& centroid = *it;
+                return {centroid[0], centroid[1], centroid[2]};
+            };
+        }
+    }
+
+    bool gridFromFile() const
+    {
+        return cellCentroidsFromGrid_;
     }
 
     const std::vector<int>& globalCell()
@@ -319,6 +337,11 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
 protected:
     void createGrids_()
     {
+        // Check if grid will be loaded from file
+        const auto gridFileName = Parameters::Get<Parameters::UnstructuredGridFileName>();
+        if (!gridFileName.empty()) {
+            cellCentroidsFromGrid_ = true;
+        }
         this->doCreateGrids_(this->edgeConformal(), this->eclState());
     }
 
@@ -333,7 +356,8 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
                                                     getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal ||
                                                     getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
                                                     getPropValue<TypeTag, Properties::EnableDiffusion>(),
-                                                    getPropValue<TypeTag, Properties::EnableDispersion>()));
+                                                    getPropValue<TypeTag, Properties::EnableDispersion>(),
+                                                    cellCentroidsFromGrid_));
         globalTrans_->update(false, TransmissibilityType::TransUpdateQuantities::Trans);
     }
 
@@ -364,6 +388,9 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
     // diffusivity_ abd dispersivity_. The main reason is to reduce the memory usage for rank 0
     // during parallel running.
     std::unique_ptr<TransmissibilityType> globalTrans_;
+
+    // Flag to indicate if cell centroids should be read from the grid file
+    bool cellCentroidsFromGrid_ = false;
 };
 
 } // namespace Opm

opm/simulators/flow/EclGenericWriter_impl.hpp

@@ -49,7 +49,10 @@
 #include <opm/output/eclipse/RestartValue.hpp>
 #include <opm/output/eclipse/Summary.hpp>
 
+#include <opm/models/utils/parametersystem.hpp>
+
 #include <opm/simulators/flow/EclGenericWriter.hpp>
+#include <opm/simulators/flow/FlowGenericVanguard.hpp>
 
 #if HAVE_MPI
 #include <opm/simulators/utils/MPISerializer.hpp>
@@ -246,10 +249,22 @@ EclGenericWriter(const Schedule& schedule,
     outputNnc_.resize(1);
 
     if (this->collectOnIORank_.isIORank()) {
-        this->eclIO_ = std::make_unique<EclipseIO>
-            (this->eclState_,
-             UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
-             this->schedule_, summaryConfig, "", enableEsmry);
+        // If an external unstructured grid is provided, use input grid as-is.
+        // Otherwise, build an EclipseGrid that reflects the simulation grid.
+        const auto gridFileName = Parameters::Get<Parameters::UnstructuredGridFileName>();
+        if (!gridFileName.empty()) {
+            this->eclIO_ = std::make_unique<EclipseIO>
+                (this->eclState_,
+                 eclState_.getInputGrid(),
+                 this->schedule_, summaryConfig, "", enableEsmry);
+        }
+        else {
+            assert(equilGrid != nullptr);
+            this->eclIO_ = std::make_unique<EclipseIO>
+                (this->eclState_,
+                 UgGridHelpers::createEclipseGrid(*equilGrid, eclState_.getInputGrid()),
+                 this->schedule_, summaryConfig, "", enableEsmry);
+        }
     }
 
     // create output thread if enabled and rank is I/O rank

opm/simulators/flow/FlowGenericVanguard.cpp

@@ -544,8 +544,9 @@ void FlowGenericVanguard::registerParameters_()
     // register here for the use in the tests without BlackoilModelParameters
     Parameters::Register<Parameters::UseMultisegmentWell>
         ("Use the well model for multi-segment wells instead of the one for single-segment wells");
+    Parameters::Register<Parameters::UnstructuredGridFileName>
+        ("Filename for unstructured grid input. If empty, the grid will be constructed from the ECL deck.");
 }
-
 template void FlowGenericVanguard::registerParameters_<double>();
 
 #if FLOW_INSTANTIATE_FLOAT

opm/simulators/flow/FlowGenericVanguard.hpp

@@ -88,7 +88,7 @@ struct ZoltanImbalanceTol { static constexpr Scalar value = 1.1; };
 struct ZoltanPhgEdgeSizeThreshold { static constexpr auto value = 0.35; };
 
 struct ZoltanParams { static constexpr auto value = "graph"; };
-
+struct UnstructuredGridFileName { static constexpr auto* value = ""; };
 } // namespace Opm::Parameters
 
 namespace Opm {
@@ -145,6 +145,11 @@ class FlowGenericVanguard {
      */
     static std::string canonicalDeckPath(const std::string& caseName);
 
+    bool gridFromFile() const
+    {
+        return false;
+    }
+
     /*!
      * \brief Returns the wall time required to set up the simulator before it was born.
      */

opm/simulators/flow/FlowProblem.hpp

@@ -231,7 +231,8 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
                               (energyModuleType == EnergyModules::FullyImplicitThermal ||
                                energyModuleType == EnergyModules::SequentialImplicitThermal),
                               enableDiffusion,
-                              enableDispersion)
+                              enableDispersion,
+                              simulator.vanguard().gridFromFile())
         , wellModel_(simulator, this->iterationContext())
         , aquiferModel_(simulator)
         , pffDofData_(simulator.gridView(), this->elementMapper())

opm/simulators/flow/GenericCpGridVanguard.cpp

@@ -40,6 +40,8 @@
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/LgrCollection.hpp>
 
+#include <opm/models/utils/parametersystem.hpp>
+
 #include <opm/simulators/utils/ParallelEclipseState.hpp>
 #include <opm/simulators/utils/ParallelSerialization.hpp>
 #include <opm/simulators/utils/PropsDataHandle.hpp>
@@ -465,46 +467,107 @@ doCreateGrids_(const bool edge_conformal, EclipseState& eclState)
         OpmLog::info("\nProcessing grid");
     }
 
+    const auto gridFileName = Parameters::Get<Parameters::UnstructuredGridFileName>();
+
     // --- Create grid ---
     OPM_TIMEBLOCK(createGrids);
 
+    std::vector<std::size_t> removed_cells;
+    if (gridFileName.empty()) {
 #if HAVE_MPI
-    this->grid_ = std::make_unique<Dune::CpGrid>(FlowGenericVanguard::comm());
+        this->grid_ = std::make_unique<Dune::CpGrid>(FlowGenericVanguard::comm());
 #else
-    this->grid_ = std::make_unique<Dune::CpGrid>();
+        this->grid_ = std::make_unique<Dune::CpGrid>();
 #endif
 
-    // Note: removed_cells is guaranteed to be empty on ranks other than 0.
-    auto removed_cells = this->grid_
-        ->processEclipseFormat(input_grid,
-                               &eclState,
-                               /* isPeriodic = */ false,
-                               /* flipNormals = */ false,
-                               /* clipZ = */ false,
-                               edge_conformal);
+        // Note: removed_cells is guaranteed to be empty on ranks other than 0.
+        removed_cells = this->grid_
+            ->processEclipseFormat(input_grid,
+                                   &eclState,
+                                   /* isPeriodic = */ false,
+                                   /* flipNormals = */ false,
+                                   /* clipZ = */ false,
+                                   edge_conformal);
+
+        if (isRoot) {
+            const auto& active_porv = eclState.fieldProps().porv(false);
+            const auto& unit_system = eclState.getUnits();
+            const auto& volume_unit = unit_system.name(UnitSystem::measure::volume);
+            double total_pore_volume = unit_system.from_si(UnitSystem::measure::volume,
+                                                           std::accumulate(active_porv.begin(), active_porv.end(), 0.0));
+            OpmLog::info(fmt::format("Total number of active cells: {} / total pore volume: {:0.0f} {}",
+                                     grid_->numCells(), total_pore_volume , volume_unit));
+
+            double removed_pore_volume =
+                std::accumulate(removed_cells.begin(), removed_cells.end(), 0.0,
+                                [&active_porv, &eclState](const auto acc, const auto global_index)
+                                { return acc + active_porv[eclState.getInputGrid().activeIndex(global_index)]; });
+
+            if (removed_pore_volume > 0) {
+                removed_pore_volume = unit_system.from_si(UnitSystem::measure::volume, removed_pore_volume);
+                OpmLog::info(fmt::format("Removed {} cells with a pore volume of {:0.0f} "
+                                         "{} ({:5.3f} %) due to MINPV/MINPVV",
+                                         removed_cells.size(),
+                                         removed_pore_volume,
+                                         volume_unit,
+                                         100 * removed_pore_volume / total_pore_volume));
+            }
+        }
+    }
+    else {
+#if HAVE_MPI
+        if (FlowGenericVanguard::comm().size() > 1) {
+            OpmLog::info(fmt::format("Using unstructured grid from file (MPI): {}", gridFileName));
+            this->grid_ = std::make_unique<Dune::CpGrid>(FlowGenericVanguard::comm());
+            this->grid_->readUnstructuredGridFile(gridFileName);
+
+            if (isRoot) {
+                if (!input_grid->allActive()) {
+                    OPM_THROW(std::runtime_error,
+                              fmt::format("Cannot use unstructured grid file '{}': the Eclipse input grid "
+                                          "contains inactive cells ({} active out of {} total). "
+                                          "Unstructured grid files require all cells to be active.",
+                                          gridFileName,
+                                          input_grid->getNumActive(),
+                                          input_grid->getCartesianSize()));
+                }
 
-    if (isRoot) {
-        const auto& active_porv = eclState.fieldProps().porv(false);
-        const auto& unit_system = eclState.getUnits();
-        const auto& volume_unit = unit_system.name(UnitSystem::measure::volume);
-        double total_pore_volume = unit_system.from_si(UnitSystem::measure::volume,
-                                                       std::accumulate(active_porv.begin(), active_porv.end(), 0.0));
-        OpmLog::info(fmt::format("Total number of active cells: {} / total pore volume: {:0.0f} {}",
-                                 grid_->numCells(), total_pore_volume , volume_unit));
-
-        double removed_pore_volume =
-            std::accumulate(removed_cells.begin(), removed_cells.end(), 0.0,
-                            [&active_porv, &eclState](const auto acc, const auto global_index)
-                            { return acc + active_porv[eclState.getInputGrid().activeIndex(global_index)]; });
-
-        if (removed_pore_volume > 0) {
-            removed_pore_volume = unit_system.from_si(UnitSystem::measure::volume, removed_pore_volume);
-            OpmLog::info(fmt::format("Removed {} cells with a pore volume of {:0.0f} "
-                                     "{} ({:5.3f} %) due to MINPV/MINPVV",
-                                     removed_cells.size(),
-                                     removed_pore_volume,
-                                     volume_unit,
-                                     100 * removed_pore_volume / total_pore_volume));
+                const int numCellsFile = this->grid_->numCells();
+                const int numCellsEcl  = static_cast<int>(input_grid->getNumActive());
+                if (numCellsFile != numCellsEcl) {
+                    OPM_THROW(std::runtime_error,
+                              fmt::format("Cell count mismatch: unstructured grid file '{}' has {} cells, "
+                                          "but the Eclipse input grid has {} active cells.",
+                                          gridFileName, numCellsFile, numCellsEcl));
+                }
+            }
+        }
+        else
+#endif
+        {
+            OpmLog::info(fmt::format("Using unstructured grid from file: {}", gridFileName));
+            this->grid_ = std::make_unique<Dune::CpGrid>(gridFileName);
+
+            if (isRoot) {
+                if (!input_grid->allActive()) {
+                    OPM_THROW(std::runtime_error,
+                              fmt::format("Cannot use unstructured grid file '{}': the Eclipse input grid "
+                                          "contains inactive cells ({} active out of {} total). "
+                                          "Unstructured grid files require all cells to be active.",
+                                          gridFileName,
+                                          input_grid->getNumActive(),
+                                          input_grid->getCartesianSize()));
+                }
+
+                const int numCellsFile = this->grid_->numCells();
+                const int numCellsEcl  = static_cast<int>(input_grid->getNumActive());
+                if (numCellsFile != numCellsEcl) {
+                    OPM_THROW(std::runtime_error,
+                              fmt::format("Cell count mismatch: unstructured grid file '{}' has {} cells, "
+                                          "but the Eclipse input grid has {} active cells.",
+                                          gridFileName, numCellsFile, numCellsEcl));
+                }
+            }
         }
     }