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 #1118 from GitPaean/calculating_well_potentials_for_each_timestep

Browse Source 
[WIP] Changes to improve the StandardWellsDense
This commit is contained in:
Andreas Lauser 2017-04-12 10:33:57 +02:00 committed by GitHub
commit d334fc1c70
20 changed files with 663 additions and 317 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
opm/autodiff/BlackoilModelBase.hpp
View File

@ -558,6 +558,11 @@ namespace Opm {
std::vector<double>& maxNormWell,
int nc) const;
/// Set up the group control related at the beginning of each time step
void
setupGroupControl(const ReservoirState& reservoir_state,
WellState& well_state);
double dpMaxRel() const { return param_.dp_max_rel_; }
double dbhpMaxRel() const {return param_.dbhp_max_rel_; }
double dsMax() const { return param_.ds_max_; }

83
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -764,6 +764,11 @@ typedef Eigen::Array<double,
asImpl().wellModel().computeWellConnectionPressures(state0, well_state);
}
// set up the guide rate and group control
if (asImpl().wellModel().wellCollection()->groupControlActive() && initial_assembly) {
setupGroupControl(reservoir_state, well_state);
}
// Possibly switch well controls and updating well state to
// get reasonable initial conditions for the wells
asImpl().wellModel().updateWellControls(well_state);
@ -820,12 +825,6 @@ typedef Eigen::Array<double,
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
if (param_.compute_well_potentials_) {
SolutionState state0 = state;
asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
return report;
}
@ -2551,6 +2550,78 @@ typedef Eigen::Array<double,
}
}
template <class Grid, class WellModel, class Implementation>
void
BlackoilModelBase<Grid, WellModel, Implementation>::
setupGroupControl(const ReservoirState& reservoir_state,
WellState& well_state)
{
if (asImpl().wellModel().wellCollection()->requireWellPotentials()) {
SolutionState state = asImpl().variableState(reservoir_state, well_state);
asImpl().makeConstantState(state);
asImpl().wellModel().computeWellConnectionPressures(state, well_state);
const int np = numPhases();
std::vector<ADB> b(np, ADB::null());
std::vector<ADB> mob(np, ADB::null());
const ADB& press = state.pressure;
const ADB& temp = state.temperature;
const ADB& rs = state.rs;
const ADB& rv = state.rv;
const std::vector<PhasePresence> cond = phaseCondition();
const ADB pv_mult = poroMult(press);
const ADB tr_mult = transMult(press);
const std::vector<ADB> kr = asImpl().computeRelPerm(state);
std::vector<ADB> mob_perfcells(np, ADB::null());
std::vector<ADB> b_perfcells(np, ADB::null());
for (int phase = 0; phase < np; ++phase) {
if (active_[phase]) {
const std::vector<int>& well_cells = asImpl().wellModel().wellOps().well_cells;
const ADB mu = asImpl().fluidViscosity(canph_[phase], state.canonical_phase_pressures[canph_[phase]],
temp, rs, rv, cond);
mob[phase] = tr_mult * kr[canph_[phase]] / mu;
mob_perfcells[phase] = subset(mob[phase], well_cells);
b[phase] = asImpl().fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond);
b_perfcells[phase] = subset(b[phase], well_cells);
}
}
// well potentials for each well
std::vector<double> well_potentials;
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, well_state, state, well_potentials);
asImpl().wellModel().wellCollection()->setGuideRatesWithPotentials(asImpl().wellModel().wellsPointer(),
fluid_.phaseUsage(), well_potentials);
} // end of if
applyVREPGroupControl(reservoir_state, well_state);
if (asImpl().wellModel().wellCollection()->groupControlApplied()) {
asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
} else {
asImpl().wellModel().wellCollection()->applyGroupControls();
// the well collections do not have access to Well State, so the currentControls() of Well State need to
// be updated based on the group control setup
const int nw = wells().number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellNode& well_node = asImpl().wellModel().wellCollection()->findWellNode(wells().name[w]);
if (!well_node.individualControl()) {
well_state.currentControls()[w] = well_node.groupControlIndex();
}
}
} // end of else
}
} // namespace Opm
#endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED

5
opm/autodiff/BlackoilModelEbos.hpp
View File

@ -212,9 +212,12 @@ namespace Opm {
/// \param[in, out] reservoir_state reservoir state variables
/// \param[in, out] well_state well state variables
void prepareStep(const SimulatorTimerInterface& /*timer*/,
const ReservoirState& /*reservoir_state*/,
const ReservoirState& reservoir_state,
const WellState& /* well_state */)
{
if ( wellModel().wellCollection()->havingVREPGroups() ) {
updateRateConverter(reservoir_state);
}
}

2
opm/autodiff/BlackoilModelParameters.cpp
View File

@ -54,7 +54,6 @@ namespace Opm
param.getDefault("max_single_precision_days", unit::convert::to( maxSinglePrecisionTimeStep_, unit::day) ), unit::day );
solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
compute_well_potentials_ = param.getDefault("compute_well_potentials", compute_well_potentials_);
use_update_stabilization_ = param.getDefault("use_update_stabilization", use_update_stabilization_);
deck_file_name_ = param.template get<std::string>("deck_filename");
}
@ -78,7 +77,6 @@ namespace Opm
maxSinglePrecisionTimeStep_ = unit::convert::from( 20.0, unit::day );
solve_welleq_initially_ = true;
update_equations_scaling_ = false;
compute_well_potentials_ = false;
use_update_stabilization_ = true;
}

4
opm/autodiff/BlackoilModelParameters.hpp
View File

@ -62,10 +62,6 @@ namespace Opm
/// Update scaling factors for mass balance equations
bool update_equations_scaling_;
/// Compute well potentials, needed to calculate default guide rates for group
/// controlled wells
bool compute_well_potentials_;
/// Try to detect oscillation or stagnation.
bool use_update_stabilization_;

9
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -150,10 +150,13 @@ namespace Opm {
// get reasonable initial conditions for the wells
wellModel().updateWellControls(well_state);
// enforce VREP control when necessary.
Base::applyVREPGroupControl(reservoir_state, well_state);
// TODO: I do not think the multi_segment well can handle group control yet
if (asImpl().wellModel().wellCollection()->groupControlActive()) {
// enforce VREP control when necessary.
Base::applyVREPGroupControl(reservoir_state, well_state);
asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
asImpl().wellModel().wellCollection()->updateWellTargets(well_state.wellRates());
}
// Create the primary variables.
SolutionState state = asImpl().variableState(reservoir_state, well_state);

6
opm/autodiff/BlackoilTransportModel.hpp
View File

@ -149,12 +149,6 @@ namespace Opm {
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
if (param_.compute_well_potentials_) {
SolutionState state0 = state;
asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
}
return report;
}

1
opm/autodiff/ParallelDebugOutput.hpp
View File

@ -642,7 +642,6 @@ namespace Opm
// as we get the following error otherwise
// with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
// converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructor
std::vector<double>(),
std::unordered_set<std::string>());
const Wells* wells = wells_manager.c_wells();

3
opm/autodiff/SimulatorBase.hpp
View File

@ -175,9 +175,6 @@ namespace Opm
void
outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg);
void computeWellPotentials(const Wells* wells,
const WellState& xw,
std::vector<double>& well_potentials);
void updateListEconLimited(const std::unique_ptr<Solver>& solver,
const Schedule& schedule,

25
opm/autodiff/SimulatorBase_impl.hpp
View File

@ -138,8 +138,6 @@ namespace Opm
desiredRestoreStep );
}
bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
SimulatorReport report;
SimulatorReport stepReport;
@ -178,7 +176,6 @@ namespace Opm
Opm::UgGridHelpers::beginFaceCentroids(grid_),
dynamic_list_econ_limited,
is_parallel_run_,
well_potentials,
defunct_well_names_);
const Wells* wells = wells_manager.c_wells();
WellState well_state;
@ -328,11 +325,6 @@ namespace Opm
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;
// The well potentials are only computed if they are needed
// For now thay are only used to determine default guide rates for group controlled wells
if ( is_well_potentials_computed ) {
asImpl().computeWellPotentials(wells, well_state, well_potentials);
}
asImpl().updateListEconLimited(solver, eclipse_state_->getSchedule(), timer.currentStepNum(), wells,
well_state, dynamic_list_econ_limited);
@ -491,23 +483,6 @@ namespace Opm
return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
}
template <class Implementation>
void SimulatorBase<Implementation>::computeWellPotentials(const Wells* wells,
const WellState& xw,
std::vector<double>& well_potentials)
{
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
well_potentials.clear();
well_potentials.resize(nw*np,0.0);
for (int w = 0; w < nw; ++w) {
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
for (int phase = 0; phase < np; ++phase) {
well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
}
}
}
}
template <class Implementation>
void SimulatorBase<Implementation>::computeRESV(const std::size_t step,

25
opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp
View File

@ -199,8 +199,6 @@ public:
desiredRestoreStep );
}
bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
std::vector<double> well_potentials;
DynamicListEconLimited dynamic_list_econ_limited;
SimulatorReport report;
SimulatorReport stepReport;
@ -240,7 +238,6 @@ public:
Opm::UgGridHelpers::beginFaceCentroids(grid()),
dynamic_list_econ_limited,
is_parallel_run_,
well_potentials,
defunct_well_names_ );
const Wells* wells = wells_manager.c_wells();
WellState well_state;
@ -395,11 +392,6 @@ public:
report.output_write_time += perfTimer.stop();
prev_well_state = well_state;
// The well potentials are only computed if they are needed
// For now thay are only used to determine default guide rates for group controlled wells
if ( is_well_potentials_computed ) {
computeWellPotentials(wells, well_state, well_potentials);
}
updateListEconLimited(solver, eclState().getSchedule(), timer.currentStepNum(), wells,
well_state, dynamic_list_econ_limited);
@ -605,23 +597,6 @@ protected:
}
void computeWellPotentials(const Wells* wells,
const WellState& xw,
std::vector<double>& well_potentials)
{
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
well_potentials.clear();
well_potentials.resize(nw*np,0.0);
for (int w = 0; w < nw; ++w) {
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
for (int phase = 0; phase < np; ++phase) {
well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
}
}
}
}
void updateListEconLimited(const std::unique_ptr<Solver>& solver,
const Schedule& schedule,

1
opm/autodiff/SimulatorFullyImplicitBlackoilMultiSegment_impl.hpp
View File

@ -126,7 +126,6 @@ namespace Opm
// as we get the following error otherwise
// with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
// converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructor
std::vector<double>(), // null well_potentials
Base::defunct_well_names_);
const Wells* wells = wells_manager.c_wells();
WellState well_state;

1
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp
View File

@ -446,7 +446,6 @@ namespace Opm
// with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
// converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructo
, false,
std::vector<double>(),
std::unordered_set<std::string>());
const Wells* wells = wellsmanager.c_wells();

3
opm/autodiff/StandardWells.hpp
View File

@ -152,8 +152,9 @@ namespace Opm {
void
computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const WellState& well_state,
SolutionState& state0,
WellState& well_state);
std::vector<double>& well_potentials) const;
template <class SolutionState>
void

36
opm/autodiff/StandardWellsDense.hpp
View File

@ -155,7 +155,7 @@ enum WellVariablePositions {
int numCells() const;
void resetWellControlFromState(WellState xw) const;
void resetWellControlFromState(const WellState& xw) const;
const Wells& wells() const;
@ -243,13 +243,19 @@ enum WellVariablePositions {
const double grav);
// TODO: Later we might want to change the function to only handle one well,
// the requirement for well potential calculation can be based on individual wells.
// getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
// Calculating well potentials for each well
// TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
template<typename Simulator>
void
computeWellPotentials(const Simulator& ebosSimulator,
WellState& well_state) const;
const WellState& well_state,
std::vector<double>& well_potentials) const;
// TODO: some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
template<typename Simulator>
void prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state);
WellCollection* wellCollection() const;
@ -357,6 +363,26 @@ enum WellVariablePositions {
const int well_index,
WellState& xw) const;
bool wellHasTHPConstraints(const int well_index) const;
// TODO: maybe we should provide a light version of computeWellFlux, which does not include the
// calculation of the derivatives
template<typename Simulator>
void
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const int well_index,
std::vector<double>& well_flux) const;
double mostStrictBhpFromBhpLimits(const int well_index) const;
// TODO: maybe it should be improved to be calculate general rates for THP control later
template<typename Simulator>
std::vector<double>
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const int well_index,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential) const;
};

658
opm/autodiff/StandardWellsDense_impl.hpp
View File

@ -129,16 +129,16 @@ namespace Opm {
const double dt,
WellState& well_state)
{
if (iterationIdx == 0) {
prepareTimeStep(ebosSimulator, well_state);
}
SimulatorReport report;
if ( ! localWellsActive() ) {
return report;
}
if (param_.compute_well_potentials_) {
computeWellPotentials(ebosSimulator, well_state);
}
resetWellControlFromState(well_state);
updateWellControls(well_state);
// Set the primary variables for the wells
setWellVariables(well_state);
@ -260,7 +260,7 @@ namespace Opm {
{
const int np = wells().number_of_phases;
assert (mob.size() == np);
assert (int(mob.size()) == np);
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
@ -549,7 +549,7 @@ namespace Opm {
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
void
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
resetWellControlFromState(WellState xw) const
resetWellControlFromState(const WellState& xw) const
{
const int nw = wells_->number_of_wells;
for (int w = 0; w < nw; ++w) {
@ -765,7 +765,7 @@ namespace Opm {
EvalWell well_pressure = bhp + cdp;
EvalWell drawdown = pressure - well_pressure;
// injection perforations
// producing perforations
if ( drawdown.value() > 0 ) {
//Do nothing if crossflow is not allowed
if (!allow_cf && wells().type[w] == INJECTOR) {
@ -902,6 +902,11 @@ namespace Opm {
if (!converged) {
well_state = well_state0;
// also recover the old well controls
for (int w = 0; w < nw; ++w) {
WellControls* wc = wells().ctrls[w];
well_controls_set_current(wc, well_state.currentControls()[w]);
}
}
SimulatorReport report;
@ -1286,7 +1291,11 @@ namespace Opm {
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
const double* distr = well_controls_iget_distr(wc, current);
for (int p = 0; p < np; ++p) {
F[p] /= distr[p];
if (distr[p] > 0.) { // For injection wells, there only one non-zero distr value
F[p] /= distr[p];
} else {
F[p] = 0.;
}
}
} else {
for (int p = 0; p < np; ++p) {
@ -1407,7 +1416,8 @@ namespace Opm {
const WellControls* wc = wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
// Looping over all controls until we find a THP constraint
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
int ctrl_index = 0;
for ( ; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(wc, ctrl_index) == THP) {
// the current control
const int current = well_state.currentControls()[w];
@ -1460,6 +1470,11 @@ namespace Opm {
break;
}
} // end of for loop for seaching THP constraints
// no THP constraint found
if (ctrl_index == nwc) { // not finding a THP contstraints
well_state.thp()[w] = 0.0;
}
} // end of for (int w = 0; w < nw; ++w)
}
@ -1534,12 +1549,6 @@ namespace Opm {
}
}
// upate the well targets following group controls
if (wellCollection()->groupControlActive()) {
applyVREPGroupControl(xw);
wellCollection()->updateWellTargets(xw.wellRates());
}
// the new well control indices after all the related updates,
std::vector<int> updated_control_index(nw, 0);
for (int w = 0; w < nw; ++w) {
@ -1549,11 +1558,25 @@ namespace Opm {
// checking whether control changed
wellhelpers::WellSwitchingLogger logger;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
if (updated_control_index[w] != old_control_index[w]) {
WellControls* wc = wells().ctrls[w];
logger.wellSwitched(wells().name[w],
well_controls_iget_type(wc, old_control_index[w]),
well_controls_iget_type(wc, updated_control_index[w]));
}
if (updated_control_index[w] != old_control_index[w] || well_collection_->groupControlActive()) {
updateWellStateWithTarget(wc, updated_control_index[w], w, xw);
}
}
// upate the well targets following group controls
// it will not change the control mode, only update the targets
if (wellCollection()->groupControlActive()) {
applyVREPGroupControl(xw);
wellCollection()->updateWellTargets(xw.wellRates());
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
updateWellStateWithTarget(wc, updated_control_index[w], w, xw);
}
}
@ -1703,129 +1726,57 @@ namespace Opm {
void
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
computeWellPotentials(const Simulator& ebosSimulator,
WellState& well_state) const
const WellState& well_state,
std::vector<double>& well_potentials) const
{
// number of wells and phases
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
well_potentials.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) {
// bhp needs to be determined for the well potential calculation
// There can be more than one BHP/THP constraints.
// TODO: there is an option to ignore the THP limit when calculating well potentials,
// we are not handling it for the moment, while easy to incorporate
// the bhp will be used to compute well potentials
double bhp;
// get the bhp value based on the bhp constraints
const double bhp = mostStrictBhpFromBhpLimits(w);
// type of the well, INJECTOR or PRODUCER
const WellType& well_type = wells().type[w];
// initial bhp value, making the value not usable
switch(well_type) {
case INJECTOR:
bhp = std::numeric_limits<double>::max();
break;
case PRODUCER:
bhp = -std::numeric_limits<double>::max();
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[w]);
// does the well have a THP related constraint?
const bool has_thp_control = wellHasTHPConstraints(w);
std::vector<double> potentials(np);
if ( !has_thp_control ) {
assert(std::abs(bhp) != std::numeric_limits<double>::max());
computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials);
} 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
if ( !well_state.isNewWell(w) ) {
for (int p = 0; p < np; ++p) {
// This is dangerous for new added well
// since we are not handling the initialization correctly for now
potentials[p] = well_state.wellRates()[w * np + p];
}
} else {
// We need to generate a reasonable rates to start the iteration process
computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials);
for (double& value : 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.
const double rate_safety_scaling_factor = 0.00001;
value *= rate_safety_scaling_factor;
}
}
potentials = computeWellPotentialWithTHP(ebosSimulator, w, bhp, potentials);
}
// the well controls
const WellControls* well_control = wells().ctrls[w];
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_control);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
// finding a BHP constraint
if (well_controls_iget_type(well_control, ctrl_index) == BHP) {
// get the bhp constraint value, it should always be postive assummingly
const double bhp_target = well_controls_iget_target(well_control, ctrl_index);
switch(well_type) {
case INJECTOR: // using the lower bhp contraint from Injectors
if (bhp_target < bhp) {
bhp = bhp_target;
}
break;
case PRODUCER:
if (bhp_target > bhp) {
bhp = bhp_target;
}
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[w]);
} // end of switch
}
// finding a THP constraint
if (well_controls_iget_type(well_control, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
const Opm::PhaseUsage& pu = phase_usage_;
if (active_[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(well_control, ctrl_index);
const double& thp = well_controls_iget_target(well_control, ctrl_index);
const double& alq = well_controls_iget_alq(well_control, ctrl_index);
// Calculating the BHP value based on THP
const WellType& well_type = wells().type[w];
const int first_perf = wells().well_connpos[w]; //first perforation
if (well_type == INJECTOR) {
const double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities()[first_perf], gravity_);
const double bhp_calculated = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if (bhp_calculated < bhp) {
bhp = bhp_calculated;
}
}
else if (well_type == PRODUCER) {
const double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities()[first_perf], gravity_);
const double bhp_calculated = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictest of the bhp controlls i.e. largest bhp for producers
if (bhp_calculated > bhp) {
bhp = bhp_calculated;
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
// there should be always some avaible bhp/thp constraints there
assert(std::abs(bhp) != std::numeric_limits<double>::max());
// Should we consider crossflow when calculating well potentionals?
const bool allow_cf = allow_cross_flow(w, ebosSimulator);
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const int cell_index = wells().well_cells[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0));
std::vector<EvalWell> well_potentials(np, 0.0);
std::vector<EvalWell> mob(np, 0.0);
getMobility(ebosSimulator, perf, cell_index, mob);
computeWellFlux(w, wells().WI[perf], intQuants.fluidState(), mob, bhp, wellPerforationPressureDiffs()[perf], allow_cf, well_potentials);
for(int p = 0; p < np; ++p) {
well_state.wellPotentials()[perf * np + p] = well_potentials[p].value();
}
// putting the sucessfully calculated potentials to the well_potentials
for (int p = 0; p < np; ++p) {
well_potentials[w * np + p] = std::abs(potentials[p]);
}
} // end of for (int w = 0; w < nw; ++w)
}
@ -1834,6 +1785,93 @@ namespace Opm {
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
template<typename Simulator>
void
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
prepareTimeStep(const Simulator& ebos_simulator,
WellState& well_state)
{
const int nw = wells().number_of_wells;
for (int w = 0; w < nw; ++w) {
// after restarting, the well_controls can be modified while
// the well_state still uses the old control index
// we need to synchronize these two.
resetWellControlFromState(well_state);
if (wellCollection()->groupControlActive()) {
WellControls* wc = wells().ctrls[w];
WellNode& well_node = well_collection_->findWellNode(std::string(wells().name[w]));
// handling the situation that wells do not have a valid control
// it happens the well specified with GRUP and restarting due to non-convergencing
// putting the well under group control for this situation
int ctrl_index = well_controls_get_current(wc);
const int group_control_index = well_node.groupControlIndex();
if (group_control_index >= 0 && ctrl_index < 0) {
// put well under group control
well_controls_set_current(wc, group_control_index);
well_state.currentControls()[w] = group_control_index;
}
// Final step, update whehter the well is under group control or individual control
// updated ctrl_index from the well control
ctrl_index = well_controls_get_current(wc);
if (well_node.groupControlIndex() >= 0 && ctrl_index == well_node.groupControlIndex()) {
// under group control
well_node.setIndividualControl(false);
} else {
// individual control
well_node.setIndividualControl(true);
}
}
}
if (well_collection_->groupControlActive()) {
if (well_collection_->requireWellPotentials()) {
// calculate the well potentials
setWellVariables(well_state);
computeWellConnectionPressures(ebos_simulator, well_state);
// To store well potentials for each well
std::vector<double> well_potentials;
computeWellPotentials(ebos_simulator, well_state, well_potentials);
// update/setup guide rates for each well based on the well_potentials
well_collection_->setGuideRatesWithPotentials(wellsPointer(), phase_usage_, well_potentials);
}
applyVREPGroupControl(well_state);
if (!wellCollection()->groupControlApplied()) {
wellCollection()->applyGroupControls();
} else {
wellCollection()->updateWellTargets(well_state.wellRates());
}
}
// since the controls are all updated, we should update well_state accordingly
for (int w = 0; w < nw; ++w) {
WellControls* wc = wells().ctrls[w];
const int control = well_controls_get_current(wc);
well_state.currentControls()[w] = control;
updateWellStateWithTarget(wc, control, w, well_state);
// The wells are not considered to be newly added
// for next time step
if (well_state.isNewWell(w) ) {
well_state.setNewWell(w, false);
}
} // end of for (int w = 0; w < nw; ++w)
}
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
WellCollection*
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
@ -1956,11 +1994,14 @@ namespace Opm {
computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs);
wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
// for the wells under group control, update the currentControls for the well_state
// for the wells under group control, update the control index for the well_state and well_controls
for (const WellNode* well_node : wellCollection()->getLeafNodes()) {
if (well_node->isInjector() && !well_node->individualControl()) {
const int well_index = well_node->selfIndex();
well_state.currentControls()[well_index] = well_node->groupControlIndex();
WellControls* wc = wells().ctrls[well_index];
well_controls_set_current(wc, well_node->groupControlIndex());
}
}
}
@ -2167,7 +2208,13 @@ namespace Opm {
const WellControls* wc = wells().ctrls[wellIdx];
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
const double* distr = well_controls_get_current_distr(wc);
return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
if (distr[phaseIdx] > 0.) {
return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
} else {
// TODO: not sure why return EvalWell(0.) causing problem here
// Probably due to the wrong Jacobians.
return wellVolumeFraction(wellIdx, phaseIdx);
}
}
std::vector<double> g = {1,1,0.01};
return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]);
@ -2403,9 +2450,13 @@ namespace Opm {
switch (well_controls_iget_type(wc, current)) {
case BHP:
xw.bhp()[well_index] = target;
// TODO: similar to the way below to handle THP
// we should not something related to thp here when there is thp constraint
break;
case THP: {
xw.thp()[well_index] = target;
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
@ -2453,38 +2504,58 @@ namespace Opm {
break;
}
case RESERVOIR_RATE:
// No direct change to any observable quantity at
// surface condition. In this case, use existing
// flow rates as initial conditions as reservoir
// rate acts only in aggregate.
break;
case RESERVOIR_RATE: // intentional fall-through
case SURFACE_RATE:
// assign target value as initial guess for injectors and
// single phase producers (orat, grat, wrat)
// checking the number of the phases under control
int numPhasesWithTargetsUnderThisControl = 0;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
numPhasesWithTargetsUnderThisControl += 1;
}
}
assert(numPhasesWithTargetsUnderThisControl > 0);
const WellType& well_type = wells().type[well_index];
if (well_type == INJECTOR) {
// assign target value as initial guess for injectors
// only handles single phase control at the moment
assert(numPhasesWithTargetsUnderThisControl == 1);
for (int phase = 0; phase < np; ++phase) {
const double& compi = wells().comp_frac[np * well_index + phase];
// TODO: it was commented out from the master branch already.
//if (compi > 0.0) {
xw.wellRates()[np*well_index + phase] = target * compi;
//}
}
} else if (well_type == PRODUCER) {
// only set target as initial rates for single phase
// producers. (orat, grat and wrat, and not lrat)
// lrat will result in numPhasesWithTargetsUnderThisControl == 2
int numPhasesWithTargetsUnderThisControl = 0;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
numPhasesWithTargetsUnderThisControl += 1;
if (distr[phase] > 0.) {
xw.wellRates()[np*well_index + phase] = target / distr[phase];
} else {
xw.wellRates()[np * well_index + phase] = 0.;
}
}
} else if (well_type == PRODUCER) {
// update the rates of phases under control based on the target,
// and also update rates of phases not under control to keep the rate ratio,
// assuming the mobility ratio does not change for the production wells
double original_rates_under_phase_control = 0.0;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
xw.wellRates()[np*well_index + phase] = target * distr[phase];
if (distr[phase] > 0.0) {
original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase];
}
}
if (original_rates_under_phase_control != 0.0 ) {
double scaling_factor = target / original_rates_under_phase_control;
for (int phase = 0; phase < np; ++phase) {
xw.wellRates()[np * well_index + phase] *= scaling_factor;
}
} else { // scaling factor is not well defied when original_rates_under_phase_control is zero
// separating targets equally between phases under control
const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
} else {
// this only happens for SURFACE_RATE control
xw.wellRates()[np * well_index + phase] = target_rate_divided;
}
}
}
} else {
@ -2539,16 +2610,281 @@ namespace Opm {
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
}
} else {
if (active_[ Water ]) {
xw.wellSolutions()[WFrac*nw + well_index] = wells().comp_frac[np*well_index + Water];
}
const WellType& well_type = wells().type[well_index];
if (well_type == INJECTOR) {
// only single phase injection handled
if (active_[Water]) {
if (distr[Water] > 0.0) {
xw.wellSolutions()[WFrac * nw + well_index] = 1.0;
} else {
xw.wellSolutions()[WFrac * nw + well_index] = 0.0;
}
}
if (active_[ Gas ]) {
xw.wellSolutions()[GFrac*nw + well_index] = wells().comp_frac[np*well_index + Gas];
if (active_[Gas]) {
if (distr[Gas] > 0.0) {
xw.wellSolutions()[GFrac * nw + well_index] = 1.0;
} else {
xw.wellSolutions()[GFrac * nw + well_index] = 0.0;
}
}
// TODO: it is possible to leave injector as a oil well,
// when F_w and F_g both equals to zero, not sure under what kind of circumstance
// this will happen.
} else if (well_type == PRODUCER) { // producers
if (active_[Water]) {
xw.wellSolutions()[WFrac * nw + well_index] = 1.0 / np;
}
if (active_[Gas]) {
xw.wellSolutions()[GFrac * nw + well_index] = 1.0 / np;
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
bool
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
wellHasTHPConstraints(const int well_index) const
{
const WellControls* well_control = wells().ctrls[well_index];
const int nwc = well_controls_get_num(well_control);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_control, ctrl_index) == THP) {
return true;
}
}
return false;
}
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
template <typename Simulator>
void
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
computeWellRatesWithBhp(const Simulator& ebosSimulator,
const EvalWell& bhp,
const int well_index,
std::vector<double>& well_flux) const
{
const int np = wells().number_of_phases;
well_flux.resize(np, 0.0);
const bool allow_cf = allow_cross_flow(well_index, ebosSimulator);
for (int perf = wells().well_connpos[well_index]; perf < wells().well_connpos[well_index + 1]; ++perf) {
const int cell_index = wells().well_cells[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_index, /*timeIdx=*/ 0));
// flux for each perforation
std::vector<EvalWell> cq_s(np, 0.0);
std::vector<EvalWell> mob(np, 0.0);
getMobility(ebosSimulator, perf, cell_index, mob);
computeWellFlux(well_index, wells().WI[perf], intQuants.fluidState(), mob, bhp,
wellPerforationPressureDiffs()[perf], allow_cf, cq_s);
for(int p = 0; p < np; ++p) {
well_flux[p] += cq_s[p].value();
}
}
}
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
double
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
mostStrictBhpFromBhpLimits(const int well_index) const
{
double bhp;
// type of the well, INJECTOR or PRODUCER
const WellType& well_type = wells().type[well_index];
// initial bhp value, making the value not usable
switch(well_type) {
case INJECTOR:
bhp = std::numeric_limits<double>::max();
break;
case PRODUCER:
bhp = -std::numeric_limits<double>::max();
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[well_index]);
}
// the well controls
const WellControls* well_control = wells().ctrls[well_index];
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_control);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
// finding a BHP constraint
if (well_controls_iget_type(well_control, ctrl_index) == BHP) {
// get the bhp constraint value, it should always be postive assummingly
const double bhp_target = well_controls_iget_target(well_control, ctrl_index);
switch(well_type) {
case INJECTOR: // using the lower bhp contraint from Injectors
if (bhp_target < bhp) {
bhp = bhp_target;
}
break;
case PRODUCER:
if (bhp_target > bhp) {
bhp = bhp_target;
}
break;
default:
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << wells().name[well_index]);
} // end of switch
}
}
return bhp;
}
template<typename FluidSystem, typename BlackoilIndices, typename ElementContext, typename MaterialLaw>
template <typename Simulator>
std::vector<double>
StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
computeWellPotentialWithTHP(const Simulator& ebosSimulator,
const int well_index,
const double initial_bhp, // bhp from BHP constraints
const std::vector<double>& initial_potential) const
{
// 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?
const int np = wells().number_of_phases;
assert( np == int(initial_potential.size()) );
std::vector<double> potentials = initial_potential;
std::vector<double> old_potentials = potentials; // keeping track of the old potentials
double bhp = initial_bhp;
double old_bhp = bhp;
bool converged = false;
const int max_iteration = 1000;
const double bhp_tolerance = 1000.; // 1000 pascal
int iteration = 0;
while ( !converged && iteration < max_iteration ) {
// for each iteration, we calculate the bhp based on the rates/potentials with thp constraints
// with considering the bhp value from the bhp limits. At the beginning of each iteration,
// we initialize the bhp to be the bhp value from the bhp limits. Then based on the bhp values calculated
// from the thp constraints, we decide the effective bhp value for well potential calculation.
bhp = initial_bhp;
// the well controls
const WellControls* well_control = wells().ctrls[well_index];
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_control);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_control, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
const Opm::PhaseUsage& pu = phase_usage_;
if (active_[ Water ]) {
aqua = potentials[pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = potentials[pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = potentials[pu.phase_pos[ Gas ] ];
}
const int vfp = well_controls_iget_vfp(well_control, ctrl_index);
const double thp = well_controls_iget_target(well_control, ctrl_index);
const double alq = well_controls_iget_alq(well_control, ctrl_index);
// Calculating the BHP value based on THP
// TODO: check whether it is always correct to do calculation based on the depth of the first perforation.
const int first_perf = wells().well_connpos[well_index]; //first perforation
const WellType& well_type = wells().type[well_index];
if (well_type == INJECTOR) {
const double dp = wellhelpers::computeHydrostaticCorrection(
wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities()[first_perf], gravity_);
const double bhp_calculated = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if (bhp_calculated < bhp) {
bhp = bhp_calculated;
}
}
else if (well_type == PRODUCER) {
const double dp = wellhelpers::computeHydrostaticCorrection(
wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities()[first_perf], gravity_);
const double bhp_calculated = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictest of the bhp controlls i.e. largest bhp for producers
if (bhp_calculated > bhp) {
bhp = bhp_calculated;
}
} else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
// there should be always some available bhp/thp constraints there
if (std::isinf(bhp) || std::isnan(bhp)) {
OPM_THROW(std::runtime_error, "Unvalid bhp value obtained during the potential calculation for well " << wells().name[well_index]);
}
converged = std::abs(old_bhp - bhp) < bhp_tolerance;
computeWellRatesWithBhp(ebosSimulator, bhp, well_index, potentials);
// checking whether the potentials have valid values
for (const double value : potentials) {
if (std::isinf(value) || std::isnan(value)) {
OPM_THROW(std::runtime_error, "Unvalid potential value obtained during the potential calculation for well " << wells().name[well_index]);
}
}
if (!converged) {
old_bhp = bhp;
for (int p = 0; p < np; ++p) {
// TODO: improve the interpolation, will it always be valid with the way below?
// TODO: finding better paramters, better iteration strategy for better convergence rate.
const double potential_update_damping_factor = 0.001;
potentials[p] = potential_update_damping_factor * potentials[p] + (1.0 - potential_update_damping_factor) * old_potentials[p];
old_potentials[p] = potentials[p];
}
}
++iteration;
}
if (!converged) {
OPM_THROW(std::runtime_error, "Failed in getting converged for the potential calculation for well " << wells().name[well_index]);
}
return potentials;
}
} // namespace Opm

22
opm/autodiff/StandardWells_impl.hpp
View File

@ -1002,8 +1002,9 @@ namespace Opm
void
StandardWells::computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells,
const WellState& well_state,
SolutionState& state0,
WellState& well_state)
std::vector<double>& well_potentials) const
{
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
@ -1079,17 +1080,18 @@ namespace Opm
// compute well potentials
Vector aliveWells;
std::vector<ADB> well_potentials;
computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_potentials);
std::vector<ADB> perf_potentials;
computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, perf_potentials);
// store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw];
Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
well_potentials.resize(nw * np, 0.0);
for (int p = 0; p < np; ++p) {
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) {
well_potentials[w * np + p] += perf_potentials[p].value()[perf];
}
}
}
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
}

41
opm/autodiff/WellStateFullyImplicitBlackoil.hpp
View File

@ -103,8 +103,7 @@ namespace Opm
current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
}
well_potentials_.clear();
well_potentials_.resize(nperf * np, 0.0);
is_new_well_.resize(nw, true);
// intialize wells that have been there before
// order may change so the mapping is based on the well name
@ -117,12 +116,18 @@ namespace Opm
const_iterator it = prevState.wellMap().find( name );
if( it != end )
{
// this is not a new added well
is_new_well_[w] = false;
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
// bhp
bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
// thp
thp()[ newIndex ] = prevState.thp()[ oldIndex ];
// wellrates
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
@ -159,17 +164,6 @@ namespace Opm
perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
}
}
// currentControls
const int old_control_index = prevState.currentControls()[ oldIndex ];
if (old_control_index < well_controls_get_num(wells->ctrls[w])) {
// If the set of controls have changed, this may not be identical
// to the last control, but it must be a valid control.
currentControls()[ newIndex ] = old_control_index;
} else {
assert(well_controls_get_num(wells->ctrls[w]) > 0);
currentControls()[ newIndex ] = 0;
}
}
@ -210,10 +204,6 @@ namespace Opm
std::vector<int>& currentControls() { return current_controls_; }
const std::vector<int>& currentControls() const { return current_controls_; }
/// One rate per phase and well connection.
std::vector<double>& wellPotentials() { return well_potentials_; }
const std::vector<double>& wellPotentials() const { return well_potentials_; }
data::Wells report(const PhaseUsage &pu) const override {
data::Wells res = WellState::report(pu);
@ -265,10 +255,25 @@ namespace Opm
return res;
}
bool isNewWell(const int w) const {
return is_new_well_[w];
}
void setNewWell(const int w, const bool is_new_well) {
is_new_well_[w] = is_new_well;
}
private:
std::vector<double> perfphaserates_;
std::vector<int> current_controls_;
std::vector<double> well_potentials_;
// marking whether the well is just added
// for newly added well, the current initialized rates from WellState
// will have very wrong compositions for production wells, will mostly cause
// problem with VFP interpolation
std::vector<bool> is_new_well_;
};
} // namespace Opm

47
opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp
View File

@ -57,7 +57,6 @@ namespace Opm
using BaseType :: numPhases;
using BaseType :: perfPhaseRates;
using BaseType :: currentControls;
using BaseType :: wellPotentials;
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
@ -68,51 +67,15 @@ namespace Opm
// call init on base class
BaseType :: init(wells, state, prevState);
// TODO: the reason to keep this is to avoid getting defaulted value BHP
// some facilities needed from opm-parser or opm-core
// It is a little tricky, since sometimes before applying group control, the only
// available constraints in the well_controls is the defaulted BHP value, and it
// is really not desirable to use this value to enter the Newton iterations.
setWellSolutions(pu);
setWellSolutionsFromPrevState(prevState);
}
template <class PrevState>
void setWellSolutionsFromPrevState(const PrevState& prevState)
{
// Set nw and np, or return if no wells.
if (wells_.get() == nullptr) {
return;
}
const int nw = wells_->number_of_wells;
if (nw == 0) {
return;
}
const int np = wells_->number_of_phases;
// intialize wells that have been there before
// order may change so the mapping is based on the well name
// if there are no well, do nothing in init
if( ! prevState.wellMap().empty() )
{
typedef typename WellMapType :: const_iterator const_iterator;
const_iterator end = prevState.wellMap().end();
int nw_old = prevState.bhp().size();
for (int w = 0; w < nw; ++w) {
std::string name( wells_->name[ w ] );
const_iterator it = prevState.wellMap().find( name );
if( it != end )
{
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
// wellSolutions
for( int i = 0; i < np; ++i)
{
wellSolutions()[ i*nw + newIndex ] = prevState.wellSolutions()[i * nw_old + oldIndex ];
}
}
}
}
}
/// Set wellSolutions() based on the base class members.
void setWellSolutions(const PhaseUsage& pu)

3
tests/test_multisegmentwells.cpp
View File

@ -107,12 +107,11 @@ struct SetupMSW {
Opm::UgGridHelpers::cell2Faces(grid),
Opm::UgGridHelpers::beginFaceCentroids(grid),
dummy_dynamic_list,
false
false,
// We need to pass the optionaly arguments
// as we get the following error otherwise
// with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
// converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructor
, std::vector<double>(), // null well_potentials
std::unordered_set<std::string>());
const Wells* wells = wells_manager.c_wells();