Merge pull request #554 from totto82/nnc_thpress

NNC and threshold pressure

examples/flow.cpp

@@ -413,7 +413,8 @@ try
     std::map<std::pair<int, int>, double> maxDp;
     computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
     std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
-
+    std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
+    threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
     SimulatorFullyImplicitBlackoil< Grid >  simulator(param,
                                                       grid,
                                                       geoprops,

examples/flow_multisegment.cpp

@@ -413,6 +413,8 @@ try
     std::map<std::pair<int, int>, double> maxDp;
     computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
     std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
+    std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
+    threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
 
     SimulatorFullyImplicitBlackoilMultiSegment< Grid >  simulator(param,
                                                       grid,

examples/flow_polymer.cpp

@@ -287,6 +287,8 @@ try
     std::map<std::pair<int, int>, double> maxDp;
     computeMaxDp(maxDp, deck, eclipseState, *grid->c_grid(), state, *props, gravity[2]);
     std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid(), maxDp);
+    std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geology.nnc(), maxDp);
+    threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
 
     Opm::BlackoilOutputWriter
         outputWriter(cGrid, param, eclipseState, pu,

examples/flow_solvent.cpp

@@ -383,6 +383,8 @@ try
     std::map<std::pair<int, int>, double> maxDp;
     computeMaxDp(maxDp, deck, eclipseState, grid, state, props, gravity[2]);
     std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, grid, maxDp);
+    std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipseState, geoprops.nnc(), maxDp);
+    threshold_pressures.insert(threshold_pressures.begin(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
 
     SimulatorFullyImplicitBlackoilSolvent< Grid >  simulator(param,
                                                       grid,

opm/autodiff/AutoDiffHelpers.hpp

@@ -23,7 +23,11 @@
 
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
+#include <opm/autodiff/GeoProps.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/grid/PinchProcessor.hpp>
+#include <opm/core/props/rock/RockFromDeck.hpp>
+
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -69,7 +73,7 @@ struct HelperOps
 
     /// Constructs all helper vectors and matrices.
     template<class Grid>
-    HelperOps(const Grid& grid, Opm::EclipseStateConstPtr eclState = EclipseStateConstPtr () )
+    HelperOps(const Grid& grid, const NNC& nnc = NNC())
     {
         using namespace AutoDiffGrid;
         const int nc = numCells(grid);
@@ -78,39 +82,31 @@ struct HelperOps
 
         TwoColInt nbi;
         extractInternalFaces(grid, internal_faces, nbi);
-        int num_internal = internal_faces.size();
+        const int num_internal = internal_faces.size();
-        // num_connections may also include non-neighboring connections
-        int num_connections = num_internal;
-        int numNNC = 0;
 
         // handle non-neighboring connections
-        std::shared_ptr<const NNC> nnc = eclState ? eclState->getNNC()
+        const bool has_nnc = nnc.hasNNC();
-            : std::shared_ptr<const Opm::NNC>();
+        size_t numNNC = nnc.numNNC();
-        const bool has_nnc = nnc && nnc->hasNNC();
+
+        // num_connections may also include non-neighboring connections
+        const int num_connections = num_internal + numNNC;
         if (has_nnc) {
-            numNNC = nnc->numNNC();
+            const int *cartDims = AutoDiffGrid::cartDims(grid);
-            num_connections += numNNC;
+            const int numCartesianCells = cartDims[0] * cartDims[1] * cartDims[2];
-            //std::cout << "Added " << numNNC << " NNC" <<std::endl;
-            nbi.resize(num_internal, 2);
 
             // the nnc's acts on global indicies and must be mapped to cell indicies
-            size_t cartesianSize = eclState->getEclipseGrid()->getCartesianSize();
+            std::vector<int> global2localIdx(numCartesianCells,0);
-            std::vector<int> global2localIdx(cartesianSize,0);
             for (int i = 0; i< nc; ++i) {
                 global2localIdx[ globalCell( grid )[i] ] = i;
             }
-            const std::vector<size_t>& NNC1 = nnc->nnc1();
-            const std::vector<size_t>& NNC2 = nnc->nnc2();
             nnc_cells.resize(numNNC,2);
+            nnc_trans.resize(numNNC);
             for (int i = 0; i < numNNC; ++i) {
-                nnc_cells(i,0) = global2localIdx[NNC1[i]];
+                nnc_cells(i,0) = global2localIdx[nnc.nncdata()[i].cell1];
-                nnc_cells(i,1) = global2localIdx[NNC2[i]];
+                nnc_cells(i,1) = global2localIdx[nnc.nncdata()[i].cell2];
+                // store the nnc transmissibilities for later usage.
+                nnc_trans(i) = nnc.nncdata()[i].trans;
             }
-            // store the nnc transmissibilities for later usage.
-            nnc_trans = Eigen::Map<const V>(nnc->trans().data(), numNNC);
-        } else {
-            nnc_trans.resize(0);
-            nnc_cells.resize(0, 2); // Required to have two columns.
         }
 
 
@@ -166,7 +162,6 @@ struct HelperOps
         fulldiv = fullngrad.transpose();
     }
 };
-
 // -------------------- upwinding helper class --------------------
 
 

opm/autodiff/BlackoilModelBase.hpp

@@ -298,7 +298,7 @@ namespace Opm {
         ModelParameters                 param_;
         bool use_threshold_pressure_;
         bool wells_active_;
-        V threshold_pressures_by_interior_face_;
+        V threshold_pressures_by_connection_;
 
         std::vector<ReservoirResidualQuant> rq_;
         std::vector<PhasePresence> phaseCondition_;

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -163,7 +163,7 @@ namespace detail {
         , active_(detail::activePhases(fluid.phaseUsage()))
         , canph_ (detail::active2Canonical(fluid.phaseUsage()))
         , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
-        , ops_   (grid, eclState)
+        , ops_   (grid, geo.nnc())
         , wops_  (wells_)
         , has_disgas_(has_disgas)
         , has_vapoil_(has_vapoil)
@@ -370,19 +370,27 @@ namespace detail {
     template <class Grid, class Implementation>
     void
     BlackoilModelBase<Grid, Implementation>::
-    setThresholdPressures(const std::vector<double>& threshold_pressures_by_face)
+    setThresholdPressures(const std::vector<double>& threshold_pressures)
     {
         const int num_faces = AutoDiffGrid::numFaces(grid_);
-        if (int(threshold_pressures_by_face.size()) != num_faces) {
+        const int num_nnc = geo_.nnc().numNNC();
-            OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures_by_face input, must be equal to number of faces.");
+        const int num_connections = num_faces + num_nnc;
+        if (int(threshold_pressures.size()) != num_connections) {
+            OPM_THROW(std::runtime_error, "Illegal size of threshold_pressures input, must be equal to number of faces + nncs.");
         }
         use_threshold_pressure_ = true;
         // Map to interior faces.
         const int num_ifaces = ops_.internal_faces.size();
-        threshold_pressures_by_interior_face_.resize(num_ifaces);
+        threshold_pressures_by_connection_.resize(num_ifaces + num_nnc);
         for (int ii = 0; ii < num_ifaces; ++ii) {
-            threshold_pressures_by_interior_face_[ii] = threshold_pressures_by_face[ops_.internal_faces[ii]];
+            threshold_pressures_by_connection_[ii] = threshold_pressures[ops_.internal_faces[ii]];
         }
+        // Handle NNCs
+        // Note: the nnc threshold pressures is appended after the face threshold pressures
+        for (int ii = 0; ii < num_nnc; ++ii) {
+            threshold_pressures_by_connection_[ii + num_ifaces] = threshold_pressures[ii + num_faces];
+        }
+
     }
 
 
@@ -2447,14 +2455,14 @@ namespace detail {
         // threshold, that shall have zero flow. Storing the bool
         // Array as a V (a double Array) with 1 and 0 elements, a
         // 1 where flow is allowed, a 0 where it is not.
-        const V high_potential = (dp.value().abs() >= threshold_pressures_by_interior_face_).template cast<double>();
+        const V high_potential = (dp.value().abs() >= threshold_pressures_by_connection_).template cast<double>();
 
         // Create a sparse vector that nullifies the low potential elements.
         const M keep_high_potential(high_potential.matrix().asDiagonal());
 
         // Find the current sign for the threshold modification
         const V sign_dp = sign(dp.value());
-        const V threshold_modification = sign_dp * threshold_pressures_by_interior_face_;
+        const V threshold_modification = sign_dp * threshold_pressures_by_connection_;
 
         // Modify potential and nullify where appropriate.
         dp = keep_high_potential * (dp - threshold_modification);

opm/autodiff/BlackoilSolventModel.hpp

@@ -123,7 +123,7 @@ namespace Opm {
         using Base::has_vapoil_;
         using Base::param_;
         using Base::use_threshold_pressure_;
-        using Base::threshold_pressures_by_interior_face_;
+        using Base::threshold_pressures_by_connection_;
         using Base::rq_;
         using Base::phaseCondition_;
         using Base::well_perforation_pressure_diffs_;

opm/autodiff/GeoProps.hpp

@@ -28,6 +28,8 @@
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
+#include <opm/core/grid/PinchProcessor.hpp>
 #include <opm/common/utility/platform_dependent/disable_warnings.h>
 
 #include <Eigen/Eigen>
@@ -101,10 +103,9 @@ namespace Opm
                 ntg = eclState->getDoubleGridProperty("NTG")->getData();
             }
 
-            // get grid from parser.
+            // Get grid from parser.
-
-            // Get original grid cell volume.
             EclipseGridConstPtr eclgrid = eclState->getEclipseGrid();
+
             // Pore volume.
             // New keywords MINPVF will add some PV due to OPM cpgrid process algorithm.
             // But the default behavior is to get the comparable pore volume with ECLIPSE.
@@ -125,10 +126,14 @@ namespace Opm
             // for MINPV. Note that the change does not effect the pore volume calculations
             // as the pore volume is currently defaulted to be comparable to ECLIPSE, but
             // only the transmissibility calculations.
-            minPvFillProps_(grid, eclState,ntg);
+            bool opmfil = eclgrid->getMinpvMode() == MinpvMode::ModeEnum::OpmFIL;
+            // opmfil is hardcoded to be true. i.e the volume weighting is always used
+            opmfil = true;
+            if (opmfil) {
+                minPvFillProps_(grid, eclState,ntg);
+            }
 
             // Transmissibility
-
             Vector htrans(AutoDiffGrid::numCellFaces(grid));
             Grid* ug = const_cast<Grid*>(& grid);
 
@@ -142,6 +147,38 @@ namespace Opm
             std::vector<double> mult;
             multiplyHalfIntersections_(grid, eclState, ntg, htrans, mult);
 
+            // Handle NNCs
+            if (eclState) {
+                nnc_ = *(eclState->getNNC());
+            }
+
+            // opmfil is hardcoded to be true. i.e the pinch processor is never used
+            if (~opmfil && eclgrid->isPinchActive()) {
+                const double minpv = eclgrid->getMinpvValue();
+                const double thickness = eclgrid->getPinchThresholdThickness();
+                auto transMode = eclgrid->getPinchOption();
+                auto multzMode = eclgrid->getMultzOption();
+                PinchProcessor<Grid> pinch(minpv, thickness, transMode, multzMode);
+
+                std::vector<double> htrans_copy(htrans.size());
+                std::copy_n(htrans.data(), htrans.size(), htrans_copy.begin());
+
+                std::vector<int> actnum;
+                eclgrid->exportACTNUM(actnum);
+
+                auto transMult = eclState->getTransMult();
+                std::vector<double> multz(numCells, 0.0);
+                const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
+
+                for (int i = 0; i < numCells; ++i) {
+                    multz[i] = transMult->getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
+                }
+
+                // Note the pore volume from eclState is used and not the pvol_ calculated above
+                std::vector<double> porv = eclState->getDoubleGridProperty("PORV")->getData();
+                pinch.process(grid, htrans_copy, actnum, multz, porv, nnc_);
+            }           
+
             // combine the half-face transmissibilites into the final face
             // transmissibilites.
             tpfa_trans_compute(ug, htrans.data(), trans_.data());
@@ -191,6 +228,7 @@ namespace Opm
         const double* gravity()          const { return gravity_;}
         Vector&       poreVolume()             { return pvol_   ;}
         Vector&       transmissibility()       { return trans_  ;}
+        const NNC& nnc() const { return nnc_;}
 
     private:
         template <class Grid>
@@ -218,6 +256,13 @@ namespace Opm
         double gravity_[3]; // Size 3 even if grid is 2-dim.
         bool use_local_perm_;
 
+
+        /// Non-neighboring connections
+        NNC nnc_;
+
+
+
+
     };
 
     template <class GridType>
@@ -369,7 +414,6 @@ namespace Opm
             for(auto cellFaceIter = cellFacesRange.begin(), cellFaceEnd = cellFacesRange.end();
                 cellFaceIter != cellFaceEnd; ++cellFaceIter, ++cellFaceIdx)
             {
-
                 // The index of the face in the compressed grid
                 const int faceIdx = *cellFaceIter;
 
@@ -402,9 +446,7 @@ namespace Opm
 
                 int cartesianCellIdx = AutoDiffGrid::globalCell(grid)[cellIdx];
                 auto cellCenter = eclGrid->getCellCenter(cartesianCellIdx);
-
                 for (int indx = 0; indx < dim; ++indx) {
-
                     const double Ci = Opm::UgGridHelpers::faceCentroid(grid, faceIdx)[indx] - cellCenter[indx];
                     dist += Ci*Ci;
                     cn += sgn * Ci * scaledFaceNormal[ indx ];

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -169,7 +169,7 @@ namespace Opm {
         using Base::has_vapoil_;
         using Base::param_;
         using Base::use_threshold_pressure_;
-        using Base::threshold_pressures_by_interior_face_;
+        using Base::threshold_pressures_by_connection_;
         using Base::rq_;
         using Base::phaseCondition_;
         using Base::well_perforation_pressure_diffs_;