Merge pull request #78 from dr-robertk/PR/use-correct-cellcenterdepth

Use correct cellcenterdepth

applications/ebos/eclfluxmodule.hh

@@ -210,7 +210,6 @@ protected:
         Valgrind::SetUndefined(*this);
 
         const auto& problem = elemCtx.problem();
-        const auto& grid = elemCtx.simulator().gridManager().grid();
         const auto& stencil = elemCtx.stencil(timeIdx);
         const auto& scvf = stencil.interiorFace(scvfIdx);
 
@@ -237,8 +236,8 @@ protected:
         // grids). The "good" solution would be to take the Z coordinate of the element
         // centers, but since ECL seems to like to be inconsistent on that front, it
         // needs to be done like here...
-        Scalar zIn = grid.cellCenterDepth(I);
-        Scalar zEx = grid.cellCenterDepth(J);
+        Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx_, timeIdx);
+        Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx_, timeIdx);
 
         // the distances from the DOF's depths. (i.e., the additional depth of the
         // exterior DOF)

applications/ebos/eclproblem.hh

@@ -319,6 +319,7 @@ public:
             this->gravity_[dim - 1] = 9.80665;
 
         initFluidSystem_();
+        updateElementDepths_();
         readRockParameters_();
         readMaterialParameters_();
         transmissibilities_.finishInit();
@@ -613,6 +614,19 @@ public:
         return porosity_[globalSpaceIdx];
     }
 
+    /*!
+     * \brief Returns the depth of an degree of freedom [m]
+     *
+     * For ECL problems this is defined as the average of the depth of an element and is
+     * thus slightly different from the depth of an element's centroid.
+     */
+    template <class Context>
+    Scalar dofCenterDepth(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
+    {
+        unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        return elementCenterDepth_[globalSpaceIdx];
+    }
+
     /*!
      * \copydoc BlackoilProblem::rockCompressibility
      */
@@ -791,6 +805,45 @@ private:
     static bool enableEclOutput_()
     { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
 
+    Scalar cellCenterDepth( const Element& element ) const
+    {
+        typedef typename Element :: Geometry Geometry;
+        static constexpr int zCoord = Geometry ::dimension - 1;
+        Scalar zz = 0.0;
+
+        const Geometry geometry = element.geometry();
+
+        const int corners = geometry.corners();
+        for (int i=0; i<corners; ++i)
+        {
+            zz += geometry.corner( i )[ zCoord ];
+        }
+        return zz/Scalar(corners);
+    }
+
+    void updateElementDepths_()
+    {
+        const auto& gridManager = this->simulator().gridManager();
+        const auto& gridView = gridManager.gridView();
+        const auto& elemMapper = this->elementMapper();;
+
+        int numElements = gridView.size(/*codim=*/0);
+        elementCenterDepth_.resize(numElements);
+
+        auto elemIt = gridView.template begin</*codim=*/0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const Element& element = *elemIt;
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+            const unsigned int elemIdx = elemMapper.index(element);
+#else
+            const unsigned int elemIdx = elemMapper.map(element);
+#endif
+
+            elementCenterDepth_[elemIdx] = cellCenterDepth( element );
+        }
+    }
+
     void readRockParameters_()
     {
         auto deck = this->simulator().gridManager().deck();
@@ -1186,6 +1239,7 @@ private:
     }
 
     std::vector<Scalar> porosity_;
+    std::vector<Scalar> elementCenterDepth_;
     std::vector<DimMatrix> intrinsicPermeability_;
     EclTransmissibility<TypeTag> transmissibilities_;