From eb278c3c9a405bd09ae808b5ca16921d9a93e716 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 18 Apr 2016 11:22:19 +0200 Subject: [PATCH 1/7] moving addWellFluxEq() to StandardWells Conflicts: opm/autodiff/BlackoilModelBase_impl.hpp --- opm/autodiff/BlackoilModelBase.hpp | 7 ---- opm/autodiff/BlackoilModelBase_impl.hpp | 32 ++----------------- opm/autodiff/StandardWells.hpp | 7 ++++ opm/autodiff/StandardWells_impl.hpp | 29 +++++++++++++++++ .../BlackoilPolymerModel_impl.hpp | 2 +- 5 files changed, 39 insertions(+), 38 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index d39e4e719..1e9d64795 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -410,13 +410,6 @@ namespace Opm { 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); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 32d7bee7b..255cab70d 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -901,7 +901,7 @@ 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); @@ -1069,34 +1069,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:: @@ -1174,7 +1146,7 @@ namespace detail { 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().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); converged = getWellConvergence(it); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 33a15812f..badf9941a 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -141,6 +141,13 @@ namespace Opm { 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); + + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 1ccba74f7..36a9b0795 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -778,4 +778,33 @@ 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; + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c6535023c..306b15392 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -573,7 +573,7 @@ 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); } From dcca0b0b76bb07bba490ba6cabf0f7a8501a3992 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 19 Apr 2016 11:20:18 +0200 Subject: [PATCH 2/7] moving addWellControlEq to StandardWells --- opm/autodiff/BlackoilModelBase.hpp | 5 - opm/autodiff/BlackoilModelBase_impl.hpp | 175 +----------------- opm/autodiff/BlackoilSolventModel.hpp | 2 +- opm/autodiff/StandardWells.hpp | 10 + opm/autodiff/StandardWells_impl.hpp | 171 +++++++++++++++++ .../fullyimplicit/BlackoilPolymerModel.hpp | 2 +- .../BlackoilPolymerModel_impl.hpp | 2 +- 7 files changed, 188 insertions(+), 179 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 1e9d64795..b8cf69a88 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -414,11 +414,6 @@ namespace Opm { 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 255cab70d..5453d1d87 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -903,7 +903,7 @@ namespace detail { asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().stdWells().addWellFluxEq(cq_s, state, residual_); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); - asImpl().addWellControlEq(state, well_state, aliveWells); + asImpl().stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state); } @@ -1147,7 +1147,9 @@ namespace detail { 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().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_); - asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells, + active_, vfp_properties_, gravity, residual_); converged = getWellConvergence(it); if (converged) { @@ -1171,7 +1173,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); } @@ -1214,174 +1215,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:: diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index ddc3aa40c..1e0409667 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -147,7 +147,7 @@ namespace Opm { using Base::maxResidualAllowed; // using Base::updateWellControls; using Base::computeWellConnectionPressures; - using Base::addWellControlEq; + // using Base::addWellControlEq; // using Base::computePropertiesForWellConnectionPressures; std::vector diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index badf9941a..0b1d7e591 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -147,6 +147,16 @@ namespace Opm { 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); + protected: diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 36a9b0795..49efed601 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -807,4 +807,175 @@ namespace Opm 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); + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 71f513552..2f4c1a0e7 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -203,7 +203,7 @@ namespace Opm { // using Base::updateWellControls; using Base::computeWellConnectionPressures; - using Base::addWellControlEq; + // using Base::addWellControlEq; using Base::computeRelPerm; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 306b15392..d6cb9b1cd 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -575,7 +575,7 @@ namespace Opm { stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_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_); } From 13acc8ee03452725e4c4c00b6d47d86e9f22d2a8 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 19 Apr 2016 17:09:50 +0200 Subject: [PATCH 3/7] moving computeWellConnectionPressures to StandardWells the results look okay, while the running for flow_solvent needs further investigation even the results with flow_solvent actually look okay. With two different version of computePropertiesForWellConnectionPressures, flow_solvent produces the same results. This is something needs further investigation. The current implementation requires a copy of computeWellConnectionPressure in StandardWells and StandardWellsSolvent. That means probably we need to introduce the asImpl() for the Wells. --- opm/autodiff/BlackoilModelBase.hpp | 5 +- opm/autodiff/BlackoilModelBase_impl.hpp | 54 +++++-------------- opm/autodiff/BlackoilSolventModel.hpp | 2 +- opm/autodiff/StandardWells.hpp | 52 ++++++++++-------- opm/autodiff/StandardWellsSolvent.hpp | 11 ++++ opm/autodiff/StandardWellsSolvent_impl.hpp | 36 +++++++++++++ opm/autodiff/StandardWells_impl.hpp | 36 +++++++++++++ .../fullyimplicit/BlackoilPolymerModel.hpp | 2 +- .../BlackoilPolymerModel_impl.hpp | 4 +- 9 files changed, 135 insertions(+), 67 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8cf69a88..8d60bf619 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -380,9 +380,6 @@ namespace Opm { computeAccum(const SolutionState& state, const int aix ); - void computeWellConnectionPressures(const SolutionState& state, - const WellState& xw); - void assembleMassBalanceEq(const SolutionState& state); @@ -410,6 +407,8 @@ namespace Opm { const std::vector& mob_perfcells, const std::vector& b_perfcells, WellState& well_state); + + void addWellContributionToMassBalanceEq(const std::vector& cq_s, const SolutionState& state, const WellState& xw); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 5453d1d87..74cb2e980 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -798,41 +798,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 +807,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 +817,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 +838,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); @@ -1133,6 +1103,9 @@ namespace detail { } } + // gravity + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + int it = 0; bool converged; do { @@ -1147,7 +1120,6 @@ namespace detail { 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().stdWells().addWellFluxEq(cq_s, wellSolutionState, residual_); - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); asImpl().stdWells().addWellControlEq(wellSolutionState, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); converged = getWellConvergence(it); @@ -1201,7 +1173,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) { diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 1e0409667..d3b09322e 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -146,7 +146,7 @@ namespace Opm { using Base::drMaxRel; using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; + // using Base::computeWellConnectionPressures; // using Base::addWellControlEq; // using Base::computePropertiesForWellConnectionPressures; diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 0b1d7e591..c10f5ba85 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -133,29 +133,39 @@ 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: 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); - // 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); 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 49efed601..2298d2661 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:: diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 2f4c1a0e7..e63c80c58 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -202,7 +202,7 @@ namespace Opm { using Base::maxResidualAllowed; // using Base::updateWellControls; - using Base::computeWellConnectionPressures; + // 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 d6cb9b1cd..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); From b87920e5e494e9b83269e92b9434edcade25ab15 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 21 Apr 2016 15:46:20 +0200 Subject: [PATCH 4/7] moving computeWellPotentials to StandardWells --- opm/autodiff/BlackoilModelBase.hpp | 6 -- opm/autodiff/BlackoilModelBase_impl.hpp | 114 ++---------------------- opm/autodiff/StandardWells.hpp | 12 +++ opm/autodiff/StandardWells_impl.hpp | 104 +++++++++++++++++++++ 4 files changed, 124 insertions(+), 112 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 8d60bf619..23d3437ad 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -402,12 +402,6 @@ 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 addWellContributionToMassBalanceEq(const std::vector& cq_s, const SolutionState& state, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 74cb2e980..3feae3d46 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -874,7 +874,14 @@ namespace detail { asImpl().stdWells().addWellFluxEq(cq_s, state, residual_); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().stdWells().addWellControlEq(state, well_state, aliveWells, active_, vfp_properties_, gravity, residual_); - asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state); + // 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); + } } @@ -1546,111 +1553,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/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index c10f5ba85..d094cb138 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -166,6 +166,18 @@ namespace Opm { 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); diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 2298d2661..7ca1cbc61 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -1014,4 +1014,108 @@ namespace Opm // 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); + } + + } + } From e7d00f4f99470306e23aba555aad9b10ab5833a4 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 21 Apr 2016 17:42:50 +0200 Subject: [PATCH 5/7] adding variableStateWellIndices to StandardWells to handle different types of variables later for different types of wells. --- opm/autodiff/BlackoilModelBase_impl.hpp | 3 +-- opm/autodiff/StandardWells.hpp | 5 +++++ opm/autodiff/StandardWells_impl.hpp | 12 ++++++++++++ 3 files changed, 18 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 3feae3d46..0765a0f0e 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -604,8 +604,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; } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index d094cb138..0f894a242 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -180,6 +180,11 @@ namespace Opm { WellState& well_state); + void + variableStateWellIndices(std::vector& indices, + int& next) const; + + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 7ca1cbc61..44a2afb73 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -1118,4 +1118,16 @@ namespace Opm } + + + + + void + StandardWells::variableStateWellIndices(std::vector& indices, + int& next) const + { + indices[Qs] = next++; + indices[Bhp] = next++; + } + } From 15380fd3709a171d104ebd573c67187bc14ab045 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 21 Apr 2016 18:03:16 +0200 Subject: [PATCH 6/7] moving variableWellStateIndices to StandardWells Probably it is the time to introduce SeqQs and SeqP for the multi-segment wells. --- opm/autodiff/BlackoilModelBase.hpp | 3 --- opm/autodiff/BlackoilModelBase_impl.hpp | 20 +------------------ opm/autodiff/BlackoilMultiSegmentModel.hpp | 2 +- .../BlackoilMultiSegmentModel_impl.hpp | 2 +- opm/autodiff/StandardWells.hpp | 2 ++ opm/autodiff/StandardWells_impl.hpp | 18 +++++++++++++++++ 6 files changed, 23 insertions(+), 24 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 23d3437ad..301b637ed 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -363,9 +363,6 @@ namespace Opm { std::vector variableStateIndices() const; - std::vector - variableWellStateIndices() const; - SolutionState variableStateExtractVars(const ReservoirState& x, const std::vector& indices, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 0765a0f0e..97a495c06 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -612,24 +612,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 @@ -1092,7 +1074,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); diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index 3ff7094b2..76dd8befd 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -223,7 +223,7 @@ namespace Opm { using Base::convergenceReduction; using Base::maxResidualAllowed; using Base::variableState; - using Base::variableWellStateIndices; + // using Base::variableWellStateIndices; using Base::asImpl; const std::vector& wellsMultiSegment() const { return wells_multisegment_; } diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index c87c2c7ac..3741d4a15 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); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 0f894a242..0e41ef5fa 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -185,6 +185,8 @@ namespace Opm { int& next) const; + std::vector + variableWellStateIndices() const; protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 44a2afb73..24fc8f8a6 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -1130,4 +1130,22 @@ namespace Opm 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; + } + } From 08e691e262f08c1fde3b72e4e966d5f4049332df Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 22 Apr 2016 10:51:17 +0200 Subject: [PATCH 7/7] moving variableWellStateInitials to StandardWells. --- opm/autodiff/BlackoilModelBase.hpp | 3 -- opm/autodiff/BlackoilModelBase_impl.hpp | 39 +------------------ opm/autodiff/BlackoilMultiSegmentModel.hpp | 5 +++ .../BlackoilMultiSegmentModel_impl.hpp | 22 +++++++++++ opm/autodiff/StandardWells.hpp | 5 +++ opm/autodiff/StandardWells_impl.hpp | 35 +++++++++++++++++ 6 files changed, 69 insertions(+), 40 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 301b637ed..f5341eae5 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -356,9 +356,6 @@ namespace Opm { void variableReservoirStateInitials(const ReservoirState& x, std::vector& vars0) const; - void - variableWellStateInitials(const WellState& xw, - std::vector& vars0) const; std::vector variableStateIndices() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 97a495c06..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:: @@ -1100,7 +1065,7 @@ namespace detail { // 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; diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index 76dd8befd..4d659c920 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -225,6 +225,7 @@ namespace Opm { using Base::variableState; // 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 3741d4a15..3026df908 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -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/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 0e41ef5fa..9f12a9596 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -188,6 +188,11 @@ namespace Opm { std::vector variableWellStateIndices() const; + template + void + variableWellStateInitials(const WellState& xw, + std::vector& vars0) const; + protected: bool wells_active_; const Wells* wells_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 24fc8f8a6..be11a96a6 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -1148,4 +1148,39 @@ namespace Opm 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()); + } + } + }