From aa9e28a5b83dd0c493b89d2c824be06ddfec7822 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 17 Jan 2014 20:07:51 +0100 Subject: [PATCH] Add reverse look-up mapping for region vectors Class RegionMapping<> provides an easy way of extracting the cells that belong to any identified region (e.g., as defined by EQLNUM) of the deck. --- opm/core/simulator/initStateEquil.hpp | 91 +++++++++++++++++++++ tests/test_equil.cpp | 110 ++++++++++++++++++++++++++ 2 files changed, 201 insertions(+) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index b81183cf..5df47a0c 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -27,6 +27,7 @@ #include #include +#include #include struct UnstructuredGrid; @@ -102,6 +103,96 @@ namespace Opm }; } // namespace miscibility + template < class Region = std::vector > + class RegionMapping { + public: + explicit + RegionMapping(const Region& reg) + : reg_(reg) + { + rev_.init(reg_); + } + + typedef typename Region::value_type RegionId; + typedef typename Region::size_type CellId; + typedef typename std::vector::const_iterator CellIter; + + class CellRange { + public: + CellRange(const CellIter b, + const CellIter e) + : b_(b), e_(e) + {} + + typedef CellIter iterator; + typedef CellIter const_iterator; + + iterator begin() const { return b_; } + iterator end() const { return e_; } + + private: + iterator b_; + iterator e_; + }; + + RegionId + numRegions() const { return RegionId(rev_.p.size()) - 1; } + + RegionId + region(const CellId c) const { return reg_[c]; } + + CellRange + 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]); + } + + private: + Region reg_; + + struct { + typedef typename std::vector::size_type Pos; + std::vector p; + std::vector c; + RegionId low; + + void + init(const Region& reg) + { + typedef typename Region::const_iterator CI; + const std::pair + 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; + } + + for (typename std::vector::size_type + i = 1, sz = p.size(); i < sz; ++i) { + p[0] += p[i]; + p[i] = p[0] - p[i]; + } + + assert (p[0] == + static_cast(reg.size())); + + c.resize(reg.size()); + for (CellId i = 0, nc = reg.size(); i < nc; ++i) { + c[ p[ reg[i] - low + 1 ] ++ ] = i; + } + + p[0] = 0; + } + } rev_; + }; + struct EquilRecord { struct { double depth; diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 36d1d452..80ae2652 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -25,6 +25,8 @@ #include #include +#include + #include #include @@ -201,4 +203,112 @@ BOOST_AUTO_TEST_CASE (CellSubset) BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); } +BOOST_AUTO_TEST_CASE (RegMapping) +{ + typedef std::vector PVal; + typedef std::vector PPress; + + std::shared_ptr + G(create_grid_cart3d(10, 1, 10), destroy_grid); + + Opm::parameter::ParameterGroup param; + { + using Opm::unit::kilogram; + using Opm::unit::meter; + using Opm::unit::cubic; + + std::stringstream dens; dens << 700*kilogram/cubic(meter); + param.insertParameter("rho2", dens.str()); + } + + typedef Opm::BlackoilPropertiesBasic Props; + Props props(param, G->dimensions, G->number_of_cells); + + typedef Opm::equil::DensityCalculator RhoCalc; + RhoCalc calc(props, 0); + + Opm::equil::EquilRecord record[] = + { + { + { 0 , 1e5 } , // Datum depth, pressure + { 2.5 , -0.075e5 } , // Zwoc , Pcow_woc + { 0 , 0 } // Zgoc , Pcgo_goc + } + , + { + { 5 , 1.35e5 } , // Datum depth, pressure + { 7.5 , -0.225e5 } , // Zwoc , Pcow_woc + { 5 , 0 } // Zgoc , Pcgo_goc + } + }; + + Opm::equil::EquilReg region[] = + { + Opm::equil::EquilReg(record[0], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[0], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + , + Opm::equil::EquilReg(record[1], calc, + Opm::equil::miscibility::NoMixing(), + Opm::equil::miscibility::NoMixing(), + props.phaseUsage()) + }; + + std::vector eqlnum(G->number_of_cells); + { + std::vector cells(G->number_of_cells); + std::iota(cells.begin(), cells.end(), 0); + + const int cdim[] = { 2, 1, 2 }; + int ncoarse = cdim[0]; + for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; } + + partition_unif_idx(G->dimensions, G->number_of_cells, + G->cartdims, cdim, + &cells[0], &eqlnum[0]); + } + Opm::equil::RegionMapping<> eqlmap(eqlnum); + + PPress ppress(2, PVal(G->number_of_cells, 0)); + for (int r = 0, e = eqlmap.numRegions(); r != e; ++r) + { + const Opm::equil::RegionMapping<>::CellRange& + rng = eqlmap.cells(r); + + const int rno = r; + const double grav = 10; + const PPress p = + Opm::equil::phasePressures(*G, region[rno], rng, grav); + + PVal::size_type i = 0; + for (Opm::equil::RegionMapping<>::CellRange::const_iterator + c = rng.begin(), ce = rng.end(); + c != ce; ++c, ++i) + { + assert (i < p[0].size()); + + ppress[0][*c] = p[0][i]; + ppress[1][*c] = p[1][i]; + } + } + + const int first = 0, last = G->number_of_cells - 1; + const double reltol = 1.0e-8; + BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol); + BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol); + BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol); + BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol); +} + BOOST_AUTO_TEST_SUITE_END()