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

add StandardWellGeneric

Browse Source 
avoid rebuilding this for all simulators when code
is only dependent on Scalar. specialized for double
This commit is contained in:
Arne Morten Kvarving 2021-06-03 15:34:14 +02:00
parent 6c1ca7450f
commit 10ff86af52
6 changed files with 1096 additions and 899 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -62,6 +62,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/GroupState.cpp
opm/simulators/wells/ParallelWellInfo.cpp
opm/simulators/wells/SegmentState.cpp
opm/simulators/wells/StandardWellGeneric.cpp
opm/simulators/wells/TargetCalculator.cpp
opm/simulators/wells/VFPHelpers.cpp
opm/simulators/wells/VFPProdProperties.cpp

81
opm/simulators/wells/StandardWell.hpp
View File

@ -29,9 +29,9 @@
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/StandardWellGeneric.hpp>
#include <opm/simulators/wells/VFPInjProperties.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
@ -58,7 +58,8 @@ namespace Opm
{
template<typename TypeTag>
class StandardWell: public WellInterface<TypeTag>
class StandardWell : public WellInterface<TypeTag>
, public StandardWellGeneric<GetPropType<TypeTag, Properties::Scalar>>
{
public:
@ -136,22 +137,6 @@ namespace Opm
using typename Base::BVector;
using typename Base::Eval;
// sparsity pattern for the matrices
//[A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
typedef Dune::DynamicVector<Scalar> VectorBlockWellType;
typedef Dune::BlockVector<VectorBlockWellType> BVectorWell;
// the matrix type for the diagonal matrix D
typedef Dune::DynamicMatrix<Scalar> DiagMatrixBlockWellType;
typedef Dune::BCRSMatrix <DiagMatrixBlockWellType> DiagMatWell;
// the matrix type for the non-diagonal matrix B and C^T
typedef Dune::DynamicMatrix<Scalar> OffDiagMatrixBlockWellType;
typedef Dune::BCRSMatrix<OffDiagMatrixBlockWellType> OffDiagMatWell;
typedef DenseAd::DynamicEvaluation<Scalar, numStaticWellEq + numEq + 1> EvalWell;
using Base::contiSolventEqIdx;
@ -200,9 +185,6 @@ namespace Opm
#if HAVE_CUDA || HAVE_OPENCL
/// add the contribution (C, D^-1, B matrices) of this Well to the WellContributions object
void addWellContribution(WellContributions& wellContribs) const;
/// get the number of blocks of the C and B matrices, used to allocate memory in a WellContributions object
void getNumBlocks(unsigned int& _nnzs) const;
#endif
/// using the solution x to recover the solution xw for wells and applying
@ -247,12 +229,6 @@ namespace Opm
return param_.matrix_add_well_contributions_;
}
bool doGasLiftOptimize(
const WellState& well_state,
const Simulator& ebosSimulator,
DeferredLogger& deferred_logger
) const;
virtual void gasLiftOptimizationStage1 (
WellState& well_state,
const Simulator& ebosSimulator,
@ -262,15 +238,6 @@ namespace Opm
GLiftWellStateMap &state_map
) const override;
bool checkGliftNewtonIterationIdxOk(
const Simulator& ebosSimulator,
DeferredLogger& deferred_logger
) const;
void gliftDebug(
const std::string &msg,
DeferredLogger& deferred_logger) const;
/* returns BHP */
double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
const SummaryState &summary_state,
@ -364,33 +331,13 @@ namespace Opm
using Base::ipr_a_;
using Base::ipr_b_;
using Base::changed_to_stopped_this_step_;
using typename StandardWellGeneric<Scalar>::BVectorWell;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
// densities of the fluid in each perforation
std::vector<double> perf_densities_;
// pressure drop between different perforations
std::vector<double> perf_pressure_diffs_;
// residuals of the well equations
BVectorWell resWell_;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
DiagMatWell invDuneD_;
// Wrapper for the parallel application of B for distributed wells
wellhelpers::ParallelStandardWellB<Scalar> parallelB_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
// the values for the primary varibles
// based on different solutioin strategies, the wells can have different primary variables
mutable std::vector<double> primary_variables_;
@ -401,9 +348,6 @@ namespace Opm
// the saturations in the well bore under surface conditions at the beginning of the time step
std::vector<double> F0_;
// Enable GLIFT debug mode. This will enable output of logging messages.
bool glift_debug = false;
// Optimize only wells under THP control
bool glift_optimize_only_thp_wells = true;
@ -448,8 +392,6 @@ namespace Opm
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
void computeConnectionPressureDelta();
void computeWellConnectionDensitesPressures(const Simulator& ebosSimulator,
const WellState& well_state,
const std::vector<double>& b_perf,
@ -485,8 +427,6 @@ namespace Opm
const WellState &well_state) const;
double calculateThpFromBhp(const WellState &well_state, const std::vector<double>& rates, const double bhp, DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const override;
// get the mobility for specific perforation
@ -570,19 +510,11 @@ namespace Opm
// TODO: looking for better alternative to avoid wrong-signed well rates
bool openCrossFlowAvoidSingularity(const Simulator& ebos_simulator) const;
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
const double dx);
// calculate a relaxation factor to avoid overshoot of the fractions for producers
// which might result in negative rates
static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// calculate the skin pressure based on water velocity, throughput and polymer concentration.
// throughput is used to describe the formation damage during water/polymer injection.
// calculated skin pressure will be applied to the drawdown during perforation rate calculation
@ -616,11 +548,6 @@ namespace Opm
virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
// checking the convergence of the well control equations
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
// checking convergence of extra equations, if there are any
void checkConvergenceExtraEqs(const std::vector<double>& res,
ConvergenceReport& report) const;

824
opm/simulators/wells/StandardWellGeneric.cpp Normal file
View File

@ -0,0 +1,824 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS AS.
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/StandardWellGeneric.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GasLiftOpt.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/VFPInjTable.hpp>
#include <opm/simulators/timestepping/ConvergenceReport.hpp>
#include <opm/simulators/utils/DeferredLogger.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 <fmt/format.h>
#include <stdexcept>
namespace Opm
{
template<class Scalar>
StandardWellGeneric<Scalar>::
StandardWellGeneric(int Bhp,
const ParallelWellInfo& pw_info,
const WellInterfaceGeneric& baseif)
: baseif_(baseif)
, perf_densities_(baseif_.numPerfs())
, perf_pressure_diffs_(baseif_.numPerfs())
, parallelB_(duneB_, pw_info)
, Bhp_(Bhp)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
invDuneD_.setBuildMode(DiagMatWell::row_wise);
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells)
{
double relaxation_factor = 1.0;
static constexpr int WQTotal = 0;
// For injector, we only check the total rates to avoid sign change of rates
const double original_total_rate = primary_variables[WQTotal];
const double newton_update = dwells[0][WQTotal];
const double possible_update_total_rate = primary_variables[WQTotal] - newton_update;
// 0.8 here is a experimental value, which remains to be optimized
// if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will
// always be 1.0, more thoughts might be needed.
if (original_total_rate * possible_update_total_rate < 0.) { // sign changed
relaxation_factor = std::abs(original_total_rate / newton_update) * 0.8;
}
assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
return relaxation_factor;
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
relaxationFactorFraction(const double old_value,
const double dx)
{
assert(old_value >= 0. && old_value <= 1.0);
double relaxation_factor = 1.;
// updated values without relaxation factor
const double possible_updated_value = old_value - dx;
// 0.95 is an experimental value remains to be optimized
if (possible_updated_value < 0.0) {
relaxation_factor = std::abs(old_value / dx) * 0.95;
} else if (possible_updated_value > 1.0) {
relaxation_factor = std::abs((1. - old_value) / dx) * 0.95;
}
// if possible_updated_value is between 0. and 1.0, then relaxation_factor
// remains to be one
assert(relaxation_factor >= 0. && relaxation_factor <= 1.);
return relaxation_factor;
}
template<class Scalar>
double
StandardWellGeneric<Scalar>::
calculateThpFromBhp(const WellState &well_state,
const std::vector<double>& rates,
const double bhp,
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];
// pick the density in the top layer
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, getRho(), 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_.getALQ(well_state);
const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth();
const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, getRho(), 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<class Scalar>
void
StandardWellGeneric<Scalar>::
computeConnectionPressureDelta()
{
// Algorithm:
// We'll assume the perforations are given in order from top to
// bottom for each well. By top and bottom we do not necessarily
// mean in a geometric sense (depth), but in a topological sense:
// the 'top' perforation is nearest to the surface topologically.
// Our goal is to compute a pressure delta for each perforation.
// 1. Compute pressure differences between perforations.
// dp_perf will contain the pressure difference between a
// perforation and the one above it, except for the first
// perforation for each well, for which it will be the
// difference to the reference (bhp) depth.
const int nperf = baseif_.numPerfs();
perf_pressure_diffs_.resize(nperf, 0.0);
auto z_above = baseif_.parallelWellInfo().communicateAboveValues(baseif_.refDepth(), baseif_.perfDepth());
for (int perf = 0; perf < nperf; ++perf) {
const double dz = baseif_.perfDepth()[perf] - z_above[perf];
perf_pressure_diffs_[perf] = dz * perf_densities_[perf] * baseif_.gravity();
}
// 2. Compute pressure differences to the reference point (bhp) by
// accumulating the already computed adjacent pressure
// differences, storing the result in dp_perf.
// This accumulation must be done per well.
const auto beg = perf_pressure_diffs_.begin();
const auto end = perf_pressure_diffs_.end();
baseif_.parallelWellInfo().partialSumPerfValues(beg, end);
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::
gliftDebug(const std::string &msg,
DeferredLogger& deferred_logger) const
{
if (glift_debug) {
const std::string message = fmt::format(
" GLIFT (DEBUG) : SW : Well {} : {}", baseif_.name(), msg);
deferred_logger.info(message);
}
}
template<class Scalar>
bool
StandardWellGeneric<Scalar>::
checkGliftNewtonIterationIdxOk(const int report_step_idx,
const int iteration_idx,
const Schedule& schedule,
DeferredLogger& deferred_logger ) const
{
const GasLiftOpt& glo = schedule.glo(report_step_idx);
if (glo.all_newton()) {
const int nupcol = schedule[report_step_idx].nupcol();
if (this->glift_debug) {
const std::string msg = fmt::format(
"LIFTOPT item4 == YES, it = {}, nupcol = {} --> GLIFT optimize = {}",
iteration_idx,
nupcol,
((iteration_idx <= nupcol) ? "TRUE" : "FALSE"));
gliftDebug(msg, deferred_logger);
}
return iteration_idx <= nupcol;
}
else {
if (this->glift_debug) {
const std::string msg = fmt::format(
"LIFTOPT item4 == NO, it = {} --> GLIFT optimize = {}",
iteration_idx, ((iteration_idx == 1) ? "TRUE" : "FALSE"));
gliftDebug(msg, deferred_logger);
}
return iteration_idx == 1;
}
}
/* At this point we know that the well does not have BHP control mode and
that it does have THP constraints, see computeWellPotentials().
* TODO: Currently we limit the application of gas lift optimization to wells
* operating under THP control mode, does it make sense to
* extend it to other modes?
*/
template<class Scalar>
bool
StandardWellGeneric<Scalar>::
doGasLiftOptimize(const WellState &well_state,
const int report_step_idx,
const int iteration_idx,
const Schedule& schedule,
DeferredLogger& deferred_logger) const
{
gliftDebug("checking if GLIFT should be done..", deferred_logger);
if (!well_state.gliftOptimizationEnabled()) {
gliftDebug("Optimization disabled in WellState", deferred_logger);
return false;
}
if (well_state.gliftCheckAlqOscillation(baseif_.name())) {
gliftDebug("further optimization skipped due to oscillation in ALQ",
deferred_logger);
return false;
}
if (glift_optimize_only_thp_wells) {
const int well_index = baseif_.indexOfWell();
auto control_mode = well_state.currentProductionControl(well_index);
if (control_mode != Well::ProducerCMode::THP ) {
gliftDebug("Not THP control", deferred_logger);
return false;
}
}
if (!checkGliftNewtonIterationIdxOk(report_step_idx,
iteration_idx,
schedule,
deferred_logger)) {
return false;
}
const GasLiftOpt& glo = schedule.glo(report_step_idx);
if (!glo.has_well(baseif_.name())) {
gliftDebug("Gas Lift not activated: WLIFTOPT is probably missing",
deferred_logger);
return false;
}
auto increment = glo.gaslift_increment();
// NOTE: According to the manual: LIFTOPT, item 1, :
// "Increment size for lift gas injection rate. Lift gas is
// allocated to individual wells in whole numbers of the increment
// size. If gas lift optimization is no longer required, it can be
// turned off by entering a zero or negative number."
if (increment <= 0) {
if (glift_debug) {
const std::string msg = fmt::format(
"Gas Lift switched off in LIFTOPT item 1 due to non-positive "
"value: {}", increment);
gliftDebug(msg, deferred_logger);
}
return false;
}
else {
return true;
}
}
template<class Scalar>
std::optional<double>
StandardWellGeneric<Scalar>::
computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) 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().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, getRho(), baseif_.gravity());
auto fbhp = [this, &controls, thp_limit, dp, alq_value](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, 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]);
};
// Get the flo samples, add extra samples at low rates and bhp
// limit point if necessary. Then the sign must be flipped
// since the VFP code expects that production flo values are
// negative.
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);
}
for (double& x : flo_samples) {
x = -x;
}
// 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 under 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 = 600.0 * unit::barsa;
const int max_iteration = 50;
const double flo_tolerance = 1e-6 * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<>::
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() ? high : 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 = 50;
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<>::
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<class Scalar>
std::optional<double>
StandardWellGeneric<Scalar>::
computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
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, getRho(), 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 = 50;
const double flo_tolerance = 1e-6 * std::fabs(flo_samples.back());
int iteration = 0;
try {
const double solved_bhp = RegulaFalsiBisection<>::
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 = 50;
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<>::
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<class Scalar>
void
StandardWellGeneric<Scalar>::
checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger,
const double max_residual_allowed) const
{
double control_tolerance = 0.;
using CR = ConvergenceReport;
CR::WellFailure::Type ctrltype = CR::WellFailure::Type::Invalid;
const int well_index = baseif_.indexOfWell();
if (baseif_.wellIsStopped()) {
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; // use smaller tolerance for zero control?
}
else if (baseif_.isInjector() )
{
auto current = well_state.currentInjectionControl(well_index);
switch(current) {
case Well::InjectorCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case Well::InjectorCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case Well::InjectorCMode::RATE:
case Well::InjectorCMode::RESV:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; //
break;
case Well::InjectorCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; //
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
else if (baseif_.isProducer() )
{
auto current = well_state.currentProductionControl(well_index);
switch(current) {
case Well::ProducerCMode::THP:
ctrltype = CR::WellFailure::Type::ControlTHP;
control_tolerance = 1.e4; // 0.1 bar
break;
case Well::ProducerCMode::BHP:
ctrltype = CR::WellFailure::Type::ControlBHP;
control_tolerance = 1.e3; // 0.01 bar
break;
case Well::ProducerCMode::ORAT:
case Well::ProducerCMode::WRAT:
case Well::ProducerCMode::GRAT:
case Well::ProducerCMode::LRAT:
case Well::ProducerCMode::RESV:
case Well::ProducerCMode::CRAT:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-4; // smaller tolerance for rate control
break;
case Well::ProducerCMode::GRUP:
ctrltype = CR::WellFailure::Type::ControlRate;
control_tolerance = 1.e-6; // smaller tolerance for rate control
break;
default:
OPM_DEFLOG_THROW(std::runtime_error, "Unknown well control control types for well " << baseif_.name(), deferred_logger);
}
}
const double well_control_residual = std::abs(this->resWell_[0][Bhp_]);
const int dummy_component = -1;
if (std::isnan(well_control_residual)) {
report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (well_control_residual > max_residual_allowed * 10.) {
report.setWellFailed({ctrltype, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if ( well_control_residual > control_tolerance) {
report.setWellFailed({ctrltype, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
template<class Scalar>
void
StandardWellGeneric<Scalar>::
checkConvergencePolyMW(const std::vector<double>& res,
ConvergenceReport& report,
const double maxResidualAllowed) const
{
if (baseif_.isInjector()) {
// checking the convergence of the perforation rates
const double 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 < baseif_.numPerfs(); ++perf) {
const double wat_vel_residual = res[Bhp_ + 1 + perf];
if (std::isnan(wat_vel_residual)) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (wat_vel_residual > maxResidualAllowed * 10.) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if (wat_vel_residual > wat_vel_tol) {
report.setWellFailed({wat_vel_failure_type, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
// checking the convergence of the skin pressure
const double pskin_tol = 1000.; // 1000 pascal
const auto pskin_failure_type = CR::WellFailure::Type::Pressure;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
const double pskin_residual = res[Bhp_ + 1 + perf + baseif_.numPerfs()];
if (std::isnan(pskin_residual)) {
report.setWellFailed({pskin_failure_type, CR::Severity::NotANumber, dummy_component, baseif_.name()});
} else if (pskin_residual > maxResidualAllowed * 10.) {
report.setWellFailed({pskin_failure_type, CR::Severity::TooLarge, dummy_component, baseif_.name()});
} else if (pskin_residual > pskin_tol) {
report.setWellFailed({pskin_failure_type, CR::Severity::Normal, dummy_component, baseif_.name()});
}
}
}
}
#if HAVE_CUDA || HAVE_OPENCL
template<class Scalar>
void
StandardWellGeneric<Scalar>::
getNumBlocks(unsigned int& numBlocks) const
{
numBlocks = duneB_.nonzeroes();
}
#endif
template class StandardWellGeneric<double>;
}

160
opm/simulators/wells/StandardWellGeneric.hpp Normal file
View File

@ -0,0 +1,160 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2016 - 2017 IRIS AS.
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_STANDARDWELL_GENERIC_HEADER_INCLUDED
#define OPM_STANDARDWELL_GENERIC_HEADER_INCLUDED
#include <dune/common/dynmatrix.hh>
#include <dune/common/dynvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <optional>
#include <vector>
namespace Opm
{
class ConvergenceReport;
class DeferredLogger;
class ParallelWellInfo;
class Schedule;
class SummaryState;
class WellInterfaceGeneric;
class WellState;
template<class Scalar>
class StandardWellGeneric
{
protected:
// sparsity pattern for the matrices
//[A C^T [x = [ res
// B D ] x_well] res_well]
// the vector type for the res_well and x_well
using VectorBlockWellType = Dune::DynamicVector<Scalar>;
using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
// the matrix type for the diagonal matrix D
using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
// the matrix type for the non-diagonal matrix B and C^T
using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
public:
#if HAVE_CUDA || HAVE_OPENCL
/// get the number of blocks of the C and B matrices, used to allocate memory in a WellContributions object
void getNumBlocks(unsigned int& _nnzs) const;
#endif
protected:
StandardWellGeneric(int Bhp,
const ParallelWellInfo& pw_info,
const WellInterfaceGeneric& baseif);
// calculate a relaxation factor to avoid overshoot of total rates
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
const double dx);
double calculateThpFromBhp(const WellState& well_state,
const std::vector<double>& rates,
const double bhp,
DeferredLogger& deferred_logger) const;
// checking the convergence of the well control equations
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger,
const double max_residual_allowed) const;
void checkConvergencePolyMW(const std::vector<double>& res,
ConvergenceReport& report,
const double maxResidualAllowed) const;
void computeConnectionPressureDelta();
std::optional<double> computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const;
std::optional<double> computeBhpAtThpLimitProdWithAlq(const std::function<std::vector<double>(const double)>& frates,
const SummaryState& summary_state,
DeferredLogger& deferred_logger,
double alq_value) const;
void gliftDebug(const std::string &msg,
DeferredLogger& deferred_logger) const;
bool checkGliftNewtonIterationIdxOk(const int report_step_idx,
const int iteration_idx,
const Schedule& schedule,
DeferredLogger& deferred_logger) const;
bool doGasLiftOptimize(const WellState &well_state,
const int report_step_idx,
const int iteration_idx,
const Schedule& schedule,
DeferredLogger& deferred_logger) const;
// Base interface reference
const WellInterfaceGeneric& baseif_;
// residuals of the well equations
BVectorWell resWell_;
// densities of the fluid in each perforation
std::vector<double> perf_densities_;
// pressure drop between different perforations
std::vector<double> perf_pressure_diffs_;
// Enable GLIFT debug mode. This will enable output of logging messages.
bool glift_debug = false;
// Optimize only wells under THP control
bool glift_optimize_only_thp_wells = true;
// two off-diagonal matrices
OffDiagMatWell duneB_;
OffDiagMatWell duneC_;
// diagonal matrix for the well
DiagMatWell invDuneD_;
// Wrapper for the parallel application of B for distributed wells
wellhelpers::ParallelStandardWellB<Scalar> parallelB_;
// several vector used in the matrix calculation
mutable BVectorWell Bx_;
mutable BVectorWell invDrw_;
double getRho() const { return perf_densities_[0]; }
private:
int Bhp_; // index of Bhp
};
}
#endif // OPM_STANDARDWELL_GENERIC_HEADER_INCLUDED

913
opm/simulators/wells/StandardWell_impl.hpp
View File

File diff suppressed because it is too large Load Diff

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

@ -123,16 +123,28 @@ public:
return number_of_phases_;
}
int numPerfs() const {
return number_of_perforations_;
}
double refDepth() const {
return ref_depth_;
}
double gravity() const {
return gravity_;
}
const VFPProperties* vfpProperties() const {
return vfp_properties_;
}
double gravity() const {
return gravity_;
const ParallelWellInfo& parallelWellInfo() const {
return parallel_well_info_;
}
const std::vector<double>& perfDepth() const {
return perf_depth_;
}
double getTHPConstraint(const SummaryState& summaryState) const;