diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index fa8a01dc6..3d82cf1f0 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -206,9 +206,6 @@ namespace Opm const Scalar& segment_pressure, const bool& allow_cf, std::vector& cq_s, - Scalar& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, DeferredLogger& deferred_logger) const; void computePerfRateEval(const IntensiveQuantities& int_quants, @@ -261,10 +258,15 @@ namespace Opm DeferredLogger& deferred_logger) const; void computeWellRatesWithBhp(const Simulator& ebosSimulator, - const Scalar bhp, + const Scalar& bhp, std::vector& well_flux, DeferredLogger& deferred_logger) const; + void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator, + const Scalar& bhp, + std::vector& well_flux, + DeferredLogger& deferred_logger) const; + std::vector computeWellPotentialWithTHP(const Simulator& ebos_simulator, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index bb06aa61b..5f19edbc5 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -316,23 +316,64 @@ namespace Opm { if (this->well_ecl_.isInjector()) { const auto controls = this->well_ecl_.injectionControls(ebosSimulator.vanguard().summaryState()); - computeWellRatesWithBhp(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger); + computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger); } else { const auto controls = this->well_ecl_.productionControls(ebosSimulator.vanguard().summaryState()); - computeWellRatesWithBhp(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger); + computeWellRatesWithBhpIterations(ebosSimulator, controls.bhp_limit, well_flux, deferred_logger); } } - - template void MultisegmentWell:: computeWellRatesWithBhp(const Simulator& ebosSimulator, - const Scalar bhp, + const Scalar& bhp, std::vector& well_flux, DeferredLogger& deferred_logger) const { + + const int np = this->number_of_phases_; + + well_flux.resize(np, 0.0); + const bool allow_cf = this->getAllowCrossFlow(); + const int nseg = this->numberOfSegments(); + const WellState& well_state = ebosSimulator.problem().wellModel().wellState(); + const auto& ws = well_state.well(this->indexOfWell()); + auto segments_copy = ws.segments; + segments_copy.scale_pressure(bhp); + const auto& segment_pressure = segments_copy.pressure; + for (int seg = 0; seg < nseg; ++seg) { + for (const int perf : this->segment_perforations_[seg]) { + const int cell_idx = this->well_cells_[perf]; + const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); + // flux for each perforation + std::vector mob(this->num_components_, 0.); + getMobilityScalar(ebosSimulator, perf, mob); + double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); + const double Tw = this->well_index_[perf] * trans_mult; + + const Scalar seg_pressure = segment_pressure[seg]; + std::vector cq_s(this->num_components_, 0.); + computePerfRateScalar(intQuants, mob, Tw, seg, perf, seg_pressure, + allow_cf, cq_s, deferred_logger); + + for(int p = 0; p < np; ++p) { + well_flux[this->ebosCompIdxToFlowCompIdx(p)] += cq_s[p]; + } + } + } + this->parallel_well_info_.communication().sum(well_flux.data(), well_flux.size()); + } + + + template + void + MultisegmentWell:: + computeWellRatesWithBhpIterations(const Simulator& ebosSimulator, + const Scalar& bhp, + std::vector& well_flux, + DeferredLogger& deferred_logger) const + { // creating a copy of the well itself, to avoid messing up the explicit informations // during this copy, the only information not copied properly is the well controls MultisegmentWell well_copy(*this); @@ -404,7 +445,7 @@ namespace Opm if (bhp_at_thp_limit) { const auto& controls = well.injectionControls(summary_state); const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit); - computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); deferred_logger.debug("Converged thp based potential calculation for well " + this->name() + ", at bhp = " + std::to_string(bhp)); } else { @@ -413,14 +454,14 @@ namespace Opm + this->name() + ". Instead the bhp based value is used"); const auto& controls = well.injectionControls(summary_state); const double bhp = controls.bhp_limit; - computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } } else { auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger); if (bhp_at_thp_limit) { const auto& controls = well.productionControls(summary_state); const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit); - computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); deferred_logger.debug("Converged thp based potential calculation for well " + this->name() + ", at bhp = " + std::to_string(bhp)); } else { @@ -429,7 +470,7 @@ namespace Opm + this->name() + ". Instead the bhp based value is used"); const auto& controls = well.productionControls(summary_state); const double bhp = controls.bhp_limit; - computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } } @@ -883,9 +924,6 @@ namespace Opm const Scalar& segment_pressure, const bool& allow_cf, std::vector& cq_s, - Scalar& perf_press, - double& perf_dis_gas_rate, - double& perf_vap_oil_rate, DeferredLogger& deferred_logger) const { @@ -912,6 +950,10 @@ namespace Opm cmix_s[comp_idx] = getValue(this->surfaceVolumeFraction(seg, comp_idx)); } + Scalar perf_dis_gas_rate = 0.0; + Scalar perf_vap_oil_rate = 0.0; + Scalar perf_press = 0.0; + this->computePerfRate(pressure_cell, rs, rv, @@ -1701,12 +1743,32 @@ namespace Opm return rates; }; - return this->MultisegmentWellGeneric:: + auto bhpAtLimit = this->MultisegmentWellGeneric:: computeBhpAtThpLimitProd(frates, summary_state, maxPerfPress(ebos_simulator), getRefDensity(), deferred_logger); + + if(bhpAtLimit) + return bhpAtLimit; + + auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) { + // Solver the well iterations to see if we are + // able to get a solution with an update + // solution + std::vector rates(3); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger); + return rates; + }; + + return this->MultisegmentWellGeneric:: + computeBhpAtThpLimitProd(fratesIter, + summary_state, + maxPerfPress(ebos_simulator), + getRefDensity(), + deferred_logger); + } @@ -1731,8 +1793,26 @@ namespace Opm return rates; }; + auto bhpAtLimit = this->MultisegmentWellGeneric:: + computeBhpAtThpLimitInj(frates, + summary_state, + getRefDensity(), + deferred_logger); + + if(bhpAtLimit) + return bhpAtLimit; + + auto fratesIter = [this, &ebos_simulator, &deferred_logger](const double bhp) { + // Solver the well iterations to see if we are + // able to get a solution with an update + // solution + std::vector rates(3); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, rates, deferred_logger); + return rates; + }; + return this->MultisegmentWellGeneric:: - computeBhpAtThpLimitInj(frates, summary_state, getRefDensity(), deferred_logger); + computeBhpAtThpLimitInj(fratesIter, summary_state, getRefDensity(), deferred_logger); } @@ -1783,10 +1863,7 @@ namespace Opm const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); - Scalar perf_press; - double perf_dis_gas_rate = 0.; - double perf_vap_oil_rate = 0.; - computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, perf_press, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger); + computePerfRateScalar(int_quants, mob, Tw, seg, perf, seg_pressure, allow_cf, cq_s, deferred_logger); for (int comp = 0; comp < this->num_components_; ++comp) { well_q_s[comp] += cq_s[comp]; } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 8e842ce28..68a7dbccd 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -327,7 +327,7 @@ namespace Opm double& perf_vap_oil_rate, DeferredLogger& deferred_logger) const; - void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator, + void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator, const double& bhp, std::vector& well_flux, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 13f9202a4..7863fdb73 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1673,10 +1673,10 @@ namespace Opm template void StandardWell:: - computeWellRatesWithBhpPotential(const Simulator& ebosSimulator, - const double& bhp, - std::vector& well_flux, - DeferredLogger& deferred_logger) const + computeWellRatesWithBhpIterations(const Simulator& ebosSimulator, + const double& bhp, + std::vector& well_flux, + DeferredLogger& deferred_logger) const { // iterate to get a more accurate well density @@ -1738,13 +1738,13 @@ namespace Opm auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger); if (bhp_at_thp_limit) { const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit); - computeWellRatesWithBhpPotential(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } else { deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL", "Failed in getting converged thp based potential calculation for well " + name() + ". Instead the bhp based value is used"); const double bhp = controls.bhp_limit; - computeWellRatesWithBhpPotential(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } } else { computeWellRatesWithThpAlqProd( @@ -1771,7 +1771,7 @@ namespace Opm if (bhp_at_thp_limit) { const auto& controls = this->well_ecl_.productionControls(summary_state); bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit); - computeWellRatesWithBhpPotential(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } else { deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL", @@ -1779,7 +1779,7 @@ namespace Opm + name() + ". Instead the bhp based value is used"); const auto& controls = this->well_ecl_.productionControls(summary_state); bhp = controls.bhp_limit; - computeWellRatesWithBhpPotential(ebos_simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger); } return bhp; } @@ -1890,7 +1890,7 @@ namespace Opm // get the bhp value based on the bhp constraints const double bhp = this->mostStrictBhpFromBhpLimits(summaryState); assert(std::abs(bhp) != std::numeric_limits::max()); - computeWellRatesWithBhpPotential(ebosSimulator, bhp, well_potentials, deferred_logger); + computeWellRatesWithBhpIterations(ebosSimulator, bhp, well_potentials, deferred_logger); } else { // the well has a THP related constraint well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);