OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-09 10:15:34 -06:00
Code Issues Packages Projects Releases Wiki Activity

Add well potential output to MSW

Browse Source 
The common part is moved to the well interface class
This commit also adds iterations to the standardWells
well potential calculations to improve the results
This commit is contained in:
Tor Harald Sandve 2019-04-23 13:30:12 +02:00
parent a618a1a777
commit 821794b0ad
11 changed files with 465 additions and 173 deletions
Show Stats 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);

11
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;
@ -337,9 +337,11 @@ namespace Opm {
// xw to update Well State
void recoverWellSolutionAndUpdateWellState(const BVector& x);
void updateWellControls(Opm::DeferredLogger& deferred_logger);
void updateWellControls(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger);
void updateGroupControls(Opm::DeferredLogger& deferred_logger);
void updateGroupControls(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger);
// setting the well_solutions_ based on well_state.
void updatePrimaryVariables(Opm::DeferredLogger& deferred_logger);
@ -387,7 +389,8 @@ namespace Opm {
// some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
void prepareTimeStep(Opm::DeferredLogger& deferred_logger);
void prepareTimeStep(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger);
void prepareGroupControl(Opm::DeferredLogger& deferred_logger);

186
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() ) );
}
}
@ -679,18 +681,18 @@ namespace Opm {
int exception_thrown = 0;
try {
if (iterationIdx == 0) {
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(local_deferredLogger);
}
updateWellControls(local_deferredLogger);
// Set the well primary variables based on the value of well solutions
initPrimaryVariablesEvaluation();
std::vector< Scalar > B_avg(numComponents(), Scalar() );
computeAverageFormationFactor(B_avg);
if (iterationIdx == 0) {
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(B_avg, local_deferredLogger);
}
updateWellControls(B_avg, local_deferredLogger);
// Set the well primary variables based on the value of well solutions
initPrimaryVariablesEvaluation();
if (param_.solve_welleq_initially_ && iterationIdx == 0) {
// solve the well equations as a pre-processing step
last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
@ -699,7 +701,7 @@ namespace Opm {
if (initial_step_) {
// update the explicit quantities to get the initial fluid distribution in the well correct.
calculateExplicitQuantities(local_deferredLogger);
prepareTimeStep(local_deferredLogger);
prepareTimeStep(B_avg, local_deferredLogger);
last_report_ = solveWellEq(B_avg, dt, local_deferredLogger);
initial_step_ = false;
}
@ -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;
@ -921,7 +923,7 @@ namespace Opm {
// are active wells anywhere in the global domain.
if( wellsActive() )
{
updateWellControls(deferred_logger);
updateWellControls(B_avg, deferred_logger);
initPrimaryVariablesEvaluation();
}
} catch (std::exception& e) {
@ -1025,7 +1027,8 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateWellControls(Opm::DeferredLogger& deferred_logger)
updateWellControls(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger)
{
// Even if there are no wells active locally, we cannot
// return as the DeferredLogger uses global communication.
@ -1033,10 +1036,10 @@ namespace Opm {
if( !wellsActive() ) return ;
for (const auto& well : well_container_) {
well->updateWellControl(ebosSimulator_, well_state_, deferred_logger);
well->updateWellControl(ebosSimulator_, B_avg, well_state_, deferred_logger);
}
updateGroupControls(deferred_logger);
updateGroupControls(B_avg, deferred_logger);
}
@ -1066,17 +1069,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 +1146,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]);
@ -1111,7 +1171,8 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
prepareTimeStep(Opm::DeferredLogger& deferred_logger)
prepareTimeStep(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger)
{
if ( wellCollection().havingVREPGroups() ) {
@ -1132,7 +1193,7 @@ namespace Opm {
int exception_thrown = 0;
try {
for (const auto& well : well_container_) {
well->checkWellOperability(ebosSimulator_, well_state_, deferred_logger);
well->checkWellOperability(ebosSimulator_, B_avg, well_state_, deferred_logger);
}
// since the controls are all updated, we should update well_state accordingly
for (const auto& well : well_container_) {
@ -1144,7 +1205,7 @@ namespace Opm {
if (!well->isOperable() ) continue;
if (well_state_.effectiveEventsOccurred(w) ) {
well->updateWellStateWithTarget(ebosSimulator_, well_state_, deferred_logger);
well->updateWellStateWithTarget(ebosSimulator_, B_avg, well_state_, deferred_logger);
}
// there is no new well control change input within a report step,
@ -1372,7 +1433,8 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateGroupControls(Opm::DeferredLogger& deferred_logger)
updateGroupControls(const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger)
{
if (wellCollection().groupControlActive()) {
@ -1399,7 +1461,7 @@ namespace Opm {
// TODO: we should only do the well is involved in the update group targets
for (auto& well : well_container_) {
well->updateWellStateWithTarget(ebosSimulator_, well_state_, deferred_logger);
well->updateWellStateWithTarget(ebosSimulator_, B_avg, well_state_, deferred_logger);
well->updatePrimaryVariables(well_state_, deferred_logger);
}
}

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

@ -120,8 +120,9 @@ namespace Opm
/// updating the well state based the current control mode
virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const override;
Opm::DeferredLogger& deferred_logger) override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
@ -139,6 +140,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 +307,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 +331,13 @@ namespace Opm
const int perf,
std::vector<EvalWell>& mob) const;
virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) override;
void assembleControlEq(Opm::DeferredLogger& deferred_logger) const;
void assemblePressureEq(const int seg) const;
@ -345,6 +356,7 @@ namespace Opm
// checking the operability of the well based on current reservoir condition
// it is not implemented for multisegment well yet
virtual void checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override;
@ -366,7 +378,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;

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

@ -260,8 +260,9 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
updateWellStateWithTarget(const Simulator& /* ebos_simulator */,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& /* deferred_logger */) const
Opm::DeferredLogger& /* deferred_logger */)
{
// Updating well state bas on well control
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
@ -541,24 +542,73 @@ 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());
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, 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>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger)
{
// store a copy of the well state, we don't want to update the real well state
WellState copy = ebosSimulator.problem().wellModel().wellState();
initPrimaryVariablesEvaluation();
if (iterate) {
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();
}
}
@ -743,7 +793,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 +823,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 +1011,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 +1129,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;
}
}
}
@ -1673,6 +1748,7 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
checkWellOperability(const Simulator& /* ebos_simulator */,
const std::vector<Scalar>& /*B_avg */,
const WellState& /* well_state */,
Opm::DeferredLogger& deferred_logger)
{
@ -1829,7 +1905,7 @@ namespace Opm
updateWellState(dx_well, well_state, deferred_logger, relaxation_factor);
// TODO: should we do something more if a switching of control happens
this->updateWellControl(ebosSimulator, well_state, deferred_logger);
this->updateWellControl(ebosSimulator, B_avg, well_state, deferred_logger);
initPrimaryVariablesEvaluation();
}
@ -1877,6 +1953,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 +2023,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 +2083,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)
{

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

@ -148,8 +148,9 @@ namespace Opm
Opm::DeferredLogger& deferred_logger) override;
virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const override;
Opm::DeferredLogger& deferred_logger) override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg,
@ -168,6 +169,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 +329,18 @@ 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;
virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) override;
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;
@ -373,6 +376,7 @@ namespace Opm
// update the operability status of the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
virtual void checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
) override;
@ -380,24 +384,29 @@ namespace Opm
// 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 std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
);
// check whether the well is operable under BHP limit with current reservoir condition
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under THP limit with current reservoir condition
void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
// update WellState based on IPR and associated VFP table
void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
// for a well, when all drawdown are in the wrong direction, then this well will not
// be able to produce/inject .
@ -433,7 +442,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;

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

@ -149,8 +149,9 @@ namespace Opm
Opm::DeferredLogger& deferred_logger) override;
virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const override;
Opm::DeferredLogger& deferred_logger) override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
@ -168,6 +169,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 +332,18 @@ 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;
virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) override;
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;
@ -376,29 +379,35 @@ namespace Opm
// update the operability status of the well is operable under the current reservoir condition
// mostly related to BHP limit and THP limit
virtual void checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) override;
// 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 std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under BHP limit with current reservoir condition
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger);
// check whether the well is operable under THP limit with current reservoir condition
void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
// update WellState based on IPR and associated VFP table
void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
Opm::DeferredLogger& deferred_logger);
// for a well, when all drawdown are in the wrong direction, then this well will not
// be able to produce/inject .
@ -434,7 +443,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;

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

@ -480,7 +480,7 @@ namespace Opm
void
StandardWellV<TypeTag>::
assembleWellEq(const Simulator& ebosSimulator,
const std::vector<Scalar>& /* B_avg */,
const std::vector<Scalar>& B_avg,
const double dt,
WellState& well_state,
Opm::DeferredLogger& deferred_logger
@ -488,7 +488,7 @@ namespace Opm
{
// TODO: only_wells should be put back to save some computation
checkWellOperability(ebosSimulator, well_state, deferred_logger);
checkWellOperability(ebosSimulator, B_avg, well_state, deferred_logger);
if (!this->isOperable()) return;
@ -1194,8 +1194,9 @@ namespace Opm
void
StandardWellV<TypeTag>::
updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
// number of phases
const int np = number_of_phases_;
@ -1220,7 +1221,7 @@ namespace Opm
case THP: {
// when a well can not work under THP target, it switches to BHP control
if (this->operability_status_.isOperableUnderTHPLimit() ) {
updateWellStateWithTHPTargetIPR(ebos_simulator, well_state, deferred_logger);
updateWellStateWithTHPTargetIPR(ebos_simulator, B_avg, well_state, deferred_logger);
} else { // go to BHP limit
assert(this->operability_status_.isOperableUnderBHPLimit() );
@ -1398,6 +1399,7 @@ namespace Opm
void
StandardWellV<TypeTag>::
checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger)
{
@ -1412,7 +1414,7 @@ namespace Opm
const bool old_well_operable = this->operability_status_.isOperable();
updateWellOperability(ebos_simulator, well_state, deferred_logger);
updateWellOperability(ebos_simulator, B_avg, well_state, deferred_logger);
const bool well_operable = this->operability_status_.isOperable();
@ -1431,6 +1433,7 @@ namespace Opm
void
StandardWellV<TypeTag>::
updateWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& /* well_state */,
Opm::DeferredLogger& deferred_logger)
{
@ -1439,7 +1442,7 @@ namespace Opm
updateIPR(ebos_simulator, deferred_logger);
// checking the BHP limit related
checkOperabilityUnderBHPLimitProducer(ebos_simulator, deferred_logger);
checkOperabilityUnderBHPLimitProducer(ebos_simulator, B_avg, deferred_logger);
// checking whether the well can operate under the THP constraints.
if (this->wellHasTHPConstraints()) {
@ -1454,7 +1457,9 @@ namespace Opm
template<typename TypeTag>
void
StandardWellV<TypeTag>::
checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger)
{
const double bhp_limit = mostStrictBhpFromBhpLimits(deferred_logger);
// Crude but works: default is one atmosphere.
@ -1478,7 +1483,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, B_avg, bhp_limit, /*iterate=*/false, 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);
@ -1583,7 +1588,7 @@ namespace Opm
{
const double bhp = well_state.bhp()[index_of_well_];
std::vector<double> well_rates;
computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, bhp, /*iterate=*/ false, well_rates, deferred_logger);
const double sign = (well_type_ == PRODUCER) ? -1. : 1.;
const double threshold = sign * std::numeric_limits<double>::min();
@ -1626,11 +1631,13 @@ namespace Opm
void
StandardWellV<TypeTag>::
updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
if (well_type_ == PRODUCER) {
updateWellStateWithTHPTargetIPRProducer(ebos_simulator,
B_avg,
well_state,
deferred_logger);
}
@ -1650,8 +1657,9 @@ namespace Opm
void
StandardWellV<TypeTag>::
updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
well_state.thp()[index_of_well_] = this->getTHPConstraint(deferred_logger);
@ -1667,7 +1675,7 @@ namespace Opm
initPrimaryVariablesEvaluation();
std::vector<double> rates;
computeWellRatesWithBhp(ebos_simulator, EvalWell(numWellEq_ + numEq, bhp), rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, /*iterate=*/false, rates, deferred_logger);
// TODO: double checke the obtained rates
// this is another places we might obtain negative rates
@ -1690,7 +1698,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,12 +2271,42 @@ namespace Opm
void
StandardWellV<TypeTag>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
const int np = number_of_phases_;
well_flux.resize(np, 0.0);
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
if (iterate) {
// 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();
}
const bool allow_cf = getAllowCrossFlow();
@ -2278,17 +2316,24 @@ 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) {
well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value();
}
}
// reset bhp limit
well_controls_iset_target(wc, bhp_index, orig_bhp);
well_controls_set_current(wc, orig_current);
}
@ -2299,9 +2344,10 @@ namespace Opm
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 +2411,7 @@ namespace Opm
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), potentials, deferred_logger);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, potentials, deferred_logger);
// checking whether the potentials have valid values
for (const double value : potentials) {
@ -2403,6 +2449,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 +2471,7 @@ namespace Opm
if ( !wellHasTHPConstraints() ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
computeWellRatesWithBhp(ebosSimulator, EvalWell(numWellEq_ + numEq, bhp), well_potentials, deferred_logger);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, 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 +2483,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);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, 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 +2492,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 +2908,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)
{
@ -2884,7 +2931,7 @@ namespace Opm
// we should be able to provide a better initialization
calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
updateWellOperability(ebos_simulator, well_state_copy, deferred_logger);
updateWellOperability(ebos_simulator, B_avg, well_state_copy, deferred_logger);
if ( !this->isOperable() ) {
const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
@ -2892,7 +2939,7 @@ namespace Opm
return;
}
updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger);
updateWellStateWithTarget(ebos_simulator, B_avg, well_state_copy, deferred_logger);
calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
updatePrimaryVariables(well_state_copy, deferred_logger);

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

@ -444,13 +444,13 @@ namespace Opm
void
StandardWell<TypeTag>::
assembleWellEq(const Simulator& ebosSimulator,
const std::vector<Scalar>& /* B_avg */,
const std::vector<Scalar>& B_avg,
const double dt,
WellState& well_state,
Opm::DeferredLogger& deferred_logger)
{
checkWellOperability(ebosSimulator, well_state, deferred_logger);
checkWellOperability(ebosSimulator, B_avg, well_state, deferred_logger);
if (!this->isOperable()) return;
@ -1117,8 +1117,9 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
// number of phases
@ -1144,7 +1145,7 @@ namespace Opm
case THP: {
// when a well can not work under THP target, it switches to BHP control
if (this->operability_status_.isOperableUnderTHPLimit() ) {
updateWellStateWithTHPTargetIPR(ebos_simulator, well_state, deferred_logger);
updateWellStateWithTHPTargetIPR(ebos_simulator, B_avg, well_state, deferred_logger);
} else { // go to BHP limit
assert(this->operability_status_.isOperableUnderBHPLimit() );
@ -1322,6 +1323,7 @@ namespace Opm
void
StandardWell<TypeTag>::
checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger
)
@ -1337,7 +1339,7 @@ namespace Opm
const bool old_well_operable = this->operability_status_.isOperable();
updateWellOperability(ebos_simulator, well_state, deferred_logger);
updateWellOperability(ebos_simulator, B_avg, well_state, deferred_logger);
const bool well_operable = this->operability_status_.isOperable();
@ -1356,6 +1358,7 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& /* well_state */,
Opm::DeferredLogger& deferred_logger
)
@ -1365,7 +1368,7 @@ namespace Opm
updateIPR(ebos_simulator, deferred_logger);
// checking the BHP limit related
checkOperabilityUnderBHPLimitProducer(ebos_simulator, deferred_logger);
checkOperabilityUnderBHPLimitProducer(ebos_simulator, B_avg, deferred_logger);
// checking whether the well can operate under the THP constraints.
if (this->wellHasTHPConstraints()) {
@ -1380,7 +1383,9 @@ namespace Opm
template<typename TypeTag>
void
StandardWell<TypeTag>::
checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
Opm::DeferredLogger& deferred_logger)
{
const double bhp_limit = mostStrictBhpFromBhpLimits(deferred_logger);
// Crude but works: default is one atmosphere.
@ -1404,7 +1409,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, bhp_limit, well_rates_bhp_limit, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, B_avg, bhp_limit, /*iterate=*/ false, 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);
@ -1510,7 +1515,7 @@ namespace Opm
{
const double bhp = well_state.bhp()[index_of_well_];
std::vector<double> well_rates;
computeWellRatesWithBhp(ebos_simulator, bhp, well_rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, bhp, /*iterate=*/ false, well_rates, deferred_logger);
const double sign = (well_type_ == PRODUCER) ? -1. : 1.;
const double threshold = sign * std::numeric_limits<double>::min();
@ -1553,11 +1558,13 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
if (well_type_ == PRODUCER) {
updateWellStateWithTHPTargetIPRProducer(ebos_simulator,
B_avg,
well_state,
deferred_logger);
}
@ -1577,8 +1584,9 @@ namespace Opm
void
StandardWell<TypeTag>::
updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
well_state.thp()[index_of_well_] = this->getTHPConstraint(deferred_logger);
@ -1594,7 +1602,7 @@ namespace Opm
initPrimaryVariablesEvaluation();
std::vector<double> rates;
computeWellRatesWithBhp(ebos_simulator, bhp, rates, deferred_logger);
computeWellRatesWithBhp(ebos_simulator, B_avg, bhp, /*iterate=*/ false, rates, deferred_logger);
// TODO: double checke the obtained rates
// this is another places we might obtain negative rates
@ -1617,7 +1625,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);
@ -2157,12 +2165,43 @@ namespace Opm
void
StandardWell<TypeTag>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) const
Opm::DeferredLogger& deferred_logger)
{
const int np = number_of_phases_;
well_flux.resize(np, 0.0);
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
if (iterate) {
// 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();
}
const bool allow_cf = getAllowCrossFlow();
@ -2185,6 +2224,11 @@ namespace Opm
well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value();
}
}
// reset bhp limit
well_controls_iset_target(wc, bhp_index, orig_bhp);
well_controls_set_current(wc, orig_current);
}
@ -2195,9 +2239,10 @@ namespace Opm
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 +2306,7 @@ namespace Opm
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhp(ebosSimulator, bhp, potentials, deferred_logger);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, potentials, deferred_logger);
// checking whether the potentials have valid values
for (const double value : potentials) {
@ -2303,6 +2348,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
@ -2315,6 +2361,7 @@ namespace Opm
// TODO: for computeWellPotentials, no derivative is required actually
initPrimaryVariablesEvaluation();
const int np = number_of_phases_;
well_potentials.resize(np, 0.0);
@ -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);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, 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);
computeWellRatesWithBhp(ebosSimulator, B_avg, bhp, /*iterate=*/ true, 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)
@ -2778,7 +2825,7 @@ namespace Opm
// we should be able to provide a better initialization
calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
updateWellOperability(ebos_simulator, well_state_copy, deferred_logger);
updateWellOperability(ebos_simulator, B_avg, well_state_copy, deferred_logger);
if ( !this->isOperable() ) {
const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
@ -2786,7 +2833,7 @@ namespace Opm
return;
}
updateWellStateWithTarget(ebos_simulator, well_state_copy, deferred_logger);
updateWellStateWithTarget(ebos_simulator, B_avg, well_state_copy, deferred_logger);
calculateExplicitQuantities(ebos_simulator, well_state_copy, deferred_logger);
updatePrimaryVariables(well_state_copy, deferred_logger);

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

@ -182,15 +182,18 @@ 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;
virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) const = 0;
Opm::DeferredLogger& deferred_logger) = 0;
void updateWellControl(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) /* const */;
@ -235,7 +238,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,
@ -243,7 +246,10 @@ namespace Opm
void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
virtual void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0;
virtual void checkWellOperability(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) = 0;
// whether the well is operable
bool isOperable() const;
@ -354,7 +360,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 +391,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,17 +411,25 @@ 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);
void scaleProductivityIndex(const int perfIdx, double& productivity_index, const bool new_well, Opm::DeferredLogger& deferred_logger);
virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
const std::vector<Scalar>& B_avg,
const double& bhp,
const bool iterate,
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger) = 0;
// count the number of times an output log message is created in the productivity
// index calculations
int well_productivity_index_logger_counter_;

21
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;
}
}
@ -426,6 +426,7 @@ namespace Opm
void
WellInterface<TypeTag>::
updateWellControl(const Simulator& ebos_simulator,
const std::vector<Scalar>& B_avg,
WellState& well_state,
Opm::DeferredLogger& deferred_logger) /* const */
{
@ -504,7 +505,7 @@ namespace Opm
}
if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) {
updateWellStateWithTarget(ebos_simulator, well_state, deferred_logger);
updateWellStateWithTarget(ebos_simulator, B_avg, well_state, deferred_logger);
updatePrimaryVariables(well_state, deferred_logger);
}
}
@ -934,7 +935,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 +960,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 +1167,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)
@ -1189,7 +1190,7 @@ namespace Opm
++it;
solveEqAndUpdateWellState(well_state, deferred_logger);
updateWellControl(ebosSimulator, well_state, deferred_logger);
updateWellControl(ebosSimulator, B_avg, well_state, deferred_logger);
initPrimaryVariablesEvaluation();
} while (it < max_iter);
@ -1238,7 +1239,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)
{