Merge pull request #320 from atgeirr/fix-pvtregions

Fix pvt regions
This commit is contained in:
Bård Skaflestad 2015-02-27 12:29:42 +01:00
commit 90ff9c055d
5 changed files with 79 additions and 57 deletions

View File

@ -90,7 +90,6 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
// For data that is dependant on the subgrid we simply allocate space
// and initialize with obviously bogus numbers.
cellPvtRegionIdx_.resize(number_of_cells, std::numeric_limits<int>::min());
pvtTableIdx_.resize(number_of_cells, std::numeric_limits<int>::min());
}
/// Initializes the properties.
@ -136,14 +135,6 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
}
}
// first, calculate the PVT table index for each compressed
// cell. This array is required to construct the PVT classes
// below.
Opm::extractPvtTableIndex(pvtTableIdx_,
deck,
number_of_cells,
global_cell);
const int numSamples = 0;
// Resize the property objects container
@ -333,13 +324,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pw.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.data(), T.data(), rs,
props_[phase_usage_.phase_pos[Water]]->mu(n, pvt_region_.data(), pw.data(), T.data(), rs,
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
@ -361,12 +353,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.data(), T.data(), rs.data(), &cond[0],
props_[phase_usage_.phase_pos[Oil]]->mu(n, pvt_region_.data(), po.data(), T.data(), rs.data(), &cond[0],
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
@ -384,13 +377,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), T.data(), rs,
props_[phase_usage_.phase_pos[Gas]]->mu(n, pvt_region_.data(), pg.data(), T.data(), rs,
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
@ -410,12 +404,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.data(), T.data(), rv.data(),&cond[0],
props_[phase_usage_.phase_pos[Gas]]->mu(n, pvt_region_.data(), pg.data(), T.data(), rv.data(),&cond[0],
mu.data(), dmudp.data(), dmudr.data());
return mu;
}
@ -433,13 +428,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pw.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->mu(n, &pvtTableIdx_[0], pw.value().data(), T.value().data(), rs,
props_[phase_usage_.phase_pos[Water]]->mu(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs,
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
const int num_blocks = pw.numBlocks();
@ -467,12 +463,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Oil]]->mu(n, &pvtTableIdx_[0], po.value().data(), T.value().data(), rs.value().data(),
props_[phase_usage_.phase_pos[Oil]]->mu(n, pvt_region_.data(), po.value().data(), T.value().data(), rs.value().data(),
&cond[0], mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -501,13 +498,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.value().size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
const double* rv = 0;
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv,
props_[phase_usage_.phase_pos[Gas]]->mu(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv,
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -536,12 +534,13 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.value().size() == n);
V mu(n);
V dmudp(n);
V dmudr(n);
props_[phase_usage_.phase_pos[Gas]]->mu(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv.value().data(),&cond[0],
props_[phase_usage_.phase_pos[Gas]]->mu(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv.value().data(),&cond[0],
mu.data(), dmudp.data(), dmudr.data());
ADB::M dmudp_diag = spdiag(dmudp);
@ -588,6 +587,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call bWat(): water phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pw.size() == n);
V b(n);
@ -595,7 +595,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.data(), T.data(), rs,
props_[phase_usage_.phase_pos[Water]]->b(n, pvt_region_.data(), pw.data(), T.data(), rs,
b.data(), dbdp.data(), dbdr.data());
return b;
@ -618,13 +618,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call bOil(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V b(n);
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.data(), T.data(), rs.data(), &cond[0],
props_[phase_usage_.phase_pos[Oil]]->b(n, pvt_region_.data(), po.data(), T.data(), rs.data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
return b;
@ -643,6 +644,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call bGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V b(n);
@ -650,7 +652,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), T.data(), rs,
props_[phase_usage_.phase_pos[Gas]]->b(n, pvt_region_.data(), pg.data(), T.data(), rs,
b.data(), dbdp.data(), dbdr.data());
return b;
@ -673,13 +675,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V b(n);
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.data(), T.data(), rv.data(), &cond[0],
props_[phase_usage_.phase_pos[Gas]]->b(n, pvt_region_.data(), pg.data(), T.data(), rv.data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
return b;
@ -698,6 +701,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muWat(): water phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pw.size() == n);
V b(n);
@ -705,7 +709,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V dbdr(n);
const double* rs = 0;
props_[phase_usage_.phase_pos[Water]]->b(n, &pvtTableIdx_[0], pw.value().data(), T.value().data(), rs,
props_[phase_usage_.phase_pos[Water]]->b(n, pvt_region_.data(), pw.value().data(), T.value().data(), rs,
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -734,13 +738,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muOil(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V b(n);
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Oil]]->b(n, &pvtTableIdx_[0], po.value().data(), T.value().data(), rs.value().data(),
props_[phase_usage_.phase_pos[Oil]]->b(n, pvt_region_.data(), po.value().data(), T.value().data(), rs.value().data(),
&cond[0], b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -769,6 +774,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V b(n);
@ -776,7 +782,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
V dbdr(n);
const double* rv = 0;
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv,
props_[phase_usage_.phase_pos[Gas]]->b(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv,
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -805,13 +811,14 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call muGas(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(pg.size() == n);
V b(n);
V dbdp(n);
V dbdr(n);
props_[phase_usage_.phase_pos[Gas]]->b(n, &pvtTableIdx_[0], pg.value().data(), T.value().data(), rv.value().data(), &cond[0],
props_[phase_usage_.phase_pos[Gas]]->b(n, pvt_region_.data(), pg.value().data(), T.value().data(), rv.value().data(), &cond[0],
b.data(), dbdp.data(), dbdr.data());
ADB::M dbdp_diag = spdiag(dbdp);
@ -842,10 +849,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call rsMax(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V rbub(n);
V drbubdp(n);
props_[phase_usage_.phase_pos[Oil]]->rsSat(n, &pvtTableIdx_[0], po.data(), rbub.data(), drbubdp.data());
props_[phase_usage_.phase_pos[Oil]]->rsSat(n, pvt_region_.data(), po.data(), rbub.data(), drbubdp.data());
return rbub;
}
@ -874,10 +882,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call rsMax(): oil phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V rbub(n);
V drbubdp(n);
props_[phase_usage_.phase_pos[Oil]]->rsSat(n, &pvtTableIdx_[0], po.value().data(), rbub.data(), drbubdp.data());
props_[phase_usage_.phase_pos[Oil]]->rsSat(n, pvt_region_.data(), po.value().data(), rbub.data(), drbubdp.data());
ADB::M drbubdp_diag = spdiag(drbubdp);
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
@ -914,10 +923,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call rvMax(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V rv(n);
V drvdp(n);
props_[phase_usage_.phase_pos[Gas]]->rvSat(n, &pvtTableIdx_[0], po.data(), rv.data(), drvdp.data());
props_[phase_usage_.phase_pos[Gas]]->rvSat(n, pvt_region_.data(), po.data(), rv.data(), drvdp.data());
return rv;
}
@ -946,10 +956,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
OPM_THROW(std::runtime_error, "Cannot call rvMax(): gas phase not present.");
}
const int n = cells.size();
mapPvtRegions(cells);
assert(po.size() == n);
V rv(n);
V drvdp(n);
props_[phase_usage_.phase_pos[Gas]]->rvSat(n, &pvtTableIdx_[0], po.value().data(), rv.data(), drvdp.data());
props_[phase_usage_.phase_pos[Gas]]->rvSat(n, pvt_region_.data(), po.value().data(), rv.data(), drvdp.data());
ADB::M drvdp_diag = spdiag(drvdp);
const int num_blocks = po.numBlocks();
std::vector<ADB::M> jacs(num_blocks);
@ -1241,5 +1252,20 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
}
}
// Fills pvt_region_ with cellPvtRegionIdx_[cells].
void BlackoilPropsAdFromDeck::mapPvtRegions(const std::vector<int>& cells) const
{
const int n = cells.size();
pvt_region_.resize(n);
for (int ii = 0; ii < n; ++ii) {
pvt_region_[ii] = cellPvtRegionIdx_[cells[ii]];
}
}
} // namespace Opm

View File

@ -455,6 +455,9 @@ namespace Opm
const std::vector<int>& cells,
const double vap) const;
// Fills pvt_region_ with cellPvtRegionIdx_[cells].
void mapPvtRegions(const std::vector<int>& cells) const;
RockFromDeck rock_;
// This has to be a shared pointer as we must
// be able to make a copy of *this in the parallel case.
@ -467,18 +470,14 @@ namespace Opm
// The PVT region which is to be used for each cell
std::vector<int> cellPvtRegionIdx_;
// Used for storing the region-per-cell array computed in calls
// to pvt functions.
mutable std::vector<int> pvt_region_;
// The PVT properties. One object per active fluid phase.
std::vector<std::shared_ptr<Opm::PvtInterface> > props_;
// The index of the PVT table which ought to be used for each
// cell. Eclipse does not seem to allow specifying fluid-phase
// specific table indices, so for the sake of simplicity, we
// don't do this either. (if it turns out that Eclipes does in
// fact support it or if it by some miracle gains this feature
// in the future, this attribute needs to become a vector of
// vectors of ints.)
std::vector<int> pvtTableIdx_;
// Densities, one std::array per PVT region.
std::vector<std::array<double, BlackoilPhases::MaxNumPhases> > densities_;
// VAPPARS

View File

@ -145,12 +145,12 @@ public:
BlackoilPropsDataHandle(const BlackoilPropsAdFromDeck& sendProps,
BlackoilPropsAdFromDeck& recvProps)
: sendProps_(sendProps), recvProps_(recvProps),
size_(2)
size_(1)
{
// satOilMax might be non empty. In this case we will need to send it, too.
if ( sendProps.satOilMax_.size()>0 )
{
recvProps_.satOilMax_.resize(recvProps_.pvtTableIdx_.size(),
recvProps_.satOilMax_.resize(recvProps_.satOilMax_.size(),
-std::numeric_limits<double>::max());
++size_;
}
@ -166,7 +166,7 @@ public:
{
if ( T::codimension == 0)
{
// We only send pvtTableIdx_, cellPvtRegionIdx_, and maybe satOilMax_
// We only send cellPvtRegionIdx_, and maybe satOilMax_
return size_;
}
else
@ -180,12 +180,9 @@ public:
assert( T::codimension == 0);
buffer.write(sendProps_.cellPvtRegionIndex()[e.index()]);
buffer.write(sendProps_.pvtTableIdx_[e.index()]);
if ( size_==2 )
{
return;
if ( size_ > 1 ) {
buffer.write(sendProps_.satOilMax_[e.index()]);
}
buffer.write(sendProps_.satOilMax_[e.index()]);
}
template<class B, class T>
void scatter(B& buffer, const T& e, std::size_t size)
@ -195,12 +192,10 @@ public:
double val;
buffer.read(val);
recvProps_.cellPvtRegionIdx_[e.index()]=val;
buffer.read(val);
recvProps_.pvtTableIdx_[e.index()]=val;
if ( size_==2 )
return;
buffer.read(val);
recvProps_.satOilMax_[e.index()]=val;
if ( size_ > 1 ) {
buffer.read(val);
recvProps_.satOilMax_[e.index()]=val;
}
}
bool contains(int dim, int codim)
{

View File

@ -49,8 +49,9 @@ PVCDO
/
PVDG
-- Pg Bg(Pg) mug
1 1 1
-- Pg Bg(Pg) mug
1 1 1
800 0.99999999 1
/
SWOF

View File

@ -94,7 +94,7 @@ BOOST_FIXTURE_TEST_CASE(Construction, TestFixture<SetupSimple>)
}
BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture<SetupSimple>)
BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture<SetupSimple>)
{
// Immiscible and incompressible two-phase fluid
typedef std::vector<int> Region;
@ -106,11 +106,11 @@ BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture<SetupSimple>)
RCvrt cvrt(ad_props, reg);
Opm::BlackoilState x;
x.init(*grid.c_grid(), 2);
x.init(*grid.c_grid(), 3);
cvrt.defineState(x);
std::vector<double> qs{1.0e3, 1.0e1};
std::vector<double> qs{1.0e3, 1.0e1, 1.0e-1};
std::vector<double> coeff(qs.size(), 0.0);
// Immiscible and incompressible: All coefficients are one (1),
@ -118,4 +118,5 @@ BOOST_FIXTURE_TEST_CASE(TwoPhaseII, TestFixture<SetupSimple>)
cvrt.calcCoeff(qs, 0, coeff);
BOOST_CHECK_CLOSE(coeff[0], 1.0, 1.0e-6);
BOOST_CHECK_CLOSE(coeff[1], 1.0, 1.0e-6);
BOOST_CHECK_CLOSE(coeff[2], 1.0, 1.0e-6);
}