From 242a500f52f027a4e115d5ae12c3c9e394eaeb9d Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 5 Jan 2016 19:32:36 +0100 Subject: [PATCH] ebos: replace the macro mess to select the grid type by specialized grid manager classes i.e., there is now a base class for the EclGridManagers and one class for each type of grid (Dune::ALUGrid, Dune::PolyhedralGrid and Dune::CpGrid). Selecting the concrete grid type is now done by deriving the EclProblem's type tag from the type tag of the respective grid manager. --- applications/ebos/eclalugridmanager.hh | 194 +++++++++++++++++++ applications/ebos/eclbasegridmanager.hh | 230 +++++++++++++++++++++++ applications/ebos/eclcpgridmanager.hh | 172 +++++++++++++++++ applications/ebos/eclproblem.hh | 37 +++- applications/ebos/ecltransmissibility.hh | 2 - 5 files changed, 624 insertions(+), 11 deletions(-) create mode 100644 applications/ebos/eclalugridmanager.hh create mode 100644 applications/ebos/eclbasegridmanager.hh create mode 100644 applications/ebos/eclcpgridmanager.hh diff --git a/applications/ebos/eclalugridmanager.hh b/applications/ebos/eclalugridmanager.hh new file mode 100644 index 000000000..073d513a5 --- /dev/null +++ b/applications/ebos/eclalugridmanager.hh @@ -0,0 +1,194 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 . +*/ +/*! + * \file + * \copydoc Ewoms::EclAluGridManager + */ +#ifndef EWOMS_ECL_ALU_GRID_MANAGER_HH +#define EWOMS_ECL_ALU_GRID_MANAGER_HH + +#include "eclbasegridmanager.hh" +#include "alucartesianindexmapper.hh" + +#include +#include +#include + +namespace Ewoms { +template +class EclAluGridManager; + +namespace Properties { +NEW_TYPE_TAG(EclAluGridManager, INHERITS_FROM(EclBaseGridManager)); + +// declare the properties +SET_TYPE_PROP(EclAluGridManager, GridManager, Ewoms::EclAluGridManager); +SET_TYPE_PROP(EclAluGridManager, Grid, Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>); +SET_TYPE_PROP(EclAluGridManager, EquilGrid, Dune::CpGrid); +} // namespace Properties + +/*! + * \ingroup EclBlackOilSimulator + * + * \brief Helper class for grid instantiation of ECL file-format using problems. + * + * This class uses Dune::ALUGrid as the simulation grid. + */ +template +class EclAluGridManager : public EclBaseGridManager +{ + friend class EclBaseGridManager; + typedef EclBaseGridManager ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + +public: + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid; + +private: + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef Grid* GridPointer; + typedef EquilGrid* EquilGridPointer; + typedef Ewoms::AluCartesianIndexMapper CartesianIndexMapper; + typedef CartesianIndexMapper* CartesianIndexMapperPointer; + + static const int dimension = Grid::dimension; + +public: + /*! + * \brief Inherit the constructors from the base class. + */ + using EclBaseGridManager::EclBaseGridManager; + + ~EclAluGridManager() + { + delete cartesianIndexMapper_; + delete grid_; + delete equilGrid_; + } + + /*! + * \brief Return a reference to the simulation grid. + */ + Grid& grid() + { return *grid_; } + + /*! + * \brief Return a reference to the simulation grid. + */ + const Grid& grid() const + { return *grid_; } + + /*! + * \brief Returns a refefence to the grid which should be used by the EQUIL + * initialization code. + * + * The EQUIL keyword is used to specify the initial condition of the reservoir in + * hydrostatic equilibrium. Since the code which does this is not accepting arbitrary + * DUNE grids (the code is part of the opm-core module), this is not necessarily the + * same as the grid which is used for the actual simulation. + */ + const EquilGrid& equilGrid() const + { return *equilGrid_; } + + /*! + * \brief Indicates that the initial condition has been computed and the memory used + * by the EQUIL grid can be released. + * + * Depending on the implementation, subsequent accesses to the EQUIL grid lead to + * crashes. + */ + void releaseEquilGrid() + { + delete equilGrid_; + equilGrid_ = 0; + } + + /*! + * \brief Distribute the simulation grid over multiple processes + * + * (For parallel simulation runs.) + */ + void loadBalance() + { + auto gridView = grid().leafGridView(); + auto dataHandle = cartesianIndexMapper_->dataHandle(gridView); + grid().loadBalance(*dataHandle); + + // communicate non-interior cells values + grid().communicate(*dataHandle, + Dune::InteriorBorder_All_Interface, + Dune::ForwardCommunication ); + } + + /*! + * \brief Returns the object which maps a global element index of the simulation grid + * to the corresponding element index of the logically Cartesian index. + */ + const CartesianIndexMapper& cartesianIndexMapper() const + { return *cartesianIndexMapper_; } + +protected: + void createGrids_() + { + std::vector porv = this->eclState()->getDoubleGridProperty("PORV")->getData(); + + // we use separate grid objects: one for the calculation of the initial condition + // via EQUIL and one for the actual simulation. The reason is that the EQUIL code + // cannot cope with arbitrary Dune grids and is also allergic to distributed + // grids. + + ///// + // create the EQUIL grid + ///// + equilGrid_ = new EquilGrid(); + equilGrid_->processEclipseFormat(this->eclState()->getEclipseGrid(), + /*isPeriodic=*/false, + /*flipNormals=*/false, + /*clipZ=*/false, + porv); + + ///// + // create the simulation grid + ///// + cartesianCellId_ = equilGrid_->globalCell(); + + for (unsigned i = 0; i < dimension; ++i) + cartesianDimension_[i] = equilGrid_->logicalCartesianSize()[i]; + + Dune::FromToGridFactory factory; + grid_ = factory.convert(*equilGrid_, cartesianCellId_); + + cartesianIndexMapper_ = + new CartesianIndexMapper(*grid_, cartesianDimension_, cartesianCellId_); + } + + GridPointer grid_; + EquilGridPointer equilGrid_; + std::vector cartesianCellId_; + std::array cartesianDimension_; + CartesianIndexMapperPointer cartesianIndexMapper_; +}; + +} // namespace Ewoms + +#endif diff --git a/applications/ebos/eclbasegridmanager.hh b/applications/ebos/eclbasegridmanager.hh new file mode 100644 index 000000000..8deeb9eb8 --- /dev/null +++ b/applications/ebos/eclbasegridmanager.hh @@ -0,0 +1,230 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 . +*/ +/*! + * \file + * \copydoc Ewoms::EclBaseGridManager + */ +#ifndef EWOMS_ECL_BASE_GRID_MANAGER_HH +#define EWOMS_ECL_BASE_GRID_MANAGER_HH + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Ewoms { +template +class EclBaseGridManager; + +namespace Properties { +NEW_TYPE_TAG(EclBaseGridManager); + +// declare the properties required by the for the ecl grid manager +NEW_PROP_TAG(Grid); +NEW_PROP_TAG(EquilGrid); +NEW_PROP_TAG(Scalar); +NEW_PROP_TAG(EclDeckFileName); + +SET_STRING_PROP(EclBaseGridManager, EclDeckFileName, "ECLDECK.DATA"); +} // namespace Properties + +/*! + * \ingroup EclBlackOilSimulator + * + * \brief Helper class for grid instantiation of ECL file-format using problems. + */ +template +class EclBaseGridManager : public BaseGridManager +{ + typedef BaseGridManager ParentType; + typedef typename GET_PROP_TYPE(TypeTag, GridManager) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + +public: + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + +protected: + static const int dimension = Grid::dimension; + +public: + /*! + * \brief Register all run-time parameters for the grid manager. + */ + static void registerParameters() + { + EWOMS_REGISTER_PARAM(TypeTag, std::string, EclDeckFileName, + "The name of the file which contains the ECL deck to be simulated"); + } + + /*! + * \brief Create the grid for problem data files which use the ECL file format. + * + * This is the file format used by the commercial ECLiPSE simulator. Usually it uses + * a cornerpoint description of the grid. + */ + EclBaseGridManager(Simulator &simulator) + : ParentType(simulator) + { + std::string fileName = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName); + + // compute the base name of the input file name + const char directorySeparator = '/'; + long int i; + for (i = fileName.size(); i >= 0; -- i) + if (fileName[i] == directorySeparator) + break; + std::string baseName = fileName.substr(i + 1, fileName.size()); + + // remove the extension from the input file + for (i = baseName.size(); i >= 0; -- i) + if (baseName[i] == '.') + break; + std::string rawCaseName; + if (i < 0) + rawCaseName = baseName; + else + rawCaseName = baseName.substr(0, i); + + // transform the result to ALL_UPPERCASE + caseName_ = ""; + for (size_t i = 0; i < rawCaseName.size(); ++i) + caseName_ += std::toupper(rawCaseName[i]); + + Opm::ParserPtr parser(new Opm::Parser()); + typedef std::pair ParseModePair; + typedef std::vector ParseModePairs; + ParseModePairs tmp; + tmp.push_back(ParseModePair(Opm::ParseMode::PARSE_RANDOM_SLASH , Opm::InputError::IGNORE)); + Opm::ParseMode parseMode(tmp); + std::cout << "Reading the deck file '" << fileName << "'" << std::endl; + deck_ = parser->parseFile(fileName , parseMode); + eclState_.reset(new Opm::EclipseState(deck_, parseMode)); + + asImp_().createGrids_(); + + asImp_().finalizeInit_(); + } + + /*! + * \brief Return a pointer to the parsed ECL deck + */ + Opm::DeckConstPtr deck() const + { return deck_; } + + /*! + * \brief Return a pointer to the internalized ECL deck + */ + Opm::EclipseStateConstPtr eclState() const + { return eclState_; } + + /*! + * \brief Return a pointer to the internalized schedule of the ECL deck + */ + Opm::ScheduleConstPtr schedule() const + { return eclState_->getSchedule(); } + + /*! + * \brief Return a pointer to the EclipseGrid object + * + * The EclipseGrid class is provided by the opm-parser module and is used to + * internalize the cornerpoint grid representation and, amongst others, can be used + * to write EGRID files (which tends to be difficult with a plain Dune::CpGrid) + */ + Opm::EclipseGridConstPtr eclGrid() const + { return eclState_->getEclipseGrid(); } + + /*! + * \brief Returns the name of the case. + * + * i.e., the all-uppercase version of the file name from which the + * deck is loaded with the ".DATA" suffix removed. + */ + const std::string& caseName() const + { return caseName_; } + + /*! + * \brief Returns the number of logically Cartesian cells in each direction + */ + const std::array& cartesianDimensions() const + { return asImp_().cartesianIndexMapper().cartesianDimensions(); } + + /*! + * \brief Returns the overall number of cells of the logically Cartesian grid + */ + int cartesianSize() const + { return asImp_().cartesianIndexMapper().cartesianSize(); } + + /*! + * \brief Returns the Cartesian cell id for identifaction with ECL data + */ + unsigned cartesianIndex(unsigned compressedCellIdx) const + { return asImp_().cartesianIndexMapper().cartesianIndex(compressedCellIdx); } + + /*! + * \brief Return the index of the cells in the logical Cartesian grid + */ + unsigned cartesianIndex(const std::array& coords) const + { + unsigned cartIndex = coords[0]; + int factor = cartesianDimensions()[0]; + for (unsigned i = 1; i < dimension; ++i) { + cartIndex += coords[i]*factor; + factor *= cartesianDimensions()[i]; + } + + return cartIndex; + } + + /*! + * \brief Extract Cartesian index triplet (i,j,k) of an active cell. + * + * \param [in] cellIdx Active cell index. + * \param [out] ijk Cartesian index triplet + */ + void cartesianCoordinate(unsigned cellIdx, std::array& ijk) const + { return asImp_().cartesianIndexMapper().cartesianCoordinate(cellIdx, ijk); } + +private: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + + std::string caseName_; + Opm::DeckConstPtr deck_; + Opm::EclipseStateConstPtr eclState_; + +}; + +} // namespace Ewoms + +#endif diff --git a/applications/ebos/eclcpgridmanager.hh b/applications/ebos/eclcpgridmanager.hh new file mode 100644 index 000000000..427bae14d --- /dev/null +++ b/applications/ebos/eclcpgridmanager.hh @@ -0,0 +1,172 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 . +*/ +/*! + * \file + * \copydoc Ewoms::EclCpGridManager + */ +#ifndef EWOMS_ECL_CP_GRID_MANAGER_HH +#define EWOMS_ECL_CP_GRID_MANAGER_HH + +#include "eclbasegridmanager.hh" + +#include + +namespace Ewoms { +template +class EclCpGridManager; + +namespace Properties { +NEW_TYPE_TAG(EclCpGridManager, INHERITS_FROM(EclBaseGridManager)); + +// declare the properties +SET_TYPE_PROP(EclCpGridManager, GridManager, Ewoms::EclCpGridManager); +SET_TYPE_PROP(EclCpGridManager, Grid, Dune::CpGrid); +SET_TYPE_PROP(EclCpGridManager, EquilGrid, typename GET_PROP_TYPE(TypeTag, Grid)); +} // namespace Properties + +/*! + * \ingroup EclBlackOilSimulator + * + * \brief Helper class for grid instantiation of ECL file-format using problems. + * + * This class uses Dune::CpGrid as the simulation grid. + */ +template +class EclCpGridManager : public EclBaseGridManager +{ + friend class EclBaseGridManager; + typedef EclBaseGridManager ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + +public: + typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid; + typedef typename GET_PROP_TYPE(TypeTag, EquilGrid) EquilGrid; + +private: + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef Grid* GridPointer; + typedef EquilGrid* EquilGridPointer; + typedef Dune::CartesianIndexMapper CartesianIndexMapper; + typedef CartesianIndexMapper* CartesianIndexMapperPointer; + +public: + /*! + * \brief Inherit the constructors from the base class. + */ + using EclBaseGridManager::EclBaseGridManager; + + ~EclCpGridManager() + { + delete cartesianIndexMapper_; + delete grid_; + delete equilGrid_; + } + + /*! + * \brief Return a reference to the simulation grid. + */ + Grid& grid() + { return *grid_; } + + /*! + * \brief Return a reference to the simulation grid. + */ + const Grid& grid() const + { return *grid_; } + + /*! + * \brief Returns a refefence to the grid which should be used by the EQUIL + * initialization code. + * + * The EQUIL keyword is used to specify the initial condition of the reservoir in + * hydrostatic equilibrium. Since the code which does this is not accepting arbitrary + * DUNE grids (the code is part of the opm-core module), this is not necessarily the + * same as the grid which is used for the actual simulation. + */ + const EquilGrid& equilGrid() const + { return *equilGrid_; } + + /*! + * \brief Indicates that the initial condition has been computed and the memory used + * by the EQUIL grid can be released. + * + * Depending on the implementation, subsequent accesses to the EQUIL grid lead to + * crashes. + */ + void releaseEquilGrid() + { + delete equilGrid_; + equilGrid_ = 0; + } + + /*! + * \brief Distribute the simulation grid over multiple processes + * + * (For parallel simulation runs.) + */ + void loadBalance() + { + // distribute the grid and switch to the distributed view + grid_->loadBalance(); + grid_->switchToDistributedView(); + cartesianIndexMapper_ = new CartesianIndexMapper(*grid_); + } + + /*! + * \brief Returns the object which maps a global element index of the simulation grid + * to the corresponding element index of the logically Cartesian index. + */ + const CartesianIndexMapper& cartesianIndexMapper() const + { return *cartesianIndexMapper_; } + +protected: + void createGrids_() + { + std::vector porv = this->eclState()->getDoubleGridProperty("PORV")->getData(); + + grid_ = new Dune::CpGrid(); + grid_->processEclipseFormat(this->eclState()->getEclipseGrid(), + /*isPeriodic=*/false, + /*flipNormals=*/false, + /*clipZ=*/false, + porv); + + // we use separate grid objects: one for the calculation of the initial condition + // via EQUIL and one for the actual simulation. The reason is that the EQUIL code + // is allergic to distributed grids and the simulation grid is distributed before + // the initial condition is calculated. + equilGrid_ = new Dune::CpGrid(); + equilGrid_->processEclipseFormat(this->eclState()->getEclipseGrid(), + /*isPeriodic=*/false, + /*flipNormals=*/false, + /*clipZ=*/false, + porv); + } + + GridPointer grid_; + EquilGridPointer equilGrid_; + CartesianIndexMapperPointer cartesianIndexMapper_; +}; + +} // namespace Ewoms + +#endif diff --git a/applications/ebos/eclproblem.hh b/applications/ebos/eclproblem.hh index 8b286fdb7..9ba820f55 100644 --- a/applications/ebos/eclproblem.hh +++ b/applications/ebos/eclproblem.hh @@ -26,9 +26,26 @@ #ifndef EWOMS_ECL_PROBLEM_HH #define EWOMS_ECL_PROBLEM_HH -#include +// make sure that the EBOS_USE_ALUGRID macro. using the preprocessor for this is slightly +// hacky... +#if EBOS_USE_ALUGRID +//#define DISABLE_ALUGRID_SFC_ORDERING 1 +#if !HAVE_DUNE_ALUGRID || !DUNE_VERSION_NEWER(DUNE_ALUGRID, 2,4) +#warning "ALUGrid was indicated to be used for the ECL black oil simulator, but this " +#warning "requires the presence of dune-alugrid >= 2.4. Falling back to Dune::CpGrid" +#undef EBOS_USE_ALUGRID +#define EBOS_USE_ALUGRID 0 +#endif +#else +#define EBOS_USE_ALUGRID 0 +#endif -#include "eclgridmanager.hh" +#if EBOS_USE_ALUGRID +#include "eclalugridmanager.hh" +#else +#include "eclpolyhedralgridmanager.hh" +#include "eclcpgridmanager.hh" +#endif #include "eclwellmanager.hh" #include "eclequilinitializer.hh" #include "eclwriter.hh" @@ -53,9 +70,6 @@ #include #include -// for this simulator to make sense, dune-cornerpoint and opm-parser -// must be available -#include #include #include @@ -73,7 +87,12 @@ template class EclProblem; namespace Properties { -NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclGridManager, EclOutputBlackOil)); +#if EBOS_USE_ALUGRID +NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclAluGridManager, EclOutputBlackOil)); +#else +NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclCpGridManager, EclOutputBlackOil)); +//NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclPolyhedralGridManager, EclOutputBlackOil)); +#endif // Write all solutions for visualization, not just the ones for the // report steps... @@ -841,6 +860,9 @@ private: readExplicitInitialCondition_(); else readEquilInitialCondition_(); + + // release the memory of the EQUIL grid since it's no longer needed after this point + this->simulator().gridManager().releaseEquilGrid(); } void readEquilInitialCondition_() @@ -861,9 +883,6 @@ private: auto &elemFluidState = initialFluidStates_[elemIdx]; elemFluidState.assign(equilInitializer.initialFluidState(elemIdx)); } - - // release the equil grid pointer since it's no longer needed. - this->simulator().gridManager().releaseEquilGrid(); } void readExplicitInitialCondition_() diff --git a/applications/ebos/ecltransmissibility.hh b/applications/ebos/ecltransmissibility.hh index 46d453a69..0b01430b5 100644 --- a/applications/ebos/ecltransmissibility.hh +++ b/applications/ebos/ecltransmissibility.hh @@ -26,8 +26,6 @@ #ifndef EWOMS_ECL_TRANSMISSIBILITY_HH #define EWOMS_ECL_TRANSMISSIBILITY_HH -#include "eclgridmanager.hh" - #include #include #include