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 #1802 from totto82/addWellPotMSW

Browse Source 
Add well potential output to MSW
This commit is contained in:
Atgeirr Flø Rasmussen
2019-05-27 15:31:21 +02:00
committed by GitHub
parent 8ac85b892c 613473082f
commit c554fdb5fb
11 changed files with 381 additions and 108 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
opm/autodiff/BlackoilModelParametersEbos.hpp
View File

@@ -72,7 +72,7 @@ SET_BOOL_PROP(FlowModelParameters, UpdateEquationsScaling, false);
SET_BOOL_PROP(FlowModelParameters, UseUpdateStabilization, true);
SET_BOOL_PROP(FlowModelParameters, MatrixAddWellContributions, false);
SET_SCALAR_PROP(FlowModelParameters, TolerancePressureMsWells, 0.01 *1e5);
SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 2.0 *1e5);
SET_SCALAR_PROP(FlowModelParameters, MaxPressureChangeMsWells, 1e6);
SET_BOOL_PROP(FlowModelParameters, UseInnerIterationsMsWells, true);
SET_INT_PROP(FlowModelParameters, MaxInnerIterMsWells, 100);

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

@@ -281,7 +281,7 @@ namespace Opm {
std::vector<bool> is_cell_perforated_;
// create the well container
std::vector<WellInterfacePtr > createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger);
std::vector<WellInterfacePtr > createWellContainer(const int time_step, const Wells* wells, const bool allow_closing_opening_wells, Opm::DeferredLogger& deferred_logger);
WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step, Opm::DeferredLogger& deferred_logger) const;

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

@@ -306,7 +306,7 @@ namespace Opm {
wellTesting(reportStepIdx, simulationTime, local_deferredLogger);
// create the well container
well_container_ = createWellContainer(reportStepIdx, local_deferredLogger);
well_container_ = createWellContainer(reportStepIdx, wells(), /*allow_closing_opening_wells=*/true, local_deferredLogger);
// do the initialization for all the wells
// TODO: to see whether we can postpone of the intialization of the well containers to
@@ -414,7 +414,6 @@ namespace Opm {
updateWellTestState(simulationTime, wellTestState_);
// calculate the well potentials for output
// TODO: when necessary
try {
std::vector<double> well_potentials;
computeWellPotentials(well_potentials, local_deferredLogger);
@@ -513,7 +512,7 @@ namespace Opm {
template<typename TypeTag>
std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
BlackoilWellModel<TypeTag>::
createWellContainer(const int time_step, Opm::DeferredLogger& deferred_logger)
createWellContainer(const int time_step, const Wells* wells, const bool allow_closing_opening_wells, Opm::DeferredLogger& deferred_logger)
{
std::vector<WellInterfacePtr> well_container;
@@ -525,7 +524,7 @@ namespace Opm {
// With the following way, it will have the same order with wells struct
// Hopefully, it can generate the same residual history with master branch
for (int w = 0; w < nw; ++w) {
const std::string well_name = std::string(wells()->name[w]);
const std::string well_name = std::string(wells->name[w]);
// finding the location of the well in wells_ecl
const int nw_wells_ecl = wells_ecl_.size();
@@ -543,58 +542,61 @@ namespace Opm {
const Well2& well_ecl = wells_ecl_[index_well];
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit
if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) {
// TODO: more checking here, to make sure this standard more specific and complete
// maybe there is some WCON keywords will not open the well
if (well_state_.effectiveEventsOccurred(w)) {
if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) {
// The well was shut this timestep, we are most likely retrying
// a timestep without the well in question, after it caused
// repeated timestep cuts. It should therefore not be opened,
// even if it was new or received new targets this report step.
well_state_.setEffectiveEventsOccurred(w, false);
} else {
wellTestState_.openWell(well_name);
if (allow_closing_opening_wells) {
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit
if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL ) ) {
// TODO: more checking here, to make sure this standard more specific and complete
// maybe there is some WCON keywords will not open the well
if (well_state_.effectiveEventsOccurred(w)) {
if (wellTestState_.lastTestTime(well_name) == ebosSimulator_.time()) {
// The well was shut this timestep, we are most likely retrying
// a timestep without the well in question, after it caused
// repeated timestep cuts. It should therefore not be opened,
// even if it was new or received new targets this report step.
well_state_.setEffectiveEventsOccurred(w, false);
} else {
wellTestState_.openWell(well_name);
}
}
}
}
// TODO: should we do this for all kinds of closing reasons?
// something like wellTestState_.hasWell(well_name)?
if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) ||
wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) {
if( well_ecl.getAutomaticShutIn() ) {
// shut wells are not added to the well container
// TODO: make a function from well_state side to handle the following
well_state_.thp()[w] = 0.;
well_state_.bhp()[w] = 0.;
const int np = numPhases();
for (int p = 0; p < np; ++p) {
well_state_.wellRates()[np * w + p] = 0.;
// TODO: should we do this for all kinds of closing reasons?
// something like wellTestState_.hasWell(well_name)?
if ( wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) ||
wellTestState_.hasWell(well_name, WellTestConfig::Reason::PHYSICAL) ) {
if( well_ecl.getAutomaticShutIn() ) {
// shut wells are not added to the well container
// TODO: make a function from well_state side to handle the following
well_state_.thp()[w] = 0.;
well_state_.bhp()[w] = 0.;
const int np = numPhases();
for (int p = 0; p < np; ++p) {
well_state_.wellRates()[np * w + p] = 0.;
}
continue;
} else {
// close wells are added to the container but marked as closed
struct WellControls* well_controls = wells->ctrls[w];
well_controls_stop_well(well_controls);
}
continue;
} else {
// close wells are added to the container but marked as closed
struct WellControls* well_controls = wells()->ctrls[w];
well_controls_stop_well(well_controls);
}
}
// Use the pvtRegionIdx from the top cell
const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]];
const int well_cell_top = wells->well_cells[wells->well_connpos[w]];
const int pvtreg = pvt_region_idx_[well_cell_top];
if ( !well_ecl.isMultiSegment() || !param_.use_multisegment_well_) {
if ( GET_PROP_VALUE(TypeTag, EnablePolymerMW) && well_ecl.isInjector() ) {
well_container.emplace_back(new StandardWellV<TypeTag>(well_ecl, time_step, wells(),
well_container.emplace_back(new StandardWellV<TypeTag>(well_ecl, time_step, wells,
param_, *rateConverter_, pvtreg, numComponents() ) );
} else {
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells(),
well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells,
param_, *rateConverter_, pvtreg, numComponents() ) );
}
} else {
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells(),
well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells,
param_, *rateConverter_, pvtreg, numComponents() ) );
}
}
@@ -708,7 +710,7 @@ namespace Opm {
// reservoir state, will tihs be a better place to inialize the explict information?
}
assembleWellEq(B_avg, dt, local_deferredLogger);
assembleWellEq(B_avg, dt, local_deferredLogger);
} catch (std::exception& e) {
exception_thrown = 1;
@@ -1066,17 +1068,74 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
computeWellPotentials(std::vector<double>& well_potentials, Opm::DeferredLogger& deferred_logger)
{
Opm::DeferredLogger local_deferredLogger;
// number of wells and phases
const int nw = numWells();
const int np = numPhases();
well_potentials.resize(nw * np, 0.0);
const int reportStepIdx = ebosSimulator_.episodeIndex();
const double invalid_alq = -1e100;
const double invalid_vfp = -2147483647;
auto well_state_copy = well_state_;
const Wells* local_wells = clone_wells(wells());
std::vector<WellInterfacePtr> well_container_copy = createWellContainer(reportStepIdx, local_wells, /*allow_closing_opening_wells=*/ false, deferred_logger);
// average B factors are required for the convergence checking of well equations
// Note: this must be done on all processes, even those with
// no wells needing testing, otherwise we will have locking.
std::vector< Scalar > B_avg(numComponents(), Scalar() );
computeAverageFormationFactor(B_avg);
const Opm::SummaryConfig& summaryConfig = ebosSimulator_.vanguard().summaryConfig();
int exception_thrown = 0;
try {
for (const auto& well : well_container_) {
for (const auto& well : well_container_copy) {
// Only compute the well potential when asked for
well->init(&phase_usage_, depth_, gravity_, number_of_cells_);
WellControls* wc = well->wellControls();
well_controls_clear(wc);
well_controls_assert_number_of_phases( wc , np);
if (well->wellType() == INJECTOR) {
const auto& injectionProperties = well->wellEcl()->getInjectionProperties();
if (injectionProperties.hasInjectionControl(WellInjector::THP)) {
const double thp_limit = injectionProperties.THPLimit;
const int vfp_number = injectionProperties.VFPTableNumber;
well_controls_add_new(THP, thp_limit, invalid_alq, vfp_number, NULL, wc);
}
// we always have a bhp limit
const double bhp_limit = injectionProperties.BHPLimit;
well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
} else {
const auto& productionProperties = well->wellEcl()->getProductionProperties();
if (productionProperties.hasProductionControl(WellProducer::THP)) {
const double thp_limit = productionProperties.THPLimit;
const double alq_value = productionProperties.ALQValue;
const int vfp_number = productionProperties.VFPTableNumber;
well_controls_add_new(THP, thp_limit, alq_value, vfp_number, NULL, wc);
}
// we always have a bhp limit
const double bhp_limit = productionProperties.BHPLimit;
well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
well->setVFPProperties(vfp_properties_.get());
}
if (has_polymer_)
{
const Grid& grid = ebosSimulator_.vanguard().grid();
if (PolymerModule::hasPlyshlog() || GET_PROP_VALUE(TypeTag, EnablePolymerMW) ) {
well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_, deferred_logger);
}
}
const bool needed_for_output = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) ||
summaryConfig.hasSummaryKey( "WOPI:" + well->name()) ||
summaryConfig.hasSummaryKey( "WGPI:" + well->name())) && well->wellType() == INJECTOR) ||
@@ -1086,7 +1145,7 @@ namespace Opm {
if (needed_for_output || wellCollection().requireWellPotentials())
{
std::vector<double> potentials;
well->computeWellPotentials(ebosSimulator_, well_state_, potentials, deferred_logger);
well->computeWellPotentials(ebosSimulator_, B_avg, well_state_copy, potentials, deferred_logger);
// putting the sucessfully calculated potentials to the well_potentials
for (int p = 0; p < np; ++p) {
well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]);

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

@@ -139,6 +139,7 @@ namespace Opm
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) override;
@@ -305,6 +306,8 @@ namespace Opm
const bool& allow_cf,
std::vector<EvalWell>& cq_s,
EvalWell& perf_press,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const;
// convert a Eval from reservoir to contain the derivative related to wells
@@ -327,6 +330,12 @@ namespace Opm
const int perf,
std::vector<EvalWell>& mob) const;
void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger);
void assembleControlEq(Opm::DeferredLogger& deferred_logger) const;
void assemblePressureEq(const int seg) const;
@@ -366,7 +375,7 @@ namespace Opm
WellState& well_state,
Opm::DeferredLogger& deferred_logger);
virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override;

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

@@ -541,24 +541,83 @@ namespace Opm
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellPotentials(const Simulator& /* ebosSimulator */,
const WellState& /* well_state */,
computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger)
{
const std::string msg = std::string("Well potential calculation is not supported for multisegment wells \n")
+ "A well potential of zero is returned for output purposes. \n"
+ "If you need well potential to set the guide rate for group controled wells \n"
+ "you will have to change the " + name() + " well to a standard well \n";
deferred_logger.warning("WELL_POTENTIAL_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg);
const int np = number_of_phases_;
well_potentials.resize(np, 0.0);
updatePrimaryVariables(well_state, deferred_logger);
// initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
// TODO: for computeWellPotentials, no derivative is required actually
initPrimaryVariablesEvaluation();
// get the bhp value based on the bhp constraints
const double bhp = Base::mostStrictBhpFromBhpLimits(deferred_logger);
// does the well have a THP related constraint?
if ( !Base::wellHasTHPConstraints() ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
} else {
const std::string msg = std::string("Well potential calculation is not supported for thp controlled multisegment wells \n")
+ "A well potential of zero is returned for output purposes. \n"
+ "If you need well potential computed from thp to set the guide rate for group controled wells \n"
+ "you will have to change the " + name() + " well to a standard well \n";
deferred_logger.warning("WELL_POTENTIAL_FOR_THP_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg);
return;
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger)
{
WellControls* wc = well_controls_;
const int bhp_index = Base::getControlIndex(BHP);
const double orig_bhp = well_controls_iget_target(wc, bhp_index);
const auto orig_current = well_controls_get_current(wc);
well_controls_iset_target(wc, bhp_index, bhp);
well_controls_set_current(wc, bhp_index);
// store a copy of the well state, we don't want to update the real well state
WellState copy = ebosSimulator.problem().wellModel().wellState();
initPrimaryVariablesEvaluation();
const double dt = ebosSimulator.timeStepSize();
// iterate to get a solution that satisfies the bhp potential.
iterateWellEquations(ebosSimulator, B_avg, dt, copy, deferred_logger);
// compute the potential and store in the flux vector.
const int np = number_of_phases_;
well_flux.resize(np, 0.0);
for(int compIdx = 0; compIdx < num_components_; ++compIdx) {
const EvalWell rate = getSegmentRate(0, compIdx);
well_flux[ebosCompIdxToFlowCompIdx(compIdx)] += rate.value();
}
// reset bhp limit
well_controls_iset_target(wc, bhp_index, orig_bhp);
well_controls_set_current(wc, orig_current);
}
@@ -743,7 +802,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
updateWellState(const BVectorWell& dwells,
WellState& well_state,
Opm::DeferredLogger& /* deferred_logger */,
Opm::DeferredLogger& deferred_logger,
const double relaxation_factor) const
{
const double dFLimit = param_.dwell_fraction_max_;
@@ -773,7 +832,7 @@ namespace Opm
const int sign = dwells[seg][SPres] > 0.? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]), relaxation_factor * max_pressure_change);
// const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
primary_variables_[seg][SPres] = old_primary_variables[seg][SPres] - dx_limited;
primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
}
// update the total rate // TODO: should we have a limitation of the total rate change?
@@ -961,6 +1020,8 @@ namespace Opm
const bool& allow_cf,
std::vector<EvalWell>& cq_s,
EvalWell& perf_press,
double& perf_dis_gas_rate,
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const
{
@@ -1077,6 +1138,29 @@ namespace Opm
cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
}
} // end for injection perforations
// calculating the perforation solution gas rate and solution oil rates
if (well_type_ == PRODUCER) {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// TODO: the formulations here remain to be tested with cases with strong crossflow through production wells
// s means standard condition, r means reservoir condition
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// d = 1.0 - rs * rv
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - rv.value() * rs.value();
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_vap_oil_rate = rv.value() * (cq_s[gasCompIdx].value() - rs.value() * cq_s[oilCompIdx].value()) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = rs.value() * (cq_s[oilCompIdx].value() - rv.value() * cq_s[gasCompIdx].value()) / d;
}
}
}
@@ -1877,6 +1961,9 @@ namespace Opm
duneD_ = 0.0;
resWell_ = 0.0;
well_state.wellVaporizedOilRates()[index_of_well_] = 0.;
well_state.wellDissolvedGasRates()[index_of_well_] = 0.;
// for the black oil cases, there will be four equations,
// the first three of them are the mass balance equations, the last one is the pressure equations.
//
@@ -1944,7 +2031,16 @@ namespace Opm
getMobility(ebosSimulator, perf, mob);
std::vector<EvalWell> cq_s(num_components_, 0.0);
EvalWell perf_press;
computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, deferred_logger);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRatePressure(int_quants, mob, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
// updating the solution gas rate and solution oil rate
if (well_type_ == PRODUCER) {
well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate;
well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate;
}
// store the perf pressure and rates
const int rate_start_offset = (first_perf_ + perf) * number_of_phases_;
@@ -1995,7 +2091,7 @@ namespace Opm
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
wellTestingPhysical(Simulator& /* simulator */, const std::vector<double>& /* B_avg */,
wellTestingPhysical(const Simulator& /* simulator */, const std::vector<double>& /* B_avg */,
const double /* simulation_time */, const int /* report_step */,
WellState& /* well_state */, WellTestState& /* welltest_state */, Opm::DeferredLogger& deferred_logger)
{

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

@@ -168,6 +168,7 @@ namespace Opm
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) /* const */ override;
@@ -327,17 +328,22 @@ namespace Opm
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const;
// TODO: maybe we should provide a light version of computePerfRate, which does not include the
// calculation of the derivatives
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const;
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger);
std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const;
@@ -433,7 +439,7 @@ namespace Opm
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
WellState& well_state, WellTestState& welltest_state,
Opm::DeferredLogger& deferred_logger) override;

20
opm/simulators/wells/StandardWellV.hpp
View File

@@ -168,6 +168,7 @@ namespace Opm
/// computing the well potentials for group control
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) /* const */ override;
@@ -330,17 +331,22 @@ namespace Opm
double& perf_vap_oil_rate,
Opm::DeferredLogger& deferred_logger) const;
// TODO: maybe we should provide a light version of computePerfRate, which does not include the
// calculation of the derivatives
void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const;
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger);
std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const;
@@ -434,7 +440,7 @@ namespace Opm
static double relaxationFactorRate(const std::vector<double>& primary_variables,
const BVectorWell& dwells);
virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) override;

71
opm/simulators/wells/StandardWellV_impl.hpp
View File

@@ -1478,7 +1478,7 @@ namespace Opm
// option 2: stick with the above IPR curve
// we use IPR here
std::vector<double> well_rates_bhp_limit;
computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp_limit), well_rates_bhp_limit, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger);
const double thp_limit = this->getTHPConstraint(deferred_logger);
@@ -1667,7 +1667,7 @@ namespace Opm
initPrimaryVariablesEvaluation();
std::vector<double> rates;
computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp), rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
// TODO: double checke the obtained rates
// this is another places we might obtain negative rates
@@ -1690,7 +1690,7 @@ namespace Opm
calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const
{
const double thp_target = this->getTHPConstraint(deferred_logger);
const double thp_control_index = this->getTHPControlIndex();
const double thp_control_index = this->getControlIndex(THP);
const int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index);
const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
@@ -2263,10 +2263,11 @@ namespace Opm
void
StandardWellV<TypeTag>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const
{
const int np = number_of_phases_;
well_flux.resize(np, 0.0);
@@ -2278,11 +2279,13 @@ namespace Opm
// flux for each perforation
std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.});
getMobility(ebosSimulator, perf, mob, deferred_logger);
//double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
//const double Tw = well_index_[perf] * trans_mult;
std::vector<EvalWell> cq_s(num_components_, {numWellEq_ + numEq, 0.});
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, bhp, perf, allow_cf,
computePerfRate(intQuants, mob, EvalWell(numWellEq_ + numEq, bhp), perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
for(int p = 0; p < np; ++p) {
@@ -2293,15 +2296,60 @@ namespace Opm
template<typename TypeTag>
void
StandardWellV<TypeTag>::
computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger)
{
WellControls* wc = well_controls_;
const int bhp_index = Base::getControlIndex(BHP);
const double orig_bhp = well_controls_iget_target(wc, bhp_index);
const auto orig_current = well_controls_get_current(wc);
well_controls_iset_target(wc, bhp_index, bhp);
well_controls_set_current(wc, bhp_index);
// iterate to get a more accurate well density
// create a copy of the well_state to use. If the operability checking is sucessful, we use this one
// to replace the original one
WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
well_state_copy.currentControls()[index_of_well_] = bhp_index;
bool converged = this->solveWellEqUntilConverged(ebosSimulator, B_avg, well_state_copy, deferred_logger);
if (!converged) {
const std::string msg = " well " + name() + " did not get converged during well potential calculations "
"returning zero values for the potential";
deferred_logger.debug(msg);
return;
}
updatePrimaryVariables(well_state_copy, deferred_logger);
computeWellConnectionPressures(ebosSimulator, well_state_copy);
initPrimaryVariablesEvaluation();
computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
// reset bhp limit
well_controls_iset_target(wc, bhp_index, orig_bhp);
well_controls_set_current(wc, orig_current);
}
template<typename TypeTag>
std::vector<double>
StandardWellV<TypeTag>::
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
// TODO: pay attention to the situation that finally the potential is calculated based on the bhp control
// TODO: should we consider the bhp constraints during the iterative process?
@@ -2365,7 +2413,7 @@ namespace Opm
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, potentials, deferred_logger);
// checking whether the potentials have valid values
for (const double value : potentials) {
@@ -2403,6 +2451,7 @@ namespace Opm
void
StandardWellV<TypeTag>::
computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) // const
@@ -2424,7 +2473,7 @@ namespace Opm
if ( !wellHasTHPConstraints() ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
} else {
// the well has a THP related constraint
// checking whether a well is newly added, it only happens at the beginning of the report step
@@ -2436,7 +2485,7 @@ namespace Opm
}
} else {
// We need to generate a reasonable rates to start the iteration process
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
for (double& value : well_potentials) {
// make the value a little safer in case the BHP limits are default ones
// TODO: a better way should be a better rescaling based on the investigation of the VFP table.
@@ -2445,7 +2494,7 @@ namespace Opm
}
}
well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials, deferred_logger);
well_potentials = computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
}
}
@@ -2861,7 +2910,7 @@ namespace Opm
template<typename TypeTag>
void
StandardWellV<TypeTag>::
wellTestingPhysical(Simulator& ebos_simulator, const std::vector<double>& B_avg,
wellTestingPhysical(const Simulator& ebos_simulator, const std::vector<double>& B_avg,
const double /* simulation_time */ , const int /* report_step */,
WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger)
{

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

@@ -1617,7 +1617,7 @@ namespace Opm
calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const
{
const double thp_target = this->getTHPConstraint(deferred_logger);
const double thp_control_index = this->getTHPControlIndex();
const double thp_control_index = this->getControlIndex(THP);
const int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index);
const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
@@ -2152,15 +2152,15 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const
{
const int np = number_of_phases_;
well_flux.resize(np, 0.0);
@@ -2178,7 +2178,7 @@ namespace Opm
std::vector<EvalWell> cq_s(num_components_, 0.0);
double perf_dis_gas_rate = 0.;
double perf_vap_oil_rate = 0.;
computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
computePerfRate(intQuants, mob, EvalWell(bhp), Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
for(int p = 0; p < np; ++p) {
@@ -2189,15 +2189,61 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger)
{
WellControls* wc = well_controls_;
const int bhp_index = Base::getControlIndex(BHP);
const double orig_bhp = well_controls_iget_target(wc, bhp_index);
const auto orig_current = well_controls_get_current(wc);
well_controls_iset_target(wc, bhp_index, bhp);
well_controls_set_current(wc, bhp_index);
// iterate to get a more accurate well density
// create a copy of the well_state to use. If the operability checking is sucessful, we use this one
// to replace the original one
WellState well_state_copy = ebosSimulator.problem().wellModel().wellState();
well_state_copy.currentControls()[index_of_well_] = bhp_index;
bool converged = this->solveWellEqUntilConverged(ebosSimulator, B_avg, well_state_copy, deferred_logger);
if (!converged) {
const std::string msg = " well " + name() + " did not get converged during well potential calculations "
"returning zero values for the potential";
deferred_logger.debug(msg);
return;
}
updatePrimaryVariables(well_state_copy, deferred_logger);
computeWellConnectionPressures(ebosSimulator, well_state_copy);
initPrimaryVariablesEvaluation();
computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
// reset bhp limit
well_controls_iset_target(wc, bhp_index, orig_bhp);
well_controls_set_current(wc, orig_current);
}
template<typename TypeTag>
std::vector<double>
StandardWell<TypeTag>::
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
// TODO: pay attention to the situation that finally the potential is calculated based on the bhp control
// TODO: should we consider the bhp constraints during the iterative process?
@@ -2261,7 +2307,7 @@ namespace Opm
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhp(ebosSimulator, bhp, potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, potentials, deferred_logger);
// checking whether the potentials have valid values
for (const double value : potentials) {
@@ -2303,6 +2349,7 @@ namespace Opm
void
StandardWell<TypeTag>::
computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) // const
@@ -2324,7 +2371,7 @@ namespace Opm
// does the well have a THP related constraint?
if ( !wellHasTHPConstraints() ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
} else {
// the well has a THP related constraint
// checking whether a well is newly added, it only happens at the beginning of the report step
@@ -2336,7 +2383,7 @@ namespace Opm
}
} else {
// We need to generate a reasonable rates to start the iteration process
computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials, deferred_logger);
computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
for (double& value : well_potentials) {
// make the value a little safer in case the BHP limits are default ones
// TODO: a better way should be a better rescaling based on the investigation of the VFP table.
@@ -2345,7 +2392,7 @@ namespace Opm
}
}
well_potentials = computeWellPotentialWithTHP(ebosSimulator, bhp, well_potentials, deferred_logger);
well_potentials = computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
}
}
@@ -2754,7 +2801,7 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
wellTestingPhysical(Simulator& ebos_simulator, const std::vector<double>& B_avg,
wellTestingPhysical(const Simulator& ebos_simulator, const std::vector<double>& B_avg,
const double /* simulation_time */, const int /* report_step */,
WellState& well_state, WellTestState& welltest_state,
Opm::DeferredLogger& deferred_logger)

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

@@ -182,6 +182,7 @@ namespace Opm
// TODO: before we decide to put more information under mutable, this function is not const
virtual void computeWellPotentials(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
std::vector<double>& well_potentials,
Opm::DeferredLogger& deferred_logger) = 0;
@@ -235,7 +236,7 @@ namespace Opm
// TODO: theoretically, it should be a const function
// Simulator is not const is because that assembleWellEq is non-const Simulator
void wellTesting(Simulator& simulator, const std::vector<double>& Bavg,
void wellTesting(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
const WellTestConfig::Reason testing_reason,
/* const */ WellState& well_state, WellTestState& welltest_state,
@@ -354,7 +355,7 @@ namespace Opm
double getTHPConstraint(Opm::DeferredLogger& deferred_logger) const;
int getTHPControlIndex() const;
int getControlIndex(const WellControlType& type) const;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
@@ -385,11 +386,11 @@ namespace Opm
OperabilityStatus operability_status_;
void wellTestingEconomic(Simulator& simulator, const std::vector<double>& B_avg,
void wellTestingEconomic(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
const WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger);
virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
virtual void wellTestingPhysical(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger) = 0;
@@ -405,11 +406,11 @@ namespace Opm
WellTestState& well_test_state,
Opm::DeferredLogger& deferred_logger) const;
void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state,
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state,
const std::vector<double>& B_avg,
Opm::DeferredLogger& deferred_logger);
bool solveWellEqUntilConverged(Simulator& ebosSimulator,
bool solveWellEqUntilConverged(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger);

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

@@ -380,7 +380,7 @@ namespace Opm
WellInterface<TypeTag>::
wellHasTHPConstraints() const
{
return getTHPControlIndex() >= 0;
return getControlIndex(THP) >= 0;
}
@@ -391,7 +391,7 @@ namespace Opm
WellInterface<TypeTag>::
getTHPConstraint(Opm::DeferredLogger& deferred_logger) const
{
const int thp_control_index = getTHPControlIndex();
const int thp_control_index = getControlIndex(THP);
if (thp_control_index < 0) {
OPM_DEFLOG_THROW(std::runtime_error, " there is no THP constraint/limit for well " << name()
@@ -407,11 +407,11 @@ namespace Opm
template<typename TypeTag>
int
WellInterface<TypeTag>::
getTHPControlIndex() const
getControlIndex(const WellControlType& type) 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) {
if (well_controls_iget_type(well_controls_, ctrl_index) == type) {
return ctrl_index;
}
}
@@ -934,7 +934,7 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::
wellTesting(Simulator& simulator, const std::vector<double>& B_avg,
wellTesting(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
const WellTestConfig::Reason testing_reason,
/* const */ WellState& well_state,
@@ -959,7 +959,7 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::
wellTestingEconomic(Simulator& simulator, const std::vector<double>& B_avg,
wellTestingEconomic(const Simulator& simulator, const std::vector<double>& B_avg,
const double simulation_time, const int report_step,
const WellState& well_state, WellTestState& welltest_state, Opm::DeferredLogger& deferred_logger)
{
@@ -1166,7 +1166,7 @@ namespace Opm
template<typename TypeTag>
bool
WellInterface<TypeTag>::
solveWellEqUntilConverged(Simulator& ebosSimulator,
solveWellEqUntilConverged(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger)
@@ -1238,7 +1238,7 @@ namespace Opm
template<typename TypeTag>
void
WellInterface<TypeTag>::
solveWellForTesting(Simulator& ebosSimulator, WellState& well_state,
solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state,
const std::vector<double>& B_avg,
Opm::DeferredLogger& deferred_logger)
{