OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5354 from akva2/various_well_helpers_template_scalar

Browse Source 
Various well helpers: template Scalar type
This commit is contained in:
Bård Skaflestad
2024-05-21 13:14:50 +02:00
committed by GitHub
parent fd4ac412a6 8e000d08ce
commit 8a4e78e7c4
18 changed files with 686 additions and 590 deletions
Download Patch File Download Diff File Expand all files Collapse all files

59
opm/simulators/flow/ActionHandler.cpp
View File

@@ -112,12 +112,14 @@ namespace {
namespace Opm {
ActionHandler::ActionHandler(EclipseState& ecl_state,
Schedule& schedule,
Action::State& actionState,
SummaryState& summaryState,
BlackoilWellModelGeneric<double>& wellModel,
Parallel::Communication comm)
template<class Scalar>
ActionHandler<Scalar>::
ActionHandler(EclipseState& ecl_state,
Schedule& schedule,
Action::State& actionState,
SummaryState& summaryState,
BlackoilWellModelGeneric<Scalar>& wellModel,
Parallel::Communication comm)
: ecl_state_(ecl_state)
, schedule_(schedule)
, actionState_(actionState)
@@ -126,9 +128,11 @@ ActionHandler::ActionHandler(EclipseState& ecl_state,
, comm_(comm)
{}
void ActionHandler::applyActions(const int reportStep,
const double sim_time,
const TransFunc& transUp)
template<class Scalar>
void ActionHandler<Scalar>::
applyActions(const int reportStep,
const double sim_time,
const TransFunc& transUp)
{
OPM_TIMEBLOCK(applyActions);
const auto& actions = schedule_[reportStep].actions();
@@ -184,10 +188,12 @@ void ActionHandler::applyActions(const int reportStep,
}
}
void ActionHandler::applySimulatorUpdate(const int report_step,
const SimulatorUpdate& sim_update,
bool& commit_wellstate,
const TransFunc& updateTrans)
template<class Scalar>
void ActionHandler<Scalar>::
applySimulatorUpdate(const int report_step,
const SimulatorUpdate& sim_update,
bool& commit_wellstate,
const TransFunc& updateTrans)
{
OPM_TIMEBLOCK(applySimulatorUpdate);
@@ -207,12 +213,13 @@ void ActionHandler::applySimulatorUpdate(const int report_step,
}
}
std::unordered_map<std::string, double>
ActionHandler::fetchWellPI(const int reportStep,
const Action::ActionX& action,
const std::vector<std::string>& matching_wells) const
template<class Scalar>
std::unordered_map<std::string, Scalar>
ActionHandler<Scalar>::
fetchWellPI(const int reportStep,
const Action::ActionX& action,
const std::vector<std::string>& matching_wells) const
{
auto wellpi_wells = action.wellpi_wells(WellMatcher(schedule_[reportStep].well_order(),
schedule_[reportStep].wlist_manager()),
matching_wells);
@@ -221,7 +228,7 @@ ActionHandler::fetchWellPI(const int reportStep,
return {};
const auto num_wells = schedule_[reportStep].well_order().size();
std::vector<double> wellpi_vector(num_wells);
std::vector<Scalar> wellpi_vector(num_wells);
for (const auto& wname : wellpi_wells) {
if (this->wellModel_.hasWell(wname)) {
const auto& well = schedule_.getWell( wname, reportStep );
@@ -230,10 +237,10 @@ ActionHandler::fetchWellPI(const int reportStep,
}
if (comm_.size() > 1) {
std::vector<double> wellpi_buffer(num_wells * comm_.size());
std::vector<Scalar> wellpi_buffer(num_wells * comm_.size());
comm_.gather( wellpi_vector.data(), wellpi_buffer.data(), num_wells, 0 );
if (comm_.rank() == 0) {
for (int rank=1; rank < comm_.size(); rank++) {
for (int rank = 1; rank < comm_.size(); rank++) {
for (std::size_t well_index=0; well_index < num_wells; well_index++) {
const auto global_index = rank*num_wells + well_index;
const auto value = wellpi_buffer[global_index];
@@ -245,7 +252,7 @@ ActionHandler::fetchWellPI(const int reportStep,
comm_.broadcast(wellpi_vector.data(), wellpi_vector.size(), 0);
}
std::unordered_map<std::string, double> wellpi;
std::unordered_map<std::string, Scalar> wellpi;
for (const auto& wname : wellpi_wells) {
const auto& well = schedule_.getWell( wname, reportStep );
wellpi[wname] = wellpi_vector[ well.seqIndex() ];
@@ -253,8 +260,10 @@ ActionHandler::fetchWellPI(const int reportStep,
return wellpi;
}
void ActionHandler::evalUDQAssignments(const unsigned episodeIdx,
UDQState& udq_state)
template<class Scalar>
void ActionHandler<Scalar>::
evalUDQAssignments(const unsigned episodeIdx,
UDQState& udq_state)
{
const auto& udq = schedule_[episodeIdx].udq();
@@ -266,4 +275,6 @@ void ActionHandler::evalUDQAssignments(const unsigned episodeIdx,
udq_state);
}
template class ActionHandler<double>;
} // namespace Opm

7
opm/simulators/flow/ActionHandler.hpp
View File

@@ -46,6 +46,7 @@ class SummaryState;
class UDQState;
//! \brief Class handling Action support in simulator
template<class Scalar>
class ActionHandler
{
public:
@@ -56,7 +57,7 @@ public:
Schedule& schedule,
Action::State& actionState,
SummaryState& summaryState,
BlackoilWellModelGeneric<double>& wellModel,
BlackoilWellModelGeneric<Scalar>& wellModel,
Parallel::Communication comm);
void applyActions(int reportStep,
@@ -79,7 +80,7 @@ public:
bool& commit_wellstate,
const TransFunc& updateTrans);
std::unordered_map<std::string, double>
std::unordered_map<std::string, Scalar>
fetchWellPI(int reportStep,
const Action::ActionX& action,
const std::vector<std::string>& matching_wells) const;
@@ -88,7 +89,7 @@ public:
Schedule& schedule_;
Action::State& actionState_;
SummaryState& summaryState_;
BlackoilWellModelGeneric<double>& wellModel_;
BlackoilWellModelGeneric<Scalar>& wellModel_;
Parallel::Communication comm_;
};

2
opm/simulators/flow/FlowProblem.hpp
View File

@@ -2800,7 +2800,7 @@ private:
PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
TracerModel tracerModel_;
ActionHandler actionHandler_;
ActionHandler<Scalar> actionHandler_;
template<class T>
struct BCData

71
opm/simulators/wells/WellAssemble.cpp
View File

@@ -37,7 +37,6 @@
#include <opm/simulators/wells/WellState.hpp>
#include <cassert>
#include <cmath>
#include <stdexcept>
namespace Opm
@@ -66,7 +65,7 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
{
const auto current = well_state.well(well_.indexOfWell()).production_cmode;
const auto& pu = well_.phaseUsage();
const double efficiencyFactor = well_.wellEcl().getEfficiencyFactor();
const Scalar efficiencyFactor = well_.wellEcl().getEfficiencyFactor();
switch (current) {
case Well::ProducerCMode::ORAT: {
@@ -102,7 +101,7 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
case Well::ProducerCMode::RESV: {
auto total_rate = rates[0]; // To get the correct type only.
total_rate = 0.0;
std::vector<double> convert_coeff(well_.numPhases(), 1.0);
std::vector<Scalar> convert_coeff(well_.numPhases(), 1.0);
well_.rateConverter().calcCoeff(/*fipreg*/ 0, well_.pvtRegionIdx(), well_state.well(well_.indexOfWell()).surface_rates, convert_coeff);
for (int phase = 0; phase < 3; ++phase) {
if (pu.phase_used[phase]) {
@@ -113,7 +112,7 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
if (controls.prediction_mode) {
control_eq = total_rate - controls.resv_rate;
} else {
std::vector<double> hrates(well_.numPhases(), 0.);
std::vector<Scalar> hrates(well_.numPhases(), 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
hrates[pu.phase_pos[Water]] = controls.water_rate;
}
@@ -123,9 +122,9 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
hrates[pu.phase_pos[Gas]] = controls.gas_rate;
}
std::vector<double> hrates_resv(well_.numPhases(), 0.);
std::vector<Scalar> hrates_resv(well_.numPhases(), 0.);
well_.rateConverter().calcReservoirVoidageRates(/*fipreg*/ 0, well_.pvtRegionIdx(), hrates, hrates_resv);
double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
Scalar target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
control_eq = total_rate - target;
}
break;
@@ -153,7 +152,9 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
active_rates[pu.phase_pos[canonical_phase]] = rates[canonical_phase];
}
}
auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
auto rCoeff = [this, &group_state](const RegionId id, const int region,
const std::optional<std::string>& prod_gname,
std::vector<Scalar>& coeff)
{
if (prod_gname)
well_.rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
@@ -161,15 +162,17 @@ assembleControlEqProd(const WellState<Scalar>& well_state,
well_.rateConverter().calcCoeff(id, region, coeff);
};
WellGroupControls(well_).getGroupProductionControl(group, well_state,
group_state,
schedule,
summaryState,
bhp, active_rates,
rCoeff,
efficiencyFactor,
control_eq,
deferred_logger);
WellGroupControls(well_).getGroupProductionControl(group,
well_state,
group_state,
schedule,
summaryState,
bhp,
active_rates,
rCoeff,
efficiencyFactor,
control_eq,
deferred_logger);
break;
}
case Well::ProducerCMode::CMODE_UNDEFINED: {
@@ -203,7 +206,7 @@ assembleControlEqInj(const WellState<Scalar>& well_state,
auto current = well_state.well(well_.indexOfWell()).injection_cmode;
const InjectorType injectorType = controls.injector_type;
const auto& pu = well_.phaseUsage();
const double efficiencyFactor = well_.wellEcl().getEfficiencyFactor();
const Scalar efficiencyFactor = well_.wellEcl().getEfficiencyFactor();
switch (current) {
case Well::InjectorCMode::RATE: {
@@ -211,10 +214,10 @@ assembleControlEqInj(const WellState<Scalar>& well_state,
break;
}
case Well::InjectorCMode::RESV: {
std::vector<double> convert_coeff(well_.numPhases(), 1.0);
std::vector<Scalar> convert_coeff(well_.numPhases(), 1.0);
well_.rateConverter().calcInjCoeff(/*fipreg*/ 0, well_.pvtRegionIdx(), convert_coeff);
double coeff;
Scalar coeff;
switch (injectorType) {
case InjectorType::WATER: {
@@ -247,26 +250,29 @@ assembleControlEqInj(const WellState<Scalar>& well_state,
case Well::InjectorCMode::GRUP: {
assert(well_.wellEcl().isAvailableForGroupControl());
const auto& group = schedule.getGroup(well_.wellEcl().groupName(), well_.currentStep());
auto rCoeff = [this, &group_state](const RegionId id, const int region, const std::optional<std::string>& prod_gname, std::vector<double>& coeff)
auto rCoeff = [this, &group_state](const RegionId id, const int region,
const std::optional<std::string>& prod_gname,
std::vector<Scalar>& coeff)
{
if(prod_gname) {
well_.rateConverter().calcCoeff(id, region, group_state.production_rates(*prod_gname), coeff);
well_.rateConverter().calcCoeff(id, region,
group_state.production_rates(*prod_gname), coeff);
} else {
well_.rateConverter().calcInjCoeff(id, region, coeff);
}
};
WellGroupControls(well_).getGroupInjectionControl(group,
well_state,
group_state,
schedule,
summaryState,
injectorType,
bhp,
injection_rate,
rCoeff,
efficiencyFactor,
control_eq,
deferred_logger);
well_state,
group_state,
schedule,
summaryState,
injectorType,
bhp,
injection_rate,
rCoeff,
efficiencyFactor,
control_eq,
deferred_logger);
break;
}
case Well::InjectorCMode::CMODE_UNDEFINED: {
@@ -318,4 +324,5 @@ INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,8u>)
INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,9u>)
INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,10u>)
INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,11u>)
} // namespace Opm

433
opm/simulators/wells/WellBhpThpCalculator.cpp
View File

@@ -48,11 +48,11 @@
static constexpr bool extraBhpAtThpLimitOutput = false;
static constexpr bool extraThpFromBhpOutput = false;
namespace Opm
{
namespace Opm {
bool
WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) const
template<class Scalar>
bool WellBhpThpCalculator<Scalar>::
wellHasTHPConstraints(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
@@ -70,7 +70,9 @@ WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) co
return false;
}
double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState) const
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
getTHPConstraint(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
@@ -86,7 +88,9 @@ double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState)
return 0.0;
}
double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
{
const auto& well_ecl = well_.wellEcl();
if (well_ecl.isInjector()) {
@@ -102,12 +106,14 @@ double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summ
return 0.0;
}
double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
const std::optional<double>& alq,
const double thp_limit,
DeferredLogger& deferred_logger) const
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
calculateThpFromBhp(const std::vector<Scalar>& rates,
const Scalar bhp,
const Scalar rho,
const std::optional<Scalar>& alq,
const Scalar thp_limit,
DeferredLogger& deferred_logger) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@@ -115,31 +121,31 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
const double aqua = rates[Water];
const double liquid = rates[Oil];
const double vapour = rates[Gas];
const Scalar aqua = rates[Water];
const Scalar liquid = rates[Oil];
const Scalar vapour = rates[Gas];
// pick the density in the top layer
double thp = 0.0;
Scalar thp = 0.0;
const int table_id = well_.wellEcl().vfp_table_number();
if (well_.isInjector()) {
assert(!alq.has_value());
const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
const Scalar vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
auto thp_func =
[this, table_id, aqua, liquid, vapour, dp](
const double bhp_value, const double pressure_loss) {
const Scalar bhp_value, const Scalar pressure_loss) {
return this->well_.vfpProperties()->getInj()->thp(
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss);
};
thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger);
}
else if (well_.isProducer()) {
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
auto thp_func =
[this, table_id, aqua, liquid, vapour, dp, &alq]
(const double bhp_value, const double pressure_loss) {
(const Scalar bhp_value, const Scalar pressure_loss) {
return this->well_.vfpProperties()->getProd()->thp(
table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value());
};
@@ -151,36 +157,35 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector<double>& rate
return thp;
}
double
WellBhpThpCalculator::
findThpFromBhpIteratively(
const std::function<double(const double, const double)>& thp_func,
const double bhp,
const double thp_limit,
const double dp,
DeferredLogger& deferred_logger) const
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
const Scalar bhp,
const Scalar thp_limit,
const Scalar dp,
DeferredLogger& deferred_logger) const
{
auto pressure_loss = getVfpBhpAdjustment(bhp + dp, thp_limit);
auto thp = thp_func(bhp, pressure_loss);
const double tolerance = 1e-5 * unit::barsa;
const Scalar tolerance = 1e-5 * unit::barsa;
bool do_iterate = true;
int it = 1;
int max_iterations = 50;
while(do_iterate) {
while (do_iterate) {
if (it > max_iterations) {
break;
}
double thp_prev = thp;
Scalar thp_prev = thp;
pressure_loss = getVfpBhpAdjustment(bhp + dp - pressure_loss, thp_prev);
thp = thp_func(bhp, pressure_loss);
auto error = std::fabs(thp-thp_prev);
auto error = std::fabs(thp - thp_prev);
if (extraThpFromBhpOutput) {
const std::string msg = fmt::format(
"findThpFromBhpIteratively(): iteration {}, thp = {}, bhp = {}, "
"pressure_loss = {}, error = {}", it, thp, bhp+dp-pressure_loss, pressure_loss, error);
deferred_logger.debug(msg);
}
if (std::fabs(thp-thp_prev) < tolerance) {
if (std::fabs(thp - thp_prev) < tolerance) {
break;
}
it++;
@@ -188,14 +193,15 @@ findThpFromBhpIteratively(
return thp;
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
template<class Scalar>
std::optional<Scalar>
WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitProd(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
const double thp_limit,
const Scalar maxPerfPress,
const Scalar rho,
const Scalar alq_value,
const Scalar thp_limit,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
@@ -221,34 +227,43 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
// Make the fbhp() function.
const auto& controls = well_.wellEcl().productionControls(summary_state);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
const Scalar vfp_ref_depth = table.getDatumDepth();
const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
vfp_ref_depth,
rho,
well_.gravity());
auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<double>& rates) {
auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector<Scalar>& rates) {
assert(rates.size() == 3);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number,
well_.indexOfWell());
const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number,
well_.indexOfWell());
const bool use_vfpexp = well_.useVfpExplicit();
const double bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
rates[Water],
rates[Oil],
rates[Gas],
thp_limit,
alq_value,
wfr,
gfr,
use_vfpexp);
const Scalar bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number,
rates[Water],
rates[Oil],
rates[Gas],
thp_limit,
alq_value,
wfr,
gfr,
use_vfpexp);
return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit);
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
auto flo = [&table](const std::vector<Scalar>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Find the bhp-point where production becomes nonzero.
auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); };
auto bhp_max = this->bhpMax(fflo, controls.bhp_limit, maxPerfPress, table.getFloAxis().front(), deferred_logger);
auto fflo = [&flo, &frates](Scalar bhp) { return flo(frates(bhp)); };
auto bhp_max = this->bhpMax(fflo,
controls.bhp_limit,
maxPerfPress,
table.getFloAxis().front(),
deferred_logger);
// could not solve for the bhp-point, we could not continue to find the bhp
if (!bhp_max.has_value()) {
@@ -257,16 +272,17 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
"find bhp-point where production becomes non-zero for well " + well_.name());
return std::nullopt;
}
const std::array<double, 2> range {controls.bhp_limit, *bhp_max};
const std::array<Scalar, 2> range {controls.bhp_limit, *bhp_max};
return this->computeBhpAtThpLimit(frates, fbhp, range, deferred_logger);
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
template<class Scalar>
std::optional<Scalar>
WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const Scalar rho,
const Scalar flo_rel_tol,
const int max_iteration,
const bool throwOnError,
DeferredLogger& deferred_logger) const
@@ -282,13 +298,15 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
}
}
void WellBhpThpCalculator::updateThp(const double rho,
const bool stop_or_zero_rate_target,
const std::function<double()>& alq_value,
const std::array<unsigned,3>& active,
WellState<double>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
template<class Scalar>
void WellBhpThpCalculator<Scalar>::
updateThp(const Scalar rho,
const bool stop_or_zero_rate_target,
const std::function<Scalar()>& alq_value,
const std::array<unsigned,3>& active,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
@@ -308,7 +326,7 @@ void WellBhpThpCalculator::updateThp(const double rho,
}
// the well is under other control types, we calculate the thp based on bhp and rates
std::vector<double> rates(3, 0.0);
std::vector<Scalar> rates(3, 0.0);
const PhaseUsage& pu = well_.phaseUsage();
if (active[Water]) {
@@ -320,18 +338,19 @@ void WellBhpThpCalculator::updateThp(const double rho,
if (active[Gas]) {
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
const std::optional<double> alq = this->well_.isProducer() ? std::optional<double>(alq_value()) : std::nullopt;
const double thp_limit = well_.getTHPConstraint(summary_state);
const std::optional<Scalar> alq = this->well_.isProducer() ? std::optional<Scalar>(alq_value()) : std::nullopt;
const Scalar thp_limit = well_.getTHPConstraint(summary_state);
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, thp_limit, deferred_logger);
}
template<class Scalar>
template<class EvalWell>
EvalWell WellBhpThpCalculator::
calculateBhpFromThp(const WellState<double>& well_state,
EvalWell WellBhpThpCalculator<Scalar>::
calculateBhpFromThp(const WellState<Scalar>& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
const Scalar rho,
DeferredLogger& deferred_logger) const
{
// TODO: when well is under THP control, the BHP is dependent on the rates,
@@ -349,8 +368,8 @@ calculateBhpFromThp(const WellState<double>& well_state,
const EvalWell aqua = rates[Water];
const EvalWell liquid = rates[Oil];
const EvalWell vapour = rates[Gas];
const double thp_limit = well_.getTHPConstraint(summaryState);
double vfp_ref_depth;
const Scalar thp_limit = well_.getTHPConstraint(summaryState);
Scalar vfp_ref_depth;
EvalWell bhp_tab;
if (well_.isInjector() )
{
@@ -375,57 +394,62 @@ calculateBhpFromThp(const WellState<double>& well_state,
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger);
}
double bhp_tab_double_value;
if constexpr (std::is_same_v<EvalWell, double>) {
Scalar bhp_tab_double_value;
if constexpr (std::is_same_v<EvalWell, Scalar>) {
bhp_tab_double_value = bhp_tab;
}
else { // EvalWell and bhp_tab is of type DenseAd::Evaluation<double,...,...>
bhp_tab_double_value = bhp_tab.value();
}
const auto bhp_adjustment = getVfpBhpAdjustment(bhp_tab_double_value, thp_limit);
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(
well_.refDepth(), vfp_ref_depth, rho, well_.gravity());
const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
vfp_ref_depth,
rho,
well_.gravity());
return bhp_tab - dp_hydro + bhp_adjustment;
}
double
WellBhpThpCalculator::
calculateMinimumBhpFromThp(const WellState<double>& well_state,
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
calculateMinimumBhpFromThp(const WellState<Scalar>& well_state,
const Well& well,
const SummaryState& summaryState,
const double rho) const
const Scalar rho) const
{
assert(well_.isProducer()); // only producers can go here for now
const double thp_limit = well_.getTHPConstraint(summaryState);
const Scalar thp_limit = well_.getTHPConstraint(summaryState);
const auto& controls = well.productionControls(summaryState);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
const double bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number,
const Scalar bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number,
thp_limit, wfr, gfr,
well_.getALQ(well_state));
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit);
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity());
return bhp_min - dp_hydro + bhp_adjustment;
}
double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const
template<class Scalar>
Scalar WellBhpThpCalculator<Scalar>::
getVfpBhpAdjustment(const Scalar bhp_tab, const Scalar thp_limit) const
{
return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit);
}
template<class Scalar>
template<class ErrorPolicy>
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates,
std::optional<Scalar>
WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimitInjImpl(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const Scalar rho,
const Scalar flo_rel_tol,
const int max_iteration,
DeferredLogger& deferred_logger) const
{
@@ -475,12 +499,12 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// Make the fbhp() function.
const auto& controls = well_.wellEcl().injectionControls(summary_state);
const auto& table = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = this->getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
const Scalar vfp_ref_depth = table.getDatumDepth();
const Scalar thp_limit = this->getTHPConstraint(summary_state);
const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(),
vfp_ref_depth,
rho, well_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<Scalar>& rates) {
assert(rates.size() == 3);
const auto bhp = well_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit);
@@ -488,25 +512,25 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
auto flo = [&table](const std::vector<Scalar>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary.
std::vector<double> flo_samples = table.getFloAxis();
std::vector<Scalar> flo_samples = table.getFloAxis();
if (flo_samples[0] > 0.0) {
const double f0 = flo_samples[0];
const Scalar f0 = flo_samples[0];
flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 });
}
const double flo_bhp_limit = flo(frates(controls.bhp_limit));
const Scalar flo_bhp_limit = flo(frates(controls.bhp_limit));
if (flo_samples.back() < flo_bhp_limit) {
flo_samples.push_back(flo_bhp_limit);
}
// Find bhp values for inflow relation corresponding to flo samples.
std::vector<double> bhp_samples;
for (double flo_sample : flo_samples) {
std::vector<Scalar> bhp_samples;
for (Scalar flo_sample : flo_samples) {
if (flo_sample > flo_bhp_limit) {
// We would have to go over the bhp limit to obtain a
// flow of this magnitude. We associate all such flows
@@ -516,16 +540,16 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
bhp_samples.push_back(controls.bhp_limit);
break;
}
auto eq = [&flo, &frates, flo_sample](double bhp) {
auto eq = [&flo, &frates, flo_sample](Scalar bhp) {
return flo(frates(bhp)) - flo_sample;
};
// TODO: replace hardcoded low/high limits.
const double low = 10.0 * unit::barsa;
const double high = 800.0 * unit::barsa;
const double flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back());
const Scalar low = 10.0 * unit::barsa;
const Scalar high = 800.0 * unit::barsa;
const Scalar flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
const Scalar solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp);
}
@@ -539,7 +563,7 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// Find bhp values for VFP relation corresponding to flo samples.
const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size()
std::vector<double> fbhp_samples(num_samples);
std::vector<Scalar> fbhp_samples(num_samples);
for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
}
@@ -564,8 +588,8 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
// We only look at the valid
int sign_change_index = -1;
for (int ii = 0; ii < num_samples - 1; ++ii) {
const double curr = fbhp_samples[ii] - bhp_samples[ii];
const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
const Scalar curr = fbhp_samples[ii] - bhp_samples[ii];
const Scalar next = fbhp_samples[ii + 1] - bhp_samples[ii + 1];
if (curr * next < 0.0) {
// Sign change in the [ii, ii + 1] interval.
sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution.
@@ -578,13 +602,13 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
}
// Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double bhp) {
auto eq = [&fbhp, &frates](Scalar bhp) {
return fbhp(frates(bhp)) - bhp;
};
// TODO: replace hardcoded low/high limits.
const double low = bhp_samples[sign_change_index + 1];
const double high = bhp_samples[sign_change_index];
const double bhp_tolerance = 0.01 * unit::barsa;
const Scalar low = bhp_samples[sign_change_index + 1];
const Scalar high = bhp_samples[sign_change_index];
const Scalar bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
if (low == high) {
// We are in the high flow regime where the bhp_samples
@@ -595,7 +619,7 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
return std::nullopt;
}
try {
const double solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
const Scalar solved_bhp = RegulaFalsiBisection<ErrorPolicy>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
if constexpr (extraBhpAtThpLimitOutput) {
OpmLog::debug("***** " + well_.name() + " solved_bhp = " + std::to_string(solved_bhp)
@@ -610,20 +634,21 @@ computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double
}
}
std::optional<double>
WellBhpThpCalculator::
bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
template<class Scalar>
std::optional<Scalar>
WellBhpThpCalculator<Scalar>::
bhpMax(const std::function<Scalar(const Scalar)>& fflo,
const Scalar bhp_limit,
const Scalar maxPerfPress,
const Scalar vfp_flo_front,
DeferredLogger& deferred_logger) const
{
// Find the bhp-point where production becomes nonzero.
double bhp_max = 0.0;
double low = bhp_limit;
double high = maxPerfPress + 1.0 * unit::barsa;
double f_low = fflo(low);
double f_high = fflo(high);
Scalar bhp_max = 0.0;
Scalar low = bhp_limit;
Scalar high = maxPerfPress + 1.0 * unit::barsa;
Scalar f_low = fflo(low);
Scalar f_high = fflo(high);
if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) +
@@ -633,7 +658,7 @@ bhpMax(const std::function<double(const double)>& fflo,
}
int adjustments = 0;
const int max_adjustments = 10;
const double adjust_amount = 5.0 * unit::barsa;
const Scalar adjust_amount = 5.0 * unit::barsa;
while (f_low * f_high > 0.0 && adjustments < max_adjustments) {
// Same sign, adjust high to see if we can flip it.
high += adjust_amount;
@@ -656,12 +681,12 @@ bhpMax(const std::function<double(const double)>& fflo,
} else {
// Bisect to find a bhp point where we produce, but
// not a large amount ('eps' below).
const double eps = 0.1 * std::fabs(vfp_flo_front);
const Scalar eps = 0.1 * std::fabs(vfp_flo_front);
const int maxit = 50;
int it = 0;
while (std::fabs(f_low) > eps && it < maxit) {
const double curr = 0.5*(low + high);
const double f_curr = fflo(curr);
const Scalar curr = 0.5*(low + high);
const Scalar f_curr = fflo(curr);
if (f_curr * f_low > 0.0) {
low = curr;
f_low = f_curr;
@@ -690,11 +715,12 @@ bhpMax(const std::function<double(const double)>& fflo,
return bhp_max;
}
std::optional<double>
WellBhpThpCalculator::
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
template<class Scalar>
std::optional<Scalar>
WellBhpThpCalculator<Scalar>::
computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const std::function<Scalar(const std::vector<Scalar>)>& fbhp,
const std::array<Scalar, 2>& range,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
@@ -714,15 +740,16 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
// highest rate) should be returned.
// Define the equation we want to solve.
auto eq = [&fbhp, &frates](double bhp) {
auto eq = [&fbhp, &frates](Scalar bhp) {
return fbhp(frates(bhp)) - bhp;
};
// Find appropriate brackets for the solution.
std::optional<double> approximate_solution;
double low, high;
std::optional<Scalar> approximate_solution;
Scalar low, high;
// trying to use bisect way to locate a bracket
bool finding_bracket = this->bisectBracket(eq, range, low, high, approximate_solution, deferred_logger);
bool finding_bracket = this->bisectBracket(eq, range, low, high,
approximate_solution, deferred_logger);
// based on the origional design, if an approximate solution is suggested, we use this value directly
// in the long run, we might change it
@@ -744,10 +771,10 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
// Solve for the proper solution in the given interval.
const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
const Scalar bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ThrowOnError>::
const Scalar solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp;
}
@@ -758,21 +785,21 @@ computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& fra
}
}
bool
WellBhpThpCalculator::
bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
template<class Scalar>
bool WellBhpThpCalculator<Scalar>::
bisectBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const
{
bool finding_bracket = false;
low = range[0];
high = range[1];
double eq_high = eq(high);
double eq_low = eq(low);
const double eq_bhplimit = eq_low;
Scalar eq_high = eq(high);
Scalar eq_low = eq(low);
const Scalar eq_bhplimit = eq_low;
if constexpr (extraBhpAtThpLimitOutput) {
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() +
" low = " + std::to_string(low) +
@@ -783,12 +810,12 @@ bisectBracket(const std::function<double(const double)>& eq,
if (eq_low * eq_high > 0.0) {
// Failed to bracket the zero.
// If this is due to having two solutions, bisect until bracketed.
double abs_low = std::fabs(eq_low);
double abs_high = std::fabs(eq_high);
Scalar abs_low = std::fabs(eq_low);
Scalar abs_high = std::fabs(eq_high);
int bracket_attempts = 0;
const int max_bracket_attempts = 20;
double interval = high - low;
const double min_interval = 1.0 * unit::barsa;
Scalar interval = high - low;
const Scalar min_interval = 1.0 * unit::barsa;
while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) {
if (abs_high < abs_low) {
low = 0.5 * (low + high);
@@ -815,7 +842,7 @@ bisectBracket(const std::function<double(const double)>& eq,
}
} else { // eq_low * eq_high > 0.0
// Still failed bracketing!
const double limit = 0.1 * unit::barsa;
const Scalar limit = 0.1 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) {
// Return the least bad solution if less off than 0.1 bar.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
@@ -834,20 +861,20 @@ bisectBracket(const std::function<double(const double)>& eq,
return finding_bracket;
}
bool
WellBhpThpCalculator::
bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
template<class Scalar>
bool WellBhpThpCalculator<Scalar>::
bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
DeferredLogger& deferred_logger)
{
bool bracket_found = false;
low = range[0];
high = range[1];
const int sample_number = 200;
const double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high = 0.0;
const Scalar interval = (high - low) / sample_number;
Scalar eq_low = eq(low);
Scalar eq_high = 0.0;
for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i;
eq_high = eq(high);
@@ -866,10 +893,11 @@ bruteForceBracket(const std::function<double(const double)>& eq,
return bracket_found;
}
bool WellBhpThpCalculator::
isStableSolution(const WellState<double>& well_state,
template<class Scalar>
bool WellBhpThpCalculator<Scalar>::
isStableSolution(const WellState<Scalar>& well_state,
const Well& well,
const std::vector<double>& rates,
const std::vector<Scalar>& rates,
const SummaryState& summaryState) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@@ -879,10 +907,10 @@ isStableSolution(const WellState<double>& well_state,
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
const double aqua = rates[Water];
const double liquid = rates[Oil];
const double vapour = rates[Gas];
const double thp = well_.getTHPConstraint(summaryState);
const Scalar aqua = rates[Water];
const Scalar liquid = rates[Oil];
const Scalar vapour = rates[Gas];
const Scalar thp = well_.getTHPConstraint(summaryState);
const auto& controls = well.productionControls(summaryState);
const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
@@ -897,27 +925,28 @@ isStableSolution(const WellState<double>& well_state,
return true;
} else { // maybe check if ipr is available
const auto ipr = getFloIPR(well_state, well, summaryState);
return bhp.dflo + 1/ipr.second >= 0;
return bhp.dflo + 1.0 / ipr.second >= 0;
}
}
std::optional<double> WellBhpThpCalculator::
estimateStableBhp(const WellState<double>& well_state,
template<class Scalar>
std::optional<Scalar> WellBhpThpCalculator<Scalar>::
estimateStableBhp(const WellState<Scalar>& well_state,
const Well& well,
const std::vector<double>& rates,
const double rho,
const std::vector<Scalar>& rates,
const Scalar rho,
const SummaryState& summaryState) const
{
// Given a *converged* well_state with ipr, estimate bhp of the stable solution
const auto& controls = well.productionControls(summaryState);
const double thp = well_.getTHPConstraint(summaryState);
const Scalar thp = well_.getTHPConstraint(summaryState);
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double aqua = rates[BlackoilPhases::Aqua];
const double liquid = rates[BlackoilPhases::Liquid];
const double vapour = rates[BlackoilPhases::Vapour];
double flo = detail::getFlo(table, aqua, liquid, vapour);
double wfr, gfr;
const Scalar aqua = rates[BlackoilPhases::Aqua];
const Scalar liquid = rates[BlackoilPhases::Liquid];
const Scalar vapour = rates[BlackoilPhases::Vapour];
Scalar flo = detail::getFlo(table, aqua, liquid, vapour);
Scalar wfr, gfr;
if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) {
wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell());
gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell());
@@ -928,11 +957,11 @@ estimateStableBhp(const WellState<double>& well_state,
auto ipr = getFloIPR(well_state, well, summaryState);
const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth,
rho, well_.gravity());
auto bhp_adjusted = [this, &thp, &dp_hydro](const double bhp) {
auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) {
return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp);
};
const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted);
@@ -942,10 +971,11 @@ estimateStableBhp(const WellState<double>& well_state,
} else {
return std::nullopt;
}
}
}
std::pair<double, double> WellBhpThpCalculator::
getFloIPR(const WellState<double>& well_state,
template<class Scalar>
std::pair<Scalar, Scalar> WellBhpThpCalculator<Scalar>::
getFloIPR(const WellState<Scalar>& well_state,
const Well& well,
const SummaryState& summary_state) const
{
@@ -954,21 +984,23 @@ getFloIPR(const WellState<double>& well_state,
const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const auto& pu = well_.phaseUsage();
const auto& ipr_a = well_state.well(well_.indexOfWell()).implicit_ipr_a;
const double& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const double& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const double& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
const Scalar& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const Scalar& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const Scalar& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
const auto& ipr_b = well_state.well(well_.indexOfWell()).implicit_ipr_b;
const double& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const double& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const double& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
const Scalar& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0;
const Scalar& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0;
const Scalar& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0;
// The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type),
// but we can equally use it to pick/add the corresponding ipr_a, ipr_b
return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a),
detail::getFlo(table, aqua_b, liquid_b, vapour_b));
}
template class WellBhpThpCalculator<double>;
#define INSTANCE(...) \
template __VA_ARGS__ WellBhpThpCalculator:: \
template __VA_ARGS__ WellBhpThpCalculator<double>:: \
calculateBhpFromThp<__VA_ARGS__>(const WellState<double>&, \
const std::vector<__VA_ARGS__>&, \
const Well&, \
@@ -993,4 +1025,5 @@ INSTANCE(DenseAd::Evaluation<double,-1,8u>)
INSTANCE(DenseAd::Evaluation<double,-1,9u>)
INSTANCE(DenseAd::Evaluation<double,-1,10u>)
INSTANCE(DenseAd::Evaluation<double,-1,11u>)
} // namespace Opm

127
opm/simulators/wells/WellBhpThpCalculator.hpp
View File

@@ -26,11 +26,9 @@
#include <functional>
#include <optional>
#include <string>
#include <vector>
namespace Opm
{
namespace Opm {
class DeferredLogger;
class SummaryState;
@@ -39,136 +37,137 @@ template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing BHP limits.
template<class Scalar>
class WellBhpThpCalculator {
public:
//! \brief Constructor sets reference to well.
WellBhpThpCalculator(const WellInterfaceGeneric<double>& well) : well_(well) {}
WellBhpThpCalculator(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
//! \brief Checks if well has THP constraints.
bool wellHasTHPConstraints(const SummaryState& summaryState) const;
//! \brief Get THP constraint for well.
double getTHPConstraint(const SummaryState& summaryState) const;
Scalar getTHPConstraint(const SummaryState& summaryState) const;
//! \brief Obtain the most strict BHP from BHP limits.
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
Scalar mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
//! \brief Calculates THP from BHP.
double calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
const std::optional<double>& alq,
const double thp_limit,
Scalar calculateThpFromBhp(const std::vector<Scalar>& rates,
const Scalar bhp,
const Scalar rho,
const std::optional<Scalar>& alq,
const Scalar thp_limit,
DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for a producer.
std::optional<double>
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
std::optional<Scalar>
computeBhpAtThpLimitProd(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
const double alq_value,
const double thp_limit,
const Scalar maxPerfPress,
const Scalar rho,
const Scalar alq_value,
const Scalar thp_limit,
DeferredLogger& deferred_logger) const;
//! \brief Compute BHP from THP limit for an injector.
std::optional<double>
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
std::optional<Scalar>
computeBhpAtThpLimitInj(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const Scalar rho,
const Scalar flo_rel_tol,
const int max_iteration,
const bool throwOnError,
DeferredLogger& deferred_logger) const;
//! \brief Update THP.
void updateThp(const double rho,
void updateThp(const Scalar rho,
const bool stop_or_zero_rate_target,
const std::function<double()>& alq_value,
const std::function<Scalar()>& alq_value,
const std::array<unsigned,3>& active,
WellState<double>& well_state,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const;
template<class EvalWell>
EvalWell calculateBhpFromThp(const WellState<double>& well_state,
EvalWell calculateBhpFromThp(const WellState<Scalar>& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
const Scalar rho,
DeferredLogger& deferred_logger) const;
double calculateMinimumBhpFromThp(const WellState<double>& well_state,
const Well& well,
const SummaryState& summaryState,
const double rho) const;
Scalar calculateMinimumBhpFromThp(const WellState<Scalar>& well_state,
const Well& well,
const SummaryState& summaryState,
const Scalar rho) const;
bool isStableSolution(const WellState<double>& well_state,
bool isStableSolution(const WellState<Scalar>& well_state,
const Well& well,
const std::vector<double>& rates,
const std::vector<Scalar>& rates,
const SummaryState& summaryState) const;
std::optional<double>
estimateStableBhp (const WellState<double>& well_state,
std::optional<Scalar>
estimateStableBhp (const WellState<Scalar>& well_state,
const Well& well,
const std::vector<double>& rates,
const double rho,
const std::vector<Scalar>& rates,
const Scalar rho,
const SummaryState& summaryState) const;
std::pair<double, double>
getFloIPR(const WellState<double>& well_state,
std::pair<Scalar, Scalar>
getFloIPR(const WellState<Scalar>& well_state,
const Well& well,
const SummaryState& summary_state) const;
private:
//! \brief Compute BHP from THP limit for an injector - implementation.
template<class ErrorPolicy>
std::optional<double>
computeBhpAtThpLimitInjImpl(const std::function<std::vector<double>(const double)>& frates,
std::optional<Scalar>
computeBhpAtThpLimitInjImpl(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const SummaryState& summary_state,
const double rho,
const double flo_rel_tol,
const Scalar rho,
const Scalar flo_rel_tol,
const int max_iteration,
DeferredLogger& deferred_logger) const;
//! \brief Calculate max BHP.
std::optional<double>
bhpMax(const std::function<double(const double)>& fflo,
const double bhp_limit,
const double maxPerfPress,
const double vfp_flo_front,
std::optional<Scalar>
bhpMax(const std::function<Scalar(const Scalar)>& fflo,
const Scalar bhp_limit,
const Scalar maxPerfPress,
const Scalar vfp_flo_front,
DeferredLogger& deferred_logger) const;
//! \brief Common code for finding BHP from THP limit for producers/injectors.
std::optional<double>
computeBhpAtThpLimit(const std::function<std::vector<double>(const double)>& frates,
const std::function<double(const std::vector<double>)>& fbhp,
const std::array<double, 2>& range,
std::optional<Scalar>
computeBhpAtThpLimit(const std::function<std::vector<Scalar>(const Scalar)>& frates,
const std::function<Scalar(const std::vector<Scalar>)>& fbhp,
const std::array<Scalar, 2>& range,
DeferredLogger& deferred_logger) const;
//! \brief Get pressure adjustment to the bhp calculated from VFP table
double getVfpBhpAdjustment(const double bph_tab, const double thp_limit) const;
Scalar getVfpBhpAdjustment(const Scalar bph_tab, const Scalar thp_limit) const;
//! \brief Find limits using bisection.
bool bisectBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
bool bisectBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const;
//! \brief Find limits using brute-force solver.
static bool bruteForceBracket(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
DeferredLogger& deferred_logger);
double findThpFromBhpIteratively(const std::function<double(const double, const double)>& thp_func,
const double bhp,
const double thp_limit,
const double dp,
Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
const Scalar bhp,
const Scalar thp_limit,
const Scalar dp,
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
};
}

58
opm/simulators/wells/WellConstraints.cpp
View File

@@ -31,11 +31,11 @@
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
namespace Opm
{
namespace Opm {
bool WellConstraints::
checkIndividualConstraints(SingleWellState<double>& ws,
template<class Scalar>
bool WellConstraints<Scalar>::
checkIndividualConstraints(SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
const RateConvFunc& calcReservoirVoidageRates,
bool& thp_limit_violated_but_not_switched,
@@ -69,8 +69,9 @@ checkIndividualConstraints(SingleWellState<double>& ws,
return false;
}
Well::InjectorCMode WellConstraints::
activeInjectionConstraint(const SingleWellState<double>& ws,
template<class Scalar>
Well::InjectorCMode WellConstraints<Scalar>::
activeInjectionConstraint(const SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
bool& thp_limit_violated_but_not_switched,
DeferredLogger& deferred_logger,
@@ -84,7 +85,7 @@ activeInjectionConstraint(const SingleWellState<double>& ws,
if (controls.hasControl(Well::InjectorCMode::BHP) && currentControl != Well::InjectorCMode::BHP)
{
const auto& bhp = controls.bhp_limit;
double current_bhp = ws.bhp;
Scalar current_bhp = ws.bhp;
if (bhp < current_bhp)
return Well::InjectorCMode::BHP;
}
@@ -92,7 +93,7 @@ activeInjectionConstraint(const SingleWellState<double>& ws,
if (controls.hasControl(Well::InjectorCMode::RATE) && currentControl != Well::InjectorCMode::RATE)
{
InjectorType injectorType = controls.injector_type;
double current_rate = 0.0;
Scalar current_rate = 0.0;
switch (injectorType) {
case InjectorType::WATER:
@@ -120,14 +121,14 @@ activeInjectionConstraint(const SingleWellState<double>& ws,
if (controls.hasControl(Well::InjectorCMode::RESV) && currentControl != Well::InjectorCMode::RESV)
{
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
Scalar current_rate = 0.0;
if (pu.phase_used[BlackoilPhases::Aqua])
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
if (pu.phase_used[BlackoilPhases::Liquid])
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
if (pu.phase_used[BlackoilPhases::Vapour])
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.reservoir_rate < current_rate)
@@ -140,7 +141,7 @@ activeInjectionConstraint(const SingleWellState<double>& ws,
if (controls.hasControl(Well::InjectorCMode::THP) && currentControl != Well::InjectorCMode::THP)
{
const auto& thp = well_.getTHPConstraint(summaryState);
double current_thp = ws.thp;
Scalar current_thp = ws.thp;
if (thp < current_thp) {
bool rate_less_than_potential = true;
for (int p = 0; p < well_.numPhases(); ++p) {
@@ -166,8 +167,9 @@ activeInjectionConstraint(const SingleWellState<double>& ws,
return currentControl;
}
Well::ProducerCMode WellConstraints::
activeProductionConstraint(const SingleWellState<double>& ws,
template<class Scalar>
Well::ProducerCMode WellConstraints<Scalar>::
activeProductionConstraint(const SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
const RateConvFunc& calcReservoirVoidageRates,
bool& thp_limit_violated_but_not_switched,
@@ -179,37 +181,37 @@ activeProductionConstraint(const SingleWellState<double>& ws,
const auto currentControl = ws.production_cmode;
if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP) {
const double bhp_limit = controls.bhp_limit;
double current_bhp = ws.bhp;
const Scalar bhp_limit = controls.bhp_limit;
Scalar current_bhp = ws.bhp;
if (bhp_limit > current_bhp)
return Well::ProducerCMode::BHP;
}
if (controls.hasControl(Well::ProducerCMode::ORAT) && currentControl != Well::ProducerCMode::ORAT) {
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
Scalar current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
if (controls.oil_rate < current_rate)
return Well::ProducerCMode::ORAT;
}
if (controls.hasControl(Well::ProducerCMode::WRAT) && currentControl != Well::ProducerCMode::WRAT) {
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
Scalar current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (controls.water_rate < current_rate)
return Well::ProducerCMode::WRAT;
}
if (controls.hasControl(Well::ProducerCMode::GRAT) && currentControl != Well::ProducerCMode::GRAT) {
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]];
Scalar current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]];
if (controls.gas_rate < current_rate)
return Well::ProducerCMode::GRAT;
}
if (controls.hasControl(Well::ProducerCMode::LRAT) && currentControl != Well::ProducerCMode::LRAT) {
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
Scalar current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
bool skip = false;
if (controls.liquid_rate == controls.oil_rate) {
const double current_water_rate = ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
const Scalar current_water_rate = ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (std::abs(current_water_rate) < 1e-12) {
skip = true;
deferred_logger.debug("LRAT_ORAT_WELL", "Well " + well_.name() + " The LRAT target is equal the ORAT target and the water rate is zero, skip checking LRAT");
@@ -220,7 +222,7 @@ activeProductionConstraint(const SingleWellState<double>& ws,
}
if (controls.hasControl(Well::ProducerCMode::RESV) && currentControl != Well::ProducerCMode::RESV) {
double current_rate = 0.0;
Scalar current_rate = 0.0;
if (pu.phase_used[BlackoilPhases::Aqua])
current_rate -= ws.reservoir_rates[pu.phase_pos[BlackoilPhases::Aqua]];
@@ -237,7 +239,7 @@ activeProductionConstraint(const SingleWellState<double>& ws,
const int fipreg = 0; // not considering the region for now
const int np = well_.numPhases();
std::vector<double> surface_rates(np, 0.0);
std::vector<Scalar> surface_rates(np, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua])
surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = controls.water_rate;
if (pu.phase_used[BlackoilPhases::Liquid])
@@ -245,10 +247,10 @@ activeProductionConstraint(const SingleWellState<double>& ws,
if (pu.phase_used[BlackoilPhases::Vapour])
surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.gas_rate;
std::vector<double> voidage_rates(np, 0.0);
std::vector<Scalar> voidage_rates(np, 0.0);
calcReservoirVoidageRates(fipreg, well_.pvtRegionIdx(), surface_rates, voidage_rates);
double resv_rate = 0.0;
Scalar resv_rate = 0.0;
for (int p = 0; p < np; ++p)
resv_rate += voidage_rates[p];
@@ -259,7 +261,7 @@ activeProductionConstraint(const SingleWellState<double>& ws,
if (well_.wellHasTHPConstraints(summaryState) && currentControl != Well::ProducerCMode::THP) {
const auto& thp = well_.getTHPConstraint(summaryState);
double current_thp = ws.thp;
Scalar current_thp = ws.thp;
if (thp > current_thp && !ws.trivial_target) {
// If WVFPEXP item 4 is set to YES1 or YES2
// switching to THP is prevented if the well will
@@ -291,4 +293,6 @@ activeProductionConstraint(const SingleWellState<double>& ws,
return currentControl;
}
template class WellConstraints<double>;
} // namespace Opm

16
opm/simulators/wells/WellConstraints.hpp
View File

@@ -27,7 +27,6 @@
#include <opm/input/eclipse/Schedule/Well/Well.hpp>
#include <functional>
#include <utility>
#include <vector>
#include <optional>
@@ -44,17 +43,18 @@ enum class WellInjectorCMode;
enum class WellProducerCMode;
//! \brief Class for computing well group constraints.
template<class Scalar>
class WellConstraints {
public:
//! \brief Constructor sets reference to well.
WellConstraints(const WellInterfaceGeneric<double>& well) : well_(well) {}
WellConstraints(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId, const int,
const std::vector<double>&,
std::vector<double>&)>;
const std::vector<Scalar>&,
std::vector<Scalar>&)>;
bool
checkIndividualConstraints(SingleWellState<double>& ws,
checkIndividualConstraints(SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
const RateConvFunc& calcReservoirVoidageRates,
bool& thp_limit_violated_but_not_switched,
@@ -64,21 +64,21 @@ public:
private:
WellInjectorCMode
activeInjectionConstraint(const SingleWellState<double>& ws,
activeInjectionConstraint(const SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
bool& thp_limit_violated_but_not_switched,
DeferredLogger& deferred_logger,
const std::optional<Well::InjectionControls>& inj_controls = std::nullopt) const;
WellProducerCMode
activeProductionConstraint(const SingleWellState<double>& ws,
activeProductionConstraint(const SingleWellState<Scalar>& ws,
const SummaryState& summaryState,
const RateConvFunc& calcReservoirVoidageRates,
bool& thp_limit_violated_but_not_switched,
DeferredLogger& deferred_logger,
const std::optional<Well::ProductionControls>& prod_controls = std::nullopt) const;
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
};
}

30
opm/simulators/wells/WellConvergence.cpp
View File

@@ -30,18 +30,18 @@
#include <cmath>
#include <stdexcept>
namespace Opm
{
namespace Opm {
void WellConvergence::
checkConvergenceControlEq(const WellState<double>& well_state,
template<class Scalar>
void WellConvergence<Scalar>::
checkConvergenceControlEq(const WellState<Scalar>& well_state,
const Tolerances& tolerances,
const double well_control_residual,
const Scalar well_control_residual,
const bool well_is_stopped,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
Scalar control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
@@ -120,21 +120,21 @@ checkConvergenceControlEq(const WellState<double>& well_state,
}
}
void
WellConvergence::
checkConvergencePolyMW(const std::vector<double>& res,
template<class Scalar>
void WellConvergence<Scalar>::
checkConvergencePolyMW(const std::vector<Scalar>& res,
const int Bhp,
const double maxResidualAllowed,
const Scalar maxResidualAllowed,
ConvergenceReport& report) const
{
if (well_.isInjector()) {
// checking the convergence of the perforation rates
const double wat_vel_tol = 1.e-8;
const Scalar wat_vel_tol = 1.e-8;
const int dummy_component = -1;
using CR = ConvergenceReport;
const auto wat_vel_failure_type = CR::WellFailure::Type::MassBalance;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
const double wat_vel_residual = res[Bhp + 1 + perf];
const Scalar wat_vel_residual = res[Bhp + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, well_.name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
@@ -145,10 +145,10 @@ checkConvergencePolyMW(const std::vector<double>& res,
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const Scalar pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < well_.numPerfs(); ++perf) {
const double pskin_residual = res[Bhp + 1 + perf + well_.numPerfs()];
const Scalar pskin_residual = res[Bhp + 1 + perf + well_.numPerfs()];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, well_.name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
@@ -160,4 +160,6 @@ checkConvergencePolyMW(const std::vector<double>& res,
}
}
template class WellConvergence<double>;
}

23
opm/simulators/wells/WellConvergence.hpp
View File

@@ -33,36 +33,37 @@ class DeferredLogger;
template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
template<class Scalar>
class WellConvergence
{
public:
WellConvergence(const WellInterfaceGeneric<double>& well)
WellConvergence(const WellInterfaceGeneric<Scalar>& well)
: well_(well)
{}
struct Tolerances {
double bhp; //!< Tolerance for bhp controlled well
double thp; //!< Tolerance for thp controlled well
double rates; //!< Tolerance for a rate controlled well
double grup; //!< Tolerance for a grup controlled well
double max_residual_allowed; //!< Max residual allowd
Scalar bhp; //!< Tolerance for bhp controlled well
Scalar thp; //!< Tolerance for thp controlled well
Scalar rates; //!< Tolerance for a rate controlled well
Scalar grup; //!< Tolerance for a grup controlled well
Scalar max_residual_allowed; //!< Max residual allowd
};
// checking the convergence of the well control equations
void checkConvergenceControlEq(const WellState<double>& well_state,
void checkConvergenceControlEq(const WellState<Scalar>& well_state,
const Tolerances& tolerances,
const double well_control_residual,
const Scalar well_control_residual,
const bool well_is_stopped,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
void checkConvergencePolyMW(const std::vector<double>& res,
void checkConvergencePolyMW(const std::vector<Scalar>& res,
const int Bhp,
const double maxResidualAllowed,
const Scalar maxResidualAllowed,
ConvergenceReport& report) const;
private:
const WellInterfaceGeneric<double>& well_;
const WellInterfaceGeneric<Scalar>& well_;
};
}

47
opm/simulators/wells/WellGroupConstraints.cpp
View File

@@ -33,12 +33,13 @@
namespace Opm
{
std::pair<bool, double>
WellGroupConstraints::
template<class Scalar>
std::pair<bool, Scalar>
WellGroupConstraints<Scalar>::
checkGroupConstraintsInj(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const double efficiencyFactor,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Scalar efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
@@ -69,12 +70,12 @@ checkGroupConstraintsInj(const Group& group,
}
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
std::vector<Scalar> resv_coeff(well_.phaseUsage().num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(well_.indexOfWell());
// Call check for the well's injection phase.
return WellGroupHelpers<double>::checkGroupConstraintsInj(well_.name(),
return WellGroupHelpers<Scalar>::checkGroupConstraintsInj(well_.name(),
well_.wellEcl().groupName(),
group,
well_state,
@@ -91,23 +92,24 @@ checkGroupConstraintsInj(const Group& group,
deferred_logger);
}
std::pair<bool, double>
WellGroupConstraints::
template<class Scalar>
std::pair<bool, Scalar>
WellGroupConstraints<Scalar>::
checkGroupConstraintsProd(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const double efficiencyFactor,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Scalar efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
DeferredLogger& deferred_logger) const
{
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
std::vector<Scalar> resv_coeff(well_.phaseUsage().num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(well_.indexOfWell());
return WellGroupHelpers<double>::checkGroupConstraintsProd(well_.name(),
return WellGroupHelpers<Scalar>::checkGroupConstraintsProd(well_.name(),
well_.wellEcl().groupName(),
group,
well_state,
@@ -123,9 +125,10 @@ checkGroupConstraintsProd(const Group& group,
deferred_logger);
}
bool WellGroupConstraints::
checkGroupConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
template<class Scalar>
bool WellGroupConstraints<Scalar>::
checkGroupConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
@@ -146,8 +149,8 @@ checkGroupConstraints(WellState<double>& well_state,
// check, skipping over only the single group parent whose
// control is the active one for the well (if any).
const auto& group = schedule.getGroup(well.groupName(), well_.currentStep());
const double efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, double> group_constraint =
const Scalar efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, Scalar> group_constraint =
this->checkGroupConstraintsInj(group, well_state,
group_state, efficiencyFactor,
schedule, summaryState,
@@ -177,8 +180,8 @@ checkGroupConstraints(WellState<double>& well_state,
// check, skipping over only the single group parent whose
// control is the active one for the well (if any).
const auto& group = schedule.getGroup(well.groupName(), well_.currentStep());
const double efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, double> group_constraint =
const Scalar efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, Scalar> group_constraint =
this->checkGroupConstraintsProd(group, well_state,
group_state, efficiencyFactor,
schedule, summaryState,
@@ -199,4 +202,6 @@ checkGroupConstraints(WellState<double>& well_state,
return false;
}
template class WellGroupConstraints<double>;
} // namespace Opm

27
opm/simulators/wells/WellGroupConstraints.hpp
View File

@@ -44,45 +44,46 @@ template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing well group constraints.
template<class Scalar>
class WellGroupConstraints {
public:
//! \brief Constructor sets reference to well.
WellGroupConstraints(const WellInterfaceGeneric<double>& well) : well_(well) {}
WellGroupConstraints(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId,
const int,
const std::optional<std::string>&,
std::vector<double>&)>;
std::vector<Scalar>&)>;
bool checkGroupConstraints(WellState<double>& well_state,
const GroupState<double>& group_state,
bool checkGroupConstraints(WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
DeferredLogger& deferred_logger) const;
private:
std::pair<bool, double>
std::pair<bool, Scalar>
checkGroupConstraintsInj(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const double efficiencyFactor,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Scalar efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
DeferredLogger& deferred_logger) const;
std::pair<bool, double>
std::pair<bool, Scalar>
checkGroupConstraintsProd(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const double efficiencyFactor,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Scalar efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
};
}

123
opm/simulators/wells/WellGroupControls.cpp
View File

@@ -41,21 +41,21 @@
#include <cstddef>
#include <cassert>
namespace Opm
{
namespace Opm {
template<class Scalar>
template<class EvalWell>
void WellGroupControls::
void WellGroupControls<Scalar>::
getGroupInjectionControl(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
EvalWell& control_eq,
DeferredLogger& deferred_logger) const
{
@@ -128,10 +128,10 @@ getGroupInjectionControl(const Group& group,
// This is the group containing the control we will check against.
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(pu.num_phases, 1.0);
std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), std::nullopt, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
double sales_target = 0;
Scalar sales_target = 0;
if (schedule[well_.currentStep()].gconsale().has(group.name())) {
const auto& gconsale = schedule[well_.currentStep()].gconsale().get(group.name(), summaryState);
sales_target = gconsale.sales_target;
@@ -163,7 +163,7 @@ getGroupInjectionControl(const Group& group,
};
auto localReduction = [&](const std::string& group_name) {
const std::vector<double>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
const std::vector<Scalar>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
return tcalc.calcModeRateFromRates(groupTargetReductions);
};
@@ -171,13 +171,12 @@ getGroupInjectionControl(const Group& group,
if (!group.has_gpmaint_control(injectionPhase, currentGroupControl))
ctrl = group.injectionControls(injectionPhase, summaryState);
const double orig_target = tcalc.groupTarget(ctrl,
deferred_logger);
const auto chain = WellGroupHelpers<double>::groupChainTopBot(well_.name(), group.name(),
const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<Scalar>::groupChainTopBot(well_.name(), group.name(),
schedule, well_.currentStep());
// Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
const std::size_t num_ancestors = chain.size() - 1;
double target = orig_target;
Scalar target = orig_target;
for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
if ((ii == 0) || well_.guideRate()->has(chain[ii], injectionPhase)) {
// Apply local reductions only at the control level
@@ -188,22 +187,23 @@ getGroupInjectionControl(const Group& group,
target *= localFraction(chain[ii+1]);
}
// Avoid negative target rates coming from too large local reductions.
const double target_rate = std::max(0.0, target / efficiencyFactor);
const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor);
const auto current_rate = injection_rate;
control_eq = current_rate - target_rate;
}
std::optional<double>
WellGroupControls::
template<class Scalar>
std::optional<Scalar>
WellGroupControls<Scalar>::
getGroupInjectionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const
{
// Setting some defaults to silence warnings below.
@@ -262,10 +262,10 @@ getGroupInjectionTargetRate(const Group& group,
// This is the group containing the control we will check against.
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(pu.num_phases, 1.0);
std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), std::nullopt, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
double sales_target = 0;
Scalar sales_target = 0;
if (schedule[well_.currentStep()].gconsale().has(group.name())) {
const auto& gconsale = schedule[well_.currentStep()].gconsale().get(group.name(), summaryState);
sales_target = gconsale.sales_target;
@@ -296,7 +296,7 @@ getGroupInjectionTargetRate(const Group& group,
};
auto localReduction = [&](const std::string& group_name) {
const std::vector<double>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
const std::vector<Scalar>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
return tcalc.calcModeRateFromRates(groupTargetReductions);
};
@@ -304,15 +304,15 @@ getGroupInjectionTargetRate(const Group& group,
if (!group.has_gpmaint_control(injectionPhase, currentGroupControl))
ctrl = group.injectionControls(injectionPhase, summaryState);
const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<double>::groupChainTopBot(well_.name(),
const auto chain = WellGroupHelpers<Scalar>::groupChainTopBot(well_.name(),
group.name(),
schedule,
well_.currentStep());
// Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
const std::size_t num_ancestors = chain.size() - 1;
double target = orig_target;
Scalar target = orig_target;
for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
if ((ii == 0) || well_.guideRate()->has(chain[ii], injectionPhase)) {
// Apply local reductions only at the control level
@@ -322,21 +322,23 @@ getGroupInjectionTargetRate(const Group& group,
}
target *= localFraction(chain[ii+1]);
}
return std::max(0.0, target / efficiencyFactor);
return std::max(Scalar(0.0), target / efficiencyFactor);
}
template<class Scalar>
template<class EvalWell>
void WellGroupControls::getGroupProductionControl(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
const RateConvFunc& rateConverter,
double efficiencyFactor,
EvalWell& control_eq,
DeferredLogger& deferred_logger) const
void WellGroupControls<Scalar>::
getGroupProductionControl(const Group& group,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
const RateConvFunc& rateConverter,
Scalar efficiencyFactor,
EvalWell& control_eq,
DeferredLogger& deferred_logger) const
{
const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
if (currentGroupControl == Group::ProductionCMode::FLD ||
@@ -379,12 +381,12 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
// This is the group containing the control we will check against.
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
std::vector<Scalar> resv_coeff(well_.phaseUsage().num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
// gconsale may adjust the grat target.
// the adjusted rates is send to the targetCalculator
double gratTargetFromSales = 0.0;
Scalar gratTargetFromSales = 0.0;
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
@@ -411,7 +413,7 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
};
auto localReduction = [&](const std::string& group_name) {
const std::vector<double>& groupTargetReductions = group_state.production_reduction_rates(group_name);
const std::vector<Scalar>& groupTargetReductions = group_state.production_reduction_rates(group_name);
return tcalc.calcModeRateFromRates(groupTargetReductions);
};
@@ -419,12 +421,12 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
if (!group.has_gpmaint_control(currentGroupControl))
ctrl = group.productionControls(summaryState);
const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<double>::groupChainTopBot(well_.name(), group.name(),
const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<Scalar>::groupChainTopBot(well_.name(), group.name(),
schedule, well_.currentStep());
// Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
const std::size_t num_ancestors = chain.size() - 1;
double target = orig_target;
Scalar target = orig_target;
for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
if ((ii == 0) || well_.guideRate()->has(chain[ii])) {
// Apply local reductions only at the control level
@@ -435,19 +437,20 @@ void WellGroupControls::getGroupProductionControl(const Group& group,
target *= localFraction(chain[ii+1]);
}
// Avoid negative target rates coming from too large local reductions.
const double target_rate = std::max(0.0, target / efficiencyFactor);
const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor);
const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
control_eq = current_rate - target_rate;
}
double WellGroupControls::
template<class Scalar>
Scalar WellGroupControls<Scalar>::
getGroupProductionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const
{
const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
@@ -476,12 +479,12 @@ getGroupProductionTargetRate(const Group& group,
// This is the group containing the control we will check against.
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(well_.phaseUsage().num_phases, 1.0);
std::vector<Scalar> resv_coeff(well_.phaseUsage().num_phases, 1.0);
rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
// gconsale may adjust the grat target.
// the adjusted rates is send to the targetCalculator
double gratTargetFromSales = 0.0;
Scalar gratTargetFromSales = 0.0;
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
@@ -508,7 +511,7 @@ getGroupProductionTargetRate(const Group& group,
};
auto localReduction = [&](const std::string& group_name) {
const std::vector<double>& groupTargetReductions = group_state.production_reduction_rates(group_name);
const std::vector<Scalar>& groupTargetReductions = group_state.production_reduction_rates(group_name);
return tcalc.calcModeRateFromRates(groupTargetReductions);
};
@@ -516,12 +519,12 @@ getGroupProductionTargetRate(const Group& group,
if (!group.has_gpmaint_control(currentGroupControl))
ctrl = group.productionControls(summaryState);
const double orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<double>::groupChainTopBot(well_.name(), group.name(),
const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger);
const auto chain = WellGroupHelpers<Scalar>::groupChainTopBot(well_.name(), group.name(),
schedule, well_.currentStep());
// Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
const std::size_t num_ancestors = chain.size() - 1;
double target = orig_target;
Scalar target = orig_target;
for (std::size_t ii = 0; ii < num_ancestors; ++ii) {
if ((ii == 0) || well_.guideRate()->has(chain[ii])) {
// Apply local reductions only at the control level
@@ -532,22 +535,24 @@ getGroupProductionTargetRate(const Group& group,
target *= localFraction(chain[ii+1]);
}
// Avoid negative target rates coming from too large local reductions.
const double target_rate = std::max(0.0, target / efficiencyFactor);
const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor);
const auto& ws = well_state.well(well_.indexOfWell());
const auto& rates = ws.surface_rates;
const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
double scale = 1.0;
Scalar scale = 1.0;
if (target_rate == 0.0) {
return 0.0;
}
if (current_rate > 1e-14)
scale = target_rate/current_rate;
scale = target_rate / current_rate;
return scale;
}
template class WellGroupControls<double>;
#define INSTANCE(...) \
template void WellGroupControls:: \
template void WellGroupControls<double>:: \
getGroupInjectionControl<__VA_ARGS__>(const Group&, \
const WellState<double>&, \
const GroupState<double>&, \
@@ -560,7 +565,7 @@ getGroupInjectionControl<__VA_ARGS__>(const Group&, \
double efficiencyFactor, \
__VA_ARGS__& control_eq, \
DeferredLogger& deferred_logger) const; \
template void WellGroupControls:: \
template void WellGroupControls<double>:: \
getGroupProductionControl<__VA_ARGS__>(const Group&, \
const WellState<double>&, \
const GroupState<double>&, \

36
opm/simulators/wells/WellGroupControls.hpp
View File

@@ -43,62 +43,64 @@ template<class Scalar> class WellInterfaceGeneric;
template<class Scalar> class WellState;
//! \brief Class for computing well group controls.
template<class Scalar>
class WellGroupControls {
public:
//! \brief Constructor sets reference to well.
WellGroupControls(const WellInterfaceGeneric<double>& well) : well_(well) {}
WellGroupControls(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
using RateConvFunc = std::function<void(const RegionId, const int, const std::optional<std::string>&, std::vector<double>&)>;
using RateConvFunc = std::function<void(const RegionId, const int,
const std::optional<std::string>&, std::vector<Scalar>&)>;
template<class EvalWell>
void getGroupInjectionControl(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
EvalWell& control_eq,
DeferredLogger& deferred_logger) const;
std::optional<double>
std::optional<Scalar>
getGroupInjectionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const;
template<class EvalWell>
void getGroupProductionControl(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
EvalWell& control_eq,
DeferredLogger& deferred_logger) const;
double getGroupProductionTargetRate(const Group& group,
const WellState<double>& well_state,
const GroupState<double>& group_state,
Scalar getGroupProductionTargetRate(const Group& group,
const WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const RateConvFunc& rateConverter,
double efficiencyFactor,
Scalar efficiencyFactor,
DeferredLogger& deferred_logger) const;
private:
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
};
}

16
opm/simulators/wells/WellHelpers.cpp
View File

@@ -122,11 +122,12 @@ mmv (const X& x, Y& y) const
}
}
double computeHydrostaticCorrection(const double well_ref_depth, const double vfp_ref_depth,
const double rho, const double gravity)
template<class Scalar>
Scalar computeHydrostaticCorrection(const Scalar well_ref_depth, const Scalar vfp_ref_depth,
const Scalar rho, const Scalar gravity)
{
const double dh = vfp_ref_depth - well_ref_depth;
const double dp = rho * gravity * dh;
const Scalar dh = vfp_ref_depth - well_ref_depth;
const Scalar dp = rho * gravity * dh;
return dp;
}
@@ -200,7 +201,6 @@ bool rateControlWithZeroProdTarget(const WellProductionControls& controls,
}
}
bool rateControlWithZeroInjTarget(const WellInjectionControls& controls,
const WellInjectorCMode mode)
{
@@ -214,7 +214,6 @@ bool rateControlWithZeroInjTarget(const WellInjectionControls& controls,
}
}
template class ParallelStandardWellB<double>;
template<int Dim> using Vec = Dune::BlockVector<Dune::FieldVector<double,Dim>>;
@@ -239,5 +238,10 @@ template void sumDistributedWellEntries<double,Comm>(Dune::DynamicMatrix<double>
using DMatrix = Dune::DynamicMatrix<double>;
template DMatrix transposeDenseDynMatrix<DMatrix>(const DMatrix&);
template double computeHydrostaticCorrection<double>(const double,
const double,
const double,
const double);
} // namespace wellhelpers
} // namespace Opm

11
opm/simulators/wells/WellHelpers.hpp
View File

@@ -26,8 +26,6 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/common/dynmatrix.hh>
#include <array>
namespace Opm {
class ParallelWellInfo;
@@ -70,10 +68,11 @@ private:
const ParallelWellInfo& parallel_well_info_;
};
double computeHydrostaticCorrection(const double well_ref_depth,
const double vfp_ref_depth,
const double rho, const double gravity);
template<class Scalar>
Scalar computeHydrostaticCorrection(const Scalar well_ref_depth,
const Scalar vfp_ref_depth,
const Scalar rho,
const Scalar gravity);
/// \brief Sums entries of the diagonal Matrix for distributed wells
template<typename Scalar, typename Comm>

163
opm/simulators/wells/WellTest.cpp
View File

@@ -32,55 +32,58 @@
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
namespace Opm
{
namespace Opm {
template<class Scalar>
template<class RatioFunc>
bool WellTest::checkMaxRatioLimitWell(const SingleWellState<double>& ws,
const double max_ratio_limit,
const RatioFunc& ratioFunc) const
bool WellTest<Scalar>::
checkMaxRatioLimitWell(const SingleWellState<Scalar>& ws,
const Scalar max_ratio_limit,
const RatioFunc& ratioFunc) const
{
const int np = well_.numPhases();
std::vector<double> well_rates(np, 0.0);
std::vector<Scalar> well_rates(np, 0.0);
for (int p = 0; p < np; ++p) {
well_rates[p] = ws.surface_rates[p];
}
const double well_ratio = ratioFunc(well_rates, well_.phaseUsage());
const Scalar well_ratio = ratioFunc(well_rates, well_.phaseUsage());
return (well_ratio > max_ratio_limit);
}
template<class Scalar>
template<class RatioFunc>
void WellTest::checkMaxRatioLimitCompletions(const SingleWellState<double>& ws,
const double max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const
void WellTest<Scalar>::
checkMaxRatioLimitCompletions(const SingleWellState<Scalar>& ws,
const Scalar max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const
{
int worst_offending_completion = RatioLimitCheckReport::INVALIDCOMPLETION;
// the maximum water cut value of the completions
// it is used to identify the most offending completion
double max_ratio_completion = 0;
Scalar max_ratio_completion = 0;
const int np = well_.numPhases();
const auto& perf_data = ws.perf_data;
const auto& perf_phase_rates = perf_data.phase_rates;
// look for the worst_offending_completion
for (const auto& completion : well_.getCompletions()) {
std::vector<double> completion_rates(np, 0.0);
std::vector<Scalar> completion_rates(np, 0.0);
// looping through the connections associated with the completion
const std::vector<int>& conns = completion.second;
for (const int c : conns) {
for (int p = 0; p < np; ++p) {
const double connection_rate = perf_phase_rates[c * np + p];
const Scalar connection_rate = perf_phase_rates[c * np + p];
completion_rates[p] += connection_rate;
}
} // end of for (const int c : conns)
well_.parallelWellInfo().communication().sum(completion_rates.data(), completion_rates.size());
const double ratio_completion = ratioFunc(completion_rates, well_.phaseUsage());
const Scalar ratio_completion = ratioFunc(completion_rates, well_.phaseUsage());
if (ratio_completion > max_ratio_completion) {
worst_offending_completion = completion.first;
@@ -88,7 +91,7 @@ void WellTest::checkMaxRatioLimitCompletions(const SingleWellState<double>& ws,
}
} // end of for (const auto& completion : completions_)
const double violation_extent = max_ratio_completion / max_ratio_limit;
const Scalar violation_extent = max_ratio_completion / max_ratio_limit;
if (violation_extent > report.violation_extent) {
report.worst_offending_completion = worst_offending_completion;
@@ -96,27 +99,30 @@ void WellTest::checkMaxRatioLimitCompletions(const SingleWellState<double>& ws,
}
}
void WellTest::checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
RatioLimitCheckReport& report) const
template<class Scalar>
void WellTest<Scalar>::
checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const
{
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// function to calculate gor based on rates
auto gor = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double oil_rate = -rates[pu.phase_pos[Oil]];
const double gas_rate = -rates[pu.phase_pos[Gas]];
auto gor = [](const std::vector<Scalar>& rates,
const PhaseUsage& pu)
{
const Scalar oil_rate = -rates[pu.phase_pos[Oil]];
const Scalar gas_rate = -rates[pu.phase_pos[Gas]];
if (gas_rate <= 0.)
return 0.;
return Scalar{0};
else if (oil_rate <= 0.)
return 1.e100; // big value to mark it as violated
return Scalar{1e30}; // big value to mark it as violated
else
return (gas_rate / oil_rate);
};
const double max_gor_limit = econ_production_limits.maxGasOilRatio();
const Scalar max_gor_limit = econ_production_limits.maxGasOilRatio();
assert(max_gor_limit > 0.);
const bool gor_limit_violated = this->checkMaxRatioLimitWell(ws, max_gor_limit, gor);
@@ -127,28 +133,30 @@ void WellTest::checkMaxGORLimit(const WellEconProductionLimits& econ_production_
}
}
void WellTest::checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
RatioLimitCheckReport& report) const
template<class Scalar>
void WellTest<Scalar>::
checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Water = BlackoilPhases::Aqua;
// function to calculate wgr based on rates
auto wgr = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double water_rate = -rates[pu.phase_pos[Water]];
const double gas_rate = -rates[pu.phase_pos[Gas]];
auto wgr = [](const std::vector<Scalar>& rates,
const PhaseUsage& pu)
{
const Scalar water_rate = -rates[pu.phase_pos[Water]];
const Scalar gas_rate = -rates[pu.phase_pos[Gas]];
if (water_rate <= 0.)
return 0.;
return Scalar{0};
else if (gas_rate <= 0.)
return 1.e100; // big value to mark it as violated
return Scalar{1e30}; // big value to mark it as violated
else
return (water_rate / gas_rate);
};
const double max_wgr_limit = econ_production_limits.maxWaterGasRatio();
const Scalar max_wgr_limit = econ_production_limits.maxWaterGasRatio();
assert(max_wgr_limit > 0.);
const bool wgr_limit_violated = this->checkMaxRatioLimitWell(ws, max_wgr_limit, wgr);
@@ -159,31 +167,34 @@ void WellTest::checkMaxWGRLimit(const WellEconProductionLimits& econ_production_
}
}
void WellTest::checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
RatioLimitCheckReport& report) const
template<class Scalar>
void WellTest<Scalar>::
checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const
{
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
// function to calculate water cut based on rates
auto waterCut = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double oil_rate = -rates[pu.phase_pos[Oil]];
const double water_rate = -rates[pu.phase_pos[Water]];
const double liquid_rate = oil_rate + water_rate;
auto waterCut = [](const std::vector<Scalar>& rates,
const PhaseUsage& pu)
{
const Scalar oil_rate = -rates[pu.phase_pos[Oil]];
const Scalar water_rate = -rates[pu.phase_pos[Water]];
const Scalar liquid_rate = oil_rate + water_rate;
if (liquid_rate <= 0.)
return 0.;
return Scalar{0};
else if (water_rate < 0)
return 0.;
return Scalar{0};
else if (oil_rate < 0)
return 1.;
return Scalar{1};
else
return (water_rate / liquid_rate);
};
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
const Scalar max_water_cut_limit = econ_production_limits.maxWaterCut();
assert(max_water_cut_limit > 0.);
const bool watercut_limit_violated =
@@ -196,9 +207,11 @@ void WellTest::checkMaxWaterCutLimit(const WellEconProductionLimits& econ_produc
}
}
bool WellTest::checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const std::vector<double>& rates_or_potentials,
DeferredLogger& deferred_logger) const
template<class Scalar>
bool WellTest<Scalar>::
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const std::vector<Scalar>& rates_or_potentials,
DeferredLogger& deferred_logger) const
{
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
@@ -207,26 +220,26 @@ bool WellTest::checkRateEconLimits(const WellEconProductionLimits& econ_producti
const PhaseUsage& pu = well_.phaseUsage();
if (econ_production_limits.onMinOilRate()) {
const double oil_rate = rates_or_potentials[pu.phase_pos[Oil]];
const double min_oil_rate = econ_production_limits.minOilRate();
const Scalar oil_rate = rates_or_potentials[pu.phase_pos[Oil]];
const Scalar min_oil_rate = econ_production_limits.minOilRate();
if (std::abs(oil_rate) < min_oil_rate) {
return true;
}
}
if (econ_production_limits.onMinGasRate() ) {
const double gas_rate = rates_or_potentials[pu.phase_pos[Gas]];
const double min_gas_rate = econ_production_limits.minGasRate();
const Scalar gas_rate = rates_or_potentials[pu.phase_pos[Gas]];
const Scalar min_gas_rate = econ_production_limits.minGasRate();
if (std::abs(gas_rate) < min_gas_rate) {
return true;
}
}
if (econ_production_limits.onMinLiquidRate() ) {
const double oil_rate = rates_or_potentials[pu.phase_pos[Oil]];
const double water_rate = rates_or_potentials[pu.phase_pos[Water]];
const double liquid_rate = oil_rate + water_rate;
const double min_liquid_rate = econ_production_limits.minLiquidRate();
const Scalar oil_rate = rates_or_potentials[pu.phase_pos[Oil]];
const Scalar water_rate = rates_or_potentials[pu.phase_pos[Water]];
const Scalar liquid_rate = oil_rate + water_rate;
const Scalar min_liquid_rate = econ_production_limits.minLiquidRate();
if (std::abs(liquid_rate) < min_liquid_rate) {
return true;
}
@@ -239,9 +252,11 @@ bool WellTest::checkRateEconLimits(const WellEconProductionLimits& econ_producti
return false;
}
WellTest::RatioLimitCheckReport WellTest::
template<class Scalar>
typename WellTest<Scalar>::RatioLimitCheckReport
WellTest<Scalar>::
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
const SingleWellState<Scalar>& ws,
DeferredLogger& deferred_logger) const
{
// TODO: not sure how to define the worst-offending completion when more than one
@@ -289,11 +304,13 @@ checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
return report;
}
void WellTest::updateWellTestStateEconomic(const SingleWellState<double>& ws,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const
template<class Scalar>
void WellTest<Scalar>::
updateWellTestStateEconomic(const SingleWellState<Scalar>& ws,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const
{
if (well_.wellIsStopped())
return;
@@ -440,10 +457,12 @@ void WellTest::updateWellTestStateEconomic(const SingleWellState<double>& ws,
}
}
void WellTest::updateWellTestStatePhysical(const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const
template<class Scalar>
void WellTest<Scalar>::
updateWellTestStatePhysical(const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const
{
if (well_test_state.well_is_closed(well_.name())) {
// Already closed, do nothing.
@@ -458,4 +477,6 @@ void WellTest::updateWellTestStatePhysical(const double simulation_time,
}
}
template class WellTest<double>;
} // namespace Opm

27
opm/simulators/wells/WellTest.hpp
View File

@@ -38,12 +38,13 @@ template<class Scalar> class WellInterfaceGeneric;
class WellTestState;
//! \brief Class for performing well tests.
template<class Scalar>
class WellTest {
public:
//! \brief Constructor sets reference to well.
WellTest(const WellInterfaceGeneric<double>& well) : well_(well) {}
WellTest(const WellInterfaceGeneric<Scalar>& well) : well_(well) {}
void updateWellTestStateEconomic(const SingleWellState<double>& ws,
void updateWellTestStateEconomic(const SingleWellState<Scalar>& ws,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
@@ -59,43 +60,43 @@ private:
static constexpr int INVALIDCOMPLETION = std::numeric_limits<int>::max();
bool ratio_limit_violated = false;
int worst_offending_completion = INVALIDCOMPLETION;
double violation_extent = 0.0;
Scalar violation_extent = 0.0;
};
void checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const;
void checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const;
void checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
const SingleWellState<Scalar>& ws,
RatioLimitCheckReport& report) const;
template<class RatioFunc>
bool checkMaxRatioLimitWell(const SingleWellState<double>& ws,
const double max_ratio_limit,
bool checkMaxRatioLimitWell(const SingleWellState<Scalar>& ws,
const Scalar max_ratio_limit,
const RatioFunc& ratioFunc) const;
template<class RatioFunc>
void checkMaxRatioLimitCompletions(const SingleWellState<double>& ws,
const double max_ratio_limit,
void checkMaxRatioLimitCompletions(const SingleWellState<Scalar>& ws,
const Scalar max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const std::vector<double>& rates_or_potentials,
const std::vector<Scalar>& rates_or_potentials,
DeferredLogger& deferred_logger) const;
RatioLimitCheckReport
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const SingleWellState<double>& ws,
const SingleWellState<Scalar>& ws,
DeferredLogger& deferred_logger) const;
const WellInterfaceGeneric<double>& well_; //!< Reference to well interface
const WellInterfaceGeneric<Scalar>& well_; //!< Reference to well interface
};
}