refactoring the interface of computeWellPotentials()

to reduce the cost of makeConstantState when not calculating the
potentials.
This commit is contained in:
Kai Bao 2016-05-11 13:13:08 +02:00
parent 19a256dce0
commit 846ff890de
3 changed files with 85 additions and 97 deletions

View File

@ -751,8 +751,6 @@ namespace detail {
{ {
using namespace Opm::AutoDiffGrid; using namespace Opm::AutoDiffGrid;
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
// If we have VFP tables, we need the well connection // If we have VFP tables, we need the well connection
// pressures for the "simple" hydrostatic correction // pressures for the "simple" hydrostatic correction
// between well depth and vfp table depth. // between well depth and vfp table depth.
@ -812,11 +810,11 @@ namespace detail {
asImpl().wellModel().addWellFluxEq(cq_s, state, residual_); asImpl().wellModel().addWellFluxEq(cq_s, state, residual_);
asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_); asImpl().wellModel().addWellControlEq(state, well_state, aliveWells, residual_);
{
if (param_.compute_well_potentials_) {
SolutionState state0 = state; SolutionState state0 = state;
asImpl().makeConstantState(state0); asImpl().makeConstantState(state0);
asImpl().wellModel().computeWellPotentials(state0, mob_perfcells, b_perfcells, vfp_properties_, asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state);
param_.compute_well_potentials_, gravity, well_state);
} }
} }

View File

@ -138,12 +138,9 @@ namespace Opm {
// state0 is non-constant, while it will not be used outside of the function // state0 is non-constant, while it will not be used outside of the function
template <class SolutionState, class WellState> template <class SolutionState, class WellState>
void void
computeWellPotentials(SolutionState& state0, computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells, const std::vector<ADB>& b_perfcells,
const VFPProperties& vfp_properties, SolutionState& state0,
const bool compute_well_potentials,
const bool gravity,
WellState& well_state); WellState& well_state);

View File

@ -1030,16 +1030,11 @@ namespace Opm
template <class SolutionState, class WellState> template <class SolutionState, class WellState>
void void
StandardWells::computeWellPotentials(SolutionState& state0, StandardWells::computeWellPotentials(const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& mob_perfcells,
const std::vector<ADB>& b_perfcells, const std::vector<ADB>& b_perfcells,
const VFPProperties& vfp_properties, SolutionState& state0,
const bool compute_well_potentials,
const bool gravity,
WellState& well_state) WellState& well_state)
{ {
//only compute well potentials if they are needed
if (compute_well_potentials) {
const int nw = wells().number_of_wells; const int nw = wells().number_of_wells;
const int np = wells().number_of_phases; const int np = wells().number_of_phases;
const Opm::PhaseUsage& pu = fluid_->phaseUsage(); const Opm::PhaseUsage& pu = fluid_->phaseUsage();
@ -1081,9 +1076,9 @@ namespace Opm
if (well_type == INJECTOR) { if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity); wellPerforationDensities(), gravity_);
const double bhp = vfp_properties.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictest of the bhp controlls i.e. smallest bhp for injectors // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) { if ( bhp < bhps[w]) {
bhps[w] = bhp; bhps[w] = bhp;
@ -1091,10 +1086,10 @@ namespace Opm
} }
else if (well_type == PRODUCER) { else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection( double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(), wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
wellPerforationDensities(), gravity); wellPerforationDensities(), gravity_);
const double bhp = vfp_properties.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictest of the bhp controlls i.e. largest bhp for producers // apply the strictest of the bhp controlls i.e. largest bhp for producers
if ( bhp > bhps[w]) { if ( bhp > bhps[w]) {
bhps[w] = bhp; bhps[w] = bhp;
@ -1119,15 +1114,13 @@ namespace Opm
// store well potentials in the well state // store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase // transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw]; const int nperf = wells().well_connpos[nw];
V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np); Vector cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) { for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np); cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
} }
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
} }
}
@ -1177,12 +1170,12 @@ namespace Opm
// The transpose() below switches the ordering. // The transpose() below switches the ordering.
const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose(); const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
const V qs = Eigen::Map<const V>(wrates.data(), nw*np); const Vector qs = Eigen::Map<const V>(wrates.data(), nw*np);
vars0.push_back(qs); vars0.push_back(qs);
// Initial well bottom-hole pressure. // Initial well bottom-hole pressure.
assert (not xw.bhp().empty()); assert (not xw.bhp().empty());
const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size()); const Vector bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
vars0.push_back(bhp); vars0.push_back(bhp);
} }
else else