ebos: implement threshold pressures

the implementation also supports to calculate the threshold pressures
based on the initial condition (and it should work even for
MPI-parallel runs).

applications/ebos/eclfluxmodule.hh

@@ -137,6 +137,15 @@ public:
                   "The ECL transmissibility module does not provide explicit potential gradients");
     }
 
+    /*!
+     * \brief Return the gravity corrected pressure difference between the interior and
+     *        the exterior of a face.
+     *
+     * \param phaseIdx The index of the fluid phase
+     */
+    const Evaluation& pressureDifferential(unsigned phaseIdx) const
+    { return pressureDifferential_[phaseIdx]; }
+
     /*!
      * \brief Return the filter velocity of a fluid phase at the face's integration point
      *        [m/s]
@@ -235,6 +244,19 @@ protected:
 
             pressureDifferential_[phaseIdx] = pressureExterior - pressureInterior;
 
+            // apply the threshold pressure for the intersection. note that the concept
+            // of threshold pressure is a quite big hack that only makes sense for ECL
+            // datasets. (and even there its physical justification is quite
+            // questionable IMO.)
+            if (std::abs(Toolbox::value(pressureDifferential_[phaseIdx])) > thpres_)
+                pressureDifferential_[phaseIdx] -=
+                    Ewoms::signum(pressureDifferential_[phaseIdx])*thpres_;
+            else {
+                pressureDifferential_[phaseIdx] = 0.0;
+                volumeFlux_[phaseIdx] = 0.0;
+                continue;
+            }
+
             // this is slightly hacky because in the automatic differentiation case, it
             // only works for the element centered finite volume method. for ebos this
             // does not matter, though.
@@ -266,6 +288,9 @@ protected:
     // the area of the face between the DOFs [m^2]
     Scalar faceArea_;
 
+    // threshold pressure [Pa]
+    Scalar thpres_;
+
     // the volumetric flux of all phases [m^3/s]
     Evaluation volumeFlux_[numPhases];
 

applications/ebos/eclproblem.hh

@@ -52,6 +52,7 @@
 #include "eclsummarywriter.hh"
 #include "ecloutputblackoilmodule.hh"
 #include "ecltransmissibility.hh"
+#include "eclthresholdpressure.hh"
 #include "ecldummygradientcalculator.hh"
 #include "eclfluxmodule.hh"
 #include "ecldeckunits.hh"
@@ -137,7 +138,7 @@ SET_BOOL_PROP(EclBaseProblem, EnableWriteAllSolutions, false);
 //
 // By default, stop it after the universe will probably have stopped
 // to exist. (the ECL problem will finish the simulation explicitly
-// after it simulated the last episode specified in the deck->)
+// after it simulated the last episode specified in the deck.)
 SET_SCALAR_PROP(EclBaseProblem, EndTime, 1e100);
 
 // The default for the initial time step size of the simulation [s].
@@ -259,6 +260,7 @@ public:
     EclProblem(Simulator &simulator)
         : ParentType(simulator)
         , transmissibilities_(simulator)
+        , thresholdPressures_(simulator)
         , wellManager_(simulator)
         , deckUnits_(simulator)
         , eclWriter_(simulator)
@@ -528,11 +530,17 @@ public:
     { return intrinsicPermeability_[globalElemIdx]; }
 
     /*!
-     * \copydoc FvBaseMultiPhaseProblem::transmissibility
+     * \copydoc BlackOilBaseProblem::transmissibility
      */
     Scalar transmissibility(unsigned elem1Idx, unsigned elem2Idx) const
     { return transmissibilities_.transmissibility(elem1Idx, elem2Idx); }
 
+    /*!
+     * \copydoc BlackOilBaseProblem::thresholdPressure
+     */
+    Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const
+    { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
+
     /*!
      * \copydoc FvBaseMultiPhaseProblem::porosity
      */
@@ -679,6 +687,11 @@ public:
 
         // update the data required for capillary pressure hysteresis
         updateHysteresis_();
+
+        // let the object for threshold pressures initialize itself. this is done only at
+        // this point, because determining the threshold pressures may require to access
+        // the initial solution.
+        thresholdPressures_.finishInit();
     }
 
     /*!
@@ -767,7 +780,7 @@ private:
         ////////////////////////////////
         // permeability
 
-        // read the intrinsic permeabilities from the eclState-> Note that all arrays
+        // read the intrinsic permeabilities from the eclState. Note that all arrays
         // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
         // dune-cornerpoint grid object might remove a few elements...
         if (eclState->hasDeckDoubleGridProperty("PERMX")) {
@@ -934,11 +947,11 @@ private:
 
         initialFluidStates_.resize(numDof);
 
-        const std::vector<double> &waterSaturationData =
+        const std::vector<double>& waterSaturationData =
             deck->getKeyword("SWAT").getSIDoubleData();
-        const std::vector<double> &gasSaturationData =
+        const std::vector<double>& gasSaturationData =
             deck->getKeyword("SGAS").getSIDoubleData();
-        const std::vector<double> &pressureData =
+        const std::vector<double>& pressureData =
             deck->getKeyword("PRESSURE").getSIDoubleData();
         const std::vector<double> *rsData = 0;
         if (enableDisgas)
@@ -1104,6 +1117,8 @@ private:
 
     std::shared_ptr<EclMaterialLawManager> materialLawManager_;
 
+    EclThresholdPressure<TypeTag> thresholdPressures_;
+
     std::vector<unsigned short> rockTableIdx_;
     std::vector<RockParams> rockParams_;
 

applications/ebos/eclthresholdpressure.hh

@@ -0,0 +1,299 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright (C) 2015 by Andreas Lauser
+
+  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 2 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/>.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Ewoms::EclThresholdPressure
+ */
+#ifndef EWOMS_ECL_THRESHOLD_PRESSURE_HH
+#define EWOMS_ECL_THRESHOLD_PRESSURE_HH
+
+#include <ewoms/common/propertysystem.hh>
+
+#include <opm/material/localad/Evaluation.hpp>
+#include <opm/material/localad/Math.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
+
+#include <dune/grid/common/gridenums.hh>
+#include <dune/common/version.hh>
+
+#include <array>
+#include <vector>
+#include <unordered_map>
+
+namespace Ewoms {
+namespace Properties {
+NEW_PROP_TAG(Simulator);
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(Evaluation);
+NEW_PROP_TAG(ElementContext);
+NEW_PROP_TAG(FluidSystem);
+}
+
+/*!
+ * \ingroup EclBlackOilSimulator
+ *
+ * \brief This class calculates the threshold pressure for grid faces according to the
+ *        Eclipse Reference Manual.
+ *
+ * If the difference of the pressure potential between two cells is below the threshold
+ * pressure, the pressure potential difference is assumed to be zero, if it is larger
+ * than the threshold pressure, it is reduced by the threshold pressure.
+ */
+template <class TypeTag>
+class EclThresholdPressure
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+    enum { numPhases = FluidSystem::numPhases };
+
+public:
+    EclThresholdPressure(const Simulator& simulator)
+        : simulator_(simulator)
+    {
+        enableThresholdPressure_ = false;
+    }
+
+    /*!
+     * \brief Actually compute the threshold pressures over a face as a pre-compute step.
+     */
+    void finishInit()
+    {
+        const auto& gridView = simulator_.gridView();
+        auto deck = simulator_.gridManager().deck();
+
+        unsigned numElements = gridView.size(/*codim=*/0);
+
+        // this code assumes that the DOFs are the elements. (i.e., an
+        // ECFV spatial discretization with TPFA). if you try to use
+        // it with something else, you're currently out of luck,
+        // sorry!
+        assert((int) simulator_.model().numGridDof() == numElements);
+
+        const auto& gridManager = simulator_.gridManager();
+        Opm::EclipseStateConstPtr eclState = gridManager.eclState();
+        Opm::SimulationConfigConstPtr simConfig = eclState->getSimulationConfig();
+
+        enableThresholdPressure_ = simConfig->hasThresholdPressure();
+        if (!enableThresholdPressure_)
+            return;
+
+        numEquilRegions_ =
+            deck->getKeyword("EQLDIMS").getRecord(0).getItem("NTEQUL").template get<int>(0);
+        if (numEquilRegions_ > 0xff) {
+            // make sure that the index of an equilibration region can be stored in a
+            // single byte
+            OPM_THROW(std::runtime_error,
+                      "The maximum number of supported equilibration regions is 255!");
+        }
+
+        // allocate the array which specifies the threshold pressures
+        thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+        thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+
+        // internalize the data specified using the EQLNUM keyword
+        const std::vector<int>& equilRegionData =
+            eclState->getIntGridProperty("EQLNUM")->getData();
+        elemEquilRegion_.resize(numElements, 0);
+        for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
+            int cartElemIdx = gridManager.cartesianIndex(elemIdx);
+
+            // ECL uses Fortran-style indices but we want C-style ones!
+            elemEquilRegion_[elemIdx] = equilRegionData[cartElemIdx] - 1;
+        }
+
+        computeDefaultThresholdPressures_();
+        applyExplicitThresholdPressures_();
+    }
+
+    /*!
+     * \brief Returns the theshold pressure [Pa] for the intersection between two elements.
+     *
+     * This is tailor made for the E100 threshold pressure mechanism and it is thus quite
+     * a hack: First of all threshold pressures in general are unphysical, and second,
+     * they should be different for the fluid phase but are not. Anyway, this seems to be
+     * E100's way of doing things, so we do it the same way.
+     */
+    Scalar thresholdPressure(int elemIdx1, int elemIdx2) const
+    {
+        if (!enableThresholdPressure_)
+            return 0.0;
+
+        unsigned short equilRegion1Idx = elemEquilRegion_[elemIdx1];
+        unsigned short equilRegion2Idx = elemEquilRegion_[elemIdx2];
+
+        if (equilRegion1Idx == equilRegion2Idx)
+            return 0.0;
+
+        return thpres_[equilRegion1Idx*numEquilRegions_ + equilRegion2Idx];
+    }
+
+private:
+    // compute the defaults of the threshold pressures using the initial condition
+    void computeDefaultThresholdPressures_()
+    {
+        const auto& gridManager = simulator_.gridManager();
+        const auto& gridView = gridManager.gridView();
+
+        typedef Opm::MathToolbox<Evaluation> Toolbox;
+        // loop over the whole grid and compute the maximum gravity adjusted pressure
+        // difference between two EQUIL regions.
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+        ElementContext elemCtx(simulator_);
+        for (; elemIt != elemEndIt; ++elemIt) {
+
+            const auto& elem = *elemIt;
+            if (elem.partitionType() != Dune::InteriorEntity)
+                continue;
+
+            elemCtx.updateAll(elem);
+            const auto &stencil = elemCtx.stencil(/*timeIdx=*/0);
+
+            for (int scvfIdx = 0; scvfIdx < stencil.numInteriorFaces(); ++ scvfIdx) {
+                const auto &face = stencil.interiorFace(scvfIdx);
+
+                unsigned i = face.interiorIndex();
+                unsigned j = face.exteriorIndex();
+
+                unsigned insideElemIdx = elemCtx.globalSpaceIndex(i, /*timeIdx=*/0);
+                unsigned outsideElemIdx = elemCtx.globalSpaceIndex(j, /*timeIdx=*/0);
+
+                unsigned equilRegionInside = elemEquilRegion_[insideElemIdx];
+                unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx];
+
+                if (equilRegionInside == equilRegionOutside)
+                    continue;
+
+                // determine the maximum difference of the pressure of any phase over the
+                // intersection
+                Scalar pth = 0;
+                const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0);
+                for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    pth = std::max(pth, std::abs(Toolbox::value(extQuants.pressureDifferential(phaseIdx))));
+
+                int offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
+                int offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
+                thpresDefault_[offset1] = std::max(thpres_[offset1], pth);
+                thpresDefault_[offset2] = std::max(thpres_[offset2], pth);
+            }
+        }
+
+        // make sure that the threshold pressures is consistent for parallel
+        // runs. (i.e. take the maximum of all processes)
+        for (int i = 0; i < thpres_.size(); ++i)
+            thpres_[i] = gridView.comm().max(thpres_[i]);
+    }
+
+    // internalize the threshold pressures which where explicitly specified via the
+    // THPRES keyword.
+    void applyExplicitThresholdPressures_()
+    {
+        const auto& gridManager = simulator_.gridManager();
+        const auto& gridView = gridManager.gridView();
+        const auto& elementMapper = simulator_.model().elementMapper();
+        const auto& eclState = simulator_.gridManager().eclState();
+        Opm::SimulationConfigConstPtr simConfig = eclState->getSimulationConfig();
+        Opm::ThresholdPressureConstPtr thpres = simConfig->getThresholdPressure();
+
+        // set the threshold pressures for all EQUIL region boundaries which have a
+        // intersection in the grid
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+
+            const auto& elem = *elemIt;
+            if (elem.partitionType() != Dune::InteriorEntity)
+                continue;
+
+            auto isIt = gridView.ibegin(elem);
+            const auto& isEndIt = gridView.iend(elem);
+            for (; isIt != isEndIt; ++ isIt) {
+                // store intersection, this might be costly
+                const auto& intersection = *isIt;
+
+                // ignore boundary intersections for now (TODO?)
+                if (intersection.boundary())
+                    continue;
+
+                const auto& inside = intersection.inside();
+                const auto& outside = intersection.outside();
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+                unsigned insideElemIdx = elementMapper.index(inside);
+                unsigned outsideElemIdx = elementMapper.index(outside);
+#else
+                unsigned insideElemIdx = elementMapper.map(*inside);
+                unsigned outsideElemIdx = elementMapper.map(*outside);
+#endif
+
+                unsigned equilRegionInside = elemEquilRegion_[insideElemIdx];
+                unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx];
+                if (thpres->hasRegionBarrier(equilRegionInside + 1, equilRegionOutside + 1) ||
+                    thpres->hasRegionBarrier(equilRegionOutside + 1, equilRegionInside + 1))
+                {
+                    Scalar pth = 0.0;
+                    if (thpres->hasThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1)) {
+                        // threshold pressure explicitly specified
+                        pth = thpres->getThresholdPressure(equilRegionInside + 1, equilRegionOutside + 1);
+                    }
+                    else if (thpres->hasThresholdPressure(equilRegionOutside + 1, equilRegionInside + 1)) {
+                        // threshold pressure explicitly specified
+                        pth = thpres->getThresholdPressure(equilRegionOutside + 1, equilRegionInside + 1);
+                    }
+                    else {
+                        // take the threshold pressure from the initial condition
+                        unsigned offset = equilRegionInside*numEquilRegions_ + equilRegionOutside;
+                        pth = thpresDefault_[offset];
+                    }
+
+                    unsigned offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside;
+                    unsigned offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;
+
+                    thpres_[offset1] = pth;
+                    thpres_[offset2] = pth;
+                }
+            }
+        }
+    }
+
+    const Simulator& simulator_;
+
+    std::vector<Scalar> thpresDefault_;
+    std::vector<Scalar> thpres_;
+    unsigned numEquilRegions_;
+    std::vector<unsigned char> elemEquilRegion_;
+
+    bool enableThresholdPressure_;
+};
+
+} // namespace Ewoms
+
+#endif