changed assessing safe THP range

This commit is contained in:
Paul 2024-03-24 20:47:46 +01:00
parent 393c70a83e
commit 6e76602e8f
6 changed files with 106 additions and 95 deletions

View File

@ -466,7 +466,7 @@ template<class Scalar> class WellContributions;
const double dt,
DeferredLogger& local_deferredLogger);
void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger);
void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger);
/// Update rank's notion of intersecting wells and their
/// associate solution variables.

View File

@ -49,12 +49,9 @@
#include <opm/simulators/utils/MPIPacker.hpp>
#include <opm/simulators/utils/phaseUsageFromDeck.hpp>
<<<<<<< HEAD
#if COMPILE_BDA_BRIDGE
#include <opm/simulators/linalg/bda/WellContributions.hpp>
#endif
=======
>>>>>>> f92db77d3 (moved common thp calculation to updateWellControls)
#if HAVE_MPI
#include <opm/simulators/utils/MPISerializer.hpp>
@ -1258,8 +1255,10 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger)
{
const int reportStepIdx = this->ebosSimulator_.episodeIndex();
const auto& network = schedule()[reportStepIdx].network();
const int reportStepIdx = this->simulator_.episodeIndex();
const auto& network = this->schedule()[reportStepIdx].network();
const auto& balance = this->schedule()[reportStepIdx].network_balance();
const Scalar thp_tolerance = balance.thp_tolerance();
if (!network.active()) {
return;
@ -1271,75 +1270,87 @@ namespace Opm {
for (const std::string& nodeName : network.node_names()) {
const bool has_choke = network.node(nodeName).as_choke();
if (has_choke) {
const auto& summary_state = this->ebosSimulator_.vanguard().summaryState();
const Group& group = schedule().getGroup(nodeName, reportStepIdx);
const auto& summary_state = this->simulator_.vanguard().summaryState();
const Group& group = this->schedule().getGroup(nodeName, reportStepIdx);
const auto ctrl = group.productionControls(summary_state);
const auto cmode = ctrl.cmode;
const auto pu = this->phase_usage_;
//TODO: Auto choke combined with RESV control is not supported
std::vector<double> resv_coeff(pu.num_phases, 1.0);
double gratTargetFromSales = 0.0;
std::vector<Scalar> resv_coeff(pu.num_phases, 1.0);
Scalar gratTargetFromSales = 0.0;
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
gratTargetFromSales, nodeName, group_state,
group.has_gpmaint_control(cmode));
const double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
auto mismatch = [&] (auto well_group_thp) {
double well_group_rate(0.0);
double rate(0.0);
Scalar well_group_rate(0.0);
Scalar rate(0.0);
for (auto& well : this->well_container_) {
std::string well_name = well->name();
auto& ws = well_state.well(well_name);
if (group.hasWell(well_name)) {
well->setDynamicThpLimit(well_group_thp);
const Well& well_ecl = wells_ecl_[well->indexOfWell()];
const Well& well_ecl = this->wells_ecl_[well->indexOfWell()];
const auto inj_controls = Well::InjectionControls(0);
const auto prod_controls = well_ecl.productionControls(summary_state);
well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false);
well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false);
rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
well_group_rate += rate;
}
}
return well_group_rate - orig_target;
return (well_group_rate - orig_target)/orig_target;
};
double thp(0.0);
double min_thp(1.0E8);
double max_thp(0.0);
for (auto& well : this->well_container_) {
std::string well_name = well->name();
if (group.hasWell(well_name)) {
auto& ws = well_state.well(well_name);
thp = ws.thp;
min_thp = std::min(min_thp, thp);
max_thp = std::max(max_thp, thp);
Scalar min_thp, max_thp;
std::array<Scalar, 2> range_initial;
//Find an initial bracket
if (!this->well_group_thp_calc_.has_value()){
// Retrieve the terminal pressure of the associated root of the manifold group
std::string node_name = nodeName;
while (!network.node(node_name).terminal_pressure().has_value()) {
auto branch = network.uptree_branch(node_name).value();
node_name = branch.uptree_node();
}
min_thp = network.node(node_name).terminal_pressure().value();
WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp);
// Narrow down the bracket
Scalar low1, high1;
std::array<Scalar, 2> range = {0.9*min_thp, 1.1*max_thp};
std::optional<Scalar> appr_sol;
WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger);
min_thp = low1;
max_thp = high1;
range_initial = {min_thp, max_thp};
}
const auto upbranch = network.uptree_branch(nodeName);
const auto it = node_pressures_.find((*upbranch).uptree_node());
const double nodal_pressure = it->second;
double well_group_thp = nodal_pressure;
const auto it = this->node_pressures_.find((*upbranch).uptree_node());
const Scalar nodal_pressure = it->second;
Scalar well_group_thp = nodal_pressure;
if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) {
std::array<double, 2> range;
// The bracket is based on the initial bracket or on a range based on a previous calculated common group thp
std::array<Scalar, 2> range;
this->well_group_thp_calc_.has_value() ?
range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : range = {0.9 * min_thp, 1.1 * max_thp};
range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} :
range = range_initial;
double low, high;
std::optional<double> approximate_solution;
const double tolerance1 = 1.0E-3;
Scalar low, high;
std::optional<Scalar> approximate_solution;
const Scalar tolerance1 = thp_tolerance;
local_deferredLogger.debug("Using brute force search to bracket the common THP");
const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
const bool finding_bracket = WellBhpThpCalculator<Scalar>::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger);
if (approximate_solution.has_value()) {
this->well_group_thp_calc_ = *approximate_solution;
local_deferredLogger.debug("Approximate common THP value found: " + std::to_string(this->well_group_thp_calc_.value()));
} else if (finding_bracket) {
const double tolerance2 = 1.0E-3;
const Scalar tolerance2 = thp_tolerance;
const int max_iteration_solve = 100;
int iteration = 0;
this->well_group_thp_calc_= RegulaFalsiBisection<ThrowOnError>::
@ -1348,11 +1359,12 @@ namespace Opm {
"iteration = " + std::to_string(iteration));
local_deferredLogger.debug("Common THP value = " + std::to_string(this->well_group_thp_calc_.value()));
} else {
local_deferredLogger.warning("Common THP solve failed due to bracketing failure");
this->well_group_thp_calc_ = {};
local_deferredLogger.debug("Common THP solve failed due to bracketing failure");
}
}
well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure);
this->well_group_thp_calc_.has_value() ?
well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure) : well_group_thp = nodal_pressure;
for (auto& well : this->well_container_) {
std::string well_name = well->name();
@ -1365,7 +1377,7 @@ namespace Opm {
}
}
// Use the common THP in computeNetworkPressures
// Use the common group THP in computeNetworkPressures
group_state.update_well_group_thp(nodeName, well_group_thp);
}
}
@ -2031,7 +2043,7 @@ namespace Opm {
// network related
bool more_network_update = false;
if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
const double dt = this->ebosSimulator_.timeStepSize();
const double dt = this->simulator_.timeStepSize();
// Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
computeWellGroupThp(dt, deferred_logger);
const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx);

View File

@ -101,8 +101,9 @@ GroupState<Scalar>::production_rates(const std::string& gname) const
}
//-------------------------------------------------------------------------
template<class Scalar>
bool GroupState<Scalar>::
void GroupState<Scalar>::
GroupState::update_well_group_thp(const std::string& gname, const double& thp)
{
this->group_thp[gname] = thp;
@ -119,6 +120,8 @@ GroupState::well_group_thp(const std::string& gname) const
return group_iter->second;
}
//-------------------------------------------------------------------------
template<class Scalar>
bool GroupState<Scalar>::
has_production_reduction_rates(const std::string& gname) const

View File

@ -875,15 +875,9 @@ bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
low = range[0];
high = range[1];
const int sample_number = 200;
<<<<<<< HEAD
const Scalar interval = (high - low) / sample_number;
Scalar eq_low = eq(low);
Scalar eq_high = 0.0;
=======
const double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high = 0.0;
>>>>>>> 3ce8d7eec (moved common thp calculation to updateWellControls)
for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i;
eq_high = eq(high);
@ -1011,22 +1005,23 @@ getFloIPR(const WellState<Scalar>& well_state,
detail::getFlo(table, aqua_b, liquid_b, vapour_b));
}
template<class Scalar>
bool
WellBhpThpCalculator::
bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
const double& limit,
WellBhpThpCalculator<Scalar>::
bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
std::optional<Scalar>& approximate_solution,
const Scalar& limit,
DeferredLogger& deferred_logger)
{
bool bracket_found = false;
low = range[0];
high = range[1];
const int sample_number = 300;
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);
@ -1050,6 +1045,30 @@ bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
return bracket_found;
}
template<class Scalar>
bool
WellBhpThpCalculator<Scalar>::
bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
Scalar& min_thp, Scalar& max_thp)
{
bool bracket_found = false;
const int sample_number = 1000;
const Scalar interval = 1E5;
Scalar eq_low = eq(min_thp);
Scalar eq_high = 0.0;
for (int i = 0; i < sample_number + 1; ++i) {
max_thp = min_thp + interval * i;
eq_high = eq(max_thp);
if (eq_high * eq_low <= 0.) {
bracket_found = true;
min_thp = max_thp - interval;
break;
}
eq_low = eq_high;
}
return bracket_found;
}
template class WellBhpThpCalculator<double>;
#define INSTANCE(...) \

View File

@ -120,19 +120,23 @@ public:
const SummaryState& summary_state) 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);
//! \brief Find limits using brute-force solver.
static bool bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
std::optional<double>& approximate_solution,
const double& limit,
static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
std::optional<Scalar>& approximate_solution,
const Scalar& limit,
DeferredLogger& deferred_logger);
//! \brief Find limits using brute-force solver.
static bool bruteForceBracketCommonTHP(const std::function<Scalar(const Scalar)>& eq,
Scalar& min_thp, Scalar& max_thp);
private:
//! \brief Compute BHP from THP limit for an injector - implementation.
template<class ErrorPolicy>
@ -169,12 +173,6 @@ private:
std::optional<Scalar>& approximate_solution,
DeferredLogger& deferred_logger) const;
//! \brief Find limits using brute-force solver.
static bool bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
const std::array<Scalar, 2>& range,
Scalar& low, Scalar& high,
DeferredLogger& deferred_logger);
Scalar findThpFromBhpIteratively(const std::function<Scalar(const Scalar, const Scalar)>& thp_func,
const Scalar bhp,
const Scalar thp_limit,

View File

@ -374,26 +374,15 @@ public:
void updateConnectionTransmissibilityFactor(const Simulator& simulator,
SingleWellState<Scalar>& ws) const;
SingleWellState<double>& ws) const;
virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
const double dt,
const WellInjectionControls& inj_controls,
const WellProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false) = 0;
bool solveWellWithTHPConstraint(const Simulator& simulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
protected:
// simulation parameters
const ModelParameters& param_;
@ -448,16 +437,6 @@ protected:
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger) = 0;
virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
const double dt,
const WellInjectionControls& inj_controls,
const WellProductionControls& prod_controls,
WellState<Scalar>& well_state,
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger,
const bool fixed_control = false,
const bool fixed_status = false) = 0;
virtual void updateIPRImplicit(const Simulator& simulator,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) = 0;