From 2fa90d24f6db187375be19ab858acd5905442f93 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 11 Dec 2019 22:04:40 +0100 Subject: [PATCH] [bugfix][ebos] Make compile with PolyhedralGrid again. --- ebos/eclpolyhedralgridvanguard.hh | 26 +++++++++- ebos/eclwriter.hh | 51 ++++++++++--------- .../ExtractParallelGridInformationToISTL.hpp | 5 ++ opm/simulators/linalg/ISTLSolverEbos.hpp | 16 +++--- .../linalg/findOverlapRowsAndColumns.hpp | 30 ++++++----- .../wells/BlackoilWellModel_impl.hpp | 9 ++-- 6 files changed, 86 insertions(+), 51 deletions(-) diff --git a/ebos/eclpolyhedralgridvanguard.hh b/ebos/eclpolyhedralgridvanguard.hh index 4ca53dfff..8a187637a 100644 --- a/ebos/eclpolyhedralgridvanguard.hh +++ b/ebos/eclpolyhedralgridvanguard.hh @@ -28,6 +28,7 @@ #define EWOMS_ECL_POLYHEDRAL_GRID_VANGUARD_HH #include "eclbasevanguard.hh" +#include "ecltransmissibility.hh" #include @@ -78,7 +79,8 @@ private: public: EclPolyhedralGridVanguard(Simulator& simulator) - : EclBaseVanguard(simulator) + : EclBaseVanguard(simulator), + simulator_( simulator ) { this->callImplementationInit(); } @@ -147,6 +149,24 @@ public: const CartesianIndexMapper& equilCartesianIndexMapper() const { return *cartesianIndexMapper_; } + + /*! + * \brief Free the memory occupied by the global transmissibility object. + * + * After writing the initial solution, this array should not be necessary anymore. + */ + void releaseGlobalTransmissibilities() + { + } + + std::unordered_set defunctWellNames() const + { return defunctWellNames_; } + + const EclTransmissibility& globalTransmissibility() const + { + return simulator_.problem().eclTransmissibilities(); + } + protected: void createGrids_() { @@ -162,8 +182,12 @@ protected: // not handling the removal of completions for this type of grid yet. } + Simulator& simulator_; + GridPointer grid_; CartesianIndexMapperPointer cartesianIndexMapper_; + + std::unordered_set defunctWellNames_; }; } // namespace Opm diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh index fdd71aaa8..513e6b368 100644 --- a/ebos/eclwriter.hh +++ b/ebos/eclwriter.hh @@ -148,6 +148,7 @@ class EclWriter typedef typename GET_PROP_TYPE(TypeTag, Vanguard) Vanguard; typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid; typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; @@ -180,10 +181,8 @@ public: , eclOutputModule_(simulator, collectToIORank_) { if (collectToIORank_.isIORank()) { - globalGrid_ = simulator_.vanguard().grid(); - globalGrid_.switchToGlobalView(); eclIO_.reset(new Opm::EclipseIO(simulator_.vanguard().eclState(), - Opm::UgGridHelpers::createEclipseGrid(globalGrid_, simulator_.vanguard().eclState().getInputGrid()), + Opm::UgGridHelpers::createEclipseGrid(globalGrid(), simulator_.vanguard().eclState().getInputGrid()), simulator_.vanguard().schedule(), simulator_.vanguard().summaryConfig())); } @@ -206,14 +205,19 @@ public: return *eclIO_; } + const EquilGrid& globalGrid() const + { + return simulator_.vanguard().equilGrid(); + } + void writeInit() { if (collectToIORank_.isIORank()) { std::map > integerVectors; if (collectToIORank_.isParallel()) integerVectors.emplace("MPI_RANK", collectToIORank_.globalRanks()); - auto cartMap = Opm::cartesianToCompressed(globalGrid_.size(0), - Opm::UgGridHelpers::globalCell(globalGrid_)); + auto cartMap = Opm::cartesianToCompressed(globalGrid().size(0), + Opm::UgGridHelpers::globalCell(globalGrid())); eclIO_->writeInitial(computeTrans_(cartMap), integerVectors, exportNncStructure_(cartMap)); } } @@ -274,13 +278,13 @@ public: std::map miscSummaryData; std::map> regionData; eclOutputModule_.outputFipLog(miscSummaryData, regionData, isSubStep); - + bool forceDisableProdOutput = false; bool forceDisableInjOutput = false; bool forceDisableCumOutput = false; eclOutputModule_.outputProdLog(reportStepNum, isSubStep, forceDisableProdOutput); eclOutputModule_.outputInjLog(reportStepNum, isSubStep, forceDisableInjOutput); - eclOutputModule_.outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput); + eclOutputModule_.outputCumLog(reportStepNum, isSubStep, forceDisableCumOutput); std::vector buffer; if (collectToIORank_.isIORank()) { @@ -476,7 +480,7 @@ private: Opm::data::Solution computeTrans_(const std::unordered_map& cartesianToActive) const { - const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper(); + const auto& cartMapper = simulator_.vanguard().equilCartesianIndexMapper(); const auto& cartDims = cartMapper.cartesianDimensions(); const int globalSize = cartDims[0]*cartDims[1]*cartDims[2]; @@ -490,8 +494,8 @@ private: tranz.data[0] = 0.0; } - typedef typename Grid :: LeafGridView GlobalGridView; - const GlobalGridView& globalGridView = globalGrid_.leafGridView(); + typedef typename EquilGrid :: LeafGridView GlobalGridView; + const GlobalGridView& globalGridView = globalGrid().leafGridView(); #if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); @@ -500,7 +504,6 @@ private: ElementMapper globalElemMapper(globalGridView); #endif - const auto& cartesianCellIdx = globalGrid_.globalCell(); const EclTransmissibility* globalTrans; if (!collectToIORank_.isParallel()) @@ -538,9 +541,12 @@ private: continue; // we only need to handle each connection once, thank you. // Ordering of compressed and uncompressed index should be the same - assert(cartesianCellIdx[c1] <= cartesianCellIdx[c2]); - int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]); - int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]); + const int cartIdx1 = cartMapper.cartesianIndex( c1 ); + const int cartIdx2 = cartMapper.cartesianIndex( c2 ); + // Ordering of compressed and uncompressed index should be the same + assert(cartIdx1 <= cartIdx2); + int gc1 = std::min(cartIdx1, cartIdx2); + int gc2 = std::max(cartIdx1, cartIdx2); if (gc2 - gc1 == 1) { tranx.data[gc1] = globalTrans->transmissibility(c1, c2); @@ -595,8 +601,8 @@ private: // Checked when writing NNC transmissibilities from the simulation. std::sort(nncData.begin(), nncData.end(), nncCompare); - typedef typename Grid :: LeafGridView GlobalGridView; - const GlobalGridView& globalGridView = globalGrid_.leafGridView(); + typedef typename EquilGrid :: LeafGridView GlobalGridView; + const GlobalGridView& globalGridView = globalGrid().leafGridView(); #if DUNE_VERSION_NEWER(DUNE_GRID, 2,6) typedef Dune::MultipleCodimMultipleGeomTypeMapper ElementMapper; ElementMapper globalElemMapper(globalGridView, Dune::mcmgElementLayout()); @@ -621,7 +627,10 @@ private: globalTrans = &(simulator_.vanguard().globalTransmissibility()); } - auto cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); + // Cartesian index mapper for the serial I/O grid + const auto& equilCartMapper = simulator_.vanguard().equilCartesianIndexMapper(); + + const auto& cartDims = simulator_.vanguard().cartesianIndexMapper().cartesianDimensions(); auto elemIt = globalGridView.template begin(); const auto& elemEndIt = globalGridView.template end(); for (; elemIt != elemEndIt; ++ elemIt) { @@ -641,11 +650,8 @@ private: if (c1 > c2) continue; // we only need to handle each connection once, thank you. - // TODO (?): use the cartesian index mapper to make this code work - // with grids other than Dune::CpGrid. The problem is that we need - // the a mapper for the sequential grid, not for the distributed one. - std::size_t cc1 = globalGrid_.globalCell()[c1]; - std::size_t cc2 = globalGrid_.globalCell()[c2]; + std::size_t cc1 = equilCartMapper.cartesianIndex( c1 ); //globalIOGrid_.globalCell()[c1]; + std::size_t cc2 = equilCartMapper.cartesianIndex( c2 ); //globalIOGrid_.globalCell()[c2]; if ( cc2 < cc1 ) std::swap(cc1, cc2); @@ -733,7 +739,6 @@ private: CollectDataToIORankType collectToIORank_; EclOutputBlackOilModule eclOutputModule_; std::unique_ptr eclIO_; - Grid globalGrid_; std::unique_ptr taskletRunner_; Scalar restartTimeStepSize_; diff --git a/opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp b/opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp index 5d095e67b..68974aae6 100644 --- a/opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp +++ b/opm/simulators/linalg/ExtractParallelGridInformationToISTL.hpp @@ -39,6 +39,11 @@ namespace Opm void extractParallelGridInformationToISTL(const Dune::CpGrid& grid, boost::any& anyComm); +// Grid is not CpGrid --> do nothing. +template +void extractParallelGridInformationToISTL(const Grid&, boost::any&) +{} + } // end namespace Opm #endif //defined(HAVE_OPM_GRID) #endif // OPM_EXTRACTPARALLELGRIDINFORMATIONTOISTL_HEADER_INCLUDED diff --git a/opm/simulators/linalg/ISTLSolverEbos.hpp b/opm/simulators/linalg/ISTLSolverEbos.hpp index 64b11c9a2..56b45dfd1 100644 --- a/opm/simulators/linalg/ISTLSolverEbos.hpp +++ b/opm/simulators/linalg/ISTLSolverEbos.hpp @@ -240,7 +240,7 @@ protected: { parameters_.template init(); extractParallelGridInformationToISTL(simulator_.vanguard().grid(), parallelInformation_); - auto gridForConn = simulator_.vanguard().grid(); + const auto& gridForConn = simulator_.vanguard().grid(); const auto wellsForConn = simulator_.vanguard().schedule().getWellsatEnd(); const bool useWellConn = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions); @@ -648,21 +648,21 @@ protected: /// Create sparsity pattern of matrix without off-diagonal ghost entries. void noGhostAdjacency() { - auto grid = simulator_.vanguard().grid(); + const auto& grid = simulator_.vanguard().grid(); typedef typename Matrix::size_type size_type; - size_type numCells = grid.numCells(); + size_type numCells = grid.size( 0 ); noGhostMat_.reset(new Matrix(numCells, numCells, Matrix::random)); std::vector> pattern; - pattern.resize(grid.numCells()); + pattern.resize(numCells); - auto lid = grid.localIdSet(); + const auto& lid = grid.localIdSet(); const auto& gridView = grid.leafGridView(); auto elemIt = gridView.template begin<0>(); const auto& elemEndIt = gridView.template end<0>(); //Loop over cells - for (; elemIt != elemEndIt; ++elemIt) + for (; elemIt != elemEndIt; ++elemIt) { const auto& elem = *elemIt; size_type idx = lid.id(elem); @@ -679,7 +679,7 @@ protected: } else { auto isend = gridView.iend(elem); - for (auto is = gridView.ibegin(elem); is!=isend; ++is) + for (auto is = gridView.ibegin(elem); is!=isend; ++is) { //check if face has neighbor if (is->neighbor()) @@ -725,7 +725,7 @@ protected: /// Copy interior rows to noghost matrix void copyJacToNoGhost(const Matrix& jac, Matrix& ng) { - //Loop over precalculated interior rows. + //Loop over precalculated interior rows. for (auto row = interiorRows_.begin(); row != interiorRows_.end(); row++ ) { //Copy row diff --git a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp index c1b69e121..676e5278b 100644 --- a/opm/simulators/linalg/findOverlapRowsAndColumns.hpp +++ b/opm/simulators/linalg/findOverlapRowsAndColumns.hpp @@ -30,7 +30,7 @@ namespace detail { /// \brief Find cell IDs for wells contained in local grid. /// - /// Cell IDs of wells stored in a graph, so it can be used to create + /// Cell IDs of wells stored in a graph, so it can be used to create /// an adjacency pattern. Only relevant when the UseWellContribusion option is set to true /// \tparam The type of the DUNE grid. /// \tparam Well vector type @@ -41,17 +41,19 @@ namespace detail template void setWellConnections(const Grid& grid, const W& wells, bool useWellConn, std::vector>& wellGraph) { - if ( grid.comm().size() > 1) + if ( grid.comm().size() > 1) { - wellGraph.resize(grid.numCells()); + Dune::CartesianIndexMapper< Grid > cartMapper( grid ); + const int numCells = cartMapper.compressedSize(); // grid.numCells() + wellGraph.resize(numCells); if (useWellConn) { - const auto& cpgdim = grid.logicalCartesianSize(); + const auto& cpgdim = cartMapper.cartesianDimensions(); std::vector cart(cpgdim[0]*cpgdim[1]*cpgdim[2], -1); - for( int i=0; i < grid.numCells(); ++i ) - cart[grid.globalCell()[i]] = i; + for( int i=0; i < numCells; ++i ) + cart[ cartMapper.cartesianIndex( i ) ] = i; Dune::cpgrid::WellConnections well_indices; well_indices.init(wells, cpgdim, cart); @@ -66,7 +68,7 @@ namespace detail wellGraph[*perf].insert(*perf2); wellGraph[*perf2].insert(*perf); } - } + } } } } @@ -81,27 +83,27 @@ namespace detail /// \param overlapRows List where overlap rows are stored. /// \param interiorRows List where overlap rows are stored. template - void findOverlapAndInterior(const Grid& grid, std::vector& overlapRows, + void findOverlapAndInterior(const Grid& grid, std::vector& overlapRows, std::vector& interiorRows) { //only relevant in parallel case. - if ( grid.comm().size() > 1) + if ( grid.comm().size() > 1) { //Numbering of cells - auto lid = grid.localIdSet(); + const auto& lid = grid.localIdSet(); const auto& gridView = grid.leafGridView(); auto elemIt = gridView.template begin<0>(); const auto& elemEndIt = gridView.template end<0>(); //loop over cells in mesh - for (; elemIt != elemEndIt; ++elemIt) - { + for (; elemIt != elemEndIt; ++elemIt) + { const auto& elem = *elemIt; int lcell = lid.id(elem); - + if (elem.partitionType() != Dune::InteriorEntity) - { + { //add row to list overlapRows.push_back(lcell); } else { diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 68872997a..7044e78ce 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -293,7 +293,7 @@ namespace Opm { // Compute reservoir volumes for RESV controls. rateConverter_.reset(new RateConverterType (phase_usage_, - std::vector(number_of_cells_, 0))); + std::vector(number_of_cells_, 0))); rateConverter_->template defineState(ebosSimulator_); // update VFP properties @@ -555,7 +555,7 @@ namespace Opm { // for ecl compatible restart the current controls are not written const auto& ioCfg = eclState().getIOConfig(); - const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST(); + const auto ecl_compatible_rst = ioCfg.getEclCompatibleRST(); if (true || ecl_compatible_rst) { // always set the control from the schedule for (int w = 0; w :: actionOnBrokenConstraints(const Group& group, const Group::ExceedAction& exceed_action, const Group::ProductionCMode& newControl, const int reportStepIdx, Opm::DeferredLogger& deferred_logger) { - auto& well_state = well_state_; + auto& well_state = well_state_; const Group::ProductionCMode& oldControl = well_state.currentProductionGroupControl(group.name()); std::ostringstream ss;