OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-28 02:00:59 -06:00
Code Issues Packages Projects Releases Wiki Activity

Implement RESV control mode for producers

Browse Source 
This commit introduces support for the RESV control mode of
prediction (WCONPROD) and history-matching (WCONHIST) alike.  The
implementation uses class SurfaceToReservoirVoidage<> to compute
coefficients that convert component rates at surface conditions
(i.e., the primary degrees of freedom in the well residual) to phase
rates at reservoir condition.  The resulting coefficients can be
entered directly into system matrix of the linearised residual.

Note: We abuse the "distr" mechanism of struct WellControls to store
the conversion coefficients.  This may require refactorisation and
clarification at a later stage.  In the meantime, it allows for
transparent assembly of well equations--irrespective of surface- or
reservoir (voidage) rates.

Note: We do not yet support injectors controlled by total reservoir
voidage rate--either in history-matching (WCONINJH) or
prediction-scenario capacity.
This commit is contained in:
Bård Skaflestad 2014-07-07 16:51:48 +02:00
parent 23520be41a
commit 0f663dfe9f
2 changed files with 296 additions and 43 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp
View File

@ -944,37 +944,47 @@ namespace {
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
bool broken = false;
switch (well_type) {
case INJECTOR:
{
switch (ctrl_type) {
case BHP:
return bhp.value()[well] > target;
broken = bhp.value()[well] > target;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) > target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) > target;
}
break;
}
break;
case PRODUCER:
{
switch (ctrl_type) {
case BHP:
return bhp.value()[well] < target;
broken = bhp.value()[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
return rateToCompare(well_phase_flow_rate, well, num_phases, distr) < target;
case RESERVOIR_RATE:
// Intentional fall-through.
default:
OPM_THROW(std::logic_error, "Can only handle BHP and SURFACE_RATE controls.");
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) < target;
}
break;
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
return broken;
}
} // anonymous namespace
@ -1011,14 +1021,6 @@ namespace {
// inequality constraint, and therefore skipped.
continue;
}
if (well_controls_iget_type(wc, ctrl_index) == RESERVOIR_RATE) {
// We cannot handle this yet.
#ifdef OPM_VERBOSE
std::cout << "Warning: a RESERVOIR_RATE well control exists for well "
<< wells_.name[w] << ", but will never be checked." << std::endl;
#endif
continue;
}
if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
@ -1040,6 +1042,18 @@ namespace {
xw.bhp()[w] = target;
bhp_changed = true;
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.
//
// Just record the fact that we need to recompute
// the 'well_phase_flow_rate'.
rates_changed = true;
break;
case SURFACE_RATE:
for (int phase = 0; phase < np; ++phase) {
if (distr[phase] > 0.0) {
@ -1048,9 +1062,6 @@ namespace {
}
rates_changed = true;
break;
default:
OPM_THROW(std::logic_error, "Programmer error: should not have switched "
"to anything but BHP or SURFACE_RATE.");
}
}
}
@ -1078,13 +1089,11 @@ namespace {
const WellStateFullyImplicitBlackoil& xw,
const V& aliveWells)
{
// Handling BHP and SURFACE_RATE wells.
const int np = wells_.number_of_phases;
const int nw = wells_.number_of_wells;
V bhp_targets(nw);
V rate_targets(nw);
V bhp_targets = V::Zero(nw);
V rate_targets = V::Zero(nw);
M rate_distr(nw, np*nw);
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells_.ctrls[w];
@ -1092,20 +1101,33 @@ namespace {
// the current control set in the Wells struct, which
// is instead treated as a default.
const int current = xw.currentControls()[w];
if (well_controls_iget_type(wc, current) == BHP) {
bhp_targets[w] = well_controls_iget_target(wc, current);
rate_targets[w] = -1e100;
} else if (well_controls_iget_type(wc, current) == SURFACE_RATE) {
bhp_targets[w] = -1e100;
rate_targets[w] = well_controls_iget_target(wc, current);
{
const double * distr = well_controls_iget_distr(wc, current);
for (int phase = 0; phase < np; ++phase) {
rate_distr.insert(w, phase*nw + w) = distr[phase];
}
switch (well_controls_iget_type(wc, current)) {
case BHP:
{
bhp_targets (w) = well_controls_iget_target(wc, current);
rate_targets(w) = -1e100;
}
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
{
// RESERVOIR and SURFACE rates look the same, from a
// high-level point of view, in the system of
// simultaneous linear equations.
const double* const distr =
well_controls_iget_distr(wc, current);
for (int p = 0; p < np; ++p) {
rate_distr.insert(w, p*nw + w) = distr[p];
}
} else {
OPM_THROW(std::runtime_error, "Can only handle BHP and SURFACE_RATE type controls.");
bhp_targets (w) = -1.0e100;
rate_targets(w) = well_controls_iget_target(wc, current);
}
break;
}
}
const ADB bhp_residual = state.bhp - bhp_targets;

233
opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp
View File

@ -26,8 +26,11 @@
#include <opm/autodiff/FullyImplicitBlackoilSolver.hpp>
#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
#include <opm/autodiff/RateConverter.hpp>
#include <opm/core/grid.h>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
@ -43,15 +46,25 @@
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
#include <boost/filesystem.hpp>
#include <boost/lexical_cast.hpp>
#include <algorithm>
#include <cstddef>
#include <cassert>
#include <functional>
#include <memory>
#include <numeric>
#include <fstream>
#include <iostream>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm
{
@ -78,6 +91,10 @@ namespace Opm
private:
// Data.
typedef RateConverter::
SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
std::vector<int> > RateConverterType;
// Parameters for output.
bool output_;
bool output_vtk_;
@ -94,7 +111,13 @@ namespace Opm
const DerivedGeology &geo_;
FullyImplicitBlackoilSolver<Grid> solver_;
// Misc. data
RateConverterType rateConverter_;
std::vector<int> allcells_;
void
computeRESV(const std::size_t step,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw);
};
@ -209,6 +232,7 @@ namespace Opm
gravity_(gravity),
geo_(geo),
solver_(param, grid_, props_, geo_, rock_comp_props, *wells_, linsolver, has_disgas, has_vapoil)
, rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0))
{
// For output.
output_ = param.getDefault("output", true);
@ -280,6 +304,8 @@ namespace Opm
SimulatorReport sreport;
{
computeRESV(timer.currentStepNum(), state, well_state);
// Run solver.
solver_timer.start();
solver_.step(timer.currentStepLength(), state, well_state);
@ -328,5 +354,210 @@ namespace Opm
return report;
}
namespace SimFIBODetails {
typedef std::unordered_map<std::string, WellConstPtr> WellMap;
inline WellMap
mapWells(const std::vector<WellConstPtr>& wells)
{
WellMap wmap;
for (std::vector<WellConstPtr>::const_iterator
w = wells.begin(), e = wells.end();
w != e; ++w)
{
wmap.insert(std::make_pair((*w)->name(), *w));
}
return wmap;
}
inline int
resv_control(const WellControls* ctrl)
{
int i, n = well_controls_get_num(ctrl);
bool match = false;
for (i = 0; (! match) && (i < n); ++i) {
match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
}
if (! match) { i = 0; }
return i - 1; // -1 if no match, undo final "++" otherwise
}
inline bool
is_resv_prod(const Wells& wells,
const int w)
{
return ((wells.type[w] == PRODUCER) &&
(0 <= resv_control(wells.ctrls[w])));
}
inline bool
is_resv_prod(const WellMap& wmap,
const std::string& name,
const std::size_t step)
{
bool match = false;
WellMap::const_iterator i = wmap.find(name);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
match = (wp->isProducer(step) &&
wp->getProductionProperties(step)
.hasProductionControl(WellProducer::RESV));
}
return match;
}
inline std::vector<int>
resvProducers(const Wells& wells,
const std::size_t step,
const WellMap& wmap)
{
std::vector<int> resv_prod;
for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) {
if (is_resv_prod(wells, w) ||
((wells.name[w] != 0) &&
is_resv_prod(wmap, wells.name[w], step)))
{
resv_prod.push_back(w);
}
}
return resv_prod;
}
inline void
historyRates(const PhaseUsage& pu,
const WellProductionProperties& p,
std::vector<double>& rates)
{
assert (! p.predictionMode);
assert (rates.size() ==
std::vector<double>::size_type(pu.num_phases));
if (pu.phase_used[ BlackoilPhases::Aqua ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Aqua ];
rates[i] = p.WaterRate;
}
if (pu.phase_used[ BlackoilPhases::Liquid ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Liquid ];
rates[i] = p.OilRate;
}
if (pu.phase_used[ BlackoilPhases::Vapour ]) {
const std::vector<double>::size_type
i = pu.phase_pos[ BlackoilPhases::Vapour ];
rates[i] = p.GasRate;
}
}
} // namespace SimFIBODetails
template <class T>
void
SimulatorFullyImplicitBlackoil<T>::
Impl::computeRESV(const std::size_t step,
const BlackoilState& x,
WellStateFullyImplicitBlackoil& xw)
{
typedef SimFIBODetails::WellMap WellMap;
const std::vector<WellConstPtr>& w_ecl = schedule_->getWells(step);
const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
const std::vector<int>& resv_prod =
SimFIBODetails::resvProducers(*wells_, step, wmap);
if (! resv_prod.empty()) {
const PhaseUsage& pu = props_.phaseUsage();
const std::vector<double>::size_type np = props_.numPhases();
rateConverter_.defineState(x);
std::vector<double> distr (np);
std::vector<double> hrates(np);
std::vector<double> prates(np);
for (std::vector<int>::const_iterator
rp = resv_prod.begin(), e = resv_prod.end();
rp != e; ++rp)
{
WellControls* ctrl = wells_->ctrls[*rp];
// RESV control mode, all wells
{
const int rctrl = SimFIBODetails::resv_control(ctrl);
if (0 <= rctrl) {
const std::vector<double>::size_type off = (*rp) * np;
// Convert to positive rates to avoid issues
// in coefficient calculations.
std::transform(xw.wellRates().begin() + (off + 0*np),
xw.wellRates().begin() + (off + 1*np),
prates.begin(), std::negate<double>());
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(prates, fipreg, distr);
well_controls_iset_distr(ctrl, rctrl, & distr[0]);
}
}
// RESV control, WCONHIST wells. A bit of duplicate
// work, regrettably.
if (wells_->name[*rp] != 0) {
WellMap::const_iterator i = wmap.find(wells_->name[*rp]);
if (i != wmap.end()) {
WellConstPtr wp = i->second;
const WellProductionProperties& p =
wp->getProductionProperties(step);
if (! p.predictionMode) {
// History matching (WCONHIST/RESV)
SimFIBODetails::historyRates(pu, p, hrates);
const int fipreg = 0; // Hack. Ignore FIP regions.
rateConverter_.calcCoeff(hrates, fipreg, distr);
// WCONHIST/RESV target is sum of all
// observed phase rates translated to
// reservoir conditions. Recall sign
// convention: Negative for producers.
const double target =
- std::inner_product(distr.begin(), distr.end(),
hrates.begin(), 0.0);
well_controls_clear(ctrl);
well_controls_assert_number_of_phases(ctrl, int(np));
const int ok =
well_controls_add_new(RESERVOIR_RATE, target,
& distr[0], ctrl);
if (ok != 0) {
xw.currentControls()[*rp] = 0;
well_controls_set_current(ctrl, 0);
}
}
}
}
}
}
}
} // namespace Opm