Skip to content

Supporting Infrastructure for Geo-Mechanical Simulations #942

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
bska wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Supporting Infrastructure for Geo-Mechanical Simulations #942

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
95ad6db
-- fixed out of bound nnc
hnil May 20, 2025
e395612
Changed point communicator to be used for all points. Seems to fix mi…
hnil Aug 12, 2025
0342bde
Commit on maybe working grid lableling and adding of boundary in for …
hnil Aug 19, 2025
c16c5b3
small change
hnil Aug 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions opm/grid/CpGrid.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1221,6 +1221,9 @@ namespace Dune
/// of the vertex.
int faceVertex(int face, int local_index) const;


/// \brief Get maps from vertexToCell out[0] and cellToVertex out[1];
std::array< std::vector< std::set<int> >, 2 > vertexCell() const;
/// \brief Get vertical position of cell center ("zcorn" average).
/// \brief cell_index The index of the specific cell.
double cellCenterDepth(int cell_index) const;
Expand Down
36 changes: 27 additions & 9 deletions opm/grid/MinpvProcessor.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,7 @@ MinpvProcessor::process(const std::vector<double>& thickness,
bool c_active = actnum.empty() || actnum[c];
bool c_thin = (thickness[c] <= z_tolerance);
bool c_thin_inactive = !c_active && c_thin;
bool c_low_pv_active = pv[c] < minpvv[c] && c_active;
bool c_low_pv_active = (pv[c] < minpvv[c] && c_active) || (c_thin && c_active);

if (c_low_pv_active || c_thin_inactive) {
std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
Expand Down Expand Up @@ -168,7 +168,7 @@ MinpvProcessor::process(const std::vector<double>& thickness,
bool active = actnum.empty() || actnum[c_below];
bool thin = (thickness[c_below] <= z_tolerance);
bool thin_inactive = !active && thin;
bool low_pv_active = pv[c_below] < minpvv[c_below] && active;
bool low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active);


while ( (thin_inactive || low_pv_active) && kk_iter < dims_[2] )
Expand Down Expand Up @@ -225,13 +225,26 @@ MinpvProcessor::process(const std::vector<double>& thickness,
active = actnum.empty() || actnum[c_below];
thin = (thickness[c_below] <= z_tolerance);
thin_inactive = (!actnum.empty() && !actnum[c_below]) && thin;
low_pv_active = pv[c_below] < minpvv[c_below] && active;
low_pv_active = ((pv[c_below] < minpvv[c_below]) && active) || (thin && active);
}

// create nnc if false or merge the cells if true
if (mergeMinPVCells && c_low_pv_active) {
if (mergeMinPVCells){// && c_low_pv_active) {
// try to make a topological connected grid
// in Flow this is currently called only for edge_conformal grids
// however zcorn inbetween is not modified to make zcorn sorted
// Set lower k coordinates of cell below to upper cells's coordinates.
// i.e fill the void using the cell below
if (kk==0 || kk_iter == dims_[2]) {
kk = kk_iter;
continue;
}
//bottom cell not active, hence no nnc is created
if (!actnum.empty() && !actnum[c_below]) {
kk = kk_iter;
continue;
}

std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk_iter, zcorn);
for (int count = 0; count < 4; ++count) {
cz_below[count] = cz[count];
Expand Down Expand Up @@ -260,15 +273,15 @@ MinpvProcessor::process(const std::vector<double>& thickness,
auto above_active = actnum.empty() || actnum[c_above];
auto above_inactive = !actnum.empty() && !actnum[c_above];
auto above_thin = thickness[c_above] < z_tolerance;
auto above_small_pv = pv[c_above] < minpvv[c_above];
auto above_small_pv = (pv[c_above] < minpvv[c_above]) || above_thin;

if ((above_inactive && above_thin) || (above_active && above_small_pv
&& (!pinchNOGAP || above_thin) ) ) {
for (k_above = kk - 2; k_above > 0; --k_above) {
c_above = ii + dims_[0] * (jj + dims_[1] * (k_above));
above_active = actnum.empty() || actnum[c_above];
above_inactive = !actnum.empty() && !actnum[c_above];
auto above_significant_pv = pv[c_above] > minpvv[c_above];
auto above_significant_pv = !((pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance));
auto above_broad = thickness[c_above] > z_tolerance;

// \todo if condition seems wrong and should be the negation of above?
Expand All @@ -293,25 +306,30 @@ MinpvProcessor::process(const std::vector<double>& thickness,
option4ALLZero = option4ALLZero || (!permz.empty() && permz[c_above] == 0.0) || multz(c_above) == 0.0;
nnc_allowed = nnc_allowed && (computeGap(cz_above, cz_below) < max_gap) && (!pinchOption4ALL || !option4ALLZero) ;

//bool
above_small_pv = (pv[c_above] < minpvv[c_above]) || (thickness[c_above] < z_tolerance);
bool below_small_pv = (pv[c_below] < minpvv[c_below]) || (thickness[c_below] < z_tolerance);
if ( nnc_allowed &&
(actnum.empty() || (actnum[c_above] && actnum[c_below])) &&
pv[c_above] > minpvv[c_above] && pv[c_below] > minpvv[c_below]) {
!(above_small_pv) && !(below_small_pv) ){
result.add_nnc(c_above, c_below);
}
kk = kk_iter;
}
}
else
{
if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && pv[c] > minpvv[c] &&
if (kk < dims_[2] - 1 && (actnum.empty() || actnum[c]) && !((pv[c] < minpvv[c]) || (thickness[c] < z_tolerance)) &&
multz(c) != 0.0)
{
// Check whether there is a gap to the neighbor below whose thickness is less
// than MAX_GAP. In that case we need to create an NNC if there is a gap between the two cells.
int kk_below = kk + 1;
int c_below = ii + dims_[0] * (jj + dims_[1] * kk_below);

if ((actnum.empty() || actnum[c_below]) && pv[c_below] > minpvv[c_below])
if ( (actnum.empty() || actnum[c_below])
&&
!((pv[c_below] < minpvv[c_below]) || (thickness[c_below] < z_tolerance) ) )
{
// Check MAX_GAP threshold
std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
Expand Down
58 changes: 39 additions & 19 deletions opm/grid/common/GridPartitioning.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -260,6 +260,7 @@ void addOverlapCornerCell(const CpGrid& grid,
std::vector<std::set<int> >& cell_overlap,
int level)
{
assert(false);
bool validLevel = (level>-1) && (level <= grid.maxLevel());
const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet();
int my_index = ix.index(from);
Expand Down Expand Up @@ -291,6 +292,7 @@ void addOverlapCornerCell(const CpGrid& grid,
std::vector<std::tuple<int,int,char>>& exportList,
int level)
{
assert(false);
// Add corner cells to the overlap layer. Example of a subdomain of a 4x4 grid
// with and without corner cells in the overlap is given below. Note that the
// corner cell is not needed for cell centered finite volume schemes.
Expand Down Expand Up @@ -334,6 +336,7 @@ void addOverlapLayer(const CpGrid& grid,
int recursion_deps,
int level)
{
assert(false);
bool validLevel = (level>-1) && (level <= grid.maxLevel());
const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet();

Expand Down Expand Up @@ -390,6 +393,7 @@ void addOverlapLayer(const CpGrid& grid,
bool all,
int level)
{
assert(false);
cell_overlap.resize(cell_part.size());
bool validLevel = (level>-1) && (level <= grid.maxLevel());
const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet();
Expand Down Expand Up @@ -420,6 +424,7 @@ void addOverlapLayer(const CpGrid& grid,
const std::vector<int>& cell_part,
std::vector<std::tuple<int,int,char>>& exportList,
bool addCornerCells,
const std::array<std::vector< std::set<int> >,2>& vertex_cell_maps,
int recursion_deps,
int level)
{
Expand Down Expand Up @@ -447,29 +452,43 @@ void addOverlapLayer(const CpGrid& grid,
cell_part,
exportList,
addCornerCells,
vertex_cell_maps,
recursion_deps-1,
level);

}
else if (addCornerCells) {
// Add cells to the overlap that just share a corner with e.
auto iit2 = validLevel? iit->outside().ilevelbegin() : iit->outside().ileafbegin();
const auto& endIit2 = validLevel? iit->outside().ilevelend() : iit->outside().ileafend();

for (; iit2 != endIit2; ++iit2) {
if ( iit2->neighbor() )
{
int nb_index2 = ix.index(iit2->outside());
if( cell_part[nb_index2]!=owner ) {
addOverlapCornerCell(grid,
owner,
e,
iit2->outside(),
cell_part,
exportList,
level);
auto vertices = vertex_cell_maps[1][index]; //could probably used the subindices but do not trust them here
for(auto vertex: vertices){
std::set<int> cells = vertex_cell_maps[0][vertex];
for(auto& cell: cells){
if(cell_part[cell] != owner){
// cell is the neighbor to owner
exportList.emplace_back(cell, owner, AttributeSet::copy);
// do this since it is allways done probably an ERROR for
exportList.emplace_back(index, cell_part[cell], AttributeSet::copy);
}
}
}
// Add cells to the overlap that just share a corner with e.
// auto iit2 = validLevel? iit->outside().ilevelbegin() : iit->outside().ileafbegin();
// const auto& endIit2 = validLevel? iit->outside().ilevelend() : iit->outside().ileafend();

// for (; iit2 != endIit2; ++iit2) {
// if ( iit2->neighbor() )
// {
// int nb_index2 = ix.index(iit2->outside());
// if( cell_part[nb_index2]!=owner ) {
// addOverlapCornerCell(grid,
// owner,
// e,
// iit2->outside(),
// cell_part,
// exportList,
// level);
// }
// }
// }
}
}
}
Expand All @@ -487,6 +506,7 @@ void addOverlapLayerNoZeroTrans(const CpGrid& grid,
const double* trans,
int level)
{
assert(false);
using AttributeSet = Dune::OwnerOverlapCopyAttributeSet::AttributeSet;
bool validLevel = (level>-1) && (level <= grid.maxLevel());
const auto& ix = validLevel? grid.levelIndexSet(level) : grid.leafIndexSet();
Expand Down Expand Up @@ -571,16 +591,16 @@ int addOverlapLayer([[maybe_unused]] const CpGrid& grid,

auto it = validLevel? grid.template lbegin<0>(level) : grid.template leafbegin<0>();
const auto& endIt = validLevel? grid.template lend<0>(level) : grid.template leafend<0>();

const std::array<std::vector<std::set<int>>,2> vertex_cell_maps = grid.vertexCell();
for (; it != endIt; ++it) {
int index = ix.index(*it);
auto owner = cell_part[index];
exportProcs.insert(std::make_pair(owner, 0));
if ( trans ) {
if ( trans && false) {
addOverlapLayerNoZeroTrans(grid, index, *it, owner, cell_part, exportList, addCornerCells, layers-1, trans, level);
}
else {
addOverlapLayer(grid, index, *it, owner, cell_part, exportList, addCornerCells, layers-1, level);
addOverlapLayer(grid, index, *it, owner, cell_part, exportList, addCornerCells, vertex_cell_maps, layers-1, level);
}
}
// remove multiple entries
Expand Down
3 changes: 2 additions & 1 deletion opm/grid/common/GridPartitioning.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,8 @@ namespace Dune
bool addCornerCells,
const double* trans,
int layers = 1,
int level = -1);
int level = -1
);

namespace cpgrid
{
Expand Down
25 changes: 25 additions & 0 deletions opm/grid/cpgrid/CpGrid.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -390,6 +390,7 @@ CpGrid::scatterGrid(EdgeWeightMethod method,
comm().barrier();

// first create the overlap
auto vertex_cell_map = this->vertexCell();
auto noImportedOwner = addOverlapLayer(*this,
computedCellPart,
exportList,
Expand Down Expand Up @@ -1176,6 +1177,30 @@ int CpGrid::cellFace(int cell, int local_index, int level) const
: current_data_->back()->cell_to_face_[cpgrid::EntityRep<0>(cell, true)][local_index].index();
}

std::array<std::vector<std::set<int>>,2> CpGrid::vertexCell() const{
//int level = -1;
int nc = numCells();//(level)
int nv = this->numVertices();
std::vector<std::set<int>> vertex_cell(nv);
std::vector<std::set<int>> cell_vertex(nc);
for(int cell=0; cell< nc; ++cell){ //cells
int nlf = numCellFaces(cell);//level);
for(int lf=0; lf < nlf; ++lf){//local faces of cell
int face = this->cellFace(cell,lf);//level);
int nlv = numFaceVertices(face);//numFaceVertices(face,level);
for(int lv=0; lv < nlv; ++lv){//local nodes of face
int vertex = faceVertex(face,lv);
vertex_cell[vertex].insert(cell);
cell_vertex[cell].insert(vertex);
}
}
}
std::array<std::vector<std::set<int>>,2> out;
out[0] = vertex_cell;
out[1] = cell_vertex;
return out;
}

const cpgrid::OrientedEntityTable<0,1>::row_type CpGrid::cellFaceRow(int cell) const
{
return current_data_->back()->cell_to_face_[cpgrid::EntityRep<0>(cell, true)];
Expand Down
Loading