opm-simulators/opm/autodiff/WellInterface_impl.hpp

692 lines
19 KiB
C++

/*
Copyright 2017 SINTEF ICT, Applied Mathematics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
namespace Opm
{
template<typename TypeTag>
WellInterface<TypeTag>::
WellInterface(const Well* well, const int time_step, const Wells* wells)
: well_ecl_(well)
, current_step_(time_step)
{
// TODO: trying to use wells struct as little as possible here, be prepared to
// remove the wells struct in future
const std::string& well_name = well->name();
// looking for the location of the well in the wells struct
int index_well;
for (index_well = 0; index_well < wells->number_of_wells; ++index_well) {
if (well_name == std::string(wells->name[index_well])) {
break;
}
}
// should not enter the constructor if the well does not exist in the wells struct
// here, just another assertion.
assert(index_well != wells->number_of_wells);
name_ = well_name;
index_of_well_ = index_well;
well_type_ = wells->type[index_well];
allow_cf_ = wells->allow_cf[index_well];
number_of_phases_ = wells->number_of_phases;
// copying the comp_frac
{
comp_frac_.resize(number_of_phases_);
const int index_begin = index_well * number_of_phases_;
std::copy(wells->comp_frac + index_begin,
wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
}
well_controls_ = wells->ctrls[index_well];
ref_depth_ = wells->depth_ref[index_well];
// perforations related
{
const int perf_index_begin = wells->well_connpos[index_well];
const int perf_index_end = wells->well_connpos[index_well + 1];
number_of_perforations_ = perf_index_end - perf_index_begin;
first_perf_ = perf_index_begin;
well_cell_.resize(number_of_perforations_);
std::copy(wells->well_cells + perf_index_begin,
wells->well_cells + perf_index_end,
well_cell_.begin() );
well_index_.resize(number_of_perforations_);
std::copy(wells->WI + perf_index_begin,
wells->WI + perf_index_end,
well_index_.begin() );
saturation_table_number_.resize(number_of_perforations_);
std::copy(wells->sat_table_id + perf_index_begin,
wells->sat_table_id + perf_index_end,
saturation_table_number_.begin() );
// TODO: not sure about the processing of depth for perforations here
// Will revisit here later. There are different ways and the definition for different wells
// can be different, it is possible that we need to remove this from the WellInterface
perf_depth_.resize(number_of_perforations_, 0.);
const auto& completion_set = well->getCompletions(time_step);
for (int i = 0; i < number_of_perforations_; ++i) {
perf_depth_[i] = completion_set.get(i).getCenterDepth();
}
}
well_efficiency_factor_ = 1.0;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
init(const PhaseUsage* phase_usage_arg,
const std::vector<bool>* active_arg,
const VFPProperties* vfp_properties_arg,
const double gravity_arg,
const int /* num_cells */)
{
phase_usage_ = phase_usage_arg;
active_ = active_arg;
vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg;
}
template<typename TypeTag>
const std::string&
WellInterface<TypeTag>::
name() const
{
return name_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
indexOfWell() const
{
return index_of_well_;
}
template<typename TypeTag>
WellType
WellInterface<TypeTag>::
wellType() const
{
return well_type_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
numberOfPhases() const
{
return number_of_phases_;
}
template<typename TypeTag>
const std::vector<double>&
WellInterface<TypeTag>::
compFrac() const
{
return comp_frac_;
}
template<typename TypeTag>
WellControls*
WellInterface<TypeTag>::
wellControls() const
{
return well_controls_;
}
template<typename TypeTag>
const std::vector<int>&
WellInterface<TypeTag>::
saturationTableNumber() const
{
return saturation_table_number_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
numberOfPerforations() const
{
return number_of_perforations_;
}
template<typename TypeTag>
const std::vector<double>&
WellInterface<TypeTag>::
wellIndex() const
{
return well_index_;
}
template<typename TypeTag>
const std::vector<double>&
WellInterface<TypeTag>::
perfDepth() const
{
return perf_depth_;
}
template<typename TypeTag>
const std::vector<int>&
WellInterface<TypeTag>::
wellCells() const
{
return well_cell_;
}
template<typename TypeTag>
const std::vector<bool>&
WellInterface<TypeTag>::
active() const
{
assert(active_);
return *active_;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
allowCrossFlow() const
{
return allow_cf_;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
setWellEfficiencyFactor(const double efficiency_factor)
{
well_efficiency_factor_ = efficiency_factor;
}
template<typename TypeTag>
const PhaseUsage&
WellInterface<TypeTag>::
phaseUsage() const
{
assert(phase_usage_);
return *phase_usage_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowPhaseToEbosCompIdx( const int phaseIdx ) const
{
const int phaseToComp[ 3 ] = { FluidSystem::waterCompIdx, FluidSystem::oilCompIdx, FluidSystem::gasCompIdx};
if (phaseIdx > 2 )
return phaseIdx;
return phaseToComp[ phaseIdx ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowToEbosPvIdx( const int flowPv ) const
{
const int flowToEbos[ 3 ] = {
BlackoilIndices::pressureSwitchIdx,
BlackoilIndices::waterSaturationIdx,
BlackoilIndices::compositionSwitchIdx
};
if (flowPv > 2 )
return flowPv;
return flowToEbos[ flowPv ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
{
assert(phaseIdx < 3);
const int flowToEbos[ 3 ] = { FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx };
return flowToEbos[ phaseIdx ];
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
numPhases() const
{
return number_of_phases_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
numComponents() const
{
if (numPhases() == 2) {
return 2;
}
int numComp = FluidSystem::numComponents;
if (has_solvent) {
numComp ++;
}
return numComp;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
wsolvent() const
{
if (!has_solvent) {
return 0.0;
}
WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
if (injection.injectorType == WellInjector::GAS) {
double solvent_fraction = well_ecl_->getSolventFraction(current_step_);
return solvent_fraction;
}
assert(false);
return 0.0;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
wpolymer() const
{
if (!has_polymer) {
return 0.0;
}
WellInjectionProperties injection = well_ecl_->getInjectionProperties(current_step_);
WellPolymerProperties polymer = well_ecl_->getPolymerProperties(current_step_);
if (injection.injectorType == WellInjector::WATER) {
const double polymer_injection_concentration = polymer.m_polymerConcentration;
return polymer_injection_concentration;
}
assert(false); // TODO: find a more logical way to handle this situation
return 0.0;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
mostStrictBhpFromBhpLimits() const
{
double bhp;
// initial bhp value, making the value not usable
switch( well_type_ ) {
case INJECTOR:
bhp = std::numeric_limits<double>::max();
break;
case PRODUCER:
bhp = -std::numeric_limits<double>::max();
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
}
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
// finding a BHP constraint
if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
// get the bhp constraint value, it should always be postive assummingly
const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
switch(well_type_) {
case INJECTOR: // using the lower bhp contraint from Injectors
if (bhp_target < bhp) {
bhp = bhp_target;
}
break;
case PRODUCER:
if (bhp_target > bhp) {
bhp = bhp_target;
}
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name());
} // end of switch
}
}
return bhp;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
wellHasTHPConstraints() const
{
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_controls_, ctrl_index) == THP) {
return true;
}
}
return false;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
const Opm::PhaseUsage& pu = *phase_usage_;
const int np = numberOfPhases();
if (econ_production_limits.onMinOilRate()) {
assert(active()[Oil]);
const double oil_rate = well_state.wellRates()[index_of_well_ * np + 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(active()[Gas]);
const double gas_rate = well_state.wellRates()[index_of_well_ * np + 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(active()[Oil]);
assert(active()[Water]);
const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[index_of_well_ * np + 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()) {
OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
}
return false;
}
template<typename TypeTag>
typename WellInterface<TypeTag>::RatioCheckTuple
WellInterface<TypeTag>::
checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
bool water_cut_limit_violated = false;
int worst_offending_connection = INVALIDCONNECTION;
bool last_connection = false;
double violation_extent = -1.0;
const int np = numberOfPhases();
const Opm::PhaseUsage& pu = *phase_usage_;
const int well_number = index_of_well_;
assert(active()[Oil]);
assert(active()[Water]);
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate;
double water_cut;
if (std::abs(liquid_rate) != 0.) {
water_cut = water_rate / liquid_rate;
} else {
water_cut = 0.0;
}
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
if (water_cut > max_water_cut_limit) {
water_cut_limit_violated = true;
}
if (water_cut_limit_violated) {
// need to handle the worst_offending_connection
const int perf_start = first_perf_;
const int perf_number = number_of_perforations_;
std::vector<double> water_cut_perf(perf_number);
for (int perf = 0; perf < perf_number; ++perf) {
const int i_perf = perf_start + perf;
const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
if (std::abs(liquid_perf_rate) != 0.) {
water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
} else {
water_cut_perf[perf] = 0.;
}
}
last_connection = (perf_number == 1);
if (last_connection) {
worst_offending_connection = 0;
violation_extent = water_cut_perf[0] / max_water_cut_limit;
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
double max_water_cut_perf = 0.;
for (int perf = 0; perf < perf_number; ++perf) {
if (water_cut_perf[perf] > max_water_cut_perf) {
worst_offending_connection = perf;
max_water_cut_perf = water_cut_perf[perf];
}
}
assert(max_water_cut_perf != 0.);
assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
violation_extent = max_water_cut_perf / max_water_cut_limit;
}
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
template<typename TypeTag>
typename WellInterface<TypeTag>::RatioCheckTuple
WellInterface<TypeTag>::
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
const WellState& well_state) const
{
// TODO: not sure how to define the worst-offending connection 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 connection separately.
// Among the worst-offending connections, we use the one has the biggest violation
// extent.
bool any_limit_violated = false;
bool last_connection = false;
int worst_offending_connection = INVALIDCONNECTION;
double violation_extent = -1.0;
if (econ_production_limits.onMaxWaterCut()) {
const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state);
bool water_cut_violated = std::get<0>(water_cut_return);
if (water_cut_violated) {
any_limit_violated = true;
const double violation_extent_water_cut = std::get<3>(water_cut_return);
if (violation_extent_water_cut > violation_extent) {
violation_extent = violation_extent_water_cut;
worst_offending_connection = std::get<2>(water_cut_return);
last_connection = std::get<1>(water_cut_return);
}
}
}
if (econ_production_limits.onMaxGasOilRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
}
if (econ_production_limits.onMaxWaterGasRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
}
if (econ_production_limits.onMaxGasLiquidRatio()) {
OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
}
if (any_limit_violated) {
assert(worst_offending_connection >=0);
assert(violation_extent > 1.);
}
return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
}
}