Merge pull request #520 from totto82/surface_dens

Support for non-constant surface densities
This commit is contained in:
Atgeirr Flø Rasmussen 2015-11-13 11:40:30 +01:00
commit 2a0051655c
7 changed files with 61 additions and 39 deletions

View File

@ -782,12 +782,13 @@ namespace detail {
std::vector<double> perf_depth(pdepth.data(), pdepth.data() + nperf);
// Surface density.
DataBlock surf_dens(nperf, pu.num_phases);
for (int phase = 0; phase < pu.num_phases; ++ phase) {
surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]);
// The compute density segment wants the surface densities as
// an np * number of wells cells array
V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
std::vector<double> surf_dens_perf(rho.data(), rho.data() + nperf * pu.num_phases);
// Gravity
double grav = detail::getGravity(geo_.gravity(), dimensions(grid_));
@ -2769,15 +2770,14 @@ namespace detail {
const ADB& rs,
const ADB& rv) const
{
const double* rhos = fluid_.surfaceDensity();
ADB rho = rhos[phase] * b;
const V& rhos = fluid_.surfaceDensity(phase, cells_);
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
ADB rho = rhos * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.
rho += rhos[Gas] * rs * b;
rho += fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * rs * b;
}
if (phase == Gas && active_[Oil]) {
// It is correct to index into rhos with canonical phase indices.
rho += rhos[Oil] * rv * b;
rho += fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * rv * b;
}
return rho;
}

View File

@ -349,11 +349,20 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
const double* BlackoilPropsAdFromDeck::surfaceDensity(const int cellIdx) const
/// \param[in] phaseIdx
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n density values for phase given by phaseIdx.
V BlackoilPropsAdFromDeck::surfaceDensity(const int phaseIdx, const Cells& cells) const
{
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
return &densities_[pvtRegionIdx][0];
assert( !(phaseIdx > numPhases()));
const int n = cells.size();
V rhos = V::Zero(n);
for (int cellIdx = 0; cellIdx < n; ++cellIdx) {
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
const auto* rho = &densities_[pvtRegionIdx][0];
rhos[cellIdx] = rho[phaseIdx];
}
return rhos;
}

View File

@ -187,8 +187,10 @@ namespace Opm
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
const double* surfaceDensity(const int cellIdx = 0) const;
/// \param[in] phaseIdx
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n density values for phase given by phaseIdx.
V surfaceDensity(const int phaseIdx , const Cells& cells) const;
// ------ Viscosity ------

View File

@ -88,8 +88,10 @@ namespace Opm
// ------ Density ------
/// Densities of stock components at surface conditions.
/// \return Array of 3 density values.
virtual const double* surfaceDensity(int regionIdx = 0) const = 0;
/// \param[in] phaseIdx
/// \param[in] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of n density values for phase given by phaseIdx.
virtual V surfaceDensity(const int PhaseIdx, const Cells& cells) const = 0;
// ------ Viscosity ------

View File

@ -340,14 +340,6 @@ namespace Opm {
// Use cell values for the temperature as the wells don't knows its temperature yet.
const ADB perf_temp = subset(state.temperature, well_cells);
// Surface density.
const PhaseUsage& pu = fluid_.phaseUsage();
//std::vector<double> surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases);
DataBlock surf_dens(nperf, pu.num_phases);
for (int phase = 0; phase < pu.num_phases; ++ phase) {
surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]);
}
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
// and cell values for rs, rv, phase condition and temperature.
@ -358,6 +350,7 @@ namespace Opm {
perf_cond[perf] = pc[well_cells[perf]];
}
const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases);
std::vector<double> rsmax_perf(nperf, 0.0);
std::vector<double> rvmax_perf(nperf, 0.0);
@ -374,10 +367,18 @@ namespace Opm {
const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
}
V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
if ( phase != pu.phase_pos[BlackoilPhases::Vapour]) {
continue; // the gas surface density is added after the solvent is accounted for.
}
surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells);
V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value();
V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells);
if (has_solvent_) {
const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
// A weighted sum of the b-factors of gas and solvent are used.
@ -408,11 +409,11 @@ namespace Opm {
bg = bg * (ones - F_solvent);
bg = bg + F_solvent * bs;
const V& rhog = surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]);
const V& rhos = solvent_props_.solventSurfaceDensity(well_cells);
surf_dens.col(pu.phase_pos[BlackoilPhases::Vapour]) = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases);
const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
@ -420,7 +421,7 @@ namespace Opm {
// b and surf_dens_perf is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
std::vector<double> surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases);
std::vector<double> surf_dens_perf(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases);
// Extract well connection depths.
const V depth = cellCentroidsZToEigen(grid_);

View File

@ -623,9 +623,9 @@ namespace {
V ImpesTPFAAD::fluidRho(const int phase, const V& p, const V& T, const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
V rho = fluid_.surfaceDensity(phase, cells);
V b = fluidFvf(phase, p, T, cells);
V rho = V::Constant(p.size(), 1, rhos[phase]) * b;
rho = rho * b;
return rho;
}
@ -635,9 +635,9 @@ namespace {
ADB ImpesTPFAAD::fluidRho(const int phase, const ADB& p, const ADB& T, const std::vector<int>& cells) const
{
const double* rhos = fluid_.surfaceDensity();
const V& rhos = fluid_.surfaceDensity(phase, cells);
ADB b = fluidFvf(phase, p, T, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
ADB rho = rhos * b;
return rho;
}

View File

@ -115,16 +115,24 @@ BOOST_FIXTURE_TEST_CASE(SubgridConstruction, TestFixtureAd<SetupSimple>)
BOOST_FIXTURE_TEST_CASE(SurfaceDensity, TestFixture<SetupSimple>)
{
const double* rho0AD = boprops_ad.surfaceDensity();
const Opm::BlackoilPropsAdFromDeck::Cells cells(1, 0);
typedef Opm::BlackoilPropsAdFromDeck::V V;
enum { Water = Opm::BlackoilPropsAdFromDeck::Water };
BOOST_CHECK_EQUAL(rho0AD[ Water ], 1000.0);
V rho0AD_Water = boprops_ad.surfaceDensity(Water, cells);
BOOST_REQUIRE_EQUAL(rho0AD_Water.size(), cells.size());
BOOST_CHECK_EQUAL(rho0AD_Water[0], 1000.0);
enum { Oil = Opm::BlackoilPropsAdFromDeck::Oil };
BOOST_CHECK_EQUAL(rho0AD[ Oil ], 800.0);
V rho0AD_Oil = boprops_ad.surfaceDensity(Oil, cells);
BOOST_REQUIRE_EQUAL(rho0AD_Oil.size(), cells.size());
BOOST_CHECK_EQUAL(rho0AD_Oil[0], 800.0);
enum { Gas = Opm::BlackoilPropsAdFromDeck::Gas };
BOOST_CHECK_EQUAL(rho0AD[ Gas ], 1.0);
V rho0AD_Gas = boprops_ad.surfaceDensity(Gas, cells);
BOOST_REQUIRE_EQUAL(rho0AD_Gas.size(), cells.size());
BOOST_CHECK_EQUAL(rho0AD_Gas[0], 1.0);
}