Merge pull request #2866 from bska/activate-welpi-scaling

Initial Implementation of WELPI Feature
This commit is contained in:
Kai Bao 2020-12-04 20:44:42 +01:00 committed by GitHub
commit 7ecb2e4813
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 522 additions and 58 deletions

View File

@ -791,7 +791,6 @@ namespace MissingFeatures {
"WELEVNT",
"WELMOVEL"
"WELOPENL"
"WELPI",
"WELPRI",
"WELSOMIN",
"WELSPECL",

View File

@ -34,6 +34,7 @@
#include <functional>
#include <map>
#include <memory>
#include <optional>
#include <string>
#include <tuple>
#include <unordered_map>
@ -215,16 +216,8 @@ 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);
this->assignShutConnections(wsrpt);
return wsrpt;
}
@ -296,6 +289,9 @@ namespace Opm {
void initializeWellProdIndCalculators();
void initializeWellPerfData();
void initializeWellState(const int timeStepIdx,
const int globalNumWells,
const SummaryState& summaryState);
// create the well container
std::vector<WellInterfacePtr > createWellContainer(const int time_step);
@ -331,6 +327,9 @@ namespace Opm {
bool report_step_starts_;
bool glift_debug = false;
bool alternative_well_rate_init_;
std::optional<int> last_run_wellpi_{};
std::unique_ptr<RateConverterType> rateConverter_;
std::unique_ptr<VFPProperties<VFPInjProperties,VFPProdProperties>> vfp_properties_;
@ -482,6 +481,10 @@ 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,
const Schedule& sched,
std::map<std::string, data::GroupData>& gvalues) const;

View File

@ -23,8 +23,10 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <utility>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <algorithm>
#include <utility>
#include <fmt/format.h>
namespace Opm {
@ -261,51 +263,16 @@ namespace Opm {
wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells);
local_parallel_well_info_ = createLocalParallelWellInfo(wells_ecl_);
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<double> cellPressures(nc, 0.0);
ElementContext elemCtx(ebosSimulator_);
const auto& gridView = ebosSimulator_.vanguard().gridView();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
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_);
@ -344,9 +311,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_;
}
@ -679,6 +650,51 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
initializeWellState(const int timeStepIdx,
const int globalNumWells,
const SummaryState& summaryState)
{
std::vector<double> cellPressures(this->local_num_cells_, 0.0);
ElementContext elemCtx(ebosSimulator_);
const auto& gridView = ebosSimulator_.vanguard().gridView();
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
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<typename TypeTag>
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>::
@ -2563,6 +2579,127 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
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<int>(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<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
@ -2615,6 +2752,60 @@ namespace Opm {
}
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
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 <typename TypeTag>
void
BlackoilWellModel<TypeTag>::
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 <typename TypeTag>
void
BlackoilWellModel<TypeTag>::

View File

@ -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,

View File

@ -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.
///

View File

@ -27,12 +27,13 @@
#include <opm/simulators/wells/PerforationData.hpp>
#include <array>
#include <map>
#include <memory>
#include <string>
#include <vector>
#include <cassert>
#include <cstddef>
#include <map>
#include <memory>
#include <stdexcept>
#include <string>
#include <vector>
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<PerforationData>& 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<double>& bhp() { return bhp_; }
const std::vector<double>& bhp() const { return bhp_; }

View File

@ -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()
@ -269,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

View File

@ -30,17 +30,23 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/parser/eclipse/Units/Units.hpp>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
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<bool> 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<bool> 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<bool> 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<bool> 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<bool> 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<bool> 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<bool> 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<bool> 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