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

clean up and improvements according reviewer comments

Browse Source
This commit is contained in:
Paul 2024-02-19 22:09:33 +01:00
parent 1ddf675cfd
commit 393c70a83e
6 changed files with 35 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -413,7 +413,7 @@ template<class Scalar> class WellContributions;
Scalar gravity_{};
std::vector<Scalar> depth_{};
bool alternative_well_rate_init_{};
double well_group_thp_prev_{0.0};
std::optional<double> well_group_thp_calc_;
std::unique_ptr<RateConverterType> rateConverter_{};
std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};

1
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@ -1363,6 +1363,7 @@ updateNetworkPressures(const int reportStepIdx)
}
for (auto& well : well_container_generic_) {
// Producers only, since we so far only support the
// "extended" network model (properties defined by
// BRANPROP and NODEPROP) which only applies to producers.

107
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

@ -1181,7 +1181,6 @@ namespace Opm {
if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) {
local_deferredLogger.info(" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations ");
}
const bool relax_network_balance = network_update_iteration >= iteration_to_relax;
std::tie(do_network_update, well_group_control_changed) =
updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger);
@ -1274,18 +1273,19 @@ namespace Opm {
if (has_choke) {
const auto& summary_state = this->ebosSimulator_.vanguard().summaryState();
const Group& group = schedule().getGroup(nodeName, reportStepIdx);
std::optional<Group::ProductionControls> ctrl = group.productionControls(summary_state);
const auto cmode = ctrl->cmode;
const auto ctrl = group.productionControls(summary_state);
const auto cmode = ctrl.cmode;
const auto pu = this->phase_usage_;
//PJPE: conversion factor res<-> surface rates TODO: to be calculated
//TODO: Auto choke combined with RESV control is not supported
std::vector<double> resv_coeff(pu.num_phases, 1.0);
// rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
double gratTargetFromSales = 0.0; //PJPE: check this
double 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,
gratTargetFromSales, nodeName, group_state,
group.has_gpmaint_control(cmode));
double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
//PJPE: TODO: include efficiency factor
const double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger);
auto mismatch = [&] (auto well_group_thp) {
double well_group_rate(0.0);
@ -1294,32 +1294,18 @@ namespace Opm {
std::string well_name = well->name();
auto& ws = well_state.well(well_name);
if (group.hasWell(well_name)) {
rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
std::cout << "(1) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl;
well->setDynamicThpLimit(well_group_thp);
const Well& well_ecl = 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);
// May be used instead if the wells of the sub-sea manifold has only THP constraints
// well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger);
rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
if (ws.production_cmode == Opm::WellProducerCMode::THP)
ws.thp = well_group_thp;
std::cout << "(2) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl;
well_group_rate += rate;
}
}
double diff = well_group_rate - orig_target; // (well_group_rate - orig_target)/orig_target;
std::cout << "well_group_thp: " << well_group_thp << " mismatch: " << diff << "\n" << std::endl;
return diff;
return well_group_rate - orig_target;
};
double well_group_rate(0.0);
double rate(0.0);
double thp(0.0);
double min_thp(1.0E8);
double max_thp(0.0);
@ -1327,8 +1313,6 @@ namespace Opm {
std::string well_name = well->name();
if (group.hasWell(well_name)) {
auto& ws = well_state.well(well_name);
rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
well_group_rate += rate;
thp = ws.thp;
min_thp = std::min(min_thp, thp);
max_thp = std::max(max_thp, thp);
@ -1340,83 +1324,46 @@ namespace Opm {
const double nodal_pressure = it->second;
double well_group_thp = nodal_pressure;
if ((reportStepIdx <= 1) || (well_group_rate > 0.9*orig_target)) {
std::array<double, 2> range; // = {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here?
if (well_group_thp_prev_ == 0) {
range = {0.9*min_thp, 1.1*max_thp};
}
else {
range = {0.9* this->well_group_thp_prev_, 1.1*this->well_group_thp_prev_};
}
if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) {
std::array<double, 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};
double low, high;
// For testing
// low = 3.42767e+06;// 3.32833e+06;
// // high = 3.42767e+06;// 3.32833e+06;
// // std::cout << "low: " << low << " high: " << high << std::endl;
// double low_value_test = mismatch(low);
// low_value_test = mismatch(low);
// low_value_test = mismatch(low);
// std::cout << "low_value: " << low_value_test << " high_value: " << high_value_test << std::endl;
std::optional<double> approximate_solution;
const double limit = 1.0E-3;// 1.0E-4;
local_deferredLogger.debug(" Using brute force search to bracket the common THP ");
std::cout << " Using brute force search to bracket the common THP " << std::endl;
const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, limit, local_deferredLogger);
double low_value = mismatch(low);
double high_value = mismatch(high);
const double tolerance1 = 1.0E-3;
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);
if (approximate_solution.has_value()) {
well_group_thp = *approximate_solution;
double check = mismatch(well_group_thp);
local_deferredLogger.debug("Approximate value found: mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp));
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) {
// PJPE: TODO what settings to take here?
std::cout << "Found bracket: [" << low << ", " << high << "]; mismatch: [" << low_value << ", "<< high_value << "]" << std::endl;
const double tolerance = 1.0E-3;
const double tolerance2 = 1.0E-3;
const int max_iteration_solve = 100;
int iteration = 0;
well_group_thp = RegulaFalsiBisection<ThrowOnError>::
solve(mismatch, low, high, max_iteration_solve, tolerance, iteration);
local_deferredLogger.debug( " bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " +
"mismatch = [" + std::to_string(low_value) + ", " + std::to_string(high_value) + "], " +
this->well_group_thp_calc_= RegulaFalsiBisection<ThrowOnError>::
solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration);
local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " +
"iteration = " + std::to_string(iteration));
double check = mismatch(well_group_thp);
local_deferredLogger.debug(" mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp));
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");
well_group_thp = nodal_pressure;
}
}
std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl;
} else {
well_group_thp = nodal_pressure;
}
well_group_thp = std::max(well_group_thp, nodal_pressure);
well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure);
for (auto& well : this->well_container_) {
std::string well_name = well->name();
if (group.hasWell(well_name)) {
auto& ws = well_state.well(well_name);
std::cout << "Cmode of Well: " << well_name << " is " << ws.production_cmode <<std::endl;
if (ws.production_cmode == Opm::WellProducerCMode::THP){
if (ws.production_cmode == Opm::WellProducerCMode::THP) {
well->setDynamicThpLimit(well_group_thp);
ws.thp = well_group_thp;
}
// const Well& well_ecl = wells_ecl_[well->indexOfWell()];
// const auto inj_controls = Well::InjectionControls(0);
// const auto prod_controls = well_ecl.productionControls(summary_state);
// // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger);
// well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state_copy, group_state, local_deferredLogger);
if (ws.production_cmode == Opm::WellProducerCMode::THP)
ws.thp = well_group_thp;
}
}
well_group_thp_prev_ = well_group_thp;
// Use the common THP in computeNetworkPressures
group_state.update_well_group_thp(nodeName, well_group_thp);
@ -2555,11 +2502,11 @@ namespace Opm {
const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates";
deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
}
}
// If we're using local well solves that include control switches, they also update
// operability, so reset before main iterations begin
well->resetWellOperability();
}
}
updatePrimaryVariables(deferred_logger);
// Actually do the pre-step network rebalance, using the updated well states and initial solutions

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

@ -1021,10 +1021,9 @@ bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
DeferredLogger& deferred_logger)
{
bool bracket_found = false;
// bool approximate_solution_found = false;
low = range[0];
high = range[1];
const int sample_number = 300;//10000; //300; //10000;
const int sample_number = 300;
const double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high = 0.0;
@ -1032,7 +1031,6 @@ bruteForceBracketCommonTHP(const std::function<double(const double)>& eq,
high = range[0] + interval * i;
eq_high = eq(high);
if ( (std::fabs(eq_high) < limit)) {
// approximate_solution_found = true;
approximate_solution = high;
break;
}