From aaca907f77530fb5fc4b1b8c51bb8517f01e275e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 19 Oct 2020 23:06:13 +0200 Subject: [PATCH 1/7] Productivity Index Calculator: Add Reinitialization operation This commit adds a new member function WellProdIndexCalculator::reInit(const Well& well) which reinitializes the internal arrays in the same way as the constructor. This is needed to ensure that the PI calculation device is synchronised in the case of CTF rescaling-e.g., as a result of WELPI. --- .../wells/WellProdIndexCalculator.cpp | 5 + .../wells/WellProdIndexCalculator.hpp | 10 + tests/test_wellprodindexcalculator.cpp | 200 +++++++++++++++++- 3 files changed, 213 insertions(+), 2 deletions(-) diff --git a/opm/simulators/wells/WellProdIndexCalculator.cpp b/opm/simulators/wells/WellProdIndexCalculator.cpp index 2a96a7c9b..f515f0cc3 100644 --- a/opm/simulators/wells/WellProdIndexCalculator.cpp +++ b/opm/simulators/wells/WellProdIndexCalculator.cpp @@ -96,6 +96,11 @@ Opm::WellProdIndexCalculator::WellProdIndexCalculator(const Well& well) : standardConnFactors_{ calculateStandardConnFactors(well) } {} +void Opm::WellProdIndexCalculator::reInit(const Well& well) +{ + this->standardConnFactors_ = calculateStandardConnFactors(well); +} + double Opm::WellProdIndexCalculator:: connectionProdIndStandard(const std::size_t connIdx, diff --git a/opm/simulators/wells/WellProdIndexCalculator.hpp b/opm/simulators/wells/WellProdIndexCalculator.hpp index a3b0d5a09..671f17bce 100644 --- a/opm/simulators/wells/WellProdIndexCalculator.hpp +++ b/opm/simulators/wells/WellProdIndexCalculator.hpp @@ -41,6 +41,16 @@ namespace Opm { /// per-connection static data. explicit WellProdIndexCalculator(const Well& well); + /// Reinitialization operation + /// + /// Needed to repopulate the internal data members in case of + /// changes to the Well's properties, e.g., as a result of the + /// Well's CTFs being rescaled due to WELPI. + /// + /// \param[in] well Individual well for which to collect + /// per-connection static data. + void reInit(const Well& well); + /// Compute connection-level steady-state productivity index value /// using dynamic phase mobility. /// diff --git a/tests/test_wellprodindexcalculator.cpp b/tests/test_wellprodindexcalculator.cpp index 9953f0c68..4720b683f 100644 --- a/tests/test_wellprodindexcalculator.cpp +++ b/tests/test_wellprodindexcalculator.cpp @@ -30,17 +30,23 @@ #include #include #include +#include #include #include #include #include +#include namespace { + double liquid_PI_unit() + { + return Opm::UnitSystem::newMETRIC().to_si(Opm::UnitSystem::measure::liquid_productivity_index, 1.0); + } + double cp_rm3_per_db() { - return Opm::prefix::centi*Opm::unit::Poise * Opm::unit::cubic(Opm::unit::meter) - / (Opm::unit::day * Opm::unit::barsa); + return Opm::UnitSystem::newMETRIC().to_si(Opm::UnitSystem::measure::transmissibility, 1.0); } std::string drainRadDefaulted() @@ -568,3 +574,193 @@ BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF) } BOOST_AUTO_TEST_SUITE_END() // WellLevel + +// =========================================================================== + +BOOST_AUTO_TEST_SUITE(Re_Init_Connection_Level) + +BOOST_AUTO_TEST_CASE(allDefaulted_SameCF) +{ + auto well = createWell(drainRadDefaulted(), noSkinFactor_SameCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 2.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(allDefaulted_DifferentCF) +{ + auto well = createWell(drainRadDefaulted(), noSkinFactor_DifferentCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 2.0), expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 1.0), expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 0.5), expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(defaultedDRad_Skin2_SameCF) +{ + auto well = createWell(drainRadDefaulted(), skin2_SameCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 2.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(defaultedDRad_skin421_DifferentCF) +{ + auto well = createWell(drainRadDefaulted(), skin421_DifferentCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 2.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 1.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 0.5), 1.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(logarithmic_SameCF) +{ + auto well = createWell(explicitDrainRad(), noSkinFactor_SameCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.5 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 2.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(logarithmic_DifferentCF) +{ + auto well = createWell(explicitDrainRad(), noSkinFactor_DifferentCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.25 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 4.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(logarithmic_Skin2_SameCF) +{ + auto well = createWell(explicitDrainRad(), skin2_SameCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), 0.75 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.5 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), 3.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_CASE(logarithmic_skin421_DifferentCF) +{ + auto well = createWell(explicitDrainRad(), skin421_DifferentCF()); + auto wpiCalc = Opm::WellProdIndexCalculator { well }; + + well.updateWellProductivityIndex(2.0); + const auto scalingFactor = well.getWellPIScalingFactor(1.0*liquid_PI_unit()); + + BOOST_CHECK_CLOSE(scalingFactor, 2.0, 1.0e-10); + + std::vector scalingApplicable; + well.applyWellProdIndexScaling(scalingFactor, scalingApplicable); + + wpiCalc.reInit(well); + + BOOST_REQUIRE_EQUAL(wpiCalc.numConnections(), std::size_t{3}); + + const auto expectCF = 200*cp_rm3_per_db(); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(0, 1.0), (5.0 / 6.0) * 0.5 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(1, 2.0), 1.5 * 1.0 * expectCF, 1.0e-10); + BOOST_CHECK_CLOSE(wpiCalc.connectionProdIndStandard(2, 4.0), (8.0 / 3.0) * 2.0 * expectCF, 1.0e-10); +} + +BOOST_AUTO_TEST_SUITE_END() // Re_Init_Connection_Level From c800c5376dd05aa5bcf497604d24b7b98c92c4b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 19 Oct 2020 23:58:54 +0200 Subject: [PATCH 2/7] Add Special-Purpose Operation to Reset WellState CTFs This commit adds a new member function WellState::resetConnectionTransFactors which overwrites the transmissibility factor of 'well_perf_data_' pertaining to a particular well. This is to keep the values in sync following a rescaling operation such as WELPI. --- opm/simulators/wells/WellState.hpp | 57 +++++++++++++++++-- .../wells/WellStateFullyImplicitBlackoil.hpp | 1 + 2 files changed, 54 insertions(+), 4 deletions(-) diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 7380a7ad9..108fc71e1 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -27,12 +27,13 @@ #include #include -#include -#include -#include -#include #include #include +#include +#include +#include +#include +#include namespace Opm { @@ -96,6 +97,54 @@ namespace Opm } } + /// Special purpose method to support dynamically rescaling a well's + /// CTFs through WELPI. + /// + /// \param[in] well_index Process-local linear index of single well. + /// Must be in the range 0..numWells()-1. + /// + /// \param[in] well_perf_data New perforation data. Only + /// PerforationData::connection_transmissibility_factor actually + /// used (overwrites existing internal values). + void resetConnectionTransFactors(const int well_index, + const std::vector& well_perf_data) + { + if (this->well_perf_data_[well_index].size() != well_perf_data.size()) { + throw std::invalid_argument { + "Size mismatch for perforation data in well " + + std::to_string(well_index) + }; + } + + auto connID = std::size_t{0}; + auto dst = this->well_perf_data_[well_index].begin(); + for (const auto& src : well_perf_data) { + if (dst->cell_index != src.cell_index) { + throw std::invalid_argument { + "Cell index mismatch in connection " + + std::to_string(connID) + + " of well " + + std::to_string(well_index) + }; + } + + if (dst->satnum_id != src.satnum_id) { + throw std::invalid_argument { + "Saturation function table mismatch in connection " + + std::to_string(connID) + + " of well " + + std::to_string(well_index) + }; + } + + dst->connection_transmissibility_factor = + src.connection_transmissibility_factor; + + ++dst; + ++connID; + } + } + /// One bhp pressure per well. std::vector& bhp() { return bhp_; } const std::vector& bhp() const { return bhp_; } diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 10c0ec25f..b658a40dc 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -64,6 +64,7 @@ namespace Opm using BaseType :: wellMap; using BaseType :: numWells; using BaseType :: numPhases; + using BaseType :: resetConnectionTransFactors; /// Allocate and initialize if wells is non-null. Also tries /// to give useful initial values to the bhp(), wellRates() From 4e9e60a71b39e314065e3cec4c24adf568f8ab40 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 20 Oct 2020 00:10:27 +0200 Subject: [PATCH 3/7] Refactor Well Guide Rate Summary Assignment In anticipation of adding more steps later. --- opm/simulators/wells/BlackoilWellModel.hpp | 12 ++-------- .../wells/BlackoilWellModel_impl.hpp | 23 +++++++++++++++++++ 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 43b079a0c..80713e911 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -214,16 +214,7 @@ namespace Opm { { auto wsrpt = well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid())); - for (const auto& well : this->wells_ecl_) { - auto xwPos = wsrpt.find(well.name()); - if (xwPos == wsrpt.end()) { // No well results. Unexpected. - continue; - } - - auto& grval = xwPos->second.guide_rates; - grval.clear(); - grval += this->getGuideRateValues(well); - } + this->assignWellGuideRates(wsrpt); return wsrpt; } @@ -473,6 +464,7 @@ namespace Opm { void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent); + void assignWellGuideRates(data::Wells& wsrpt) const; void assignGroupValues(const int reportStepIdx, const Schedule& sched, std::map& gvalues) const; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 32cda46d1..79c99e14e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -2591,6 +2591,29 @@ namespace Opm { } } + + + + + template + void + BlackoilWellModel:: + assignWellGuideRates(data::Wells& wsrpt) const + { + for (const auto& well : this->wells_ecl_) { + auto xwPos = wsrpt.find(well.name()); + if (xwPos == wsrpt.end()) { // No well results. Unexpected. + continue; + } + + xwPos->second.guide_rates = this->getGuideRateValues(well); + } + } + + + + + template void BlackoilWellModel:: From e780d107abdcf209048477247586faaadcc9ff7c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 20 Oct 2020 00:16:52 +0200 Subject: [PATCH 4/7] Split Well Stat Initialization Out to Helper Function Mostly to reduce the complexity of the implementation of beginReportStep() and to enable easier reordering of the stages. --- opm/simulators/wells/BlackoilWellModel.hpp | 3 + .../wells/BlackoilWellModel_impl.hpp | 87 +++++++++++-------- 2 files changed, 52 insertions(+), 38 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 80713e911..ad4167ad5 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -283,6 +283,9 @@ namespace Opm { void initializeWellProdIndCalculators(); void initializeWellPerfData(); + void initializeWellState(const int timeStepIdx, + const int globalNumWells, + const SummaryState& summaryState); // create the well container std::vector createWellContainer(const int time_step); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 79c99e14e..ceef0b4aa 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -239,50 +239,16 @@ namespace Opm { wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells); this->initializeWellProdIndCalculators(); - initializeWellPerfData(); + // The well state initialize bhp with the cell pressure in the top cell. + // We must therefore provide it with updated cell pressures + this->initializeWellPerfData(); + this->initializeWellState(timeStepIdx, globalNumWells, summaryState); // Wells are active if they are active wells on at least // one process. wells_active_ = localWellsActive() ? 1 : 0; wells_active_ = grid.comm().max(wells_active_); - // The well state initialize bhp with the cell pressure in the top cell. - // We must therefore provide it with updated cell pressures - size_t nc = local_num_cells_; - std::vector cellPressures(nc, 0.0); - ElementContext elemCtx(ebosSimulator_); - const auto& gridView = ebosSimulator_.vanguard().gridView(); - const auto& elemEndIt = gridView.template end(); - for (auto elemIt = gridView.template begin(); - elemIt != elemEndIt; - ++elemIt) - { - const auto& elem = *elemIt; - if (elem.partitionType() != Dune::InteriorEntity) { - continue; - } - elemCtx.updatePrimaryStencil(elem); - elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); - - const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& fs = intQuants.fluidState(); - // copy of get perfpressure in Standard well - // exept for value - double perf_pressure = 0.0; - if (Indices::oilEnabled) { - perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); - } else { - if (Indices::waterEnabled) { - perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); - } else { - perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); - } - } - cellPressures[cellIdx] = perf_pressure; - } - well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_, well_perf_data_, summaryState, globalNumWells); - // handling MS well related if (param_.use_multisegment_well_&& anyMSWellOpenLocal()) { // if we use MultisegmentWell model well_state_.initWellStateMSWell(wells_ecl_, phase_usage_, &previous_well_state_); @@ -656,6 +622,51 @@ namespace Opm { + template + void + BlackoilWellModel:: + initializeWellState(const int timeStepIdx, + const int globalNumWells, + const SummaryState& summaryState) + { + std::vector cellPressures(this->local_num_cells_, 0.0); + ElementContext elemCtx(ebosSimulator_); + + const auto& gridView = ebosSimulator_.vanguard().gridView(); + const auto& elemEndIt = gridView.template end(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + if (elemIt->partitionType() != Dune::InteriorEntity) { + continue; + } + + elemCtx.updatePrimaryStencil(*elemIt); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + + const auto& fs = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0).fluidState(); + + // copy of get perfpressure in Standard well except for value + double& perf_pressure = cellPressures[elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0)]; + if (Indices::oilEnabled) { + perf_pressure = fs.pressure(FluidSystem::oilPhaseIdx).value(); + } else if (Indices::waterEnabled) { + perf_pressure = fs.pressure(FluidSystem::waterPhaseIdx).value(); + } else { + perf_pressure = fs.pressure(FluidSystem::gasPhaseIdx).value(); + } + } + + well_state_.init(cellPressures, schedule(), wells_ecl_, timeStepIdx, + &previous_well_state_, phase_usage_, well_perf_data_, + summaryState, globalNumWells); + } + + + + + template std::vector::WellInterfacePtr > BlackoilWellModel:: From f6130861df6354b5f2d7666fb3e3e5103e88a50c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 29 Oct 2020 23:16:31 +0100 Subject: [PATCH 5/7] Well Data: Report CTFs for Shut Connections Needed to ensure that the CTFs reported to the summary files will be correct in the presence of WELPI even for shut connections. --- opm/simulators/wells/BlackoilWellModel.hpp | 2 ++ .../wells/BlackoilWellModel_impl.hpp | 33 ++++++++++++++++++- 2 files changed, 34 insertions(+), 1 deletion(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index ad4167ad5..0e7e11bf0 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -215,6 +215,7 @@ namespace Opm { auto wsrpt = well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid())); this->assignWellGuideRates(wsrpt); + this->assignShutConnections(wsrpt); return wsrpt; } @@ -468,6 +469,7 @@ namespace Opm { void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent); void assignWellGuideRates(data::Wells& wsrpt) const; + void assignShutConnections(data::Wells& wsrpt) const; void assignGroupValues(const int reportStepIdx, const Schedule& sched, std::map& gvalues) const; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index ceef0b4aa..65f72a3ab 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -23,8 +23,8 @@ #include #include -#include #include +#include #include namespace Opm { @@ -2625,6 +2625,37 @@ namespace Opm { + template + void + BlackoilWellModel:: + assignShutConnections(data::Wells& wsrpt) const + { + for (const auto& well : this->wells_ecl_) { + auto xwPos = wsrpt.find(well.name()); + if (xwPos == wsrpt.end()) { // No well results. Unexpected. + continue; + } + + auto& xcon = xwPos->second.connections; + for (const auto& conn : well.getConnections()) { + if (conn.state() != Connection::State::SHUT) { + continue; + } + + auto& xc = xcon.emplace_back(); + xc.index = conn.global_index(); + xc.pressure = xc.reservoir_rate = 0.0; + + xc.effective_Kh = conn.Kh(); + xc.trans_factor = conn.CF(); + } + } + } + + + + + template void BlackoilWellModel:: From 319c240336cb246eb7c89a0ff408166a2d7ab7c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 29 Oct 2020 23:30:09 +0100 Subject: [PATCH 6/7] Run WELPI Scaling At Beginning of Report Step This commit implements the WELPI feature. We calculate new PI/II values for all wells in the event of a WELPI request and use those values for well-specific WELPI request, to calculate CTF scaling factors. We then apply those factors to all subsequent editions of the well provided the connection factors are eligible for WELPI-based rescaling. If we trigger a rescaling event we also reset the WellState's internal copies of the CTFs and reinitialize the Well PI calculators to ensure the rescaling takes effect immediately. Since we rely on PI values being available at the end of each time step we must also take care to forward those values from the WellState of one report step to the WellState of the next report step. Finally, take care not to redo a WELPI scaling if we've already performed the scaling operation and restart a report step. This, in turn, happens if WELPI is requested on the first report step. --- opm/simulators/wells/BlackoilWellModel.hpp | 6 + .../wells/BlackoilWellModel_impl.hpp | 130 +++++++++++++++++- .../wells/WellStateFullyImplicitBlackoil.hpp | 10 ++ 3 files changed, 144 insertions(+), 2 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 0e7e11bf0..f31b2ba6d 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -322,6 +323,9 @@ namespace Opm { bool report_step_starts_; bool glift_debug = false; bool alternative_well_rate_init_; + + std::optional last_run_wellpi_{}; + std::unique_ptr rateConverter_; std::unique_ptr> vfp_properties_; @@ -468,6 +472,8 @@ namespace Opm { void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent); + void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger); + void assignWellGuideRates(data::Wells& wsrpt) const; void assignShutConnections(data::Wells& wsrpt) const; void assignGroupValues(const int reportStepIdx, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 65f72a3ab..1ff34b7ec 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -23,6 +23,8 @@ #include #include +#include + #include #include #include @@ -238,7 +240,6 @@ namespace Opm { // Make wells_ecl_ contain only this partition's non-shut wells. wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells); - this->initializeWellProdIndCalculators(); // The well state initialize bhp with the cell pressure in the top cell. // We must therefore provide it with updated cell pressures this->initializeWellPerfData(); @@ -287,9 +288,13 @@ namespace Opm { schedule().getVFPInjTables(timeStepIdx), schedule().getVFPProdTables(timeStepIdx)) ); + this->initializeWellProdIndCalculators(); + if (this->schedule().getEvents().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX, timeStepIdx)) { + this->runWellPIScaling(timeStepIdx, local_deferredLogger); + } + // update the previous well state. This is used to restart failed steps. previous_well_state_ = well_state_; - } @@ -2550,6 +2555,127 @@ namespace Opm { + + + template + void + BlackoilWellModel:: + runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger) + { + if (this->last_run_wellpi_.has_value() && (*this->last_run_wellpi_ == timeStepIdx)) { + // We've already run WELPI scaling for this report step. Most + // common for the very first report step. Don't redo WELPI scaling. + return; + } + + auto hasWellPIEvent = [this, timeStepIdx](const int well_index) -> bool + { + return this->schedule() + .hasWellGroupEvent(this->wells_ecl_[well_index].name(), + ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX, + timeStepIdx); + }; + + auto getWellPI = [this](const int well_index) -> double + { + const auto& pu = this->phase_usage_; + const auto np = this->numPhases(); + const auto* pi = &this->well_state_.productivityIndex()[np*well_index + 0]; + + const auto preferred = this->wells_ecl_[well_index].getPreferredPhase(); + switch (preferred) { // Should really have LIQUID = OIL + WATER here too... + case Phase::WATER: + return pu.phase_used[BlackoilPhases::PhaseIndex::Aqua] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Aqua]] + : 0.0; + + case Phase::OIL: + return pu.phase_used[BlackoilPhases::PhaseIndex::Liquid] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Liquid]] + : 0.0; + + case Phase::GAS: + return pu.phase_used[BlackoilPhases::PhaseIndex::Vapour] + ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Vapour]] + : 0.0; + + default: + throw std::invalid_argument { + "Unsupported preferred phase " + + std::to_string(static_cast(preferred)) + }; + } + }; + + auto getWellPIScalingFactor = [this](const int well_index, + const double newWellPI) -> double + { + return this->wells_ecl_[well_index].getWellPIScalingFactor(newWellPI); + }; + + auto rescaleWellPI = + [this, timeStepIdx](const int well_index, + const double scalingFactor) -> void + { + { + const auto& wname = this->wells_ecl_[well_index].name(); + auto& schedule = this->ebosSimulator_.vanguard().schedule(); // Mutable + + schedule.applyWellProdIndexScaling(wname, timeStepIdx, scalingFactor); + this->wells_ecl_[well_index] = schedule.getWell(wname, timeStepIdx); + } + + const auto& well = this->wells_ecl_[well_index]; + auto& pd = this->well_perf_data_[well_index]; + auto pdIter = pd.begin(); + for (const auto& conn : well.getConnections()) { + if (conn.state() != Connection::State::SHUT) { + pdIter->connection_transmissibility_factor = conn.CF(); + ++pdIter; + } + } + + this->well_state_.resetConnectionTransFactors(well_index, pd); + this->prod_index_calc_[well_index].reInit(well); + }; + + // Minimal well setup to compute PI/II values + { + auto saved_previous_well_state = this->previous_well_state_; + this->previous_well_state_ = this->well_state_; + + well_container_ = createWellContainer(timeStepIdx); + + for (auto& well : well_container_) { + well->init(&phase_usage_, depth_, gravity_, local_num_cells_); + } + + std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false); + for (auto& well : well_container_) { + well->updatePerforatedCell(is_cell_perforated_); + } + + this->calculateProductivityIndexValues(local_deferredLogger); + + this->previous_well_state_ = std::move(saved_previous_well_state); + } + + const auto nw = this->numLocalWells(); + for (auto wellID = 0*nw; wellID < nw; ++wellID) { + if (hasWellPIEvent(wellID)) { + const auto newWellPI = getWellPI(wellID); + const auto scalingFactor = getWellPIScalingFactor(wellID, newWellPI); + rescaleWellPI(wellID, scalingFactor); + } + } + + this->last_run_wellpi_ = timeStepIdx; + } + + + + + template void BlackoilWellModel:: diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index b658a40dc..194e95b08 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -270,6 +270,16 @@ namespace Opm } } } + + // Productivity index. + { + auto* thisWellPI = &this ->productivityIndex()[newIndex*np + 0]; + const auto* thatWellPI = &prevState->productivityIndex()[oldIndex*np + 0]; + + for (int p = 0; p < np; ++p) { + thisWellPI[p] = thatWellPI[p]; + } + } } // If in the new step, there is no THP related target/limit anymore, its thp value should be From d06b7f9ec2544b41ff368d94398e353c6ac1053e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 30 Oct 2020 19:13:16 +0100 Subject: [PATCH 7/7] Record Basic Support for WELPI Remove 'WELPI' keyword from MissingFeatures. --- opm/simulators/flow/MissingFeatures.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/simulators/flow/MissingFeatures.cpp b/opm/simulators/flow/MissingFeatures.cpp index d17650d7d..19b530bf3 100644 --- a/opm/simulators/flow/MissingFeatures.cpp +++ b/opm/simulators/flow/MissingFeatures.cpp @@ -791,7 +791,6 @@ namespace MissingFeatures { "WELEVNT", "WELMOVEL" "WELOPENL" - "WELPI", "WELPRI", "WELSOMIN", "WELSPECL",