OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

RegionMapping<>: Support arbitrary region IDs

Browse Source 
This commit introduces a new public method, activeRegions(), that
retrieves those region IDs that contain at least one active cell.
We furthermore extend the cells() method to support lookup of
arbitrary region IDs.  Non-active region IDs produce empty cell
ranges.

Intended use case is

    for (const auto& reg : rmap.activeRegions()) {
        const auto& c = rmap.cells(reg);

        // use c
    }
This commit is contained in:
Bård Skaflestad 2015-09-02 10:36:07 +02:00
parent a439f7bd88
commit c010d4c3e2
4 changed files with 167 additions and 109 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
opm/core/simulator/initStateEquil.hpp
View File

@ -399,13 +399,10 @@ namespace Opm
const Grid& G ,
const double grav)
{
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
{
const typename RMap::CellRange cells = reg.cells(r);
for (const auto& r : reg.activeRegions()) {
const auto& cells = reg.cells(r);
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],

130
opm/core/utility/RegionMapping.hpp
View File

@ -20,6 +20,9 @@
#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED
#define OPM_REGIONMAPPING_HEADER_INCLUDED
#include <boost/range.hpp>
#include <unordered_map>
#include <vector>
namespace Opm
@ -47,7 +50,7 @@ namespace Opm
*/
explicit
RegionMapping(const Region& reg)
: reg_(reg)
: reg_(reg)
{
rev_.init(reg_);
}
@ -70,58 +73,7 @@ namespace Opm
*/
typedef typename std::vector<CellId>::const_iterator CellIter;
/**
* Range of cells. Result from reverse (region-to-cell)
* mapping.
*/
class CellRange {
public:
/**
* Constructor.
*
* \param[in] b Beginning of range.
* \param[in] e One past end of range.
*/
CellRange(const CellIter b,
const CellIter e)
: b_(b), e_(e)
{}
/**
* Read-only iterator on cell ranges.
*/
typedef CellIter const_iterator;
/**
* Size type for this range.
*/
typedef typename std::vector<CellId>::size_type size_type;
/**
* Beginning of cell range.
*/
const_iterator begin() const { return b_; }
/**
* One past end of cell range.
*/
const_iterator end() const { return e_; }
/**
* Number of elements in the range.
*/
size_type size() const { return e_ - b_; }
private:
const_iterator b_;
const_iterator e_;
};
/**
* Number of declared regions in cell-to-region mapping.
*/
RegionId
numRegions() const { return RegionId(rev_.p.size()) - 1; }
typedef boost::iterator_range<CellIter> Range;
/**
* Compute region number of given active cell.
@ -132,17 +84,33 @@ namespace Opm
RegionId
region(const CellId c) const { return reg_[c]; }
const std::vector<RegionId>&
activeRegions() const
{
return rev_.active;
}
/**
* Extract active cells in particular region.
*
* \param[in] r Region number
* \returns Range of active cells in region @c r.
*
* \return Range of active cells in region @c r. Empty if @c r is
* not an active region.
*/
CellRange
Range
cells(const RegionId r) const {
const RegionId i = r - rev_.low;
return CellRange(rev_.c.begin() + rev_.p[i + 0],
rev_.c.begin() + rev_.p[i + 1]);
const auto id = rev_.binid.find(r);
if (id == rev_.binid.end()) {
// Region 'r' not an active region. Return empty.
return Range(rev_.c.end(), rev_.c.end());
}
const auto i = id->second;
return Range(rev_.c.begin() + rev_.p[i + 0],
rev_.c.begin() + rev_.p[i + 1]);
}
private:
@ -156,9 +124,12 @@ namespace Opm
*/
struct {
typedef typename std::vector<CellId>::size_type Pos;
std::unordered_map<RegionId, Pos> binid;
std::vector<RegionId> active;
std::vector<Pos> p; /**< Region start pointers */
std::vector<CellId> c; /**< Region cells */
RegionId low; /**< Smallest region number */
/**
* Compute reverse mapping. Standard linear insertion
@ -167,32 +138,37 @@ namespace Opm
void
init(const Region& reg)
{
typedef typename Region::const_iterator CI;
const std::pair<CI,CI>
m = std::minmax_element(reg.begin(), reg.end());
low = *m.first;
const typename Region::size_type
n = *m.second - low + 1;
p.resize(n + 1); std::fill(p.begin(), p.end(), Pos(0));
for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
p[ reg[i] - low + 1 ] += 1;
binid.clear();
for (const auto& r : reg) {
++binid[r];
}
for (typename std::vector<Pos>::size_type
i = 1, sz = p.size(); i < sz; ++i) {
p .clear(); p.emplace_back(0);
active.clear();
{
Pos n = 0;
for (auto& id : binid) {
active.push_back(id.first);
p .push_back(id.second);
id.second = n++;
}
}
for (decltype(p.size()) i = 1, sz = p.size(); i < sz; ++i) {
p[0] += p[i];
p[i] = p[0] - p[i];
}
assert (p[0] ==
static_cast<typename Region::size_type>(reg.size()));
assert (p[0] == static_cast<Pos>(reg.size()));
c.resize(reg.size());
for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
c[ p[ reg[i] - low + 1 ] ++ ] = i;
{
CellId i = 0;
for (const auto& r : reg) {
auto& pos = p[ binid[r] + 1 ];
c[ pos++ ] = i++;
}
}
p[0] = 0;

17
tests/test_equil.cpp
View File

@ -307,10 +307,8 @@ BOOST_AUTO_TEST_CASE (RegMapping)
Opm::RegionMapping<> eqlmap(eqlnum);
PPress ppress(2, PVal(G->number_of_cells, 0));
for (int r = 0, e = eqlmap.numRegions(); r != e; ++r)
{
const Opm::RegionMapping<>::CellRange&
rng = eqlmap.cells(r);
for (const auto& r : eqlmap.activeRegions()) {
const auto& rng = eqlmap.cells(r);
const int rno = r;
const double grav = 10;
@ -318,14 +316,13 @@ BOOST_AUTO_TEST_CASE (RegMapping)
Opm::Equil::phasePressures(*G, region[rno], rng, grav);
PVal::size_type i = 0;
for (Opm::RegionMapping<>::CellRange::const_iterator
c = rng.begin(), ce = rng.end();
c != ce; ++c, ++i)
{
for (const auto& c : rng) {
assert (i < p[0].size());
ppress[0][*c] = p[0][i];
ppress[1][*c] = p[1][i];
ppress[0][c] = p[0][i];
ppress[1][c] = p[1][i];
++i;
}
}

120
tests/test_regionmapping.cpp
View File

@ -26,42 +26,130 @@
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE UnitsTest
#define BOOST_TEST_MODULE RegionMapping
#include <opm/core/utility/platform_dependent/disable_warnings.h>
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <opm/core/utility/platform_dependent/reenable_warnings.h>
/* --- our own headers --- */
#include <opm/core/utility/RegionMapping.hpp>
#include <algorithm>
#include <map>
BOOST_AUTO_TEST_SUITE ()
BOOST_AUTO_TEST_CASE (RegionMapping)
BOOST_AUTO_TEST_CASE (Forward)
{
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
for (size_t i = 0; i < regions.size(); ++i) {
for (decltype(regions.size()) i = 0, n = regions.size(); i < n; ++i) {
BOOST_CHECK_EQUAL(rm.region(i), regions[i]);
}
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
std::vector< std::vector<int> > region_cells = { { 0, 2, 4 }, { 7 }, { 3 }, { 1 }, { 6, 8 }, { 5 } };
BOOST_REQUIRE_EQUAL(rm.numRegions(), region_ids.size());
for (size_t i = 0; i < region_ids.size(); ++i) {
auto cells = rm.cells(region_ids[i]);
BOOST_REQUIRE_EQUAL(std::distance(cells.begin(), cells.end()), region_cells[i].size());
size_t count = 0;
for (auto iter = cells.begin(); iter != cells.end(); ++iter, ++count) {
BOOST_CHECK_EQUAL(*iter, region_cells[i][count]);
}
}
}
BOOST_AUTO_TEST_CASE (ActiveRegions)
{
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
auto active = [&region_ids](const int reg)
{
auto b = region_ids.begin();
auto e = region_ids.end();
return std::find(b, e, reg) != e;
};
BOOST_CHECK_EQUAL(rm.activeRegions().size(), region_ids.size());
for (const auto& reg : rm.activeRegions()) {
BOOST_CHECK(active(reg));
}
}
BOOST_AUTO_TEST_CASE (Consecutive)
{
using RegionCells = std::map<int, std::vector<int>>;
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
RegionCells region_cells;
{
using VT = RegionCells::value_type;
region_cells.insert(VT(2, { 0, 2, 4 }));
region_cells.insert(VT(3, { 7 }));
region_cells.insert(VT(4, { 3 }));
region_cells.insert(VT(5, { 1 }));
region_cells.insert(VT(6, { 6, 8 }));
region_cells.insert(VT(7, { 5 }));
}
for (const auto& reg : region_ids) {
const auto& cells = rm.cells(reg);
const auto& expect = region_cells[reg];
BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
expect.begin(), expect.end());
}
// Verify that there are no cells in unused regions 0 and 1.
for (const auto& r : { 0, 1 }) {
BOOST_CHECK(rm.cells(r).empty());
}
}
BOOST_AUTO_TEST_CASE (NonConsecutive)
{
using RegionCells = std::map<int, std::vector<int>>;
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 4, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
std::vector<int> region_ids = { 2, 3, 4, 6, 7 };
RegionCells region_cells;
{
using VT = RegionCells::value_type;
region_cells.insert(VT(2, { 0, 2, 4 }));
region_cells.insert(VT(3, { 7 }));
region_cells.insert(VT(4, { 1, 3 }));
region_cells.insert(VT(6, { 6, 8 }));
region_cells.insert(VT(7, { 5 }));
}
for (const auto& reg : region_ids) {
const auto& cells = rm.cells(reg);
const auto& expect = region_cells[reg];
BOOST_CHECK_EQUAL_COLLECTIONS(cells .begin(), cells .end(),
expect.begin(), expect.end());
}
// Verify that there are no cells in unused regions 0, 1, and 5.
for (const auto& r : { 0, 1, 5 }) {
BOOST_CHECK(rm.cells(r).empty());
}
}
BOOST_AUTO_TEST_SUITE_END()