ebos: use the transmissibilities as edge weights for load balancing

this makes creating the grid a bit slower because the
transmissibilities need to be calculated twice: once for the
sequential grid and once for the distributed one. while corresponds to
the way `flow_legacy` does the load balancing and it should allow
better results, this does not seem to be the case for the Norne deck
if ZOLTAN is not available:

After loadbalancing process 3 has 4413 cells.
After loadbalancing process 2 has 12390 cells.
After loadbalancing process 0 has 13629 cells.
After loadbalancing process 1 has 21253 cells.

i.e., process 1 is responsible for almost 5 as many cells as process
3.

ebos/eclcpgridmanager.hh

@@ -28,9 +28,12 @@
 #define EWOMS_ECL_CP_GRID_MANAGER_HH
 
 #include "eclbasegridmanager.hh"
+#include "ecltransmissibility.hh"
 
 #include <dune/grid/CpGrid.hpp>
 
+#include <dune/common/version.hh>
+
 namespace Ewoms {
 template <class TypeTag>
 class EclCpGridManager;
@@ -59,6 +62,7 @@ class EclCpGridManager : public EclBaseGridManager<TypeTag>
 
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
 
 public:
     typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
@@ -135,9 +139,55 @@ public:
         MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
         MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
 
-        // distribute the grid and switch to the distributed view.
+        if (mpiSize > 1) {
-        grid_->loadBalance();
+            // the CpGrid's loadBalance() method likes to have the transmissibilities as
-        grid_->switchToDistributedView();
+            // its edge weights. since this is (kind of) a layering violation and
+            // transmissibilities are relatively expensive to compute, we only do it if
+            // more than a single process is involved in the simulation.
+            cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);
+            EclTransmissibility<TypeTag> eclTrans(*this);
+            eclTrans.update();
+
+            // convert to transmissibility for faces
+            // TODO: grid_->numFaces() is not generic. use grid_->size(1) instead? (might
+            // not work)
+            const auto& gridView = grid_->leafGridView();
+            unsigned numFaces = grid_->numFaces();
+            std::vector<double> faceTrans(numFaces, 0.0);
+            ElementMapper elemMapper(this->gridView());
+            auto elemIt = gridView.template begin</*codim=*/0>();
+            const auto& elemEndIt = gridView.template end</*codim=*/0>();
+            for (; elemIt != elemEndIt; ++ elemIt) {
+                const auto& elem = *elemIt;
+                auto isIt = gridView.ibegin(elem);
+                const auto& isEndIt = gridView.iend(elem);
+                for (; isIt != isEndIt; ++ isIt) {
+                    const auto& is = *isIt;
+                    if (!is.neighbor())
+                        continue;
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
+                    unsigned I = elemMapper.index(is.inside());
+                    unsigned J = elemMapper.index(is.outside());
+#else
+                    unsigned I = elemMapper.map(is.inside());
+                    unsigned J = elemMapper.map(is.outside());
+#endif
+
+                    // FIXME (?): this is not portable!
+                    unsigned faceIdx = is.id();
+
+                    faceTrans[faceIdx] = eclTrans.transmissibility(I, J);
+                }
+            }
+
+            //distribute the grid and switch to the distributed view.
+            grid_->loadBalance(&this->eclState(), faceTrans.data());
+            grid_->switchToDistributedView();
+
+            delete cartesianIndexMapper_;
+            cartesianIndexMapper_ = nullptr;
+        }
 #endif
 
         cartesianIndexMapper_ = new CartesianIndexMapper(*grid_);

ebos/eclproblem.hh

@@ -301,7 +301,7 @@ public:
      */
     EclProblem(Simulator& simulator)
         : ParentType(simulator)
-        , transmissibilities_(simulator)
+        , transmissibilities_(simulator.gridManager())
         , thresholdPressures_(simulator)
         , wellManager_(simulator)
         , deckUnits_(simulator)
@@ -597,7 +597,7 @@ public:
                                            unsigned timeIdx) const
     {
         unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
-        return intrinsicPermeability_[globalSpaceIdx];
+        return transmissibilities_.permeability(globalSpaceIdx);
     }
 
     /*!
@@ -607,7 +607,7 @@ public:
      * Its main (only?) usage is the ECL transmissibility calculation code...
      */
     const DimMatrix& intrinsicPermeability(unsigned globalElemIdx) const
-    { return intrinsicPermeability_[globalElemIdx]; }
+    { return transmissibilities_.permeability(globalElemIdx); }
 
     /*!
      * \copydoc BlackOilBaseProblem::transmissibility
@@ -930,48 +930,10 @@ private:
         const auto& gridManager = this->simulator().gridManager();
         const auto& deck = gridManager.deck();
         const auto& eclState = gridManager.eclState();
-        const auto& props = eclState.get3DProperties();
-
-        size_t numDof = this->model().numGridDof();
 
         // the PVT region number
         updatePvtnum_();
 
-        ////////////////////////////////
-        // intrinsic permeability
-
-        intrinsicPermeability_.resize(numDof);
-
-        // read the intrinsic permeabilities from the eclState. Note that all arrays
-        // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
-        // opm-grid CpGrid object might remove a few elements...
-        if (props.hasDeckDoubleGridProperty("PERMX")) {
-                const std::vector<double>& permxData =
-                props.getDoubleGridProperty("PERMX").getData();
-            std::vector<double> permyData(permxData);
-            if (props.hasDeckDoubleGridProperty("PERMY"))
-                permyData = props.getDoubleGridProperty("PERMY").getData();
-            std::vector<double> permzData(permxData);
-            if (props.hasDeckDoubleGridProperty("PERMZ"))
-                permzData = props.getDoubleGridProperty("PERMZ").getData();
-
-            for (size_t dofIdx = 0; dofIdx < numDof; ++ dofIdx) {
-                unsigned cartesianElemIdx = gridManager.cartesianIndex(dofIdx);
-                intrinsicPermeability_[dofIdx] = 0.0;
-                intrinsicPermeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
-                intrinsicPermeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
-                intrinsicPermeability_[dofIdx][2][2] = permzData[cartesianElemIdx];
-            }
-
-            // for now we don't care about non-diagonal entries
-        }
-        else
-            OPM_THROW(std::logic_error,
-                      "Can't read the intrinsic permeability from the ecl state. "
-                      "(The PERM{X,Y,Z} keywords are missing)");
-        ////////////////////////////////
-
-
         ////////////////////////////////
         // porosity
         updatePorosity_();
@@ -979,6 +941,7 @@ private:
 
         ////////////////////////////////
         // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
+        size_t numDof = this->model().numGridDof();
         std::vector<int> compressedToCartesianElemIdx(numDof);
         for (unsigned elemIdx = 0; elemIdx < numDof; ++elemIdx)
             compressedToCartesianElemIdx[elemIdx] = gridManager.cartesianIndex(elemIdx);
@@ -1350,7 +1313,6 @@ private:
 
     std::vector<Scalar> porosity_;
     std::vector<Scalar> elementCenterDepth_;
-    std::vector<DimMatrix> intrinsicPermeability_;
     EclTransmissibility<TypeTag> transmissibilities_;
 
     std::shared_ptr<EclMaterialLawManager> materialLawManager_;

ebos/ecltransmissibility.hh

@@ -49,7 +49,8 @@ namespace Ewoms {
 namespace Properties {
 NEW_PROP_TAG(GridView);
 NEW_PROP_TAG(Scalar);
-NEW_PROP_TAG(Simulator);
+NEW_PROP_TAG(GridManager);
+NEW_PROP_TAG(ElementMapper);
 }
 
 /*!
@@ -63,7 +64,8 @@ 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 GET_PROP_TYPE(TypeTag, GridManager) GridManager;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
     typedef typename GridView::Intersection Intersection;
 
     // Grid and world dimension
@@ -73,8 +75,8 @@ class EclTransmissibility
     typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
 
 public:
-    EclTransmissibility(const Simulator& simulator)
+    EclTransmissibility(const GridManager& gridManager)
-        : simulator_(simulator)
+        : gridManager_(gridManager)
     {}
 
     /*!
@@ -96,25 +98,19 @@ public:
 
     void update()
     {
-        const auto& problem = simulator_.problem();
+        const auto& gridView = gridManager_.gridView();
-        const auto& gridManager = simulator_.gridManager();
+        const auto& cartMapper = gridManager_.cartesianIndexMapper();
-        const auto& gridView = simulator_.gridView();
+        const auto& eclState = gridManager_.eclState();
-        const auto& elementMapper = simulator_.model().elementMapper();
-        const auto& cartMapper = gridManager.cartesianIndexMapper();
-        const auto& eclState = gridManager.eclState();
         const auto& eclGrid = eclState.getInputGrid();
-        const auto& transMult = eclState.getTransMult();
+        auto& transMult = eclState.getTransMult();
+        ElementMapper elemMapper(gridView);
 
         const std::vector<double>& ntg =
             eclState.get3DProperties().getDoubleGridProperty("NTG").getData();
 
-        unsigned numElements = elementMapper.size();
+        unsigned numElements = elemMapper.size();
 
-        // this code assumes that the DOFs are the elements. (i.e., an
+        extractPermeability_(eclState);
-        // 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;
@@ -127,9 +123,9 @@ public:
         for (; elemIt != elemEndIt; ++elemIt) {
             const auto& elem = *elemIt;
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
-            unsigned elemIdx = elementMapper.index(elem);
+            unsigned elemIdx = elemMapper.index(elem);
 #else
-            unsigned elemIdx = elementMapper.map(elem);
+            unsigned elemIdx = elemMapper.map(elem);
 #endif
 
             // compute the axis specific "centroids" used for the transmissibilities. for
@@ -168,11 +164,11 @@ public:
                 const auto& inside = intersection.inside();
                 const auto& outside = intersection.outside();
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2,4)
-                unsigned insideElemIdx = elementMapper.index(inside);
+                unsigned insideElemIdx = elemMapper.index(inside);
-                unsigned outsideElemIdx = elementMapper.index(outside);
+                unsigned outsideElemIdx = elemMapper.index(outside);
 #else
-                unsigned insideElemIdx = elementMapper.map(*inside);
+                unsigned insideElemIdx = elemMapper.map(*inside);
-                unsigned outsideElemIdx = elementMapper.map(*outside);
+                unsigned outsideElemIdx = elemMapper.map(*outside);
 #endif
 
                 // we only need to calculate a face's transmissibility
@@ -198,7 +194,7 @@ public:
                                                   intersection.indexInInside(),
                                                   insideElemIdx,
                                                   axisCentroids),
-                                  problem.intrinsicPermeability(insideElemIdx));
+                                  permeability_[insideElemIdx]);
                 computeHalfTrans_(halfTrans2,
                                   intersection,
                                   outsideFaceIdx,
@@ -206,7 +202,7 @@ public:
                                                   intersection.indexInOutside(),
                                                   outsideElemIdx,
                                                   axisCentroids),
-                                  problem.intrinsicPermeability(outsideElemIdx));
+                                  permeability_[outsideElemIdx]);
 
                 applyNtg_(halfTrans1, insideFaceIdx, insideCartElemIdx, ntg);
                 applyNtg_(halfTrans2, outsideFaceIdx, outsideCartElemIdx, ntg);
@@ -257,10 +253,50 @@ public:
         }
     }
 
+    const DimMatrix& permeability(unsigned elemIdx) const
+    { return permeability_[elemIdx]; }
+
     Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
     { return trans_.at(isId_(elemIdx1, elemIdx2)); }
 
 private:
+    void extractPermeability_(const Opm::EclipseState& eclState)
+    {
+        const auto& props = gridManager_.eclState().get3DProperties();
+
+        unsigned numElem = gridManager_.gridView().size(/*codim=*/0);
+        permeability_.resize(numElem);
+
+        // read the intrinsic permeabilities from the eclState. Note that all arrays
+        // provided by eclState are one-per-cell of "uncompressed" grid, whereas the
+        // simulation grid might remove a few elements. (e.g. because it is distributed
+        // over several processes.)
+        if (props.hasDeckDoubleGridProperty("PERMX")) {
+            const std::vector<double>& permxData =
+                props.getDoubleGridProperty("PERMX").getData();
+            std::vector<double> permyData(permxData);
+            if (props.hasDeckDoubleGridProperty("PERMY"))
+                permyData = props.getDoubleGridProperty("PERMY").getData();
+            std::vector<double> permzData(permxData);
+            if (props.hasDeckDoubleGridProperty("PERMZ"))
+                permzData = props.getDoubleGridProperty("PERMZ").getData();
+
+            for (size_t dofIdx = 0; dofIdx < numElem; ++ dofIdx) {
+                unsigned cartesianElemIdx = gridManager_.cartesianIndex(dofIdx);
+                permeability_[dofIdx] = 0.0;
+                permeability_[dofIdx][0][0] = permxData[cartesianElemIdx];
+                permeability_[dofIdx][1][1] = permyData[cartesianElemIdx];
+                permeability_[dofIdx][2][2] = permzData[cartesianElemIdx];
+            }
+
+            // for now we don't care about non-diagonal entries
+        }
+        else
+            OPM_THROW(std::logic_error,
+                      "Can't read the intrinsic permeability from the ecl state. "
+                      "(The PERM{X,Y,Z} keywords are missing)");
+    }
+
     std::uint64_t isId_(unsigned elemIdx1, unsigned elemIdx2) const
     {
         static const unsigned elemIdxShift = 32; // bits
@@ -362,7 +398,8 @@ private:
         }
     }
 
-    const Simulator& simulator_;
+    const GridManager& gridManager_;
+    std::vector<DimMatrix> permeability_;
     std::unordered_map<std::uint64_t, Scalar> trans_;
 };