Support for non constant surface densities

The surface density function returns one value pr cell to allow for non
constant surface densities.
This commit is contained in:
Tor Harald Sandve
2015-11-10 12:49:19 +01:00
parent 91af12034c
commit a46b64adcd
7 changed files with 59 additions and 31 deletions
+8 -6
View File
@@ -782,12 +782,14 @@ 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]]);
const std::vector<V> rhos = fluid_.surfaceDensity(well_cells);
// The compute density segment wants the surface densities as
// an np * number of wells cells array
V rho = superset(rhos[0], Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
rho += superset(rhos[phase], 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_));
@@ -2690,7 +2692,7 @@ namespace detail {
const ADB& rs,
const ADB& rv) const
{
const double* rhos = fluid_.surfaceDensity();
std::vector<V> rhos = fluid_.surfaceDensity(cells_);
ADB rho = rhos[phase] * b;
if (phase == Oil && active_[Gas]) {
// It is correct to index into rhos with canonical phase indices.
+18 -4
View File
@@ -349,11 +349,25 @@ 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] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of number of phases with n density values each.
std::vector<V> BlackoilPropsAdFromDeck::surfaceDensity(const Cells& cells) const
{
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
return &densities_[pvtRegionIdx][0];
const int n = cells.size();
std::vector<V> rhos(BlackoilPhases::MaxNumPhases);
for (size_t phaseIdx = 0; phaseIdx < rhos.size(); ++phaseIdx) {
rhos[phaseIdx] = V::Zero(n);
}
for (int cellIdx = 0; cellIdx < n; ++cellIdx) {
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
const double* rho = &densities_[pvtRegionIdx][0];
for (size_t phaseIdx = 0; phaseIdx < rhos.size(); ++phaseIdx) {
rhos[phaseIdx][cellIdx] = rho[phaseIdx];
}
}
return rhos;
}
+3 -2
View File
@@ -187,8 +187,9 @@ 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] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of number of phases with n density values each.
std::vector<V> surfaceDensity(const Cells& cells) const;
// ------ Viscosity ------
+3 -2
View File
@@ -88,8 +88,9 @@ 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] cells Array of n cell indices to be associated with the pressure values.
/// \return Array of number of phases with n density values each.
virtual std::vector<V> surfaceDensity(const Cells& cells) const = 0;
// ------ Viscosity ------
+10 -9
View File
@@ -341,12 +341,7 @@ namespace Opm {
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]]);
}
std::vector<V> surf_dens = fluid_.surfaceDensity(well_cells);
// Compute b, rsmax, rvmax values for perforations.
// Evaluate the properties using average well block pressures
@@ -358,6 +353,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);
@@ -408,9 +404,9 @@ 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& rhog = surf_dens[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);
surf_dens[pu.phase_pos[BlackoilPhases::Vapour]] = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
}
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
@@ -419,8 +415,13 @@ namespace Opm {
}
// b and surf_dens_perf is row major, so can just copy data.
V surf_dens_copy = superset(surf_dens[0], Span(nperf, pu.num_phases, 0), nperf*pu.num_phases);
for (int phase = 1; phase < pu.num_phases; ++phase) {
surf_dens_copy += superset(surf_dens[phase], Span(nperf, pu.num_phases, phase), nperf*pu.num_phases);
}
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_);
+4 -4
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();
std::vector<V> rhos = fluid_.surfaceDensity(cells);
V b = fluidFvf(phase, p, T, cells);
V rho = V::Constant(p.size(), 1, rhos[phase]) * b;
V rho = rhos[phase] * 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();
std::vector<V> rhos = fluid_.surfaceDensity(cells);
ADB b = fluidFvf(phase, p, T, cells);
ADB rho = V::Constant(p.size(), 1, rhos[phase]) * b;
ADB rho = rhos[phase] * b;
return rho;
}
+13 -4
View File
@@ -115,16 +115,25 @@ 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;
std::vector<V> rho0AD = boprops_ad.surfaceDensity(cells);
BOOST_REQUIRE_EQUAL(rho0AD.size(), 3);
enum { Water = Opm::BlackoilPropsAdFromDeck::Water };
BOOST_CHECK_EQUAL(rho0AD[ Water ], 1000.0);
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);
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);
BOOST_REQUIRE_EQUAL(rho0AD[Gas].size(), cells.size());
BOOST_CHECK_EQUAL(rho0AD[ Gas ][0], 1.0);
}