Skip to content

Commit

Permalink
Modify global_cell_ for lgrs
Browse files Browse the repository at this point in the history
  • Loading branch information
@aritorto
aritorto committed Sep 25, 2024
1 parent fc42f67 commit 694547c
Show file tree
Hide file tree
Showing 10 changed files with 478 additions and 388 deletions.
26 changes: 22 additions & 4 deletions opm/grid/CpGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -309,7 +309,7 @@ namespace Dune
/// those dealing with permeability fields from the input deck
/// from whence the current CpGrid was constructed.
const std::vector<int>& globalCell() const;

/// @brief Returns either data_ or distributed_data_(if non empty).
const std::vector<std::shared_ptr<Dune::cpgrid::CpGridData>>& currentData() const;

Expand Down Expand Up @@ -889,10 +889,9 @@ namespace Dune
const std::vector<std::array<int,3>>& cells_per_dim_vec,
const int& preAdaptMaxLevel) const;

/// @brief Define the cells, cell_to_point_, global_cell_, cell_to_face_, face_to_cell_, for the leaf grid view (or adapted grid).
/// @brief Define the cells, cell_to_point_, cell_to_face_, face_to_cell_, for the leaf grid view (or adapted grid).
void populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
std::vector<std::array<int,8>>& adapted_cell_to_point,
std::vector<int>& adapted_global_cell,
const int& cell_count,
cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
Expand Down Expand Up @@ -922,7 +921,6 @@ namespace Dune
/* Leaf grid View Cells argumemts */
Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
std::vector<std::array<int,8>>& adapted_cell_to_point,
std::vector<int>& adapted_global_cell,
const int& cell_count,
cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
Expand Down Expand Up @@ -950,6 +948,26 @@ namespace Dune
const int& preAdaptMaxLevel,
const int& newLevels);


/// @brief For refined level grids created based on startIJK and endIJK values, compute the "local ijk/Cartesian index" within the LGR.
///
/// It's confusing that this "localIJK" is stored in CpGridData member global_cell_. Potential explanation: level zero grid is also called
/// "GLOBAL" grid.
/// Example: a level zero grid with dimension 4x3x3, an LGR with startIJK = {1,2,2}, endIJK = {3,3,3}, and cells_per_dim = {2,2,2}.
/// Then the dimension of the LGR is (3-1)*2 x (3-2)*2 x (3-2)* 2 = 4x2x2 = 16. Therefore the global_cell_lgr minimim value should be 0,
/// the maximum should be 15.
/// To invoke this method, each refined level grid must have 1. logical_cartesian_size_, 2. cell_to_idxInParentCell_, and 3. cells_per_dim_
/// already populated.
///
/// @param [in] level Grid index where LGR is stored
/// @param [out] global_cell_lgr
void computeGlobalCellLgr(const int& level, const std::array<int,3>& startIJK, std::vector<int>& global_cell_lgr);

/// @brief For a leaf grid with with LGRs, we assign the global_cell_ values of either the parent cell or the equivalent cell from
/// level zero.
/// For nested refinement, we lookup the oldest ancestor, from level zero.
void computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf);

/// @brief Get the ijk index of a refined corner, given its corner index of a single-cell-refinement.
///
/// Given a single-cell, we refine it in {nx, ny, nz} refined children cells (per direction). Then, this single-cell-refinement
Expand Down
2 changes: 1 addition & 1 deletion opm/grid/cpgrid/CartesianIndexMapper.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ namespace Dune

int compressedLevelZeroSize() const
{

return (*grid_.currentData()[0]).size(0);
}

Expand Down
86 changes: 75 additions & 11 deletions opm/grid/cpgrid/CpGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -630,7 +630,60 @@ const std::vector<int>& CpGrid::globalCell() const
{
// Temporary. For a grid with LGRs, we set the globalCell() of the as the one for level 0.
// Goal: CartesianIndexMapper well-defined for CpGrid LeafView with LGRs.
return currentData().back() -> global_cell_;
return currentData().back() /*current_view_data_ */-> global_cell_;
}

void CpGrid::computeGlobalCellLgr(const int& level, const std::array<int,3>& startIJK, std::vector<int>& global_cell_lgr)
{
assert(level);
for (const auto& element : elements(levelGridView(level))) {
// Element belogns to an LGR, therefore has a father. Get IJK of the father in the level grid the father was born.
// For CARFIN, parent cells belong to level 0.
std::array<int,3> parentIJK = {0,0,0};
currentData()[element.father().level()]->getIJK(element.father().index(), parentIJK);
// Each parent cell has been refined in cells_per_dim[0]*cells_per_dim[1]*cells_per_dim[2] child cells.
// element has certain 'position' inside its parent cell that can be described with 'IJK' indices, let's denote them by ijk,
// where 0<= i < cells_per_dim[0], 0<= j < cells_per_dim[1], 0<= k < cells_per_dim[2].
const auto& cells_per_dim = currentData()[level]->cells_per_dim_;
//
// Refined cell (here 'element') has "index in parent cell": k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0] + i
// and it's stored in cell_to_idxInParentCell_.
auto idx_in_parent_cell = currentData()[level]-> cell_to_idxInParentCell_[element.index()]; // deep copy
// Find ijk.
std::array<int,3> childIJK = {0,0,0};

childIJK[0] = idx_in_parent_cell % cells_per_dim[0];
idx_in_parent_cell -= childIJK[0]; // k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0]
idx_in_parent_cell /= cells_per_dim[0]; // k*cells_per_dim[1] + j
childIJK[1] = idx_in_parent_cell % cells_per_dim[1];
idx_in_parent_cell -= childIJK[1]; // k*cells_per_dim[1]
childIJK[2] = idx_in_parent_cell /cells_per_dim[1];
// The corresponding lgrIJK can be computed as follows:
const std::array<int,3>& lgrIJK = { ( (parentIJK[0] - startIJK[0])*cells_per_dim[0] ) + childIJK[0], // Shift parent index according to the startIJK of the LGR.
( (parentIJK[1] - startIJK[1])*cells_per_dim[1] ) + childIJK[1],
( (parentIJK[2] - startIJK[2])*cells_per_dim[2] ) + childIJK[2] };
// Dimensions of the "patch of cells" formed when providing startIJK and endIJK for an LGR
const auto& lgr_logical_cartesian_size = currentData()[level]->logical_cartesian_size_;
global_cell_lgr[element.index()] = (lgrIJK[2]*lgr_logical_cartesian_size[0]*lgr_logical_cartesian_size[1]) + (lgrIJK[1]*lgr_logical_cartesian_size[0]) + lgrIJK[0];
}
}

void CpGrid::computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf)
{
for (const auto& element: elements(leafGridView()))
{
// When refine via CpGrid::addLgrsUpdateGridView(/*...*/), level-grid to lookup global_cell_ is equal to level-zero-grid
// In the context of allowed nested refinement, we lookup for the oldest ancestor, also belonging to level-zero-grid.
auto ancestor = element.getOrigin();
int origin_in_level_zero = ancestor.index();
while (ancestor.level()>0){
ancestor = ancestor.getOrigin();
origin_in_level_zero = ancestor.getOrigin().index();
}
assert(ancestor.level()==0);
assert(origin_in_level_zero < currentData().front()->size(0));
global_cell_leaf[element.index()] = currentData().front()-> global_cell_[origin_in_level_zero];
}
}

void CpGrid::getIJK(const int c, std::array<int,3>& ijk) const
Expand Down Expand Up @@ -1809,7 +1862,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
cornerInMarkedElemWithEquivRefinedCorner,
cells_per_dim_vec);

std::vector<int> adapted_global_cell(cell_count, 0);
updateLeafGridViewGeometries( /* Leaf grid View Corners arguments */
adapted_corners,
corner_count,
Expand All @@ -1822,7 +1874,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
/* Leaf grid View Cells argumemts */
adapted_cells,
adapted_cell_to_point,
adapted_global_cell,
cell_count,
adapted_cell_to_face,
adapted_face_to_cell,
Expand Down Expand Up @@ -1895,7 +1946,6 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
(*data[refinedLevelGridIdx]).child_to_parent_cells_ = refined_child_to_parent_cells_vec[level];
(*data[refinedLevelGridIdx]).cell_to_idxInParentCell_ = refined_cell_to_idxInParentCell_vec[level];
(*data[refinedLevelGridIdx]).level_to_leaf_cells_ = refined_level_to_leaf_cells_vec[level];
(*data[refinedLevelGridIdx]).global_cell_.swap(refined_global_cell_vec[level]);
(*data[refinedLevelGridIdx]).index_set_ = std::make_unique<cpgrid::IndexSet>(data[refinedLevelGridIdx]->size(0),
data[refinedLevelGridIdx]->size(3));
// Determine the amount of cells per direction, per parent cell, of the corresponding LGR.
Expand All @@ -1912,6 +1962,7 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
}
else {
(*data[refinedLevelGridIdx]).logical_cartesian_size_ = (*data[0]).logical_cartesian_size_;
(*data[refinedLevelGridIdx]).global_cell_.swap(refined_global_cell_vec[level]);
}
// One alternative definition for logical_cartesian_size_ in the case where the marked elements for refinement do not form a block of cells,
// therefore, are not associated with the keyword CARFIN, is to imagine that we put all the marked elements one next to the other, along
Expand All @@ -1927,14 +1978,31 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
(*data[levels + preAdaptMaxLevel +1]).child_to_parent_cells_ = adapted_child_to_parent_cells;
(*data[levels + preAdaptMaxLevel +1]).cell_to_idxInParentCell_ = adapted_cell_to_idxInParentCell;
(*data[levels + preAdaptMaxLevel +1]).leaf_to_level_cells_ = leaf_to_level_cells;
(*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(adapted_global_cell);
(*data[levels + preAdaptMaxLevel +1]).index_set_ = std::make_unique<cpgrid::IndexSet>(data[levels + preAdaptMaxLevel +1]->size(0),
data[levels + preAdaptMaxLevel +1]->size(3));
(*data[levels + preAdaptMaxLevel +1]).logical_cartesian_size_ = (*data[0]).logical_cartesian_size_;

// Update the leaf grid view
current_view_data_ = data.back().get();

// When the refinement is determined by startIJK and endIJK values, the LGR has a (local) Cartesian size.
// Therefore, each refined cell belonging to the LGR can be associated with a (local) IJK and its (local) Cartesian index.
// If the LGR has NXxNYxNZ dimension, then the Cartesian indices take values
// k*NN*NY + j*NX + i, where i<NX, j<Ny, k<NZ.
// This index is stored in <refined-level-grid>.global_cell_[ refined cell index (~element.index()) ] = k*NN*NY + j*NX + i.
if (isCARFIN) {
for (int level = 0; level < levels; ++level) {
const int refinedLevelGridIdx = level + preAdaptMaxLevel +1;
std::vector<int> global_cell_lgr(data[refinedLevelGridIdx]->size(0));
computeGlobalCellLgr(refinedLevelGridIdx, startIJK_vec[level], global_cell_lgr);
(*data[refinedLevelGridIdx]).global_cell_.swap(global_cell_lgr);
}
}

std::vector<int> global_cell_leaf( data[levels + preAdaptMaxLevel +1]->size(0));
computeGlobalCellLeafGridViewWithLgrs(global_cell_leaf);
(*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(global_cell_leaf);

updateCornerHistoryLevels(cornerInMarkedElemWithEquivRefinedCorner,
elemLgrAndElemLgrCorner_to_refinedLevelAndRefinedCorner,
adaptedCorner_to_elemLgrAndElemLgrCorner,
Expand Down Expand Up @@ -2112,6 +2180,7 @@ void CpGrid::addLgrsUpdateLeafView(const std::vector<std::array<int,3>>& cells_p
auto globalActiveLgrs = comm().sum(non_empty_lgrs);
if(globalActiveLgrs == 0) {
Opm::OpmLog::warning("All the LGRs contain only inactive cells.\n");
// OPM_THROW(std::logic_error, "All the LGRs contain only inactive cells.\n");
}

preAdapt();
Expand Down Expand Up @@ -2529,7 +2598,7 @@ CpGrid::defineChildToParentAndIdxInParentCell(const std::map<std::array<int,2>,s
const auto& element = Dune::cpgrid::Entity<0>(*current_view_data_, elemLgr, true); // elemLgr == parent cell index in starting grid.
refined_child_to_parent_cells_vec[shiftedLevel][cell] = {element.level(), element.getEquivLevelElem().index()};
refined_cell_to_idxInParentCell_vec[shiftedLevel][cell] = elemLgrCell;
}
}
}

return std::make_tuple<std::vector<std::vector<std::array<int,2>>>, std::vector<std::vector<int>>,
Expand Down Expand Up @@ -3237,7 +3306,6 @@ void CpGrid::populateRefinedFaces(std::vector<Dune::cpgrid::EntityVariableBase<c

void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
std::vector<std::array<int,8>>& adapted_cell_to_point,
std::vector<int>& adapted_global_cell,
const int& cell_count,
cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
Expand Down Expand Up @@ -3265,8 +3333,6 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
const auto& [elemLgr, elemLgrCell] = adaptedCell_to_elemLgrAndElemLgrCell.at(cell);
const auto& elemLgrCellEntity = Dune::cpgrid::EntityRep<0>(elemLgrCell, true);

adapted_global_cell[cell] = current_view_data_->global_cell_[(elemLgr == -1) ? elemLgrCell : elemLgr];

// Auxiliary cell_to_face
std::vector<cpgrid::EntityRep<1>> aux_cell_to_face;

Expand Down Expand Up @@ -3574,7 +3640,6 @@ void CpGrid::updateLeafGridViewGeometries( /* Leaf grid View Corners arguments *
/* Leaf grid View Cells argumemts */
Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<3,3>>& adapted_cells,
std::vector<std::array<int,8>>& adapted_cell_to_point,
std::vector<int>& adapted_global_cell,
const int& cell_count,
cpgrid::OrientedEntityTable<0,1>& adapted_cell_to_face,
cpgrid::OrientedEntityTable<1,0>& adapted_face_to_cell,
Expand Down Expand Up @@ -3617,7 +3682,6 @@ void CpGrid::updateLeafGridViewGeometries( /* Leaf grid View Corners arguments *
// --- Adapted cells ---
populateLeafGridCells(adapted_cells,
adapted_cell_to_point,
adapted_global_cell,
cell_count,
adapted_cell_to_face,
adapted_face_to_cell,
Expand Down
15 changes: 8 additions & 7 deletions opm/grid/cpgrid/CpGridData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1710,16 +1710,17 @@ void CpGridData::computeCommunicationInterfaces([[maybe_unused]] int noExistingP
#endif
}

std::array<Dune::FieldVector<double,3>,8> CpGridData::getReferenceRefinedCorners(int idxInParentCell, const std::array<int,3>& cells_per_dim) const
std::array<Dune::FieldVector<double,3>,8> CpGridData::getReferenceRefinedCorners(const int& idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const
{
// Refined cells in parent cell: k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0] + i
std::array<int,3> ijk = {0,0,0};
ijk[0] = idxInParentCell % cells_per_dim[0];
idxInParentCell -= ijk[0]; // k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0]
idxInParentCell /= cells_per_dim[0]; // k*cells_per_dim[1] + j
ijk[1] = idxInParentCell % cells_per_dim[1];
idxInParentCell -= ijk[1]; // k*cells_per_dim[1]
ijk[2] = idxInParentCell /cells_per_dim[1];
int idx_copy = idx_in_parent_cell;
ijk[0] = idx_copy % cells_per_dim[0];
idx_copy -= ijk[0]; // k*cells_per_dim[0]*cells_per_dim[1] + j*cells_per_dim[0]
idx_copy /= cells_per_dim[0]; // k*cells_per_dim[1] + j
ijk[1] = idx_copy % cells_per_dim[1];
idx_copy -= ijk[1]; // k*cells_per_dim[1]
ijk[2] = idx_copy /cells_per_dim[1];

std::array<Dune::FieldVector<double,3>,8> corners_in_parent_reference_elem = { // corner '0'
{{ double(ijk[0])/cells_per_dim[0], double(ijk[1])/cells_per_dim[1], double(ijk[2])/cells_per_dim[2] },
Expand Down
2 changes: 1 addition & 1 deletion opm/grid/cpgrid/CpGridData.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -381,7 +381,7 @@ class CpGridData
void postAdapt();

private:
std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(int idxInParentCell, const std::array<int,3>& cells_per_dim) const;
std::array<Dune::FieldVector<double,3>,8> getReferenceRefinedCorners(const int& idx_in_parent_cell, const std::array<int,3>& cells_per_dim) const;

/// @brief Compute amount of cells in each direction of a patch of cells. (Cartesian grid required).
///
Expand Down
14 changes: 6 additions & 8 deletions opm/grid/cpgrid/Entity.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -529,11 +529,10 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
if (!(this->hasFather())){
OPM_THROW(std::logic_error, "Entity has no father.");
}
if (pgrid_ -> cell_to_idxInParentCell_[this->index()] !=-1) {
int idxInParentCell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
assert(idxInParentCell>-1);
const int& idx_in_parent_cell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
if (idx_in_parent_cell !=-1) {
const auto& cells_per_dim = (*(pgrid_ -> level_data_ptr_))[this->level()] -> cells_per_dim_;
const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idxInParentCell, cells_per_dim);
const auto& auxArr = pgrid_ -> getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
FieldVector<double, 3> corners_in_father_reference_elem_temp[8] =
{ auxArr[0], auxArr[1], auxArr[2], auxArr[3], auxArr[4], auxArr[5], auxArr[6], auxArr[7]};
auto in_father_reference_elem_corners = std::make_shared<EntityVariable<cpgrid::Geometry<0, 3>, 3>>();
Expand Down Expand Up @@ -614,10 +613,9 @@ Dune::cpgrid::Entity<0> Dune::cpgrid::Entity<codim>::getEquivLevelElem() const
template<int codim>
int Dune::cpgrid::Entity<codim>::getLevelCartesianIdx() const
{
const auto entityLevel = this -> level();
const auto level = (*(pgrid_ -> level_data_ptr_))[entityLevel].get();
const auto& elemInLevel = this->getLevelElem(); // throws when the entity does not belong to the leaf grid view.
return level -> global_cell_[elemInLevel.index()];
const auto& level_data = (*(pgrid_ -> level_data_ptr_))[level()].get();
// getLevelElem() throws when the entity does not belong to the leaf grid view.
return level_data -> global_cell_[getLevelElem().index()];
}

} // namespace cpgrid
Expand Down
Loading

0 comments on commit 694547c

Please sign in to comment.