From e8b460fe109697a51528b6c16e03a55853a30eb3 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Tue, 9 Aug 2016 12:40:33 +0200
Subject: [PATCH] implement MINPV treatment

this is done my modifying the NTG values. it is just a
straight-forward port of the corresponding functionality from
opm-simulators.
---
 applications/ebos/ecltransmissibility.hh | 42 +++++++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/applications/ebos/ecltransmissibility.hh b/applications/ebos/ecltransmissibility.hh
index 0d087bf7b..c94c83894 100644
--- a/applications/ebos/ecltransmissibility.hh
+++ b/applications/ebos/ecltransmissibility.hh
@@ -100,9 +100,11 @@ public:
         const auto eclGrid = eclState->getInputGrid();
         const auto& transMult = eclState->getTransMult();
 
-        const std::vector<double>& ntg =
+        std::vector<double> ntg =
             eclState->get3DProperties().getDoubleGridProperty("NTG").getData();
 
+        applyMinPv_(ntg);
+
         unsigned numElements = elementMapper.size();
 
         // this code assumes that the DOFs are the elements. (i.e., an
@@ -329,6 +331,44 @@ private:
         }
     }
 
+    void applyMinPv_(std::vector<double> &ntg)
+    {
+        const auto& model = simulator_.model();
+        const auto& gridManager = simulator_.gridManager();
+        const auto& eclState = gridManager.eclState();
+        const Opm::EclipseGridConstPtr eclGrid = gridManager.eclGrid();
+        const auto& porv = eclState->get3DProperties().getDoubleGridProperty("PORV").getData();
+
+        Scalar minpvValue = eclGrid->getMinpvValue();
+        if (minpvValue <= 0.0)
+            return;
+
+        int numCells = model.numGridDof();
+        for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+            const int nx = eclGrid->getNX();
+            const int ny = eclGrid->getNY();
+            const int cartesianCellIdx = gridManager.cartesianIndex(cellIdx);
+
+            const Scalar cellVolume = eclGrid->getCellVolume(cartesianCellIdx);
+            ntg[cartesianCellIdx] *= cellVolume;
+            Scalar totalCellVolume = cellVolume;
+
+            // Average properties as long as there exist cells above which exhibit a pore
+            // volume less than the MINPV threshold
+            int cartesianCellIdxAbove = cartesianCellIdx - nx*ny;
+            while (cartesianCellIdxAbove >= 0
+                   && porv[cartesianCellIdxAbove] > 0.0
+                   && porv[cartesianCellIdxAbove] < minpvValue)
+            {
+                const double cellVolumeAbove = eclGrid->getCellVolume(cartesianCellIdxAbove);
+                totalCellVolume += cellVolumeAbove;
+                ntg[cartesianCellIdx] += ntg[cartesianCellIdxAbove]*cellVolumeAbove;
+                cartesianCellIdxAbove -= nx*ny;
+            }
+            ntg[cartesianCellIdx] /= totalCellVolume;
+        }
+    }
+
     void applyNtg_(Scalar &trans, unsigned faceIdx, unsigned cartElemIdx,
                    const std::vector<double>& ntg) const
     {