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

moved common thp calculation to updateWellControls

Browse Source
This commit is contained in:
Paul
2024-01-22 13:47:06 +01:00
parent 6ddf5dd01b
commit 1ddf675cfd
8 changed files with 185 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -413,6 +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::unique_ptr<RateConverterType> rateConverter_{};
std::map<std::string, std::unique_ptr<AverageRegionalPressureType>> regionalAveragePressureCalculator_{};

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

@@ -49,9 +49,12 @@
#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>
@@ -1166,8 +1169,6 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger)
{
// PJPE: calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES)
computeWellGroupThp(dt,local_deferredLogger);
// not necessarily that we always need to update once of the network solutions
bool do_network_update = true;
bool well_group_control_changed = false;
@@ -1250,7 +1251,7 @@ namespace Opm {
return {more_network_update, well_group_control_changed};
}
//PJPE: This function is to be used for well groups in an extended network that act as a subsea manifold
// This function is to be used for well groups in an extended network that act as a subsea manifold
// The wells of such group should have a common THP and total phase rate(s) obeying (if possible)
// the well group constraint set by GCONPROD
template <typename TypeTag>
@@ -1276,8 +1277,9 @@ namespace Opm {
std::optional<Group::ProductionControls> 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
//PJPE: conversion factor res<-> surface rates TODO: to be calculated
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
WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff,
gratTargetFromSales, nodeName, group_state,
@@ -1286,92 +1288,116 @@ namespace Opm {
//PJPE: TODO: include efficiency factor
auto mismatch = [&] (auto well_group_thp) {
// auto well_state_copy = this->wellState();
double well_group_rate(0.0);
double 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)) {
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);
// well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger);
auto& ws = well_state.well(well_name);
double rate = -tcalc.calcModeRateFromRates(ws.surface_rates);
std::cout << "well: " << well_name << "rate: " << rate << std::endl;
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;
}
}
// this->wellState() = well_state;
return well_group_rate - orig_target;
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;
};
double well_group_rate(0.0);
double rate(0.0);
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);
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);
}
}
const auto upbranch = network.uptree_branch(nodeName);
const auto it = node_pressures_.find((*upbranch).uptree_node());
const double nodal_pressure = it->second;
const std::array<double, 2> range {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here?
double low, high;
low = 3527000.0;
high = 3576666.7;
std::cout << "low: " << low << " high: " << high << std::endl;
double low_v = mismatch(low);
std::cout << "1: low_value: " << low_v << std::endl;
low_v = mismatch(low);
high_v = mismatch(high);
std::cout << "2: low_value: " << low_v << " high_value: " << high_v << std::endl;
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_};
}
local_deferredLogger.debug(" Trying the brute force search to bracket the common THP ");
const bool finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger);
std::cout << "low: " << low << " high: " << high << std::endl;
double low_value = mismatch(low);
double high_value = mismatch(high);
std::cout << "low_value: " << mismatch(low) << " high_value: " << mismatch(high) << std::endl;
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;
double well_group_thp = 0.0;
if ((low_value * high_value >= 0.0) && (low >= nodal_pressure) ) {
// In case no bracket is found
const double limit = 0.0003; // PJPE what limit to be used?
double abs_low = std::fabs(low_value);
double abs_high = std::fabs(high_value);
if (std::min(abs_low, abs_high) < limit) {
well_group_thp = abs_low < abs_high ? low : high;
local_deferredLogger.warning("Common THP not solved precisely");
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);
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));
} 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 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) + "], " +
"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));
} else {
local_deferredLogger.warning("Common THP solve failed due to bracketing failure");
well_group_thp = nodal_pressure;
}
} else if (high <= nodal_pressure) {
local_deferredLogger.warning("Common THP set to nodal pressure");
well_group_thp = nodal_pressure;
std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl;
} else {
// PJPE: TODO what settings to take here?
const double tolerance = 1E-4;
const int max_iteration_solve = 20;
int iteration = 0;
well_group_thp = RegulaFalsiBisection<ThrowOnError>::
solve(mismatch, low, high, max_iteration_solve, tolerance, iteration);
std::cout << "Bracket = [" << low << ", " << high <<"] " << ", iteration = " << iteration << std::endl;
double check = mismatch(well_group_thp);
std::cout << "mismatch = " << check << ", well_group_thp = " << well_group_thp << std::endl;
well_group_thp = nodal_pressure;
}
// well_group_thp = std::max(well_group_thp, nodal_pressure);
std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl;
well_group_thp = std::max(well_group_thp, nodal_pressure);
// if (well_group_thp < 2.5E6){
// double min = 24.0E5;
// double max = 34.0E5;
// double pr(0);
// for (int i = 0; i < 1000; ++i){
// // pr = (i/10.0 + 1) * 1.0E5 ;
// pr = min + i * (max-min)/1000.0;
// // std::cout << " thp = " << pr[i] << ", mismatch = " << mismatch(pr[i]) << std::endl;
// std::cout << pr << ", " << mismatch(pr) << std::endl;
// }
// return;
// }
for (auto& well : this->well_container_) {
std::string well_name = well->name();
if (group.hasWell(well_name)) {
@@ -1380,16 +1406,19 @@ namespace Opm {
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, group_state, local_deferredLogger);
}
// 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;
}
}
//PJPE: use the common THP in computeNetworkPressures
well_group_thp_prev_ = well_group_thp;
// Use the common THP in computeNetworkPressures
group_state.update_well_group_thp(nodeName, well_group_thp);
}
}
@@ -2055,6 +2084,9 @@ namespace Opm {
// network related
bool more_network_update = false;
if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) {
const double dt = this->ebosSimulator_.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);
const Scalar network_imbalance = comm.max(local_network_imbalance);
const auto& balance = this->schedule()[episodeIdx].network_balance();

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

@@ -875,9 +875,15 @@ 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);
@@ -889,11 +895,6 @@ bruteForceBracket(const std::function<Scalar(const Scalar)>& eq,
eq_low = eq_high;
}
if (bracket_found) {
// std::cout << "BFB: low: " << low << " high: " << high << std::endl;
// std::cout << "BFB: low_value: " << eq_low << " high_value: " << eq_high << std::endl;
// double low_value = eq(low);
// double high_value = eq(high);
// std::cout << "CHK: low_value: " << low_value<< " high_value: " << high_value << std::endl;
deferred_logger.debug(
" brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) +
" high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high));
@@ -1010,6 +1011,47 @@ getFloIPR(const WellState<Scalar>& well_state,
detail::getFlo(table, aqua_b, liquid_b, vapour_b));
}
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,
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 double interval = (high - low) / sample_number;
double eq_low = eq(low);
double eq_high = 0.0;
for (int i = 0; i < sample_number + 1; ++i) {
high = range[0] + interval * i;
eq_high = eq(high);
if ( (std::fabs(eq_high) < limit)) {
// approximate_solution_found = true;
approximate_solution = high;
break;
}
if (eq_high * eq_low <= 0.) {
bracket_found = true;
break;
}
low = high;
eq_low = eq_high;
}
if (bracket_found) {
deferred_logger.debug(
" brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) +
" high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high));
}
return bracket_found;
}
template class WellBhpThpCalculator<double>;
#define INSTANCE(...) \

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

@@ -117,18 +117,21 @@ public:
std::pair<Scalar, Scalar>
getFloIPR(const WellState<Scalar>& well_state,
const Well& well,
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,
DeferredLogger& deferred_logger);
const SummaryState& summary_state) const;
//! \brief Find limits using brute-force solver.
static bool bruteForceBracketHigh(const std::function<double(const double)>& eq,
const std::array<double, 2>& range,
double& low, double& high,
DeferredLogger& deferred_logger);
//! \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,
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,
DeferredLogger& deferred_logger);
private:
//! \brief Compute BHP from THP limit for an injector - implementation.

2
opm/simulators/wells/WellGroupHelpers.cpp
View File

@@ -935,7 +935,7 @@ computeNetworkPressures(const Network::ExtNetwork& network,
OpmLog::debug(oss.str());
#endif
} else {
// Table number specified as 9999 in the deck, no pressure loss.
// Table number specified as 9999 in the deck, no pressure loss.
if (network.node(node).as_choke()){
// Node pressure is set to the common THP of the wells.
// The choke pressure must be non-negative therefore the node pressure of

13
opm/simulators/wells/WellInterface.hpp
View File

@@ -382,7 +382,18 @@ public:
const WellProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
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_;

7
opm/simulators/wells/WellInterfaceGeneric.cpp
View File

@@ -487,6 +487,13 @@ WellInterfaceGeneric<Scalar>::getDynamicThpLimit() const
return dynamic_thp_limit_;
}
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
setDynamicThpLimit(const std::optional<Scalar> thp_limit)
{
dynamic_thp_limit_ = thp_limit;
}
template<class Scalar>
void WellInterfaceGeneric<Scalar>::
updatePerforatedCell(std::vector<bool>& is_cell_perforated)

1
opm/simulators/wells/WellInterfaceGeneric.hpp
View File

@@ -103,6 +103,7 @@ public:
void setWsolvent(const Scalar wsolvent);
void setDynamicThpLimit(const Scalar thp_limit);
std::optional<Scalar> getDynamicThpLimit() const;
void setDynamicThpLimit(const std::optional<Scalar> thp_limit);
void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
/// Returns true if the well has one or more THP limits/constraints.