diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index d39e4e719..f5341eae5 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -356,16 +356,10 @@ namespace Opm { void variableReservoirStateInitials(const ReservoirState& x, std::vector& vars0) const; - void - variableWellStateInitials(const WellState& xw, - std::vector& vars0) const; std::vector variableStateIndices() const; - std::vector - variableWellStateIndices() const; - SolutionState variableStateExtractVars(const ReservoirState& x, const std::vector& indices, @@ -380,9 +374,6 @@ namespace Opm { computeAccum(const SolutionState& state, const int aix ); - void computeWellConnectionPressures(const SolutionState& state, - const WellState& xw); - void assembleMassBalanceEq(const SolutionState& state); @@ -405,27 +396,11 @@ namespace Opm { SolutionState& state, WellState& well_state); - void - computeWellPotentials(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - WellState& well_state); - - - void - addWellFluxEq(const std::vector& cq_s, - const SolutionState& state); - void addWellContributionToMassBalanceEq(const std::vector& cq_s, const SolutionState& state, const WellState& xw); - void - addWellControlEq(const SolutionState& state, - const WellState& xw, - const V& aliveWells); - bool getWellConvergence(const int iteration); bool isVFPActive() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 32d7bee7b..6ecfa2576 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -507,7 +507,7 @@ namespace detail { // and bhp and Q for the wells vars0.reserve(np + 1); variableReservoirStateInitials(x, vars0); - asImpl().variableWellStateInitials(xw, vars0); + asImpl().stdWells().variableWellStateInitials(xw, vars0); return vars0; } @@ -554,41 +554,6 @@ namespace detail { - template - void - BlackoilModelBase:: - variableWellStateInitials(const WellState& xw, std::vector& vars0) const - { - // Initial well rates. - if ( stdWells().localWellsActive() ) - { - // 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; - - // 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); - 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()); - vars0.push_back(bhp); - } - else - { - // push null states for qs and bhp - vars0.push_back(V()); - vars0.push_back(V()); - } - } - - - - - template std::vector BlackoilModelBase:: @@ -604,8 +569,7 @@ namespace detail { if (active_[Gas]) { indices[Xvar] = next++; } - indices[Qs] = next++; - indices[Bhp] = next++; + asImpl().stdWells().variableStateWellIndices(indices, next); assert(next == fluid_.numPhases() + 2); return indices; } @@ -613,24 +577,6 @@ namespace detail { - template - std::vector - BlackoilModelBase:: - variableWellStateIndices() const - { - // Black oil model standard is 5 equation. - // For the pure well solve, only the well equations are picked. - std::vector indices(5, -1); - int next = 0; - indices[Qs] = next++; - indices[Bhp] = next++; - assert(next == 2); - return indices; - } - - - - template typename BlackoilModelBase::SolutionState @@ -798,41 +744,6 @@ namespace detail { - template - void - BlackoilModelBase:: - computeWellConnectionPressures(const SolutionState& state, - const WellState& xw) - { - if( ! localWellsActive() ) return ; - - using namespace Opm::AutoDiffGrid; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. - std::vector b_perf; - std::vector rsmax_perf; - std::vector rvmax_perf; - std::vector surf_dens_perf; - asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - // Extract well connection depths. - const typename WellModel::Vector depth = cellCentroidsZToEigen(grid_); - const typename WellModel::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); - - // Gravity - const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - - asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav); - - } - - - - - template void BlackoilModelBase:: @@ -842,6 +753,9 @@ namespace detail { { using namespace Opm::AutoDiffGrid; + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = cellCentroidsZToEigen(grid_); + // If we have VFP tables, we need the well connection // pressures for the "simple" hydrostatic correction // between well depth and vfp table depth. @@ -849,14 +763,15 @@ namespace detail { SolutionState state = asImpl().variableState(reservoir_state, well_state); SolutionState state0 = state; asImpl().makeConstantState(state0); - asImpl().computeWellConnectionPressures(state0, well_state); + // asImpl().computeWellConnectionPressures(state0, well_state); + // Extract well connection depths. + asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells // asImpl().updateWellControls(well_state); // asImpl().stdWells().updateWellControls(well_state); - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); // Create the primary variables. @@ -869,7 +784,8 @@ namespace detail { // Compute initial accumulation contributions // and well connection pressures. asImpl().computeAccum(state0, 0); - asImpl().computeWellConnectionPressures(state0, well_state); + // asImpl().computeWellConnectionPressures(state0, well_state); + asImpl().stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // OPM_AD_DISKVAL(state.pressure); @@ -901,10 +817,17 @@ namespace detail { std::vector cq_s; asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); - asImpl().addWellFluxEq(cq_s, state); + asImpl().stdWells().addWellFluxEq(cq_s, state, residual_); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); - asImpl().addWellControlEq(state, well_state, aliveWells); - asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state); + asImpl().stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); + // asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state); + { + SolutionState state0 = state; + asImpl().makeConstantState(state0); + asImpl().stdWells().computeWellPotentials(state0, mob_perfcells, b_perfcells, + fluid_.phaseUsage(), active_, vfp_properties_, + param_.compute_well_potentials_, gravity, well_state); + } } @@ -1069,34 +992,6 @@ namespace detail { - template - void - BlackoilModelBase:: - addWellFluxEq(const std::vector& cq_s, - const SolutionState& state) - { - if( !asImpl().localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - ADB qs = state.qs; - for (int phase = 0; phase < np; ++phase) { - qs -= superset(stdWells().wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); - - } - - residual_.well_flux_eq = qs; - } - - - - - template bool BlackoilModelBase:: @@ -1144,7 +1039,7 @@ namespace detail { V aliveWells; const int np = wells().number_of_phases; std::vector cq_s(np, ADB::null()); - std::vector indices = variableWellStateIndices(); + std::vector indices = asImpl().stdWells().variableWellStateIndices(); SolutionState state0 = state; WellState well_state0 = well_state; asImpl().makeConstantState(state0); @@ -1161,21 +1056,25 @@ namespace detail { } } + // gravity + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + int it = 0; bool converged; do { // bhp and Q for the wells std::vector vars0; vars0.reserve(2); - asImpl().variableWellStateInitials(well_state, vars0); + asImpl().stdWells().variableWellStateInitials(well_state, vars0); std::vector vars = ADB::variables(vars0); SolutionState wellSolutionState = state0; asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); - asImpl().addWellFluxEq(cq_s, wellSolutionState); - asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); + asImpl().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_); + asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells, + active_, vfp_properties_, gravity, residual_); converged = getWellConvergence(it); if (converged) { @@ -1199,7 +1098,6 @@ namespace detail { const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); assert(dx.size() == total_residual_v.size()); // asImpl().updateWellState(dx.array(), well_state); - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); asImpl().stdWells().updateWellState(dx.array(), gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); asImpl().stdWells(). updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); } @@ -1228,7 +1126,9 @@ namespace detail { std::vector old_derivs = state.qs.derivative(); state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); } - asImpl().computeWellConnectionPressures(state, well_state); + // asImpl().computeWellConnectionPressures(state, well_state); + const ADB::V depth = Opm::AutoDiffGrid::cellCentroidsZToEigen(grid_); + asImpl().stdWells().computeWellConnectionPressures(state, well_state, fluid_, active_, phaseCondition(), depth, gravity); } if (!converged) { @@ -1242,174 +1142,6 @@ namespace detail { - template - void - BlackoilModelBase:: - addWellControlEq(const SolutionState& state, - const WellState& xw, - const V& aliveWells) - { - if( ! localWellsActive() ) return; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - ADB aqua = ADB::constant(ADB::V::Zero(nw)); - ADB liquid = ADB::constant(ADB::V::Zero(nw)); - ADB vapour = ADB::constant(ADB::V::Zero(nw)); - - if (active_[Water]) { - aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); - } - if (active_[Oil]) { - liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); - } - if (active_[Gas]) { - vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); - } - - //THP calculation variables - std::vector inj_table_id(nw, -1); - std::vector prod_table_id(nw, -1); - ADB::V thp_inj_target_v = ADB::V::Zero(nw); - ADB::V thp_prod_target_v = ADB::V::Zero(nw); - ADB::V alq_v = ADB::V::Zero(nw); - - //Hydrostatic correction variables - ADB::V rho_v = ADB::V::Zero(nw); - ADB::V vfp_ref_depth_v = ADB::V::Zero(nw); - - //Target vars - ADB::V bhp_targets = ADB::V::Zero(nw); - ADB::V rate_targets = ADB::V::Zero(nw); - Eigen::SparseMatrix rate_distr(nw, np*nw); - - //Selection variables - std::vector bhp_elems; - std::vector thp_inj_elems; - std::vector thp_prod_elems; - std::vector rate_elems; - - //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]; - - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - const int current = xw.currentControls()[w]; - - switch (well_controls_iget_type(wc, current)) { - case BHP: - { - bhp_elems.push_back(w); - bhp_targets(w) = well_controls_iget_target(wc, current); - rate_targets(w) = -1e100; - } - break; - - case THP: - { - const int perf = wells().well_connpos[w]; - rho_v[w] = asImpl().stdWells().wellPerforationDensities()[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]; - if (well_type == INJECTOR) { - inj_table_id[w] = table_id; - thp_inj_target_v[w] = target; - alq_v[w] = -1e100; - - vfp_ref_depth_v[w] = vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(); - - thp_inj_elems.push_back(w); - } - else if (well_type == PRODUCER) { - prod_table_id[w] = table_id; - thp_prod_target_v[w] = target; - alq_v[w] = well_controls_iget_alq(wc, current); - - vfp_ref_depth_v[w] = vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(); - - thp_prod_elems.push_back(w); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well"); - } - bhp_targets(w) = -1e100; - rate_targets(w) = -1e100; - } - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - { - rate_elems.push_back(w); - // 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]; - } - - bhp_targets(w) = -1.0e100; - rate_targets(w) = well_controls_iget_target(wc, current); - } - break; - } - } - - //Calculate BHP target from THP - const ADB thp_inj_target = ADB::constant(thp_inj_target_v); - const ADB thp_prod_target = ADB::constant(thp_prod_target_v); - const ADB alq = ADB::constant(alq_v); - const ADB bhp_from_thp_inj = vfp_properties_.getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target); - const ADB bhp_from_thp_prod = vfp_properties_.getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq); - - //Perform hydrostatic correction to computed targets - double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const ADB::V dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), 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); - - //Calculate residuals - const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj; - const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod; - const ADB bhp_residual = state.bhp - bhp_targets; - const ADB rate_residual = rate_distr * state.qs - rate_targets; - - //Select the right residual for each well - residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) + - superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) + - superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) + - superset(subset(rate_residual, rate_elems), rate_elems, nw); - - // For wells that are dead (not flowing), and therefore not communicating - // with the reservoir, we set the equation to be equal to the well's total - // flow. This will be a solution only if the target rate is also zero. - Eigen::SparseMatrix rate_summer(nw, np*nw); - for (int w = 0; w < nw; ++w) { - for (int phase = 0; phase < np; ++phase) { - rate_summer.insert(w, phase*nw + w) = 1.0; - } - } - Selector alive_selector(aliveWells, Selector::NotEqualZero); - residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs); - // OPM_AD_DUMP(residual_.well_eq); - } - - - - - template V BlackoilModelBase:: @@ -1767,111 +1499,6 @@ namespace detail { - template - void - BlackoilModelBase:: - computeWellPotentials(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - WellState& well_state) - { - //only compute well potentials if they are needed - if (param_.compute_well_potentials_) { - const int nw = wells().number_of_wells; - const int np = wells().number_of_phases; - const Opm::PhaseUsage pu = fluid_.phaseUsage(); - V bhps = V::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]; - - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - - if (well_type == INJECTOR) { - double dp = wellhelpers::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - stdWells().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(), - stdWells().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 - SolutionState state0 = state; - asImpl().makeConstantState(state0); - state0.bhp = ADB::constant(bhps); - - // compute well potentials - V aliveWells; - std::vector well_potentials; - asImpl().stdWells().computeWellFlux(state0, fluid_.phaseUsage(), active_, 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); - } - - } - - - - - template std::vector BlackoilModelBase:: diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index 3ff7094b2..4d659c920 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -223,8 +223,9 @@ namespace Opm { using Base::convergenceReduction; using Base::maxResidualAllowed; using Base::variableState; - using Base::variableWellStateIndices; + // using Base::variableWellStateIndices; using Base::asImpl; + using Base::variableReservoirStateInitials; const std::vector& wellsMultiSegment() const { return wells_multisegment_; } @@ -234,6 +235,10 @@ namespace Opm { void updateWellState(const V& dwells, WellState& well_state); + std::vector + variableStateInitials(const ReservoirState& x, + const WellState& xw) const; + void variableWellStateInitials(const WellState& xw, std::vector& vars0) const; diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index c87c2c7ac..3026df908 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -1559,7 +1559,7 @@ namespace Opm { V aliveWells; const int np = wells().number_of_phases; std::vector cq_s(np, ADB::null()); - std::vector indices = variableWellStateIndices(); + std::vector indices = stdWells().variableWellStateIndices(); SolutionState state0 = state; WellState well_state0 = well_state; makeConstantState(state0); @@ -1653,6 +1653,28 @@ namespace Opm { } + + + + template + std::vector + BlackoilMultiSegmentModel:: + variableStateInitials(const ReservoirState& x, + const WellState& xw) const + { + assert(active_[ Oil ]); + + const int np = x.numPhases(); + + std::vector vars0; + // p, Sw and Rs, Rv or Sg is used as primary depending on solution conditions + // and bhp and Q for the wells + vars0.reserve(np + 1); + variableReservoirStateInitials(x, vars0); + variableWellStateInitials(xw, vars0); + return vars0; + } + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index ddc3aa40c..d3b09322e 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -146,8 +146,8 @@ namespace Opm { using Base::drMaxRel; using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; - using Base::addWellControlEq; + // using Base::computeWellConnectionPressures; + // using Base::addWellControlEq; // using Base::computePropertiesForWellConnectionPressures; std::vector diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 33a15812f..9f12a9596 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -133,14 +133,65 @@ namespace Opm { const VFPProperties& vfp_properties, WellState& well_state); - template - void updateWellControls(const Opm::PhaseUsage& pu, - const double gravity, - const VFPProperties& vfp_properties, - const bool terminal_output, - const std::vector& active, - WellState& xw) const; + template + void updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const; + // TODO: should LinearisedBlackoilResidual also be a template class? + template + void addWellFluxEq(const std::vector& cq_s, + const SolutionState& state, + LinearisedBlackoilResidual& residual); + + // TODO: some parameters, like gravity, maybe it is better to put in the member list + template + void addWellControlEq(const SolutionState& state, + const WellState& xw, + const Vector& aliveWells, + const std::vector active, + const VFPProperties& vfp_properties, + const double gravity, + LinearisedBlackoilResidual& residual); + + template + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity); + + // state0 is non-constant, while it will not be used outside of the function + template + void + computeWellPotentials(SolutionState& state0, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + const Opm::PhaseUsage& pu, + const std::vector active, + const VFPProperties& vfp_properties, + const bool compute_well_potentials, + const bool gravity, + WellState& well_state); + + + void + variableStateWellIndices(std::vector& indices, + int& next) const; + + + std::vector + variableWellStateIndices() const; + + template + void + variableWellStateInitials(const WellState& xw, + std::vector& vars0) const; protected: bool wells_active_; diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp index cc491d323..5801aadb0 100644 --- a/opm/autodiff/StandardWellsSolvent.hpp +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -34,6 +34,7 @@ namespace Opm { public: using Base = StandardWells; + using Base::computeWellConnectionDensitesPressures; // --------- Public methods --------- explicit StandardWellsSolvent(const Wells* wells); @@ -63,6 +64,16 @@ namespace Opm { const std::vector& active, std::vector& mob_perfcells, std::vector& b_perfcells) const; + + template + void computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity); + protected: const SolventPropsAdFromDeck* solvent_props_; int solvent_pos_; diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp index 3529e8708..fa5babbac 100644 --- a/opm/autodiff/StandardWellsSolvent_impl.hpp +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -197,6 +197,42 @@ namespace Opm + template + void + StandardWellsSolvent:: + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition, + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Extract well connection depths + // TODO: depth_perf should be a member of the StandardWells class + const Vector pdepth = subset(depth, wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + + + + template void StandardWellsSolvent:: diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 1ccba74f7..be11a96a6 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -289,6 +289,42 @@ namespace Opm + template + void + StandardWells:: + computeWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& phaseCondition, + const Vector& depth, + const double gravity) + { + if( ! localWellsActive() ) return ; + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + std::vector b_perf; + std::vector rsmax_perf; + std::vector rvmax_perf; + std::vector surf_dens_perf; + computePropertiesForWellConnectionPressures(state, xw, fluid, active, phaseCondition, + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + // Extract well connection depths + // TODO: depth_perf should be a member of the StandardWells class + const Vector pdepth = subset(depth, wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); + + computeWellConnectionDensitesPressures(xw, fluid, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, gravity); + + } + + + + + template void StandardWells:: @@ -778,4 +814,373 @@ namespace Opm } + + + + + template + void + StandardWells:: + addWellFluxEq(const std::vector& cq_s, + const SolutionState& state, + LinearisedBlackoilResidual& residual) + { + if( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + ADB qs = state.qs; + for (int phase = 0; phase < np; ++phase) { + qs -= superset(wellOps().p2w * cq_s[phase], Span(nw, 1, phase*nw), nw*np); + + } + + residual.well_flux_eq = qs; + } + + + + + + template + void + StandardWells::addWellControlEq(const SolutionState& state, + const WellState& xw, + const Vector& aliveWells, + const std::vector active, + const VFPProperties& vfp_properties, + const double gravity, + LinearisedBlackoilResidual& residual) + { + if( ! localWellsActive() ) return; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + ADB aqua = ADB::constant(Vector::Zero(nw)); + ADB liquid = ADB::constant(Vector::Zero(nw)); + ADB vapour = ADB::constant(Vector::Zero(nw)); + + if (active[Water]) { + aqua += subset(state.qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); + } + if (active[Oil]) { + liquid += subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); + } + if (active[Gas]) { + vapour += subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); + } + + //THP calculation variables + std::vector inj_table_id(nw, -1); + std::vector prod_table_id(nw, -1); + Vector thp_inj_target_v = Vector::Zero(nw); + Vector thp_prod_target_v = Vector::Zero(nw); + Vector alq_v = Vector::Zero(nw); + + //Hydrostatic correction variables + Vector rho_v = Vector::Zero(nw); + Vector vfp_ref_depth_v = Vector::Zero(nw); + + //Target vars + Vector bhp_targets = Vector::Zero(nw); + Vector rate_targets = Vector::Zero(nw); + Eigen::SparseMatrix rate_distr(nw, np*nw); + + //Selection variables + std::vector bhp_elems; + std::vector thp_inj_elems; + std::vector thp_prod_elems; + std::vector rate_elems; + + //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]; + + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + const int current = xw.currentControls()[w]; + + switch (well_controls_iget_type(wc, current)) { + case BHP: + { + bhp_elems.push_back(w); + bhp_targets(w) = well_controls_iget_target(wc, current); + rate_targets(w) = -1e100; + } + break; + + case THP: + { + const int perf = wells().well_connpos[w]; + rho_v[w] = wellPerforationDensities()[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]; + if (well_type == INJECTOR) { + inj_table_id[w] = table_id; + thp_inj_target_v[w] = target; + alq_v[w] = -1e100; + + vfp_ref_depth_v[w] = vfp_properties.getInj()->getTable(table_id)->getDatumDepth(); + + thp_inj_elems.push_back(w); + } + else if (well_type == PRODUCER) { + prod_table_id[w] = table_id; + thp_prod_target_v[w] = target; + alq_v[w] = well_controls_iget_alq(wc, current); + + vfp_ref_depth_v[w] = vfp_properties.getProd()->getTable(table_id)->getDatumDepth(); + + thp_prod_elems.push_back(w); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER type well"); + } + bhp_targets(w) = -1e100; + rate_targets(w) = -1e100; + } + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + { + rate_elems.push_back(w); + // 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]; + } + + bhp_targets(w) = -1.0e100; + rate_targets(w) = well_controls_iget_target(wc, current); + } + break; + } + } + + //Calculate BHP target from THP + const ADB thp_inj_target = ADB::constant(thp_inj_target_v); + const ADB thp_prod_target = ADB::constant(thp_prod_target_v); + const ADB alq = ADB::constant(alq_v); + const ADB bhp_from_thp_inj = vfp_properties.getInj()->bhp(inj_table_id, aqua, liquid, vapour, thp_inj_target); + const ADB bhp_from_thp_prod = vfp_properties.getProd()->bhp(prod_table_id, aqua, liquid, vapour, thp_prod_target, alq); + + //Perform hydrostatic correction to computed targets + // double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const Vector dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, wellPerforationDensities(), 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); + + //Calculate residuals + const ADB thp_inj_residual = state.bhp - bhp_from_thp_inj + dp_inj; + const ADB thp_prod_residual = state.bhp - bhp_from_thp_prod + dp_prod; + const ADB bhp_residual = state.bhp - bhp_targets; + const ADB rate_residual = rate_distr * state.qs - rate_targets; + + //Select the right residual for each well + residual.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, nw) + + superset(subset(thp_inj_residual, thp_inj_elems), thp_inj_elems, nw) + + superset(subset(thp_prod_residual, thp_prod_elems), thp_prod_elems, nw) + + superset(subset(rate_residual, rate_elems), rate_elems, nw); + + // For wells that are dead (not flowing), and therefore not communicating + // with the reservoir, we set the equation to be equal to the well's total + // flow. This will be a solution only if the target rate is also zero. + Eigen::SparseMatrix rate_summer(nw, np*nw); + for (int w = 0; w < nw; ++w) { + for (int phase = 0; phase < np; ++phase) { + rate_summer.insert(w, phase*nw + w) = 1.0; + } + } + Selector alive_selector(aliveWells, Selector::NotEqualZero); + residual.well_eq = alive_selector.select(residual.well_eq, rate_summer * state.qs); + // OPM_AD_DUMP(residual_.well_eq); + } + + + + + + template + void + StandardWells::computeWellPotentials(SolutionState& state0, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + const Opm::PhaseUsage& pu, + const std::vector active, + const VFPProperties& vfp_properties, + const bool compute_well_potentials, + const bool gravity, + 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; + 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"); + } + } + } + + } + + // use bhp limit from control + state0.bhp = ADB::constant(bhps); + + // compute well potentials + Vector aliveWells; + std::vector well_potentials; + computeWellFlux(state0, pu, active, 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); + } + + } + + + + + + void + StandardWells::variableStateWellIndices(std::vector& indices, + int& next) const + { + indices[Qs] = next++; + indices[Bhp] = next++; + } + + + + + + std::vector + StandardWells::variableWellStateIndices() const + { + // Black oil model standard is 5 equation. + // For the pure well solve, only the well equations are picked. + std::vector indices(5, -1); + int next = 0; + + variableStateWellIndices(indices, next); + + assert(next == 2); + return indices; + } + + + + + + template + void + StandardWells::variableWellStateInitials(const WellState& xw, + std::vector& vars0) const + { + // Initial well rates. + if ( localWellsActive() ) + { + // 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; + + // 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); + 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()); + vars0.push_back(bhp); + } + else + { + // push null states for qs and bhp + vars0.push_back(V()); + vars0.push_back(V()); + } + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 71f513552..e63c80c58 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -202,8 +202,8 @@ namespace Opm { using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; - using Base::addWellControlEq; + // using Base::computeWellConnectionPressures; + // using Base::addWellControlEq; using Base::computeRelPerm; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c6535023c..fce5c27b1 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -499,6 +499,7 @@ namespace Opm { // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + const V depth = cellCentroidsZToEigen(grid_); // updateWellControls(well_state); stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); @@ -512,7 +513,8 @@ namespace Opm { // Compute initial accumulation contributions // and well connection pressures. computeAccum(state0, 0); - computeWellConnectionPressures(state0, well_state); + // computeWellConnectionPressures(state0, well_state); + stdWells().computeWellConnectionPressures(state0, well_state, fluid_, active_, phaseCondition(), depth, gravity); } // OPM_AD_DISKVAL(state.pressure); @@ -573,9 +575,9 @@ namespace Opm { stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); - Base::addWellFluxEq(cq_s, state); + stdWells().addWellFluxEq(cq_s, state, residual_); addWellContributionToMassBalanceEq(cq_s, state, well_state); - addWellControlEq(state, well_state, aliveWells); + stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); }