[DualPorosity] (4/5) Apply the dual-continuum policy and twin co-location

CMakeLists_files.cmake

@@ -491,6 +491,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_milu.cpp
   tests/test_multmatrixtransposed.cpp
   tests/test_networkpressure.cpp
+  tests/test_dualporosity_trans.cpp
   tests/test_nonnc.cpp
   tests/test_norne_pvt.cpp
   tests/test_OilSatfuncConsistencyChecks.cpp

opm/simulators/flow/CpGridVanguard.hpp

@@ -27,6 +27,8 @@
 #ifndef OPM_CPGRID_VANGUARD_HPP
 #define OPM_CPGRID_VANGUARD_HPP
 
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/TimingMacros.hpp>
 
 #include <opm/models/common/multiphasebaseproperties.hh>
@@ -40,6 +42,7 @@
 #include <functional>
 #include <memory>
 #include <stdexcept>
+#include <string>
 #include <tuple>
 #include <vector>
 
@@ -236,6 +239,18 @@ class CpGridVanguard : public FlowBaseVanguard<TypeTag>
     void loadBalance()
     {
 #if HAVE_MPI
+        // Dual-continuum runs are supported in serial only.  The check is collective -- the
+        // RUNSPEC flags are distributed to every rank -- and it is raised before load balancing
+        // so the run stops with a diagnosable message instead of failing later inside the
+        // partitioner.
+        if (this->eclState().runspec().porosityModel().dualContinuum() && this->comm().size() > 1) {
+            const std::string msg =
+                "Dual-continuum runs (DUALPORO/DUALPERM) are supported in serial only. "
+                "Rerun on a single process.";
+            OpmLog::error(msg);
+            OPM_THROW(std::runtime_error, msg);
+        }
+
         if (const auto& extPFile = this->externalPartitionFile();
             !extPFile.empty() && (extPFile != "none"))
         {

opm/simulators/flow/EclGenericWriter.hpp

@@ -231,6 +231,8 @@ class EclGenericWriter
     /// Returns true if either of the two connected cells belongs to a numerical aquifer.
     bool isNumAquConn_(const std::size_t cartIdx1, const std::size_t cartIdx2) const;
 
+    bool isDualPorosityTwin_(const std::size_t cartIdx1, const std::size_t cartIdx2) const;
+
     /// Create LevelCartesianIndexMapper.
     ///
     /// For CpGrid, the LevelCartesianIndexMapper constructor takes the grid

opm/simulators/flow/EclGenericWriter_impl.hpp

@@ -562,6 +562,13 @@ computeTrans_(const std::vector<std::unordered_map<int,int>>&  levelCartToLevelC
                 continue;
             }
 
+            if (isDualPorosityTwin_(originCartIdxIn, originCartIdxOut)) {
+                // The matrix-fracture coupling is always an NNC for the
+                // purpose of file output — the twin cells are neighbours
+                // only in the doubled Cartesian bookkeeping.
+                continue;
+            }
+
             const auto minLevelCartIdx = std::min(levelCartIdxIn, levelCartIdxOut);
             const auto maxLevelCartIdx = std::max(levelCartIdxIn, levelCartIdxOut);
 
@@ -608,6 +615,22 @@ isCartesianNeighbour_(const std::array<int,3>& levelCartDims,
         || (diff == (levelCartDims[0] * levelCartDims[1]));
 }
 
+template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
+bool
+EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
+isDualPorosityTwin_(const std::size_t cartIdx1, const std::size_t cartIdx2) const
+{
+    if (! this->eclState_.runspec().porosityModel().dualContinuum()) {
+        return false;
+    }
+
+    const auto& inputGrid = this->eclState_.getInputGrid();
+    const auto lo = std::min(cartIdx1, cartIdx2);
+    const auto hi = std::max(cartIdx1, cartIdx2);
+    return inputGrid.isFractureCell(hi) && !inputGrid.isFractureCell(lo)
+        && inputGrid.matrixTwin(hi) == lo;
+}
+
 template<class Grid, class EquilGrid, class GridView, class ElementMapper, class Scalar>
 bool
 EclGenericWriter<Grid,EquilGrid,GridView,ElementMapper,Scalar>::
@@ -854,7 +877,8 @@ exportNncStructure_(const std::vector<std::unordered_map<int,int>>& levelCartToL
         assert (entry.cell2 >= entry.cell1);
 
         if (! isCartesianNeighbour_(level0CartDims, entry.cell1, entry.cell2) ||
-            isNumAquConn_(entry.cell1, entry.cell2))
+            isNumAquConn_(entry.cell1, entry.cell2) ||
+            isDualPorosityTwin_(entry.cell1, entry.cell2))
         {
             bool foundNncEdit = false;
             auto trans = entry.trans;

opm/simulators/flow/FlowBaseVanguard.hpp

@@ -363,19 +363,32 @@ class FlowBaseVanguard : public BaseVanguard<TypeTag>,
         ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
 
         const auto num_aqu_cells = this->allAquiferCells();
+        const bool dualPorosity = this->eclState().runspec().porosityModel().dualContinuum();
+        // Bind the input grid only when it is needed: it is available on the I/O rank only,
+        // so an unconditional binding would make every parallel run -- dual continuum or not --
+        // reach ParallelEclipseState::getInputGrid() and throw.  Dual-continuum runs are serial
+        // (see CpGridVanguard::loadBalance), so the guarded binding is always on rank 0.
+        const auto* inputGrid = dualPorosity ? &this->eclState().getInputGrid() : nullptr;
 
         for(const auto& element : elements(this->gridView())) {
             const unsigned int elemIdx = elemMapper.index(element);
             cellCenterDepth_[elemIdx] = cellCenterDepth(element);
+            const unsigned int global_index = cartesianIndex(elemIdx);
 
             if (!num_aqu_cells.empty()) {
-                const unsigned int global_index = cartesianIndex(elemIdx);
                 const auto search = num_aqu_cells.find(global_index);
                 if (search != num_aqu_cells.end()) {
                     // updating the cell depth using aquifer cell depth
                     cellCenterDepth_[elemIdx] = search->second->depth;
                 }
             }
+
+            // Dual porosity: fracture cells are co-located with their matrix
+            // twins — the input grid carries the twin's depth (the geometric
+            // stacking of the fracture half is bookkeeping only).
+            if (dualPorosity && inputGrid->isFractureCell(global_index)) {
+                cellCenterDepth_[elemIdx] = inputGrid->getCellDepth(global_index);
+            }
         }
     }
     void updateCellThickness_()

opm/simulators/flow/Transmissibility.hpp

@@ -271,6 +271,8 @@ class Transmissibility {
 
     void extractPermeability_(const std::function<unsigned int(unsigned int)>& map);
 
+    void applyDualPorosityPermScaling_(const std::function<unsigned int(unsigned int)>& map);
+
     void extractPorosity_();
 
     void extractDispersion_();

opm/simulators/flow/Transmissibility_impl.hpp

@@ -38,6 +38,7 @@
 #include <opm/grid/utility/ElementChunks.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/FieldPropsManager.hpp>
 #include <opm/input/eclipse/EclipseState/Grid/TransMult.hpp>
@@ -198,6 +199,17 @@ update(bool global, const TransUpdateQuantities update_quantities,
     // whether only update the permeability related transmissibility
     const bool onlyTrans = (update_quantities == TransUpdateQuantities::Trans);
     const auto& cartDims = cartMapper_.cartesianDimensions();
+    const auto& porosityModel = eclState_.runspec().porosityModel();
+    const bool dualPorosity = porosityModel.dualContinuum();
+    const bool dualPermeability = porosityModel.dualPermeability();
+    // Twin classification is arithmetic on the global Cartesian index: the fracture half is
+    // the upper half of the index range.  Deriving it from the Cartesian dimensions keeps this
+    // rank-local -- EclipseState::getInputGrid() is available on the I/O rank only.
+    const std::size_t matrixCellCount = EclipseGrid::matrixCellCount(cartDims);
+    const auto isFractureCell = [matrixCellCount](const std::size_t cartIdx)
+    {
+        return cartIdx >= matrixCellCount;
+    };
     const auto& transMult = eclState_.getTransMult();
     const auto& comm = gridView_.comm();
     ElementMapper elemMapper(gridView_, Dune::mcmgElementLayout());
@@ -546,6 +558,22 @@ update(bool global, const TransUpdateQuantities update_quantities,
                                                            faceIdToDir(inside.faceIdx));
                 }
 
+                // Dual-continuum runs: the matrix and fracture halves never
+                // connect through grid faces (their coupling comes exclusively
+                // through the input NNCs).  Matrix-matrix faces carry flow
+                // only in dual-permeability runs; in single-permeability dual
+                // porosity the matrix half has no internal flow.
+                if (dualPorosity) {
+                    const bool insideFracture  = isFractureCell(inside.cartElemIdx);
+                    const bool outsideFracture = isFractureCell(outside.cartElemIdx);
+                    if (insideFracture != outsideFracture) {
+                        trans = 0.0;
+                    }
+                    else if (!insideFracture && !dualPermeability) {
+                        trans = 0.0;
+                    }
+                }
+
                 transMap.insert_or_assign(details::isId(inside.elemIdx, outside.elemIdx), trans);
 
                 // update the "thermal half transmissibility" for the intersection
@@ -678,12 +706,48 @@ extractPermeability_()
 
         // for now we don't care about non-diagonal entries
 
+        this->applyDualPorosityPermScaling_([](const unsigned int i) { return i; });
     }
     else
         throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
                                "(The PERM{X,Y,Z} keywords are missing)");
 }
 
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
+applyDualPorosityPermScaling_(const std::function<unsigned int(unsigned int)>& map)
+{
+    // Dual porosity: the effective fracture permeability is scaled by the
+    // fracture porosity unless the run disables that scaling. Matrix cells
+    // are untouched, and so is the matrix-fracture coupling transmissibility
+    // (it is computed from the matrix permeability upstream and arrives here
+    // as an input NNC). The rule itself is Runspec's -- the well connection
+    // factors apply the same one, and spelling it out separately here is how
+    // the two last diverged.
+    if (!eclState_.runspec().porosityModel().fracturePermeabilityScalingActive())
+        return;
+
+    const auto& fp = eclState_.fieldProps();
+    const std::vector<double>& poroData = this->lookUpData_.assignFieldPropsDoubleOnLeaf(fp, "PORO");
+
+    // Classify twins from the Cartesian dimensions rather than from the input grid:
+    // the same arithmetic EclipseGrid uses, but available on every process. This file
+    // asked the question two different ways -- the face policy already derives it
+    // locally -- and the input-grid form is the pattern that broke every parallel run.
+    const std::size_t matrixCellCount =
+        EclipseGrid::matrixCellCount(cartMapper_.cartesianDimensions());
+
+    // The porosity must be read through the same element-to-input mapping
+    // the permeability extraction used, so reordered grids scale the right
+    // cells.
+    for (std::size_t elemIdx = 0; elemIdx < permeability_.size(); ++elemIdx) {
+        const auto inputDofIdx = map(static_cast<unsigned int>(elemIdx));
+        if (static_cast<std::size_t>(cartMapper_.cartesianIndex(elemIdx)) >= matrixCellCount) {
+            permeability_[elemIdx] *= poroData[inputDofIdx];
+        }
+    }
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 extractPermeability_(const std::function<unsigned int(unsigned int)>& map)
@@ -725,6 +789,8 @@ extractPermeability_(const std::function<unsigned int(unsigned int)>& map)
         }
 
         // for now we don't care about non-diagonal entries
+
+        this->applyDualPorosityPermScaling_(map);
     }
     else {
         throw std::logic_error("Can't read the intrinsic permeability from the ecl state. "
@@ -781,6 +847,16 @@ void Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 removeNonCartesianTransmissibilities_(bool removeAll)
 {
     const auto& cartDims = cartMapper_.cartesianDimensions();
+
+    // A dual-continuum coupling is a connection between a cell and its twin, exactly
+    // half the Cartesian index range apart. It is the physics of the run, not a sparse
+    // non-neighbour connection the deck happened to add, so it must survive both the
+    // threshold prune and a blanket removal: a tight matrix with a small shape factor
+    // produces a legitimately small transmissibility, and zeroing it would strand the
+    // matrix continuum silently while the run completed and the material balance closed.
+    const bool dualPorosity = eclState_.runspec().porosityModel().dualContinuum();
+    const std::size_t twinGap = EclipseGrid::matrixCellCount(cartDims);
+
     for (auto&& trans: trans_) {
         //either remove all NNC transmissibilities or those less than the threshold (by default 1e-6 in the deck's unit system)
         if (removeAll || trans.second < transmissibilityThreshold_) {
@@ -796,6 +872,11 @@ removeNonCartesianTransmissibilities_(bool removeAll)
                 continue;
             }
 
+            // the matrix-fracture coupling, kept for the reason above
+            if (dualPorosity && (static_cast<std::size_t>(gc2 - gc1) == twinGap)) {
+                continue;
+            }
+
             trans.second = 0.0;
         }
     }
@@ -1191,6 +1272,24 @@ applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompr
         }
 
         if (low == -1 || high == -1) {
+            // A dual-continuum coupling must never be silently discarded.  In a parallel run a
+            // cell missing from this rank's map is inactive OR owned by another rank, and both
+            // arrive here -- so this is the path by which a partition that separates a twin pair
+            // drops its coupling.  Measured on a 3x3x2 case: none split at two processes, at
+            // least four of nine at four processes.  Dual-continuum runs are refused before load
+            // balancing for exactly this reason; if that guard is ever lifted, this must be an
+            // error rather than a warning.
+            if (eclState_.runspec().porosityModel().dualContinuum()) {
+                const auto& cd = cartMapper_.cartesianDimensions();
+                const std::size_t half =
+                    (static_cast<std::size_t>(cd[0]) * cd[1] * cd[2]) / 2;
+                if ((c2 > c1 ? c2 - c1 : c1 - c2) == half) {
+                    OPM_THROW(std::runtime_error,
+                              "Dual-continuum coupling between cells " + std::to_string(c1) +
+                              " and " + std::to_string(c2) + " cannot be built: one of the two "
+                              "cells is inactive or not owned by this process.");
+                }
+            }
             // Discard the NNC if it is between active cell and inactive cell
             std::ostringstream sstr;
             sstr << "NNC between active and inactive cells ("

opm/simulators/flow/equil/InitStateEquil.hpp

@@ -737,7 +737,8 @@ class InitialStateComputer
     template <class RMap>
     void updateInitialSaltSaturation_(const EclipseState& eclState, const RMap& reg);
 
-    void updateCellProps_(const GridView& gridView,
+    void updateCellProps_(const EclipseState& eclipseState,
+                          const GridView& gridView,
                           const NumericalAquifers& aquifer);
 
     void applyNumericalAquifers_(const GridView& gridView,

opm/simulators/flow/equil/InitStateEquil_impl.hpp

@@ -1514,7 +1514,7 @@ InitialStateComputer(MaterialLawManager& materialLawManager,
     // Querry cell depth, cell top-bottom.
     // numerical aquifer cells might be specified with different depths.
     const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
-    updateCellProps_(gridView, num_aquifers);
+    updateCellProps_(eclipseState, gridView, num_aquifers);
 
     // Get the equilibration records.
     const std::vector<EquilRecord> rec = getEquil(eclipseState);
@@ -1863,7 +1863,8 @@ void InitialStateComputer<FluidSystem,
                           GridView,
                           ElementMapper,
                           CartesianIndexMapper>::
-updateCellProps_(const GridView& gridView,
+updateCellProps_(const EclipseState& eclipseState,
+                 const GridView& gridView,
                  const NumericalAquifers& aquifer)
 {
     ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
@@ -1877,6 +1878,9 @@ updateCellProps_(const GridView& gridView,
     auto elemIt = gridView.template begin</*codim=*/0>();
     const auto& elemEndIt = gridView.template end</*codim=*/0>();
     const auto num_aqu_cells = aquifer.allAquiferCells();
+    const bool dualPorosity = eclipseState.runspec().porosityModel().dualContinuum();
+    // Bound only under dual continuum: getInputGrid() is I/O-rank only (see FlowBaseVanguard).
+    const auto* dpInputGrid = dualPorosity ? &eclipseState.getInputGrid() : nullptr;
     for (; elemIt != elemEndIt; ++elemIt) {
         const Element& element = *elemIt;
         const unsigned int elemIdx = elemMapper.index(element);
@@ -1898,6 +1902,19 @@ updateCellProps_(const GridView& gridView,
                 cellZMinMax_[elemIdx].second += depth_change_num_aqu;
             }
         }
+
+        // Dual porosity: equilibrate the fracture cell at its matrix twin's
+        // depth (the input grid carries it) — the geometric stacking of the
+        // fracture half is bookkeeping only.
+        if (dualPorosity && dpInputGrid->isFractureCell(cartIx)) {
+            const Scalar depth_change_dp =
+                dpInputGrid->getCellDepth(cartIx) - cellCenterDepth_[elemIdx];
+            cellCenterDepth_[elemIdx] += depth_change_dp;
+            cellZSpan_[elemIdx].first += depth_change_dp;
+            cellZSpan_[elemIdx].second += depth_change_dp;
+            cellZMinMax_[elemIdx].first += depth_change_dp;
+            cellZMinMax_[elemIdx].second += depth_change_dp;
+        }
     }
 }
 

tests/test_ParallelSerialization.cpp

@@ -268,6 +268,7 @@ TEST_FOR_TYPE_NAMED(Network::Node, NetworkNode)
 TEST_FOR_TYPE(OilVaporizationProperties)
 TEST_FOR_TYPE(PAvg)
 TEST_FOR_TYPE(Phases)
+TEST_FOR_TYPE(PorosityModel)
 TEST_FOR_TYPE(PlymwinjTable)
 TEST_FOR_TYPE(PlyshlogTable)
 TEST_FOR_TYPE(PvcdoTable)