Skip to content

Commit

Permalink
Merge pull request #342 from totto82/fix_minpv_pinch_PR
Browse files Browse the repository at this point in the history 
Add support for nnc minpv option
  • Loading branch information
@atgeirr
atgeirr authored Sep 21, 2018
2 parents 471be7e + 3f9916e commit f54ae2f
Show file tree
Hide file tree
Showing 7 changed files with 203 additions and 64 deletions.
20 changes: 12 additions & 8 deletions opm/grid/GridManager.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -133,7 +133,6 @@ namespace Opm
const std::vector<double>& poreVolumes)
{
struct grdecl g;
int cells_modified = 0;
std::vector<int> actnum;
std::vector<double> coord;
std::vector<double> zcorn;
Expand All @@ -156,10 +155,16 @@ namespace Opm
if (!poreVolumes.empty() && (inputGrid.getMinpvMode() != MinpvMode::ModeEnum::Inactive)) {
MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
const double minpv_value = inputGrid.getMinpvValue();
// Currently the pinchProcessor is not used and only opmfil is supported
//bool opmfil = inputGrid.getMinpvMode() == MinpvMode::OpmFIL;
bool opmfil = true;
cells_modified = mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
std::vector<double> thickness(cartGridSize);
for (size_t i = 0; i < cartGridSize; ++i) {
thickness[i] = inputGrid.getCellThicknes(i);
}

// The legacy code only supports the opmfil option
bool opmfil = true; //inputGrid.getMinpvMode() == MinpvMode::OpmFIL;
const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
mp.process(thickness, z_tolerance, poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
}

const double z_tolerance = inputGrid.isPinchActive() ? inputGrid.getPinchThresholdThickness() : 0.0;
Expand All @@ -168,9 +173,8 @@ namespace Opm
OPM_THROW(std::runtime_error, "Failed to construct grid.");
}

if (cells_modified > 0) {
attach_zcorn_copy( ug_ , zcorn.data() );
}
attach_zcorn_copy( ug_ , zcorn.data() );

}


Expand Down
104 changes: 88 additions & 16 deletions opm/grid/MinpvProcessor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,9 @@


#include <opm/grid/utility/ErrorMacros.hpp>
#include <array>
#include <opm/grid/utility/OpmParserIncludes.hpp>

#include <array>

namespace Opm
{
Expand All @@ -38,6 +39,8 @@ namespace Opm
/// \param[in] nz logical cartesian number of cells in K-direction
MinpvProcessor(const int nx, const int ny, const int nz);
/// Change zcorn so that it respects the minpv property.
/// \param[in] thickness thickness of the cell
/// \param[in] z_tolerance cells with thickness below z_tolerance will be bypassed in the minpv process.
/// \param[in] pv pore volumes of all logical cartesian cells
/// \param[in] minpv minimum pore volume to accept a cell
/// \param[in] actnum active cells, inactive cells are not considered
Expand All @@ -48,7 +51,10 @@ namespace Opm
/// will have the zcorn numbers changed so they are zero-thickness. Any
/// cell below will be changed to include the deleted volume if mergeMinPCCells is true
/// els the volume will be lost
int process(const std::vector<double>& pv,
std::map<int,int> process(
const std::vector<double>& thickness,
const double z_tolerance,
const std::vector<double>& pv,
const double minpv,
const std::vector<int>& actnum,
const bool mergeMinPVCells,
Expand All @@ -68,7 +74,10 @@ namespace Opm



inline int MinpvProcessor::process(const std::vector<double>& pv,
inline std::map<int,int> MinpvProcessor::process(
const std::vector<double>& thickness,
const double z_tolerance,
const std::vector<double>& pv,
const double minpv,
const std::vector<int>& actnum,
const bool mergeMinPVCells,
Expand All @@ -81,11 +90,21 @@ namespace Opm
// pv[c] is less than minpv.
// 3. If below the minpv threshold, move the lower four
// zcorn associated with the cell c to coincide with
// the upper four (so it becomes degenerate). If mergeMinPVcells
// the upper four (so it becomes degenerate).
// 4. Look for the next active cell by skipping
// inactive cells with thickness below the z_tolerance.
// 5. If mergeMinPVcells:
// is true, the higher four zcorn associated with the cell below
// is moved to these values (so it gains the deleted volume).
// is false, a nnc is created between the cell above the removed
// cell and the cell below it. Note that the connection is only
// created if the cell below and above are active
// Inactive cells with thickness below z_tolerance and cells with porv<minpv
// are bypassed.


int cells_modified = 0;
// return a list of the non-neighbor connection.
std::map<int,int> nnc;

// Check for sane input sizes.
const size_t log_size = dims_[0] * dims_[1] * dims_[2];
Expand All @@ -103,31 +122,84 @@ namespace Opm
const int c = ii + dims_[0] * (jj + dims_[1] * kk);
if (pv[c] < minpv && (actnum.empty() || actnum[c])) {
// Move deeper (higher k) coordinates to lower k coordinates.
// i.e remove the cell
std::array<double, 8> cz = getCellZcorn(ii, jj, kk, zcorn);
for (int count = 0; count < 4; ++count) {
cz[count + 4] = cz[count];
}
setCellZcorn(ii, jj, kk, cz, zcorn);

// Find the next cell
int kk_iter = kk + 1;
if (kk_iter == dims_[2]) // we are at the end of the pillar.
continue;

int c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
// bypass inactive cells with thickness less then the tolerance
while ( ((actnum.empty() || !actnum[c_below]) && (thickness[c_below] <= z_tolerance)) ){
// move these cell to the posistion of the first cell to make the
// coordinates strictly sorted
setCellZcorn(ii, jj, kk_iter, cz, zcorn);
kk_iter ++;
if (kk_iter == dims_[2])
break;

c_below = ii + dims_[0] * (jj + dims_[1] * (kk_iter));
}

if (kk_iter == dims_[2]) // we have come to the end of the pillar.
continue;

// optionally add removed volume to the cell below.
if (mergeMinPVCells) {
// Check if there is a cell below.
if (pv[c] > 0.0 && kk < dims_[2] - 1) {
// Set lower k coordinates of cell below to upper cells's coordinates.
std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk + 1, zcorn);
for (int count = 0; count < 4; ++count) {
cz_below[count] = cz[count];
// create nnc if false or merge the cells if true
if (!mergeMinPVCells) {

// We are at the top, so no nnc is created.
if (kk == 0)
continue;

int c_above = ii + dims_[0] * (jj + dims_[1] * (kk - 1));

// Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
if (((actnum.empty() || !actnum[c_above]) && thickness[c_above] < z_tolerance) || ((actnum.empty() || actnum[c_above]) && pv[c_above] < minpv) ) {
for (int topk = kk - 2; topk > 0; --topk) {
c_above = ii + dims_[0] * (jj + dims_[1] * (topk));
if ( ((actnum.empty() || actnum[c_above]) && pv[c_above] > minpv) || ((actnum.empty() || !actnum[c_above]) && thickness[c_above] > z_tolerance)) {
break;
}
}
setCellZcorn(ii, jj, kk + 1, cz_below, zcorn);
}

// Bypass inactive cells with thickness below tolerance and active cells with volume below minpv
if (((actnum.empty() || (!actnum[c_below])) && thickness[c_below] < z_tolerance) || ((actnum.empty() || actnum[c_below]) && pv[c_below] < minpv) ) {
for (int botk = kk_iter + 1; botk < dims_[2]; ++botk) {
c_below = ii + dims_[0] * (jj + dims_[1] * (botk));
if ( ((actnum.empty() || actnum[c_below]) && pv[c_below] > minpv) || ((actnum.empty() || !actnum[c_below]) && thickness[c_below] > z_tolerance)) {
break;
}
}
}

// Add a connection if the cell above and below is active and has porv > minpv
if ((actnum.empty() || (actnum[c_above] && actnum[c_below])) && pv[c_above] > minpv && pv[c_below] > minpv) {
nnc.insert(std::make_pair(c_above, c_below));
}
} else {

// Set lower k coordinates of cell below to upper cells's coordinates.
// i.e fill the void using the cell below
std::array<double, 8> cz_below = getCellZcorn(ii, jj, kk_iter, zcorn);
for (int count = 0; count < 4; ++count) {
cz_below[count] = cz[count];
}
setCellZcorn(ii, jj, kk_iter, cz_below, zcorn);
}
++cells_modified;

}
}

}
}
return cells_modified;
return nnc;
}


Expand Down
6 changes: 4 additions & 2 deletions opm/grid/cpgrid/CpGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -248,7 +248,8 @@ CpGrid::scatterGrid(const std::vector<const cpgrid::OpmWellType *> * wells,
g.coord = &coord[0];
g.zcorn = &zcorn[0];
g.actnum = &actnum[0];
current_view_data_->processEclipseFormat(g, 0.0, false, false);
std::map<int,int> nnc;
current_view_data_->processEclipseFormat(g, nnc, 0.0, false, false);
}

void CpGrid::readSintefLegacyFormat(const std::string& grid_prefix)
Expand Down Expand Up @@ -276,7 +277,8 @@ CpGrid::scatterGrid(const std::vector<const cpgrid::OpmWellType *> * wells,
void CpGrid::processEclipseFormat(const grdecl& input_data, double z_tolerance,
bool remove_ij_boundary, bool turn_normals)
{
current_view_data_->processEclipseFormat(input_data, z_tolerance, remove_ij_boundary, turn_normals);
std::map<int,int> nnc;
current_view_data_->processEclipseFormat(input_data, nnc, z_tolerance, remove_ij_boundary, turn_normals);
}

} // namespace Dune
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 @@ -199,7 +199,7 @@ class CpGridData
/// \param z_tolerance points along a pillar that are closer together in z
/// coordinate than this parameter, will be replaced by a single point.
/// \param remove_ij_boundary if true, will remove (i, j) boundaries. Used internally.
void processEclipseFormat(const grdecl& input_data, double z_tolerance, bool remove_ij_boundary, bool turn_normals = false);
void processEclipseFormat(const grdecl& input_data, std::map<int,int>& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals = false);


/// @brief
Expand Down
56 changes: 43 additions & 13 deletions opm/grid/cpgrid/processEclipseFormat.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ namespace Dune
void removeOuterCellLayer(processed_grid& grid);
// void removeUnusedNodes(processed_grid& grid); // NOTE: not deleted, see comment at definition.
void buildTopo(const processed_grid& output,
const std::map<int,int>& nnc,
std::vector<int>& global_cell,
cpgrid::OrientedEntityTable<0, 1>& c2f,
cpgrid::OrientedEntityTable<1, 0>& f2c,
Expand Down Expand Up @@ -125,15 +126,21 @@ namespace cpgrid
g.zcorn = &zcornData[0];

g.actnum = actnumData.empty() ? nullptr : &actnumData[0];
std::map<int,int> nnc;

// Possibly process MINPV
if (!poreVolume.empty() && (ecl_grid.getMinpvMode() != Opm::MinpvMode::ModeEnum::Inactive)) {
Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
// Currently the pinchProcessor is not used and only opmfil is supported
//bool opmfil = ecl_grid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
bool opmfil = true;
size_t cells_modified = mp.process(poreVolume, ecl_grid.getMinpvValue(), actnumData, opmfil, zcornData.data());
if (cells_modified > 0) {
// Currently PINCH is always assumed to be active
bool opmfil = ecl_grid.getMinpvMode() == Opm::MinpvMode::OpmFIL;
const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
std::vector<double> thickness(cartGridSize);
for (size_t i = 0; i < cartGridSize; ++i) {
thickness[i] = ecl_grid.getCellThicknes(i);
}
const double z_tolerance = ecl_grid.isPinchActive() ? ecl_grid.getPinchThresholdThickness() : 0.0;
nnc = mp.process(thickness, z_tolerance, poreVolume, ecl_grid.getMinpvValue(), actnumData, opmfil, zcornData.data());
if (opmfil || nnc.size() > 0) {
this->zcorn = zcornData;
}
}
Expand Down Expand Up @@ -193,10 +200,10 @@ namespace cpgrid
grdecl new_g;
addOuterCellLayer(g, new_coord, new_zcorn, new_actnum, new_g);
// Make the grid.
processEclipseFormat(new_g, z_tolerance, true, turn_normals);
processEclipseFormat(new_g, nnc, z_tolerance, true, turn_normals);
} else {
// Make the grid.
processEclipseFormat(g, z_tolerance, false, turn_normals);
processEclipseFormat(g, nnc, z_tolerance, false, turn_normals);
}
}
#endif // #if HAVE_ECL_INPUT
Expand All @@ -205,7 +212,7 @@ namespace cpgrid


/// Read the Eclipse grid format ('.grdecl').
void CpGridData::processEclipseFormat(const grdecl& input_data, double z_tolerance, bool remove_ij_boundary, bool turn_normals)
void CpGridData::processEclipseFormat(const grdecl& input_data, std::map<int,int>& nnc, double z_tolerance, bool remove_ij_boundary, bool turn_normals)
{
// Process.
#ifdef VERBOSE
Expand All @@ -223,7 +230,7 @@ namespace cpgrid
std::cout << "Building topology." << std::endl;
#endif
std::vector<int> face_to_output_face;
buildTopo(output, global_cell_, cell_to_face_, face_to_cell_, face_to_point_, cell_to_point_, face_to_output_face);
buildTopo(output, nnc, global_cell_, cell_to_face_, face_to_cell_, face_to_point_, cell_to_point_, face_to_output_face);
std::copy(output.dimensions, output.dimensions + 3, logical_cartesian_size_.begin());

#ifdef VERBOSE
Expand Down Expand Up @@ -633,6 +640,7 @@ namespace cpgrid


void buildTopo(const processed_grid& output,
const std::map<int,int>& nnc,
std::vector<int>& global_cell,
cpgrid::OrientedEntityTable<0, 1>& c2f,
cpgrid::OrientedEntityTable<1, 0>& f2c,
Expand All @@ -644,6 +652,20 @@ namespace cpgrid
global_cell.assign(output.local_cell_index,
output.local_cell_index + output.number_of_cells);

std::vector<int> global_to_local;
if (nnc.size() > 0) {
int cart_size = 1;
int num_dims = sizeof(output.dimensions)/ sizeof(*output.dimensions);
for (int idx = 0; idx < num_dims ; ++idx) {
cart_size *= output.dimensions[idx];
}
// create the inverse map of global_cell
global_to_local.resize(cart_size, -1);
for (size_t idx = 0; idx < global_cell.size(); ++idx) {
global_to_local[global_cell[idx]] = idx;
}
}

// Build face to cell.
f2c.clear();
int nf = output.number_of_faces;
Expand All @@ -660,6 +682,18 @@ namespace cpgrid
cells[cellcount].setValue(fnc[1], false);
++cellcount;
}

if (output.face_tag[i] == TOP ) {
// add the NNC created from the minpv processs
// at the bottom of the cell
if (fnc[1] == -1 ) {
auto it = nnc.find(global_cell[fnc[0]]);
if (it != nnc.end()) {
cells[cellcount].setValue(global_to_local[it->second], false);
++cellcount;
}
}
}
// Assertation below is no longer true, due to periodic_extension etc.
// Instead, the appendRow() is put inside an if test.
// assert(cellcount == 1 || cellcount == 2);
Expand Down Expand Up @@ -718,10 +752,6 @@ namespace cpgrid
#endif
}





/// Encapsulate a vector<T>, and a permutation array used for access.
template <typename T>
class IndirectArray
Expand Down
9 changes: 8 additions & 1 deletion opm/grid/polyhedralgrid/grid.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -820,9 +820,16 @@ namespace Dune
Opm::MinpvProcessor mp(g.dims[0], g.dims[1], g.dims[2]);
const double minpv_value = eclipseGrid->getMinpvValue();
// Currently the pinchProcessor is not used and only opmfil is supported
// The polyhedralgrid only only supports the opmfil option
//bool opmfil = eclipseGrid->getMinpvMode() == Opm::MinpvMode::OpmFIL;
bool opmfil = true;
mp.process(poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
const size_t cartGridSize = g.dims[0] * g.dims[1] * g.dims[2];
std::vector<double> thickness(cartGridSize);
for (size_t i = 0; i < cartGridSize; ++i) {
thickness[i] = eclipseGrid->getCellThicknes(i);
}
const double z_tolerance = eclipseGrid->isPinchActive() ? eclipseGrid->getPinchThresholdThickness() : 0.0;
mp.process(thickness, z_tolerance, poreVolumes, minpv_value, actnum, opmfil, zcorn.data());
}

const double z_tolerance = eclipseGrid->isPinchActive() ?
Expand Down
Loading

0 comments on commit f54ae2f

Please sign in to comment.