From b875ba812244abfd13c608ddb8f2b7ddf77876cd Mon Sep 17 00:00:00 2001 From: Antonella Ritorto Date: Tue, 1 Oct 2024 16:22:46 +0200 Subject: [PATCH] LevelCartesianIndexMapper for AluGrid --- CMakeLists_files.cmake | 1 + .../flow/AluGridCartesianIndexMapper.hpp | 13 -- .../flow/AluGridLevelCartesianIndexMapper.hpp | 113 ++++++++++++++++++ opm/simulators/flow/AluGridVanguard.hpp | 10 ++ opm/simulators/flow/FlowProblem.hpp | 2 +- opm/simulators/flow/GenericCpGridVanguard.cpp | 7 ++ opm/simulators/flow/GenericCpGridVanguard.hpp | 8 ++ .../flow/PolyhedralGridVanguard.hpp | 9 ++ .../utils/satfunc/RelpermDiagnostics.cpp | 21 ++-- .../utils/satfunc/RelpermDiagnostics.hpp | 9 +- tests/test_relpermdiagnostics.cpp | 7 +- 11 files changed, 170 insertions(+), 30 deletions(-) create mode 100644 opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4c6360e9f..6825a3c2c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -780,6 +780,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/models/utils/timerguard.hh opm/simulators/flow/ActionHandler.hpp opm/simulators/flow/AluGridCartesianIndexMapper.hpp + opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp opm/simulators/flow/AluGridVanguard.hpp opm/simulators/flow/Banners.hpp opm/simulators/flow/BaseAquiferModel.hpp diff --git a/opm/simulators/flow/AluGridCartesianIndexMapper.hpp b/opm/simulators/flow/AluGridCartesianIndexMapper.hpp index fbca4c569..7984b1ea6 100644 --- a/opm/simulators/flow/AluGridCartesianIndexMapper.hpp +++ b/opm/simulators/flow/AluGridCartesianIndexMapper.hpp @@ -193,10 +193,6 @@ public: int compressedSize() const { return cartesianIndex_.size(); } - /** \brief return number of cells in the active grid. Only for unifying calls with CpGrid and PolyhedralGrid specializations. */ - int compressedLevelZeroSize() const - { return cartesianIndex_.size(); } - /** \brief return index of the cells in the logical Cartesian grid */ int cartesianIndex(const int compressedElementIndex) const { @@ -237,15 +233,6 @@ public: throw std::invalid_argument("cartesianCoordinate not implemented for dimension " + std::to_string(dimension)); } - /** \brief Only for unifying calls with CartesianIndexMapper where levels are relevant */ - void cartesianCoordinateLevel(const int compressedElementIndex, std::array& coords, int level) const - { - if (level) { - throw std::invalid_argument("Invalid level.\n"); - } - cartesianCoordinate(compressedElementIndex, coords); - } - template std::unique_ptr > dataHandle(const GridView& gridView) { diff --git a/opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp b/opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp new file mode 100644 index 000000000..04d3a7e3d --- /dev/null +++ b/opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp @@ -0,0 +1,113 @@ +//=========================================================================== +// +// File: AluGridLevelCartesianIndexMapper.hpp +// +// Created: Tue October 01 09:44:00 2024 +// +// Author(s): Antonella Ritorto +// +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + Copyright 2024 Equinor ASA. + + 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 3 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 . +*/ +#ifndef OPM_ALUGRIDLEVELCARTESIANINDEXMAPPER_HPP +#define OPM_ALUGRIDLEVELCARTESIANINDEXMAPPER_HPP + +#include +#include +#include + +#include +#include +#include + + +namespace Opm +{ + +template<> +class LevelCartesianIndexMapper> +{ + +#if HAVE_MPI + using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>; +#else + using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>; +#endif //HAVE_MPI + + public: + static const int dimension = 3 ; + + LevelCartesianIndexMapper(const Grid& grid, + const Dune::CartesianIndexMapper& cartesianIndexMapper) + : grid_{&grid} + , cartesianIndexMapper_{std::make_unique>(cartesianIndexMapper)} + {} + + const std::array& cartesianDimensions(int level) const + { + throwIfLevelPositive(level); + return cartesianIndexMapper_ ->cartesianDimensions(); + } + + int cartesianSize(int level) const + { + throwIfLevelPositive(level); + return cartesianIndexMapper_->cartesianSize(); + } + + int compressedSize(int level) const + { + throwIfLevelPositive(level); + return cartesianIndexMapper_-> compressedSize(); + } + + + int cartesianIndex( const int compressedElementIndex, const int level) const + { + throwIfLevelPositive(level);; + return cartesianIndexMapper_->cartesianIndex(compressedElementIndex); + } + + void cartesianCoordinate(const int compressedElementIndex, std::array& coords, int level) const + { + throwIfLevelPositive(level); + cartesianIndexMapper_->cartesianCoordinate(compressedElementIndex, coords); + } + + private: + const Grid* grid_; + std::unique_ptr> cartesianIndexMapper_; + + void throwIfLevelPositive(int level) const + { + if (level) { + throw std::invalid_argument("Invalid level.\n"); + } + } +}; + +} + +#endif diff --git a/opm/simulators/flow/AluGridVanguard.hpp b/opm/simulators/flow/AluGridVanguard.hpp index 5646db1b4..2f65caf5f 100644 --- a/opm/simulators/flow/AluGridVanguard.hpp +++ b/opm/simulators/flow/AluGridVanguard.hpp @@ -34,10 +34,12 @@ #include #include +#include #include #include +#include #include #include #include @@ -106,6 +108,7 @@ public: using EquilGrid = GetPropType; using GridView = GetPropType; using CartesianIndexMapper = Dune::CartesianIndexMapper; + using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper; using EquilCartesianIndexMapper = Dune::CartesianIndexMapper; using TransmissibilityType = Transmissibility; using Factory = Dune::FromToGridFactory; @@ -232,6 +235,13 @@ public: const CartesianIndexMapper& cartesianIndexMapper() const { return *cartesianIndexMapper_; } + /*! + * \brief Returns the object which maps a global element index of the simulation grid + * to the corresponding element index of the level logically Cartesian index. + */ + const LevelCartesianIndexMapper levelCartesianIndexMapper() const + { return LevelCartesianIndexMapper(*grid_, *cartesianIndexMapper_); } + /*! * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid */ diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index 95e0a449f..e96aa0b58 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -253,7 +253,7 @@ public: // temporary measure. RelpermDiagnostics relpermDiagnostics{}; relpermDiagnostics.diagnosis(simulator.vanguard().eclState(), - simulator.vanguard().cartesianIndexMapper()); + simulator.vanguard().levelCartesianIndexMapper()); } } diff --git a/opm/simulators/flow/GenericCpGridVanguard.cpp b/opm/simulators/flow/GenericCpGridVanguard.cpp index 81e00e0f6..203d72b9e 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.cpp +++ b/opm/simulators/flow/GenericCpGridVanguard.cpp @@ -604,6 +604,13 @@ GenericCpGridVanguard::cartesianIndexMapper() con return *cartesianIndexMapper_; } +template +const LevelCartesianIndexMapper +GenericCpGridVanguard::levelCartesianIndexMapper() const +{ + return LevelCartesianIndexMapper(*grid_); +} + template const Dune::CartesianIndexMapper& GenericCpGridVanguard::equilCartesianIndexMapper() const diff --git a/opm/simulators/flow/GenericCpGridVanguard.hpp b/opm/simulators/flow/GenericCpGridVanguard.hpp index 62248c745..1921720e5 100644 --- a/opm/simulators/flow/GenericCpGridVanguard.hpp +++ b/opm/simulators/flow/GenericCpGridVanguard.hpp @@ -28,6 +28,7 @@ #define OPM_GENERIC_CPGRID_VANGUARD_HPP #include +#include #include @@ -78,6 +79,7 @@ template class GenericCpGridVanguard { protected: using CartesianIndexMapper = Dune::CartesianIndexMapper; + using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper; using Element = typename GridView::template Codim<0>::Entity; public: @@ -131,6 +133,12 @@ public: */ const CartesianIndexMapper& cartesianIndexMapper() const; + /*! + * \brief Returns the object which maps a global element index of the simulation grid + * to the corresponding element index of the level logically Cartesian index. + */ + const LevelCartesianIndexMapper levelCartesianIndexMapper() const; + /*! * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid */ diff --git a/opm/simulators/flow/PolyhedralGridVanguard.hpp b/opm/simulators/flow/PolyhedralGridVanguard.hpp index cfdb698ce..aea3cb0b7 100644 --- a/opm/simulators/flow/PolyhedralGridVanguard.hpp +++ b/opm/simulators/flow/PolyhedralGridVanguard.hpp @@ -28,6 +28,7 @@ #define OPM_POLYHEDRAL_GRID_VANGUARD_HPP #include +#include #include @@ -93,6 +94,7 @@ public: using EquilGrid = GetPropType; using GridView = GetPropType; using CartesianIndexMapper = Dune::CartesianIndexMapper; + using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper; using EquilCartesianIndexMapper = Dune::CartesianIndexMapper; static constexpr int dimension = Grid::dimension; static constexpr int dimensionworld = Grid::dimensionworld; @@ -169,6 +171,13 @@ public: const CartesianIndexMapper& cartesianIndexMapper() const { return *cartesianIndexMapper_; } + /*! + * \brief Returns the object which maps a global element index of the simulation grid + * to the corresponding element index of the level logically Cartesian index. + */ + const LevelCartesianIndexMapper levelCartesianIndexMapper() const + { return LevelCartesianIndexMapper(*grid_); } + /*! * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid * diff --git a/opm/simulators/utils/satfunc/RelpermDiagnostics.cpp b/opm/simulators/utils/satfunc/RelpermDiagnostics.cpp index 111c603fe..7aceb718d 100644 --- a/opm/simulators/utils/satfunc/RelpermDiagnostics.cpp +++ b/opm/simulators/utils/satfunc/RelpermDiagnostics.cpp @@ -45,10 +45,13 @@ #include #include +#include +#include #ifdef HAVE_DUNE_ALUGRID #include #include +#include #endif // HAVE_DUNE_ALUGRID namespace Opm { @@ -822,9 +825,9 @@ namespace Opm { } } - template + template void RelpermDiagnostics::diagnosis(const EclipseState& eclState, - const CartesianIndexMapper& cartesianIndexMapper) + const LevelCartesianIndexMapper& levelCartesianIndexMapper) { OpmLog::info("\n===============Saturation Functions Diagnostics===============\n"); bool doDiagnostics = phaseCheck_(eclState); @@ -833,16 +836,16 @@ namespace Opm { satFamilyCheck_(eclState); tableCheck_(eclState); unscaledEndPointsCheck_(eclState); - scaledEndPointsCheck_(eclState, cartesianIndexMapper); + scaledEndPointsCheck_(eclState, levelCartesianIndexMapper); } - template + template void RelpermDiagnostics::scaledEndPointsCheck_(const EclipseState& eclState, - const CartesianIndexMapper& cartesianIndexMapper) + const LevelCartesianIndexMapper& levelCartesianIndexMapper) { // All end points are subject to round-off errors, checks should account for it const float tolerance = 1e-6; - const int nc = cartesianIndexMapper.compressedLevelZeroSize(); + const int nc = levelCartesianIndexMapper.compressedSize(0); const bool threepoint = eclState.runspec().endpointScaling().threepoint(); scaledEpsInfo_.resize(nc); EclEpsGridProperties epsGridProperties(eclState, false); @@ -852,7 +855,7 @@ namespace Opm { std::string cellIdx; { std::array ijk; - cartesianIndexMapper.cartesianCoordinateLevel(c, ijk, 0); + levelCartesianIndexMapper.cartesianCoordinate(c, ijk, 0); cellIdx = "(" + std::to_string(ijk[0]) + ", " + std::to_string(ijk[1]) + ", " + std::to_string(ijk[2]) + ")"; @@ -887,8 +890,8 @@ namespace Opm { } #define INSTANCE_DIAGNOSIS(...) \ - template void RelpermDiagnostics::diagnosis>(const EclipseState&, const Dune::CartesianIndexMapper<__VA_ARGS__>&); \ - template void RelpermDiagnostics::scaledEndPointsCheck_>(const EclipseState&, const Dune::CartesianIndexMapper<__VA_ARGS__>&); + template void RelpermDiagnostics::diagnosis>(const EclipseState&, const Opm::LevelCartesianIndexMapper<__VA_ARGS__>&); \ + template void RelpermDiagnostics::scaledEndPointsCheck_>(const EclipseState&, const Opm::LevelCartesianIndexMapper<__VA_ARGS__>&); INSTANCE_DIAGNOSIS(Dune::CpGrid) INSTANCE_DIAGNOSIS(Dune::PolyhedralGrid<3,3>) diff --git a/opm/simulators/utils/satfunc/RelpermDiagnostics.hpp b/opm/simulators/utils/satfunc/RelpermDiagnostics.hpp index b594d21cd..4ea30d145 100644 --- a/opm/simulators/utils/satfunc/RelpermDiagnostics.hpp +++ b/opm/simulators/utils/satfunc/RelpermDiagnostics.hpp @@ -28,6 +28,7 @@ namespace Opm { class EclipseState; + template class LevelCartesianIndexMapper; class MiscTable; class MsfnTable; class SgcwmisTable; @@ -54,9 +55,9 @@ namespace Opm { ///output if they're found. ///\param[in] eclState eclipse state. ///\param[in] grid unstructured grid. - template + template void diagnosis(const EclipseState& eclState, - const CartesianIndexMapper& cartesianIndexMapper); + const LevelCartesianIndexMapper& levelCartesianIndexMapper); private: enum FluidSystem { @@ -95,9 +96,9 @@ namespace Opm { ///Check endpoints in the saturation tables. void unscaledEndPointsCheck_(const EclipseState& eclState); - template + template void scaledEndPointsCheck_(const EclipseState& eclState, - const CartesianIndexMapper& cartesianIndexMapper); + const LevelCartesianIndexMapper& levelCartesianIndexMapper); ///For every table, need to deal with case by case. void swofTableCheck_(const SwofTable& swofTables, diff --git a/tests/test_relpermdiagnostics.cpp b/tests/test_relpermdiagnostics.cpp index ad92b7b19..d7df0ac54 100644 --- a/tests/test_relpermdiagnostics.cpp +++ b/tests/test_relpermdiagnostics.cpp @@ -41,6 +41,7 @@ #include #include +#include #include #include @@ -84,12 +85,12 @@ BOOST_AUTO_TEST_CASE(diagnosis) /*flipNormals=*/false, /*clipZ=*/false); - typedef Dune::CartesianIndexMapper CartesianIndexMapper; - CartesianIndexMapper cartesianIndexMapper = CartesianIndexMapper(grid); + typedef Opm::LevelCartesianIndexMapper LevelCartesianIndexMapper; + LevelCartesianIndexMapper levelCartesianIndexMapper = LevelCartesianIndexMapper(grid); std::shared_ptr counterLog = std::make_shared(Log::DefaultMessageTypes); OpmLog::addBackend( "COUNTERLOG" , counterLog ); RelpermDiagnostics diagnostics; - diagnostics.diagnosis(eclState, cartesianIndexMapper); + diagnostics.diagnosis(eclState, levelCartesianIndexMapper); BOOST_CHECK_EQUAL(1, counterLog->numMessages(Log::MessageType::Warning)); } BOOST_AUTO_TEST_SUITE_END()