first attempt to introduce StandardWells struct.

To replace the const Wells struct in BlackoilModelBase.
Only testing with Flow for the moment. Will update other Flow siblings
later.
This commit is contained in:
Kai Bao 2016-03-31 15:19:46 +02:00
parent b406be3839
commit 8367eaacb3
2 changed files with 103 additions and 67 deletions

View File

@ -276,13 +276,23 @@ namespace Opm {
std::vector<int> well_cells; // the set of perforated cells
};
struct StandardWells {
// keeping the underline, later they will be private members
StandardWells(const Wells* wells);
const Wells& wells() const;
bool wells_active_;
const Wells* wells_;
V well_perforation_densities_; //Density of each well perforation
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
};
// --------- Data members ---------
const Grid& grid_;
const BlackoilPropsAdInterface& fluid_;
const DerivedGeology& geo_;
const RockCompressibility* rock_comp_props_;
const Wells* wells_;
StandardWells std_wells_;
VFPProperties vfp_properties_;
const NewtonIterationBlackoilInterface& linsolver_;
// For each canonical phase -> true if active
@ -305,8 +315,6 @@ namespace Opm {
V isRs_;
V isRv_;
V isSg_;
V well_perforation_densities_; //Density of each well perforation
V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
LinearisedBlackoilResidual residual_;
@ -339,11 +347,15 @@ namespace Opm {
}
// return true if wells are available in the reservoir
bool wellsActive() const { return wells_active_; }
bool wellsActive() const { return std_wells_.wells_active_; }
// return true if wells are available on this process
bool localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; }
bool localWellsActive() const { return std_wells_.wells_ ? (std_wells_.wells_->number_of_wells > 0 ) : false; }
// return wells object
const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; }
// const Wells& wells () const { assert( bool(std_wells_.wells_ != 0) ); return *(std_wells_.wells_); }
// return the StandardWells object
StandardWells& stdWells() { return std_wells_; }
const StandardWells& stdWells() const { return std_wells_; }
int numWellVars() const;

View File

@ -159,14 +159,14 @@ namespace detail {
, fluid_ (fluid)
, geo_ (geo)
, rock_comp_props_(rock_comp_props)
, wells_ (wells_arg)
, std_wells_ (wells_arg)
, vfp_properties_(eclState->getTableManager()->getVFPInjTables(), eclState->getTableManager()->getVFPProdTables())
, linsolver_ (linsolver)
, active_(detail::activePhases(fluid.phaseUsage()))
, canph_ (detail::active2Canonical(fluid.phaseUsage()))
, cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid)))
, ops_ (grid, geo.nnc())
, wops_ (wells_)
, wops_ (wells_arg)
, has_disgas_(has_disgas)
, has_vapoil_(has_vapoil)
, param_( param )
@ -195,9 +195,9 @@ namespace detail {
// Only rank 0 does print to std::cout if terminal_output is enabled
terminal_output_ = (info.communicator().rank()==0);
}
int local_number_of_wells = wells_ ? wells_->number_of_wells : 0;
int local_number_of_wells = std_wells_.wells_ ? std_wells_.wells_->number_of_wells : 0;
int global_number_of_wells = info.communicator().sum(local_number_of_wells);
wells_active_ = ( wells_ && global_number_of_wells > 0 );
std_wells_.wells_active_ = ( std_wells_.wells_ && global_number_of_wells > 0 );
// Compute the global number of cells
std::vector<int> v( Opm::AutoDiffGrid::numCells(grid_), 1);
global_nc_ = 0;
@ -205,7 +205,7 @@ namespace detail {
}else
#endif
{
wells_active_ = ( wells_ && wells_->number_of_wells > 0 );
std_wells_.wells_active_ = ( std_wells_.wells_ && std_wells_.wells_->number_of_wells > 0 );
global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
}
}
@ -457,12 +457,36 @@ namespace detail {
template <class Grid, class Implementation>
BlackoilModelBase<Grid, Implementation>::
StandardWells::StandardWells(const Wells* wells)
: wells_(wells)
{
}
template <class Grid, class Implementation>
const Wells&
BlackoilModelBase<Grid, Implementation>::
StandardWells::wells() const
{
assert(wells_ != 0);
return *(wells_);
}
template <class Grid, class Implementation>
int
BlackoilModelBase<Grid, Implementation>::numWellVars() const
{
// For each well, we have a bhp variable, and one flux per phase.
const int nw = localWellsActive() ? wells().number_of_wells : 0;
const int nw = localWellsActive() ? stdWells().wells().number_of_wells : 0;
return (numPhases() + 1) * nw;
}
@ -582,8 +606,8 @@ namespace detail {
{
// Need to reshuffle well rates, from phase running fastest
// to wells running fastest.
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
const int np = stdWells().wells().number_of_phases;
// The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
@ -803,15 +827,15 @@ namespace detail {
std::vector<double>& rvmax_perf,
std::vector<double>& surf_dens_perf)
{
const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
const int nperf = stdWells().wells().well_connpos[stdWells().wells().number_of_wells];
const int nw = stdWells().wells().number_of_wells;
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
for (int perf = stdWells().wells().well_connpos[w]; perf < stdWells().wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == stdWells().wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
@ -891,7 +915,7 @@ namespace detail {
// 2. Compute densities
std::vector<double> cd =
WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_.phaseUsage(),
stdWells().wells(), xw, fluid_.phaseUsage(),
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
const int nperf = wells().well_connpos[wells().number_of_wells];
@ -908,11 +932,11 @@ namespace detail {
// 3. Compute pressure deltas
std::vector<double> cdp =
WellDensitySegmented::computeConnectionPressureDelta(
wells(), perf_depth, cd, grav);
stdWells().wells(), perf_depth, cd, grav);
// 4. Store the results
well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
stdWells().well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf);
stdWells().well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
}
@ -1159,15 +1183,15 @@ namespace detail {
{
if( ! localWellsActive() ) return ;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
const int nperf = stdWells().wells().well_connpos[nw];
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
V Tw = Eigen::Map<const V>(wells().WI, nperf);
V Tw = Eigen::Map<const V>(stdWells().wells().WI, nperf);
const std::vector<int>& well_cells = wops_.well_cells;
// pressure diffs computed already (once per step, not changing per iteration)
const V& cdp = well_perforation_pressure_diffs_;
const V& cdp = std_wells_.well_perforation_pressure_diffs_;
// Extract needed quantities for the perforation cells
const ADB& p_perfcells = subset(state.pressure, well_cells);
const ADB& rv_perfcells = subset(state.rv, well_cells);
@ -1197,15 +1221,15 @@ namespace detail {
const V numInjectingPerforations = (wops_.p2w * ADB::constant(selectInjectingPerforations)).value();
const V numProducingPerforations = (wops_.p2w * ADB::constant(selectProducingPerforations)).value();
for (int w = 0; w < nw; ++w) {
if (!wells().allow_cf[w]) {
for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) {
if (!stdWells().wells().allow_cf[w]) {
for (int perf = stdWells().wells().well_connpos[w] ; perf < stdWells().wells().well_connpos[w+1]; ++perf) {
// Crossflow is not allowed; reverse flow is prevented.
// At least one of the perforation must be open in order to have a meeningful
// equation to solve. For the special case where all perforations have reverse flow,
// and the target rate is non-zero all of the perforations are keept open.
if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
if (stdWells().wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) {
selectProducingPerforations[perf] = 0.0;
} else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
} else if (stdWells().wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){
selectInjectingPerforations[perf] = 0.0;
}
}
@ -1242,7 +1266,7 @@ namespace detail {
// and the well injection rates.
// compute avg. and total wellbore phase volumetric rates at standard conds
const DataBlock compi = Eigen::Map<const DataBlock>(wells().comp_frac, nw, np);
const DataBlock compi = Eigen::Map<const DataBlock>(stdWells().wells().comp_frac, nw, np);
std::vector<ADB> wbq(np, ADB::null());
ADB wbqt = ADB::constant(V::Zero(nw));
for (int phase = 0; phase < np; ++phase) {
@ -1313,9 +1337,9 @@ namespace detail {
}
// Update the perforation phase rates (used to calculate the pressure drop in the wellbore).
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int nperf = wells().well_connpos[nw];
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
const int nperf = stdWells().wells().well_connpos[nw];
V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np);
@ -1323,7 +1347,7 @@ namespace detail {
xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np);
// Update the perforation pressures.
const V& cdp = well_perforation_pressure_diffs_;
const V& cdp = std_wells_.well_perforation_pressure_diffs_;
const V perfpressure = (wops_.w2p * state.bhp.value().matrix()).array() + cdp;
xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf);
}
@ -1343,8 +1367,8 @@ namespace detail {
return;
}
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
ADB qs = state.qs;
for (int phase = 0; phase < np; ++phase) {
qs -= superset(wops_.p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np);
@ -1501,10 +1525,10 @@ namespace detail {
return false;
}
const int nw = wells().number_of_wells;
const int nw = stdWells().wells().number_of_wells;
//Loop over all wells
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const WellControls* wc = stdWells().wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
@ -1530,12 +1554,12 @@ namespace detail {
std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" };
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
#pragma omp parallel for schedule(dynamic)
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const WellControls* wc = stdWells().wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
// is instead treated as a default.
@ -1554,7 +1578,7 @@ namespace detail {
}
if (detail::constraintBroken(
xw.bhp(), xw.thp(), xw.wellRates(),
w, np, wells().type[w], wc, ctrl_index)) {
w, np, stdWells().wells().type[w], wc, ctrl_index)) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
}
@ -1563,7 +1587,7 @@ namespace detail {
// Constraint number ctrl_index was broken, switch to it.
if (terminal_output_)
{
std::cout << "Switching control mode for well " << wells().name[w]
std::cout << "Switching control mode for well " << stdWells().wells().name[w]
<< " from " << modestring[well_controls_iget_type(wc, current)]
<< " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl;
}
@ -1603,19 +1627,19 @@ namespace detail {
const double& alq = well_controls_iget_alq(wc, current);
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
const WellType& well_type = stdWells().wells().type[w];
if (well_type == INJECTOR) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
well_perforation_densities_, gravity);
stdWells().wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
stdWells().well_perforation_densities_, gravity);
xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
}
else if (well_type == PRODUCER) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
well_perforation_densities_, gravity);
stdWells().wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
std_wells_.well_perforation_densities_, gravity);
xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
}
@ -1681,7 +1705,7 @@ namespace detail {
WellState& well_state)
{
V aliveWells;
const int np = wells().number_of_phases;
const int np = stdWells().wells().number_of_phases;
std::vector<ADB> cq_s(np, ADB::null());
std::vector<int> indices = variableWellStateIndices();
SolutionState state0 = state;
@ -1746,7 +1770,7 @@ namespace detail {
if ( terminal_output_ ) {
std::cout << "well converged iter: " << it << std::endl;
}
const int nw = wells().number_of_wells;
const int nw = stdWells().wells().number_of_wells;
{
// We will set the bhp primary variable to the new ones,
// but we do not change the derivatives here.
@ -1786,8 +1810,8 @@ namespace detail {
{
if( ! localWellsActive() ) return;
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
ADB aqua = ADB::constant(ADB::V::Zero(nw));
ADB liquid = ADB::constant(ADB::V::Zero(nw));
@ -1828,7 +1852,7 @@ namespace detail {
//Run through all wells to calculate BHP/RATE targets
//and gather info about current control
for (int w = 0; w < nw; ++w) {
auto wc = wells().ctrls[w];
auto wc = stdWells().wells().ctrls[w];
// The current control in the well state overrides
// the current control set in the Wells struct, which
@ -1846,13 +1870,13 @@ namespace detail {
case THP:
{
const int perf = wells().well_connpos[w];
rho_v[w] = well_perforation_densities_[perf];
const int perf = stdWells().wells().well_connpos[w];
rho_v[w] = std_wells_.well_perforation_densities_[perf];
const int table_id = well_controls_iget_vfp(wc, current);
const double target = well_controls_iget_target(wc, current);
const WellType& well_type = wells().type[w];
const WellType& well_type = stdWells().wells().type[w];
if (well_type == INJECTOR) {
inj_table_id[w] = table_id;
thp_inj_target_v[w] = target;
@ -1910,7 +1934,7 @@ namespace detail {
//Perform hydrostatic correction to computed targets
double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, well_perforation_densities_, gravity);
const ADB::V dp_v = detail::computeHydrostaticCorrection(stdWells().wells(), vfp_ref_depth_v, std_wells_.well_perforation_densities_, gravity);
const ADB dp = ADB::constant(dp_v);
const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw);
const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw);
@ -2286,8 +2310,8 @@ namespace detail {
if( localWellsActive() )
{
const int np = wells().number_of_phases;
const int nw = wells().number_of_wells;
const int np = stdWells().wells().number_of_phases;
const int nw = stdWells().wells().number_of_wells;
// Extract parts of dwells corresponding to each part.
int varstart = 0;
@ -2323,7 +2347,7 @@ namespace detail {
//Loop over all wells
#pragma omp parallel for schedule(static)
for (int w=0; w<nw; ++w) {
const WellControls* wc = wells().ctrls[w];
const WellControls* wc = stdWells().wells().ctrls[w];
const int nwc = well_controls_get_num(wc);
//Loop over all controls until we find a THP control
//that specifies what we need...
@ -2347,18 +2371,18 @@ namespace detail {
double alq = well_controls_iget_alq(wc, ctrl_index);
int table_id = well_controls_iget_vfp(wc, ctrl_index);
const WellType& well_type = wells().type[w];
const WellType& well_type = stdWells().wells().type[w];
if (well_type == INJECTOR) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(),
well_perforation_densities_, gravity);
stdWells().wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(),
std_wells_.well_perforation_densities_, gravity);
well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
}
else if (well_type == PRODUCER) {
double dp = detail::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(),
well_perforation_densities_, gravity);
stdWells().wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(),
std_wells_.well_perforation_densities_, gravity);
well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
}