add option to use different flux modules and add one for ECL-transmissibilities to ebos

the ebos module implemenents what Eclipse calls 'NEWTRAN'
transmissibilities. Also, this commit required a few cleanups in the
velocity module infrastructure.

applications/ebos/ecldummygradientcalculator.hh

@@ -0,0 +1,102 @@
+/*
+  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::EclDummyGradientCalculator
+ */
+#ifndef EWOMS_ECL_DUMMY_GRADIENT_CALCULATOR_HH
+#define EWOMS_ECL_DUMMY_GRADIENT_CALCULATOR_HH
+
+#include <ewoms/disc/common/fvbaseproperties.hh>
+
+#include <dune/common/fvector.hh>
+
+namespace Ewoms {
+/*!
+ * \ingroup EclBlackOilSimulator
+ *
+ * \brief This is a "dummy" gradient calculator which does not do anything.
+ *
+ * The ECL blackoil simulator does not need any gradients: Volume fluxes are calculated
+ * via pressure differences instead of pressure gradients (i.e., transmissibilities
+ * instead of permeabilities), and an energy equation and molecular diffusion are not
+ * supported.
+ */
+template<class TypeTag>
+class EclDummyGradientCalculator
+{
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+
+public:
+    static void registerParameters()
+    { }
+
+    template <bool prepareValues = true, bool prepareGradients = true>
+    void prepare(const ElementContext &elemCtx, int timeIdx)
+    { }
+
+    template <class QuantityCallback, class QuantityType = Scalar>
+    QuantityType calculateValue(const ElementContext &elemCtx,
+                                int fapIdx,
+                                const QuantityCallback &quantityCallback) const
+    {
+        OPM_THROW(std::logic_error,
+                  "Generic values are not supported by the ECL black-oil simulator");
+    }
+
+    template <class QuantityCallback>
+    void calculateGradient(DimVector &quantityGrad,
+                           const ElementContext &elemCtx,
+                           int fapIdx,
+                           const QuantityCallback &quantityCallback) const
+    {
+        OPM_THROW(std::logic_error,
+                  "Generic gradients are not supported by the ECL black-oil simulator");
+    }
+
+    template <class QuantityCallback>
+    Scalar calculateBoundaryValue(const ElementContext &elemCtx,
+                                  int fapIdx,
+                                  const QuantityCallback &quantityCallback)
+    {
+        OPM_THROW(std::logic_error,
+                  "Generic boundary values are not supported by the ECL black-oil simulator");
+    }
+
+    template <class QuantityCallback>
+    void calculateBoundaryGradient(DimVector &quantityGrad,
+                                   const ElementContext &elemCtx,
+                                   int fapIdx,
+                                   const QuantityCallback &quantityCallback) const
+    {
+        OPM_THROW(std::logic_error,
+                  "Generic boundary gradients are not supported by the ECL black-oil simulator");
+    }
+};
+} // namespace Ewoms
+
+#endif

applications/ebos/eclfluxmodule.hh

@@ -0,0 +1,265 @@
+/*
+  Copyright (C) 2014 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
+ *
+ * \brief This file contains the flux module which is used for ECL problems
+ *        two-point flux approximation
+ *
+ * This is used by the ECL blackoil simulator
+ */
+#ifndef EWOMS_ECL_FLUX_MODULE_HH
+#define EWOMS_ECL_FLUX_MODULE_HH
+
+#include <ewoms/disc/common/fvbaseproperties.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+namespace Opm {
+namespace Properties {
+NEW_PROP_TAG(MaterialLaw);
+}}
+
+namespace Ewoms {
+template <class TypeTag>
+class EclTransIntensiveQuantities;
+
+template <class TypeTag>
+class EclTransExtensiveQuantities;
+
+template <class TypeTag>
+class EclTransBaseProblem;
+
+/*!
+ * \ingroup EclTransmissibility
+ * \brief Specifies a velocity module which uses the transmissibilities.
+ */
+template <class TypeTag>
+struct EclTransVelocityModule
+{
+    typedef EclTransIntensiveQuantities<TypeTag> VelocityIntensiveQuantities;
+    typedef EclTransExtensiveQuantities<TypeTag> VelocityExtensiveQuantities;
+    typedef EclTransBaseProblem<TypeTag> VelocityBaseProblem;
+
+    /*!
+     * \brief Register all run-time parameters for the velocity module.
+     */
+    static void registerParameters()
+    { }
+};
+
+/*!
+ * \ingroup EclTransmissibility
+ * \brief Provides the defaults for the parameters required by the
+ *        transmissibility based volume flux calculation.
+ */
+template <class TypeTag>
+class EclTransBaseProblem
+{ };
+
+/*!
+ * \ingroup EclTransmissibility
+ * \brief Provides the intensive quantities for the Darcy velocity module
+ */
+template <class TypeTag>
+class EclTransIntensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+protected:
+    void update_(const ElementContext &elemCtx, int dofIdx, int timeIdx)
+    { }
+};
+
+/*!
+ * \ingroup EclTransmissibility
+ * \brief Provides the ECL "velocity module"
+ */
+template <class TypeTag>
+class EclTransExtensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+
+    enum { dimWorld = GridView::dimensionworld };
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
+
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
+
+public:
+    /*!
+     * \brief Returns transmissibility for a given sub-control volume face.
+     */
+    Scalar transmissibility() const
+    { return trans_; }
+
+    /*!
+     * \brief Return the intrinsic permeability tensor at a face [m^2]
+     */
+    const DimMatrix& intrinsicPermeability() const
+    {
+        OPM_THROW(Opm::NotImplemented,
+                  "The ECL transmissibility module does not provide an explicit intrinsic permeability");
+    }
+
+    /*!
+     * \brief Return the pressure potential gradient of a fluid phase at the
+     *        face's integration point [Pa/m]
+     *
+     * \param phaseIdx The index of the fluid phase
+     */
+    const DimVector& potentialGrad(int phaseIdx) const
+    {
+        OPM_THROW(Opm::NotImplemented,
+                  "The ECL transmissibility module does not provide explicit potential gradients");
+    }
+
+    /*!
+     * \brief Return the filter velocity of a fluid phase at the
+     *        face's integration point [m/s]
+     *
+     * \param phaseIdx The index of the fluid phase
+     */
+    const DimVector& filterVelocity(int phaseIdx) const
+    {
+        OPM_THROW(Opm::NotImplemented,
+                  "The ECL transmissibility module does not provide explicit filter velocities");
+    }
+
+    /*!
+     * \brief Return the volume flux of a fluid phase at the face's integration point
+     *        \f$[m^3/s / m^2]\f$
+     *
+     * This is the fluid volume of a phase per second and per square meter of face
+     * area.
+     *
+     * \param phaseIdx The index of the fluid phase
+     */
+    Scalar volumeFlux(int phaseIdx) const
+    { return - pressureDifferential_[phaseIdx]*mobility_[phaseIdx] * trans_/faceArea_; }
+
+protected:
+    /*!
+     * \brief Returns the local index of the degree of freedom in which is
+     *        in upstream direction.
+     *
+     * i.e., the DOF which exhibits a higher effective pressure for
+     * the given phase.
+     */
+    int upstreamIndex_(int phaseIdx) const
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+        return (pressureDifferential_[phaseIdx] >= 0)?exteriorDofIdx_:interiorDofIdx_;
+    }
+
+    /*!
+     * \brief Returns the local index of the degree of freedom in which is
+     *        in downstream direction.
+     *
+     * i.e., the DOF which exhibits a lower effective pressure for the
+     * given phase.
+     */
+    int downstreamIndex_(int phaseIdx) const
+    {
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+        return (pressureDifferential_[phaseIdx] >= 0)?interiorDofIdx_:exteriorDofIdx_;
+    }
+
+    /*!
+     * \brief Update the required gradients for interior faces
+     */
+    void calculateGradients_(const ElementContext &elemCtx, int scvfIdx, int timeIdx)
+    {
+        Valgrind::SetUndefined(*this);
+
+        const auto& problem = elemCtx.problem();
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        interiorDofIdx_ = scvf.interiorIndex();
+        exteriorDofIdx_ = scvf.exteriorIndex();
+        assert(interiorDofIdx_ != exteriorDofIdx_);
+
+        trans_ = problem.transmissibility(stencil.globalSpaceIndex(interiorDofIdx_),
+                                          stencil.globalSpaceIndex(exteriorDofIdx_));
+        faceArea_ = scvf.area();
+
+        // estimate the gravity correction: for performance reasons we use a simplified
+        // approach for this flux module that assumes that gravity is constant and always
+        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
+        Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
+
+        const auto &intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx_, timeIdx);
+        const auto &intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx_, timeIdx);
+
+        Scalar zIn = elemCtx.pos(interiorDofIdx_, timeIdx)[dimWorld - 1];
+        Scalar zEx = elemCtx.pos(exteriorDofIdx_, timeIdx)[dimWorld - 1];
+        Scalar zFace = scvf.integrationPos()[dimWorld - 1];
+
+        Scalar distZIn = zIn - zFace;
+        Scalar distZEx = zEx - zFace;
+
+        for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
+            // calculate the hydrostatic pressures at the face's integration point
+            Scalar rhoIn = intQuantsIn.fluidState().density(phaseIdx);
+            Scalar rhoEx = intQuantsEx.fluidState().density(phaseIdx);
+
+            Scalar pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
+            Scalar pressureExterior = intQuantsEx.fluidState().pressure(phaseIdx);
+
+            pressureInterior += - rhoIn*(g*distZIn);
+            pressureExterior += - rhoEx*(g*distZEx);
+
+            pressureDifferential_[phaseIdx] = pressureExterior - pressureInterior;
+
+            const auto& up = elemCtx.intensiveQuantities(upstreamIndex_(phaseIdx), timeIdx);
+            mobility_[phaseIdx] = up.mobility(phaseIdx);
+        }
+    }
+
+    /*!
+     * \brief Update the velocities for all fluid phases on the interior faces of the context
+     */
+    void calculateVelocities_(const ElementContext &elemCtx, int scvfIdx, int timeIdx)
+    { }
+
+    // the local indices of the interior and exterior degrees of freedom
+    int interiorDofIdx_;
+    int exteriorDofIdx_;
+
+    // transmissibility [m^3 s]
+    Scalar trans_;
+
+    // the area of the face between the DOFs [m^2]
+    Scalar faceArea_;
+
+    // the mobility of all phases [1 / (Pa s)]
+    Scalar mobility_[numPhases];
+
+    // the difference in effective pressure between the two degrees of
+    // freedom [Pa]
+    Scalar pressureDifferential_[numPhases];
+};
+
+} // namespace Ewoms
+
+#endif

applications/ebos/eclproblem.hh

@@ -29,6 +29,9 @@
 #include "eclwriter.hh"
 #include "eclsummarywriter.hh"
 #include "ecloutputblackoilmodule.hh"
+#include "ecltransmissibility.hh"
+#include "ecldummygradientcalculator.hh"
+#include "eclfluxmodule.hh"
 
 #include <ewoms/models/blackoil/blackoilmodel.hh>
 #include <ewoms/disc/ecfv/ecfvdiscretization.hh>
@@ -154,6 +157,12 @@ SET_BOOL_PROP(EclBaseProblem, EnableEclSummaryOutput, true);
 // decent speedup...
 SET_BOOL_PROP(EclBaseProblem, EnableIntensiveQuantityCache, true);
 
+// Use the "velocity module" which uses the Eclipse "NEWTRAN" transmissibilities
+SET_TYPE_PROP(EclBaseProblem, VelocityModule, Ewoms::EclTransVelocityModule<TypeTag>);
+
+// Use the dummy gradient calculator in order not to do unnecessary work.
+SET_TYPE_PROP(EclBaseProblem, GradientCalculator, Ewoms::EclDummyGradientCalculator<TypeTag>);
+
 // The default name of the data file to load
 SET_STRING_PROP(EclBaseProblem, GridFile, "data/ecl.DATA");
 }} // namespace Properties, Opm
@@ -230,6 +239,7 @@ public:
      */
     EclProblem(Simulator &simulator)
         : ParentType(simulator)
+        , transmissibilities_(simulator)
         , wellManager_(simulator)
         , eclWriter_(simulator)
         , summaryWriter_(simulator)
@@ -253,9 +263,15 @@ public:
         // (z coodinates represent depth, not height.)
         this->gravity_[dim - 1] *= -1;
 
+        // the "NOGRAV" keyword from Frontsim disables gravity...
+        const auto& deck = simulator.gridManager().deck();
+        if (deck->hasKeyword("NOGRAV"))
+            this->gravity_ = 0.0;
+
         initFluidSystem_();
         readRockParameters_();
         readMaterialParameters_();
+        transmissibilities_.finishInit();
         readInitialCondition_();
 
         // initialize the wells. Note that this needs to be done after initializing the
@@ -401,27 +417,19 @@ public:
     }
 
     /*!
-     * \copydoc FvBaseMultiPhaseProblem::intersectionIntrinsicPermeability
+     * \brief This method returns the intrinsic permeability tensor
+     *        given a global element index.
+     *
+     * Its main (only?) usage is the ECL transmissibility calculation code...
      */
-    template <class Context>
+    const DimMatrix &intrinsicPermeability(int globalElemIdx) const
-    void intersectionIntrinsicPermeability(DimMatrix &result,
+    { return intrinsicPermeability_[globalElemIdx]; }
-                                           const Context &context,
-                                           int localIntersectionIdx, int timeIdx) const
-    {
-        // calculate the intersection index
-        const auto &scvf = context.stencil(timeIdx).interiorFace(localIntersectionIdx);
 
-        int numElements = this->model().numGridDof();
+    /*!
-
+     * \copydoc FvBaseMultiPhaseProblem::transmissibility
-        size_t interiorElemIdx = context.globalSpaceIndex(scvf.interiorIndex(), timeIdx);
+     */
-        size_t exteriorElemIdx = context.globalSpaceIndex(scvf.exteriorIndex(), timeIdx);
+    Scalar transmissibility(int elem1Idx, int elem2Idx) const
-
+    { return transmissibilities_.transmissibility(elem1Idx, elem2Idx); }
-        size_t elem1Idx = std::min(interiorElemIdx, exteriorElemIdx);
-        size_t elem2Idx = std::max(interiorElemIdx, exteriorElemIdx);
-
-        size_t globalIntersectionIdx = elem1Idx*numElements + elem2Idx;
-        result = intersectionIntrinsicPermeability_.at(globalIntersectionIdx);
-    }
 
     /*!
      * \copydoc FvBaseMultiPhaseProblem::porosity
@@ -668,20 +676,6 @@ private:
             OPM_THROW(std::logic_error,
                       "Can't read the intrinsic permeability from the ecl state. "
                       "(The PERM{X,Y,Z} keywords are missing)");
-
-        // apply the NTG keyword to the X and Y permeabilities
-        if (eclState->hasDoubleGridProperty("NTG")) {
-            const std::vector<double> &ntgData =
-                eclState->getDoubleGridProperty("NTG")->getData();
-
-            for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
-                int cartesianElemIdx = grid.globalCell()[dofIdx];
-                intrinsicPermeability_[dofIdx][0][0] *= ntgData[cartesianElemIdx];
-                intrinsicPermeability_[dofIdx][1][1] *= ntgData[cartesianElemIdx];
-            }
-        }
-
-        computeFaceIntrinsicPermeabilities_();
         ////////////////////////////////
 
 
@@ -963,141 +957,9 @@ private:
         }
     }
 
-    void computeFaceIntrinsicPermeabilities_()
-    {
-        auto eclState = this->simulator().gridManager().eclState();
-        const auto &grid = this->simulator().gridManager().grid();
-
-        int numElements = this->gridView().size(/*codim=*/0);
-
-        std::vector<double> multx(numElements, 1.0);
-        std::vector<double> multy(numElements, 1.0);
-        std::vector<double> multz(numElements, 1.0);
-        std::vector<double> multxMinus(numElements, 1.0);
-        std::vector<double> multyMinus(numElements, 1.0);
-        std::vector<double> multzMinus(numElements, 1.0);
-
-        // retrieve the transmissibility multiplier keywords. Note that we use them as
-        // permeability multipliers...
-        if (eclState->hasDoubleGridProperty("MULTX"))
-            multx = eclState->getDoubleGridProperty("MULTX")->getData();
-        if (eclState->hasDoubleGridProperty("MULTX-"))
-            multxMinus = eclState->getDoubleGridProperty("MULTX-")->getData();
-        if (eclState->hasDoubleGridProperty("MULTY"))
-            multy = eclState->getDoubleGridProperty("MULTY")->getData();
-        if (eclState->hasDoubleGridProperty("MULTY-"))
-            multyMinus = eclState->getDoubleGridProperty("MULTY-")->getData();
-        if (eclState->hasDoubleGridProperty("MULTZ"))
-            multz = eclState->getDoubleGridProperty("MULTZ")->getData();
-        if (eclState->hasDoubleGridProperty("MULTZ-"))
-            multzMinus = eclState->getDoubleGridProperty("MULTZ-")->getData();
-
-        // making this specific to clang or gcc > 4.7 is slightly hacky, but this call is
-        // only an optimization anyway...
-#if defined __clang__ || (__GNUC__ > 4 && __GNUC_MINOR__ >= 7)
-        // resize the hash map to a appropriate size for a conforming 3D grid
-        float maxLoadFactor = intersectionIntrinsicPermeability_.max_load_factor();
-
-        intersectionIntrinsicPermeability_.reserve(numElements * 6 / maxLoadFactor * 1.05 );
-#endif
-
-        auto elemIt = this->gridView().template begin</*codim=*/0>();
-        const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
-        for (; elemIt != elemEndIt; ++elemIt) {
-            auto intersectIt = elemIt->ileafbegin();
-            const auto &intersectEndIt = elemIt->ileafend();
-            for (; intersectIt != intersectEndIt; ++intersectIt) {
-                if (!intersectIt->neighbor())
-                    // skip boundary intersections...
-                    continue;
-
-                // calculate the "intersection index"
-#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
-                size_t interiorElemIdx = this->elementMapper().index(intersectIt->inside());
-                size_t exteriorElemIdx = this->elementMapper().index(intersectIt->outside());
-#else
-                size_t interiorElemIdx = this->elementMapper().map(intersectIt->inside());
-                size_t exteriorElemIdx = this->elementMapper().map(intersectIt->outside());
-#endif
-
-                size_t elem1Idx = std::min(interiorElemIdx, exteriorElemIdx);
-                size_t elem2Idx = std::max(interiorElemIdx, exteriorElemIdx);
-
-                size_t intersectIdx = elem1Idx*numElements + elem2Idx;
-
-                // do nothing if this intersection was already seen "from the other side"
-                if (intersectionIntrinsicPermeability_.count(intersectIdx) > 0)
-                    continue;
-
-                auto K1 = intrinsicPermeability_[interiorElemIdx];
-                auto K2 = intrinsicPermeability_[exteriorElemIdx];
-
-                int interiorElemCartIdx = grid.globalCell()[interiorElemIdx];
-                int exteriorElemCartIdx = grid.globalCell()[exteriorElemIdx];
-
-                // local index of the face of the interior element which contains the
-                // intersection
-                int insideFaceIdx = intersectIt->indexInInside();
-
-                // take the transmissibility multipliers into account (i.e., the
-                // MULT[XYZ]-? keywords)
-                if (insideFaceIdx == 1) { // right
-                    K1 *= multx[interiorElemCartIdx];
-                    K2 *= multxMinus[exteriorElemCartIdx];
-                }
-                else if (insideFaceIdx == 0) { // left
-                    K1 *= multxMinus[interiorElemCartIdx];
-                    K2 *= multx[exteriorElemCartIdx];
-                }
-
-                else if (insideFaceIdx == 3) { // back
-                    K1 *= multy[interiorElemCartIdx];
-                    K2 *= multyMinus[exteriorElemCartIdx];
-                }
-                else if (insideFaceIdx == 2) { // front
-                    K1 *= multyMinus[interiorElemCartIdx];
-                    K2 *= multy[exteriorElemCartIdx];
-                }
-
-                else if (insideFaceIdx == 5) { // top
-                    K1 *= multz[interiorElemCartIdx];
-                    K2 *= multzMinus[exteriorElemCartIdx];
-                }
-                else if (insideFaceIdx == 4) { // bottom
-                    K1 *= multzMinus[interiorElemCartIdx];
-                    K2 *= multz[exteriorElemCartIdx];
-                }
-
-                // element-wise harmonic average
-                auto &K = intersectionIntrinsicPermeability_[intersectIdx];
-                K = 0.0;
-                for (int i = 0; i < dimWorld; ++i)
-                    for (int j = 0; j < dimWorld; ++j)
-                        K[i][j] = Opm::utils::harmonicAverage(K1[i][j], K2[i][j]);
-            }
-        }
-    }
-
     std::vector<Scalar> porosity_;
     std::vector<DimMatrix> intrinsicPermeability_;
-
+    EclTransmissibility<TypeTag> transmissibilities_;
-    // the intrinsic permeabilities for interior faces. since grids may be
-    // non-conforming, and there does not seem to be a mapper for interfaces in DUNE,
-    // these transmissibilities are accessed via the (elementIndex1, elementIndex2) pairs
-    // of the interfaces where
-    //
-    // elementIndex1 = min(interiorElementIndex, exteriorElementIndex)
-    //
-    // and
-    //
-    // elementIndex2 = max(interiorElementIndex, exteriorElementIndex)
-    //
-    // To make this perform better, this is first mingled into a single index using
-    //
-    // intersectionIndex = elementIndex1*numElements + elementIndex2
-    //
-    // as the index for the hash map.
-    std::unordered_map<size_t, DimMatrix> intersectionIntrinsicPermeability_;
 
     std::vector<unsigned short> materialParamTableIdx_;
     std::vector<MaterialLawParams> materialParams_;

applications/ebos/ecltransmissibility.hh

@@ -0,0 +1,330 @@
+/*
+  Copyright (C) 2014 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::EclTransmissibility
+ */
+#ifndef EWOMS_ECL_TRANSMISSIBILITY_HH
+#define EWOMS_ECL_TRANSMISSIBILITY_HH
+
+#include "eclgridmanager.hh"
+
+#include <dune/common/version.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <array>
+#include <vector>
+#include <unordered_map>
+
+namespace Ewoms {
+/*!
+ * \ingroup EclBlackOilSimulator
+ *
+ * \brief This class calculates the transmissibilites for grid faces according to the
+ *        Eclipse Technical Description.
+ */
+template <class TypeTag>
+class EclTransmissibility
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GridView::Intersection Intersection;
+
+    // Grid and world dimension
+    enum { dim = GridView::dimension };
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+
+public:
+    EclTransmissibility(const Simulator& simulator)
+        : simulator_(simulator)
+    {}
+
+    /*!
+     * \brief Actually compute the transmissibilty over a face as a pre-compute step.
+     *
+     * This code actually uses the direction specific "centroids" of
+     * each element. These "centroids" are _not_ the identical
+     * barycenter of the element, but the middle of the centers of the
+     * faces of the logical Cartesian cells, i.e., the centers of the
+     * faces of the reference elements. We do things this way because
+     * the barycenter of the element can be located outside of the
+     * element for sufficiently "ugly" (i.e., thin and "non-flat")
+     * elements which in turn leads to quite wrong
+     * permeabilities. This approach is probably not always correct
+     * either but at least it seems to be much better.
+     */
+    void finishInit()
+    {
+        const auto& elementMapper = simulator_.model().elementMapper();
+        const auto& gridView = simulator_.gridView();
+        const auto& problem = simulator_.problem();
+
+        int numElements = elementMapper.size();
+
+        // 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(simulator_.model().numGridDof() == numElements);
+
+        // calculate the axis specific centroids of all elements
+        std::array<std::vector<DimVector>, dimWorld> axisCentroids;
+
+        for (int dimIdx = 0; dimIdx < dimWorld; ++dimIdx)
+            axisCentroids[dimIdx].resize(numElements);
+
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+            int elemIdx = elementMapper.index(elemIt);
+#else
+            int elemIdx = elementMapper.map(elemIt);
+#endif
+
+            // get the geometry of the current element
+            const auto& geom = elemIt->geometry();
+
+            // compute the axis specific "centroids" used for the
+            // transmissibilities
+            for (int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+                DimVector x0Local(0.5);
+                DimVector x1Local(0.5);
+
+                x0Local[dimIdx] = 0.0;
+                x1Local[dimIdx] = 1.0;
+
+                DimVector x = geom.global(x0Local);
+                x += geom.global(x1Local);
+                x /= 2;
+
+                axisCentroids[dimIdx][elemIdx] = x;
+            }
+        }
+
+        Opm::EclipseStateConstPtr eclState = simulator_.gridManager().eclState();
+        const std::vector<double>& multx =
+            eclState->getDoubleGridProperty("MULTX")->getData();
+        const std::vector<double>& multy =
+            eclState->getDoubleGridProperty("MULTY")->getData();
+        const std::vector<double>& multz =
+            eclState->getDoubleGridProperty("MULTZ")->getData();
+        const std::vector<double>& multxMinus =
+            eclState->getDoubleGridProperty("MULTX-")->getData();
+        const std::vector<double>& multyMinus =
+            eclState->getDoubleGridProperty("MULTY-")->getData();
+        const std::vector<double>& multzMinus =
+            eclState->getDoubleGridProperty("MULTZ-")->getData();
+
+        const std::vector<double>& ntg =
+            eclState->getDoubleGridProperty("NTG")->getData();
+
+        // reserving some space in the hashmap upfront saves quite a
+        // bit of time because resizes are costly for hashmaps and
+        // there would be quite a few of them if we would not have a
+        // rough idea of how large the final map will be (the rough
+        // idea is a conforming Cartesian grid)
+        trans_.reserve(numElements*3*1.05);
+
+        // compute the transmissibilities for all intersections
+        elemIt = gridView.template begin</*codim=*/ 0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            auto isIt = elemIt->ileafbegin();
+            const auto& isEndIt = elemIt->ileafend();
+            for (; isIt != isEndIt; ++ isIt) {
+                // ignore boundary intersections for now (TODO?)
+                if (isIt->boundary())
+                    continue;
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+                int insideElemIdx = elementMapper.index(*isIt->inside());
+                int outsideElemIdx = elementMapper.index(*isIt->outside());
+#else
+                int insideElemIdx = elementMapper.map(*isIt->inside());
+                int outsideElemIdx = elementMapper.map(*isIt->outside());
+#endif
+
+                // we only need to calculate a face's transmissibility
+                // once...
+                if (insideElemIdx > outsideElemIdx)
+                    continue;
+
+                // local indices of the faces of the inside and
+                // outside elements which contain the intersection
+                int insideFaceIdx = isIt->indexInInside();
+                int outsideFaceIdx = isIt->indexInOutside();
+
+                Scalar halfTrans1;
+                Scalar halfTrans2;
+
+                computeHalfTrans_(halfTrans1,
+                                  *isIt,
+                                  insideFaceIdx,
+                                  distanceVector_(*isIt,
+                                                  isIt->indexInInside(),
+                                                  insideElemIdx,
+                                                  axisCentroids),
+                                  problem.intrinsicPermeability(insideElemIdx));
+                computeHalfTrans_(halfTrans2,
+                                  *isIt,
+                                  outsideFaceIdx,
+                                  distanceVector_(*isIt,
+                                                  isIt->indexInOutside(),
+                                                  outsideElemIdx,
+                                                  axisCentroids),
+                                  problem.intrinsicPermeability(outsideElemIdx));
+
+                applyNtg_(halfTrans1, insideFaceIdx, insideElemIdx, ntg);
+                applyNtg_(halfTrans2, outsideFaceIdx, outsideElemIdx, ntg);
+
+                // convert half transmissibilities to full face
+                // transmissibilities using the harmonic mean
+                Scalar trans = 1.0 / (1.0/halfTrans1 + 1.0/halfTrans2);
+
+                // apply the full face transmissibility multipliers
+                // for the inside ...
+                applyMultipliers_(trans, insideFaceIdx, insideElemIdx,
+                                  multx, multxMinus,
+                                  multy, multyMinus,
+                                  multz, multzMinus);
+                // ... and outside elements
+                applyMultipliers_(trans, outsideFaceIdx, outsideElemIdx,
+                                  multx, multxMinus,
+                                  multy, multyMinus,
+                                  multz, multzMinus);
+
+                trans_[isId_(insideElemIdx, outsideElemIdx)] = trans;
+            }
+        }
+    }
+
+    Scalar transmissibility(int elemIdx1, int elemIdx2) const
+    { return trans_.at(isId_(elemIdx1, elemIdx2)); }
+
+private:
+    std::uint64_t isId_(int elemIdx1, int elemIdx2) const
+    {
+        static const int elemIdxShift = 32; // bits
+
+        int elemAIdx = std::min(elemIdx1, elemIdx2);
+        std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
+
+        return (elemBIdx<<elemIdxShift) + elemAIdx;
+    }
+
+    void computeHalfTrans_(Scalar& halfTrans,
+                           const Intersection& is,
+                           int faceIdx, // in the reference element that contains the intersection
+                           const DimVector& distance,
+                           const DimMatrix& perm) const
+    {
+        int dimIdx = faceIdx/2;
+        assert(dimIdx < dimWorld);
+        halfTrans = perm[dimIdx][dimIdx];
+        halfTrans *= is.geometry().volume();
+        halfTrans *= std::abs<Scalar>(is.centerUnitOuterNormal()*distance);
+        halfTrans /= distance*distance;
+    }
+
+    DimVector distanceVector_(const Intersection& is,
+                              int faceIdx, // in the reference element that contains the intersection
+                              int elemIdx,
+                              const std::array<std::vector<DimVector>, dimWorld>& axisCentroids) const
+    {
+        int dimIdx = faceIdx/2;
+        assert(dimIdx < dimWorld);
+        DimVector x = is.geometry().center();
+        x -= axisCentroids[dimIdx][elemIdx];
+
+        return x;
+    }
+
+    void applyMultipliers_(Scalar &trans, int faceIdx, int elemIdx,
+                           const std::vector<Scalar>& multx,
+                           const std::vector<Scalar>& multxMinus,
+                           const std::vector<Scalar>& multy,
+                           const std::vector<Scalar>& multyMinus,
+                           const std::vector<Scalar>& multz,
+                           const std::vector<Scalar>& multzMinus) const
+    {
+        // apply multiplyer for the transmissibility of the face. (the
+        // face index is the index of the reference-element face which
+        // contains the intersection of interest.)
+        switch (faceIdx) {
+        case 0: // left
+            trans *= multxMinus[elemIdx];
+            break;
+        case 1: // right
+            trans *= multx[elemIdx];
+            break;
+
+        case 2: // front
+            trans *= multyMinus[elemIdx];
+            break;
+        case 3: // back
+            trans *= multy[elemIdx];
+            break;
+
+        case 4: // bottom
+            trans *= multzMinus[elemIdx];
+            break;
+        case 5: // top
+            trans *= multz[elemIdx];
+            break;
+        }
+    }
+
+    void applyNtg_(Scalar &trans, int faceIdx, int elemIdx,
+                   const std::vector<Scalar>& ntg) const
+    {
+        // apply multiplyer for the transmissibility of the face. (the
+        // face index is the index of the reference-element face which
+        // contains the intersection of interest.)
+        switch (faceIdx) {
+        case 0: // left
+            trans *= ntg[elemIdx];
+            break;
+        case 1: // right
+            trans *= ntg[elemIdx];
+            break;
+
+        case 2: // front
+            trans *= ntg[elemIdx];
+            break;
+        case 3: // back
+            trans *= ntg[elemIdx];
+            break;
+
+            // NTG does not apply to top and bottom faces
+        }
+    }
+
+    const Simulator& simulator_;
+    std::unordered_map<std::uint64_t, Scalar> trans_;
+};
+
+} // namespace Ewoms
+
+#endif