From 9c431d1921b2d52b4d96c4366ba779df93ce91e2 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Tue, 20 Feb 2024 11:14:54 +0100 Subject: [PATCH] WellBhpThpCalculator: template Scalar type --- opm/simulators/wells/WellBhpThpCalculator.cpp | 433 ++++++++++-------- opm/simulators/wells/WellBhpThpCalculator.hpp | 127 +++-- 2 files changed, 296 insertions(+), 264 deletions(-) diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 9d2689c80..14a4cabaf 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -48,11 +48,11 @@ static constexpr bool extraBhpAtThpLimitOutput = false; static constexpr bool extraThpFromBhpOutput = false; -namespace Opm -{ +namespace Opm { -bool -WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) const +template +bool WellBhpThpCalculator:: +wellHasTHPConstraints(const SummaryState& summaryState) const { const auto& well_ecl = well_.wellEcl(); if (well_ecl.isInjector()) { @@ -70,7 +70,9 @@ WellBhpThpCalculator::wellHasTHPConstraints(const SummaryState& summaryState) co return false; } -double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState) const +template +Scalar WellBhpThpCalculator:: +getTHPConstraint(const SummaryState& summaryState) const { const auto& well_ecl = well_.wellEcl(); if (well_ecl.isInjector()) { @@ -86,7 +88,9 @@ double WellBhpThpCalculator::getTHPConstraint(const SummaryState& summaryState) return 0.0; } -double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const +template +Scalar WellBhpThpCalculator:: +mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const { const auto& well_ecl = well_.wellEcl(); if (well_ecl.isInjector()) { @@ -102,12 +106,14 @@ double WellBhpThpCalculator::mostStrictBhpFromBhpLimits(const SummaryState& summ return 0.0; } -double WellBhpThpCalculator::calculateThpFromBhp(const std::vector& rates, - const double bhp, - const double rho, - const std::optional& alq, - const double thp_limit, - DeferredLogger& deferred_logger) const +template +Scalar WellBhpThpCalculator:: +calculateThpFromBhp(const std::vector& rates, + const Scalar bhp, + const Scalar rho, + const std::optional& alq, + const Scalar thp_limit, + DeferredLogger& deferred_logger) const { assert(int(rates.size()) == 3); // the vfp related only supports three phases now. @@ -115,31 +121,31 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector& rate 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]; + const Scalar aqua = rates[Water]; + const Scalar liquid = rates[Oil]; + const Scalar vapour = rates[Gas]; // pick the density in the top layer - double thp = 0.0; + Scalar thp = 0.0; const int table_id = well_.wellEcl().vfp_table_number(); if (well_.isInjector()) { assert(!alq.has_value()); - const double vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth(); - const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); + const Scalar vfp_ref_depth = well_.vfpProperties()->getInj()->getTable(table_id).getDatumDepth(); + const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); auto thp_func = [this, table_id, aqua, liquid, vapour, dp]( - const double bhp_value, const double pressure_loss) { + const Scalar bhp_value, const Scalar pressure_loss) { return this->well_.vfpProperties()->getInj()->thp( table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss); }; thp = findThpFromBhpIteratively(thp_func, bhp, thp_limit, dp, deferred_logger); } else if (well_.isProducer()) { - const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth(); - const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); + const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(table_id).getDatumDepth(); + const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); auto thp_func = [this, table_id, aqua, liquid, vapour, dp, &alq] - (const double bhp_value, const double pressure_loss) { + (const Scalar bhp_value, const Scalar pressure_loss) { return this->well_.vfpProperties()->getProd()->thp( table_id, aqua, liquid, vapour, bhp_value + dp - pressure_loss, alq.value()); }; @@ -151,36 +157,35 @@ double WellBhpThpCalculator::calculateThpFromBhp(const std::vector& rate return thp; } -double -WellBhpThpCalculator:: -findThpFromBhpIteratively( - const std::function& thp_func, - const double bhp, - const double thp_limit, - const double dp, - DeferredLogger& deferred_logger) const +template +Scalar WellBhpThpCalculator:: +findThpFromBhpIteratively(const std::function& thp_func, + const Scalar bhp, + const Scalar thp_limit, + const Scalar dp, + DeferredLogger& deferred_logger) const { auto pressure_loss = getVfpBhpAdjustment(bhp + dp, thp_limit); auto thp = thp_func(bhp, pressure_loss); - const double tolerance = 1e-5 * unit::barsa; + const Scalar tolerance = 1e-5 * unit::barsa; bool do_iterate = true; int it = 1; int max_iterations = 50; - while(do_iterate) { + while (do_iterate) { if (it > max_iterations) { break; } - double thp_prev = thp; + Scalar thp_prev = thp; pressure_loss = getVfpBhpAdjustment(bhp + dp - pressure_loss, thp_prev); thp = thp_func(bhp, pressure_loss); - auto error = std::fabs(thp-thp_prev); + auto error = std::fabs(thp - thp_prev); if (extraThpFromBhpOutput) { const std::string msg = fmt::format( "findThpFromBhpIteratively(): iteration {}, thp = {}, bhp = {}, " "pressure_loss = {}, error = {}", it, thp, bhp+dp-pressure_loss, pressure_loss, error); deferred_logger.debug(msg); } - if (std::fabs(thp-thp_prev) < tolerance) { + if (std::fabs(thp - thp_prev) < tolerance) { break; } it++; @@ -188,14 +193,15 @@ findThpFromBhpIteratively( return thp; } -std::optional -WellBhpThpCalculator:: -computeBhpAtThpLimitProd(const std::function(const double)>& frates, +template +std::optional +WellBhpThpCalculator:: +computeBhpAtThpLimitProd(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double maxPerfPress, - const double rho, - const double alq_value, - const double thp_limit, + const Scalar maxPerfPress, + const Scalar rho, + const Scalar alq_value, + const Scalar thp_limit, DeferredLogger& deferred_logger) const { // Given a VFP function returning bhp as a function of phase @@ -221,34 +227,43 @@ computeBhpAtThpLimitProd(const std::function(const double)>& // Make the fbhp() function. const auto& controls = well_.wellEcl().productionControls(summary_state); const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); - const double vfp_ref_depth = table.getDatumDepth(); - const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); + const Scalar vfp_ref_depth = table.getDatumDepth(); + const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), + vfp_ref_depth, + rho, + well_.gravity()); - auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp, alq_value](const std::vector& rates) { assert(rates.size() == 3); - const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); - const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); + const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, + well_.indexOfWell()); + const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, + well_.indexOfWell()); const bool use_vfpexp = well_.useVfpExplicit(); - const double bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, - rates[Water], - rates[Oil], - rates[Gas], - thp_limit, - alq_value, - wfr, - gfr, - use_vfpexp); + const Scalar bhp = well_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, + rates[Water], + rates[Oil], + rates[Gas], + thp_limit, + alq_value, + wfr, + gfr, + use_vfpexp); return bhp - dp + getVfpBhpAdjustment(bhp, thp_limit); }; // Make the flo() function. - auto flo = [&table](const std::vector& rates) { + 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. - auto fflo = [&flo, &frates](double bhp) { return flo(frates(bhp)); }; - auto bhp_max = this->bhpMax(fflo, controls.bhp_limit, maxPerfPress, table.getFloAxis().front(), deferred_logger); + auto fflo = [&flo, &frates](Scalar bhp) { return flo(frates(bhp)); }; + auto bhp_max = this->bhpMax(fflo, + controls.bhp_limit, + maxPerfPress, + table.getFloAxis().front(), + deferred_logger); // could not solve for the bhp-point, we could not continue to find the bhp if (!bhp_max.has_value()) { @@ -257,16 +272,17 @@ computeBhpAtThpLimitProd(const std::function(const double)>& "find bhp-point where production becomes non-zero for well " + well_.name()); return std::nullopt; } - const std::array range {controls.bhp_limit, *bhp_max}; + const std::array range {controls.bhp_limit, *bhp_max}; return this->computeBhpAtThpLimit(frates, fbhp, range, deferred_logger); } -std::optional -WellBhpThpCalculator:: -computeBhpAtThpLimitInj(const std::function(const double)>& frates, +template +std::optional +WellBhpThpCalculator:: +computeBhpAtThpLimitInj(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double rho, - const double flo_rel_tol, + const Scalar rho, + const Scalar flo_rel_tol, const int max_iteration, const bool throwOnError, DeferredLogger& deferred_logger) const @@ -282,13 +298,15 @@ computeBhpAtThpLimitInj(const std::function(const double)>& } } -void WellBhpThpCalculator::updateThp(const double rho, - const bool stop_or_zero_rate_target, - const std::function& alq_value, - const std::array& active, - WellState& well_state, - const SummaryState& summary_state, - DeferredLogger& deferred_logger) const +template +void WellBhpThpCalculator:: +updateThp(const Scalar rho, + const bool stop_or_zero_rate_target, + const std::function& alq_value, + const std::array& active, + WellState& well_state, + const SummaryState& summary_state, + DeferredLogger& deferred_logger) const { static constexpr int Gas = BlackoilPhases::Vapour; static constexpr int Oil = BlackoilPhases::Liquid; @@ -308,7 +326,7 @@ void WellBhpThpCalculator::updateThp(const double rho, } // the well is under other control types, we calculate the thp based on bhp and rates - std::vector rates(3, 0.0); + std::vector rates(3, 0.0); const PhaseUsage& pu = well_.phaseUsage(); if (active[Water]) { @@ -320,18 +338,19 @@ void WellBhpThpCalculator::updateThp(const double rho, if (active[Gas]) { rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ]; } - const std::optional alq = this->well_.isProducer() ? std::optional(alq_value()) : std::nullopt; - const double thp_limit = well_.getTHPConstraint(summary_state); + const std::optional alq = this->well_.isProducer() ? std::optional(alq_value()) : std::nullopt; + const Scalar thp_limit = well_.getTHPConstraint(summary_state); ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, alq, thp_limit, deferred_logger); } +template template -EvalWell WellBhpThpCalculator:: -calculateBhpFromThp(const WellState& well_state, +EvalWell WellBhpThpCalculator:: +calculateBhpFromThp(const WellState& well_state, const std::vector& rates, const Well& well, const SummaryState& summaryState, - const double rho, + const Scalar rho, DeferredLogger& deferred_logger) const { // TODO: when well is under THP control, the BHP is dependent on the rates, @@ -349,8 +368,8 @@ calculateBhpFromThp(const WellState& well_state, const EvalWell aqua = rates[Water]; const EvalWell liquid = rates[Oil]; const EvalWell vapour = rates[Gas]; - const double thp_limit = well_.getTHPConstraint(summaryState); - double vfp_ref_depth; + const Scalar thp_limit = well_.getTHPConstraint(summaryState); + Scalar vfp_ref_depth; EvalWell bhp_tab; if (well_.isInjector() ) { @@ -375,57 +394,62 @@ calculateBhpFromThp(const WellState& well_state, else { OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + well_.name(), deferred_logger); } - double bhp_tab_double_value; - if constexpr (std::is_same_v) { + Scalar bhp_tab_double_value; + if constexpr (std::is_same_v) { bhp_tab_double_value = bhp_tab; } else { // EvalWell and bhp_tab is of type DenseAd::Evaluation bhp_tab_double_value = bhp_tab.value(); } const auto bhp_adjustment = getVfpBhpAdjustment(bhp_tab_double_value, thp_limit); - const double dp_hydro = wellhelpers::computeHydrostaticCorrection( - well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); + const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), + vfp_ref_depth, + rho, + well_.gravity()); return bhp_tab - dp_hydro + bhp_adjustment; } -double -WellBhpThpCalculator:: -calculateMinimumBhpFromThp(const WellState& well_state, +template +Scalar WellBhpThpCalculator:: +calculateMinimumBhpFromThp(const WellState& well_state, const Well& well, const SummaryState& summaryState, - const double rho) const + const Scalar rho) const { assert(well_.isProducer()); // only producers can go here for now - const double thp_limit = well_.getTHPConstraint(summaryState); + const Scalar thp_limit = well_.getTHPConstraint(summaryState); const auto& controls = well.productionControls(summaryState); const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); const auto& gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); - const double bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number, + const Scalar bhp_min = well_.vfpProperties()->getProd()->minimumBHP(controls.vfp_table_number, thp_limit, wfr, gfr, well_.getALQ(well_state)); - const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); + const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); const auto bhp_adjustment = getVfpBhpAdjustment(bhp_min, thp_limit); - const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, + const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); return bhp_min - dp_hydro + bhp_adjustment; } -double WellBhpThpCalculator::getVfpBhpAdjustment(const double bhp_tab, const double thp_limit) const +template +Scalar WellBhpThpCalculator:: +getVfpBhpAdjustment(const Scalar bhp_tab, const Scalar thp_limit) const { return well_.wellEcl().getWVFPDP().getPressureLoss(bhp_tab, thp_limit); } +template template -std::optional -WellBhpThpCalculator:: -computeBhpAtThpLimitInjImpl(const std::function(const double)>& frates, +std::optional +WellBhpThpCalculator:: +computeBhpAtThpLimitInjImpl(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double rho, - const double flo_rel_tol, + const Scalar rho, + const Scalar flo_rel_tol, const int max_iteration, DeferredLogger& deferred_logger) const { @@ -475,12 +499,12 @@ computeBhpAtThpLimitInjImpl(const std::function(const double // Make the fbhp() function. const auto& controls = well_.wellEcl().injectionControls(summary_state); const auto& table = well_.vfpProperties()->getInj()->getTable(controls.vfp_table_number); - const double vfp_ref_depth = table.getDatumDepth(); - const double thp_limit = this->getTHPConstraint(summary_state); - const double dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), + const Scalar vfp_ref_depth = table.getDatumDepth(); + const Scalar thp_limit = this->getTHPConstraint(summary_state); + const Scalar dp = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); - auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { + auto fbhp = [this, &controls, thp_limit, dp](const std::vector& rates) { assert(rates.size() == 3); const auto bhp = well_.vfpProperties()->getInj() ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit); @@ -488,25 +512,25 @@ computeBhpAtThpLimitInjImpl(const std::function(const double }; // Make the flo() function. - auto flo = [&table](const std::vector& rates) { + 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(); + std::vector flo_samples = table.getFloAxis(); if (flo_samples[0] > 0.0) { - const double f0 = flo_samples[0]; + const Scalar 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)); + const Scalar 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) { + std::vector bhp_samples; + for (Scalar 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 @@ -516,16 +540,16 @@ computeBhpAtThpLimitInjImpl(const std::function(const double bhp_samples.push_back(controls.bhp_limit); break; } - auto eq = [&flo, &frates, flo_sample](double bhp) { + auto eq = [&flo, &frates, flo_sample](Scalar 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 double flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back()); + const Scalar low = 10.0 * unit::barsa; + const Scalar high = 800.0 * unit::barsa; + const Scalar flo_tolerance = flo_rel_tol * std::fabs(flo_samples.back()); int iteration = 0; try { - const double solved_bhp = RegulaFalsiBisection:: + const Scalar solved_bhp = RegulaFalsiBisection:: solve(eq, low, high, max_iteration, flo_tolerance, iteration); bhp_samples.push_back(solved_bhp); } @@ -539,7 +563,7 @@ computeBhpAtThpLimitInjImpl(const std::function(const double // 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); + std::vector fbhp_samples(num_samples); for (int ii = 0; ii < num_samples; ++ii) { fbhp_samples[ii] = fbhp(frates(bhp_samples[ii])); } @@ -564,8 +588,8 @@ computeBhpAtThpLimitInjImpl(const std::function(const double // 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]; + const Scalar curr = fbhp_samples[ii] - bhp_samples[ii]; + const Scalar 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. @@ -578,13 +602,13 @@ computeBhpAtThpLimitInjImpl(const std::function(const double } // Solve for the proper solution in the given interval. - auto eq = [&fbhp, &frates](double bhp) { + auto eq = [&fbhp, &frates](Scalar 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 double bhp_tolerance = 0.01 * unit::barsa; + const Scalar low = bhp_samples[sign_change_index + 1]; + const Scalar high = bhp_samples[sign_change_index]; + const Scalar bhp_tolerance = 0.01 * unit::barsa; int iteration = 0; if (low == high) { // We are in the high flow regime where the bhp_samples @@ -595,7 +619,7 @@ computeBhpAtThpLimitInjImpl(const std::function(const double return std::nullopt; } try { - const double solved_bhp = RegulaFalsiBisection:: + const Scalar solved_bhp = RegulaFalsiBisection:: solve(eq, low, high, max_iteration, bhp_tolerance, iteration); if constexpr (extraBhpAtThpLimitOutput) { OpmLog::debug("***** " + well_.name() + " solved_bhp = " + std::to_string(solved_bhp) @@ -610,20 +634,21 @@ computeBhpAtThpLimitInjImpl(const std::function(const double } } -std::optional -WellBhpThpCalculator:: -bhpMax(const std::function& fflo, - const double bhp_limit, - const double maxPerfPress, - const double vfp_flo_front, +template +std::optional +WellBhpThpCalculator:: +bhpMax(const std::function& fflo, + const Scalar bhp_limit, + const Scalar maxPerfPress, + const Scalar vfp_flo_front, DeferredLogger& deferred_logger) const { // Find the bhp-point where production becomes nonzero. - double bhp_max = 0.0; - double low = bhp_limit; - double high = maxPerfPress + 1.0 * unit::barsa; - double f_low = fflo(low); - double f_high = fflo(high); + Scalar bhp_max = 0.0; + Scalar low = bhp_limit; + Scalar high = maxPerfPress + 1.0 * unit::barsa; + Scalar f_low = fflo(low); + Scalar f_high = fflo(high); if constexpr (extraBhpAtThpLimitOutput) { deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() + " low = " + std::to_string(low) + @@ -633,7 +658,7 @@ bhpMax(const std::function& fflo, } int adjustments = 0; const int max_adjustments = 10; - const double adjust_amount = 5.0 * unit::barsa; + const Scalar 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; @@ -656,12 +681,12 @@ bhpMax(const std::function& fflo, } else { // Bisect to find a bhp point where we produce, but // not a large amount ('eps' below). - const double eps = 0.1 * std::fabs(vfp_flo_front); + const Scalar eps = 0.1 * std::fabs(vfp_flo_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); + const Scalar curr = 0.5*(low + high); + const Scalar f_curr = fflo(curr); if (f_curr * f_low > 0.0) { low = curr; f_low = f_curr; @@ -690,11 +715,12 @@ bhpMax(const std::function& fflo, return bhp_max; } -std::optional -WellBhpThpCalculator:: -computeBhpAtThpLimit(const std::function(const double)>& frates, - const std::function)>& fbhp, - const std::array& range, +template +std::optional +WellBhpThpCalculator:: +computeBhpAtThpLimit(const std::function(const Scalar)>& frates, + const std::function)>& fbhp, + const std::array& range, DeferredLogger& deferred_logger) const { // Given a VFP function returning bhp as a function of phase @@ -714,15 +740,16 @@ computeBhpAtThpLimit(const std::function(const double)>& fra // highest rate) should be returned. // Define the equation we want to solve. - auto eq = [&fbhp, &frates](double bhp) { + auto eq = [&fbhp, &frates](Scalar bhp) { return fbhp(frates(bhp)) - bhp; }; // Find appropriate brackets for the solution. - std::optional approximate_solution; - double low, high; + std::optional approximate_solution; + Scalar low, high; // trying to use bisect way to locate a bracket - bool finding_bracket = this->bisectBracket(eq, range, low, high, approximate_solution, deferred_logger); + bool finding_bracket = this->bisectBracket(eq, range, low, high, + approximate_solution, deferred_logger); // based on the origional design, if an approximate solution is suggested, we use this value directly // in the long run, we might change it @@ -744,10 +771,10 @@ computeBhpAtThpLimit(const std::function(const double)>& fra // Solve for the proper solution in the given interval. const int max_iteration = 100; - const double bhp_tolerance = 0.01 * unit::barsa; + const Scalar bhp_tolerance = 0.01 * unit::barsa; int iteration = 0; try { - const double solved_bhp = RegulaFalsiBisection:: + const Scalar solved_bhp = RegulaFalsiBisection:: solve(eq, low, high, max_iteration, bhp_tolerance, iteration); return solved_bhp; } @@ -758,21 +785,21 @@ computeBhpAtThpLimit(const std::function(const double)>& fra } } -bool -WellBhpThpCalculator:: -bisectBracket(const std::function& eq, - const std::array& range, - double& low, double& high, - std::optional& approximate_solution, +template +bool WellBhpThpCalculator:: +bisectBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, + std::optional& approximate_solution, DeferredLogger& deferred_logger) const { bool finding_bracket = false; low = range[0]; high = range[1]; - double eq_high = eq(high); - double eq_low = eq(low); - const double eq_bhplimit = eq_low; + Scalar eq_high = eq(high); + Scalar eq_low = eq(low); + const Scalar eq_bhplimit = eq_low; if constexpr (extraBhpAtThpLimitOutput) { deferred_logger.debug("computeBhpAtThpLimitProd(): well = " + well_.name() + " low = " + std::to_string(low) + @@ -783,12 +810,12 @@ bisectBracket(const std::function& eq, 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); + Scalar abs_low = std::fabs(eq_low); + Scalar 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; + Scalar interval = high - low; + const Scalar 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); @@ -815,7 +842,7 @@ bisectBracket(const std::function& eq, } } else { // eq_low * eq_high > 0.0 // Still failed bracketing! - const double limit = 0.1 * unit::barsa; + const Scalar limit = 0.1 * unit::barsa; if (std::min(abs_low, abs_high) < limit) { // Return the least bad solution if less off than 0.1 bar. deferred_logger.warning("FAILED_ROBUST_BHP_THP_SOLVE_BRACKETING_FAILURE", @@ -834,20 +861,20 @@ bisectBracket(const std::function& eq, return finding_bracket; } -bool -WellBhpThpCalculator:: -bruteForceBracket(const std::function& eq, - const std::array& range, - double& low, double& high, +template +bool WellBhpThpCalculator:: +bruteForceBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, DeferredLogger& deferred_logger) { bool bracket_found = false; low = range[0]; high = range[1]; const int sample_number = 200; - const double interval = (high - low) / sample_number; - double eq_low = eq(low); - double eq_high = 0.0; + const Scalar interval = (high - low) / sample_number; + Scalar eq_low = eq(low); + Scalar eq_high = 0.0; for (int i = 0; i < sample_number + 1; ++i) { high = range[0] + interval * i; eq_high = eq(high); @@ -866,10 +893,11 @@ bruteForceBracket(const std::function& eq, return bracket_found; } -bool WellBhpThpCalculator:: -isStableSolution(const WellState& well_state, +template +bool WellBhpThpCalculator:: +isStableSolution(const WellState& well_state, const Well& well, - const std::vector& rates, + const std::vector& rates, const SummaryState& summaryState) const { assert(int(rates.size()) == 3); // the vfp related only supports three phases now. @@ -879,10 +907,10 @@ isStableSolution(const WellState& well_state, static constexpr int Oil = BlackoilPhases::Liquid; static constexpr int Water = BlackoilPhases::Aqua; - const double aqua = rates[Water]; - const double liquid = rates[Oil]; - const double vapour = rates[Gas]; - const double thp = well_.getTHPConstraint(summaryState); + const Scalar aqua = rates[Water]; + const Scalar liquid = rates[Oil]; + const Scalar vapour = rates[Gas]; + const Scalar thp = well_.getTHPConstraint(summaryState); const auto& controls = well.productionControls(summaryState); const auto& wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); @@ -897,27 +925,28 @@ isStableSolution(const WellState& well_state, return true; } else { // maybe check if ipr is available const auto ipr = getFloIPR(well_state, well, summaryState); - return bhp.dflo + 1/ipr.second >= 0; + return bhp.dflo + 1.0 / ipr.second >= 0; } } -std::optional WellBhpThpCalculator:: -estimateStableBhp(const WellState& well_state, +template +std::optional WellBhpThpCalculator:: +estimateStableBhp(const WellState& well_state, const Well& well, - const std::vector& rates, - const double rho, + const std::vector& rates, + const Scalar rho, const SummaryState& summaryState) const { // Given a *converged* well_state with ipr, estimate bhp of the stable solution const auto& controls = well.productionControls(summaryState); - const double thp = well_.getTHPConstraint(summaryState); + const Scalar thp = well_.getTHPConstraint(summaryState); const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); - const double aqua = rates[BlackoilPhases::Aqua]; - const double liquid = rates[BlackoilPhases::Liquid]; - const double vapour = rates[BlackoilPhases::Vapour]; - double flo = detail::getFlo(table, aqua, liquid, vapour); - double wfr, gfr; + const Scalar aqua = rates[BlackoilPhases::Aqua]; + const Scalar liquid = rates[BlackoilPhases::Liquid]; + const Scalar vapour = rates[BlackoilPhases::Vapour]; + Scalar flo = detail::getFlo(table, aqua, liquid, vapour); + Scalar wfr, gfr; if (well_.useVfpExplicit() || -flo < table.getFloAxis().front()) { wfr = well_.vfpProperties()->getExplicitWFR(controls.vfp_table_number, well_.indexOfWell()); gfr = well_.vfpProperties()->getExplicitGFR(controls.vfp_table_number, well_.indexOfWell()); @@ -928,11 +957,11 @@ estimateStableBhp(const WellState& well_state, auto ipr = getFloIPR(well_state, well, summaryState); - const double vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); + const Scalar vfp_ref_depth = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth(); - const double dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, + const Scalar dp_hydro = wellhelpers::computeHydrostaticCorrection(well_.refDepth(), vfp_ref_depth, rho, well_.gravity()); - auto bhp_adjusted = [this, &thp, &dp_hydro](const double bhp) { + auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) { return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp); }; const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted); @@ -942,10 +971,11 @@ estimateStableBhp(const WellState& well_state, } else { return std::nullopt; } -} +} -std::pair WellBhpThpCalculator:: -getFloIPR(const WellState& well_state, +template +std::pair WellBhpThpCalculator:: +getFloIPR(const WellState& well_state, const Well& well, const SummaryState& summary_state) const { @@ -954,21 +984,23 @@ getFloIPR(const WellState& well_state, const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const auto& pu = well_.phaseUsage(); const auto& ipr_a = well_state.well(well_.indexOfWell()).implicit_ipr_a; - const double& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; - const double& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; - const double& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; + const Scalar& aqua_a = pu.phase_used[BlackoilPhases::Aqua]? ipr_a[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; + const Scalar& liquid_a = pu.phase_used[BlackoilPhases::Liquid]? ipr_a[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; + const Scalar& vapour_a = pu.phase_used[BlackoilPhases::Vapour]? ipr_a[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; const auto& ipr_b = well_state.well(well_.indexOfWell()).implicit_ipr_b; - const double& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; - const double& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; - const double& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; + const Scalar& aqua_b = pu.phase_used[BlackoilPhases::Aqua]? ipr_b[pu.phase_pos[BlackoilPhases::Aqua]] : 0.0; + const Scalar& liquid_b = pu.phase_used[BlackoilPhases::Liquid]? ipr_b[pu.phase_pos[BlackoilPhases::Liquid]] : 0.0; + const Scalar& vapour_b = pu.phase_used[BlackoilPhases::Vapour]? ipr_b[pu.phase_pos[BlackoilPhases::Vapour]] : 0.0; // The getFlo helper is indended to pick one or add two of the phase rates (depending on FLO-type), // but we can equally use it to pick/add the corresponding ipr_a, ipr_b return std::make_pair(detail::getFlo(table, aqua_a, liquid_a, vapour_a), detail::getFlo(table, aqua_b, liquid_b, vapour_b)); } +template class WellBhpThpCalculator; + #define INSTANCE(...) \ -template __VA_ARGS__ WellBhpThpCalculator:: \ +template __VA_ARGS__ WellBhpThpCalculator:: \ calculateBhpFromThp<__VA_ARGS__>(const WellState&, \ const std::vector<__VA_ARGS__>&, \ const Well&, \ @@ -993,4 +1025,5 @@ INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) + } // namespace Opm diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index fdf0502ca..10c059e01 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -26,11 +26,9 @@ #include #include -#include #include -namespace Opm -{ +namespace Opm { class DeferredLogger; class SummaryState; @@ -39,136 +37,137 @@ template class WellInterfaceGeneric; template class WellState; //! \brief Class for computing BHP limits. +template class WellBhpThpCalculator { public: //! \brief Constructor sets reference to well. - WellBhpThpCalculator(const WellInterfaceGeneric& well) : well_(well) {} + WellBhpThpCalculator(const WellInterfaceGeneric& well) : well_(well) {} //! \brief Checks if well has THP constraints. bool wellHasTHPConstraints(const SummaryState& summaryState) const; //! \brief Get THP constraint for well. - double getTHPConstraint(const SummaryState& summaryState) const; + Scalar getTHPConstraint(const SummaryState& summaryState) const; //! \brief Obtain the most strict BHP from BHP limits. - double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; + Scalar mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const; //! \brief Calculates THP from BHP. - double calculateThpFromBhp(const std::vector& rates, - const double bhp, - const double rho, - const std::optional& alq, - const double thp_limit, + Scalar calculateThpFromBhp(const std::vector& rates, + const Scalar bhp, + const Scalar rho, + const std::optional& alq, + const Scalar thp_limit, DeferredLogger& deferred_logger) const; //! \brief Compute BHP from THP limit for a producer. - std::optional - computeBhpAtThpLimitProd(const std::function(const double)>& frates, + std::optional + computeBhpAtThpLimitProd(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double maxPerfPress, - const double rho, - const double alq_value, - const double thp_limit, + const Scalar maxPerfPress, + const Scalar rho, + const Scalar alq_value, + const Scalar thp_limit, DeferredLogger& deferred_logger) const; //! \brief Compute BHP from THP limit for an injector. - std::optional - computeBhpAtThpLimitInj(const std::function(const double)>& frates, + std::optional + computeBhpAtThpLimitInj(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double rho, - const double flo_rel_tol, + const Scalar rho, + const Scalar flo_rel_tol, const int max_iteration, const bool throwOnError, DeferredLogger& deferred_logger) const; //! \brief Update THP. - void updateThp(const double rho, + void updateThp(const Scalar rho, const bool stop_or_zero_rate_target, - const std::function& alq_value, + const std::function& alq_value, const std::array& active, - WellState& well_state, + WellState& well_state, const SummaryState& summary_state, DeferredLogger& deferred_logger) const; template - EvalWell calculateBhpFromThp(const WellState& well_state, + EvalWell calculateBhpFromThp(const WellState& well_state, const std::vector& rates, const Well& well, const SummaryState& summaryState, - const double rho, + const Scalar rho, DeferredLogger& deferred_logger) const; - double calculateMinimumBhpFromThp(const WellState& well_state, - const Well& well, - const SummaryState& summaryState, - const double rho) const; + Scalar calculateMinimumBhpFromThp(const WellState& well_state, + const Well& well, + const SummaryState& summaryState, + const Scalar rho) const; - bool isStableSolution(const WellState& well_state, + bool isStableSolution(const WellState& well_state, const Well& well, - const std::vector& rates, + const std::vector& rates, const SummaryState& summaryState) const; - std::optional - estimateStableBhp (const WellState& well_state, + std::optional + estimateStableBhp (const WellState& well_state, const Well& well, - const std::vector& rates, - const double rho, + const std::vector& rates, + const Scalar rho, const SummaryState& summaryState) const; - std::pair - getFloIPR(const WellState& well_state, + std::pair + getFloIPR(const WellState& well_state, const Well& well, const SummaryState& summary_state) const; private: //! \brief Compute BHP from THP limit for an injector - implementation. template - std::optional - computeBhpAtThpLimitInjImpl(const std::function(const double)>& frates, + std::optional + computeBhpAtThpLimitInjImpl(const std::function(const Scalar)>& frates, const SummaryState& summary_state, - const double rho, - const double flo_rel_tol, + const Scalar rho, + const Scalar flo_rel_tol, const int max_iteration, DeferredLogger& deferred_logger) const; //! \brief Calculate max BHP. - std::optional - bhpMax(const std::function& fflo, - const double bhp_limit, - const double maxPerfPress, - const double vfp_flo_front, + std::optional + bhpMax(const std::function& fflo, + const Scalar bhp_limit, + const Scalar maxPerfPress, + const Scalar vfp_flo_front, DeferredLogger& deferred_logger) const; //! \brief Common code for finding BHP from THP limit for producers/injectors. - std::optional - computeBhpAtThpLimit(const std::function(const double)>& frates, - const std::function)>& fbhp, - const std::array& range, + std::optional + computeBhpAtThpLimit(const std::function(const Scalar)>& frates, + const std::function)>& fbhp, + const std::array& range, DeferredLogger& deferred_logger) const; //! \brief Get pressure adjustment to the bhp calculated from VFP table - double getVfpBhpAdjustment(const double bph_tab, const double thp_limit) const; + Scalar getVfpBhpAdjustment(const Scalar bph_tab, const Scalar thp_limit) const; //! \brief Find limits using bisection. - bool bisectBracket(const std::function& eq, - const std::array& range, - double& low, double& high, - std::optional& approximate_solution, + bool bisectBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, + std::optional& approximate_solution, DeferredLogger& deferred_logger) const; //! \brief Find limits using brute-force solver. - static bool bruteForceBracket(const std::function& eq, - const std::array& range, - double& low, double& high, + static bool bruteForceBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, DeferredLogger& deferred_logger); - double findThpFromBhpIteratively(const std::function& thp_func, - const double bhp, - const double thp_limit, - const double dp, + Scalar findThpFromBhpIteratively(const std::function& thp_func, + const Scalar bhp, + const Scalar thp_limit, + const Scalar dp, DeferredLogger& deferred_logger) const; - const WellInterfaceGeneric& well_; //!< Reference to well interface + const WellInterfaceGeneric& well_; //!< Reference to well interface }; }