Skip to content

Commit

Permalink
GlobalCellLevel to Leaf Index defined and tested
Browse files Browse the repository at this point in the history
  • Loading branch information
aritorto committed Sep 26, 2024
1 parent c70338f commit f1f25d0
Show file tree
Hide file tree
Showing 8 changed files with 247 additions and 21 deletions.
1 change: 1 addition & 0 deletions CMakeLists_files.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -122,6 +122,7 @@ if(Boost_VERSION_STRING VERSION_GREATER 1.53)
tests/cpgrid/geometry_test.cpp
tests/cpgrid/grid_lgr_test.cpp
tests/cpgrid/inactiveCell_lgr_test.cpp
tests/cpgrid/lgr_cartesian_idx_test.cpp
tests/cpgrid/lookUpCellCentroid_cpgrid_test.cpp
tests/cpgrid/lookupdataCpGrid_test.cpp
tests/cpgrid/shifted_cart_test.cpp
Expand Down
5 changes: 5 additions & 0 deletions opm/grid/CpGrid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -968,6 +968,11 @@ namespace Dune
/// For nested refinement, we lookup the oldest ancestor, from level zero.
void computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell_leaf);

/// @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
/// on the leaf grid view.
/// Notice that cells that vanished and do not appear on the leaf grid view will not be considered.
void mapGlobalCellLevelToLeafIndexSet();

/// @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
20 changes: 14 additions & 6 deletions opm/grid/common/CartesianIndexMapper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -43,12 +43,6 @@ namespace Dune
return 0;
}

/** \brief return number of cells in the active level zero grid. Only relevant for CpGrid specialization. */
int compressedLevelZeroSize() const
{
return 0;
}

/** \brief return index of the cells in the logical Cartesian grid */
int cartesianIndex( const int /* compressedElementIndex */) const
{
Expand All @@ -60,11 +54,25 @@ namespace Dune
{
}

/// Only relevant for CpGrid specialization.
/** \brief return number of cells in the active level zero grid. Only relevant for CpGrid specialization. */
int compressedLevelZeroSize() const
{
return 0;
}

/** \brief return Cartesian coordinate, i.e. IJK, for a given cell. Only relevant for CpGrid specialization.*/
void cartesianCoordinateLevel(const int /* compressedElementIndexOnLevel */,
std::array<int,dimension>& /* coordsOnLevel */, int /*level*/) const
{
}

/** \brief return index of the cells in the logical Cartesian grid, for refined level grids. Only relevant for CpGrid specialization. */
int cartesianIndexLevel( const int /* compressedElementIndex */ , const int level) const
{
return 0;
}
/// Only relevant for CpGrid specialization. END
};

} // end namespace Opm
Expand Down
24 changes: 17 additions & 7 deletions opm/grid/cpgrid/CartesianIndexMapper.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,12 +50,6 @@ namespace Dune
return grid_.globalCell().size();
}

int compressedLevelZeroSize() const
{

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

int cartesianIndex( const int compressedElementIndex ) const
{
assert( compressedElementIndex >= 0 && compressedElementIndex < compressedSize() );
Expand All @@ -67,13 +61,29 @@ namespace Dune
grid_.getIJK( compressedElementIndex, coords );
}

/// Additional methods realted to LGRs
int compressedLevelZeroSize() const
{
return (*grid_.currentData()[0]).size(0);
}

int cartesianIndexLevel( const int compressedElementIndex, const int level) const
{
if ((level < 0) || (level > grid_.maxLevel())) {
throw std::invalid_argument("Invalid level.\n");
}
assert( compressedElementIndex >= 0 && compressedElementIndex < grid_.currentData()[level]->size(0) );
return grid_.currentData()[level]->globalCell()[compressedElementIndex];
}

void cartesianCoordinateLevel(const int compressedElementIndexOnLevel, std::array<int,dimension>& coordsOnLevel, int level) const
{
if ((level < 0) || (level > grid_.maxLevel())) {
throw std::invalid_argument("Invalid level.\n");
}
(*grid_.currentData()[level]).getIJK( compressedElementIndexOnLevel, coordsOnLevel);
grid_.currentData()[level]->getIJK( compressedElementIndexOnLevel, coordsOnLevel);
}
/// Additional methods realted to LGRs. END
};

} // end namespace Opm
Expand Down
10 changes: 10 additions & 0 deletions opm/grid/cpgrid/CpGrid.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -686,6 +686,14 @@ void CpGrid::computeGlobalCellLeafGridViewWithLgrs(std::vector<int>& global_cell
}
}

void CpGrid::mapGlobalCellLevelToLeafIndexSet()
{
for (const auto& element : elements(leafGridView())) {
const auto& global_cell_level = currentData()[element.level()]->globalCell()[element.getEquivLevelElem().index()];
(*currentData()[element.level()]).globalCellLevel_to_leafIdx_[global_cell_level] = element.index();
}
}

void CpGrid::getIJK(const int c, std::array<int,3>& ijk) const
{
current_view_data_->getIJK(c, ijk);
Expand Down Expand Up @@ -2003,6 +2011,8 @@ bool CpGrid::adapt(const std::vector<std::array<int,3>>& cells_per_dim_vec,
computeGlobalCellLeafGridViewWithLgrs(global_cell_leaf);
(*data[levels + preAdaptMaxLevel +1]).global_cell_.swap(global_cell_leaf);

mapGlobalCellLevelToLeafIndexSet();

updateCornerHistoryLevels(cornerInMarkedElemWithEquivRefinedCorner,
elemLgrAndElemLgrCorner_to_refinedLevelAndRefinedCorner,
adaptedCorner_to_elemLgrAndElemLgrCorner,
Expand Down
21 changes: 20 additions & 1 deletion opm/grid/cpgrid/CpGridData.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,11 @@ class CpGridData
ijk[2] = gc / logical_cartesian_size_[1];
}

const std::vector<int>& globalCell() const
{
return global_cell_;
}

/// @brief Determine if a finite amount of patches (of cells) are disjoint, namely, they do not share any corner nor face.
///
/// @param [in] startIJK_vec Vector of Cartesian triplet indices where each patch starts.
Expand Down Expand Up @@ -609,6 +614,16 @@ class CpGridData
return logical_cartesian_size_;
}

/// Get map from global_cell_ to leaf index set
///
/// Only relevant for CpGrid with LGRs. The maps of
/// all refined level grids is defined in CpGrid::adapt(/*..*/),
/// via CpGrid::mapGlobalCellLevelToLeafIndexSet().
const std::unordered_map<int,int>& getGlobalCellLevelToLeafIndexSet() const
{
return globalCellLevel_to_leafIdx_;
}

/// \brief Redistribute a global grid.
///
/// The whole grid must be available on all processors.
Expand Down Expand Up @@ -852,7 +867,11 @@ class CpGridData
/** Level-grid or Leaf-grid cell to parent cell and refined-cell-in-parent-cell index (number between zero and total amount
of children per parent (cells_per_dim[0]_*cells_per_dim_[1]*cells_per_dim_[2])). Entry is -1 when cell has no father. */
std::vector<int> cell_to_idxInParentCell_;

// Only relevant for refined level grids.
/** @brief Map from global_cell_lgr to leaf index set. Only relevant for refined level grids. */
std::unordered_map<int,int> globalCellLevel_to_leafIdx_;




/// \brief Object for collective communication operations.
Expand Down
21 changes: 14 additions & 7 deletions opm/grid/polyhedralgrid/cartesianindexmapper.hh
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@ namespace Dune

const std::array<int, dimension>& cartesianDimensions() const
{
return grid_.logicalCartesianSize();
return grid_.logicalCartesianSize();
}

int cartesianSize() const
Expand All @@ -44,12 +44,6 @@ namespace Dune
return grid_.size( 0 );
}

// Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
int compressedLevelZeroSize() const
{
return grid_.size( 0 );
}

int cartesianIndex( const int compressedElementIndex ) const
{
assert( compressedElementIndex >= 0 && compressedElementIndex < compressedSize() );
Expand All @@ -73,6 +67,19 @@ namespace Dune
coords[ 0 ] = gc ;
}

// Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
int compressedLevelZeroSize() const
{
return grid_.size( 0 );
}
// Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
int cartesianIndexLevel( const int compressedElementIndex, const int level ) const
{
if (level) {
throw std::invalid_argument("Invalid level.\n");
}
return cartesianIndex(compressedElementIndex);
}
// Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant.
void cartesianCoordinateLevel(const int compressedElementIndexOnLevel, std::array<int,dimension>& coordsOnLevel, int level) const
{
Expand Down
166 changes: 166 additions & 0 deletions tests/cpgrid/lgr_cartesian_idx_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
//===========================================================================
//
// File: lgr_cartesian_idx_test.cpp
//
// Created: Sep 26 2024 09:15
//
// Author(s): Antonella Ritorto <[email protected]>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"

#define BOOST_TEST_MODULE LGRTests
#include <boost/test/unit_test.hpp>
#include <boost/version.hpp>
#if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 71
#include <boost/test/floating_point_comparison.hpp>
#else
#include <boost/test/tools/floating_point_comparison.hpp>
#endif
#include <opm/grid/CpGrid.hpp>
#include <opm/grid/cpgrid/CpGridData.hpp>
#include <opm/grid/cpgrid/DefaultGeometryPolicy.hpp>
#include <opm/grid/cpgrid/Entity.hpp>
#include <opm/grid/cpgrid/EntityRep.hpp>
#include <opm/grid/cpgrid/Geometry.hpp>
#include <opm/grid/LookUpData.hh>

#include <dune/grid/common/mcmgmapper.hh>

#include <cassert>
#include <sstream>
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <map>

struct Fixture
{
Fixture()
{
int m_argc = boost::unit_test::framework::master_test_suite().argc;
char** m_argv = boost::unit_test::framework::master_test_suite().argv;
Dune::MPIHelper::instance(m_argc, m_argv);
Opm::OpmLog::setupSimpleDefaultLogging();
}

static int rank()
{
int m_argc = boost::unit_test::framework::master_test_suite().argc;
char** m_argv = boost::unit_test::framework::master_test_suite().argv;
return Dune::MPIHelper::instance(m_argc, m_argv).rank();
}
};

BOOST_GLOBAL_FIXTURE(Fixture);

void checkGlobalCellLgr(Dune::CpGrid& grid)
{
const Dune::CartesianIndexMapper<Dune::CpGrid> mapper{grid};

for (const auto& element : elements(grid.leafGridView()))
{
// How to get the Cartesian Index of a cell on the leaf grid view.
// 1. The "default" Cartesian Index of a cell on the leaf (which is a mixed grid with coarse and refined cells) is equal to the Cartesian Index of:
// (a) the equivalent cell in level zero (same center and volume, cell not involved in any refinement).
// (b) the parent cell in level zero when the leaf cell is a refined one (it has a equivalent cell that belongs to a refined level grid).
// This default Cartesian Index for leaf cells coincide with the globalCell values.
const auto& cartesian_idx_from_leaf_elem = mapper.cartesianIndex( element.index() );
const auto& global_cell_idx_leaf = grid.globalCell()[element.index()]; // parent cell index when the cell is a refined one.
BOOST_CHECK_EQUAL(cartesian_idx_from_leaf_elem, global_cell_idx_leaf);
// global_ijk represents the ijk values of the equivalent cell on the level zero, or parent cell if the leaf cell is a refined one.
// Notice that all the refined cells of a same parent cell will get the same global_cell_ value and the same global_ijk (since they
// inherit the parent cell value).
std::array<int,3> global_ijk = {0,0,0};
mapper.cartesianCoordinate( element.index(), global_ijk);

// How to get the Level Cartesian Index of a cell on the leaf grid view.
// Each LGR can be seen as a Cartesian Grid itself, with its own (local) logical_cartesian_size and its own (local) Cartesian indices.
// Given a leaf cell, via the CartesianIndexMapper and its method cartesianIndexLevel(...), we get the local-Cartesian-index.
const auto& cartesian_idx_from_level_elem = mapper.cartesianIndexLevel( element.getEquivLevelElem().index(), element.level() );
const auto& global_cell_idx_level = grid.currentData()[element.level()]->globalCell()[element.getEquivLevelElem().index()];
BOOST_CHECK_EQUAL(cartesian_idx_from_level_elem, global_cell_idx_level);
// local_ijk represents the ijk values of the equivalent cell on the level its was born.
std::array<int,3> local_ijk = {0,0,0};
mapper.cartesianCoordinateLevel( element.getEquivLevelElem().index(), local_ijk, element.level() );

// For leaf cells that were not involved in any refinement, global_ijk and local_ijk must coincide.
if(element.level()==0)
{
BOOST_CHECK_EQUAL( global_ijk[0], local_ijk[0]);
BOOST_CHECK_EQUAL( global_ijk[1], local_ijk[1]);
BOOST_CHECK_EQUAL( global_ijk[2], local_ijk[2]);
}
}

for (int level = 0; level < grid.maxLevel(); ++level)
{
const auto& globalCellLevel_to_leafIdx = grid.currentData()[level]->getGlobalCellLevelToLeafIndexSet();
for (const auto& element : elements(grid.levelGridView(level)))
{
if(element.isLeaf()) {
const auto& global_cell_level = grid.currentData()[element.level()]->globalCell()[element.index()];
globalCellLevel_to_leafIdx.at(global_cell_level);
}
else {
BOOST_CHECK_THROW( globalCellLevel_to_leafIdx.at(element.index()), std::out_of_range);
}
}
}
}

BOOST_AUTO_TEST_CASE(refine_one_cell)
{
// Create a grid
Dune::CpGrid grid;
const std::array<double, 3> cell_sizes = {1.0, 1.0, 1.0};
const std::array<int, 3> grid_dim = {4,3,3};
grid.createCartesian(grid_dim, cell_sizes);

const std::array<int, 3> cells_per_dim = {2,2,2};
const std::array<int, 3> startIJK = {1,0,1};
const std::array<int, 3> endIJK = {2,1,2};// Single Cell! (with index 13 in level zero grid), LGR dimensions 2x2x2
const std::string lgr_name = {"LGR1"};
grid.addLgrsUpdateLeafView({cells_per_dim}, {startIJK}, {endIJK}, {lgr_name});

checkGlobalCellLgr(grid);
}

BOOST_AUTO_TEST_CASE(three_lgrs)
{
// Create a grid
Dune::CpGrid grid;
const std::array<double, 3> cell_sizes = {1.0, 1.0, 1.0};
const std::array<int, 3> grid_dim = {4,3,3};
grid.createCartesian(grid_dim, cell_sizes);

const std::vector<std::array<int,3>> cells_per_dim_vec = {{2,2,2}, {3,3,3}, {4,4,4}};
const std::vector<std::array<int,3>> startIJK_vec = {{0,0,0}, {0,0,2}, {3,2,2}};
const std::vector<std::array<int,3>> endIJK_vec = {{2,1,1}, {1,1,3}, {4,3,3}};
const std::vector<std::string> lgr_name_vec = {"LGR1", "LGR2", "LGR3"};
grid.addLgrsUpdateLeafView(cells_per_dim_vec, startIJK_vec, endIJK_vec, lgr_name_vec);

checkGlobalCellLgr(grid);
}

0 comments on commit f1f25d0

Please sign in to comment.