Merge pull request #892 from qilicun/pinchprocessor

PinchProcessor

CMakeLists_files.cmake

@@ -196,6 +196,7 @@ list (APPEND TEST_SOURCE_FILES
 	tests/test_wellcollection.cpp
 	tests/test_timer.cpp
 	tests/test_minpvprocessor.cpp
+	tests/test_pinchprocessor.cpp
 	tests/test_gridutilities.cpp
 	tests/test_anisotropiceikonal.cpp
 	tests/test_stoppedwells.cpp
@@ -224,6 +225,7 @@ list (APPEND TEST_DATA_FILES
 	tests/testBlackoilState1.DATA
 	tests/testBlackoilState2.DATA
   tests/testBlackoilState3.DATA
+  tests/testPinch1.DATA
 	tests/wells_manager_data.data
 	tests/wells_manager_data_expanded.data
 	tests/wells_manager_data_wellSTOP.data
@@ -281,6 +283,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/grid/GridManager.hpp
 	opm/core/grid/GridUtilities.hpp
 	opm/core/grid/MinpvProcessor.hpp
+	opm/core/grid/PinchProcessor.hpp
 	opm/core/grid/cart_grid.h
 	opm/core/grid/cornerpoint_grid.h
 	opm/core/grid/cpgpreprocess/facetopology.h

opm/core/grid/PinchProcessor.hpp

@@ -0,0 +1,509 @@
+/*
+  Copyright 2015 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_PINCHPROCESSOR_HEADER_INCLUDED
+#define OPM_PINCHPROCESSOR_HEADER_INCLUDED
+
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/grid/GridHelpers.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/FaceDir.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/PinchMode.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <array>
+#include <iostream>
+#include <algorithm>
+#include <unordered_map>
+#include <limits>
+
+namespace Opm
+{
+
+    template <class Grid>
+    class PinchProcessor
+    {
+    public:
+        /// \brief Create a Pinch processor.
+        /// \param[in]    minpvValue   value in MINPV keyword
+        /// \param[in]    thickness    item 2 in PINCH keyword
+        /// \param[in]    transMode    item 4 in PINCH keyword
+        /// \param[in]    multzMode    item 5 in PINCH keyword
+        PinchProcessor(const double minpvValue,
+                       const double thickness,
+                       const PinchMode::ModeEnum transMode,
+                       const PinchMode::ModeEnum multzMode);
+        /// Generate NNCs for cells which pv is less than MINPV.
+        /// \param[in]    Grid    cpgrid or unstructured grid
+        /// \param[in]    htrans  half cell transmissibility, size is number of cellfaces.
+        /// \param[in]    multz   Z+ transmissibility multiplier for all active cells
+        /// \param[in]    pv      pore volume for all the cartesian cells
+        /// \param[in]    dz      dz for all the cartesian cells
+        /// \param[in]    nnc     non-neighbor connection class
+        /// Algorithm:
+        /// 1. Mark all the cells which pv and dz less than minpvValue and thickness.
+        /// 2. Find out proper pinchouts column and associate top and bottom cells.
+        /// 3. Compute transmissibility for nncs.
+        /// 4. Apply multz due to different multz options.
+        void process(const Grid& grid,
+                     const std::vector<double>& htrans,
+                     const std::vector<int>& actnum,
+                     const std::vector<double>& multz,
+                     const std::vector<double>& pv,
+                     const std::vector<double>& dz,
+                     NNC& nnc);
+
+    private:
+        double minpvValue_;
+        double thickness_;
+        PinchMode::ModeEnum transMode_;
+        PinchMode::ModeEnum multzMode_;
+        
+        /// Mark minpved cells.
+        std::vector<int> getMinpvCells_(const Grid& grid,
+                                        const std::vector<int>& actnum,
+                                        const std::vector<double>& pv,
+                                        const std::vector<double>& dz);
+
+        /// Get the interface for two cells.
+        int interface_(const Grid& grid,
+                       const int cellIdx1,
+                       const int cellIdx2);
+
+        /// Get the proper face for one cell.
+        int interface_(const Grid& grid,
+                       const int cellIdx,
+                       const Opm::FaceDir::DirEnum& faceDir);
+
+        /// Get pinchouts column.
+        std::vector<std::vector<int> >
+        getPinchoutsColumn_(const Grid& grid,
+                            const std::vector<int>& actnum,
+                            const std::vector<double>& pv,
+                            const std::vector<double>& dz);
+        
+        /// Get global cell index.
+        int getGlobalIndex_(const int i, const int j, const int k, const int* dims);
+
+        /// Get cartesian index.
+        std::array<int, 3> getCartIndex_(const int idx,
+                                         const int* dims);
+
+        /// Compute transmissibility for nnc.
+        std::vector<double> transCompute_(const Grid& grid,
+                                          const std::vector<double>& htrans,
+                                          const std::vector<int>& pinCells,
+                                          const std::vector<int>& pinFaces,
+                                          const std::vector<double>& multz);
+
+        /// Get map between half-trans index and the pair of face index and cell index.
+        std::vector<int> getHfIdxMap_(const Grid& grid);
+        
+        /// Get active cell index.
+        int getActiveCellIdx_(const Grid& grid,
+                              const int globalIdx);
+
+        /// Item 4 in PINCH keyword. 
+        void transTopbot_(const Grid& grid,
+                          const std::vector<double>& htrans,
+                          const std::vector<int>& actnum,
+                          const std::vector<double>& multz,
+                          const std::vector<double>& pv,
+                          const std::vector<double>& dz,
+                          NNC& nnc);
+        
+        /// Item 5 in PINCH keyword.
+        std::unordered_multimap<int, double> multzOptions_(const Grid& grid,
+                                                           const std::vector<int>& pinCells,
+                                                           const std::vector<int>& pinFaces,
+                                                           const std::vector<double>& multz,
+                                                           const std::vector<std::vector<int> >& seg);
+
+        /// Apply multz vector to face transmissibility.
+        void applyMultz_(std::vector<double>& trans,
+                         const std::unordered_multimap<int, double>& multzmap);
+
+    };
+
+
+
+    template <class Grid>
+    inline PinchProcessor<Grid>::PinchProcessor(const double minpv,
+                                                const double thickness,
+                                                const PinchMode::ModeEnum transMode,
+                                                const PinchMode::ModeEnum multzMode)
+    {
+        minpvValue_ = minpv;
+        thickness_  = thickness;
+        transMode_  = transMode;
+        multzMode_  = multzMode;
+    }
+
+    
+
+    template <class Grid>
+    inline int PinchProcessor<Grid>::getGlobalIndex_(const int i, const int j, const int k, const int* dims)
+    {
+        return i + dims[0] * (j + dims[1] * k);
+    }
+
+
+
+    template <class Grid>
+    inline std::array<int, 3> PinchProcessor<Grid>::getCartIndex_(const int idx,
+                                                                  const int* dims)
+    {
+        std::array<int, 3> ijk;
+        ijk[0] = (idx % dims[0]);
+        ijk[1] = ((idx / dims[0]) % dims[1]);
+        ijk[2] = ((idx / dims[0]) / dims[1]);
+
+        return ijk;
+    }
+
+
+
+    template<class Grid>
+    inline int PinchProcessor<Grid>::interface_(const Grid& grid,
+                                                const int cellIdx1,
+                                                const int cellIdx2)
+    {
+        const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
+        int commonFace = -1;
+        const int actCellIdx1 = getActiveCellIdx_(grid, cellIdx1);
+        const int actCellIdx2 = getActiveCellIdx_(grid, cellIdx2);
+        const auto cellFacesRange1 = cell_faces[actCellIdx1];
+        const auto cellFacesRange2 = cell_faces[actCellIdx2];
+        for (const auto& f1 : cellFacesRange1) {
+            for (const auto& f2 : cellFacesRange2) {
+                if (f1 == f2) {
+                    commonFace = f1;
+                    break;
+                }
+            }
+        }
+
+        if (commonFace == -1) {
+            const auto dims = Opm::UgGridHelpers::cartDims(grid);
+            const auto ijk1 = getCartIndex_(cellIdx1, dims);
+            const auto ijk2 = getCartIndex_(cellIdx2, dims);
+
+            OPM_THROW(std::logic_error, "Couldn't find the common face for cell " 
+                      << cellIdx1<< "("<<ijk1[0]<<","<<ijk1[1]<<","<<ijk1[2]<<")"
+                      << " and " << cellIdx2<<"("<<ijk2[0]<<","<<ijk2[1]<<","<<ijk2[2]<<")");
+        }
+
+        return commonFace;
+    }
+
+
+
+    template<class Grid>
+    inline int PinchProcessor<Grid>::interface_(const Grid& grid,
+                                                const int cellIdx,
+                                                const Opm::FaceDir::DirEnum& faceDir)
+    {
+        const auto actCellIdx = getActiveCellIdx_(grid, cellIdx);
+        const auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
+        const auto cellFacesRange = cell_faces[actCellIdx];
+        int faceIdx = -1;
+        for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
+            int tag = Opm::UgGridHelpers::faceTag(grid, cellFaceIter);
+            if ( (faceDir == Opm::FaceDir::ZMinus && tag == 4) || (faceDir == Opm::FaceDir::ZPlus && tag == 5) ) {
+                faceIdx = *cellFaceIter;
+            }
+        }
+        
+        if (faceIdx == -1) {
+            OPM_THROW(std::logic_error, "Couldn't find the face for cell ." << cellIdx);
+        }
+        
+        return faceIdx;
+    }
+
+
+
+    template<class Grid>
+    inline std::vector<int> PinchProcessor<Grid>::getMinpvCells_(const Grid& grid,
+                                                                 const std::vector<int>& actnum,
+                                                                 const std::vector<double>& pv,
+                                                                 const std::vector<double>& dz)
+    {
+        std::vector<int> minpvCells(pv.size(), 0);
+        for (int idx = 0; idx < static_cast<int>(pv.size()); ++idx) {
+            if (actnum[idx]) {
+                if (pv[idx] < minpvValue_) {
+                    minpvCells[idx] = 1;
+                }
+            }
+        }     
+        return minpvCells;
+    }
+
+
+
+    template<class Grid>
+    inline std::vector<int> PinchProcessor<Grid>::getHfIdxMap_(const Grid& grid)
+    {
+        std::vector<int> hf_ix(2*Opm::UgGridHelpers::numFaces(grid), -1);
+        const auto& f2c = Opm::UgGridHelpers::faceCells(grid);
+        const auto& cf = Opm::UgGridHelpers::cell2Faces(grid);
+        
+        for (int c = 0, i = 0; c < Opm::UgGridHelpers::numCells(grid); ++c) {
+            for (const auto& f: cf[c]) {
+                const auto off = 0 + (f2c(f, 0) != c);
+                hf_ix[2*f + off] = i++;
+            }
+        }
+        return hf_ix;
+    }
+
+
+
+    template<class Grid>
+    inline int PinchProcessor<Grid>::getActiveCellIdx_(const Grid& grid,
+                                                       const int globalIdx)
+    {
+        const int nc = Opm::UgGridHelpers::numCells(grid);
+        const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
+        int idx = -1;
+        for (int i = 0; i < nc; ++i) {
+            if (global_cell[i] == globalIdx) {
+                idx = i;
+                break;
+            }
+        }
+        return idx;
+    }
+
+
+
+    template<class Grid>
+    inline std::vector<double> PinchProcessor<Grid>::transCompute_(const Grid& grid,
+                                                                   const std::vector<double>& htrans,
+                                                                   const std::vector<int>& pinCells,
+                                                                   const std::vector<int>& pinFaces,
+                                                                   const std::vector<double>& multz)       
+    {
+        const int nc = Opm::UgGridHelpers::numCells(grid);
+        const int nf = Opm::UgGridHelpers::numFaces(grid);
+        std::vector<double> trans(nf, 0);
+        int cellFaceIdx = 0;
+        auto cell_faces = Opm::UgGridHelpers::cell2Faces(grid);
+        const auto& hfmap = getHfIdxMap_(grid); 
+        const auto& f2c = Opm::UgGridHelpers::faceCells(grid);
+        for (int cellIdx = 0; cellIdx < nc; ++cellIdx) {
+            auto cellFacesRange = cell_faces[cellIdx];
+            for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter, ++cellFaceIdx) {
+                const int faceIdx = *cellFaceIter;
+                const auto pos = std::find(pinFaces.begin(), pinFaces.end(), faceIdx);
+                if (pos == pinFaces.end()) {
+                    trans[faceIdx] += 1. / htrans[cellFaceIdx];
+                } else {
+                    const int idx1 = std::distance(std::begin(pinFaces), pos);
+                    int idx2;
+                    if (idx1 % 2 == 0) {
+                        idx2 = idx1 + 1;
+                    } else {
+                        idx2 = idx1 - 1;
+                    }
+                    const int f1 = hfmap[2*pinFaces[idx1] + (f2c(pinFaces[idx1], 0) != getActiveCellIdx_(grid, pinCells[idx1]))];
+                    const int f2 = hfmap[2*pinFaces[idx2] + (f2c(pinFaces[idx2], 0) != getActiveCellIdx_(grid, pinCells[idx2]))];
+                    trans[faceIdx] = (1. / htrans[f1] + 1. / htrans[f2]);
+                    trans[pinFaces[idx2]] = trans[faceIdx];
+                }
+            }
+        }
+
+        for (auto f = 0; f < nf; ++f) {
+            trans[f] = 1. / trans[f];
+        }
+
+        return trans;
+    }
+
+
+
+    template<class Grid>
+    inline std::vector<std::vector<int>> PinchProcessor<Grid>::getPinchoutsColumn_(const Grid& grid,
+                                                                                   const std::vector<int>& actnum,
+                                                                                   const std::vector<double>& pv,
+                                                                                   const std::vector<double>& dz)
+    {
+        const int* dims = Opm::UgGridHelpers::cartDims(grid);
+        std::vector<int> minpvCells = getMinpvCells_(grid, actnum, pv, dz);
+        std::vector<std::vector<int>> segment;
+        for (int z = 0; z < dims[2]; ++z) {
+            for (int y = 0; y < dims[1]; ++y) {
+                for (int x = 0; x < dims[0]; ++x) {
+                    const int c = getGlobalIndex_(x, y, z, dims);
+                    std::vector<int> seg;
+                    if (minpvCells[c]) {
+                        seg.push_back(c);
+                        minpvCells[c] = 0;
+                        for (int zz = z+1; zz < dims[2]; ++zz) {
+                            const int cc =  getGlobalIndex_(x, y, zz, dims);
+                            if (minpvCells[cc]) {
+                                seg.push_back(cc);
+                                minpvCells[cc] = 0;
+                            } else {
+                                break;
+                            }
+                        }
+                        segment.push_back(seg);
+                    }
+                }
+            }
+        }
+
+        return segment;
+    }
+
+
+
+    template<class Grid>
+    inline void PinchProcessor<Grid>::transTopbot_(const Grid& grid,
+                                                   const std::vector<double>& htrans,
+                                                   const std::vector<int>& actnum,
+                                                   const std::vector<double>& multz,
+                                                   const std::vector<double>& pv,
+                                                   const std::vector<double>& dz,
+                                                   NNC& nnc)
+    {
+        const int* dims = Opm::UgGridHelpers::cartDims(grid);
+        std::vector<int> pinFaces;
+        std::vector<int> pinCells;
+        std::vector<std::vector<int> > newSeg;
+        auto minpvSeg = getPinchoutsColumn_(grid, actnum, pv, dz);
+        for (auto& seg : minpvSeg) {
+            std::array<int, 3> ijk1 = getCartIndex_(seg.front(), dims);
+            std::array<int, 3> ijk2 = getCartIndex_(seg.back(), dims);
+            auto tmp = seg;
+            if ((ijk1[2]-1) >= 0 && (ijk2[2]+1) < dims[2]) {
+                int topCell = getGlobalIndex_(ijk1[0], ijk1[1], ijk1[2]-1, dims);
+                int botCell = getGlobalIndex_(ijk2[0], ijk2[1], ijk2[2]+1, dims);
+                /// for any segments, we need to find the active top and bottom cells.
+                /// if the original segment's top and bottom is inactive, we need to lookup
+                /// the column until they're found otherwise just ignore this segment.
+                if (!actnum[topCell]) {
+                    seg.insert(seg.begin(), topCell);
+                    for (int topk = ijk1[2]-2; topk > 0; --topk) {
+                        topCell = getGlobalIndex_(ijk1[0], ijk1[1], topk, dims);
+                        if (actnum[topCell]) {
+                            break;
+                        } else {
+                            auto it = seg.begin();
+                            seg.insert(it, topCell);
+                        }
+                    }
+                    pinFaces.push_back(interface_(grid, topCell, Opm::FaceDir::ZPlus));
+                } else {
+                    pinFaces.push_back(interface_(grid, topCell, seg.front()));
+                }
+                tmp.insert(tmp.begin(), topCell);
+                newSeg.push_back(tmp);
+                pinCells.push_back(topCell);
+                if (!actnum[botCell]) {
+                    seg.push_back(botCell);
+                    for (int botk = ijk2[2]+2; botk < dims[2]; ++botk) {
+                        botCell = getGlobalIndex_(ijk2[0], ijk2[1], botk, dims);
+                        if (actnum[botCell]) {
+                            break;
+                        } else {
+                            seg.push_back(botCell);
+                        }
+                    }
+                    pinFaces.push_back(interface_(grid, botCell, Opm::FaceDir::ZMinus));
+                } else {
+                    pinFaces.push_back(interface_(grid, seg.back(), botCell));
+                }
+                pinCells.push_back(botCell);
+            }
+        }
+
+        auto faceTrans = transCompute_(grid, htrans, pinCells, pinFaces, multz);
+        auto multzmap = multzOptions_(grid, pinCells, pinFaces, multz, newSeg);
+        applyMultz_(faceTrans, multzmap);
+        for (int i = 0; i < static_cast<int>(pinCells.size())/2; ++i) {
+            nnc.addNNC(static_cast<int>(pinCells[2*i]), static_cast<int>(pinCells[2*i+1]), faceTrans[pinFaces[2*i]]);
+        }
+    }
+        
+
+ 
+
+    template<class Grid>
+    inline std::unordered_multimap<int, double> PinchProcessor<Grid>::multzOptions_(const Grid& grid,
+                                                                                    const std::vector<int>& pinCells,
+                                                                                    const std::vector<int>& pinFaces,
+                                                                                    const std::vector<double>& multz,
+                                                                                    const std::vector<std::vector<int> >& segs)
+    {
+        std::unordered_multimap<int, double> multzmap;
+        if (multzMode_ == PinchMode::ModeEnum::TOP) {
+            for (int i = 0; i < static_cast<int>(pinFaces.size())/2; ++i) {
+                multzmap.insert(std::make_pair(pinFaces[2*i], multz[getActiveCellIdx_(grid, pinCells[2*i])]));
+                multzmap.insert(std::make_pair(pinFaces[2*i+1],multz[getActiveCellIdx_(grid, pinCells[2*i])]));
+            }
+        } else if (multzMode_ == PinchMode::ModeEnum::ALL) {
+            for (auto& seg : segs) {
+                //find the min multz in seg cells.
+                auto multzValue = std::numeric_limits<double>::max();
+                for (auto& cellIdx : seg) {
+                    auto activeIdx = getActiveCellIdx_(grid, cellIdx);
+                    if (activeIdx != -1) {
+                        multzValue = std::min(multzValue, multz[activeIdx]);
+                    }
+                }
+                //find the right face.
+                auto index = std::distance(std::begin(pinCells), std::find(pinCells.begin(), pinCells.end(), seg.front()));
+                multzmap.insert(std::make_pair(pinFaces[index], multzValue));
+                multzmap.insert(std::make_pair(pinFaces[index+1], multzValue));
+            }
+        }
+
+        return multzmap;
+    }
+
+
+
+    template<class Grid>
+    inline void PinchProcessor<Grid>::applyMultz_(std::vector<double>& trans,
+                                                  const std::unordered_multimap<int, double>& multzmap)
+    {
+        for (auto& x : multzmap) {
+            trans[x.first] *= x.second;
+        }
+    }
+
+
+
+    template<class Grid>
+    inline void PinchProcessor<Grid>::process(const Grid& grid,
+                                              const std::vector<double>& htrans,
+                                              const std::vector<int>& actnum,
+                                              const std::vector<double>& multz,
+                                              const std::vector<double>& pv,
+                                              const std::vector<double>& dz,
+                                              NNC& nnc)
+    {
+        transTopbot_(grid, htrans, actnum, multz, pv, dz, nnc);
+    }
+
+} // namespace Opm
+#endif // OPM_PINCHPROCESSOR_HEADER_INCLUDED

tests/testPinch1.DATA

@@ -0,0 +1,281 @@
+---------------------------------------------------------------------------
+RUNSPEC
+---------------------------------------------------------------------------
+
+TITLE
+    MINPV and PINCH tests
+METRIC
+
+DIMENS
+   4      1      4      /
+
+OIL
+WATER
+GAS
+DISGAS
+VAPOIL
+
+EQLDIMS
+    1   100   10    1    1  /
+
+TABDIMS
+-- NTSFUN     NTPVT   NSSFUN   NPPVT   NTFIP   NRPVT
+   1           1       41      12       1      12 /
+
+WELLDIMS
+--    MAX       CONN              WELLS IN
+--   WELLS     PR WELL   GROUPS     GROUP
+    2            20                 1      2          /
+
+START
+  1 'JAN' 2014  /
+
+UNIFIN
+UNIFOUT
+
+-- Linear solver stack size
+NSTACK
+25 / -- Increased from 10 due to convergence problems
+
+-- Increase the allowable amount of problem prints
+MESSAGES
+-- LIM1 LIM2  LIM3   LIM4  LIM5 LIM6 STOP1 STOP2 STOP3 STOP4 STOP5 STOP6
+-- MESS COM   WARN   PROB  ERR  BUG  MESS   COM  WARN   PROB   ERR   BUG
+    2*      1000000 1000000  4*                 1000000 1000000   /
+
+
+---------------------------------------------------------------------------
+GRID
+---------------------------------------------------------------------------
+DXV
+    4*5 /
+DYV 
+    1*5 /
+DZV
+    4*5 /
+TOPS
+    4*1000 /
+/
+
+-- In this section, the geometry of the simulation grid and the rock
+-- permeabilities and porosities are defined.
+
+INIT
+
+--  ARRAY   VALUE     ------- BOX ------
+EQUALS
+--    DX      5   1 1 1 1 1 4 /
+--    DX      5   2 2 1 1 1 4 /
+--    DX      5   3 4 1 1 1 4 /
+--    DY      5   1 4 1 1 1 4 /
+--    DZ      5   1 4 1 1 1 4 /
+--    TOPS    1000  1  1  1  4  1  1  / top layer
+    PERMX   1000  1 4 1 1 1 4 /
+    PORO 0.3     1 1 1 1 1 4 /
+    PORO 0.0000001   2 2 1 1 2 3 /
+    PORO 0.3      2 2 1 1 1 1 /
+    PORO 0.3      2 2 1 1 4 4 /
+    PORO 0.3      3 4 1 1 1 4 /
+/
+
+EQUALS
+    MULTZ 0.1 1 4 1 1 1 1 /
+    MULTZ 0.2 1 4 1 1 2 3 /
+    MULTZ 0.05 1 4 1 1 4 4 /
+/
+COPY
+    'PERMX'    'PERMY'   /
+    'PERMX'    'PERMZ'   /
+/
+
+PINCH
+------- All default
+0.001 'GAP' 1* 'TOPBOT' /
+/
+
+MINPV
+------- All default
+0.001
+/
+RPTGRID
+ COORD=1 ZCORN=1 TRANX TRANY TRANZ ALLNNC PORV MINPV PINCH DZ NNC /
+
+--RPTINIT
+-- AREAX AREAY AREAZ MULTX MULTY MULTZ MULTX- MULTY- MULTZ- TRANX TRANY TRANZ NTG DZ /
+ 
+
+PROPS     ======
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+         0          1.    1.0000     1.20  /
+        20         40.    1.0120     1.17  /
+        40         80.    1.0255     1.14  /
+        60        120.    1.0380     1.11  /
+        80        160.    1.0510     1.08  /
+       100        200.    1.0630     1.06  /
+       120        240.    1.0750     1.03  /
+       140        280.    1.0870     1.00  /
+       160        320.    1.0985      .98  /
+       180        360.    1.1100      .95  /
+       200        400.    1.1200      .94
+                  500.    1.1189      .94  /
+/
+
+PVTG
+--  Pg     Rv        Bg       Vg
+   100   0.0001       0.010      0.1
+         0.0          0.0104     0.1 /
+   200   0.0004       0.005      0.2
+         0.0          0.0054     0.2 /
+/
+
+
+SWOF
+0.1 0.0 1.0 0.9
+0.2 0.0 0.8 0.8
+0.3 0.1 0.6 0.7
+0.4 0.2 0.4 0.6
+0.7 0.5 0.1 0.3
+0.8 0.6 0.0 0.2
+0.9 0.7 0.0 0.1
+/
+
+SGOF
+0.0 0.0 1.0 0.2
+0.1 0.0 0.7 0.4
+0.2 0.1 0.6 0.6
+0.8 0.7 0.0 2.0
+0.9 1.0 0.0 2.1
+/
+
+PVTW
+--RefPres  Bw      Comp   Vw    Cv
+   1.      1.0   4.0E-5  0.96  0.0 /
+   
+
+ROCK
+--RefPres  Comp
+   1.   5.0E-5 /
+/
+
+DENSITY
+--   Oil          Water        G
+   1000       1000           0.82   /
+
+
+
+---------------------------------------------------------------------------
+REGIONS
+---------------------------------------------------------------------------
+
+EQUALS
+--            VAL        ------- BOX ------
+   SATNUM   1  /
+   PVTNUM   1  / 
+/
+
+
+---------------------------------------------------------------------------
+SOLUTION
+---------------------------------------------------------------------------
+
+-- The solution section defines the initial state of the solution 
+-- variables (phase pressures, saturations, and gas-oil ratios)
+
+-- Equilibration data specification
+--EQUIL
+-- DEPTH   PRES   OW CONTACT
+--  (m)    (bar)     (m)
+--   1000     210      1040 /
+
+-- Initial water saturation for each grid cell
+--EQUALS
+SWAT 
+    16*0.1 /
+SGAS
+    16*0.1 /
+
+PRESSURE 
+    16*200 /
+
+---------------------------------------------------------------------------
+SUMMARY
+---------------------------------------------------------------------------
+
+ALL
+
+
+
+-- This keyword requests that summary data be written to the Summary file 
+-- only at report times. The default, if the RPTONLY keyword is not 
+-- present (and it was not requested in a restart file), is to write the 
+-- summary information at every time step.
+-- Otherwise, the summary file may grow very large.
+RPTONLY
+
+
+
+---------------------------------------------------------------------------
+SCHEDULE
+---------------------------------------------------------------------------
+
+-- The schedule section defines the operations to be simulated
+
+-- Controls on output to the RESTART file
+RPTRST
+  BASIC=3  -- restart data written every FREQ time
+  FREQ=10  -- save every 10th report time
+  PCOW
+  TRANX
+  TRANY
+  TRANZ
+/
+
+WELSPECS
+--     WELL   GROUP      LOCATION     BHP    PI
+--     NAME   NAME      I        J   DEPTH  DEFN
+       PROD    'G'      1       1    1000  'OIL'  /
+       INJE    'G'      1        1    1000  'WATER'  /  
+/
+
+COMPDAT
+--     WELL     -LOCATION- OPEN/ SAT CONN  WELL
+--     NAME      I     J     K1       K2     SHUT  TAB  FACT   DIAM
+       PROD  4    1     1       1    'OPEN'  2*       0.1        /
+       INJE  1     1     1      1    'OPEN'  2*       0.1        / 
+/
+
+
+---------------------------------------------------------------------------
+-- CONTROL: PURE WATER INJECTION
+---------------------------------------------------------------------------
+
+
+WCONINJE
+--      WELL       INJ     OPEN/   CNTL     FLOW       FLOW      BHP
+--      NAME      TYPE     SHUT    MODE     RATE       LIMIT    LIMIT
+--                                        (m^3/day)  (m^3/day)  (bar)
+   INJE     'WATER'  'OPEN'  'RESV'     1*        1.0000e+00  1000000  /           
+/
+
+WCONPROD
+--      WELL     OPEN/  CNTL     OIL  WATER   GAS  LIQU   RES     BHP
+--      NAME     SHUT   MODE    RATE   RATE  RATE  RATE  RATE    LIMIT
+--                                                     (m^3/day)  (bar)
+   PROD    'OPEN'  'RESV'   4*                  1.0000e+00  /             
+/
+
+
+TSTEP
+  1
+/
+
+
+
+---------------------------------------------------------------------------
+END
+---------------------------------------------------------------------------
+
+
+

tests/test_pinchprocessor.cpp

@@ -0,0 +1,112 @@
+/*
+  Copyright 2015 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#if HAVE_DYNAMIC_BOOST_TEST
+#define BOOST_TEST_DYN_LINK
+#endif
+
+#define NVERBOSE // Suppress own messages when throw()in 
+
+#define BOOST_TEST_MODULE PinchProcessorTest
+#include <boost/test/unit_test.hpp>
+#include <vector>
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/core/pressure/tpfa/TransTpfa.hpp>
+#include <opm/core/grid/PinchProcessor.hpp>
+#include <opm/core/props/rock/RockFromDeck.hpp>
+
+
+using namespace Opm;
+BOOST_AUTO_TEST_CASE(Processing)
+{
+    const std::string filename="../tests/testPinch1.DATA";
+    Opm::ParserPtr parser(new Opm::Parser());
+    Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+    Opm::DeckConstPtr deck = parser->parseFile(filename, parseMode);
+    std::shared_ptr<EclipseState> eclstate (new Opm::EclipseState(deck, parseMode));
+    std::vector<double> porv = eclstate->getDoubleGridProperty("PORV")->getData();
+    EclipseGridConstPtr eclgrid = eclstate->getEclipseGrid();
+    Opm::GridManager gridM(eclgrid);
+    typedef UnstructuredGrid Grid;
+    const Grid& grid = *(gridM.c_grid());
+    const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
+    const int* cart_dims = Opm::UgGridHelpers::cartDims(grid);
+    const int nc = Opm::UgGridHelpers::numCells(grid);
+    Opm::RockFromDeck rock;
+    rock.init(eclstate, nc, global_cell, cart_dims);
+
+    const double minpv = eclgrid->getMinpvValue();
+    BOOST_CHECK_EQUAL(minpv, 0.001);
+
+    const double thickness = eclgrid->getPinchThresholdThickness();
+    BOOST_CHECK_EQUAL(thickness, 0.001);
+
+    auto transMode = eclgrid->getPinchOption();
+    BOOST_CHECK_EQUAL(transMode, PinchMode::ModeEnum::TOPBOT);
+
+    auto multzMode = eclgrid->getMultzOption();
+    BOOST_CHECK_EQUAL(multzMode, PinchMode::ModeEnum::TOP);
+
+    PinchProcessor<Grid> pinch(minpv, thickness, transMode, multzMode);
+    std::vector<int> actnum;
+    eclgrid->exportACTNUM(actnum);
+
+    if (actnum.empty() && (nc == eclgrid->getCartesianSize())) {
+        actnum.assign(nc, 1);
+    }
+    std::vector<double> htrans(Opm::UgGridHelpers::numCellFaces(grid));
+    Grid* ug = const_cast<Grid*>(& grid);
+    tpfa_htrans_compute(ug, rock.permeability(), htrans.data());
+    auto transMult = eclstate->getTransMult();
+    std::vector<double> multz(nc, 0.0);
+    std::vector<double> dz(porv.size(), 0.0);
+    for (int i = 0; i < nc; ++i) {
+        multz[i] = transMult->getMultiplier(global_cell[i], Opm::FaceDir::ZPlus);
+    }
+    Opm::NNC nnc(deck, eclgrid);
+    pinch.process(grid, htrans, actnum, multz, porv, dz, nnc);
+    auto nnc1 = nnc.nnc1();
+    auto nnc2 = nnc.nnc2();
+    auto trans = nnc.trans();
+
+    BOOST_CHECK(nnc.hasNNC());
+    BOOST_CHECK_EQUAL(nnc.numNNC(), 1);
+   
+    auto nnc1_index = 1 + cart_dims[0] * (0 + cart_dims[1] * 0);
+    auto nnc2_index = 1 + cart_dims[0] * (0 + cart_dims[1] * 3);
+    BOOST_CHECK_EQUAL(nnc1[0], nnc1_index);
+    BOOST_CHECK_EQUAL(nnc2[0], nnc2_index);
+    double factor = Opm::prefix::centi*Opm::unit::Poise 
+        * Opm::unit::cubic(Opm::unit::meter)
+        / Opm::unit::day
+        / Opm::unit::barsa;
+    for (auto& tran : trans) {
+        tran = unit::convert::to(tran, factor);
+    }
+    BOOST_CHECK(std::fabs(trans[0]-4.26350022) < 1e-3);
+}