diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index db00f77ce..3952697a2 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -751,8 +751,6 @@ namespace detail { { using namespace Opm::AutoDiffGrid; - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - // If we have VFP tables, we need the well connection // pressures for the "simple" hydrostatic correction // between well depth and vfp table depth. @@ -812,11 +810,11 @@ namespace detail { asImpl().wellModel().addWellFluxEq(cq_s, state, residual_); 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(state0, mob_perfcells, b_perfcells, vfp_properties_, - param_.compute_well_potentials_, gravity, well_state); + asImpl().wellModel().computeWellPotentials(mob_perfcells, b_perfcells, state0, well_state); } } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index e6bc74a9a..0968d4b2d 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -138,12 +138,9 @@ namespace Opm { // state0 is non-constant, while it will not be used outside of the function template void - computeWellPotentials(SolutionState& state0, - const std::vector& mob_perfcells, + computeWellPotentials(const std::vector& mob_perfcells, const std::vector& b_perfcells, - const VFPProperties& vfp_properties, - const bool compute_well_potentials, - const bool gravity, + SolutionState& state0, WellState& well_state); diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 3e8ba73a3..ce0f33f33 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -1030,102 +1030,95 @@ namespace Opm template void - StandardWells::computeWellPotentials(SolutionState& state0, - const std::vector& mob_perfcells, + StandardWells::computeWellPotentials(const std::vector& mob_perfcells, const std::vector& b_perfcells, - const VFPProperties& vfp_properties, - const bool compute_well_potentials, - const bool gravity, + SolutionState& state0, WellState& well_state) { - //only compute well potentials if they are needed - if (compute_well_potentials) { - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - const Opm::PhaseUsage& pu = fluid_->phaseUsage(); + const int nw = wells().number_of_wells; + const int np = wells().number_of_phases; + const Opm::PhaseUsage& pu = fluid_->phaseUsage(); - Vector bhps = Vector::Zero(nw); - for (int w = 0; w < nw; ++w) { - const WellControls* ctrl = wells().ctrls[w]; - const int nwc = well_controls_get_num(ctrl); - //Loop over all controls until we find a BHP control - //or a THP control that specifies what we need. - //Pick the value that gives the most restrictive flow - for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + Vector bhps = Vector::Zero(nw); + for (int w = 0; w < nw; ++w) { + const WellControls* ctrl = wells().ctrls[w]; + const int nwc = well_controls_get_num(ctrl); + //Loop over all controls until we find a BHP control + //or a THP control that specifies what we need. + //Pick the value that gives the most restrictive flow + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { - bhps[w] = well_controls_iget_target(ctrl, ctrl_index); - } - - if (well_controls_iget_type(ctrl, ctrl_index) == THP) { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - 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(ctrl, ctrl_index); - const double& thp = well_controls_iget_target(ctrl, ctrl_index); - const double& alq = well_controls_iget_alq(ctrl, ctrl_index); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity); - 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 - if ( bhp < bhps[w]) { - bhps[w] = bhp; - } - } - else if (well_type == PRODUCER) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(), - wellPerforationDensities(), gravity); - - 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 - if ( bhp > bhps[w]) { - bhps[w] = bhp; - } - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - } + if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { + bhps[w] = well_controls_iget_target(ctrl, ctrl_index); } + if (well_controls_iget_type(ctrl, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + 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(ctrl, ctrl_index); + const double& thp = well_controls_iget_target(ctrl, ctrl_index); + const double& alq = well_controls_iget_alq(ctrl, ctrl_index); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + 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 + if ( bhp < bhps[w]) { + bhps[w] = bhp; + } + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity_); + + 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 + if ( bhp > bhps[w]) { + bhps[w] = bhp; + } + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + } } - // use bhp limit from control - state0.bhp = ADB::constant(bhps); - - // compute well potentials - Vector aliveWells; - std::vector well_potentials; - computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_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]; - V 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_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); } + // use bhp limit from control + state0.bhp = ADB::constant(bhps); + + // compute well potentials + Vector aliveWells; + std::vector well_potentials; + computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, well_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_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); } @@ -1177,12 +1170,12 @@ namespace Opm // The transpose() below switches the ordering. const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); - const V qs = Eigen::Map(wrates.data(), nw*np); + const Vector qs = Eigen::Map(wrates.data(), nw*np); vars0.push_back(qs); // Initial well bottom-hole pressure. assert (not xw.bhp().empty()); - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + const Vector bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); vars0.push_back(bhp); } else