From 457d270bdf6970dfe4e520fa4c72c38af756fb48 Mon Sep 17 00:00:00 2001 From: Markus Blatt Date: Thu, 30 Sep 2021 11:50:22 +0200 Subject: [PATCH] [refactor] Move cell centroid lookup magic to vanguard for reuse. In addition to transmissibilities We will need this for the tracers, too, and should prevent code duplication. --- ebos/eclalugridvanguard.hh | 18 ++++++++++--- ebos/eclbasevanguard.hh | 42 +++++++++++++++++++++++-------- ebos/eclcpgridvanguard.hh | 14 ++++++++++- ebos/eclgenericcpgridvanguard.cc | 14 +++-------- ebos/eclpolyhedralgridvanguard.hh | 17 ++++++++++--- ebos/ecltransmissibility.cc | 12 ++------- ebos/ecltransmissibility.hh | 5 ++-- 7 files changed, 83 insertions(+), 39 deletions(-) diff --git a/ebos/eclalugridvanguard.hh b/ebos/eclalugridvanguard.hh index 666b2260c..b847db073 100644 --- a/ebos/eclalugridvanguard.hh +++ b/ebos/eclalugridvanguard.hh @@ -102,7 +102,6 @@ public: ~EclAluGridVanguard() { - delete cartesianIndexMapper_; delete equilCartesianIndexMapper_; delete grid_; delete equilGrid_; @@ -178,6 +177,18 @@ public: const EquilCartesianIndexMapper& equilCartesianIndexMapper() const { return *equilCartesianIndexMapper_; } + /*! + * \brief Get function to query cell centroids for a distributed grid. + * + * Currently this only non-empty for a loadbalanced CpGrid. + * It is a function return the centroid for the given element + * index. + */ + std::function(int)> + cellCentroids() const + { + return this->cellCentroids_(cartesianIndexMapper_.get()); + } protected: void createGrids_() { @@ -210,7 +221,8 @@ protected: grid_ = factory.convert(*equilGrid_, cartesianCellId_); cartesianIndexMapper_ = - new CartesianIndexMapper(*grid_, cartesianDimension_, cartesianCellId_); + std::make_unique(*grid_, cartesianDimension_, + cartesianCellId_); } void filterConnections_() @@ -222,7 +234,7 @@ protected: EquilGrid* equilGrid_; std::vector cartesianCellId_; std::array cartesianDimension_; - CartesianIndexMapper* cartesianIndexMapper_; + std::unique_ptr cartesianIndexMapper_; EquilCartesianIndexMapper* equilCartesianIndexMapper_; }; diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh index ec6ab0d4b..fe2611446 100644 --- a/ebos/eclbasevanguard.hh +++ b/ebos/eclbasevanguard.hh @@ -36,6 +36,7 @@ #include #include #include +#include #include #include @@ -185,6 +186,7 @@ public: protected: static const int dimension = Grid::dimension; + static const int dimensionworld = Grid::dimensionworld; using Element = typename GridView::template Codim<0>::Entity; using CartesianIndexMapper = Dune::CartesianIndexMapper; @@ -329,16 +331,6 @@ public: { return asImp_().equilCartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); } - /*! - * \brief Get the cell centroids for a distributed grid. - * - * Currently this only non-empty for a loadbalanced CpGrid. - */ - const std::vector& cellCentroids() const - { - return centroids_; - } - /*! * \brief Returns the depth of a degree of freedom [m] * @@ -389,6 +381,36 @@ public: } protected: + /*! + * \brief Get function to query cell centroids for a distributed grid. + * + * Currently this only non-empty for a loadbalanced CpGrid. + * It is a function return the centroid for the given element + * index. + * \param cartMapper The cartesian index mapper for lookup of + * cartesian indices + */ + template + std::function(int)> + cellCentroids_(const CartMapper& cartMapper) const + { + return [this, cartMapper](int elemIdx) { + const auto& centroids = this->centroids_; + auto rank = this->gridView().comm().rank(); + std::array centroid; + if (rank == 0) { + unsigned cartesianCellIdx = cartMapper.cartesianIndex(elemIdx); + centroid = this->eclState().getInputGrid().getCellCenter(cartesianCellIdx); + } else + { + std::copy(centroids.begin() + elemIdx * dimensionworld, + centroids.begin() + (elemIdx + 1) * dimensionworld, + centroid.begin()); + } + return centroid; + }; + } + void callImplementationInit() { asImp_().createGrids_(); diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh index 508308c97..c35fed0bd 100644 --- a/ebos/eclcpgridvanguard.hh +++ b/ebos/eclcpgridvanguard.hh @@ -90,6 +90,7 @@ public: using EquilGrid = GetPropType; using GridView = GetPropType; using TransmissibilityType = EclTransmissibility; + static const int dimensionworld = Grid::dimensionworld; private: typedef Dune::CartesianIndexMapper CartesianIndexMapper; @@ -138,7 +139,6 @@ public: this->eclState(), this->parallelWells_); #endif - this->allocCartMapper(); this->updateGridView_(); this->updateCartesianToCompressedMapping_(); this->updateCellDepths_(); @@ -149,6 +149,18 @@ public: #endif } + /*! + * \brief Get function to query cell centroids for a distributed grid. + * + * Currently this only non-empty for a loadbalanced CpGrid. + * It is a function return the centroid for the given element + * index. + */ + std::function(int)> + cellCentroids() const + { + return this->cellCentroids_(this->cartesianIndexMapper()); + } protected: void createGrids_() diff --git a/ebos/eclgenericcpgridvanguard.cc b/ebos/eclgenericcpgridvanguard.cc index 61c3db5ba..640feaefd 100644 --- a/ebos/eclgenericcpgridvanguard.cc +++ b/ebos/eclgenericcpgridvanguard.cc @@ -92,7 +92,6 @@ void EclGenericCpGridVanguard::doLoadBalance_(Dun // 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_.reset(new CartesianIndexMapper(*grid_)); if (grid_->size(0)) { this->allocTrans(); @@ -175,8 +174,6 @@ void EclGenericCpGridVanguard::doLoadBalance_(Dun } grid_->switchToDistributedView(); - cartesianIndexMapper_.reset(); - // Calling Schedule::filterConnections would remove any perforated // cells that exist only on other ranks even in the case of distributed wells // But we need all connections to figure out the first cell of a well (e.g. for @@ -210,12 +207,6 @@ void EclGenericCpGridVanguard::distributeFieldPro } #endif -template -void EclGenericCpGridVanguard::allocCartMapper() -{ - this->cartesianIndexMapper_.reset(new CartesianIndexMapper(this->grid())); -} - template void EclGenericCpGridVanguard::doCreateGrids_(EclipseState& eclState) { @@ -258,6 +249,9 @@ void EclGenericCpGridVanguard::doCreateGrids_(Ecl } } + + cartesianIndexMapper_ = std::make_unique(*grid_); + #if HAVE_MPI { const bool has_numerical_aquifer = eclState.aquifer().hasNumericalAquifer(); @@ -285,7 +279,7 @@ void EclGenericCpGridVanguard::doCreateGrids_(Ecl if (mpiRank == 0) { equilGrid_.reset(new Dune::CpGrid(*grid_)); - equilCartesianIndexMapper_.reset(new CartesianIndexMapper(*equilGrid_)); + equilCartesianIndexMapper_ = std::make_unique(*equilGrid_); std::vector actnum = UgGridHelpers::createACTNUM(*grid_); auto &field_props = eclState.fieldProps(); diff --git a/ebos/eclpolyhedralgridvanguard.hh b/ebos/eclpolyhedralgridvanguard.hh index b02757b86..5ec066573 100644 --- a/ebos/eclpolyhedralgridvanguard.hh +++ b/ebos/eclpolyhedralgridvanguard.hh @@ -92,7 +92,7 @@ private: typedef Grid* GridPointer; typedef EquilGrid* EquilGridPointer; typedef Dune::CartesianIndexMapper CartesianIndexMapper; - typedef CartesianIndexMapper* CartesianIndexMapperPointer; + typedef std::unique_ptr CartesianIndexMapperPointer; public: EclPolyhedralGridVanguard(Simulator& simulator) @@ -104,7 +104,6 @@ public: ~EclPolyhedralGridVanguard() { - delete cartesianIndexMapper_; delete grid_; } @@ -184,11 +183,23 @@ public: return simulator_.problem().eclTransmissibilities(); } + /*! + * \brief Get function to query cell centroids for a distributed grid. + * + * Currently this only non-empty for a loadbalanced CpGrid. + * It is a function return the centroid for the given element + * index. + */ + std::function(int)> + cellCentroids() const + { + return this->cellCentroids_(this->cartesianIndexMapper()); + } protected: void createGrids_() { grid_ = new Grid(this->eclState().getInputGrid(), this->eclState().fieldProps().porv(true)); - cartesianIndexMapper_ = new CartesianIndexMapper(*grid_); + cartesianIndexMapper_ = std::make_unique(*grid_); this->updateCartesianToCompressedMapping_(); this->updateCellDepths_(); } diff --git a/ebos/ecltransmissibility.cc b/ebos/ecltransmissibility.cc index 5ee307ba8..28cd1ba66 100644 --- a/ebos/ecltransmissibility.cc +++ b/ebos/ecltransmissibility.cc @@ -88,7 +88,7 @@ EclTransmissibility(const EclipseState& eclState, const GridView& gridView, const Dune::CartesianIndexMapper& cartMapper, const Grid& grid, - const std::vector& centroids, + std::function(int)> centroids, bool enableEnergy, bool enableDiffusivity) : eclState_(eclState) @@ -174,15 +174,7 @@ update(bool global) // compute the axis specific "centroids" used for the transmissibilities. for // consistency with the flow simulator, we use the element centers as // computed by opm-parser's Opm::EclipseGrid class for all axes. - std::array centroid; - if (gridView_.comm().rank() == 0) { - const auto& eclGrid = eclState_.getInputGrid(); - unsigned cartesianCellIdx = cartMapper_.cartesianIndex(elemIdx); - centroid = eclGrid.getCellCenter(cartesianCellIdx); - } else - std::copy(centroids_.begin() + centroidIdx * dimWorld, - centroids_.begin() + (centroidIdx + 1) * dimWorld, - centroid.begin()); + std::array centroid = centroids_(elemIdx); for (unsigned axisIdx = 0; axisIdx < dimWorld; ++axisIdx) for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh index bf1c51a20..7b6ffc033 100644 --- a/ebos/ecltransmissibility.hh +++ b/ebos/ecltransmissibility.hh @@ -38,6 +38,7 @@ #include #include #include +#include namespace Opm { @@ -58,7 +59,7 @@ public: const GridView& gridView, const Dune::CartesianIndexMapper& cartMapper, const Grid& grid, - const std::vector& centroids, + std::function(int)> centroids, bool enableEnergy, bool enableDiffusivity); @@ -236,7 +237,7 @@ protected: const GridView& gridView_; const Dune::CartesianIndexMapper& cartMapper_; const Grid& grid_; - const std::vector& centroids_; + std::function(int)> centroids_; Scalar transmissibilityThreshold_; std::map, Scalar> transBoundary_; std::map, Scalar> thermalHalfTransBoundary_;