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

Merge pull request #3240 from akva2/wellinterface_fluidsystem_split

Browse Source 
Split out FluidSystem dependent code from WellInterface
This commit is contained in:
Bård Skaflestad
2021-05-27 21:09:07 +02:00
committed by GitHub
parent 21899ab3bb 0f7967f86a
commit b16eaf72f1
6 changed files with 1526 additions and 1454 deletions
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/VFPProdProperties.cpp
opm/simulators/wells/VFPInjProperties.cpp
opm/simulators/wells/WellGroupHelpers.cpp
opm/simulators/wells/WellInterfaceFluidSystem.cpp
opm/simulators/wells/WellInterfaceGeneric.cpp
opm/simulators/wells/WellProdIndexCalculator.cpp
opm/simulators/wells/WellState.cpp

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

@@ -34,7 +34,6 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
@@ -60,451 +59,357 @@ namespace Opm {
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
#include <string>
#include <memory>
#include <vector>
#include <array>
#include <cassert>
#include <memory>
#include <string>
#include <vector>
namespace Opm
{
template<typename TypeTag>
class WellInterface : public WellInterfaceGeneric
{
public:
using ModelParameters = BlackoilModelParametersEbos<TypeTag>;
static const int Water = BlackoilPhases::Aqua;
static const int Oil = BlackoilPhases::Liquid;
static const int Gas = BlackoilPhases::Vapour;
using Grid = GetPropType<TypeTag, Properties::Grid>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using GasLiftSingleWell = ::Opm::GasLiftSingleWell<TypeTag>;
using GLiftOptWells = typename BlackoilWellModel<TypeTag>::GLiftOptWells;
using GLiftProdWells = typename BlackoilWellModel<TypeTag>::GLiftProdWells;
using GLiftWellStateMap =
typename BlackoilWellModel<TypeTag>::GLiftWellStateMap;
static const int numEq = Indices::numEq;
static const int numPhases = Indices::numPhases;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
using MatrixBlockType = Dune::FieldMatrix<Scalar, numEq, numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>;
static constexpr bool has_solvent = getPropValue<TypeTag, Properties::EnableSolvent>();
static constexpr bool has_zFraction = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool has_polymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool has_energy = getPropValue<TypeTag, Properties::EnableEnergy>();
static const bool has_temperature = getPropValue<TypeTag, Properties::EnableTemperature>();
// flag for polymer molecular weight related
static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
// index for the polymer molecular weight continuity equation
static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
// For the conversion between the surface volume rate and reservoir voidage rate
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
static const bool compositionSwitchEnabled = Indices::gasEnabled;
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
has_temperature,
has_energy,
compositionSwitchEnabled,
has_brine,
Indices::numPhases >;
/// Constructor
WellInterface(const Well& well,
const ParallelWellInfo& pw_info,
const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
/// Virtual destructor
virtual ~WellInterface() = default;
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg);
virtual void initPrimaryVariablesEvaluation() const = 0;
virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, DeferredLogger& deferred_logger, const bool relax_tolerance = false) const = 0;
virtual void solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger) = 0;
void assembleWellEq(const Simulator& ebosSimulator,
const double dt,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
virtual void gasLiftOptimizationStage1 (
WellState& well_state,
const Simulator& ebosSimulator,
DeferredLogger& deferred_logger,
GLiftProdWells& prod_wells,
GLiftOptWells& glift_wells,
GLiftWellStateMap& state_map
) const = 0;
void updateWellTestState(const WellState& well_state,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const;
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state,
DeferredLogger& deferred_logger) const = 0;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const = 0;
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const = 0;
// TODO: before we decide to put more information under mutable, this function is not const
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials,
DeferredLogger& deferred_logger) = 0;
void updateWellStateWithTarget(const Simulator& ebos_simulator,
WellState& well_state,
DeferredLogger& deferred_logger) const;
enum class IndividualOrGroup { Individual, Group, Both };
bool updateWellControl(const Simulator& ebos_simulator,
const IndividualOrGroup iog,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) /* const */;
virtual void updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const = 0;
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state,
DeferredLogger& deferred_logger) = 0; // should be const?
virtual void updateProductivityIndex(const Simulator& ebosSimulator,
const WellProdIndexCalculator& wellPICalc,
WellState& well_state,
DeferredLogger& deferred_logger) const = 0;
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const
{
return false;
}
// updating the voidage rates in well_state when requested
void calculateReservoirRates(WellState& well_state) const;
// Add well contributions to matrix
virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
void addCellRates(RateVector& rates, int cellIdx) const;
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;
template <class EvalWell>
Eval restrictEval(const EvalWell& in) const
{
Eval out = 0.0;
out.setValue(in.value());
for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
out.setDerivative(eqIdx, in.derivative(eqIdx));
}
return out;
}
// TODO: theoretically, it should be a const function
// Simulator is not const is because that assembleWellEq is non-const Simulator
void wellTesting(const Simulator& simulator,
const double simulation_time, const int report_step,
const WellTestConfig::Reason testing_reason,
/* const */ WellState& well_state, const GroupState& group_state, WellTestState& welltest_state,
DeferredLogger& deferred_logger);
void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
// check whether the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
void updateWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
DeferredLogger& deferred_logger);
// update perforation water throughput based on solved water rate
virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
/// Compute well rates based on current reservoir conditions and well variables.
/// Used in updateWellStateRates().
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const = 0;
/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
/// to the rates returned by computeCurrentWellRates().
void updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const;
void solveWellEquation(const Simulator& ebosSimulator,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
virtual bool useInnerIterations() const = 0;
protected:
// to indicate a invalid completion
static const int INVALIDCOMPLETION = INT_MAX;
// simulation parameters
const ModelParameters& param_;
// For the conversion between the surface volume rate and resrevoir voidage rate
const RateConverterType& rateConverter_;
std::vector<RateVector> connectionRates_;
std::vector< Scalar > B_avg_;
bool changed_to_stopped_this_step_ = false;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
double wpolymer() const;
double wfoam() const;
double wsalt() const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const double * rates_or_potentials,
DeferredLogger& deferred_logger) const;
template <class ValueType>
ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const = 0;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
struct RatioLimitCheckReport;
void checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report,
DeferredLogger& deferred_logger) const;
template <typename RatioFunc>
bool checkMaxRatioLimitWell(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc) const;
template <typename RatioFunc>
void checkMaxRatioLimitCompletions(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const;
double scalingFactor(const int comp_idx) const;
std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const std::vector<double>& potentials) const;
// check whether the well is operable under BHP limit with current reservoir condition
virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) =0;
// check whether the well is operable under THP limit with current reservoir condition
virtual void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger) =0;
virtual void updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const=0;
void wellTestingEconomic(const Simulator& simulator,
const double simulation_time, const WellState& well_state, const GroupState& group_state,
WellTestState& welltest_state, DeferredLogger& deferred_logger);
void wellTestingPhysical(const Simulator& simulator,
const double simulation_time, const int report_step,
WellState& well_state,
const GroupState& group_state,
WellTestState& welltest_state, DeferredLogger& deferred_logger);
virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
// iterate well equations with the specified control until converged
virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
bool iterateWellEquations(const Simulator& ebosSimulator,
const double dt,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
void updateWellTestStateEconomic(const WellState& well_state,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const;
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
DeferredLogger& deferred_logger);
bool checkConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
bool checkIndividualConstraints(WellState& well_state,
const SummaryState& summaryState) const;
bool checkGroupConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
template<typename TypeTag>
class WellInterface : public WellInterfaceFluidSystem<GetPropType<TypeTag, Properties::FluidSystem>>
{
public:
using ModelParameters = BlackoilModelParametersEbos<TypeTag>;
using Grid = GetPropType<TypeTag, Properties::Grid>;
using Simulator = GetPropType<TypeTag, Properties::Simulator>;
using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
using Indices = GetPropType<TypeTag, Properties::Indices>;
using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
using RateVector = GetPropType<TypeTag, Properties::RateVector>;
using GasLiftSingleWell = ::Opm::GasLiftSingleWell<TypeTag>;
using GLiftOptWells = typename BlackoilWellModel<TypeTag>::GLiftOptWells;
using GLiftProdWells = typename BlackoilWellModel<TypeTag>::GLiftProdWells;
using GLiftWellStateMap =
typename BlackoilWellModel<TypeTag>::GLiftWellStateMap;
static const int numEq = Indices::numEq;
static const int numPhases = Indices::numPhases;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using VectorBlockType = Dune::FieldVector<Scalar, numEq>;
using MatrixBlockType = Dune::FieldMatrix<Scalar, numEq, numEq>;
using BVector = Dune::BlockVector<VectorBlockType>;
using Eval = DenseAd::Evaluation<Scalar, /*size=*/numEq>;
using RateConverterType =
typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType;
using WellInterfaceFluidSystem<FluidSystem>::Gas;
using WellInterfaceFluidSystem<FluidSystem>::Oil;
using WellInterfaceFluidSystem<FluidSystem>::Water;
using RatioLimitCheckReport = typename WellInterfaceFluidSystem<FluidSystem>::RatioLimitCheckReport;
static constexpr bool has_solvent = getPropValue<TypeTag, Properties::EnableSolvent>();
static constexpr bool has_zFraction = getPropValue<TypeTag, Properties::EnableExtbo>();
static constexpr bool has_polymer = getPropValue<TypeTag, Properties::EnablePolymer>();
static constexpr bool has_energy = getPropValue<TypeTag, Properties::EnableEnergy>();
static const bool has_temperature = getPropValue<TypeTag, Properties::EnableTemperature>();
// flag for polymer molecular weight related
static constexpr bool has_polymermw = getPropValue<TypeTag, Properties::EnablePolymerMW>();
static constexpr bool has_foam = getPropValue<TypeTag, Properties::EnableFoam>();
static constexpr bool has_brine = getPropValue<TypeTag, Properties::EnableBrine>();
static const int contiSolventEqIdx = Indices::contiSolventEqIdx;
static const int contiZfracEqIdx = Indices::contiZfracEqIdx;
static const int contiPolymerEqIdx = Indices::contiPolymerEqIdx;
// index for the polymer molecular weight continuity equation
static const int contiPolymerMWEqIdx = Indices::contiPolymerMWEqIdx;
static const int contiFoamEqIdx = Indices::contiFoamEqIdx;
static const int contiBrineEqIdx = Indices::contiBrineEqIdx;
// For the conversion between the surface volume rate and reservoir voidage rate
static const bool compositionSwitchEnabled = Indices::gasEnabled;
using FluidState = BlackOilFluidState<Eval,
FluidSystem,
has_temperature,
has_energy,
compositionSwitchEnabled,
has_brine,
Indices::numPhases >;
/// Constructor
WellInterface(const Well& well,
const ParallelWellInfo& pw_info,
const int time_step,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
/// Virtual destructor
virtual ~WellInterface() = default;
virtual void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg);
virtual void initPrimaryVariablesEvaluation() const = 0;
virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, DeferredLogger& deferred_logger, const bool relax_tolerance = false) const = 0;
virtual void solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger) = 0;
void assembleWellEq(const Simulator& ebosSimulator,
const double dt,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
virtual void gasLiftOptimizationStage1 (
WellState& well_state,
const Simulator& ebosSimulator,
DeferredLogger& deferred_logger,
GLiftProdWells& prod_wells,
GLiftOptWells& glift_wells,
GLiftWellStateMap& state_map
) const = 0;
/// using the solution x to recover the solution xw for wells and applying
/// xw to update Well State
virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
WellState& well_state,
DeferredLogger& deferred_logger) const = 0;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const = 0;
/// r = r - C D^-1 Rw
virtual void apply(BVector& r) const = 0;
// TODO: before we decide to put more information under mutable, this function is not const
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const WellState& well_state,
std::vector<double>& well_potentials,
DeferredLogger& deferred_logger) = 0;
void updateWellStateWithTarget(const Simulator& ebos_simulator,
WellState& well_state,
DeferredLogger& deferred_logger) const;
std::pair<bool, double> checkGroupConstraintsProd(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
enum class IndividualOrGroup { Individual, Group, Both };
bool updateWellControl(const Simulator& ebos_simulator,
const IndividualOrGroup iog,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) /* const */;
std::pair<bool, double> checkGroupConstraintsInj(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
virtual void updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const = 0;
template <class EvalWell>
void getGroupInjectionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
EvalWell& control_eq,
double efficiencyFactor,
DeferredLogger& deferred_logger);
virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
const WellState& well_state,
DeferredLogger& deferred_logger) = 0; // should be const?
template <class EvalWell>
void getGroupProductionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
EvalWell& control_eq,
double efficiencyFactor);
virtual void updateProductivityIndex(const Simulator& ebosSimulator,
const WellProdIndexCalculator& wellPICalc,
WellState& well_state,
DeferredLogger& deferred_logger) const = 0;
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqInj(const WellState& well_state,
/// \brief Wether the Jacobian will also have well contributions in it.
virtual bool jacobianContainsWellContributions() const
{
return false;
}
// Add well contributions to matrix
virtual void addWellContributions(SparseMatrixAdapter&) const = 0;
void addCellRates(RateVector& rates, int cellIdx) const;
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;
template <class EvalWell>
Eval restrictEval(const EvalWell& in) const
{
Eval out = 0.0;
out.setValue(in.value());
for(int eqIdx = 0; eqIdx < numEq;++eqIdx) {
out.setDerivative(eqIdx, in.derivative(eqIdx));
}
return out;
}
// TODO: theoretically, it should be a const function
// Simulator is not const is because that assembleWellEq is non-const Simulator
void wellTesting(const Simulator& simulator,
const double simulation_time, const int report_step,
const WellTestConfig::Reason testing_reason,
/* const */ WellState& well_state, const GroupState& group_state, WellTestState& welltest_state,
DeferredLogger& deferred_logger);
void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
// check whether the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
void updateWellOperability(const Simulator& ebos_simulator,
const WellState& well_state,
DeferredLogger& deferred_logger);
// update perforation water throughput based on solved water rate
virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
/// Compute well rates based on current reservoir conditions and well variables.
/// Used in updateWellStateRates().
virtual std::vector<double> computeCurrentWellRates(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger) const = 0;
/// Modify the well_state's rates if there is only one nonzero rate.
/// If so, that rate is kept as is, but the others are set proportionally
/// to the rates returned by computeCurrentWellRates().
void updateWellStateRates(const Simulator& ebosSimulator,
WellState& well_state,
DeferredLogger& deferred_logger) const;
void solveWellEquation(const Simulator& ebosSimulator,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
virtual bool useInnerIterations() const = 0;
protected:
// simulation parameters
const ModelParameters& param_;
std::vector<RateVector> connectionRates_;
std::vector< Scalar > B_avg_;
bool changed_to_stopped_this_step_ = false;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
double wpolymer() const;
double wfoam() const;
double wsalt() const;
template <class ValueType>
ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const = 0;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
double scalingFactor(const int comp_idx) const;
std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const std::vector<double>& potentials) const;
// check whether the well is operable under BHP limit with current reservoir condition
virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) =0;
// check whether the well is operable under THP limit with current reservoir condition
virtual void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger) =0;
virtual void updateIPR(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const=0;
void wellTestingEconomic(const Simulator& simulator,
const double simulation_time, const WellState& well_state, const GroupState& group_state,
WellTestState& welltest_state, DeferredLogger& deferred_logger);
void wellTestingPhysical(const Simulator& simulator,
const double simulation_time, const int report_step,
WellState& well_state,
const GroupState& group_state,
WellTestState& welltest_state, DeferredLogger& deferred_logger);
virtual void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
// iterate well equations with the specified control until converged
virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger) = 0;
bool iterateWellEquations(const Simulator& ebosSimulator,
const double dt,
WellState& well_state,
const GroupState& group_state,
DeferredLogger& deferred_logger);
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
DeferredLogger& deferred_logger);
template <class EvalWell>
void getGroupInjectionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::InjectionControls& controls,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
double efficiencyFactor,
DeferredLogger& deferred_logger);
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqProd(const WellState& well_state,
template <class EvalWell>
void getGroupProductionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::ProductionControls& controls,
const EvalWell& bhp,
const std::vector<EvalWell>& rates, // Always 3 canonical rates.
BhpFromThpFunc bhp_from_thp,
const std::vector<EvalWell>& rates,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
};
double efficiencyFactor);
template<typename TypeTag>
struct
WellInterface<TypeTag>::
RatioLimitCheckReport{
bool ratio_limit_violated = false;
int worst_offending_completion = INVALIDCOMPLETION;
double violation_extent = 0.0;
};
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqInj(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::InjectionControls& controls,
const EvalWell& bhp,
const EvalWell& injection_rate,
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqProd(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::ProductionControls& controls,
const EvalWell& bhp,
const std::vector<EvalWell>& rates, // Always 3 canonical rates.
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
};
}

906
opm/simulators/wells/WellInterfaceFluidSystem.cpp Normal file
View File

@@ -0,0 +1,906 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2018 IRIS
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/WellInterfaceFluidSystem.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <ebos/eclalternativeblackoilindices.hh>
#include <cassert>
#include <cmath>
namespace Opm
{
template<class FluidSystem>
WellInterfaceFluidSystem<FluidSystem>::
WellInterfaceFluidSystem(const Well& well,
const ParallelWellInfo& parallel_well_info,
const int time_step,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data)
: WellInterfaceGeneric(well, parallel_well_info, time_step,
pvtRegionIdx, num_components, num_phases,
index_of_well, first_perf_index, perf_data)
, rateConverter_(rate_converter)
{
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
calculateReservoirRates(WellState& well_state) const
{
const int fipreg = 0; // not considering the region for now
const int np = number_of_phases_;
std::vector<double> surface_rates(np, 0.0);
for (int p = 0; p < np; ++p) {
surface_rates[p] = well_state.wellRates(index_of_well_)[p];
}
std::vector<double> voidage_rates(np, 0.0);
rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegionIdx_, surface_rates, voidage_rates);
for (int p = 0; p < np; ++p) {
well_state.wellReservoirRates(index_of_well_)[p] = voidage_rates[p];
}
}
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkIndividualConstraints(WellState& well_state,
const SummaryState& summaryState) const
{
const auto& well = well_ecl_;
const PhaseUsage& pu = phaseUsage();
const int well_index = index_of_well_;
if (well.isInjector()) {
const auto controls = well.injectionControls(summaryState);
auto currentControl = well_state.currentInjectionControl(well_index);
if (controls.hasControl(Well::InjectorCMode::BHP) && currentControl != Well::InjectorCMode::BHP)
{
const auto& bhp = controls.bhp_limit;
double current_bhp = well_state.bhp(well_index);
if (bhp < current_bhp) {
well_state.currentInjectionControl(well_index, Well::InjectorCMode::BHP);
return true;
}
}
if (controls.hasControl(Well::InjectorCMode::RATE) && currentControl != Well::InjectorCMode::RATE)
{
InjectorType injectorType = controls.injector_type;
double current_rate = 0.0;
switch (injectorType) {
case InjectorType::WATER:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
break;
}
case InjectorType::OIL:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
break;
}
case InjectorType::GAS:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
break;
}
default:
throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
}
if (controls.surface_rate < current_rate) {
well_state.currentInjectionControl(well_index, Well::InjectorCMode::RATE);
return true;
}
}
if (controls.hasControl(Well::InjectorCMode::RESV) && currentControl != Well::InjectorCMode::RESV)
{
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.reservoir_rate < current_rate) {
currentControl = Well::InjectorCMode::RESV;
return true;
}
}
if (controls.hasControl(Well::InjectorCMode::THP) && currentControl != Well::InjectorCMode::THP)
{
const auto& thp = getTHPConstraint(summaryState);
double current_thp = well_state.thp(well_index);
if (thp < current_thp) {
currentControl = Well::InjectorCMode::THP;
return true;
}
}
}
if (well.isProducer( )) {
const auto controls = well.productionControls(summaryState);
auto currentControl = well_state.currentProductionControl(well_index);
if (controls.hasControl(Well::ProducerCMode::BHP) && currentControl != Well::ProducerCMode::BHP )
{
const double bhp = controls.bhp_limit;
double current_bhp = well_state.bhp(well_index);
if (bhp > current_bhp) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::BHP);
return true;
}
}
if (controls.hasControl(Well::ProducerCMode::ORAT) && currentControl != Well::ProducerCMode::ORAT) {
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
if (controls.oil_rate < current_rate ) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::ORAT);
return true;
}
}
if (controls.hasControl(Well::ProducerCMode::WRAT) && currentControl != Well::ProducerCMode::WRAT ) {
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
if (controls.water_rate < current_rate ) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::WRAT);
return true;
}
}
if (controls.hasControl(Well::ProducerCMode::GRAT) && currentControl != Well::ProducerCMode::GRAT ) {
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.gas_rate < current_rate ) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::GRAT);
return true;
}
}
if (controls.hasControl(Well::ProducerCMode::LRAT) && currentControl != Well::ProducerCMode::LRAT) {
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
current_rate -= well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
if (controls.liquid_rate < current_rate ) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::LRAT);
return true;
}
}
if (controls.hasControl(Well::ProducerCMode::RESV) && currentControl != Well::ProducerCMode::RESV ) {
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
current_rate -= well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
current_rate -= well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
current_rate -= well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.prediction_mode && controls.resv_rate < current_rate) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::RESV);
return true;
}
if (!controls.prediction_mode) {
const int fipreg = 0; // not considering the region for now
const int np = number_of_phases_;
std::vector<double> surface_rates(np, 0.0);
if( pu.phase_used[BlackoilPhases::Aqua] )
surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = controls.water_rate;
if( pu.phase_used[BlackoilPhases::Liquid] )
surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.oil_rate;
if( pu.phase_used[BlackoilPhases::Vapour] )
surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.gas_rate;
std::vector<double> voidage_rates(np, 0.0);
rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegionIdx_, surface_rates, voidage_rates);
double resv_rate = 0.0;
for (int p = 0; p < np; ++p) {
resv_rate += voidage_rates[p];
}
if (resv_rate < current_rate) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::RESV);
return true;
}
}
}
if (controls.hasControl(Well::ProducerCMode::THP) && currentControl != Well::ProducerCMode::THP)
{
const auto& thp = getTHPConstraint(summaryState);
double current_thp = well_state.thp(well_index);
if (thp > current_thp) {
well_state.currentProductionControl(well_index, Well::ProducerCMode::THP);
return true;
}
}
}
return false;
}
template <typename FluidSystem>
std::pair<bool, double>
WellInterfaceFluidSystem<FluidSystem>::
checkGroupConstraintsInj(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
// Translate injector type from control to Phase.
const auto& well_controls = this->well_ecl_.injectionControls(summaryState);
auto injectorType = well_controls.injector_type;
Phase injectionPhase;
switch (injectorType) {
case InjectorType::WATER:
{
injectionPhase = Phase::WATER;
break;
}
case InjectorType::OIL:
{
injectionPhase = Phase::OIL;
break;
}
case InjectorType::GAS:
{
injectionPhase = Phase::GAS;
break;
}
default:
throw("Expected WATER, OIL or GAS as type for injector " + name());
}
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(phaseUsage().num_phases, 1.0);
rateConverter_.calcCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
// Call check for the well's injection phase.
return WellGroupHelpers::checkGroupConstraintsInj(name(),
well_ecl_.groupName(),
group,
well_state,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
injectionPhase,
phaseUsage(),
efficiencyFactor,
schedule,
summaryState,
resv_coeff,
deferred_logger);
}
template <typename FluidSystem>
std::pair<bool, double>
WellInterfaceFluidSystem<FluidSystem>::
checkGroupConstraintsProd(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
// Make conversion factors for RESV <-> surface rates.
std::vector<double> resv_coeff(this->phaseUsage().num_phases, 1.0);
rateConverter_.calcCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
return WellGroupHelpers::checkGroupConstraintsProd(name(),
well_ecl_.groupName(),
group,
well_state,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
phaseUsage(),
efficiencyFactor,
schedule,
summaryState,
resv_coeff,
deferred_logger);
}
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkGroupConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
const auto& well = well_ecl_;
const int well_index = index_of_well_;
if (well.isInjector()) {
auto currentControl = well_state.currentInjectionControl(well_index);
if (currentControl != Well::InjectorCMode::GRUP) {
// This checks only the first encountered group limit,
// in theory there could be several, and then we should
// test all but the one currently applied. At that point,
// this if-statement should be removed and we should always
// check, skipping over only the single group parent whose
// control is the active one for the well (if any).
const auto& group = schedule.getGroup( well.groupName(), current_step_ );
const double efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, double> group_constraint =
checkGroupConstraintsInj(group, well_state, group_state, efficiencyFactor,
schedule, summaryState, deferred_logger);
// If a group constraint was broken, we set the current well control to
// be GRUP.
if (group_constraint.first) {
well_state.currentInjectionControl(index_of_well_, Well::InjectorCMode::GRUP);
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
}
}
return group_constraint.first;
}
}
if (well.isProducer( )) {
auto currentControl = well_state.currentProductionControl(well_index);
if (currentControl != Well::ProducerCMode::GRUP) {
// This checks only the first encountered group limit,
// in theory there could be several, and then we should
// test all but the one currently applied. At that point,
// this if-statement should be removed and we should always
// check, skipping over only the single group parent whose
// control is the active one for the well (if any).
const auto& group = schedule.getGroup( well.groupName(), current_step_ );
const double efficiencyFactor = well.getEfficiencyFactor();
const std::pair<bool, double> group_constraint =
checkGroupConstraintsProd(group, well_state, group_state, efficiencyFactor,
schedule, summaryState, deferred_logger);
// If a group constraint was broken, we set the current well control to
// be GRUP.
if (group_constraint.first) {
well_state.currentProductionControl(index_of_well_, Well::ProducerCMode::GRUP);
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
}
}
return group_constraint.first;
}
}
return false;
}
template <typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const
{
const bool ind_broken = checkIndividualConstraints(well_state, summaryState);
if (ind_broken) {
return true;
} else {
return checkGroupConstraints(well_state, group_state, schedule, summaryState, deferred_logger);
}
}
template<typename FluidSystem>
bool
WellInterfaceFluidSystem<FluidSystem>::
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const double* rates_or_potentials,
DeferredLogger& deferred_logger) const
{
const PhaseUsage& pu = phaseUsage();
if (econ_production_limits.onMinOilRate()) {
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
const double oil_rate = rates_or_potentials[pu.phase_pos[ Oil ] ];
const double min_oil_rate = econ_production_limits.minOilRate();
if (std::abs(oil_rate) < min_oil_rate) {
return true;
}
}
if (econ_production_limits.onMinGasRate() ) {
assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
const double gas_rate = rates_or_potentials[pu.phase_pos[ Gas ] ];
const double min_gas_rate = econ_production_limits.minGasRate();
if (std::abs(gas_rate) < min_gas_rate) {
return true;
}
}
if (econ_production_limits.onMinLiquidRate() ) {
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
const double oil_rate = rates_or_potentials[pu.phase_pos[ Oil ] ];
const double water_rate = rates_or_potentials[pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate;
const double min_liquid_rate = econ_production_limits.minLiquidRate();
if (std::abs(liquid_rate) < min_liquid_rate) {
return true;
}
}
if (econ_production_limits.onMinReservoirFluidRate()) {
deferred_logger.warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
}
return false;
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const
{
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
// function to calculate water cut based on rates
auto waterCut = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double oil_rate = rates[pu.phase_pos[Oil]];
const double water_rate = rates[pu.phase_pos[Water]];
// both rate should be in the same direction
assert(oil_rate * water_rate >= 0.);
const double liquid_rate = oil_rate + water_rate;
if (liquid_rate != 0.) {
return (water_rate / liquid_rate);
} else {
return 0.;
}
};
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
assert(max_water_cut_limit > 0.);
const bool watercut_limit_violated = checkMaxRatioLimitWell(well_state, max_water_cut_limit, waterCut);
if (watercut_limit_violated) {
report.ratio_limit_violated = true;
checkMaxRatioLimitCompletions(well_state, max_water_cut_limit, waterCut, report);
}
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const
{
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
// function to calculate gor based on rates
auto gor = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double oil_rate = rates[pu.phase_pos[Oil]];
const double gas_rate = rates[pu.phase_pos[Gas]];
// both rate should be in the same direction
assert(oil_rate * gas_rate >= 0.);
double gas_oil_ratio = 0.;
if (oil_rate != 0.) {
gas_oil_ratio = gas_rate / oil_rate;
} else {
if (gas_rate != 0.) {
gas_oil_ratio = 1.e100; // big value to mark it as violated
} else {
gas_oil_ratio = 0.0;
}
}
return gas_oil_ratio;
};
const double max_gor_limit = econ_production_limits.maxGasOilRatio();
assert(max_gor_limit > 0.);
const bool gor_limit_violated = checkMaxRatioLimitWell(well_state, max_gor_limit, gor);
if (gor_limit_violated) {
report.ratio_limit_violated = true;
checkMaxRatioLimitCompletions(well_state, max_gor_limit, gor, report);
}
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const
{
assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
// function to calculate wgr based on rates
auto wgr = [](const std::vector<double>& rates,
const PhaseUsage& pu) {
const double water_rate = rates[pu.phase_pos[Water]];
const double gas_rate = rates[pu.phase_pos[Gas]];
// both rate should be in the same direction
assert(water_rate * gas_rate >= 0.);
double water_gas_ratio = 0.;
if (gas_rate != 0.) {
water_gas_ratio = water_rate / gas_rate;
} else {
if (water_rate != 0.) {
water_gas_ratio = 1.e100; // big value to mark it as violated
} else {
water_gas_ratio = 0.0;
}
}
return water_gas_ratio;
};
const double max_wgr_limit = econ_production_limits.maxWaterGasRatio();
assert(max_wgr_limit > 0.);
const bool wgr_limit_violated = checkMaxRatioLimitWell(well_state, max_wgr_limit, wgr);
if (wgr_limit_violated) {
report.ratio_limit_violated = true;
checkMaxRatioLimitCompletions(well_state, max_wgr_limit, wgr, report);
}
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report,
DeferredLogger& deferred_logger) const
{
// TODO: not sure how to define the worst-offending completion when more than one
// ratio related limit is violated.
// The defintion used here is that we define the violation extent based on the
// ratio between the value and the corresponding limit.
// For each violated limit, we decide the worst-offending completion separately.
// Among the worst-offending completions, we use the one has the biggest violation
// extent.
if (econ_production_limits.onMaxWaterCut()) {
checkMaxWaterCutLimit(econ_production_limits, well_state, report);
}
if (econ_production_limits.onMaxGasOilRatio()) {
checkMaxGORLimit(econ_production_limits, well_state, report);
}
if (econ_production_limits.onMaxWaterGasRatio()) {
checkMaxWGRLimit(econ_production_limits, well_state, report);
}
if (econ_production_limits.onMaxGasLiquidRatio()) {
deferred_logger.warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
}
if (report.ratio_limit_violated) {
assert(report.worst_offending_completion != INVALIDCOMPLETION);
assert(report.violation_extent > 1.);
}
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
updateWellTestStateEconomic(const WellState& well_state,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const
{
if (this->wellIsStopped())
return;
const WellEconProductionLimits& econ_production_limits = well_ecl_.getEconLimits();
// if no limit is effective here, then continue to the next well
if ( !econ_production_limits.onAnyEffectiveLimit() ) {
return;
}
// flag to check if the mim oil/gas rate limit is violated
bool rate_limit_violated = false;
const auto& quantity_limit = econ_production_limits.quantityLimit();
const int np = number_of_phases_;
if (econ_production_limits.onAnyRateLimit()) {
if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN)
rate_limit_violated = checkRateEconLimits(econ_production_limits, &well_state.wellPotentials()[index_of_well_ * np], deferred_logger);
else {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger);
}
}
if (rate_limit_violated) {
if (econ_production_limits.endRun()) {
const std::string warning_message = std::string("ending run after well closed due to economic limits")
+ std::string("is not supported yet \n")
+ std::string("the program will keep running after ") + name()
+ std::string(" is closed");
deferred_logger.warning("NOT_SUPPORTING_ENDRUN", warning_message);
}
if (econ_production_limits.validFollowonWell()) {
deferred_logger.warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
}
well_test_state.closeWell(name(), WellTestConfig::Reason::ECONOMIC, simulation_time);
if (write_message_to_opmlog) {
if (this->well_ecl_.getAutomaticShutIn()) {
const std::string msg = std::string("well ") + name() + std::string(" will be shut due to rate economic limit");
deferred_logger.info(msg);
} else {
const std::string msg = std::string("well ") + name() + std::string(" will be stopped due to rate economic limit");
deferred_logger.info(msg);
}
}
// the well is closed, not need to check other limits
return;
}
if ( !econ_production_limits.onAnyRatioLimit() ) {
// there is no need to check the ratio limits
return;
}
// checking for ratio related limits, mostly all kinds of ratio.
RatioLimitCheckReport ratio_report;
checkRatioEconLimits(econ_production_limits, well_state, ratio_report, deferred_logger);
if (ratio_report.ratio_limit_violated) {
const auto workover = econ_production_limits.workover();
switch (workover) {
case WellEconProductionLimits::EconWorkover::CON:
{
const int worst_offending_completion = ratio_report.worst_offending_completion;
well_test_state.addClosedCompletion(name(), worst_offending_completion, simulation_time);
if (write_message_to_opmlog) {
if (worst_offending_completion < 0) {
const std::string msg = std::string("Connection ") + std::to_string(- worst_offending_completion)
+ std::string(" for well ") + name() + std::string(" will be closed due to economic limit");
deferred_logger.info(msg);
} else {
const std::string msg = std::string("Completion ") + std::to_string(worst_offending_completion)
+ std::string(" for well ") + name() + std::string(" will be closed due to economic limit");
deferred_logger.info(msg);
}
}
bool allCompletionsClosed = true;
const auto& connections = well_ecl_.getConnections();
for (const auto& connection : connections) {
if (connection.state() == Connection::State::OPEN
&& !well_test_state.hasCompletion(name(), connection.complnum())) {
allCompletionsClosed = false;
}
}
if (allCompletionsClosed) {
well_test_state.closeWell(name(), WellTestConfig::Reason::ECONOMIC, simulation_time);
if (write_message_to_opmlog) {
if (this->well_ecl_.getAutomaticShutIn()) {
const std::string msg = name() + std::string(" will be shut due to last completion closed");
deferred_logger.info(msg);
} else {
const std::string msg = name() + std::string(" will be stopped due to last completion closed");
deferred_logger.info(msg);
}
}
}
break;
}
case WellEconProductionLimits::EconWorkover::WELL:
{
well_test_state.closeWell(name(), WellTestConfig::Reason::ECONOMIC, simulation_time);
if (write_message_to_opmlog) {
if (well_ecl_.getAutomaticShutIn()) {
// tell the control that the well is closed
const std::string msg = name() + std::string(" will be shut due to ratio economic limit");
deferred_logger.info(msg);
} else {
const std::string msg = name() + std::string(" will be stopped due to ratio economic limit");
deferred_logger.info(msg);
}
}
break;
}
case WellEconProductionLimits::EconWorkover::NONE:
break;
default:
{
deferred_logger.warning("NOT_SUPPORTED_WORKOVER_TYPE",
"not supporting workover type " + WellEconProductionLimits::EconWorkover2String(workover) );
}
}
}
}
template<typename FluidSystem>
void
WellInterfaceFluidSystem<FluidSystem>::
updateWellTestState(const WellState& well_state,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const
{
// currently, we only updateWellTestState for producers
if (this->isInjector()) {
return;
}
// Based on current understanding, only under prediction mode, we need to shut well due to various
// reasons or limits. With more knowlage or testing cases later, this might need to be corrected.
if (!underPredictionMode() ) {
return;
}
// updating well test state based on physical (THP/BHP) limits.
updateWellTestStatePhysical(well_state, simulationTime, writeMessageToOPMLog, wellTestState, deferred_logger);
// updating well test state based on Economic limits.
updateWellTestStateEconomic(well_state, simulationTime, writeMessageToOPMLog, wellTestState, deferred_logger);
// TODO: well can be shut/closed due to other reasons
}
template<typename FluidSystem>
template <typename RatioFunc>
void WellInterfaceFluidSystem<FluidSystem>::
checkMaxRatioLimitCompletions(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const
{
int worst_offending_completion = INVALIDCOMPLETION;
// the maximum water cut value of the completions
// it is used to identify the most offending completion
double max_ratio_completion = 0;
const int np = number_of_phases_;
const auto * perf_phase_rates = &well_state.perfPhaseRates()[first_perf_ * np];
// look for the worst_offending_completion
for (const auto& completion : completions_) {
std::vector<double> completion_rates(np, 0.0);
// looping through the connections associated with the completion
const std::vector<int>& conns = completion.second;
for (const int c : conns) {
for (int p = 0; p < np; ++p) {
const double connection_rate = perf_phase_rates[c * np + p];
completion_rates[p] += connection_rate;
}
} // end of for (const int c : conns)
parallel_well_info_.communication().sum(completion_rates.data(), completion_rates.size());
const double ratio_completion = ratioFunc(completion_rates, phaseUsage());
if (ratio_completion > max_ratio_completion) {
worst_offending_completion = completion.first;
max_ratio_completion = ratio_completion;
}
} // end of for (const auto& completion : completions_)
assert(max_ratio_completion > max_ratio_limit);
assert(worst_offending_completion != INVALIDCOMPLETION);
const double violation_extent = max_ratio_completion / max_ratio_limit;
assert(violation_extent > 1.0);
if (violation_extent > report.violation_extent) {
report.worst_offending_completion = worst_offending_completion;
report.violation_extent = violation_extent;
}
}
template<typename FluidSystem>
template<typename RatioFunc>
bool WellInterfaceFluidSystem<FluidSystem>::
checkMaxRatioLimitWell(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc) const
{
const int np = number_of_phases_;
std::vector<double> well_rates(np, 0.0);
for (int p = 0; p < np; ++p) {
well_rates[p] = well_state.wellRates(index_of_well_)[p];
}
const double well_ratio = ratioFunc(well_rates, phaseUsage());
return (well_ratio > max_ratio_limit);
}
template class WellInterfaceFluidSystem<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;
template class WellInterfaceFluidSystem<BlackOilFluidSystem<double,EclAlternativeBlackOilIndexTraits>>;
} // namespace Opm

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

@@ -0,0 +1,160 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
Copyright 2019 Norce
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_WELLINTERFACE_FLUID_SYSTEM_HEADER_INCLUDED
#define OPM_WELLINTERFACE_FLUID_SYSTEM_HEADER_INCLUDED
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <limits>
namespace Opm
{
namespace RateConverter
{
template <class FluidSystem, class Region> class SurfaceToReservoirVoidage;
}
class Group;
class GroupState;
class Schedule;
class WellState;
template<class FluidSystem>
class WellInterfaceFluidSystem : public WellInterfaceGeneric {
public:
void updateWellTestState(const WellState& well_state,
const double& simulationTime,
const bool& writeMessageToOPMLog,
WellTestState& wellTestState,
DeferredLogger& deferred_logger) const;
protected:
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<FluidSystem, std::vector<int>>;
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
// to indicate a invalid completion
static constexpr int INVALIDCOMPLETION = std::numeric_limits<int>::max();
WellInterfaceFluidSystem(const Well& well,
const ParallelWellInfo& parallel_well_info,
const int time_step,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
// updating the voidage rates in well_state when requested
void calculateReservoirRates(WellState& well_state) const;
bool checkIndividualConstraints(WellState& well_state,
const SummaryState& summaryState) const;
std::pair<bool, double> checkGroupConstraintsInj(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
std::pair<bool, double> checkGroupConstraintsProd(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const double efficiencyFactor,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
bool checkGroupConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
bool checkConstraints(WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
DeferredLogger& deferred_logger) const;
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const double* rates_or_potentials,
Opm::DeferredLogger& deferred_logger) const;
struct RatioLimitCheckReport{
bool ratio_limit_violated = false;
int worst_offending_completion = INVALIDCOMPLETION;
double violation_extent = 0.0;
};
void checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkMaxGORLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkMaxWGRLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report) const;
void checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state,
RatioLimitCheckReport& report,
DeferredLogger& deferred_logger) const;
void updateWellTestStateEconomic(const WellState& well_state,
const double simulation_time,
const bool write_message_to_opmlog,
WellTestState& well_test_state,
DeferredLogger& deferred_logger) const;
// For the conversion between the surface volume rate and reservoir voidage rate
const RateConverterType& rateConverter_;
private:
template <typename RatioFunc>
void checkMaxRatioLimitCompletions(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc,
RatioLimitCheckReport& report) const;
template<typename RatioFunc>
bool checkMaxRatioLimitWell(const WellState& well_state,
const double max_ratio_limit,
const RatioFunc& ratioFunc) const;
};
}
#endif // OPM_WELLINTERFACE_FLUID_SYSTEM_HEADER_INCLUDED

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

@@ -382,5 +382,4 @@ double WellInterfaceGeneric::getALQ(const WellState& well_state) const
return well_state.getALQ(name());
}
} // namespace Opm

1167
opm/simulators/wells/WellInterface_impl.hpp
View File

File diff suppressed because it is too large Load Diff