From 7bdd91d78f9a01589e18dd5d84e785086729f2b3 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 21 Sep 2015 10:57:12 +0200 Subject: [PATCH] Allow for different surface densities in well perforations The surface density input in well_perforation_densities() in WellDensitySegmented.hpp is changed from one value pr. phase to one value pr phase and perforation. This allow for different densities in different perforation. The test is changed accordingly. --- opm/autodiff/BlackoilModelBase_impl.hpp | 11 +++++++++-- opm/autodiff/WellDensitySegmented.cpp | 11 ++++++++--- opm/autodiff/WellDensitySegmented.hpp | 4 ++-- tests/test_welldensitysegmented.cpp | 11 ++++++++++- 4 files changed, 29 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 9bf85ee78..9d63708fe 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -777,8 +777,15 @@ namespace detail { const V depth = cellCentroidsZToEigen(grid_); const V pdepth = subset(depth, well_cells); std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + // Surface density. - std::vector 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 surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases); + // Gravity double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); @@ -786,7 +793,7 @@ namespace detail { std::vector cd = WellDensitySegmented::computeConnectionDensities( wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens); + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // 3. Compute pressure deltas std::vector cdp = diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp index ae43c9532..8844fd292 100644 --- a/opm/autodiff/WellDensitySegmented.cpp +++ b/opm/autodiff/WellDensitySegmented.cpp @@ -36,7 +36,7 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens) + const std::vector& surf_dens_perf) { // Verify that we have consistent input. const int np = wells.number_of_phases; @@ -45,8 +45,8 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells, if (wells.number_of_phases != phase_usage.num_phases) { OPM_THROW(std::logic_error, "Inconsistent input: wells vs. phase_usage."); } - if (surf_dens.size() != size_t(wells.number_of_phases)) { - OPM_THROW(std::logic_error, "Inconsistent input: surf_dens vs. phase_usage."); + if (nperf*np != int(surf_dens_perf.size())) { + OPM_THROW(std::logic_error, "Inconsistent input: wells vs. surf_dens."); } if (nperf*np != int(wstate.perfPhaseRates().size())) { OPM_THROW(std::logic_error, "Inconsistent input: wells vs. wstate."); @@ -94,6 +94,7 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells, const int oilpos = phase_usage.phase_pos[BlackoilPhases::Liquid]; std::vector mix(np); std::vector x(np); + std::vector surf_dens(np); std::vector dens(nperf); for (int w = 0; w < nw; ++w) { for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) { @@ -130,6 +131,10 @@ Opm::WellDensitySegmented::computeConnectionDensities(const Wells& wells, for (int phase = 0; phase < np; ++phase) { volrat += x[phase] / b_perf[perf*np + phase]; } + for (int phase = 0; phase < np; ++phase) { + surf_dens[phase] = surf_dens_perf[perf*np + phase]; + } + // Compute segment density. dens[perf] = std::inner_product(surf_dens.begin(), surf_dens.end(), mix.begin(), 0.0) / volrat; } diff --git a/opm/autodiff/WellDensitySegmented.hpp b/opm/autodiff/WellDensitySegmented.hpp index cf476d6ac..162faa628 100644 --- a/opm/autodiff/WellDensitySegmented.hpp +++ b/opm/autodiff/WellDensitySegmented.hpp @@ -47,14 +47,14 @@ namespace Opm /// \param[in] b_perf inverse ('little b') formation volume factor, size NP, P values per perforation /// \param[in] rsmax_perf saturation point for rs (gas in oil) at each perforation, size N /// \param[in] rvmax_perf saturation point for rv (oil in gas) at each perforation, size N - /// \param[in] surf_dens surface densities for active components, size P + /// \param[in] surf_dens surface densities for active components, size NP, P values per perforation static std::vector computeConnectionDensities(const Wells& wells, const WellStateFullyImplicitBlackoil& wstate, const PhaseUsage& phase_usage, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens); + const std::vector& surf_dens_perf); /// Compute pressure deltas. /// Notation: N = number of perforations, P = number of phases. diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp index c3f79b512..5d58ab05a 100644 --- a/tests/test_welldensitysegmented.cpp +++ b/tests/test_welldensitysegmented.cpp @@ -88,7 +88,16 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas) const std::vector rsmax_perf = { 50, 50, 50, 50, 50, 50, 50, 50, 50, 50 }; const std::vector rvmax_perf = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 }; const std::vector z_perf = { 10, 30, 50, 70, 90, 10, 30, 50, 70, 90 }; - const std::vector surf_dens = { 1000.0, 800.0, 10.0 }; + const std::vector surf_dens = { 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0, + 1000.0, 800.0, 10.0}; const double gravity = Opm::unit::gravity; std::vector cd =