diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 2b863da81..00af3c86e 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -260,6 +260,7 @@ namespace Opm { V isRs_; V isRv_; V isSg_; + V well_perforation_densities_; //Density of each well perforation V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation. LinearisedBlackoilResidual residual_; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index d6bab10fb..c94420f24 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -657,7 +657,19 @@ namespace detail { } } - + namespace detail { + double getGravity(const double* g, const int dim) { + double grav = 0.0; + if (g) { + // Guard against gravity in anything but last dimension. + for (int dd = 0; dd < dim - 1; ++dd) { + assert(g[dd] == 0.0); + } + grav = g[dim - 1]; + } + return grav; + } + } @@ -731,22 +743,21 @@ namespace detail { // Surface density. std::vector surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases); // Gravity - double grav = 0.0; - const double* g = geo_.gravity(); - const int dim = dimensions(grid_); - if (g) { - // Guard against gravity in anything but last dimension. - for (int dd = 0; dd < dim - 1; ++dd) { - assert(g[dd] == 0.0); - } - grav = g[dim - 1]; - } + double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - // 2. Compute pressure deltas, and store the results. - std::vector cdp = WellDensitySegmented - ::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, perf_depth, - surf_dens, grav); + // 2. Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid_.phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens); + + // 3. Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), perf_depth, cd, grav); + + // 4. Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf); well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); } @@ -1166,6 +1177,25 @@ namespace detail { } // namespace detail + namespace detail { + double computeHydrostaticCorrection(const Wells& wells, const int w, const double vfp_ref_depth, + const ADB::V& well_perforation_densities, const double gravity) { + //For the initial iteration, we have no perforation densities. + if (well_perforation_densities.size() > w) { + const double well_ref_depth = wells.depth_ref[w]; + const double dh = vfp_ref_depth - well_ref_depth; + const int perf = wells.well_connpos[w]; + const double rho = well_perforation_densities[perf]; + const double dp = rho*gravity*dh; + + return dp; + } + else { + return 0.0; + } + } + + } //Namespace @@ -1218,6 +1248,9 @@ namespace detail { current = xw.currentControls()[w]; } + //Get gravity for THP hydrostatic corrrection + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + // Updating well state and primary variables. // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE const double target = well_controls_iget_target(wc, current); @@ -1248,11 +1281,24 @@ namespace detail { //Set *BHP* target by calculating bhp from THP const WellType& well_type = wells().type[w]; - if (well_type == PRODUCER) { - xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq); + + //Gather variables for hydrostatic reconstruction + + if (well_type == INJECTOR) { + double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(); + double dp = detail::computeHydrostaticCorrection( + wells(), w, vfp_ref_depth, + well_perforation_densities_, gravity); + + xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } - else if (well_type == INJECTOR) { - xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp); + else if (well_type == PRODUCER) { + double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(); + double dp = detail::computeHydrostaticCorrection( + wells(), w, vfp_ref_depth, + well_perforation_densities_, gravity); + + xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } else { OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); @@ -1386,10 +1432,14 @@ namespace detail { //THP calculation variables std::vector inj_table_id(nw, -1); std::vector prod_table_id(nw, -1); - ADB::V thp_inj_v = ADB::V::Zero(nw); - ADB::V thp_prod_v = ADB::V::Zero(nw); + ADB::V thp_inj_target_v = ADB::V::Zero(nw); + ADB::V thp_prod_target_v = ADB::V::Zero(nw); ADB::V alq_v = ADB::V::Zero(nw); + //Hydrostatic correction variables + ADB::V rho_v = ADB::V::Zero(nw); + ADB::V vfp_ref_depth_v = ADB::V::Zero(nw); + //Target vars ADB::V bhp_targets = ADB::V::Zero(nw); ADB::V rate_targets = ADB::V::Zero(nw); @@ -1422,19 +1472,29 @@ namespace detail { case THP: { + const int perf = wells().well_connpos[w]; + rho_v[w] = well_perforation_densities_[perf]; + + const int table_id = well_controls_iget_vfp(wc, current); + const double target = well_controls_iget_target(wc, current); + const WellType& well_type = wells().type[w]; if (well_type == INJECTOR) { - inj_table_id[w] = well_controls_iget_vfp(wc, current); - thp_inj_v[w] = well_controls_iget_target(wc, current); + inj_table_id[w] = table_id; + thp_inj_target_v[w] = target; alq_v[w] = -1e100; + vfp_ref_depth_v[w] = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(); + thp_inj_elems.push_back(w); } else if (well_type == PRODUCER) { - prod_table_id[w] = well_controls_iget_vfp(wc, current); - thp_prod_v[w] = well_controls_iget_target(wc, current); + prod_table_id[w] = table_id; + thp_prod_target_v[w] = target; alq_v[w] = well_controls_iget_alq(wc, current); + vfp_ref_depth_v[w] = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(); + thp_prod_elems.push_back(w); } else { @@ -1468,23 +1528,33 @@ namespace detail { } //Calculate BHP target from THP - const ADB thp_inj = ADB::constant(thp_inj_v); - const ADB thp_prod = ADB::constant(thp_prod_v); + const ADB thp_inj_target = ADB::constant(thp_inj_target_v); + const ADB thp_prod_target = ADB::constant(thp_prod_target_v); const ADB alq = ADB::constant(alq_v); - const ADB thp_inj_targets = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj); - const ADB thp_prod_targets = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod, alq); + const ADB bhp_from_thp_inj = vfp_properties_->getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target); + const ADB bhp_from_thp_prod = vfp_properties_->getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq); + + //Perform hydrostatic correction to computed targets + const ADB well_ref_depth = ADB::constant(ADB::V::Map(wells().depth_ref, nw)); + const ADB vfp_ref_depth = ADB::constant(vfp_ref_depth_v); + const ADB dh = vfp_ref_depth - well_ref_depth; + const ADB rho = ADB::constant(rho_v); + double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const ADB dp = rho*gravity*dh; + const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); + const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); //Calculate residuals - const ADB thp_inj_residual = state.bhp - thp_inj_targets; - const ADB thp_prod_residual = state.bhp - thp_prod_targets; + const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj; + const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod; const ADB bhp_residual = state.bhp - bhp_targets; const ADB rate_residual = rate_distr * state.qs - rate_targets; //Select the right residual for each well - residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, bhp_residual.size()) + - superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, thp_inj_residual.size()) + - superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, thp_prod_residual.size()) + - superset(subset(rate_residual, rate_elems), rate_elems, rate_residual.size()); + residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) + + superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) + + superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) + + superset(subset(rate_residual, rate_elems), rate_elems, nw); // For wells that are dead (not flowing), and therefore not communicating // with the reservoir, we set the equation to be equal to the well's total diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index ef8a5e3ba..8c28aa7c5 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -94,7 +94,6 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, - VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, const ADB& aqua, const ADB& liquid, @@ -118,7 +117,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, //Get the table for each well std::vector well_tables(nw, NULL); for (int i=0; i= 0) { + if (table_id[i] > 0) { well_tables[i] = detail::getTable(m_tables, table_id[i]); } } @@ -226,4 +225,14 @@ double VFPInjProperties::thp(int table_id, +const VFPInjTable* VFPInjProperties::getTable(const int table_id) const { + return detail::getTable(m_tables, table_id); +} + + + + + + + } //Namespace Opm diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index 4d7f31fc8..42bbe8fba 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -126,6 +126,12 @@ public: const double& vapour, const double& bhp) const; + /** + * Returns the table associated with the ID, or throws an exception if + * the table does not exist + */ + const VFPInjTable* getTable(const int table_id) const; + private: // Map which connects the table number with the table itself std::map m_tables; diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index ecdefd548..d95718b3d 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -235,6 +235,14 @@ double VFPProdProperties::thp(int table_id, +const VFPProdTable* VFPProdProperties::getTable(const int table_id) const { + return detail::getTable(m_tables, table_id); +} + + + + + } diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index e86aa1abc..877b51de3 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -137,6 +137,11 @@ public: const double& bhp, const double& alq) const; + /** + * Returns the table associated with the ID, or throws an exception if + * the table does not exist + */ + const VFPProdTable* getTable(const int table_id) const; private: // Map which connects the table number with the table itself diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp index 669c260c8..576fb8583 100644 --- a/opm/autodiff/WellDensitySegmented.cpp +++ b/opm/autodiff/WellDensitySegmented.cpp @@ -26,16 +26,17 @@ #include + + + std::vector -Opm::WellDensitySegmented::computeConnectionPressureDelta(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& z_perf, - const std::vector& surf_dens, - const double gravity) +Opm::WellDensitySegmented::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) { // Verify that we have consistent input. const int np = wells.number_of_phases; @@ -53,9 +54,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, if (nperf*np != int(b_perf.size())) { OPM_THROW(std::logic_error, "Inconsistent input: wells vs. b_perf."); } - if (nperf != int(z_perf.size())) { - OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf."); - } if ((!rsmax_perf.empty()) || (!rvmax_perf.empty())) { // Need both oil and gas phases. if (!phase_usage.phase_used[BlackoilPhases::Liquid]) { @@ -66,14 +64,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, } } - // Algorithm: - - // We'll assume the perforations are given in order from top to - // bottom for each well. By top and bottom we do not necessarily - // mean in a geometric sense (depth), but in a topological sense: - // the 'top' perforation is nearest to the surface topologically. - // Our goal is to compute a pressure delta for each perforation. - // 1. Compute the flow (in surface volume units for each // component) exiting up the wellbore from each perforation, // taking into account flow from lower in the well, and @@ -145,7 +135,38 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, } } - // 3. Compute pressure differences between perforations. + return dens; +} + + + + + +std::vector +Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, + const std::vector& z_perf, + const std::vector& dens_perf, + const double gravity) { + const int np = wells.number_of_phases; + const int nw = wells.number_of_wells; + const int nperf = wells.well_connpos[nw]; + + if (nperf != int(z_perf.size())) { + OPM_THROW(std::logic_error, "Inconsistent input: wells vs. z_perf."); + } + if (nperf != int(dens_perf.size())) { + OPM_THROW(std::logic_error, "Inconsistent input: wells vs. dens_perf."); + } + + // Algorithm: + + // We'll assume the perforations are given in order from top to + // bottom for each well. By top and bottom we do not necessarily + // mean in a geometric sense (depth), but in a topological sense: + // the 'top' perforation is nearest to the surface topologically. + // Our goal is to compute a pressure delta for each perforation. + + // 1. Compute pressure differences between perforations. // dp_perf will contain the pressure difference between a // perforation and the one above it, except for the first // perforation for each well, for which it will be the @@ -155,11 +176,12 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w+1]; ++perf) { const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1]; const double dz = z_perf[perf] - z_above; - dp_perf[perf] = dz * dens[perf] * gravity; + dp_perf[perf] = dz * dens_perf[perf] * gravity; + //dens[wells.well_connpos[w]] } } - // 4. Compute pressure differences to the reference point (bhp) by + // 2. Compute pressure differences to the reference point (bhp) by // accumulating the already computed adjacent pressure // differences, storing the result in dp_perf. // This accumulation must be done per well. diff --git a/opm/autodiff/WellDensitySegmented.hpp b/opm/autodiff/WellDensitySegmented.hpp index 82a4494b8..cf476d6ac 100644 --- a/opm/autodiff/WellDensitySegmented.hpp +++ b/opm/autodiff/WellDensitySegmented.hpp @@ -39,7 +39,7 @@ namespace Opm class WellDensitySegmented { public: - /// Compute pressure deltas. + /// Compute well segment densities /// Notation: N = number of perforations, P = number of phases. /// \param[in] wells struct with static well info /// \param[in] wstate dynamic well solution information, only perfRates() is used @@ -47,17 +47,24 @@ 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] z_perf depth values for each perforation, size N /// \param[in] surf_dens surface densities for active components, size P + 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); + + /// Compute pressure deltas. + /// Notation: N = number of perforations, P = number of phases. + /// \param[in] wells struct with static well info + /// \param[in] z_perf depth values for each perforation, size N + /// \param[in] dens_perf densities for each perforation, size N (typically computed using computeConnectionDensities) /// \param[in] gravity gravity acceleration constant static std::vector computeConnectionPressureDelta(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& z_perf, - const std::vector& surf_dens, + const std::vector& dens_perf, const double gravity); }; diff --git a/tests/test_welldensitysegmented.cpp b/tests/test_welldensitysegmented.cpp index 0ab4e38a5..c3f79b512 100644 --- a/tests/test_welldensitysegmented.cpp +++ b/tests/test_welldensitysegmented.cpp @@ -90,7 +90,15 @@ BOOST_AUTO_TEST_CASE(TestPressureDeltas) 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 double gravity = Opm::unit::gravity; - const std::vector dp = WellDensitySegmented::computeConnectionPressureDelta(*wells, wellstate, pu, b_perf, rsmax_perf, rvmax_perf, z_perf, surf_dens, gravity); + + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + *wells, wellstate, pu, + b_perf, rsmax_perf, rvmax_perf, surf_dens); + std::vector dp = + WellDensitySegmented::computeConnectionPressureDelta( + *wells, z_perf, cd, gravity); + const std::vector answer = { 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity, 20e3*gravity, 62e3*gravity, 106e3*gravity, 152e3*gravity, 200e3*gravity }; BOOST_REQUIRE_EQUAL(dp.size(), answer.size());