OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-27 09:40:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3351 from akva2/msw_generic

Browse Source 
Add MultisegmentWellGeneric
This commit is contained in:
Bård Skaflestad 2021-06-08 13:06:01 +02:00 committed by GitHub
commit 8ad44d046d
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 852 additions and 713 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -60,6 +60,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/GasLiftStage2.cpp
opm/simulators/wells/GlobalWellInfo.cpp
opm/simulators/wells/GroupState.cpp
opm/simulators/wells/MultisegmentWellGeneric.cpp
opm/simulators/wells/ParallelWellInfo.cpp
opm/simulators/wells/SegmentState.cpp
opm/simulators/wells/StandardWellGeneric.cpp

48
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -23,6 +23,7 @@
#define OPM_MULTISEGMENTWELL_HEADER_INCLUDED
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
@ -31,7 +32,8 @@ namespace Opm
class DeferredLogger;
template<typename TypeTag>
class MultisegmentWell: public WellInterface<TypeTag>
class MultisegmentWell : public WellInterface<TypeTag>
, public MultisegmentWellGeneric<GetPropType<TypeTag, Properties::Scalar>>
{
public:
typedef WellInterface<TypeTag> Base;
@ -195,10 +197,6 @@ namespace Opm
virtual void addWellContributions(SparseMatrixAdapter& jacobian) const override;
/// number of segments for this well
/// int number_of_segments_;
int numberOfSegments() const;
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const override;
@ -225,9 +223,6 @@ namespace Opm
// multi-phase flow model
WellSegments::MultiPhaseModel multiphaseModel() const;
// get the WellSegments from the well_ecl_
const WellSegments& segmentSet() const;
// protected member variables from the Base class
using Base::well_ecl_;
using Base::vfp_properties_;
@ -267,24 +262,6 @@ namespace Opm
using Base::calculateBhpFromThp;
using Base::getALQ;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int> > segment_perforations_;
// the inlet segments for each segment. It is for convinience and efficiency reason
std::vector<std::vector<int> > segment_inlets_;
// segment number is an ID of the segment, it is specified in the deck
// get the loation of the segment with a segment number in the segmentSet
int segmentNumberToIndex(const int segment_number) const;
// TODO, the following should go to a class for computing purpose
// two off-diagonal matrices
mutable OffDiagMatWell duneB_;
@ -312,11 +289,6 @@ namespace Opm
// center depth of the grid block
std::vector<double> cell_perforation_pressure_diffs_;
// depth difference between the segment and the peforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_segment_depth_diffs_;
// the intial amount of fluids in each segment under surface condition
std::vector<std::vector<double> > segment_fluid_initial_;
@ -330,8 +302,6 @@ namespace Opm
// the mass rate of the segments
std::vector<EvalWell> segment_mass_rates_;
std::vector<double> segment_depth_diffs_;
// the upwinding segment for each segment based on the flow direction
std::vector<int> upwinding_segments_;
@ -360,10 +330,6 @@ namespace Opm
const double relaxation_factor=1.0) const;
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
// computing the accumulation term for later use in well mass equations
void computeInitialSegmentFluids(const Simulator& ebos_simulator);
@ -396,7 +362,6 @@ namespace Opm
// convert a Eval from reservoir to contain the derivative related to wells
EvalWell extendEval(const Eval& in) const;
double calculateThpFromBhp(const std::vector<double>& rates, const double bhp, DeferredLogger& deferred_logger) const;
void updateThp(WellState& well_state, DeferredLogger& deferred_logger) const;
// compute the fluid properties, such as densities, viscosities, and so on, in the segments
@ -457,10 +422,6 @@ namespace Opm
virtual double getRefDensity() const override;
bool frictionalPressureLossConsidered() const;
bool accelerationalPressureLossConsidered() const;
virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
@ -484,9 +445,6 @@ namespace Opm
std::vector<Scalar> getWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const;
void detectOscillations(const std::vector<double>& measure_history,
const int it, bool& oscillate, bool& stagnate) const;
double getResidualMeasureValue(const WellState& well_state,
const std::vector<double>& residuals,
DeferredLogger& deferred_logger) const;

643
opm/simulators/wells/MultisegmentWellGeneric.cpp Normal file
View File

@ -0,0 +1,643 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPInjTable.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/VFPHelpers.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <cassert>
#include <cmath>
#include <stdexcept>
namespace Opm
{
template<typename Scalar>
MultisegmentWellGeneric<Scalar>::
MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
: baseif_(baseif)
, segment_perforations_(numberOfSegments())
, segment_inlets_(numberOfSegments())
, segment_depth_diffs_(numberOfSegments(), 0.0)
, perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
// for other facilities needed but not available from parser, we need to process them here
// initialize the segment_perforations_ and update perforation_segment_depth_diffs_
const WellConnections& completion_set = baseif_.wellEcl().getConnections();
// index of the perforation within wells struct
// there might be some perforations not active, which causes the number of the perforations in
// well_ecl_ and wells struct different
// the current implementation is a temporary solution for now, it should be corrected from the parser
// side
int i_perf_wells = 0;
baseif.perfDepth().resize(baseif_.numPerfs(), 0.);
for (size_t perf = 0; perf < completion_set.size(); ++perf) {
const Connection& connection = completion_set.get(perf);
if (connection.state() == Connection::State::OPEN) {
const int segment_index = segmentNumberToIndex(connection.segment());
segment_perforations_[segment_index].push_back(i_perf_wells);
baseif.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segmentSet()[segment_index].depth();
perforation_segment_depth_diffs_[i_perf_wells] = baseif.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
}
// initialize the segment_inlets_
for (int seg = 0; seg < numberOfSegments(); ++seg) {
const Segment& segment = segmentSet()[seg];
const int segment_number = segment.segmentNumber();
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) {
const int segment_index = segmentNumberToIndex(segment_number);
const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
segment_inlets_[outlet_segment_index].push_back(segment_index);
}
}
// calculating the depth difference between the segment and its oulet_segments
// for the top segment, we will make its zero unless we find other purpose to use this value
for (int seg = 1; seg < numberOfSegments(); ++seg) {
const double segment_depth = segmentSet()[seg].depth();
const int outlet_segment_number = segmentSet()[seg].outletSegment();
const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
const double outlet_depth = outlet_segment.depth();
segment_depth_diffs_[seg] = segment_depth - outlet_depth;
}
}
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
scaleSegmentRatesWithWellRates(WellState& well_state) const
{
auto& segments = well_state.segments(baseif_.indexOfWell());
auto& segment_rates = segments.rates;
for (int phase = 0; phase < baseif_.numPhases(); ++phase) {
const double unscaled_top_seg_rate = segment_rates[phase];
const double well_phase_rate = well_state.wellRates(baseif_.indexOfWell())[phase];
if (std::abs(unscaled_top_seg_rate) > 1e-12)
{
for (int seg = 0; seg < numberOfSegments(); ++seg) {
segment_rates[baseif_.numPhases()*seg + phase] *= well_phase_rate/unscaled_top_seg_rate;
}
} else {
// for newly opened wells, the unscaled rate top segment rate is zero
// and we need to initialize the segment rates differently
double sumTw = 0;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
sumTw += baseif_.wellIndex()[perf];
}
std::vector<double> perforation_rates(baseif_.numPhases() * baseif_.numPerfs(),0.0);
const double perf_phaserate_scaled = well_state.wellRates(baseif_.indexOfWell())[phase] / sumTw;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
perforation_rates[baseif_.numPhases()* perf + phase] = baseif_.wellIndex()[perf] * perf_phaserate_scaled;
}
std::vector<double> rates;
WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, baseif_.numPhases(), 0, rates);
std::copy(rates.begin(), rates.end(), segment_rates.begin());
}
}
}
template <typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
scaleSegmentPressuresWithBhp(WellState& well_state) const
{
auto& segments = well_state.segments(baseif_.indexOfWell());
auto bhp = well_state.bhp(baseif_.indexOfWell());
segments.scale_pressure(bhp);
}
template<typename Scalar>
const WellSegments&
MultisegmentWellGeneric<Scalar>::
segmentSet() const
{
return baseif_.wellEcl().getSegments();
}
template <typename Scalar>
int
MultisegmentWellGeneric<Scalar>::
numberOfSegments() const
{
return segmentSet().size();
}
template <typename Scalar>
WellSegments::CompPressureDrop
MultisegmentWellGeneric<Scalar>::
compPressureDrop() const
{
return segmentSet().compPressureDrop();
}
template<typename Scalar>
int
MultisegmentWellGeneric<Scalar>::
segmentNumberToIndex(const int segment_number) const
{
return segmentSet().segmentNumberToIndex(segment_number);
}
template<typename Scalar>
double
MultisegmentWellGeneric<Scalar>::
calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
DeferredLogger& deferred_logger) const
{
assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
static constexpr int Water = BlackoilPhases::Aqua;
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];
double thp = 0.0;
if (baseif_.isInjector()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
thp = baseif_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
}
else if (baseif_.isProducer()) {
const int table_id = baseif_.wellEcl().vfp_table_number();
const double alq = baseif_.wellEcl().alq_value();
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
thp = baseif_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
}
else {
OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
}
return thp;
}
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
detectOscillations(const std::vector<double>& measure_history,
const int it,
bool& oscillate,
bool& stagnate) const
{
if ( it < 2 ) {
oscillate = false;
stagnate = false;
return;
}
stagnate = true;
const double F0 = measure_history[it];
const double F1 = measure_history[it - 1];
const double F2 = measure_history[it - 2];
const double d1 = std::abs((F0 - F2) / F0);
const double d2 = std::abs((F0 - F1) / F0);
const double oscillaton_rel_tol = 0.2;
oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2);
const double stagnation_rel_tol = 1.e-2;
stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol;
}
template<typename Scalar>
std::optional<double>
MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) is returned.
//
// In order to detect these situations, we will find piecewise
// linear approximations both to the inverse of the frates
// function and to the fbhp function.
//
// We first take the FLO sample points of the VFP curve, and
// find the corresponding bhp values by solving the equation:
// flo(frates(bhp)) - flo_sample = 0
// for bhp, for each flo_sample. The resulting (flo_sample,
// bhp_sample) values give a piecewise linear approximation to
// the true inverse inflow function, at the same flo values as
// the VFP data.
//
// Then we extract a piecewise linear approximation from the
// multilinear fbhp() by evaluating it at the flo_sample
// points, with fractions given by the frates(bhp_sample)
// values.
//
// When we have both piecewise linear curves defined on the
// same flo_sample points, it is easy to distinguish between
// the 0, 1 or 2 solution cases, and obtain the right interval
// in which to solve for the solution we want (with highest
// flow in case of 2 solutions).
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = baseif_.wellEcl().injectionControls(summary_state);
const auto& table = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = baseif_.getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
assert(rates.size() == 3);
return baseif_.vfpProperties()->getInj()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& 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();
if (flo_samples[0] > 0.0) {
const double 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));
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) {
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
// with simply the bhp limit. The first one
// encountered is considered valid, the rest not. They
// are therefore skipped.
bhp_samples.push_back(controls.bhp_limit);
break;
}
auto eq = [&flo, &frates, flo_sample](double 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 int max_iteration = 100;
const double flo_tolerance = 0.05 * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
solve(eq, low, high, max_iteration, flo_tolerance, iteration);
bhp_samples.push_back(solved_bhp);
}
catch (...) {
// Use previous value (or max value if at start) if we failed.
bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back());
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES",
"Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + baseif_.name());
}
}
// 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);
for (int ii = 0; ii < num_samples; ++ii) {
fbhp_samples[ii] = fbhp(frates(bhp_samples[ii]));
}
// #define EXTRA_THP_DEBUGGING
#ifdef EXTRA_THP_DEBUGGING
std::string dbgmsg;
dbgmsg += "flo: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(flo_samples[ii]);
}
dbgmsg += "\nbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(bhp_samples[ii]);
}
dbgmsg += "\nfbhp: ";
for (int ii = 0; ii < num_samples; ++ii) {
dbgmsg += " " + std::to_string(fbhp_samples[ii]);
}
OpmLog::debug(dbgmsg);
#endif // EXTRA_THP_DEBUGGING
// Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve.
// 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];
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.
}
}
// Handle the no solution case.
if (sign_change_index == -1) {
return std::optional<double>();
}
// Solve for the proper solution in the given interval.
auto eq = [&fbhp, &frates](double 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 int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
if (low == high) {
// We are in the high flow regime where the bhp_samples
// are all equal to the bhp_limit.
assert(low == controls.bhp_limit);
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::optional<double>();
}
try {
const double solved_bhp = RegulaFalsiBisection<WarnAndContinueOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
#ifdef EXTRA_THP_DEBUGGING
OpmLog::debug("***** " + name() + " solved_bhp = " + std::to_string(solved_bhp)
+ " flo_bhp_limit = " + std::to_string(flo_bhp_limit));
#endif // EXTRA_THP_DEBUGGING
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::optional<double>();
}
}
template<typename Scalar>
std::optional<double>
MultisegmentWellGeneric<Scalar>::
computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
DeferredLogger& deferred_logger) const
{
// Given a VFP function returning bhp as a function of phase
// rates and thp:
// fbhp(rates, thp),
// a function extracting the particular flow rate used for VFP
// lookups:
// flo(rates)
// and the inflow function (assuming the reservoir is fixed):
// frates(bhp)
// we want to solve the equation:
// fbhp(frates(bhp, thplimit)) - bhp = 0
// for bhp.
//
// This may result in 0, 1 or 2 solutions. If two solutions,
// the one corresponding to the lowest bhp (and therefore
// highest rate) should be returned.
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// Make the fbhp() function.
const auto& controls = baseif_.wellEcl().productionControls(summary_state);
const auto& table = baseif_.vfpProperties()->getProd()->getTable(controls.vfp_table_number);
const double vfp_ref_depth = table.getDatumDepth();
const double thp_limit = baseif_.getTHPConstraint(summary_state);
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
assert(rates.size() == 3);
return baseif_.vfpProperties()->getProd()
->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp;
};
// Make the flo() function.
auto flo = [&table](const std::vector<double>& rates) {
return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]);
};
// Find the bhp-point where production becomes nonzero.
double bhp_max = 0.0;
{
auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); };
double low = controls.bhp_limit;
double high = maxPerfPress + 1.0 * unit::barsa;
double f_low = fflo(low);
double f_high = fflo(high);
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high));
int adjustments = 0;
const int max_adjustments = 10;
const double 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;
f_high = fflo(high);
++adjustments;
}
if (f_low * f_high > 0.0) {
if (f_low > 0.0) {
// Even at the BHP limit, we are injecting.
// There will be no solution here, return an
// empty optional.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE",
"Robust bhp(thp) solve failed due to inoperability for well " + baseif_.name());
return std::optional<double>();
} else {
// Still producing, even at high bhp.
assert(f_high < 0.0);
bhp_max = high;
}
} else {
// Bisect to find a bhp point where we produce, but
// not a large amount ('eps' below).
const double eps = 0.1 * std::fabs(table.getFloAxis().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);
if (f_curr * f_low > 0.0) {
low = curr;
f_low = f_curr;
} else {
high = curr;
f_high = f_curr;
}
++it;
}
bhp_max = low;
}
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" f(low) = " + std::to_string(f_low) +
" f(high) = " + std::to_string(f_high) +
" bhp_max = " + std::to_string(bhp_max));
}
// Define the equation we want to solve.
auto eq = [&fbhp, &frates](double bhp) {
return fbhp(frates(bhp)) - bhp;
};
// Find appropriate brackets for the solution.
double low = controls.bhp_limit;
double high = bhp_max;
{
double eq_high = eq(high);
double eq_low = eq(low);
const double eq_bhplimit = eq_low;
deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() +
" low = " + std::to_string(low) +
" high = " + std::to_string(high) +
" eq(low) = " + std::to_string(eq_low) +
" eq(high) = " + std::to_string(eq_high));
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);
int bracket_attempts = 0;
const int max_bracket_attempts = 20;
double interval = high - low;
const double 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);
eq_low = eq(low);
abs_low = std::fabs(eq_low);
} else {
high = 0.5 * (low + high);
eq_high = eq(high);
abs_high = std::fabs(eq_high);
}
++bracket_attempts;
}
if (eq_low * eq_high > 0.0) {
// Still failed bracketing!
const double limit = 3.0 * unit::barsa;
if (std::min(abs_low, abs_high) < limit) {
// Return the least bad solution if less off than 3 bar.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
"Robust bhp(thp) not solved precisely for well " + baseif_.name());
return abs_low < abs_high ? low : high;
} else {
// Return failure.
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE",
"Robust bhp(thp) solve failed due to bracketing failure for well " + baseif_.name());
return std::optional<double>();
}
}
}
// We have a bracket!
// Now, see if (bhplimit, low) is a bracket in addition to (low, high).
// If so, that is the bracket we shall use, choosing the solution with the
// highest flow.
if (eq_low * eq_bhplimit <= 0.0) {
high = low;
low = controls.bhp_limit;
}
}
// Solve for the proper solution in the given interval.
const int max_iteration = 100;
const double bhp_tolerance = 0.01 * unit::barsa;
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<ThrowOnError>::
solve(eq, low, high, max_iteration, bhp_tolerance, iteration);
return solved_bhp;
}
catch (...) {
deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE",
"Robust bhp(thp) solve failed for well " + baseif_.name());
return std::optional<double>();
}
}
template<typename Scalar>
bool
MultisegmentWellGeneric<Scalar>::
frictionalPressureLossConsidered() const
{
// HF- and HFA needs to consider frictional pressure loss
return (segmentSet().compPressureDrop() != WellSegments::CompPressureDrop::H__);
}
template<typename Scalar>
bool
MultisegmentWellGeneric<Scalar>::
accelerationalPressureLossConsidered() const
{
return (segmentSet().compPressureDrop() == WellSegments::CompPressureDrop::HFA);
}
template class MultisegmentWellGeneric<double>;
} // namespace Opm

113
opm/simulators/wells/MultisegmentWellGeneric.hpp Normal file
View File

@ -0,0 +1,113 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/WellSegments.hpp>
#include <functional>
#include <optional>
#include <vector>
namespace Opm
{
class DeferredLogger;
class SummaryState;
class WellInterfaceGeneric;
class WellState;
template <typename Scalar>
class MultisegmentWellGeneric
{
protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
// get the WellSegments from the well_ecl_
const WellSegments& segmentSet() const;
// components of the pressure drop to be included
WellSegments::CompPressureDrop compPressureDrop() const;
// segment number is an ID of the segment, it is specified in the deck
// get the loation of the segment with a segment number in the segmentSet
int segmentNumberToIndex(const int segment_number) const;
/// number of segments for this well
int numberOfSegments() const;
double calculateThpFromBhp(const std::vector<double>& rates,
const double bhp,
const double rho,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double rho,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
const double maxPerfPress,
const double rho,
DeferredLogger& deferred_logger) const;
/// Detect oscillation or stagnation based on the residual measure history
void detectOscillations(const std::vector<double>& measure_history,
const int it,
bool& oscillate,
bool& stagnate) const;
bool accelerationalPressureLossConsidered() const;
bool frictionalPressureLossConsidered() const;
const WellInterfaceGeneric& baseif_;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int>> segment_perforations_;
// the inlet segments for each segment. It is for convinience and efficiency reason
std::vector<std::vector<int>> segment_inlets_;
std::vector<double> segment_depth_diffs_;
// depth difference between the segment and the perforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_segment_depth_diffs_;
};
}
#endif // OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED

752
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

File diff suppressed because it is too large Load Diff

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

@ -147,6 +147,14 @@ public:
return perf_depth_;
}
std::vector<double>& perfDepth() {
return perf_depth_;
}
const std::vector<double>& wellIndex() const {
return well_index_;
}
double getTHPConstraint(const SummaryState& summaryState) const;
double getALQ(const WellState& well_state) const;