From d36ced5448cfb21d6269b4fd26c74fc2f601889c Mon Sep 17 00:00:00 2001 From: Antonella Ritorto Date: Tue, 19 Nov 2024 15:04:03 +0100 Subject: [PATCH] Point global ids for refined level grids --- CMakeLists_files.cmake | 1 + opm/grid/CpGrid.hpp | 109 ++- opm/grid/cpgrid/CpGrid.cpp | 647 +++++++++++------- opm/grid/cpgrid/CpGridData.cpp | 35 + opm/grid/cpgrid/CpGridData.hpp | 25 +- ...ParentToChildCellToPointGlobalIdHandle.hpp | 177 +++++ .../cpgrid/addLgrsOnDistributedGrid_test.cpp | 268 ++++++-- tests/cpgrid/avoidNNCinLGRsCpGrid_test.cpp | 12 + 8 files changed, 957 insertions(+), 317 deletions(-) create mode 100644 opm/grid/cpgrid/ParentToChildCellToPointGlobalIdHandle.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index c1bb7e3a0..c7c878d31 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -202,6 +202,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/grid/LookUpData.hh opm/grid/cpgrid/OrientedEntityTable.hpp opm/grid/cpgrid/ParentToChildrenCellGlobalIdHandle.hpp + opm/grid/cpgrid/ParentToChildCellToPointGlobalIdHandle.hpp opm/grid/cpgrid/PartitionIteratorRule.hpp opm/grid/cpgrid/PartitionTypeIndicator.hpp opm/grid/cpgrid/PersistentContainer.hpp diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp index 1f9b854d2..73d310590 100644 --- a/opm/grid/CpGrid.hpp +++ b/opm/grid/CpGrid.hpp @@ -415,28 +415,12 @@ namespace Dune template cpgrid::Entity entity(const cpgrid::Entity& seed) const; - /// @brief Create a grid out of a coarse one and a refinement(LGR) of a selected block-shaped patch of cells from that coarse grid. - /// - /// Level0 refers to the coarse grid, assumed to be this-> data_[0]. Level1 refers to the LGR (stored in this->data_[1]). - /// LeafView (stored in this-> data_[2]) is built with the level0-entities which weren't involded in the - /// refinenment, together with the new born entities created in level1. - /// Old-corners and old-faces (from coarse grid) lying on the boundary of the patch, get replaced by new-born-equivalent corners - /// and new-born-faces. - /// - /// @param [in] cells_per_dim Number of (refined) cells in each direction that each parent cell should be refined to. - /// @param [in] startIJK Cartesian triplet index where the patch starts. - /// @param [in] endIJK Cartesian triplet index where the patch ends. - /// Last cell part of the lgr will be {endijk[0]-1, ... endIJK[2]-1}. - /// @param [in] lgr_name Name (std::string) for the lgr/level1 - void addLgrUpdateLeafView(const std::array& cells_per_dim, const std::array& startIJK, - const std::array& endIJK, const std::string& lgr_name); - /// @brief Create a grid out of a coarse one and (at most) 2 refinements(LGRs) of selected block-shaped disjoint patches /// of cells from that coarse grid. /// /// Level0 refers to the coarse grid, assumed to be this-> data_[0]. Level1 and level2 refer to the LGRs (stored in this->data_[1] /// data_[2]). LeafView (stored in this-> data_[3]) is built with the level0-entities which weren't involded in the - /// refinenment, together with the new born entities created in level1 and level2. + /// refinenment, together with the new born entities created in level1 and level2. /// Old-corners and old-faces (from coarse grid) lying on the boundary of the patches, get replaced by new-born-equivalent corners /// and new-born-faces. /// @@ -476,7 +460,6 @@ namespace Dune // from certain LGR Dune::cpgrid::Intersection getParentIntersectionFromLgrBoundaryFace(const Dune::cpgrid::Intersection& intersection) const; - /// --------------- Adaptivity (begin) --------------- /// @brief Mark entity for refinement (or coarsening). /// @@ -1083,12 +1066,100 @@ namespace Dune /// @param [in] cells_per_dim: Total children cells in each direction (x-,y-, and z-direction) of the single-cell-refinement. /// @param [in] faceIdxInLgr: Face index in the single-cell-refinement. /// @param [in] elemLgr_ptr: Pointer to the elemLgr single-cell-refinement grid. - /// @param [in] elemLgr: Cell index from starting grid, that has been refined into a single-cell-refinement. + /// @param [in] elemLgr: Cell index from starting grid, that has been refined into a single-cell-refinement. int getParentFaceWhereNewRefinedFaceLiesOn(const std::array& cells_per_dim, int faceIdxInLgr, const std::shared_ptr& elemLgr_ptr, int elemLgr) const; + /// --------------- Auxiliary methods to support Adaptivity (end) --------------- + /// @brief Check if there are non neighboring connections on blocks of cells selected for refinement. + /// + /// @param [in] startIJK_vec Vector of ijk values denoting the start of each block of cells selected for refinement. + /// @param [in] endIJK_vec Vector of ijk values denoting the end of each block of cells selected for refinement. + /// + /// @return True if all blocks of cells do not contain any NNCs (non-neighboring-connection). + /// False if there is a block with at least one cell with a NNC. + bool nonNNCsSelectedCellsLGR( const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec) const; + + /// @brief Mark selected elements and assign them their corresponding level. + /// + /// Given blocks of cells selected for refinement, Mark selected elements and assign them their corresponding + /// (refined) level (grid). When level zero grid is distributed before refinement, detect which LGRs are active + /// in each process. + /// + /// @param [in] startIJK_vec Vector of ijk values denoting the start of each block of cells selected for refinement. + /// @param [in] endIJK_vec Vector of ijk values denoting the end of each block of cells selected for refinement. + /// @param [out] assignRefinedLevel Assign level for the refinement of each marked cell. Example: refined element from + /// LGR1 have level 1, refined element rfom LGR2 have level 2, etc. + /// @param [out] lgr_with_at_least_one_active_cell Determine if an LGR is not empty in a given process, we set + /// lgr_with_at_least_one_active_cell[in that level] to 1 if it contains + /// at least one active cell, and to 0 otherwise. + void markElemAssignLevelDetectActiveLgrs(const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec, + std::vector& assignRefinedLevel, + std::vector& lgr_with_at_least_one_active_cell); + + /// @brief Predict minimum cell and point global ids per process. + /// + /// Predict how many new cells/points (born in refined level grids) need new globalIds, so we can assign unique + /// new ids ( and anticipate the maximum). At this point, the grid is already refined according to the LGR specification. + /// + /// @param [in] assignRefinedLevel Assign level for the refinement of each marked cell. Example: refined element from + /// LGR1 have level 1, refined element rfom LGR2 have level 2, etc. + /// @param [in] cells_per_dim_vec Total child cells in each direction (x-,y-, and z-direction) per block of cells. + /// @param [in] lgr_with_at_least_one_active_cell Determine if an LGR is not empty in a given process: + /// lgr_with_at_least_one_active_cell[level] = 1 if it contains + /// at least one active cell in the current process, and 0 otherwise. + /// @param [out] min_globalId_cell_in_proc + /// @param [out] min_globalId_point_in_proc + void predictMinCellAndPointGlobalIdPerProcess(const std::vector& assignRefinedLevel, + const std::vector>& cells_per_dim_vec, + const std::vector& lgr_with_at_least_one_active_cell, + int& min_globalId_cell_in_proc, + int& min_globalId_point_in_proc) const; + + /// @brief Assign cell global ids of new born cell from refined level grids. Assign 'candidate' point global ids + /// for points in refined level grids. + /// + /// @param [out] localToGlobal_cells_per_level Relation local element.index() to assigned cell global id. + /// @param [out] localToGlobal_points_per_level Relation local point.index() to assigned 'candidate' global id. + /// @param [in] min_globalId_cell_in_proc Minimum cell global id per process. + /// @param [in] min_globalId_point_in_proc Minimum point global id per process. + /// @param [in] cells_per_dim_vec Total child cells in each direction (x-,y-, and z-direction) per block of cells. + void assignCellIdsAndCandidatePointIds( std::vector>& localToGlobal_cells_per_level, + std::vector>& localToGlobal_points_per_level, + int min_globalId_cell_in_proc, + int min_globalId_point_in_proc, + const std::vector>& cells_per_dim_vec) const; + + /// @brief Select and re-write point global ids. + /// + /// After assigning global IDs to points in refined-level grids, a single point may have + /// "multiple unique" global IDs, one in each process to which it belongs. + /// To reduce the unnucesary id assigments, since global IDs must be distinct across the global leaf view + /// and consistent across each refined-level grid, we will rewrite the entries in + /// localToGlobal_points_per_level. Using cell_to_point_ across all refined cells through + /// communication: gathering the 8 corner points of each interior cell and scattering the + /// 8 corner points of overlapping cells, for all child cells of a parent cell in level zero grid. + /// + /// @param [out] localToGlobal_points_per_level Relation local point.index() to assigned 'candidate' global id. + /// @param [in] parent_to_children The communication step is based on level zero grid, via the relation parent-children-cells. + /// @param [in] cells_per_dim_vec Total child cells in each direction (x-,y-, and z-direction) per block of cells. + void selectWinnerPointIds(std::vector>& localToGlobal_points_per_level, + const std::vector>>& parent_to_children, + const std::vector>& cells_per_dim_vec) const; + + /// @brief For a grid whose level zero has been distributed and then locally refined, populate the cell_index_set_ of each refined level grid. + void populateCellIndexSetRefinedGrid(int level); + + /// @brief For a grid whose level zero has been distributed and then locally refined, populate the cell_index_set_ of the leaf grid view. + void populateCellIndexSetLeafGridView(); + + /// @brief For a grid whose level zero has been distributed and then locally refined, populate the global_id_set_ of the leaf grid view. + void populateLeafGlobalIdSet(); + public: /// @brief Compute for each level grid, a map from the global_cell_[ cell index in level grid ] to the leaf index of the equivalent cell diff --git a/opm/grid/cpgrid/CpGrid.cpp b/opm/grid/cpgrid/CpGrid.cpp index 0c32640b1..72d6ba927 100644 --- a/opm/grid/cpgrid/CpGrid.cpp +++ b/opm/grid/cpgrid/CpGrid.cpp @@ -51,6 +51,7 @@ #include "../CpGrid.hpp" #include "ParentToChildrenCellGlobalIdHandle.hpp" +#include "ParentToChildCellToPointGlobalIdHandle.hpp" #include #include //#include @@ -1212,6 +1213,328 @@ Dune::cpgrid::Intersection CpGrid::getParentIntersectionFromLgrBoundaryFace(cons OPM_THROW(std::invalid_argument, "Face is on the boundary of the grid"); } +bool CpGrid::nonNNCsSelectedCellsLGR( const std::vector>& startIJK_vec, const std::vector>& endIJK_vec) const +{ + // Non neighboring connections: Currently, adding LGRs whose cells have NNCs is not supported yet. + // Find out which (ACTIVE) elements belong to the block cells defined by startIJK and endIJK values + // and check if they have NNCs. + // Note: at this point, the (level zero) grid might be already distributed, before adding the LGRs. Therefore, + // we get the active index of each element and its corresponding ijk from the global grid, to be able to determine + // if the cell bolengs to certain block of cells selected for refinement (comparing element's ijk with start/endIJK values). + // It is not correct to make the comparasion with element.index() and minimum/maximum index of each block of cell, since + // element.index() is local and the block of cells are defined with global values. + for(const auto& element: elements(this->leafGridView())) { + std::array ijk; + getIJK(element.index(), ijk); + for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { + bool belongsToLevel = ( (ijk[0] >= startIJK_vec[level][0]) && (ijk[0] < endIJK_vec[level][0]) ) + && ( (ijk[1] >= startIJK_vec[level][1]) && (ijk[1] < endIJK_vec[level][1]) ) + && ( (ijk[2] >= startIJK_vec[level][2]) && (ijk[2] < endIJK_vec[level][2]) ); + if(belongsToLevel) { + // Check that the cell to be marked for refinement has no NNC (no neighbouring connections). + if (this->currentData().back()->hasNNCs({element.index()})){ + return false; + } + } + } + } + return true; +} + +void CpGrid::markElemAssignLevelDetectActiveLgrs(const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec, + std::vector& assignRefinedLevel, + std::vector& lgr_with_at_least_one_active_cell) +{ + // Find out which (ACTIVE) elements belong to the block cells defined by startIJK and endIJK values. + for(const auto& element: elements(this->leafGridView())) { + std::array ijk; + getIJK(element.index(), ijk); + for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { + bool belongsToLevel = true; + for (int c = 0; c < 3; ++c) { + belongsToLevel = belongsToLevel && ( (ijk[c] >= startIJK_vec[level][c]) && (ijk[c] < endIJK_vec[level][c]) ); + if (!belongsToLevel) + break; + } + if(belongsToLevel) { + this-> mark(1, element); + assignRefinedLevel[element.index()] = level+1; // shifted since starting grid is level 0, and refined grids levels are >= 1. + lgr_with_at_least_one_active_cell[level] = 1; + } + } + } +} + +void CpGrid::predictMinCellAndPointGlobalIdPerProcess([[maybe_unused]] const std::vector& assignRefinedLevel, + [[maybe_unused]] const std::vector>& cells_per_dim_vec, + [[maybe_unused]] const std::vector& lgr_with_at_least_one_active_cell, + [[maybe_unused]] int& min_globalId_cell_in_proc, + [[maybe_unused]] int& min_globalId_point_in_proc) const +{ +#if HAVE_MPI + // Maximum global id from level zero. (Then, new entities get global id values greater than max_globalId_levelZero). + // Recall that only cells and points are taken into account; faces are ignored (do not have any global id). + auto max_globalId_levelZero = comm().max(current_data_->front()->global_id_set_->getMaxGlobalId()); + + // Predict how many new cell ids per process are needed. + std::vector cell_ids_needed_by_proc(comm().size()); + std::size_t local_cell_ids_needed = 0; + for ( const auto& element : elements( levelGridView(0), Dune::Partitions::interior) ) { + // Get old mark (from level zero). After calling adapt, all marks are set to zero. + bool hasBeenMarked = currentData().front()->getMark(element) == 1; + if ( hasBeenMarked ) { + const auto& level = assignRefinedLevel[element.index()]; + // Shift level (to level -1) since cells_per_dim_vec stores number of subdivisions in each direction (xyz) + // per parent cell, per level, starting from level 1, ..., maxLevel. + local_cell_ids_needed += cells_per_dim_vec[level-1][0]*cells_per_dim_vec[level-1][1]*cells_per_dim_vec[level-1][2]; + } + } + comm().allgather(&local_cell_ids_needed, 1, cell_ids_needed_by_proc.data()); + + // Overestimate ('predict') how many new point ids per process are needed. + // Assign for all partition type points a 'candidate of global id' (unique in each process). + std::vector point_ids_needed_by_proc(comm().size()); + std::size_t local_point_ids_needed = 0; + for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level){ + if(lgr_with_at_least_one_active_cell[level-1]>0) { + // Amount of local_point_ids_needed might be overestimated. + for (const auto& point : vertices(levelGridView(level))){ + // If point coincides with an existing corner from level zero, then it does not need a new global id. + if ( !(*current_data_)[level]->corner_history_.empty() ) { + const auto& bornLevel_bornIdx = (*current_data_)[level]->corner_history_[point.index()]; + if (bornLevel_bornIdx[0] == -1) { // Corner is new-> it needs a new(candidate) global id + local_point_ids_needed += 1; + } + } + } + } + } + comm().allgather(&local_point_ids_needed, 1, point_ids_needed_by_proc.data()); + + auto expected_max_globalId_cell = std::accumulate(cell_ids_needed_by_proc.begin(), + cell_ids_needed_by_proc.end(), + max_globalId_levelZero + 1); + min_globalId_cell_in_proc = std::accumulate(cell_ids_needed_by_proc.begin(), + cell_ids_needed_by_proc.begin()+comm().rank(), + max_globalId_levelZero + 1); + min_globalId_point_in_proc = std::accumulate(point_ids_needed_by_proc.begin(), + point_ids_needed_by_proc.begin()+ comm().rank(), + expected_max_globalId_cell+1); +#endif +} + +void CpGrid::assignCellIdsAndCandidatePointIds( std::vector>& localToGlobal_cells_per_level, + std::vector>& localToGlobal_points_per_level, + int min_globalId_cell_in_proc, + int min_globalId_point_in_proc, + const std::vector>& cells_per_dim_vec ) const +{ + for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level) { + localToGlobal_cells_per_level[level-1].resize((*current_data_)[level]-> size(0)); + localToGlobal_points_per_level[level-1].resize((*current_data_)[level]-> size(3)); + // Notice that in general, (*current_data_)[level]-> size(0) != local owned cells/points. + + // Global ids for cells (for owned cells) + for (const auto& element : elements(levelGridView(level))) { + // At his point, partition_type_indicator_ of refined level grids is not set. However, for refined cells, + // element.partitionType() returns the partition type of the parent cell. Therefore, all child cells of + // an interior/overlap parent cell are also interior/overlap. + if (element.partitionType() == InteriorEntity) { + localToGlobal_cells_per_level[level - 1][element.index()] = min_globalId_cell_in_proc; + ++min_globalId_cell_in_proc; + } + } + for (const auto& point : vertices(levelGridView(level))) { + // Check if it coincides with a corner from level zero. In that case, no global id is needed. + const auto& bornLevel_bornIdx = (*current_data_)[level]->corner_history_[point.index()]; + if (bornLevel_bornIdx[0] != -1) { // Corner in the refined grid coincides with a corner from level 0. + // Therefore, search and assign the global id of the previous existing equivalent corner. + const auto& equivPoint = cpgrid::Entity<3>(*( (*current_data_)[bornLevel_bornIdx[0]]), bornLevel_bornIdx[1], true); + localToGlobal_points_per_level[level-1][point.index()] = current_data_->front()->global_id_set_->id( equivPoint ); + } + else { + // Assign CANDIDATE global id to (all partition type) points that do not coincide with + // any corners from level zero. + // TO DO after invoking this method: make a final decision on a unique id for points of refined level grids, + // via a communication step. + localToGlobal_points_per_level[level-1][point.index()] = min_globalId_point_in_proc; + ++min_globalId_point_in_proc; + } + } + } +} + +void CpGrid::selectWinnerPointIds([[maybe_unused]] std::vector>& localToGlobal_points_per_level, + [[maybe_unused]] const std::vector>>& parent_to_children, + [[maybe_unused]] const std::vector>& cells_per_dim_vec) const +{ +#if HAVE_MPI + // To store cell_to_point_ information of all refined level grids. + std::vector>> level_cell_to_point(cells_per_dim_vec.size()); + // To decide which "candidate" point global id wins, the rank is stored. The smallest ranks wins, + // i.e., the other non-selected candidates get rewritten with the values from the smallest (winner) rank. + std::vector> level_winning_ranks(cells_per_dim_vec.size()); + + for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level) { + + level_cell_to_point[level -1] = currentData()[level]->cell_to_point_; + // Set std::numeric_limits::max() to make sure that, during communication, the rank of the interior cell + // wins (int between 0 and comm().size()). + level_winning_ranks[level-1].resize(currentData()[level]->size(3), std::numeric_limits::max()); + + for (const auto& element : elements(levelGridView(level))) { + // For interior cells, rewrite the rank value - later used in "point global id competition". + if (element.partitionType() == InteriorEntity) { + for (const auto& corner : currentData()[level]->cell_to_point_[element.index()]){ + int rank = comm().rank(); + level_winning_ranks[level -1][corner] = rank; + } + } + } + } + ParentToChildCellToPointGlobalIdHandle parentToChildCellToPointGlobalId_handle(comm(), + parent_to_children, + level_cell_to_point, + level_winning_ranks, + localToGlobal_points_per_level); + currentData().front()->communicate(parentToChildCellToPointGlobalId_handle, + Dune::InteriorBorder_All_Interface, + Dune::ForwardCommunication ); +#endif +} + +void CpGrid::populateCellIndexSetRefinedGrid([[maybe_unused]] int level) +{ +#if HAVE_MPI + const auto& level_global_id_set = (*current_data_)[level]->global_id_set_; + auto& level_index_set = currentData()[level]->cellIndexSet(); + + level_index_set.beginResize(); + + // ParallelIndexSet::LocalIndex( elemIdx, attribute /* owner or copy */, true/false) + // The bool indicates whether the global index might exist on other processes with other attributes. + // RemoteIndices::rebuild has a Boolean as a template parameter telling the method whether to ignore this + // Boolean on the indices or not when building. + // + // For refined level grids, we check if the parent cell is fully interior. Then, its children won't be seen + // by any other process. Therefore, the boolean is set to false. + + for(const auto& element : elements(levelGridView(level))) { + if ( element.partitionType() == InteriorEntity) { + + // Check if it has an overlap neighbor + bool parentFullyInterior = true; + const auto& parent_cell = element.father(); + + for (const auto& intersection : intersections(levelGridView(parent_cell.level()), parent_cell)) { + if ( intersection.neighbor() ) { + const auto& neighborPartitionType = intersection.outside().partitionType(); + // To help detection of fully interior cells, i.e., without overlap neighbors + if (neighborPartitionType == OverlapEntity ) { + parentFullyInterior = false; + // Interior cells with overlap neighbor may appear in other processess -> true + level_index_set.add( level_global_id_set->id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), true)); + // Store it only once + break; + } + } + } + if(parentFullyInterior) { // Fully interior cells do not appear in any other process -> false + level_index_set.add( level_global_id_set->id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), false)); + } + } + else { // overlap cell + assert(element.partitionType() == OverlapEntity); + level_index_set.add( level_global_id_set->id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::copy), true)); + } + } + + level_index_set.endResize(); + + currentData()[level]->cellRemoteIndices().template rebuild(); // Do not ignore bool in ParallelIndexSet! +#endif +} + +void CpGrid::populateCellIndexSetLeafGridView() +{ +#if HAVE_MPI + auto& leaf_index_set = (*current_data_).back()->cellIndexSet(); + + leaf_index_set.beginResize(); + + // ParallelIndexSet::LocalIndex( elemIdx, attribute /* owner or copy */, true/false) + // The bool indicates whether the global index might exist on other processes with other attributes. + // RemoteIndices::rebuild has a Boolean as a template parameter telling the method whether to ignore this + // Boolean on the indices or not when building. + for(const auto& element : elements(leafGridView())) { + const auto& elemPartitionType = element.getEquivLevelElem().partitionType(); + if ( elemPartitionType == InteriorEntity) { + // Check if it has an overlap neighbor + bool isFullyInterior = true; + for (const auto& intersection : intersections(leafGridView(), element)) { + if ( intersection.neighbor() ) { + const auto& neighborPartitionType = intersection.outside().getEquivLevelElem().partitionType(); + // To help detection of fully interior cells, i.e., without overlap neighbors + if (neighborPartitionType == OverlapEntity ) { + isFullyInterior = false; + // Interior cells with overlap neighbor may appear in other processess -> false + leaf_index_set.add(globalIdSet().id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), true)); + // Store it only once + break; + } + } + } + if(isFullyInterior) { // Fully interior cells do not appear in any other process -> false + leaf_index_set.add(globalIdSet().id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), false)); + } + } + else { // overlap cell + assert(elemPartitionType == OverlapEntity); + leaf_index_set.add(globalIdSet().id(element), + ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::copy), true)); + } + } + leaf_index_set.endResize(); + + (*current_data_).back()->cellRemoteIndices().template rebuild(); + +#endif +} + +void CpGrid::populateLeafGlobalIdSet() +{ + // Global id for the cells in leaf grid view + std::vector leafCellIds(current_data_->back()->size(0)); + for(const auto& element: elements(leafGridView())){ + // Notice that for level zero cells the global_id_set_ is given, for refined level grids was defined + // under the assumption of each lgr being fully contained in the interior of a process. + // Therefore, it is not needed here to distingish between owned and overlap cells. + auto equivElem = element.getEquivLevelElem(); + leafCellIds[element.index()] = (*current_data_)[element.level()]->global_id_set_->id(equivElem); + } + + // Global id for the faces in leaf grid view. Empty vector (Entity<1> not supported for CpGrid). + std::vector leafFaceIds{}; + + // Global id for the points in leaf grid view + std::vector leafPointIds(current_data_->back()->size(3)); + for(const auto& point : vertices(leafGridView())){ + const auto& level_pointLevelIdx = current_data_->back()->corner_history_[point.index()]; + assert(level_pointLevelIdx[0] != -1); + assert(level_pointLevelIdx[1] != -1); + const auto& pointLevelEntity = cpgrid::Entity<3>(*( (*current_data_)[level_pointLevelIdx[0]]), level_pointLevelIdx[1], true); + leafPointIds[point.index()] = (*current_data_)[level_pointLevelIdx[0]]->global_id_set_->id(pointLevelEntity); + } + + current_data_->back()->global_id_set_->swap(leafCellIds, leafFaceIds, leafPointIds); +} + double CpGrid::cellCenterDepth(int cell_index) const { // Here cell center depth is computed as a raw average of cell corner depths. @@ -2048,19 +2371,13 @@ void CpGrid::postAdapt() current_view_data_ -> postAdapt(); } -void CpGrid::addLgrUpdateLeafView(const std::array& cells_per_dim, const std::array& startIJK, - const std::array& endIJK, const std::string& lgr_name) -{ - this -> addLgrsUpdateLeafView({cells_per_dim}, {startIJK}, {endIJK}, {lgr_name}); -} - - void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_per_dim_vec, const std::vector>& startIJK_vec, const std::vector>& endIJK_vec, const std::vector& lgr_name_vec) { - // For parallel run, level zero grid is stored in distributed_data_[0]. If CpGrid::scatterGrid has been invoked, then current_view_data_ == distributed_data_[0]. + // For parallel run, level zero grid is stored in distributed_data_[0]. If CpGrid::scatterGrid has been invoked, + // then current_view_data_ == distributed_data_[0]. // For serial run, level zero grid is stored in data_[0]. In this case, current_view_data_ == data_[0]. // Note: currentData() returns data_ (if grid is not distributed) or distributed_data_ otherwise. @@ -2082,37 +2399,10 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p // active/inactive cells, instead, relies on "ijk-computations". // TO DO: improve/remove. // To check "Compatibility of numbers of subdivisions of neighboring LGRs". - bool compatibleSubdivisionsHasFailed = false; - if (startIJK_vec.size() > 1) { - bool notAllowedYet = false; - for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { - for (std::size_t otherLevel = level+1; otherLevel < startIJK_vec.size(); ++otherLevel) { - const auto& sharedFaceTag = - current_view_data_->sharedFaceTag({startIJK_vec[level], startIJK_vec[otherLevel]}, {endIJK_vec[level],endIJK_vec[otherLevel]}); - if(sharedFaceTag == -1){ - break; // Go to the next "other patch" - } - if (sharedFaceTag == 0 ) { - notAllowedYet = notAllowedYet || - ((cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); - } - if (sharedFaceTag == 1) { - notAllowedYet = notAllowedYet || - ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); - } - if (sharedFaceTag == 2) { - notAllowedYet = notAllowedYet || - ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1])); - } - if (notAllowedYet){ - compatibleSubdivisionsHasFailed = true; - break; - } - } // end-otherLevel-for-loop - } // end-level-for-loop - }// end-if-patchesShareFace - compatibleSubdivisionsHasFailed = comm().max(compatibleSubdivisionsHasFailed); - if(compatibleSubdivisionsHasFailed) { + // The method compatibleSubdivision returns a bool. We convert it into an int since MPI within DUNE does not support bool directly. + int compatibleSubdivisions = current_view_data_->compatibleSubdivisions(cells_per_dim_vec, startIJK_vec, endIJK_vec); + compatibleSubdivisions = comm().min(compatibleSubdivisions); // 0 when at least one process returns false (un-compatible subdivisions). + if(!compatibleSubdivisions) { if (comm().rank()==0){ OPM_THROW(std::logic_error, "Subdivisions of neighboring LGRs sharing at least one face do not coincide. Not suppported yet."); } @@ -2120,44 +2410,27 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p OPM_THROW_NOLOG(std::logic_error, "Subdivisions of neighboring LGRs sharing at least one face do not coincide. Not suppported yet."); } } - + // Non neighboring connections: Currently, adding LGRs whose cells have NNCs is not supported yet. - // To check "Non-NNCs (non neighboring connections)" for all processes. - bool nonNNCsHasFailed = false; + // The method nonNNCsSelectedCellsLGR returns a bool. We convert it into an int since MPI within DUNE does not support bool directly. + int nonNNCs = this->nonNNCsSelectedCellsLGR(startIJK_vec, endIJK_vec); + // To check "Non-NNCs (non non neighboring connections)" for all processes. + nonNNCs = comm().min(nonNNCs); // 0 when at least one process returns false (there are NNCs on a selected cell for refinement). + if(!nonNNCs) { + OPM_THROW(std::logic_error, "NNC face on a cell containing LGR is not supported yet."); + } + // To determine if an LGR is not empty in a given process, we set // lgr_with_at_least_one_active_cell[in that level] to 1 if it contains // at least one active cell, and to 0 otherwise. std::vector lgr_with_at_least_one_active_cell(startIJK_vec.size()); // Determine the assigned level for the refinement of each marked cell std::vector assignRefinedLevel(current_view_data_->size(0)); - // Find out which (ACTIVE) elements belong to the block cells defined by startIJK and endIJK values. - for(const auto& element: elements(this->leafGridView())) { - std::array ijk; - getIJK(element.index(), ijk); - for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { - bool belongsToLevel = true; - for (int c = 0; c < 3; ++c) { - belongsToLevel = belongsToLevel && ( (ijk[c] >= startIJK_vec[level][c]) && (ijk[c] < endIJK_vec[level][c]) ); - if (!belongsToLevel) - break; - } - if(belongsToLevel) { - // Check that the cell to be marked for refinement has no NNC (no neighbouring connections). - if (current_view_data_->hasNNCs({element.index()})){ - nonNNCsHasFailed = true; - break; - } - this-> mark(1, element); - assignRefinedLevel[element.index()] = level+1; // shifted since starting grid is level 0, and refined grids levels are >= 1. - lgr_with_at_least_one_active_cell[level] = 1; - } // end-if-belongsToLevel - } // end-level-for-loop - } // end-element-for-loop - - nonNNCsHasFailed = comm().max(nonNNCsHasFailed); - if(nonNNCsHasFailed) { - OPM_THROW(std::logic_error, "NNC face on a cell containing LGR is not supported yet."); - } + + markElemAssignLevelDetectActiveLgrs(startIJK_vec, + endIJK_vec, + assignRefinedLevel, + lgr_with_at_least_one_active_cell); int non_empty_lgrs = 0; for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { @@ -2187,12 +2460,11 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p // - Define ParallelIndex for overlap cells and their neighbors if(comm().size()>1) { #if HAVE_MPI - // Maximum global id from level zero. (Then, new entities get global id values greater than max_globalId_levelZero). - // Recall that only cells and points are taken into account; faces are ignored (do not have any global id). - auto max_globalId_levelZero = comm().max(current_data_->front()->global_id_set_->getMaxGlobalId()); - + // Prediction min cell and point global ids per process + // // Predict how many new cells/points (born in refined level grids) need new globalIds, so we can assign unique // new ids ( and anticipate the maximum). + // The grid is already refined according to the LGR specification // At this point, neither cell_index_set_ nor partition_type_indicator_ are populated. // Refined level grid cells: // 1. Inherit their partition type from their parent cell (i.e., element.father().partitionType()). @@ -2203,55 +2475,19 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p // cell and face partition types are already defined, not available yet for refined level grids. // 1. Assign for all partition type points a 'candidate of global id' (unique in each process). // Except the points that coincide with a point from level zero. - // 2. To make point ids globally unique, re-write the values for points that are corners of overlap - // refined cells, via communication. + // 2. Re-write the values for points that are corners of overlap refined cells, via communication. // Under the assumption of LGRs fully-interior, no communication is needed. In the general case, communication will be used // to populate overlap cell/point global ids on the refined level grids. - - // Predict how many new cell ids per process are needed. - std::vector cell_ids_needed_by_proc(comm().size()); - std::size_t local_cell_ids_needed = 0; - for ( const auto& element : elements( levelGridView(0), Dune::Partitions::interior) ) { - // Get old mark (from level zero). After calling adapt, all marks are set to zero. - bool hasBeenMarked = currentData().front()->getMark(element) == 1; - if ( hasBeenMarked ) { - const auto& level = assignRefinedLevel[element.index()]; - // Shift level (to level -1) since cells_per_dim_vec stores number of subdivisions in each direction (xyz) - // per parent cell, per level, starting from level 1, ..., maxLevel. - local_cell_ids_needed += cells_per_dim_vec[level-1][0]*cells_per_dim_vec[level-1][1]*cells_per_dim_vec[level-1][2]; - } - } - comm().allgather(&local_cell_ids_needed, 1, cell_ids_needed_by_proc.data()); - - // Overestimate ('predict') how many new point ids per process are needed. - // Assign for all partition type points a 'candidate of global id' (unique in each process). - std::vector point_ids_needed_by_proc(comm().size()); - std::size_t local_point_ids_needed = 0; - for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level){ - if(lgr_with_at_least_one_active_cell[level-1]>0) { - // Amount of local_point_ids_needed might be overestimated. - for (const auto& point : vertices(levelGridView(level))){ - // If point coincides with an existing corner from level zero, then it does not need a new global id. - if ( !(*current_data_)[level]->corner_history_.empty() ) { - const auto& bornLevel_bornIdx = (*current_data_)[level]->corner_history_[point.index()]; - if (bornLevel_bornIdx[0] == -1) { // Corner is new-> it needs a new global id - local_point_ids_needed += 1; - } - } - } - } - } - comm().allgather(&local_point_ids_needed, 1, point_ids_needed_by_proc.data()); - - auto expected_max_globalId_cell = std::accumulate(cell_ids_needed_by_proc.begin(), - cell_ids_needed_by_proc.end(), - max_globalId_levelZero + 1); - auto min_globalId_cell_in_proc = std::accumulate(cell_ids_needed_by_proc.begin(), - cell_ids_needed_by_proc.begin()+comm().rank(), - max_globalId_levelZero + 1); - auto min_globalId_point_in_proc= std::accumulate(point_ids_needed_by_proc.begin(), - point_ids_needed_by_proc.begin()+ comm().rank(), - expected_max_globalId_cell); + /** Warning: due to the overlap layer size (equal to 1) cells that share corners or edges (not faces) with interior cells + are not included/seen by the process. This, in some cases, ends up in duplicated point ids. */ + // + int min_globalId_cell_in_proc = 0; + int min_globalId_point_in_proc = 0; + predictMinCellAndPointGlobalIdPerProcess(assignRefinedLevel, + cells_per_dim_vec, + lgr_with_at_least_one_active_cell, + min_globalId_cell_in_proc, + min_globalId_point_in_proc); // Only for level 1,2,.., maxLevel grids. // For each level, define the local-to-global maps for cells and points (for faces: empty). @@ -2262,54 +2498,50 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p // Ignore faces - empty vectors. std::vector> localToGlobal_faces_per_level(cells_per_dim_vec.size()); - for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level) { - localToGlobal_cells_per_level[level-1].resize((*current_data_)[level]-> size(0)); - localToGlobal_points_per_level[level-1].resize((*current_data_)[level]-> size(3)); - // Notice that in general, (*current_data_)[level]-> size(0) != local owned cells/points. - - // Global ids for cells (for owned cells) - for (const auto& element : elements(levelGridView(level))) { - if (element.partitionType() == InteriorEntity) { - localToGlobal_cells_per_level[level - 1][element.index()] = min_globalId_cell_in_proc; - ++min_globalId_cell_in_proc; - } - } - for (const auto& point : vertices(levelGridView(level))) { - // Check if it coincides with a corner from level zero. In that case, no global id is needed. - const auto& bornLevel_bornIdx = (*current_data_)[level]->corner_history_[point.index()]; - if (bornLevel_bornIdx[0] != -1) { // Corner in the refined grid coincides with a corner from level 0. - // Therefore, search and assign the global id of the previous existing equivalent corner. - const auto& equivPoint = cpgrid::Entity<3>(*( (*current_data_)[bornLevel_bornIdx[0]]), bornLevel_bornIdx[1], true); - localToGlobal_points_per_level[level-1][point.index()] = current_data_->front()->global_id_set_->id( equivPoint ); - } - else { - // Assign new global id to (all partition type) points that do not coincide with - // any corners from level zero. - // TO DO: Implement a final decision on a unique id for points of refined level grids. - localToGlobal_points_per_level[level-1][point.index()] = min_globalId_point_in_proc; - ++min_globalId_point_in_proc; - } - } - } + assignCellIdsAndCandidatePointIds(localToGlobal_cells_per_level, + localToGlobal_points_per_level, + min_globalId_cell_in_proc, + min_globalId_point_in_proc, + cells_per_dim_vec); + const auto& parent_to_children = current_data_->front()->parent_to_children_cells_; ParentToChildrenCellGlobalIdHandle parentToChildrenGlobalId_handle(parent_to_children, localToGlobal_cells_per_level); currentData().front()->communicate(parentToChildrenGlobalId_handle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication ); - - // After assigning global IDs to points in refined-level grids, a single point may have + + // After assigning global IDs to points in refined-level grids, a single point may have // a "unique" global ID in each local leaf grid view for every process to which it belongs. - // To ensure true uniqueness, since global IDs must be distinct across the global leaf view - // and consistent across each refined-level grid, we will rewrite the entries in - // localToGlobal_points_per_level. Correction: TO DO. + // To ensure true uniqueness, since global IDs must be distinct across the global leaf view + // and consistent across each refined-level grid, we will rewrite the entries in + // localToGlobal_points_per_level. + // + // This correction is done using cell_to_point_ across all refined cells through + // communication: gathering the 8 corner points of each interior cell and scattering the + // 8 corner points of overlapping cells, for all child cells of a parent cell in level zero grid. + // + // Design decision: Why we communicate via level zero grid instead of in each refined level grid. + // The reason is that how children are stored (the ordering) in parent_to_children_cells_ + // is always the same, accross all processes. + // Even though the ordering of the corners in cell_to_point_ is the same accross all processes, + // this may not be enough to correctly overwrite the "winner" point global ids for refined cells. + // + /** Current approach avoids duplicated point ids when + // 1. the LGR is distributed in P_{i_0}, ..., P_{i_n}, with n+1 < comm().size(), + // AND + // 2. there is no coarse cell seen by a process P with P != P_{i_j}, j = 0, ..., n. + // Otherwise, there will be duplicated point ids. + // + // Reason: neighboring cells that only share corners (not faces) are NOT considered in the + // overlap layer of the process.*/ + selectWinnerPointIds(localToGlobal_points_per_level, + parent_to_children, + cells_per_dim_vec); for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level) { - // For the general case where the LGRs might be also distributed, a communication step is needed to assign global ids - // for overlap cells and points. - // Global id set for each (refined) level grid. - if(lgr_with_at_least_one_active_cell[level-1]>0) { + if(lgr_with_at_least_one_active_cell[level-1]>0) { // Check if LGR is active in currect process. (*current_data_)[level]->global_id_set_->swap(localToGlobal_cells_per_level[level-1], localToGlobal_faces_per_level[level-1], localToGlobal_points_per_level[level-1]); @@ -2317,108 +2549,25 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p } for (std::size_t level = 1; level < cells_per_dim_vec.size()+1; ++level) { - - const auto& level_global_id_set = (*current_data_)[level]->global_id_set_; - auto& level_index_set = currentData()[level]->cellIndexSet(); - - level_index_set.beginResize(); - - for(const auto& element : elements(levelGridView(level))) { - if ( element.partitionType() == InteriorEntity) { - level_index_set.add( level_global_id_set->id(element), - ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), true)); - } - else { // overlap cell - assert(element.partitionType() == OverlapEntity); - level_index_set.add( level_global_id_set->id(element), - ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::copy), true)); - } - } - level_index_set.endResize(); - - currentData()[level]->cellRemoteIndices().template rebuild(); - + + populateCellIndexSetRefinedGrid(level); // Compute the partition type for cell currentData()[level]->computeCellPartitionType(); - // Compute the partition type for point currentData()[level]->computePointPartitionType(); - // Now we can compute the communication interface. currentData()[level]->computeCommunicationInterfaces(currentData()[level]->size(3)); - } // end-for-loop-level - //////////////////////////////// - - // Global id for the cells in leaf grid view - std::vector leafCellIds(current_data_->back()->size(0)); - for(const auto& element: elements(leafGridView())){ - // Notice that for level zero cells the global_id_set_ is given, for refined level grids was defined - // under the assumption of each lgr being fully contained in the interior of a process. - // Therefore, it is not needed here to distingish between owned and overlap cells. - auto equivElem = element.getEquivLevelElem(); - leafCellIds[element.index()] = (*current_data_)[element.level()]->global_id_set_->id(equivElem); - } - leafCellIds.shrink_to_fit(); - - // Global id for the faces in leaf grid view. Empty vector (Entity<1> not supported for CpGrid). - std::vector leafFaceIds{}; - - // Global id for the points in leaf grid view - std::vector leafPointIds(current_data_->back()->size(3)); - for(const auto& point : vertices(leafGridView())){ - const auto& level_pointLevelIdx = current_data_->back()->corner_history_[point.index()]; - assert(level_pointLevelIdx[0] != -1); - assert(level_pointLevelIdx[1] != -1); - const auto& pointLevelEntity = cpgrid::Entity<3>(*( (*current_data_)[level_pointLevelIdx[0]]), level_pointLevelIdx[1], true); - leafPointIds[point.index()] = (*current_data_)[level_pointLevelIdx[0]]->global_id_set_->id(pointLevelEntity); } - leafPointIds.shrink_to_fit(); - - current_data_->back()->global_id_set_->swap(leafCellIds, leafFaceIds, leafPointIds); + populateLeafGlobalIdSet(); this->global_id_set_ptr_ = std::make_shared(*(current_data_->back())); for (std::size_t level = 0; level < cells_per_dim_vec.size()+1; ++level) { this->global_id_set_ptr_->insertIdSet(*(*current_data_)[level]); } - auto& leaf_index_set = (*current_data_).back()->cellIndexSet(); - - leaf_index_set.beginResize(); - - for(const auto& element : elements(leafGridView())) { - const auto& elemPartitionType = element.getEquivLevelElem().partitionType(); - if ( elemPartitionType == InteriorEntity) { - // Check if it has an overlap neighbor - bool isFullyInterior = true; - for (const auto& intersection : intersections(leafGridView(), element)) { - if ( intersection.neighbor() ) { - const auto& neighborPartitionType = intersection.outside().getEquivLevelElem().partitionType(); - // To help detection of fully interior cells, i.e., without overlap neighbors - if (neighborPartitionType == OverlapEntity ) { - isFullyInterior = false; - leaf_index_set.add(globalIdSet().id(element), - ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), true)); - // Store it only once - break; - } - } - } - if(isFullyInterior) { // In case we do not need these indices, then modify/remove the assert below regarding leaf_index_set.size(). - leaf_index_set.add(globalIdSet().id(element), - ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::owner), false)); - } - } - else { // overlap cell - assert(elemPartitionType == OverlapEntity); - leaf_index_set.add(globalIdSet().id(element), - ParallelIndexSet::LocalIndex(element.index(), AttributeSet(AttributeSet::copy), true)); - } - } - leaf_index_set.endResize(); - - (*current_data_).back()->cellRemoteIndices().template rebuild(); - + populateCellIndexSetLeafGridView(); + // Compute the partition type for cell (*current_data_).back()->computeCellPartitionType(); @@ -2427,9 +2576,7 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector>& cells_p // Now we can compute the communication interface. current_data_->back()->computeCommunicationInterfaces(current_data_->back()->size(3)); - assert(static_cast(leaf_index_set.size()) == static_cast(current_data_->back()->size(0)) ); - // Alternatively, in case fully interior cells should't be added in leaf_index_set, the next assert can be used: - //assert(static_cast(leaf_index_set.size()) == static_cast(current_data_->back()->size(0) - /* fully interior cells */) ); + assert(static_cast(current_data_->back()->cellIndexSet().size()) == static_cast(current_data_->back()->size(0)) ); #endif } // end-if-comm().size()>1 diff --git a/opm/grid/cpgrid/CpGridData.cpp b/opm/grid/cpgrid/CpGridData.cpp index 001de5e49..bd5366d54 100644 --- a/opm/grid/cpgrid/CpGridData.cpp +++ b/opm/grid/cpgrid/CpGridData.cpp @@ -2239,6 +2239,41 @@ void CpGridData::validStartEndIJKs(const std::vector>& startIJ } } +bool CpGridData::compatibleSubdivisions(const std::vector>& cells_per_dim_vec, + const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec) const +{ + bool compatibleSubdivisions = true; + if (startIJK_vec.size() > 1) { + bool notAllowedYet = false; + for (std::size_t level = 0; level < startIJK_vec.size(); ++level) { + for (std::size_t otherLevel = level+1; otherLevel < startIJK_vec.size(); ++otherLevel) { + const auto& sharedFaceTag = this-> sharedFaceTag({startIJK_vec[level], startIJK_vec[otherLevel]}, {endIJK_vec[level],endIJK_vec[otherLevel]}); + if(sharedFaceTag == -1){ + break; // Go to the next "other patch" + } + if (sharedFaceTag == 0 ) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); + } + if (sharedFaceTag == 1) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][2] != cells_per_dim_vec[otherLevel][2])); + } + if (sharedFaceTag == 2) { + notAllowedYet = notAllowedYet || + ((cells_per_dim_vec[level][0] != cells_per_dim_vec[otherLevel][0]) || (cells_per_dim_vec[level][1] != cells_per_dim_vec[otherLevel][1])); + } + if (notAllowedYet){ + compatibleSubdivisions = false; + break; + } + } // end-otherLevel-for-loop + } // end-level-for-loop + }// end-if-patchesShareFace + return compatibleSubdivisions; +} + Geometry<3,3> CpGridData::cellifyPatch(const std::array& startIJK, const std::array& endIJK, const std::vector& patch_cells, DefaultGeometryPolicy& cellifiedPatch_geometry, diff --git a/opm/grid/cpgrid/CpGridData.hpp b/opm/grid/cpgrid/CpGridData.hpp index 92fd3f8c7..c4fff7de0 100644 --- a/opm/grid/cpgrid/CpGridData.hpp +++ b/opm/grid/cpgrid/CpGridData.hpp @@ -364,9 +364,9 @@ class CpGridData /// @brief Check startIJK and endIJK of each patch of cells to be refined are valid, i.e. /// startIJK and endIJK vectors have the same size and, startIJK < endIJK coordenate by coordenate. /// - /// @param [in] startIJK_vec Vector of Cartesian triplet indices where each patch starts. - /// @param [in] endIJK_vec Vector of Cartesian triplet indices where each patch ends. - /// Last cell part of the lgr will be {endIJK_vec[patch][0]-1, ..., endIJK_vec[patch][2]-1}. + /// @param [in] startIJK_vec Vector of Cartesian triplet indices where each patch starts. + /// @param [in] endIJK_vec Vector of Cartesian triplet indices where each patch ends. + /// Last cell part of the lgr will be {endIJK_vec[patch][0]-1, ..., endIJK_vec[patch][2]-1}. void validStartEndIJKs(const std::vector>& startIJK_vec, const std::vector>& endIJK_vec) const; /// @brief Check that every cell to be refined has cuboid shape. @@ -412,6 +412,25 @@ class CpGridData void postAdapt(); private: + /// @brief Check compatibility of number of subdivisions of neighboring LGRs. + /// + /// Check shared faces on boundaries of LGRs. Not optimal since the code below does not take into account + /// active/inactive cells, instead, relies on "ijk-computations". + /// + /// @param [in] cells_per_dim_vec Vector of expected subdivisions per cell, per direction, in each LGR. + /// @param [in] startIJK_vec Vector of Cartesian triplet indices where each patch starts. + /// @param [in] endIJK_vec Vector of Cartesian triplet indices where each patch ends. + /// Last cell part of the lgr will be {endIJK_vec[patch][0]-1, ..., endIJK_vec[patch][2]-1}. + /// @return True if all block of cells either do not share faces on their boundaries, or they may share faces with compatible + /// subdivisions. Example: block1 and block2 share an I_FACE, then number of subdivisions NY NZ should coincide, i.e. + /// if block1, block2 cells_per_dim values are {NX1, NY1, NZ1}, {NX2, NY2, NZ2}, respectively, then NY1 == NY2 and + /// NZ1 == Nz2. + /// False if at least two blocks share a face and their subdivions are not compatible. In the example above, + /// if NY1 != NY2 or NZ1 != NZ2. + bool compatibleSubdivisions(const std::vector>& cells_per_dim_vec, + const std::vector>& startIJK_vec, + const std::vector>& endIJK_vec) const; + std::array,8> getReferenceRefinedCorners(int idx_in_parent_cell, const std::array& cells_per_dim) const; /// @brief Compute amount of cells in each direction of a patch of cells. (Cartesian grid required). diff --git a/opm/grid/cpgrid/ParentToChildCellToPointGlobalIdHandle.hpp b/opm/grid/cpgrid/ParentToChildCellToPointGlobalIdHandle.hpp new file mode 100644 index 000000000..14711bea8 --- /dev/null +++ b/opm/grid/cpgrid/ParentToChildCellToPointGlobalIdHandle.hpp @@ -0,0 +1,177 @@ +//=========================================================================== +// +// File: ParentToChildCellToPointGlobalIdHandle.hpp +// +// Created: November 19 2024 +// +// Author(s): Antonella Ritorto +// Markus Blatt +// +// $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_PARENTTOCHILDCELLTOPOINTGLOBALIDHANDLE_HEADER +#define OPM_PARENTTOCHILDCELLTOPOINTGLOBALIDHANDLE_HEADER + +#include +#include + +#include +#include +#include + + +namespace +{ +#if HAVE_MPI + +/// \brief Handle for assignment of point global ids of refined cells. +struct ParentToChildCellToPointGlobalIdHandle { + // - The container used for gather and scatter contains "candidates" point global ids for interior elements of the refined level grids (LGRs). + // Access is done with the local index of its parent cell index and its parent cell list of children. + // level_point_global_ids[ level-1 ][ child_cell_local_index [ corner ] ] = "candidate" point global id, + // when child_cell_local_index belongs to the children_local_index_list: + // parent_to_children_[ element.index() ] = parent_to_children_[ element.index() ] = { level, children_local_index_list } + // and corner = 0, ...,7. + // To decide which "candidate" point global id wins, we use the rank. The smallest ranks wins, + // i.e., the other non-selected candidates get rewritten with the values from the smallest (winner) rank. + // - In the scatter method, the "winner" rank and the 8 point global ids of each number of children) get rewritten. + + using DataType = int; + + /// \param comm Communication + /// \param parent_to_children Map from parent index to all children, and the level they are stored. + /// parent_to_children_[ element.index() ] = { level, children_list local indices } + /// \param level_cell_to_point + /// \param level_point_global_ids A container that for the elements of a level contains all candidate point global ids. + /// \param level_winning_ranks + /// \param level_point_global_ids + ParentToChildCellToPointGlobalIdHandle(const Dune::CpGrid::Communication& comm, + const std::vector>>& parent_to_children, + const std::vector>>& level_cell_to_point, + std::vector>& level_winning_ranks, + std::vector>& level_point_global_ids) + : comm_(comm) + , parent_to_children_(parent_to_children) + , level_cell_to_point_(level_cell_to_point) + , level_winning_ranks_(level_winning_ranks) + , level_point_global_ids_(level_point_global_ids) + {} + + bool fixedSize(std::size_t, std::size_t) + { + // Not every cell has children. When they have children, the amount might vary. + return false; + } + + bool contains(std::size_t, std::size_t codim) + { + // Only communicate values attached to cells. + return codim == 0; + } + + template // T = Entity<0> + std::size_t size(const T& element) + { + // Communicate variable size: 1 (rank) + (8* amount of child cells) from an interior parent cell from level zero grid. + // Skip values that are not interior, or have no children (in that case, 'invalid' level = -1) + const auto& [level, children] = parent_to_children_[element.index()]; + // [Bug in dune-common] VariableSizeCommunicator will deadlock if a process attempts to send a message of size zero. + // This can happen if the size method returns zero for all entities that are shared with another process. + // Therefore, when skipping cells without children or for overlap cells, we set the size to 1. + if ( (element.partitionType() != Dune::InteriorEntity) || (level == -1)) + return 1; + return 1 + ( 8*children.size()); // rank + 8 "winner" point global ids per child cell + } + + // Gather global ids of child cells of a coarse interior parent cell + template // T = Entity<0> + void gather(B& buffer, const T& element) + { + // Skip values that are not interior, or have no children (in that case, 'invalid level' = -1) + const auto& [level, children] = parent_to_children_[element.index()]; + // [Bug in dune-common] VariableSizeCommunicator will deadlock if a process tries to send a message with size zero. + // To avoid this, for cells without children or for overlap cells, we set the size to 1 and write a single DataType + // value (e.g., '42'). + if ( (element.partitionType() != Dune::InteriorEntity) || (level==-1)) { + buffer.write(42); + return; + } + // Store the children's corner global ids in the buffer when the element is interior and has children. + // Write the rank first, for example via the "corner 0" of cell_to_point_ of the first child: + // First child: children[0] + // First corner of first child: level_cell_to_point_[ level -1 ][children[0]] [0] + buffer.write( comm_.rank() ); + for (const auto& child : children) + for (const auto& corner : level_cell_to_point_[level -1][child]) + buffer.write(level_point_global_ids_[level-1][corner]); + } + + // Scatter global ids of child cells of a coarse overlap parent cell + template // T = Entity<0> + void scatter(B& buffer, const T& element, std::size_t size) + { + const auto& [level, children] = parent_to_children_[element.index()]; + // Read all values to advance the pointer used by the buffer to the correct index. + // (Skip overlap-cells-without-children and interior-cells). + if ( ( (element.partitionType() == Dune::OverlapEntity) && (level==-1) ) || (element.partitionType() == Dune::InteriorEntity ) ) { + // Read all values to advance the pointer used by the buffer + // to the correct index + for (std::size_t received_int = 0; received_int < size; ++received_int) { + DataType tmp; + buffer.read(tmp); + } + } + else { // Overlap cell with children. + // Read and store the values in the correct location directly. + // The order of the children is the same on each process. + assert(children.size()>0); + assert(level>0); + // Read and store the values in the correct location directly. + DataType tmp_rank; + buffer.read(tmp_rank); + for (const auto& child : children) { + for (const auto& corner : level_cell_to_point_[level -1][child]) { + auto& min_rank = level_winning_ranks_[level-1][corner]; + // Rewrite the rank (smaller rank wins) + if (tmp_rank < min_rank) { + min_rank = tmp_rank; + auto& target_entry = level_point_global_ids_[level-1][corner]; + buffer.read(target_entry); + } else { + DataType rubbish; + buffer.read(rubbish); + } + } + } + } + } + +private: + const Dune::CpGrid::Communication& comm_; + const std::vector>>& parent_to_children_; + const std::vector>>& level_cell_to_point_; + std::vector>& level_winning_ranks_; + std::vector>& level_point_global_ids_; +}; +#endif // HAVE_MPI +} // namespace +#endif diff --git a/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp b/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp index c4e9cd8b7..d3a523ca2 100644 --- a/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp +++ b/tests/cpgrid/addLgrsOnDistributedGrid_test.cpp @@ -366,29 +366,18 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, // Check global id is not duplicated for interior cells in each refined level grid. std::vector interior_cell_global_ids; int local_interior_cell_count = 0; - interior_cell_global_ids.reserve(coarse_grid.currentData()[level]->size(0)); + interior_cell_global_ids.reserve(data[level]->size(0)); + for (const auto& element: elements(coarse_grid.levelGridView(level), Dune::Partitions::interior)){ - interior_cell_global_ids.push_back(coarse_grid.currentData()[level]->globalIdSet().id(element)); + interior_cell_global_ids.push_back(data[level]->globalIdSet().id(element)); ++local_interior_cell_count; } + auto global_level_interior_cell_count = coarse_grid.comm().sum(local_interior_cell_count); auto [all_level_interior_cell_global_ids, displ] = Opm::allGatherv(interior_cell_global_ids, coarse_grid.comm()); const std::set all_level_interior_cell_global_ids_set(all_level_interior_cell_global_ids.begin(), all_level_interior_cell_global_ids.end()); BOOST_CHECK( all_level_interior_cell_global_ids.size() == all_level_interior_cell_global_ids_set.size() ); BOOST_CHECK( global_level_interior_cell_count == static_cast(all_level_interior_cell_global_ids_set.size()) ); - - // Check global id is not duplicated for interior points in each refined level grid. - std::vector interior_point_global_ids; - int local_interior_point_count = 0; - interior_point_global_ids.reserve(coarse_grid.currentData()[level]->size(3)); - for (const auto& point : vertices(coarse_grid.levelGridView(level), Dune::Partitions::interior)){ - interior_point_global_ids.push_back(coarse_grid.currentData()[level]->globalIdSet().id(point)); - ++local_interior_point_count; - } - auto global_level_interior_point_count = coarse_grid.comm().sum(local_interior_point_count); - auto [all_level_interior_point_global_ids, displ_point] = Opm::allGatherv(interior_point_global_ids, coarse_grid.comm()); - const std::set all_level_interior_point_global_ids_set(all_level_interior_point_global_ids.begin(), all_level_interior_point_global_ids.end()); - BOOST_CHECK( static_cast(global_level_interior_point_count) == all_level_interior_point_global_ids_set.size() ); } // Check global id is not duplicated for interior cells @@ -409,20 +398,9 @@ void refinePatch_and_check(Dune::CpGrid& coarse_grid, BOOST_CHECK( static_cast(allGlobalIds_cells.size()) == global_cells_count); BOOST_CHECK( allGlobalIds_cells.size() == allGlobalIds_cells_set.size() ); - // Check global id is not duplicated for interior points - std::vector localInteriorPointIds_vec; - localInteriorPointIds_vec.reserve(data.back()->size(3)); // more than actually needed since only care about interior points - int local_interior_points_count = 0; - for (const auto& point: vertices(coarse_grid.leafGridView(), Dune::Partitions::interior)) { - localInteriorPointIds_vec.push_back(data.back()->globalIdSet().id(point)); - ++local_interior_points_count; - } - auto global_point_count = coarse_grid.comm().sum(local_interior_points_count); - auto [allGlobalIds_points, displ_point] = Opm::allGatherv(localInteriorPointIds_vec, coarse_grid.comm()); - - const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); - BOOST_CHECK( static_cast(allGlobalIds_points.size()) == global_point_count); - BOOST_CHECK( allGlobalIds_points.size() == allGlobalIds_points_set.size() ); + /** [Bug] Uniqueness of point global ids cannot be checked in general since current code sets overlap layer size equal to 1, + which in particular means that cells that share corners or edges (and not faces) with interior cells are not considered/ + seen by the process. Therefore, depending how the LGRs are distributed, there may be "multiple ids" for the same points.*/ // Local/Global id sets for level grids (level 0, 1, ..., maxLevel). For level grids, local might differ from global id. for (int level = 0; level < coarse_grid.maxLevel() +1; ++level) @@ -520,6 +498,26 @@ BOOST_AUTO_TEST_CASE(threeLgrs) // Leaf grid view total amount of cells: 10x8x8 -6 + 32 + 27 + 64 = 757 refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + // Check global id is not duplicated for points for each LGR + // LGR1 dim 4x2x4 -> 5x3x5 = 75 points + // LGR2 dim 3x3x3 -? 4x4x4 = 64 points + // LGR3 dim 4x4x4 -> 5x5x5 = 125 points + const std::vector& expected_point_ids_per_lgr = { 75, 64, 125}; + for (int lgr = 1; lgr < 4; ++lgr) { + std::vector local_point_ids; + local_point_ids.reserve(expected_point_ids_per_lgr[lgr-1]); + for (const auto& element : elements(grid.levelGridView(lgr))) { + for (int corner = 0; corner < 8; ++corner) + { + const auto& point = element.subEntity<3>(corner); + local_point_ids.push_back(grid.currentData()[lgr]->globalIdSet().id(point)); + } + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + BOOST_CHECK( static_cast(all_point_ids_set.size()) == expected_point_ids_per_lgr[lgr-1]); + } + // Check global id is not duplicated for points std::vector localPointIds_vec; localPointIds_vec.reserve(grid.currentData().back()->size(3)); @@ -567,13 +565,34 @@ BOOST_AUTO_TEST_CASE(atLeastOneLgr_per_process_attempt) // LGR4 element indices = 27, 31 in rank 3.Total 16 refined cells, 45 points (45-12 = 33 with new global id). grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); - // LGR1 dim 1x2x1 (16 refined cells) - // LGR2 dim 1x1x1 (27 refined cells) - // LGR3 dim 1x1x1 (64 refined cells) - // LGR4 dim 1x2x1 (16 refined cells) 3x5x3 + // LGR1 dim 2x4x2 (16 refined cells) + // LGR2 dim 3x3x3 (27 refined cells) + // LGR3 dim 4x4x4 (64 refined cells) + // LGR4 dim 2x4x2 (16 refined cells) // Total global ids in leaf grid view for cells: 36-(6 marked cells) + 16 + 27 + 64 + 16 = 153 // Total global ids in leaf grid view for points: 80 + 33 + 56 + 117 + 33 = 319 refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + + // Check global id is not duplicated for points for each LGR + // LGR1 dim 2x4x2 -> 3x5x3 = 45 points + // LGR2 dim 3x3x3 -> 4x4x4 = 64 points + // LGR3 dim 4x4x4 -> 5x5x5 = 125 points + // LGR4 dim 2x4x2 -> 3x5x3 = 45 points + const std::vector& expected_point_ids_per_lgr = { 45, 64, 125, 45}; + for (int lgr = 1; lgr < 5; ++lgr) { + std::vector local_point_ids; + local_point_ids.reserve(expected_point_ids_per_lgr[lgr-1]); + for (const auto& element : elements(grid.levelGridView(lgr))) { + for (int corner = 0; corner < 8; ++corner) + { + const auto& point = element.subEntity<3>(corner); + local_point_ids.push_back(grid.currentData()[lgr]->globalIdSet().id(point)); + } + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + BOOST_CHECK( static_cast(all_point_ids_set.size()) == expected_point_ids_per_lgr[lgr-1]); + } // Check global id is not duplicated for points std::vector localPointIds_vec; @@ -584,7 +603,6 @@ BOOST_AUTO_TEST_CASE(atLeastOneLgr_per_process_attempt) } auto [allGlobalIds_points, displPoint ] = Opm::allGatherv(localPointIds_vec, grid.comm()); const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); - // Total global ids in leaf grid view for points: 80 + 33 + 56 + 117 + 33 = 319 BOOST_CHECK( allGlobalIds_points_set.size() == 319 ); } @@ -623,6 +641,39 @@ BOOST_AUTO_TEST_CASE(throw_not_fully_interior_lgr) grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + + // Check global id is not duplicated for points for each LGR + // LGR1 dim 2x4x2 -> 3x5x3 = 45 points + // LGR2 dim 3x3x3 -> 4x4x4 = 64 points + // LGR3 dim 4x4x4 -> 5x5x5 = 125 points + // LGR4 dim 2x4x2 -> 3x5x3 = 45 points + const std::vector& expected_point_ids_per_lgr = { 45, 64, 125, 45}; + for (int lgr = 1; lgr < 5; ++lgr) { + std::vector local_point_ids; + local_point_ids.reserve(expected_point_ids_per_lgr[lgr-1]); + for (const auto& element : elements(grid.levelGridView(lgr))) { + for (int corner = 0; corner < 8; ++corner) + { + const auto& point = element.subEntity<3>(corner); + local_point_ids.push_back(grid.currentData()[lgr]->globalIdSet().id(point)); + } + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + BOOST_CHECK( static_cast(all_point_ids_set.size()) == expected_point_ids_per_lgr[lgr-1]); + } + + // Check global id is not duplicated for points + std::vector localPointIds_vec; + localPointIds_vec.reserve(grid.currentData().back()->size(3)); + for (const auto& point : vertices(grid.leafGridView())) { + // Notice that all partition type points are pushed back. Selecting only interior points does not bring us to the expected value. + localPointIds_vec.push_back(grid.currentData().back()->globalIdSet().id(point)); + } + auto [allGlobalIds_points, displPoint ] = Opm::allGatherv(localPointIds_vec, grid.comm()); + const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); + // Total global ids in leaf grid view for points: 80 + 33 + 56 + 117 + 33 = 319 + BOOST_CHECK( allGlobalIds_points_set.size() == 319 ); } } @@ -645,9 +696,6 @@ BOOST_AUTO_TEST_CASE(globalRefine2) BOOST_AUTO_TEST_CASE(distributed_lgr) { - // Only for testing assignment of new global ids for refined entities (cells and point belonging to - // refined level grids). - // Create a grid Dune::CpGrid grid; const std::array cell_sizes = {1.0, 1.0, 1.0}; @@ -673,19 +721,64 @@ BOOST_AUTO_TEST_CASE(distributed_lgr) const std::vector> endIJK_vec = {{3,1,1}}; const std::vector lgr_name_vec = {"LGR1"}; // LGR1 element indices = 1 (rank 0), 2 (rank 2). Total 16 refined cells, 45 points (45-12 = 33 with new global id). - // LGR1 dim 2x1x1 (16 refined cells) (45 points - only 33 new points) + // LGR1 dim 4x2x2 (16 refined cells) (45 points - only 33 new points) grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + + // Check global id is not duplicated for points for each LGR + // LGR1 dim 4x2x2 -> 5x3x3 = 45 points + std::vector local_point_ids; + local_point_ids.reserve(45); // expected_point_ids in LGR1 + for (const auto& element : elements(grid.levelGridView(1))) { + for (int corner = 0; corner < 8; ++corner) + { + const auto& point = element.subEntity<3>(corner); + local_point_ids.push_back(grid.currentData()[1]->globalIdSet().id(point)); + } + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + // Coarse cells 1 in rank 0 and 2 in rank 2 share an I_FACE where (LGR1_dim[1]+1)*(LGR1_dim[2]+ 1), + // here (2 +1)*(2 +1) = 9 points lying on, 4 of them being the 4 corners of the coarse I_FACE shared + // by cell 1 and cell 2. Leaving us with 9 - 4 = 5 potential duplicated ids. + // + // Global cell ids of cells to be refined = {1, 2} + // cell 1 = { interior in P0, overlap in P1, overlap in P2, does not exist in P3} + // cell 2 = { overlap in P0, does not exist in P1, interior in P2, overlap in P3} + // Remark: the LGR is distributed in only two of the four processes, P0 and P2. + // - P0 sees the entire LGR. + // - P1 does NOT see cell 2, but sees cell 1 + // - P2 does NOT see the LGR at all. + // - P3 does NOT see cell 1, but sees cell 2 + // P3 creates: + // +5 duplicated ids on {I_FACE, false} cell 2 (seen in P3) [equivalent face: {I_FACE, true} of unseen-in-P3 cell 1] + // That means that the "unfortunate" expected point ids count is 45 (desired value) + 5 (duplicated laying on + // shared I_FACE) = 50. + BOOST_CHECK( static_cast(all_point_ids_set.size()) == 50); + + /** Current approach avoids duplicated point ids when + // 1. the LGR is distributed in P_{i_0}, ..., P_{i_n}, with n+1 < grid.comm().size(), + // AND + // 2. there is no coarse cell seen by a process P with P != P_{i_j}, j = 0, ..., n. + // Otherwise, there will be points with multiple ids.*/ + + // Check global id is not duplicated for points + std::vector localPointIds_vec; + localPointIds_vec.reserve(grid.currentData().back()->size(3)); + for (const auto& point : vertices(grid.leafGridView())) { + // Notice that all partition type points are pushed back. Selecting only interior points does not bring us to the expected value. + localPointIds_vec.push_back(grid.currentData().back()->globalIdSet().id(point)); + } + auto [allGlobalIds_points, displPointLeaf ] = Opm::allGatherv(localPointIds_vec, grid.comm()); + const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); + // Total global ids in leaf grid view for points: 80 + (45 - 12) + 5 (undesired duplicated ids on shared I_FACE) = 118 + BOOST_CHECK( allGlobalIds_points_set.size() == 118 ); } } - BOOST_AUTO_TEST_CASE(distributed_lgr_II) { - // Only for testing assignment of new global ids for refined entities (cells and point belonging to - // refined level grids). - // Create a grid Dune::CpGrid grid; const std::array cell_sizes = {1.0, 1.0, 1.0}; @@ -714,14 +807,44 @@ BOOST_AUTO_TEST_CASE(distributed_lgr_II) grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + + // Check global id is not duplicated for points for each LGR + // LGR1 dim 6x2x2 -> 7x3x3 = 63 points + std::vector local_point_ids; + local_point_ids.reserve(63); // expected_point_ids in LGR1 + for (const auto& point : vertices(grid.levelGridView(1))) { + local_point_ids.push_back(grid.currentData()[1]->globalIdSet().id(point)); + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + BOOST_CHECK( static_cast(all_point_ids_set.size()) == 63); + // Difference with previous test case: + // Global cell ids of cells to be refined = {8, 9, 10} + // cell 8 = { interior in P0, does not exist in P1, does not exist in P2, does not exist in P3} + // cell 9 = { interior in P0, does not exist in P1, overlap in P2, does not exist in P3} + // cell 10 = { overlap in P0, does not exist in P1, interior in P2, does not exist in P3} + // Remark: the LGR is distributed only in 2 processes which are the only two processes seeing these cells. + // - P0 sees the entire LGR. + // - P1 does NOT see the LGR at all. + // - P2 does NOT see cell 8, but NO OTHER process sees cell 8 since it's fully inteior of P0. + // - P3 does NOT see the LGR at all. + + // Check global id is not duplicated for points + std::vector localPointIds_vec; + localPointIds_vec.reserve(grid.currentData().back()->size(3)); + for (const auto& point : vertices(grid.leafGridView())) { + // Notice that all partition type points are pushed back. Selecting only interior points does not bring us to the expected value. + localPointIds_vec.push_back(grid.currentData().back()->globalIdSet().id(point)); + } + auto [allGlobalIds_points, displPointLeaf ] = Opm::allGatherv(localPointIds_vec, grid.comm()); + const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); + // Total global ids in leaf grid view for points: 80 + (63 - 16) = 127 + BOOST_CHECK( allGlobalIds_points_set.size() == 127 ); } } BOOST_AUTO_TEST_CASE(distributed_in_all_ranks_lgr) { - // Only for testing assignment of new global ids for refined entities (cells and point belonging to - // refined level grids). - // Create a grid Dune::CpGrid grid; const std::array cell_sizes = {1.0, 1.0, 1.0}; @@ -753,6 +876,61 @@ BOOST_AUTO_TEST_CASE(distributed_in_all_ranks_lgr) // 64 new refined cells. 5x5x5 = 125 points (only 98 = 125 - 3x3x3 parent corners new points - new global ids). grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); refinePatch_and_check(grid, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec); + + // Check global id is not duplicated for points for each LGR + // LGR1 dim 4x4x4 -> 5x5x5 = 125 points + std::vector local_point_ids; + local_point_ids.reserve(125); // expected_point_ids in LGR1 + for (const auto& element : elements(grid.levelGridView(1))) { + for (int corner = 0; corner < 8; ++corner) + { + const auto& point = element.subEntity<3>(corner); + local_point_ids.push_back(grid.currentData()[1]->globalIdSet().id(point)); + } + } + auto [all_point_ids, displPoint ] = Opm::allGatherv(local_point_ids, grid.comm()); + const std::set all_point_ids_set(all_point_ids.begin(), all_point_ids.end()); + + /** [Bug] Uniqueness of point global ids cannot be checked in general since current code sets overlap layer size equal to 1, + which in particular means that cells that share corners or edges (and not faces) with interior cells are not considered/ + seen by the process. Therefore, depending how the LGRs are distributed, there may be "multiple ids" for the same points.*/ + // Difference with previous test case: + // Global cell ids of cells to be refined = {1,2,5,6,13,14,17,18} + // cell 1 = { interior in P0, overlap in P1, overlap in P2, does not exist in P3} + // cell 2 = { overlap in P0, does not exist in P1, interior in P2, overlap in P3} + // cell 5 = { interior in P0, overlap in P1, overlap in P2, does not exist in P3} + // cell 6 = { overlap in P0, does not exist in P1, interior in P2, does not exist in P3} + // cell 13 = { overlap in P0, interior in P1, does not exist in P2, overlap in P3} + // cell 14 = { does not exist in P0, overlap in P1, overlap in P2, interior in P3} + // cell 17 = { overlap in P0, interior in P1, overlap in P2, does not exist in P3} + // cell 18 = { does not exist in P0, overlap in P1, interior in P2, overlap in P3} + // Remark: the LGR is distributed in ALL processes BUT + // - P0 does NOT see cell 18, but sees all the others. + // - P1 does NOT see cells 2 and 6, but sees the rest of them. + // - P2 does NOT see cell 13, but sees all the others. + // - P3 does NOT see cell 1, 5 and 17, but sees all the others. + // More details: + // P0 creates +1 duplicated id (middle point edge) + // P1 creates +2 duplicated ids (middle point of an edge in cell 2, cell 6 respectively). + // P2 creates +1 duplicated id (middle point edge) + // P3 creates: + // +5 duplicated ids on {I_FACE, false} cell 2 (seen in P3) [equivalent face: {I_FACE, true} of unseen-in-P3 cell 1] + // +1 duplicated id on edge unseen-in-P3 cell 5 + // +5 duplicated ids on {J_FACE, true} cell 13 (seen in P3) [equivalent face: {J_FACE, false} of unseen-in-P3 cell 17] + // +5 duplicated ids on {I_FACE, false} cell 18 (seen in P3) [equivalent face: {I_FACE, true} of unseen-in-P3 cell 17] + BOOST_CHECK( static_cast(all_point_ids_set.size()) == 145); // Desired value: 125; 20 (=1+2+1+5+1+5+5) duplicated ids. + + // Check global id is not duplicated for points + std::vector localPointIds_vec; + localPointIds_vec.reserve(grid.currentData().back()->size(3)); + for (const auto& point : vertices(grid.leafGridView())) { + // Notice that all partition type points are pushed back. Selecting only interior points does not bring us to the expected value. + localPointIds_vec.push_back(grid.currentData().back()->globalIdSet().id(point)); + } + auto [allGlobalIds_points, displPointLeaf ] = Opm::allGatherv(localPointIds_vec, grid.comm()); + const std::set allGlobalIds_points_set(allGlobalIds_points.begin(), allGlobalIds_points.end()); + // Total global ids in leaf grid view for points: 80 + (125 - 27) = 178 desired value; +20 duplicated ids. + BOOST_CHECK( allGlobalIds_points_set.size() == 198 ); } } diff --git a/tests/cpgrid/avoidNNCinLGRsCpGrid_test.cpp b/tests/cpgrid/avoidNNCinLGRsCpGrid_test.cpp index 240d2cb24..976985aa0 100644 --- a/tests/cpgrid/avoidNNCinLGRsCpGrid_test.cpp +++ b/tests/cpgrid/avoidNNCinLGRsCpGrid_test.cpp @@ -150,6 +150,18 @@ BOOST_AUTO_TEST_CASE(NNCAtSeveralLgrs) testCase(deckString, nnc, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec, true); } +BOOST_AUTO_TEST_CASE(LgrWithNNC_and_lgrsWithoutNNC) +{ + Opm::NNC nnc; + nnc.addNNC(0, 2, 1.0); // connect cell 0 and cell 2 (both belong to LGR1). LGR2 does not have NNCs. + const std::vector> cells_per_dim_vec = {{2,2,2}, {4,4,4}}; + const std::vector> startIJK_vec = {{0,0,0},{0,0,4}}; + const std::vector> endIJK_vec = {{1,1,3}, {1,1,5}}; + // LGR1 cell indices = {0,1,2}, LGR2 cell indices = {4}. + const std::vector lgr_name_vec = {"LGR1", "LGR2"}; + testCase(deckString, nnc, cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec, true); +} + BOOST_AUTO_TEST_CASE(NNCoutsideLgrs) { Opm::NNC nnc;