From c623fba017847aabc01c2a3382fbec76875b321b Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 3 Jun 2021 15:34:14 +0200 Subject: [PATCH] add MultisegmentWellGeneric avoid rebuilding this for all simulators when code is only dependent on Scalar. instanced for double --- CMakeLists_files.cmake | 1 + opm/simulators/wells/MultisegmentWell.hpp | 48 +- .../wells/MultisegmentWellGeneric.cpp | 643 +++++++++++++++ .../wells/MultisegmentWellGeneric.hpp | 113 +++ .../wells/MultisegmentWell_impl.hpp | 752 ++---------------- opm/simulators/wells/WellInterfaceGeneric.hpp | 8 + 6 files changed, 852 insertions(+), 713 deletions(-) create mode 100644 opm/simulators/wells/MultisegmentWellGeneric.cpp create mode 100644 opm/simulators/wells/MultisegmentWellGeneric.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 79c216c53..f1f2cf053 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -60,6 +60,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/GasLiftStage2.cpp opm/simulators/wells/GlobalWellInfo.cpp opm/simulators/wells/GroupState.cpp + opm/simulators/wells/MultisegmentWellGeneric.cpp opm/simulators/wells/ParallelWellInfo.cpp opm/simulators/wells/SegmentState.cpp opm/simulators/wells/StandardWellGeneric.cpp diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 6d4d98805..d29a40532 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -23,6 +23,7 @@ #define OPM_MULTISEGMENTWELL_HEADER_INCLUDED #include +#include #include @@ -31,7 +32,8 @@ namespace Opm class DeferredLogger; template - class MultisegmentWell: public WellInterface + class MultisegmentWell : public WellInterface + , public MultisegmentWellGeneric> { public: typedef WellInterface Base; @@ -195,10 +197,6 @@ namespace Opm virtual void addWellContributions(SparseMatrixAdapter& jacobian) const override; - /// number of segments for this well - /// int number_of_segments_; - int numberOfSegments() const; - virtual std::vector computeCurrentWellRates(const Simulator& ebosSimulator, DeferredLogger& deferred_logger) const override; @@ -225,9 +223,6 @@ namespace Opm // multi-phase flow model WellSegments::MultiPhaseModel multiphaseModel() const; - // get the WellSegments from the well_ecl_ - const WellSegments& segmentSet() const; - // protected member variables from the Base class using Base::well_ecl_; using Base::vfp_properties_; @@ -267,24 +262,6 @@ namespace Opm using Base::calculateBhpFromThp; using Base::getALQ; - // TODO: trying to use the information from the Well opm-parser as much - // as possible, it will possibly be re-implemented later for efficiency reason. - - // the completions that is related to each segment - // the completions's ids are their index in the vector well_index_, well_cell_ - // This is also assuming the order of the completions in Well is the same with - // the order of the completions in wells. - // it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations - // belonging to this segment - std::vector > segment_perforations_; - - // the inlet segments for each segment. It is for convinience and efficiency reason - std::vector > segment_inlets_; - - // segment number is an ID of the segment, it is specified in the deck - // get the loation of the segment with a segment number in the segmentSet - int segmentNumberToIndex(const int segment_number) const; - // TODO, the following should go to a class for computing purpose // two off-diagonal matrices mutable OffDiagMatWell duneB_; @@ -312,11 +289,6 @@ namespace Opm // center depth of the grid block std::vector cell_perforation_pressure_diffs_; - // depth difference between the segment and the peforation - // or in another way, the depth difference between the perforation and - // the segment the perforation belongs to - std::vector perforation_segment_depth_diffs_; - // the intial amount of fluids in each segment under surface condition std::vector > segment_fluid_initial_; @@ -330,8 +302,6 @@ namespace Opm // the mass rate of the segments std::vector segment_mass_rates_; - std::vector segment_depth_diffs_; - // the upwinding segment for each segment based on the flow direction std::vector upwinding_segments_; @@ -360,10 +330,6 @@ namespace Opm const double relaxation_factor=1.0) const; - // scale the segment rates and pressure based on well rates and bhp - void scaleSegmentRatesWithWellRates(WellState& well_state) const; - void scaleSegmentPressuresWithBhp(WellState& well_state) const; - // computing the accumulation term for later use in well mass equations void computeInitialSegmentFluids(const Simulator& ebos_simulator); @@ -396,7 +362,6 @@ namespace Opm // convert a Eval from reservoir to contain the derivative related to wells EvalWell extendEval(const Eval& in) const; - double calculateThpFromBhp(const std::vector& rates, const double bhp, DeferredLogger& deferred_logger) const; void updateThp(WellState& well_state, DeferredLogger& deferred_logger) const; // compute the fluid properties, such as densities, viscosities, and so on, in the segments @@ -457,10 +422,6 @@ namespace Opm virtual double getRefDensity() const override; - bool frictionalPressureLossConsidered() const; - - bool accelerationalPressureLossConsidered() const; - virtual bool iterateWellEqWithControl(const Simulator& ebosSimulator, const double dt, const Well::InjectionControls& inj_controls, @@ -484,9 +445,6 @@ namespace Opm std::vector getWellResiduals(const std::vector& B_avg, DeferredLogger& deferred_logger) const; - void detectOscillations(const std::vector& measure_history, - const int it, bool& oscillate, bool& stagnate) const; - double getResidualMeasureValue(const WellState& well_state, const std::vector& residuals, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp new file mode 100644 index 000000000..935d4f0c2 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -0,0 +1,643 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace Opm +{ + +template +MultisegmentWellGeneric:: +MultisegmentWellGeneric(WellInterfaceGeneric& baseif) + : baseif_(baseif) + , segment_perforations_(numberOfSegments()) + , segment_inlets_(numberOfSegments()) + , segment_depth_diffs_(numberOfSegments(), 0.0) + , perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0) +{ + // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. + // for other facilities needed but not available from parser, we need to process them here + + // initialize the segment_perforations_ and update perforation_segment_depth_diffs_ + const WellConnections& completion_set = baseif_.wellEcl().getConnections(); + // index of the perforation within wells struct + // there might be some perforations not active, which causes the number of the perforations in + // well_ecl_ and wells struct different + // the current implementation is a temporary solution for now, it should be corrected from the parser + // side + int i_perf_wells = 0; + baseif.perfDepth().resize(baseif_.numPerfs(), 0.); + for (size_t perf = 0; perf < completion_set.size(); ++perf) { + const Connection& connection = completion_set.get(perf); + if (connection.state() == Connection::State::OPEN) { + const int segment_index = segmentNumberToIndex(connection.segment()); + segment_perforations_[segment_index].push_back(i_perf_wells); + baseif.perfDepth()[i_perf_wells] = connection.depth(); + const double segment_depth = segmentSet()[segment_index].depth(); + perforation_segment_depth_diffs_[i_perf_wells] = baseif.perfDepth()[i_perf_wells] - segment_depth; + i_perf_wells++; + } + } + + // initialize the segment_inlets_ + for (int seg = 0; seg < numberOfSegments(); ++seg) { + const Segment& segment = segmentSet()[seg]; + const int segment_number = segment.segmentNumber(); + const int outlet_segment_number = segment.outletSegment(); + if (outlet_segment_number > 0) { + const int segment_index = segmentNumberToIndex(segment_number); + const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number); + segment_inlets_[outlet_segment_index].push_back(segment_index); + } + } + + // calculating the depth difference between the segment and its oulet_segments + // for the top segment, we will make its zero unless we find other purpose to use this value + for (int seg = 1; seg < numberOfSegments(); ++seg) { + const double segment_depth = segmentSet()[seg].depth(); + const int outlet_segment_number = segmentSet()[seg].outletSegment(); + const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)]; + const double outlet_depth = outlet_segment.depth(); + segment_depth_diffs_[seg] = segment_depth - outlet_depth; + } +} + +template +void +MultisegmentWellGeneric:: +scaleSegmentRatesWithWellRates(WellState& well_state) const +{ + auto& segments = well_state.segments(baseif_.indexOfWell()); + auto& segment_rates = segments.rates; + for (int phase = 0; phase < baseif_.numPhases(); ++phase) { + const double unscaled_top_seg_rate = segment_rates[phase]; + const double well_phase_rate = well_state.wellRates(baseif_.indexOfWell())[phase]; + if (std::abs(unscaled_top_seg_rate) > 1e-12) + { + for (int seg = 0; seg < numberOfSegments(); ++seg) { + segment_rates[baseif_.numPhases()*seg + phase] *= well_phase_rate/unscaled_top_seg_rate; + } + } else { + // for newly opened wells, the unscaled rate top segment rate is zero + // and we need to initialize the segment rates differently + double sumTw = 0; + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + sumTw += baseif_.wellIndex()[perf]; + } + + std::vector perforation_rates(baseif_.numPhases() * baseif_.numPerfs(),0.0); + const double perf_phaserate_scaled = well_state.wellRates(baseif_.indexOfWell())[phase] / sumTw; + for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + perforation_rates[baseif_.numPhases()* perf + phase] = baseif_.wellIndex()[perf] * perf_phaserate_scaled; + } + + std::vector rates; + WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, baseif_.numPhases(), 0, rates); + std::copy(rates.begin(), rates.end(), segment_rates.begin()); + } + } +} + +template +void +MultisegmentWellGeneric:: +scaleSegmentPressuresWithBhp(WellState& well_state) const +{ + auto& segments = well_state.segments(baseif_.indexOfWell()); + auto bhp = well_state.bhp(baseif_.indexOfWell()); + segments.scale_pressure(bhp); +} + +template +const WellSegments& +MultisegmentWellGeneric:: +segmentSet() const +{ + return baseif_.wellEcl().getSegments(); +} + +template +int +MultisegmentWellGeneric:: +numberOfSegments() const +{ + return segmentSet().size(); +} + +template +WellSegments::CompPressureDrop +MultisegmentWellGeneric:: +compPressureDrop() const +{ + return segmentSet().compPressureDrop(); +} + + +template +int +MultisegmentWellGeneric:: +segmentNumberToIndex(const int segment_number) const +{ + return segmentSet().segmentNumberToIndex(segment_number); +} + +template +double +MultisegmentWellGeneric:: +calculateThpFromBhp(const std::vector& rates, + const double bhp, + const double rho, + DeferredLogger& deferred_logger) const +{ + assert(int(rates.size()) == 3); // the vfp related only supports three phases now. + + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Gas = BlackoilPhases::Vapour; + + const double aqua = rates[Water]; + const double liquid = rates[Oil]; + const double vapour = rates[Gas]; + + double thp = 0.0; + if (baseif_.isInjector()) { + const int table_id = baseif_.wellEcl().vfp_table_number(); + const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth(); + const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity()); + + thp = baseif_.vfpProperties()->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); + } + else if (baseif_.isProducer()) { + const int table_id = baseif_.wellEcl().vfp_table_number(); + const double alq = baseif_.wellEcl().alq_value(); + const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth(); + const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity()); + + thp = baseif_.vfpProperties()->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); + } + else { + OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger); + } + + return thp; +} + +template +void +MultisegmentWellGeneric:: +detectOscillations(const std::vector& measure_history, + const int it, + bool& oscillate, + bool& stagnate) const +{ + if ( it < 2 ) { + oscillate = false; + stagnate = false; + return; + } + + stagnate = true; + const double F0 = measure_history[it]; + const double F1 = measure_history[it - 1]; + const double F2 = measure_history[it - 2]; + const double d1 = std::abs((F0 - F2) / F0); + const double d2 = std::abs((F0 - F1) / F0); + + const double oscillaton_rel_tol = 0.2; + oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2); + + const double stagnation_rel_tol = 1.e-2; + stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol; +} + +template +std::optional +MultisegmentWellGeneric:: +computeBhpAtThpLimitInj(const std::function(const double)>& frates, + const SummaryState& summary_state, + const double rho, + DeferredLogger& deferred_logger) const +{ + // Given a VFP function returning bhp as a function of phase + // rates and thp: + // fbhp(rates, thp), + // a function extracting the particular flow rate used for VFP + // lookups: + // flo(rates) + // and the inflow function (assuming the reservoir is fixed): + // frates(bhp) + // we want to solve the equation: + // fbhp(frates(bhp, thplimit)) - bhp = 0 + // for bhp. + // + // This may result in 0, 1 or 2 solutions. If two solutions, + // the one corresponding to the lowest bhp (and therefore + // highest rate) is returned. + // + // In order to detect these situations, we will find piecewise + // linear approximations both to the inverse of the frates + // function and to the fbhp function. + // + // We first take the FLO sample points of the VFP curve, and + // find the corresponding bhp values by solving the equation: + // flo(frates(bhp)) - flo_sample = 0 + // for bhp, for each flo_sample. The resulting (flo_sample, + // bhp_sample) values give a piecewise linear approximation to + // the true inverse inflow function, at the same flo values as + // the VFP data. + // + // Then we extract a piecewise linear approximation from the + // multilinear fbhp() by evaluating it at the flo_sample + // points, with fractions given by the frates(bhp_sample) + // values. + // + // When we have both piecewise linear curves defined on the + // same flo_sample points, it is easy to distinguish between + // the 0, 1 or 2 solution cases, and obtain the right interval + // in which to solve for the solution we want (with highest + // flow in case of 2 solutions). + + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Gas = BlackoilPhases::Vapour; + + // Make the fbhp() function. + const auto& controls = baseif_.wellEcl().injectionControls(summary_state); + const auto& table = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number); + const double vfp_ref_depth = table.getDatumDepth(); + const double thp_limit = baseif_.getTHPConstraint(summary_state); + const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity()); + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { + assert(rates.size() == 3); + return baseif_.vfpProperties()->getInj() + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp; + }; + + // Make the flo() function. + auto flo = [&table](const std::vector& rates) { + return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); + }; + + // Get the flo samples, add extra samples at low rates and bhp + // limit point if necessary. + std::vector flo_samples = table.getFloAxis(); + if (flo_samples[0] > 0.0) { + const double f0 = flo_samples[0]; + flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 }); + } + const double flo_bhp_limit = flo(frates(controls.bhp_limit)); + if (flo_samples.back() < flo_bhp_limit) { + flo_samples.push_back(flo_bhp_limit); + } + + // Find bhp values for inflow relation corresponding to flo samples. + std::vector bhp_samples; + for (double flo_sample : flo_samples) { + if (flo_sample > flo_bhp_limit) { + // We would have to go over the bhp limit to obtain a + // flow of this magnitude. We associate all such flows + // with simply the bhp limit. The first one + // encountered is considered valid, the rest not. They + // are therefore skipped. + bhp_samples.push_back(controls.bhp_limit); + break; + } + auto eq = [&flo, &frates, flo_sample](double bhp) { + return flo(frates(bhp)) - flo_sample; + }; + // TODO: replace hardcoded low/high limits. + const double low = 10.0 * unit::barsa; + const double high = 800.0 * unit::barsa; + const int max_iteration = 100; + const double flo_tolerance = 0.05 * std::fabs(flo_samples.back()); + int iteration = 0; + try { + const double solved_bhp = RegulaFalsiBisection:: + solve(eq, low, high, max_iteration, flo_tolerance, iteration); + bhp_samples.push_back(solved_bhp); + } + catch (...) { + // Use previous value (or max value if at start) if we failed. + bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back()); + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES", + "Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + baseif_.name()); + } + } + + // Find bhp values for VFP relation corresponding to flo samples. + const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size() + std::vector fbhp_samples(num_samples); + for (int ii = 0; ii < num_samples; ++ii) { + fbhp_samples[ii] = fbhp(frates(bhp_samples[ii])); + } +// #define EXTRA_THP_DEBUGGING +#ifdef EXTRA_THP_DEBUGGING + std::string dbgmsg; + dbgmsg += "flo: "; + for (int ii = 0; ii < num_samples; ++ii) { + dbgmsg += " " + std::to_string(flo_samples[ii]); + } + dbgmsg += "\nbhp: "; + for (int ii = 0; ii < num_samples; ++ii) { + dbgmsg += " " + std::to_string(bhp_samples[ii]); + } + dbgmsg += "\nfbhp: "; + for (int ii = 0; ii < num_samples; ++ii) { + dbgmsg += " " + std::to_string(fbhp_samples[ii]); + } + OpmLog::debug(dbgmsg); +#endif // EXTRA_THP_DEBUGGING + + // Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve. + // We only look at the valid + int sign_change_index = -1; + for (int ii = 0; ii < num_samples - 1; ++ii) { + const double curr = fbhp_samples[ii] - bhp_samples[ii]; + const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1]; + if (curr * next < 0.0) { + // Sign change in the [ii, ii + 1] interval. + sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution. + } + } + + // Handle the no solution case. + if (sign_change_index == -1) { + return std::optional(); + } + + // Solve for the proper solution in the given interval. + auto eq = [&fbhp, &frates](double bhp) { + return fbhp(frates(bhp)) - bhp; + }; + // TODO: replace hardcoded low/high limits. + const double low = bhp_samples[sign_change_index + 1]; + const double high = bhp_samples[sign_change_index]; + const int max_iteration = 100; + const double bhp_tolerance = 0.01 * unit::barsa; + int iteration = 0; + if (low == high) { + // We are in the high flow regime where the bhp_samples + // are all equal to the bhp_limit. + assert(low == controls.bhp_limit); + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", + "Robust bhp(thp) solve failed for well " + baseif_.name()); + return std::optional(); + } + try { + const double solved_bhp = RegulaFalsiBisection:: + solve(eq, low, high, max_iteration, bhp_tolerance, iteration); +#ifdef EXTRA_THP_DEBUGGING + OpmLog::debug("***** " + name() + " solved_bhp = " + std::to_string(solved_bhp) + + " flo_bhp_limit = " + std::to_string(flo_bhp_limit)); +#endif // EXTRA_THP_DEBUGGING + return solved_bhp; + } + catch (...) { + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", + "Robust bhp(thp) solve failed for well " + baseif_.name()); + return std::optional(); + } +} + +template +std::optional +MultisegmentWellGeneric:: +computeBhpAtThpLimitProd(const std::function(const double)>& frates, + const SummaryState& summary_state, + const double maxPerfPress, + const double rho, + DeferredLogger& deferred_logger) const +{ + // Given a VFP function returning bhp as a function of phase + // rates and thp: + // fbhp(rates, thp), + // a function extracting the particular flow rate used for VFP + // lookups: + // flo(rates) + // and the inflow function (assuming the reservoir is fixed): + // frates(bhp) + // we want to solve the equation: + // fbhp(frates(bhp, thplimit)) - bhp = 0 + // for bhp. + // + // This may result in 0, 1 or 2 solutions. If two solutions, + // the one corresponding to the lowest bhp (and therefore + // highest rate) should be returned. + + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Gas = BlackoilPhases::Vapour; + + // Make the fbhp() function. + const auto& controls = baseif_.wellEcl().productionControls(summary_state); + const auto& table = baseif_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); + const double vfp_ref_depth = table.getDatumDepth(); + const double thp_limit = baseif_.getTHPConstraint(summary_state); + const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity()); + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { + assert(rates.size() == 3); + return baseif_.vfpProperties()->getProd() + ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp; + }; + + // Make the flo() function. + auto flo = [&table](const std::vector& rates) { + return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); + }; + + // Find the bhp-point where production becomes nonzero. + double bhp_max = 0.0; + { + auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); }; + double low = controls.bhp_limit; + double high = maxPerfPress + 1.0 * unit::barsa; + double f_low = fflo(low); + double f_high = fflo(high); + deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() + + " low = " + std::to_string(low) + + " high = " + std::to_string(high) + + " f(low) = " + std::to_string(f_low) + + " f(high) = " + std::to_string(f_high)); + int adjustments = 0; + const int max_adjustments = 10; + const double adjust_amount = 5.0 * unit::barsa; + while (f_low * f_high > 0.0 && adjustments < max_adjustments) { + // Same sign, adjust high to see if we can flip it. + high += adjust_amount; + f_high = fflo(high); + ++adjustments; + } + if (f_low * f_high > 0.0) { + if (f_low > 0.0) { + // Even at the BHP limit, we are injecting. + // There will be no solution here, return an + // empty optional. + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE", + "Robust bhp(thp) solve failed due to inoperability for well " + baseif_.name()); + return std::optional(); + } else { + // Still producing, even at high bhp. + assert(f_high < 0.0); + bhp_max = high; + } + } else { + // Bisect to find a bhp point where we produce, but + // not a large amount ('eps' below). + const double eps = 0.1 * std::fabs(table.getFloAxis().front()); + const int maxit = 50; + int it = 0; + while (std::fabs(f_low) > eps && it < maxit) { + const double curr = 0.5*(low + high); + const double f_curr = fflo(curr); + if (f_curr * f_low > 0.0) { + low = curr; + f_low = f_curr; + } else { + high = curr; + f_high = f_curr; + } + ++it; + } + bhp_max = low; + } + deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() + + " low = " + std::to_string(low) + + " high = " + std::to_string(high) + + " f(low) = " + std::to_string(f_low) + + " f(high) = " + std::to_string(f_high) + + " bhp_max = " + std::to_string(bhp_max)); + } + + // Define the equation we want to solve. + auto eq = [&fbhp, &frates](double bhp) { + return fbhp(frates(bhp)) - bhp; + }; + + // Find appropriate brackets for the solution. + double low = controls.bhp_limit; + double high = bhp_max; + { + double eq_high = eq(high); + double eq_low = eq(low); + const double eq_bhplimit = eq_low; + deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + baseif_.name() + + " low = " + std::to_string(low) + + " high = " + std::to_string(high) + + " eq(low) = " + std::to_string(eq_low) + + " eq(high) = " + std::to_string(eq_high)); + if (eq_low * eq_high > 0.0) { + // Failed to bracket the zero. + // If this is due to having two solutions, bisect until bracketed. + double abs_low = std::fabs(eq_low); + double abs_high = std::fabs(eq_high); + int bracket_attempts = 0; + const int max_bracket_attempts = 20; + double interval = high - low; + const double min_interval = 1.0 * unit::barsa; + while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) { + if (abs_high < abs_low) { + low = 0.5 * (low + high); + eq_low = eq(low); + abs_low = std::fabs(eq_low); + } else { + high = 0.5 * (low + high); + eq_high = eq(high); + abs_high = std::fabs(eq_high); + } + ++bracket_attempts; + } + if (eq_low * eq_high > 0.0) { + // Still failed bracketing! + const double limit = 3.0 * unit::barsa; + if (std::min(abs_low, abs_high) < limit) { + // Return the least bad solution if less off than 3 bar. + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", + "Robust bhp(thp) not solved precisely for well " + baseif_.name()); + return abs_low < abs_high ? low : high; + } else { + // Return failure. + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", + "Robust bhp(thp) solve failed due to bracketing failure for well " + baseif_.name()); + return std::optional(); + } + } + } + // We have a bracket! + // Now, see if (bhplimit, low) is a bracket in addition to (low, high). + // If so, that is the bracket we shall use, choosing the solution with the + // highest flow. + if (eq_low * eq_bhplimit <= 0.0) { + high = low; + low = controls.bhp_limit; + } + } + + // Solve for the proper solution in the given interval. + const int max_iteration = 100; + const double bhp_tolerance = 0.01 * unit::barsa; + int iteration = 0; + try { + const double solved_bhp = RegulaFalsiBisection:: + solve(eq, low, high, max_iteration, bhp_tolerance, iteration); + return solved_bhp; + } + catch (...) { + deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", + "Robust bhp(thp) solve failed for well " + baseif_.name()); + return std::optional(); + } +} + +template +bool +MultisegmentWellGeneric:: +frictionalPressureLossConsidered() const +{ + // HF- and HFA needs to consider frictional pressure loss + return (segmentSet().compPressureDrop() != WellSegments::CompPressureDrop::H__); +} + +template +bool +MultisegmentWellGeneric:: +accelerationalPressureLossConsidered() const +{ + return (segmentSet().compPressureDrop() == WellSegments::CompPressureDrop::HFA); +} + +template class MultisegmentWellGeneric; + +} // namespace Opm diff --git a/opm/simulators/wells/MultisegmentWellGeneric.hpp b/opm/simulators/wells/MultisegmentWellGeneric.hpp new file mode 100644 index 000000000..133a25e35 --- /dev/null +++ b/opm/simulators/wells/MultisegmentWellGeneric.hpp @@ -0,0 +1,113 @@ +/* + Copyright 2017 SINTEF Digital, Mathematics and Cybernetics. + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED +#define OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED + +#include + +#include +#include +#include + +namespace Opm +{ + +class DeferredLogger; +class SummaryState; +class WellInterfaceGeneric; +class WellState; + +template +class MultisegmentWellGeneric +{ +protected: + MultisegmentWellGeneric(WellInterfaceGeneric& baseif); + + // scale the segment rates and pressure based on well rates and bhp + void scaleSegmentRatesWithWellRates(WellState& well_state) const; + void scaleSegmentPressuresWithBhp(WellState& well_state) const; + + // get the WellSegments from the well_ecl_ + const WellSegments& segmentSet() const; + + // components of the pressure drop to be included + WellSegments::CompPressureDrop compPressureDrop() const; + + // segment number is an ID of the segment, it is specified in the deck + // get the loation of the segment with a segment number in the segmentSet + int segmentNumberToIndex(const int segment_number) const; + + /// number of segments for this well + int numberOfSegments() const; + + double calculateThpFromBhp(const std::vector& rates, + const double bhp, + const double rho, + DeferredLogger& deferred_logger) const; + + std::optional computeBhpAtThpLimitInj(const std::function(const double)>& frates, + const SummaryState& summary_state, + const double rho, + DeferredLogger& deferred_logger) const; + + std::optional computeBhpAtThpLimitProd(const std::function(const double)>& frates, + const SummaryState& summary_state, + const double maxPerfPress, + const double rho, + DeferredLogger& deferred_logger) const; + + /// Detect oscillation or stagnation based on the residual measure history + void detectOscillations(const std::vector& measure_history, + const int it, + bool& oscillate, + bool& stagnate) const; + + bool accelerationalPressureLossConsidered() const; + bool frictionalPressureLossConsidered() const; + + const WellInterfaceGeneric& baseif_; + + // TODO: trying to use the information from the Well opm-parser as much + // as possible, it will possibly be re-implemented later for efficiency reason. + + // the completions that is related to each segment + // the completions's ids are their index in the vector well_index_, well_cell_ + // This is also assuming the order of the completions in Well is the same with + // the order of the completions in wells. + // it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations + // belonging to this segment + std::vector> segment_perforations_; + + // the inlet segments for each segment. It is for convinience and efficiency reason + std::vector> segment_inlets_; + + std::vector segment_depth_diffs_; + + // depth difference between the segment and the perforation + // or in another way, the depth difference between the perforation and + // the segment the perforation belongs to + std::vector perforation_segment_depth_diffs_; +}; + +} + +#endif // OPM_MULTISEGMENTWELL_GENERIC_HEADER_INCLUDED diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 843827644..2cf0b286d 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -45,20 +45,17 @@ namespace Opm const int first_perf_index, const std::vector& perf_data) : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, first_perf_index, perf_data) - , segment_perforations_(numberOfSegments()) - , segment_inlets_(numberOfSegments()) + , MultisegmentWellGeneric(static_cast(*this)) , cell_perforation_depth_diffs_(number_of_perforations_, 0.0) , cell_perforation_pressure_diffs_(number_of_perforations_, 0.0) - , perforation_segment_depth_diffs_(number_of_perforations_, 0.0) - , segment_fluid_initial_(numberOfSegments(), std::vector(num_components_, 0.0)) - , segment_densities_(numberOfSegments(), 0.0) - , segment_viscosities_(numberOfSegments(), 0.0) - , segment_mass_rates_(numberOfSegments(), 0.0) - , segment_depth_diffs_(numberOfSegments(), 0.0) - , upwinding_segments_(numberOfSegments(), 0) - , segment_phase_fractions_(numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? - , segment_phase_viscosities_(numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? - , segment_phase_densities_(numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? + , segment_fluid_initial_(this->numberOfSegments(), std::vector(num_components_, 0.0)) + , segment_densities_(this->numberOfSegments(), 0.0) + , segment_viscosities_(this->numberOfSegments(), 0.0) + , segment_mass_rates_(this->numberOfSegments(), 0.0) + , upwinding_segments_(this->numberOfSegments(), 0) + , segment_phase_fractions_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? + , segment_phase_viscosities_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? + , segment_phase_densities_(this->numberOfSegments(), std::vector(num_components_, 0.0)) // number of phase here? { // not handling solvent or polymer for now with multisegment well if constexpr (has_solvent) { @@ -80,51 +77,6 @@ namespace Opm if constexpr (Base::has_brine) { OPM_THROW(std::runtime_error, "brine is not supported by multisegment well yet"); } - // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. - // for other facilities needed but not available from parser, we need to process them here - - // initialize the segment_perforations_ and update perforation_segment_depth_diffs_ - const WellConnections& completion_set = well_ecl_.getConnections(); - // index of the perforation within wells struct - // there might be some perforations not active, which causes the number of the perforations in - // well_ecl_ and wells struct different - // the current implementation is a temporary solution for now, it should be corrected from the parser - // side - int i_perf_wells = 0; - perf_depth_.resize(number_of_perforations_, 0.); - for (size_t perf = 0; perf < completion_set.size(); ++perf) { - const Connection& connection = completion_set.get(perf); - if (connection.state() == Connection::State::OPEN) { - const int segment_index = segmentNumberToIndex(connection.segment()); - segment_perforations_[segment_index].push_back(i_perf_wells); - perf_depth_[i_perf_wells] = connection.depth(); - const double segment_depth = segmentSet()[segment_index].depth(); - perforation_segment_depth_diffs_[i_perf_wells] = perf_depth_[i_perf_wells] - segment_depth; - i_perf_wells++; - } - } - - // initialize the segment_inlets_ - for (int seg = 0; seg < numberOfSegments(); ++seg) { - const Segment& segment = segmentSet()[seg]; - const int segment_number = segment.segmentNumber(); - const int outlet_segment_number = segment.outletSegment(); - if (outlet_segment_number > 0) { - const int segment_index = segmentNumberToIndex(segment_number); - const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number); - segment_inlets_[outlet_segment_index].push_back(segment_index); - } - } - - // calculating the depth difference between the segment and its oulet_segments - // for the top segment, we will make its zero unless we find other purpose to use this value - for (int seg = 1; seg < numberOfSegments(); ++seg) { - const double segment_depth = segmentSet()[seg].depth(); - const int outlet_segment_number = segmentSet()[seg].outletSegment(); - const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)]; - const double outlet_depth = outlet_segment.depth(); - segment_depth_diffs_[seg] = segment_depth - outlet_depth; - } } @@ -181,24 +133,24 @@ namespace Opm // calculatiing the NNZ for duneD_ // NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets) { - int nnz_d = numberOfSegments(); - for (const std::vector& inlets : segment_inlets_) { + int nnz_d = this->numberOfSegments(); + for (const std::vector& inlets : this->segment_inlets_) { nnz_d += 2 * inlets.size(); } - duneD_.setSize(numberOfSegments(), numberOfSegments(), nnz_d); + duneD_.setSize(this->numberOfSegments(), this->numberOfSegments(), nnz_d); } - duneB_.setSize(numberOfSegments(), num_cells, number_of_perforations_); - duneC_.setSize(numberOfSegments(), num_cells, number_of_perforations_); + duneB_.setSize(this->numberOfSegments(), num_cells, number_of_perforations_); + duneC_.setSize(this->numberOfSegments(), num_cells, number_of_perforations_); // we need to add the off diagonal ones for (auto row = duneD_.createbegin(), end = duneD_.createend(); row != end; ++row) { // the number of the row corrspnds to the segment now const int seg = row.index(); // adding the item related to outlet relation - const Segment& segment = segmentSet()[seg]; + const Segment& segment = this->segmentSet()[seg]; const int outlet_segment_number = segment.outletSegment(); if (outlet_segment_number > 0) { // if there is a outlet_segment - const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number); + const int outlet_segment_index = this->segmentNumberToIndex(outlet_segment_number); row.insert(outlet_segment_index); } @@ -206,7 +158,7 @@ namespace Opm row.insert(seg); // insert the item related to its inlets - for (const int& inlet : segment_inlets_[seg]) { + for (const int& inlet : this->segment_inlets_[seg]) { row.insert(inlet); } } @@ -214,7 +166,7 @@ namespace Opm // make the C matrix for (auto row = duneC_.createbegin(), end = duneC_.createend(); row != end; ++row) { // the number of the row corresponds to the segment number now. - for (const int& perf : segment_perforations_[row.index()]) { + for (const int& perf : this->segment_perforations_[row.index()]) { const int cell_idx = well_cells_[perf]; row.insert(cell_idx); } @@ -223,16 +175,16 @@ namespace Opm // make the B^T matrix for (auto row = duneB_.createbegin(), end = duneB_.createend(); row != end; ++row) { // the number of the row corresponds to the segment number now. - for (const int& perf : segment_perforations_[row.index()]) { + for (const int& perf : this->segment_perforations_[row.index()]) { const int cell_idx = well_cells_[perf]; row.insert(cell_idx); } } - resWell_.resize( numberOfSegments() ); + resWell_.resize( this->numberOfSegments() ); - primary_variables_.resize(numberOfSegments()); - primary_variables_evaluation_.resize(numberOfSegments()); + primary_variables_.resize(this->numberOfSegments()); + primary_variables_evaluation_.resize(this->numberOfSegments()); } @@ -244,7 +196,7 @@ namespace Opm MultisegmentWell:: initPrimaryVariablesEvaluation() const { - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { primary_variables_evaluation_[seg][eq_idx] = 0.0; primary_variables_evaluation_[seg][eq_idx].setValue(primary_variables_[seg][eq_idx]); @@ -265,61 +217,14 @@ namespace Opm Base::updateWellStateWithTarget(ebos_simulator, well_state, deferred_logger); // scale segment rates based on the wellRates // and segment pressure based on bhp - scaleSegmentRatesWithWellRates(well_state); - scaleSegmentPressuresWithBhp(well_state); + this->scaleSegmentRatesWithWellRates(well_state); + this->scaleSegmentPressuresWithBhp(well_state); } - template - void - MultisegmentWell:: - scaleSegmentRatesWithWellRates(WellState& well_state) const - { - auto& segments = well_state.segments(this->index_of_well_); - auto& segment_rates = segments.rates; - for (int phase = 0; phase < number_of_phases_; ++phase) { - const double unscaled_top_seg_rate = segment_rates[phase]; - const double well_phase_rate = well_state.wellRates(index_of_well_)[phase]; - if (std::abs(unscaled_top_seg_rate) > 1e-12) - { - for (int seg = 0; seg < numberOfSegments(); ++seg) { - segment_rates[this->number_of_phases_*seg + phase] *= well_phase_rate/unscaled_top_seg_rate; - } - } else { - // for newly opened wells, the unscaled rate top segment rate is zero - // and we need to initialize the segment rates differently - double sumTw = 0; - for (int perf = 0; perf < number_of_perforations_; ++perf) { - sumTw += well_index_[perf]; - } - - std::vector perforation_rates(number_of_phases_ * number_of_perforations_,0.0); - const double perf_phaserate_scaled = well_state.wellRates(index_of_well_)[phase] / sumTw; - for (int perf = 0; perf < number_of_perforations_; ++perf) { - perforation_rates[number_of_phases_ * perf + phase] = well_index_[perf] * perf_phaserate_scaled; - } - - std::vector rates; - WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, number_of_phases_, 0, rates); - std::copy(rates.begin(), rates.end(), segment_rates.begin()); - } - } - } - - template - void - MultisegmentWell:: - scaleSegmentPressuresWithBhp(WellState& well_state) const - { - auto& segments = well_state.segments(this->index_of_well_); - auto bhp = well_state.bhp(this->index_of_well_); - segments.scale_pressure(bhp); - } - - template ConvergenceReport MultisegmentWell:: @@ -328,8 +233,8 @@ namespace Opm assert(int(B_avg.size()) == num_components_); // checking if any residual is NaN or too large. The two large one is only handled for the well flux - std::vector> abs_residual(numberOfSegments(), std::vector(numWellEq, 0.0)); - for (int seg = 0; seg < numberOfSegments(); ++seg) { + std::vector> abs_residual(this->numberOfSegments(), std::vector(numWellEq, 0.0)); + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { abs_residual[seg][eq_idx] = std::abs(resWell_[seg][eq_idx]); } @@ -340,7 +245,7 @@ namespace Opm ConvergenceReport report; // TODO: the following is a little complicated, maybe can be simplified in some way? for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { if (eq_idx < num_components_) { // phase or component mass equations const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx]; if (flux_residual > maximum_residual[eq_idx]) { @@ -720,7 +625,7 @@ namespace Opm const auto& segment_pressure = segments.pressure; const PhaseUsage& pu = phaseUsage(); - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { // calculate the total rate for each segment double total_seg_rate = 0.0; // the segment pressure @@ -875,7 +780,7 @@ namespace Opm MultisegmentWell:: computeInitialSegmentFluids(const Simulator& ebos_simulator) { - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { // TODO: trying to reduce the times for the surfaceVolumeFraction calculation const double surface_volume = getSegmentSurfaceVolume(ebos_simulator, seg).value(); for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { @@ -902,7 +807,7 @@ namespace Opm const double max_pressure_change = param_.max_pressure_change_ms_wells_; const std::vector > old_primary_variables = primary_variables_; - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { if (has_wfrac_variable) { const int sign = dwells[seg][WFrac] > 0. ? 1 : -1; const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit); @@ -1072,56 +977,6 @@ namespace Opm - template - const WellSegments& - MultisegmentWell:: - segmentSet() const - { - return well_ecl_.getSegments(); - } - - - - - - template - int - MultisegmentWell:: - numberOfSegments() const - { - return segmentSet().size(); - } - - - - - - template - WellSegments::CompPressureDrop - MultisegmentWell:: - compPressureDrop() const - { - return segmentSet().compPressureDrop(); - } - - - - - - - - template - int - MultisegmentWell:: - segmentNumberToIndex(const int segment_number) const - { - return segmentSet().segmentNumberToIndex(segment_number); - } - - - - - template typename MultisegmentWell::EvalWell MultisegmentWell:: @@ -1236,7 +1091,7 @@ namespace Opm } // pressure difference between the segment and the perforation - const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf]; + const EvalWell perf_seg_press_diff = gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf]; // pressure difference between the perforation and the grid cell const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf]; @@ -1410,7 +1265,7 @@ namespace Opm surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index ); } - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { // the compostion of the components inside wellbore under surface condition std::vector mix_s(num_components_, 0.0); for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { @@ -1839,52 +1694,12 @@ namespace Opm const double bhp = well_state.bhp(index_of_well_); - well_state.update_thp(index_of_well_, calculateThpFromBhp(rates, bhp, deferred_logger)); + well_state.update_thp(index_of_well_, this->calculateThpFromBhp(rates, bhp, getRefDensity(), deferred_logger)); } - template - double - MultisegmentWell:: - calculateThpFromBhp(const std::vector& rates, - const double bhp, - DeferredLogger& deferred_logger) const - { - assert(int(rates.size()) == 3); // the vfp related only supports three phases now. - - const double aqua = rates[Water]; - const double liquid = rates[Oil]; - const double vapour = rates[Gas]; - - // pick the density in the top segment - const double rho = getRefDensity(); - - double thp = 0.0; - if (this->isInjector()) { - const int table_id = well_ecl_.vfp_table_number(); - const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id).getDatumDepth(); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - - thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp); - } - else if (this->isProducer()) { - const int table_id = well_ecl_.vfp_table_number(); - const double alq = well_ecl_.alq_value(); - const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id).getDatumDepth(); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - - thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq); - } - else { - OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger); - } - - return thp; - } - - template double MultisegmentWell:: @@ -1912,7 +1727,7 @@ namespace Opm segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value(); pressure_equation -= hydro_pressure_drop; - if (frictionalPressureLossConsidered()) { + if (this->frictionalPressureLossConsidered()) { const auto friction_pressure_drop = getFrictionPressureLoss(seg); pressure_equation -= friction_pressure_drop; segments.pressure_drop_friction[seg] = friction_pressure_drop.value(); @@ -1930,7 +1745,7 @@ namespace Opm } // contribution from the outlet segment - const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); resWell_[seg][SPres] -= outlet_pressure.value(); @@ -1938,7 +1753,7 @@ namespace Opm duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq); } - if (accelerationalPressureLossConsidered()) { + if (this->accelerationalPressureLossConsidered()) { handleAccelerationPressureLoss(seg, well_state); } } @@ -1952,7 +1767,7 @@ namespace Opm MultisegmentWell:: getHydroPressureLoss(const int seg) const { - return segment_densities_[seg] * gravity_ * segment_depth_diffs_[seg]; + return segment_densities_[seg] * gravity_ * this->segment_depth_diffs_[seg]; } @@ -1975,12 +1790,12 @@ namespace Opm density.clearDerivatives(); visc.clearDerivatives(); } - const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); - const double length = segmentSet()[seg].totalLength() - segmentSet()[outlet_segment_index].totalLength(); + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); + const double length = this->segmentSet()[seg].totalLength() - this->segmentSet()[outlet_segment_index].totalLength(); assert(length > 0.); - const double roughness = segmentSet()[seg].roughness(); - const double area = segmentSet()[seg].crossArea(); - const double diameter = segmentSet()[seg].internalDiameter(); + const double roughness = this->segmentSet()[seg].roughness(); + const double area = this->segmentSet()[seg].crossArea(); + const double diameter = this->segmentSet()[seg].internalDiameter(); const double sign = mass_rate < 0. ? 1.0 : - 1.0; @@ -1996,7 +1811,7 @@ namespace Opm MultisegmentWell:: handleAccelerationPressureLoss(const int seg, WellState& well_state) const { - const double area = segmentSet()[seg].crossArea(); + const double area = this->segmentSet()[seg].crossArea(); const EvalWell mass_rate = segment_mass_rates_[seg]; const int seg_upwind = upwinding_segments_[seg]; EvalWell density = segment_densities_[seg_upwind]; @@ -2009,10 +1824,10 @@ namespace Opm EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density); // handling the velocity head of intlet segments - for (const int inlet : segment_inlets_[seg]) { + for (const int inlet : this->segment_inlets_[seg]) { const int seg_upwind_inlet = upwinding_segments_[inlet]; - const double inlet_area = segmentSet()[inlet].crossArea(); - EvalWell inlet_density = segment_densities_[seg_upwind_inlet]; + const double inlet_area = this->segmentSet()[inlet].crossArea(); + EvalWell inlet_density = this->segment_densities_[seg_upwind_inlet]; // WARNING // We disregard the derivatives from the upwind density to make sure derivatives // wrt. to different segments dont get mixed. @@ -2142,7 +1957,7 @@ namespace Opm std::vector well_rates_bhp_limit; computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger); - const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger); + const double thp = this->calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, getRefDensity(), deferred_logger); const double thp_limit = this->getTHPConstraint(summaryState); @@ -2182,16 +1997,16 @@ namespace Opm std::fill(ipr_a_.begin(), ipr_a_.end(), 0.); std::fill(ipr_b_.begin(), ipr_b_.end(), 0.); - const int nseg = numberOfSegments(); + const int nseg = this->numberOfSegments(); double seg_bhp_press_diff = 0; double ref_depth = ref_depth_; for (int seg = 0; seg < nseg; ++seg) { // calculating the perforation rate for each perforation that belongs to this segment - const double segment_depth = segmentSet()[seg].depth(); + const double segment_depth = this->segmentSet()[seg].depth(); const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth, segment_depth, segment_densities_[seg].value(), gravity_); ref_depth = segment_depth; seg_bhp_press_diff += dp; - for (const int perf : segment_perforations_[seg]) { + for (const int perf : this->segment_perforations_[seg]) { //std::vector mob(num_components_, {numWellEq_ + numEq, 0.0}); std::vector mob(num_components_, 0.0); @@ -2203,7 +2018,7 @@ namespace Opm const auto& fs = int_quantities.fluidState(); // the pressure of the reservoir grid block the well connection is in // pressure difference between the segment and the perforation - const double perf_seg_press_diff = gravity_ * segment_densities_[seg].value() * perforation_segment_depth_diffs_[perf]; + const double perf_seg_press_diff = gravity_ * segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf]; // pressure difference between the perforation and the grid cell const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf]; const double pressure_cell = fs.pressure(FluidSystem::oilPhaseIdx).value(); @@ -2325,7 +2140,7 @@ namespace Opm auto& segments = well_state.segments(this->index_of_well_); auto& segment_rates = segments.rates; auto& segment_pressure = segments.pressure; - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { std::vector fractions(number_of_phases_, 0.0); fractions[oil_pos] = 1.0; @@ -2375,31 +2190,6 @@ namespace Opm - template - bool - MultisegmentWell:: - frictionalPressureLossConsidered() const - { - // HF- and HFA needs to consider frictional pressure loss - return (segmentSet().compPressureDrop() != WellSegments::CompPressureDrop::H__); - } - - - - - - template - bool - MultisegmentWell:: - accelerationalPressureLossConsidered() const - { - return (segmentSet().compPressureDrop() == WellSegments::CompPressureDrop::HFA); - } - - - - - template bool MultisegmentWell:: @@ -2446,7 +2236,7 @@ namespace Opm bool is_oscillate = false; bool is_stagnate = false; - detectOscillations(measure_history, it, is_oscillate, is_stagnate); + this->detectOscillations(measure_history, it, is_oscillate, is_stagnate); // TODO: maybe we should have more sophiscated strategy to recover the relaxation factor, // for example, to recover it to be bigger @@ -2558,7 +2348,7 @@ namespace Opm const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); - const int nseg = numberOfSegments(); + const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { // calculating the accumulation term @@ -2603,7 +2393,7 @@ namespace Opm // considering the contributions from the inlet segments { - for (const int inlet : segment_inlets_[seg]) { + for (const int inlet : this->segment_inlets_[seg]) { for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) { const EvalWell inlet_rate = getSegmentRateUpwinding(inlet, comp_idx) * well_efficiency_factor_; @@ -2627,7 +2417,7 @@ namespace Opm const EvalWell seg_pressure = getSegmentPressure(seg); auto& perf_rates = well_state.perfPhaseRates(this->index_of_well_); auto& perf_press_state = well_state.perfPress(this->index_of_well_); - for (const int perf : segment_perforations_[seg]) { + for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); std::vector mob(num_components_, 0.0); @@ -2699,7 +2489,7 @@ namespace Opm assemblePressureEq(const int seg, const UnitSystem& unit_system, WellState& well_state, DeferredLogger& deferred_logger) const { - switch(segmentSet()[seg].segmentType()) { + switch(this->segmentSet()[seg].segmentType()) { case Segment::SegmentType::SICD : case Segment::SegmentType::AICD : case Segment::SegmentType::VALVE : { @@ -2728,18 +2518,18 @@ namespace Opm allDrawDownWrongDirection(const Simulator& ebos_simulator) const { bool all_drawdown_wrong_direction = true; - const int nseg = numberOfSegments(); + const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { const EvalWell segment_pressure = getSegmentPressure(seg); - for (const int perf : segment_perforations_[seg]) { + for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = intQuants.fluidState(); // pressure difference between the segment and the perforation - const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * perforation_segment_depth_diffs_[perf]; + const EvalWell perf_seg_press_diff = gravity_ * segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf]; // pressure difference between the perforation and the grid cell const double cell_perf_press_diff = cell_perforation_pressure_diffs_[perf]; @@ -2769,7 +2559,7 @@ namespace Opm template void MultisegmentWell:: - updateWaterThroughput(const double dt OPM_UNUSED, WellState& well_state OPM_UNUSED) const + updateWaterThroughput(const double /*dt*/, WellState& /*well_state*/) const { } @@ -2903,7 +2693,7 @@ namespace Opm } // We increase the segment volume with a factor 10 to stabilize the system. - const double volume = segmentSet()[seg_idx].volume(); + const double volume = this->segmentSet()[seg_idx].volume(); return volume / vol_ratio; } @@ -2921,7 +2711,7 @@ namespace Opm assert(int(B_avg.size() ) == num_components_); std::vector residuals(numWellEq + 1, 0.0); - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) { double residual = 0.; if (eq_idx < num_components_) { @@ -2958,37 +2748,6 @@ namespace Opm - /// Detect oscillation or stagnation based on the residual measure history - template - void - MultisegmentWell:: - detectOscillations(const std::vector& measure_history, - const int it, bool& oscillate, bool& stagnate) const - { - if ( it < 2 ) { - oscillate = false; - stagnate = false; - return; - } - - stagnate = true; - const double F0 = measure_history[it]; - const double F1 = measure_history[it - 1]; - const double F2 = measure_history[it - 2]; - const double d1 = std::abs((F0 - F2) / F0); - const double d2 = std::abs((F0 - F1) / F0); - - const double oscillaton_rel_tol = 0.2; - oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2); - - const double stagnation_rel_tol = 1.e-2; - stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol; - } - - - - - template double MultisegmentWell:: @@ -3184,7 +2943,7 @@ namespace Opm MultisegmentWell:: updateUpwindingSegments() { - for (int seg = 0; seg < numberOfSegments(); ++seg) { + for (int seg = 0; seg < this->numberOfSegments(); ++seg) { // special treatment is needed for segment 0 if (seg == 0) { // we are not supposed to have injecting producers and producing injectors @@ -3198,7 +2957,7 @@ namespace Opm if (primary_variables_evaluation_[seg][GTotal] <= 0.) { upwinding_segments_[seg] = seg; } else { - const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); upwinding_segments_[seg] = outlet_segment_index; } } @@ -3224,7 +2983,7 @@ namespace Opm EvalWell pressure_equation = getSegmentPressure(seg); EvalWell icd_pressure_drop; - switch(segmentSet()[seg].segmentType()) { + switch(this->segmentSet()[seg].segmentType()) { case Segment::SegmentType::SICD : icd_pressure_drop = pressureDropSpiralICD(seg); break; @@ -3235,7 +2994,7 @@ namespace Opm icd_pressure_drop = pressureDropValve(seg); break; default: { - OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(segmentSet()[seg].segmentNumber()) + OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(this->segmentSet()[seg].segmentNumber()) + " for well " + name() + " is not of ICD type", deferred_logger); } } @@ -3254,7 +3013,7 @@ namespace Opm } // contribution from the outlet segment - const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment()); + const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment()); const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index); resWell_[seg][SPres] -= outlet_pressure.value(); @@ -3274,40 +3033,6 @@ namespace Opm const SummaryState& summary_state, DeferredLogger& deferred_logger) const { - // Given a VFP function returning bhp as a function of phase - // rates and thp: - // fbhp(rates, thp), - // a function extracting the particular flow rate used for VFP - // lookups: - // flo(rates) - // and the inflow function (assuming the reservoir is fixed): - // frates(bhp) - // we want to solve the equation: - // fbhp(frates(bhp, thplimit)) - bhp = 0 - // for bhp. - // - // This may result in 0, 1 or 2 solutions. If two solutions, - // the one corresponding to the lowest bhp (and therefore - // highest rate) should be returned. - - // Make the fbhp() function. - const auto& controls = well_ecl_.productionControls(summary_state); - const auto& table = vfp_properties_->getProd()->getTable(controls.vfp_table_number); - const double vfp_ref_depth = table.getDatumDepth(); - const double rho = getRefDensity(); // Use the density at the top perforation. - const double thp_limit = this->getTHPConstraint(summary_state); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { - assert(rates.size() == 3); - return this->vfp_properties_->getProd() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp; - }; - - // Make the flo() function. - auto flo = [&table](const std::vector& rates) { - return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); - }; - // Make the frates() function. auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) { // Not solving the well equations here, which means we are @@ -3320,147 +3045,12 @@ namespace Opm return rates; }; - // Find the bhp-point where production becomes nonzero. - double bhp_max = 0.0; - { - auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); }; - double low = controls.bhp_limit; - double high = maxPerfPress(ebos_simulator) + 1.0 * unit::barsa; - double f_low = fflo(low); - double f_high = fflo(high); - deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() + - " low = " + std::to_string(low) + - " high = " + std::to_string(high) + - " f(low) = " + std::to_string(f_low) + - " f(high) = " + std::to_string(f_high)); - int adjustments = 0; - const int max_adjustments = 10; - const double adjust_amount = 5.0 * unit::barsa; - while (f_low * f_high > 0.0 && adjustments < max_adjustments) { - // Same sign, adjust high to see if we can flip it. - high += adjust_amount; - f_high = fflo(high); - ++adjustments; - } - if (f_low * f_high > 0.0) { - if (f_low > 0.0) { - // Even at the BHP limit, we are injecting. - // There will be no solution here, return an - // empty optional. - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_INOPERABLE", - "Robust bhp(thp) solve failed due to inoperability for well " + name()); - return std::optional(); - } else { - // Still producing, even at high bhp. - assert(f_high < 0.0); - bhp_max = high; - } - } else { - // Bisect to find a bhp point where we produce, but - // not a large amount ('eps' below). - const double eps = 0.1 * std::fabs(table.getFloAxis().front()); - const int maxit = 50; - int it = 0; - while (std::fabs(f_low) > eps && it < maxit) { - const double curr = 0.5*(low + high); - const double f_curr = fflo(curr); - if (f_curr * f_low > 0.0) { - low = curr; - f_low = f_curr; - } else { - high = curr; - f_high = f_curr; - } - ++it; - } - bhp_max = low; - } - deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() + - " low = " + std::to_string(low) + - " high = " + std::to_string(high) + - " f(low) = " + std::to_string(f_low) + - " f(high) = " + std::to_string(f_high) + - " bhp_max = " + std::to_string(bhp_max)); - } - - // Define the equation we want to solve. - auto eq = [&fbhp, &frates](double bhp) { - return fbhp(frates(bhp)) - bhp; - }; - - // Find appropriate brackets for the solution. - double low = controls.bhp_limit; - double high = bhp_max; - { - double eq_high = eq(high); - double eq_low = eq(low); - const double eq_bhplimit = eq_low; - deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + name() + - " low = " + std::to_string(low) + - " high = " + std::to_string(high) + - " eq(low) = " + std::to_string(eq_low) + - " eq(high) = " + std::to_string(eq_high)); - if (eq_low * eq_high > 0.0) { - // Failed to bracket the zero. - // If this is due to having two solutions, bisect until bracketed. - double abs_low = std::fabs(eq_low); - double abs_high = std::fabs(eq_high); - int bracket_attempts = 0; - const int max_bracket_attempts = 20; - double interval = high - low; - const double min_interval = 1.0 * unit::barsa; - while (eq_low * eq_high > 0.0 && bracket_attempts < max_bracket_attempts && interval > min_interval) { - if (abs_high < abs_low) { - low = 0.5 * (low + high); - eq_low = eq(low); - abs_low = std::fabs(eq_low); - } else { - high = 0.5 * (low + high); - eq_high = eq(high); - abs_high = std::fabs(eq_high); - } - ++bracket_attempts; - } - if (eq_low * eq_high > 0.0) { - // Still failed bracketing! - const double limit = 3.0 * unit::barsa; - if (std::min(abs_low, abs_high) < limit) { - // Return the least bad solution if less off than 3 bar. - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", - "Robust bhp(thp) not solved precisely for well " + name()); - return abs_low < abs_high ? low : high; - } else { - // Return failure. - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", - "Robust bhp(thp) solve failed due to bracketing failure for well " + name()); - return std::optional(); - } - } - } - // We have a bracket! - // Now, see if (bhplimit, low) is a bracket in addition to (low, high). - // If so, that is the bracket we shall use, choosing the solution with the - // highest flow. - if (eq_low * eq_bhplimit <= 0.0) { - high = low; - low = controls.bhp_limit; - } - } - - // Solve for the proper solution in the given interval. - const int max_iteration = 100; - const double bhp_tolerance = 0.01 * unit::barsa; - int iteration = 0; - try { - const double solved_bhp = RegulaFalsiBisection:: - solve(eq, low, high, max_iteration, bhp_tolerance, iteration); - return solved_bhp; - } - catch (...) { - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", - "Robust bhp(thp) solve failed for well " + name()); - return std::optional(); - } + return this->MultisegmentWellGeneric:: + computeBhpAtThpLimitProd(frates, + summary_state, + maxPerfPress(ebos_simulator), + getRefDensity(), + deferred_logger); } @@ -3473,63 +3063,6 @@ namespace Opm const SummaryState& summary_state, DeferredLogger& deferred_logger) const { - // Given a VFP function returning bhp as a function of phase - // rates and thp: - // fbhp(rates, thp), - // a function extracting the particular flow rate used for VFP - // lookups: - // flo(rates) - // and the inflow function (assuming the reservoir is fixed): - // frates(bhp) - // we want to solve the equation: - // fbhp(frates(bhp, thplimit)) - bhp = 0 - // for bhp. - // - // This may result in 0, 1 or 2 solutions. If two solutions, - // the one corresponding to the lowest bhp (and therefore - // highest rate) is returned. - // - // In order to detect these situations, we will find piecewise - // linear approximations both to the inverse of the frates - // function and to the fbhp function. - // - // We first take the FLO sample points of the VFP curve, and - // find the corresponding bhp values by solving the equation: - // flo(frates(bhp)) - flo_sample = 0 - // for bhp, for each flo_sample. The resulting (flo_sample, - // bhp_sample) values give a piecewise linear approximation to - // the true inverse inflow function, at the same flo values as - // the VFP data. - // - // Then we extract a piecewise linear approximation from the - // multilinear fbhp() by evaluating it at the flo_sample - // points, with fractions given by the frates(bhp_sample) - // values. - // - // When we have both piecewise linear curves defined on the - // same flo_sample points, it is easy to distinguish between - // the 0, 1 or 2 solution cases, and obtain the right interval - // in which to solve for the solution we want (with highest - // flow in case of 2 solutions). - - // Make the fbhp() function. - const auto& controls = well_ecl_.injectionControls(summary_state); - const auto& table = vfp_properties_->getInj()->getTable(controls.vfp_table_number); - const double vfp_ref_depth = table.getDatumDepth(); - const double rho = getRefDensity(); // Use the density at the top perforation. - const double thp_limit = this->getTHPConstraint(summary_state); - const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_); - auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { - assert(rates.size() == 3); - return this->vfp_properties_->getInj() - ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit) - dp; - }; - - // Make the flo() function. - auto flo = [&table](const std::vector& rates) { - return detail::getFlo(table, rates[Water], rates[Oil], rates[Gas]); - }; - // Make the frates() function. auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) { // Not solving the well equations here, which means we are @@ -3542,126 +3075,8 @@ namespace Opm return rates; }; - // Get the flo samples, add extra samples at low rates and bhp - // limit point if necessary. - std::vector flo_samples = table.getFloAxis(); - if (flo_samples[0] > 0.0) { - const double f0 = flo_samples[0]; - flo_samples.insert(flo_samples.begin(), { f0/20.0, f0/10.0, f0/5.0, f0/2.0 }); - } - const double flo_bhp_limit = flo(frates(controls.bhp_limit)); - if (flo_samples.back() < flo_bhp_limit) { - flo_samples.push_back(flo_bhp_limit); - } - - // Find bhp values for inflow relation corresponding to flo samples. - std::vector bhp_samples; - for (double flo_sample : flo_samples) { - if (flo_sample > flo_bhp_limit) { - // We would have to go over the bhp limit to obtain a - // flow of this magnitude. We associate all such flows - // with simply the bhp limit. The first one - // encountered is considered valid, the rest not. They - // are therefore skipped. - bhp_samples.push_back(controls.bhp_limit); - break; - } - auto eq = [&flo, &frates, flo_sample](double bhp) { - return flo(frates(bhp)) - flo_sample; - }; - // TODO: replace hardcoded low/high limits. - const double low = 10.0 * unit::barsa; - const double high = 800.0 * unit::barsa; - const int max_iteration = 100; - const double flo_tolerance = 0.05 * std::fabs(flo_samples.back()); - int iteration = 0; - try { - const double solved_bhp = RegulaFalsiBisection:: - solve(eq, low, high, max_iteration, flo_tolerance, iteration); - bhp_samples.push_back(solved_bhp); - } - catch (...) { - // Use previous value (or max value if at start) if we failed. - bhp_samples.push_back(bhp_samples.empty() ? low : bhp_samples.back()); - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_EXTRACT_SAMPLES", - "Robust bhp(thp) solve failed extracting bhp values at flo samples for well " + name()); - } - } - - // Find bhp values for VFP relation corresponding to flo samples. - const int num_samples = bhp_samples.size(); // Note that this can be smaller than flo_samples.size() - std::vector fbhp_samples(num_samples); - for (int ii = 0; ii < num_samples; ++ii) { - fbhp_samples[ii] = fbhp(frates(bhp_samples[ii])); - } -// #define EXTRA_THP_DEBUGGING -#ifdef EXTRA_THP_DEBUGGING - std::string dbgmsg; - dbgmsg += "flo: "; - for (int ii = 0; ii < num_samples; ++ii) { - dbgmsg += " " + std::to_string(flo_samples[ii]); - } - dbgmsg += "\nbhp: "; - for (int ii = 0; ii < num_samples; ++ii) { - dbgmsg += " " + std::to_string(bhp_samples[ii]); - } - dbgmsg += "\nfbhp: "; - for (int ii = 0; ii < num_samples; ++ii) { - dbgmsg += " " + std::to_string(fbhp_samples[ii]); - } - OpmLog::debug(dbgmsg); -#endif // EXTRA_THP_DEBUGGING - - // Look for sign changes for the (fbhp_samples - bhp_samples) piecewise linear curve. - // We only look at the valid - int sign_change_index = -1; - for (int ii = 0; ii < num_samples - 1; ++ii) { - const double curr = fbhp_samples[ii] - bhp_samples[ii]; - const double next = fbhp_samples[ii + 1] - bhp_samples[ii + 1]; - if (curr * next < 0.0) { - // Sign change in the [ii, ii + 1] interval. - sign_change_index = ii; // May overwrite, thereby choosing the highest-flo solution. - } - } - - // Handle the no solution case. - if (sign_change_index == -1) { - return std::optional(); - } - - // Solve for the proper solution in the given interval. - auto eq = [&fbhp, &frates](double bhp) { - return fbhp(frates(bhp)) - bhp; - }; - // TODO: replace hardcoded low/high limits. - const double low = bhp_samples[sign_change_index + 1]; - const double high = bhp_samples[sign_change_index]; - const int max_iteration = 100; - const double bhp_tolerance = 0.01 * unit::barsa; - int iteration = 0; - if (low == high) { - // We are in the high flow regime where the bhp_samples - // are all equal to the bhp_limit. - assert(low == controls.bhp_limit); - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", - "Robust bhp(thp) solve failed for well " + name()); - return std::optional(); - } - try { - const double solved_bhp = RegulaFalsiBisection:: - solve(eq, low, high, max_iteration, bhp_tolerance, iteration); -#ifdef EXTRA_THP_DEBUGGING - OpmLog::debug("***** " + name() + " solved_bhp = " + std::to_string(solved_bhp) - + " flo_bhp_limit = " + std::to_string(flo_bhp_limit)); -#endif // EXTRA_THP_DEBUGGING - return solved_bhp; - } - catch (...) { - deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE", - "Robust bhp(thp) solve failed for well " + name()); - return std::optional(); - } - + return this->MultisegmentWellGeneric:: + computeBhpAtThpLimitInj(frates, summary_state, getRefDensity(), deferred_logger); } @@ -3674,9 +3089,9 @@ namespace Opm maxPerfPress(const Simulator& ebos_simulator) const { double max_pressure = 0.0; - const int nseg = numberOfSegments(); + const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { - for (const int perf : segment_perforations_[seg]) { + for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); const auto& fs = int_quants.fluidState(); @@ -3696,7 +3111,7 @@ namespace Opm MultisegmentWell:: pressureDropSpiralICD(const int seg) const { - const SICD& sicd = segmentSet()[seg].spiralICD(); + const SICD& sicd = this->segmentSet()[seg].spiralICD(); const int seg_upwind = upwinding_segments_[seg]; const std::vector& phase_fractions = segment_phase_fractions_[seg_upwind]; @@ -3862,7 +3277,7 @@ namespace Opm MultisegmentWell:: pressureDropValve(const int seg) const { - const Valve& valve = segmentSet()[seg].valve(); + const Valve& valve = this->segmentSet()[seg].valve(); const EvalWell& mass_rate = segment_mass_rates_[seg]; const int seg_upwind = upwinding_segments_[seg]; @@ -3906,11 +3321,11 @@ namespace Opm // Calculate the rates that follow from the current primary variables. std::vector well_q_s(num_components_, 0.0); const bool allow_cf = getAllowCrossFlow() || openCrossFlowAvoidSingularity(ebosSimulator); - const int nseg = numberOfSegments(); + const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { // calculating the perforation rate for each perforation that belongs to this segment const EvalWell seg_pressure = getSegmentPressure(seg); - for (const int perf : segment_perforations_[seg]) { + for (const int perf : this->segment_perforations_[seg]) { const int cell_idx = well_cells_[perf]; const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)); std::vector mob(num_components_, 0.0); @@ -4012,4 +3427,5 @@ namespace Opm const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero); connII[phase_pos] = connIICalc(mt.value() * fs.invB(flowPhaseToEbosPhaseIdx(phase_pos)).value()); } + } // namespace Opm diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 35bc1636d..0928fbe00 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -147,6 +147,14 @@ public: return perf_depth_; } + std::vector& perfDepth() { + return perf_depth_; + } + + const std::vector& wellIndex() const { + return well_index_; + } + double getTHPConstraint(const SummaryState& summaryState) const; double getALQ(const WellState& well_state) const;