small improvements to StandardWellsDense

the most significant change is that only a "PhaseUsage" object must be
passed to its constructor instead of a full "fluid object". also, the
pointers to the vector of "active" phases has been converted into a
full object. (this helps to avoid potential use-after-free errors.)
This commit is contained in:
Andreas Lauser
2017-01-11 17:02:20 +01:00
parent 4411712dce
commit b324d17003
2 changed files with 56 additions and 60 deletions

View File

@@ -195,7 +195,7 @@ namespace Opm {
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size()); const std::vector<double> pv(geo_.poreVolume().data(), geo_.poreVolume().data() + geo_.poreVolume().size());
const std::vector<double> depth(geo_.z().data(), geo_.z().data() + geo_.z().size()); const std::vector<double> depth(geo_.z().data(), geo_.z().data() + geo_.z().size());
well_model_.init(&fluid_, &active_, &vfp_properties_, gravity, depth, pv); well_model_.init(fluid_.phaseUsage(), active_, &vfp_properties_, gravity, depth, pv);
wellModel().setWellsActive( localWellsActive() ); wellModel().setWellsActive( localWellsActive() );
global_nc_ = Opm::AutoDiffGrid::numCells(grid_); global_nc_ = Opm::AutoDiffGrid::numCells(grid_);
// compute global sum of number of cells // compute global sum of number of cells

View File

@@ -44,7 +44,6 @@
#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp> #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
#include <opm/autodiff/BlackoilDetails.hpp> #include <opm/autodiff/BlackoilDetails.hpp>
#include <opm/autodiff/BlackoilModelParameters.hpp> #include <opm/autodiff/BlackoilModelParameters.hpp>
#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp> #include <opm/autodiff/WellStateFullyImplicitBlackoilDense.hpp>
#include<dune/common/fmatrix.hh> #include<dune/common/fmatrix.hh>
#include<dune/istl/bcrsmatrix.hh> #include<dune/istl/bcrsmatrix.hh>
@@ -88,9 +87,6 @@ enum WellVariablePositions {
, wells_(wells_arg) , wells_(wells_arg)
, param_(param) , param_(param)
, terminal_output_(terminal_output) , terminal_output_(terminal_output)
, fluid_(nullptr)
, active_(nullptr)
, vfp_properties_(nullptr)
, well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) , well_perforation_densities_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0)
, well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0) , well_perforation_pressure_diffs_( wells_ ? wells_arg->well_connpos[wells_arg->number_of_wells] : 0)
, wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0) , wellVariables_( wells_ ? (wells_arg->number_of_wells * wells_arg->number_of_phases) : 0)
@@ -104,8 +100,8 @@ enum WellVariablePositions {
} }
} }
void init(const BlackoilPropsAdFromDeck* fluid_arg, void init(const PhaseUsage phase_usage_arg,
const std::vector<bool>* active_arg, const std::vector<bool>& active_arg,
const VFPProperties* vfp_properties_arg, const VFPProperties* vfp_properties_arg,
const double gravity_arg, const double gravity_arg,
const std::vector<double>& depth_arg, const std::vector<double>& depth_arg,
@@ -116,7 +112,7 @@ enum WellVariablePositions {
return; return;
} }
fluid_ = fluid_arg; phase_usage_ = phase_usage_arg;
active_ = active_arg; active_ = active_arg;
vfp_properties_ = vfp_properties_arg; vfp_properties_ = vfp_properties_arg;
gravity_ = gravity_arg; gravity_ = gravity_arg;
@@ -132,7 +128,7 @@ enum WellVariablePositions {
const int nc = numCells(); const int nc = numCells();
#ifndef NDEBUG #ifndef NDEBUG
const auto pu = fluid_->phaseUsage(); const auto pu = phase_usage_;
const int np = pu.num_phases; const int np = pu.num_phases;
// assumes the gas fractions are stored after water fractions // assumes the gas fractions are stored after water fractions
@@ -263,7 +259,7 @@ enum WellVariablePositions {
// add trivial equation for 2p cases (Only support water + oil) // add trivial equation for 2p cases (Only support water + oil)
if (np == 2) { if (np == 2) {
assert(!(*active_)[ Gas ]); assert(!active_[ Gas ]);
invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0; invDuneD_[w][w][flowPhaseToEbosCompIdx(Gas)][flowToEbosPvIdx(Gas)] = 1.0;
} }
@@ -537,7 +533,7 @@ enum WellVariablePositions {
void void
computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const computeWellFlux(const int& w, const double& Tw, const intensiveQuants& intQuants, const double& cdp, const bool& allow_cf, std::vector<EvalWell>& cq_s) const
{ {
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
EvalWell bhp = getBhp(w); EvalWell bhp = getBhp(w);
const int np = wells().number_of_phases; const int np = wells().number_of_phases;
std::vector<EvalWell> cmix_s(np,0.0); std::vector<EvalWell> cmix_s(np,0.0);
@@ -574,7 +570,7 @@ enum WellVariablePositions {
cq_ps[phase] = b_perfcells_dense[phase] * cq_p; cq_ps[phase] = b_perfcells_dense[phase] * cq_p;
} }
if ((*active_)[Oil] && (*active_)[Gas]) { if (active_[Oil] && active_[Gas]) {
const int oilpos = pu.phase_pos[Oil]; const int oilpos = pu.phase_pos[Oil];
const int gaspos = pu.phase_pos[Gas]; const int gaspos = pu.phase_pos[Gas];
const EvalWell cq_psOil = cq_ps[oilpos]; const EvalWell cq_psOil = cq_ps[oilpos];
@@ -603,12 +599,12 @@ enum WellVariablePositions {
// compute volume ratio between connection at standard conditions // compute volume ratio between connection at standard conditions
EvalWell volumeRatio = 0.0; EvalWell volumeRatio = 0.0;
if ((*active_)[Water]) { if (active_[Water]) {
const int watpos = pu.phase_pos[Water]; const int watpos = pu.phase_pos[Water];
volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos]; volumeRatio += cmix_s[watpos] / b_perfcells_dense[watpos];
} }
if ((*active_)[Oil] && (*active_)[Gas]) { if (active_[Oil] && active_[Gas]) {
EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx)); EvalWell well_temperature = extendEval(fs.temperature(FluidSystem::oilPhaseIdx));
EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); EvalWell rsSatEval = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), well_temperature, well_pressure);
EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure); EvalWell rvSatEval = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), well_temperature, well_pressure);
@@ -645,11 +641,11 @@ enum WellVariablePositions {
volumeRatio += tmp_gas / b_perfcells_dense[gaspos]; volumeRatio += tmp_gas / b_perfcells_dense[gaspos];
} }
else { else {
if ((*active_)[Oil]) { if (active_[Oil]) {
const int oilpos = pu.phase_pos[Oil]; const int oilpos = pu.phase_pos[Oil];
volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos]; volumeRatio += cmix_s[oilpos] / b_perfcells_dense[oilpos];
} }
if ((*active_)[Gas]) { if (active_[Gas]) {
const int gaspos = pu.phase_pos[Gas]; const int gaspos = pu.phase_pos[Gas];
volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos]; volumeRatio += cmix_s[gaspos] / b_perfcells_dense[gaspos];
} }
@@ -846,8 +842,8 @@ enum WellVariablePositions {
{ {
const int nperf = wells().well_connpos[wells().number_of_wells]; const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells; const int nw = wells().number_of_wells;
const PhaseUsage& pu = fluid_->phaseUsage(); const PhaseUsage& pu = phase_usage_;
const int np = fluid_->numPhases(); const int np = phase_usage_.num_phases;
b_perf.resize(nperf*np); b_perf.resize(nperf*np);
surf_dens_perf.resize(nperf*np); surf_dens_perf.resize(nperf*np);
@@ -949,42 +945,42 @@ enum WellVariablePositions {
// update the second and third well variable (The flux fractions) // update the second and third well variable (The flux fractions)
std::vector<double> F(np,0.0); std::vector<double> F(np,0.0);
if ((*active_)[ Water ]) { if (active_[ Water ]) {
const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1; const int sign2 = dwells[w][flowPhaseToEbosCompIdx(WFrac)] > 0 ? 1: -1;
const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit); const double dx2_limited = sign2 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(WFrac)]),dFLimit);
well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited; well_state.wellSolutions()[WFrac*nw + w] = xvar_well_old[WFrac*nw + w] - dx2_limited;
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1; const int sign3 = dwells[w][flowPhaseToEbosCompIdx(GFrac)] > 0 ? 1: -1;
const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit); const double dx3_limited = sign3 * std::min(std::abs(dwells[w][flowPhaseToEbosCompIdx(GFrac)]),dFLimit);
well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited; well_state.wellSolutions()[GFrac*nw + w] = xvar_well_old[GFrac*nw + w] - dx3_limited;
} }
assert((*active_)[ Oil ]); assert(active_[ Oil ]);
F[Oil] = 1.0; F[Oil] = 1.0;
if ((*active_)[ Water ]) { if (active_[ Water ]) {
F[Water] = well_state.wellSolutions()[WFrac*nw + w]; F[Water] = well_state.wellSolutions()[WFrac*nw + w];
F[Oil] -= F[Water]; F[Oil] -= F[Water];
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
F[Gas] = well_state.wellSolutions()[GFrac*nw + w]; F[Gas] = well_state.wellSolutions()[GFrac*nw + w];
F[Oil] -= F[Gas]; F[Oil] -= F[Gas];
} }
if ((*active_)[ Water ]) { if (active_[ Water ]) {
if (F[Water] < 0.0) { if (F[Water] < 0.0) {
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
F[Gas] /= (1.0 - F[Water]); F[Gas] /= (1.0 - F[Water]);
} }
F[Oil] /= (1.0 - F[Water]); F[Oil] /= (1.0 - F[Water]);
F[Water] = 0.0; F[Water] = 0.0;
} }
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
if (F[Gas] < 0.0) { if (F[Gas] < 0.0) {
if ((*active_)[ Water ]) { if (active_[ Water ]) {
F[Water] /= (1.0 - F[Gas]); F[Water] /= (1.0 - F[Gas]);
} }
F[Oil] /= (1.0 - F[Gas]); F[Oil] /= (1.0 - F[Gas]);
@@ -992,19 +988,19 @@ enum WellVariablePositions {
} }
} }
if (F[Oil] < 0.0) { if (F[Oil] < 0.0) {
if ((*active_)[ Water ]) { if (active_[ Water ]) {
F[Water] /= (1.0 - F[Oil]); F[Water] /= (1.0 - F[Oil]);
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
F[Gas] /= (1.0 - F[Oil]); F[Gas] /= (1.0 - F[Oil]);
} }
F[Oil] = 0.0; F[Oil] = 0.0;
} }
if ((*active_)[ Water ]) { if (active_[ Water ]) {
well_state.wellSolutions()[WFrac*nw + w] = F[Water]; well_state.wellSolutions()[WFrac*nw + w] = F[Water];
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
well_state.wellSolutions()[GFrac*nw + w] = F[Gas]; well_state.wellSolutions()[GFrac*nw + w] = F[Gas];
} }
@@ -1056,15 +1052,15 @@ enum WellVariablePositions {
double liquid = 0.0; double liquid = 0.0;
double vapour = 0.0; double vapour = 0.0;
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
if ((*active_)[ Water ]) { if (active_[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ]; aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
} }
if ((*active_)[ Oil ]) { if (active_[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ]; liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ]; vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
} }
@@ -1203,15 +1199,15 @@ enum WellVariablePositions {
double liquid = 0.0; double liquid = 0.0;
double vapour = 0.0; double vapour = 0.0;
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
if ((*active_)[ Water ]) { if (active_[ Water ]) {
aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ];
} }
if ((*active_)[ Oil ]) { if (active_[ Oil ]) {
liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ];
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ];
} }
@@ -1325,18 +1321,18 @@ enum WellVariablePositions {
tot_well_rate += g[p] * xw.wellRates()[np*w + p]; tot_well_rate += g[p] * xw.wellRates()[np*w + p];
} }
if(std::abs(tot_well_rate) > 0) { if(std::abs(tot_well_rate) > 0) {
if ((*active_)[ Water ]) { if (active_[ Water ]) {
xw.wellSolutions()[WFrac*nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate; xw.wellSolutions()[WFrac*nw + w] = g[Water] * xw.wellRates()[np*w + Water] / tot_well_rate;
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
xw.wellSolutions()[GFrac*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ; xw.wellSolutions()[GFrac*nw + w] = g[Gas] * xw.wellRates()[np*w + Gas] / tot_well_rate ;
} }
} else { } else {
if ((*active_)[ Water ]) { if (active_[ Water ]) {
xw.wellSolutions()[WFrac*nw + w] = wells().comp_frac[np*w + Water]; xw.wellSolutions()[WFrac*nw + w] = wells().comp_frac[np*w + Water];
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
xw.wellSolutions()[GFrac*nw + w] = wells().comp_frac[np*w + Gas]; xw.wellSolutions()[GFrac*nw + w] = wells().comp_frac[np*w + Gas];
} }
} }
@@ -1466,7 +1462,7 @@ enum WellVariablePositions {
// Compute densities // Compute densities
well_perforation_densities_ = well_perforation_densities_ =
WellDensitySegmented::computeConnectionDensities( WellDensitySegmented::computeConnectionDensities(
wells(), xw, fluid_->phaseUsage(), wells(), xw, phase_usage_,
b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
// Compute pressure deltas // Compute pressure deltas
@@ -1484,8 +1480,8 @@ enum WellVariablePositions {
ModelParameters param_; ModelParameters param_;
bool terminal_output_; bool terminal_output_;
const BlackoilPropsAdFromDeck* fluid_; PhaseUsage phase_usage_;
const std::vector<bool>* active_; std::vector<bool> active_;
const VFPProperties* vfp_properties_; const VFPProperties* vfp_properties_;
double gravity_; double gravity_;
// the depth of the all the cell centers // the depth of the all the cell centers
@@ -1531,15 +1527,15 @@ enum WellVariablePositions {
EvalWell bhp = 0.0; EvalWell bhp = 0.0;
double vfp_ref_depth = 0.0; double vfp_ref_depth = 0.0;
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
if ((*active_)[ Water ]) { if (active_[ Water ]) {
aqua = getQs(wellIdx, pu.phase_pos[ Water]); aqua = getQs(wellIdx, pu.phase_pos[ Water]);
} }
if ((*active_)[ Oil ]) { if (active_[ Oil ]) {
liquid = getQs(wellIdx, pu.phase_pos[ Oil ]); liquid = getQs(wellIdx, pu.phase_pos[ Oil ]);
} }
if ((*active_)[ Gas ]) { if (active_[ Gas ]) {
vapour = getQs(wellIdx, pu.phase_pos[ Gas ]); vapour = getQs(wellIdx, pu.phase_pos[ Gas ]);
} }
if (wells().type[wellIdx] == INJECTOR) { if (wells().type[wellIdx] == INJECTOR) {
@@ -1619,11 +1615,11 @@ enum WellVariablePositions {
// Oil fraction // Oil fraction
EvalWell well_fraction = 1.0; EvalWell well_fraction = 1.0;
if ((*active_)[Water]) { if (active_[Water]) {
well_fraction -= wellVariables_[WFrac * nw + wellIdx]; well_fraction -= wellVariables_[WFrac * nw + wellIdx];
} }
if ((*active_)[Gas]) { if (active_[Gas]) {
well_fraction -= wellVariables_[GFrac * nw + wellIdx]; well_fraction -= wellVariables_[GFrac * nw + wellIdx];
} }
return well_fraction; return well_fraction;
@@ -1646,11 +1642,11 @@ enum WellVariablePositions {
const WellState& well_state, const WellState& well_state,
const int well_number) const const int well_number) const
{ {
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
const int np = well_state.numPhases(); const int np = well_state.numPhases();
if (econ_production_limits.onMinOilRate()) { if (econ_production_limits.onMinOilRate()) {
assert((*active_)[Oil]); assert(active_[Oil]);
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double min_oil_rate = econ_production_limits.minOilRate(); const double min_oil_rate = econ_production_limits.minOilRate();
if (std::abs(oil_rate) < min_oil_rate) { if (std::abs(oil_rate) < min_oil_rate) {
@@ -1659,7 +1655,7 @@ enum WellVariablePositions {
} }
if (econ_production_limits.onMinGasRate() ) { if (econ_production_limits.onMinGasRate() ) {
assert((*active_)[Gas]); assert(active_[Gas]);
const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ]; const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ];
const double min_gas_rate = econ_production_limits.minGasRate(); const double min_gas_rate = econ_production_limits.minGasRate();
if (std::abs(gas_rate) < min_gas_rate) { if (std::abs(gas_rate) < min_gas_rate) {
@@ -1668,8 +1664,8 @@ enum WellVariablePositions {
} }
if (econ_production_limits.onMinLiquidRate() ) { if (econ_production_limits.onMinLiquidRate() ) {
assert((*active_)[Oil]); assert(active_[Oil]);
assert((*active_)[Water]); assert(active_[Water]);
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
const double liquid_rate = oil_rate + water_rate; const double liquid_rate = oil_rate + water_rate;
@@ -1767,11 +1763,11 @@ enum WellVariablePositions {
double violation_extent = -1.0; double violation_extent = -1.0;
const int np = well_state.numPhases(); const int np = well_state.numPhases();
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = phase_usage_;
const int well_number = map_entry[0]; const int well_number = map_entry[0];
assert((*active_)[Oil]); assert(active_[Oil]);
assert((*active_)[Water]); assert(active_[Water]);
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ]; const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ]; const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];