OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #5375 from akva2/msw_use_scalar

Browse Source 
MultisegmentWell: use Scalar type
This commit is contained in:
Bård Skaflestad 2024-05-22 15:18:45 +02:00 committed by GitHub
commit 3544dfcabd
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
12 changed files with 309 additions and 307 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -25,8 +25,8 @@
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/MultisegmentWellEval.hpp>
namespace Opm
{
namespace Opm {
class DeferredLogger;
template<typename TypeTag>
@ -79,8 +79,8 @@ namespace Opm
const std::vector<PerforationData<Scalar>>& perf_data);
void init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const std::vector<Scalar>& depth_arg,
const Scalar gravity_arg,
const int num_cells,
const std::vector<Scalar>& B_avg,
const bool changed_to_open_this_step) override;
@ -96,7 +96,7 @@ namespace Opm
/// check whether the well equations get converged for this well
ConvergenceReport getWellConvergence(const SummaryState& summary_state,
const WellState<Scalar>& well_state,
const std::vector<double>& B_avg,
const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger,
const bool relax_tolerance) const override;
@ -115,7 +115,7 @@ namespace Opm
/// computing the well potentials for group control
void computeWellPotentials(const Simulator& simulator,
const WellState<Scalar>& well_state,
std::vector<double>& well_potentials,
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger) override;
void updatePrimaryVariables(const SummaryState& summary_state,
@ -139,7 +139,7 @@ namespace Opm
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger) const override;
double connectionDensity(const int globalConnIdx,
Scalar connectionDensity(const int globalConnIdx,
const int openConnIdx) const override;
void addWellContributions(SparseMatrixAdapter& jacobian) const override;
@ -150,26 +150,26 @@ namespace Opm
const bool use_well_weights,
const WellState<Scalar>& well_state) const override;
std::vector<double> computeCurrentWellRates(const Simulator& simulator,
std::vector<Scalar>
computeCurrentWellRates(const Simulator& simulator,
DeferredLogger& deferred_logger) const override;
std::optional<double>
std::optional<Scalar>
computeBhpAtThpLimitProdWithAlq(const Simulator& simulator,
const SummaryState& summary_state,
const double alq_value,
const Scalar alq_value,
DeferredLogger& deferred_logger) const override;
std::vector<double> getPrimaryVars() const override;
std::vector<Scalar> getPrimaryVars() const override;
int setPrimaryVars(std::vector<double>::const_iterator it) override;
int setPrimaryVars(typename std::vector<Scalar>::const_iterator it) override;
protected:
// regularize msw equation
bool regularize_;
// the intial amount of fluids in each segment under surface condition
std::vector<std::vector<double> > segment_fluid_initial_;
std::vector<std::vector<Scalar> > segment_fluid_initial_;
mutable int debug_cost_counter_ = 0;
@ -178,8 +178,7 @@ namespace Opm
const BVectorWell& dwells,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger,
const double relaxation_factor = 1.0);
const Scalar relaxation_factor = 1.0);
// computing the accumulation term for later use in well mass equations
void computeInitialSegmentFluids(const Simulator& simulator);
@ -230,31 +229,31 @@ namespace Opm
DeferredLogger& deferred_logger) const;
void computeWellRatesAtBhpLimit(const Simulator& simulator,
std::vector<double>& well_flux,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const;
void computeWellRatesWithBhp(const Simulator& simulator,
const double& bhp,
std::vector<double>& well_flux,
const Scalar& bhp,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const override;
void computeWellRatesWithBhpIterations(const Simulator& simulator,
const Scalar& bhp,
std::vector<double>& well_flux,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const override;
std::vector<double>
std::vector<Scalar>
computeWellPotentialWithTHP(const WellState<Scalar>& well_state,
const Simulator& simulator,
DeferredLogger& deferred_logger) const;
bool computeWellPotentialsImplicit(const Simulator& simulator,
std::vector<double>& well_potentials,
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const override;
Scalar getRefDensity() const override;
virtual bool iterateWellEqWithControl(const Simulator& simulator,
bool iterateWellEqWithControl(const Simulator& simulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
@ -262,7 +261,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger) override;
virtual bool iterateWellEqWithSwitching(const Simulator& simulator,
bool iterateWellEqWithSwitching(const Simulator& simulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
@ -272,7 +271,7 @@ namespace Opm
const bool fixed_control = false,
const bool fixed_status = false) override;
virtual void assembleWellEqWithoutIteration(const Simulator& simulator,
void assembleWellEqWithoutIteration(const Simulator& simulator,
const double dt,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
@ -280,7 +279,7 @@ namespace Opm
const GroupState<Scalar>& group_state,
DeferredLogger& deferred_logger) override;
virtual void updateWaterThroughput(const double dt, WellState<Scalar>& well_state) const override;
void updateWaterThroughput(const double dt, WellState<Scalar>& well_state) const override;
EvalWell getSegmentSurfaceVolume(const Simulator& simulator, const int seg_idx) const;
@ -295,20 +294,18 @@ namespace Opm
// be able to produce/inject .
bool allDrawDownWrongDirection(const Simulator& simulator) const;
std::optional<double>
std::optional<Scalar>
computeBhpAtThpLimitProd(const WellState<Scalar>& well_state,
const Simulator& simulator,
const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const;
std::optional<double>
computeBhpAtThpLimitInj(const Simulator& simulator,
std::optional<Scalar>
computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const;
double maxPerfPress(const Simulator& simulator) const;
Scalar maxPerfPress(const Simulator& simulator) const;
// check whether the well is operable under BHP limit with current reservoir condition
void checkOperabilityUnderBHPLimit(const WellState<Scalar>& well_state,

4
opm/simulators/wells/MultisegmentWellAssemble.cpp
View File

@ -89,7 +89,7 @@ assembleControlEq(const WellState<Scalar>& well_state,
const SummaryState& summaryState,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double rho,
const Scalar rho,
const PrimaryVariables& primary_variables,
Equations& eqns1,
DeferredLogger& deferred_logger) const
@ -121,7 +121,7 @@ assembleControlEq(const WellState<Scalar>& well_state,
} else if (well_.isInjector() ) {
// Find scaling factor to get injection rate,
const InjectorType injectorType = inj_controls.injector_type;
double scaling;
Scalar scaling;
const auto& pu = well_.phaseUsage();
switch (injectorType) {
case InjectorType::WATER:

2
opm/simulators/wells/MultisegmentWellAssemble.hpp
View File

@ -76,7 +76,7 @@ public:
const SummaryState& summaryState,
const Well::InjectionControls& inj_controls,
const Well::ProductionControls& prod_controls,
const double rho,
const Scalar rho,
const PrimaryVariables& primary_variables,
Equations& eqns,
DeferredLogger& deferred_logger) const;

61
opm/simulators/wells/MultisegmentWellEval.cpp
View File

@ -79,35 +79,35 @@ template<typename FluidSystem, typename Indices>
ConvergenceReport
MultisegmentWellEval<FluidSystem,Indices>::
getWellConvergence(const WellState<Scalar>& well_state,
const std::vector<double>& B_avg,
const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger,
const double max_residual_allowed,
const double tolerance_wells,
const double relaxed_inner_tolerance_flow_ms_well,
const double tolerance_pressure_ms_wells,
const double relaxed_inner_tolerance_pressure_ms_well,
const Scalar max_residual_allowed,
const Scalar tolerance_wells,
const Scalar relaxed_inner_tolerance_flow_ms_well,
const Scalar tolerance_pressure_ms_wells,
const Scalar relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance,
const bool well_is_stopped) const
{
assert(int(B_avg.size()) == baseif_.numComponents());
// checking if any residual is NaN or too large. The two large one is only handled for the well flux
std::vector<std::vector<double>> abs_residual(this->numberOfSegments(),
std::vector<double>(numWellEq, 0.0));
std::vector<std::vector<Scalar>> abs_residual(this->numberOfSegments(),
std::vector<Scalar>(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(linSys_.residual()[seg][eq_idx]);
}
}
std::vector<double> maximum_residual(numWellEq, 0.0);
std::vector<Scalar> maximum_residual(numWellEq, 0.0);
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 < this->numberOfSegments(); ++seg) {
if (eq_idx < baseif_.numComponents()) { // phase or component mass equations
const double flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
const Scalar flux_residual = B_avg[eq_idx] * abs_residual[seg][eq_idx];
if (flux_residual > maximum_residual[eq_idx]) {
maximum_residual[eq_idx] = flux_residual;
}
@ -115,7 +115,7 @@ getWellConvergence(const WellState<Scalar>& well_state,
// for the top segment (seg == 0), it is control equation, will be checked later separately
if (seg > 0) {
// Pressure equation
const double pressure_residual = abs_residual[seg][eq_idx];
const Scalar pressure_residual = abs_residual[seg][eq_idx];
if (pressure_residual > maximum_residual[eq_idx]) {
maximum_residual[eq_idx] = pressure_residual;
}
@ -127,7 +127,7 @@ getWellConvergence(const WellState<Scalar>& well_state,
using CR = ConvergenceReport;
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
if (eq_idx < baseif_.numComponents()) { // phase or component mass equations
const double flux_residual = maximum_residual[eq_idx];
const Scalar flux_residual = maximum_residual[eq_idx];
// TODO: the report can not handle the segment number yet.
if (std::isnan(flux_residual)) {
@ -140,7 +140,7 @@ getWellConvergence(const WellState<Scalar>& well_state,
report.setWellFailed({CR::WellFailure::Type::MassBalance, CR::Severity::Normal, eq_idx, baseif_.name()});
}
} else { // pressure equation
const double pressure_residual = maximum_residual[eq_idx];
const Scalar pressure_residual = maximum_residual[eq_idx];
const int dummy_component = -1;
if (std::isnan(pressure_residual)) {
report.setWellFailed({CR::WellFailure::Type::Pressure, CR::Severity::NotANumber, dummy_component, baseif_.name()});
@ -170,7 +170,7 @@ getWellConvergence(const WellState<Scalar>& well_state,
// for BHP or THP controlled wells, we need to make sure the flow direction is correct
if (!well_is_stopped && baseif_.isPressureControlled(well_state)) {
// checking the flow direction
const double sign = baseif_.isProducer() ? -1. : 1.;
const Scalar sign = baseif_.isProducer() ? -1. : 1.;
const auto weight_total_flux = this->primary_variables_.getWQTotal() * sign;
constexpr int dummy_phase = -1;
if (weight_total_flux < 0.) {
@ -207,7 +207,7 @@ assembleAccelerationPressureLoss(const int seg,
const auto& segment_set = this->segmentSet();
// Add segment head
const double seg_area = segment_set[seg].crossArea();
const Scalar seg_area = segment_set[seg].crossArea();
const EvalWell signed_velocity_head = segments_.accelerationPressureLossContribution(seg, seg_area);
segments.pressure_drop_accel[seg] = signed_velocity_head.value();
@ -225,7 +225,8 @@ assembleAccelerationPressureLoss(const int seg,
// Subtract inlet head(s), opposite signs from above
for (const int inlet : segments_.inlets(seg)) {
// area used in formula is max of areas
const double inlet_area = std::max(seg_area, segment_set[inlet].crossArea());
const Scalar inlet_area = std::max(seg_area,
static_cast<Scalar>(segment_set[inlet].crossArea()));
const EvalWell signed_velocity_head_inlet = segments_.accelerationPressureLossContribution(inlet, inlet_area);
segments.pressure_drop_accel[seg] -= signed_velocity_head_inlet.value();
@ -437,7 +438,7 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
double residual = 0.;
Scalar residual = 0.;
if (eq_idx < baseif_.numComponents()) {
residual = std::abs(linSys_.residual()[seg][eq_idx]) * B_avg[eq_idx];
} else {
@ -459,7 +460,7 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
// handling the control equation residual
{
const double control_residual = std::abs(linSys_.residual()[0][numWellEq - 1]);
const Scalar control_residual = std::abs(linSys_.residual()[0][numWellEq - 1]);
if (std::isnan(control_residual) || std::isinf(control_residual)) {
deferred_logger.debug(fmt::format("nan or inf value for control residual for well {}", baseif_.name()));
return {false, residuals};
@ -471,14 +472,14 @@ getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
}
template<typename FluidSystem, typename Indices>
double
typename MultisegmentWellEval<FluidSystem,Indices>::Scalar
MultisegmentWellEval<FluidSystem,Indices>::
getControlTolerance(const WellState<Scalar>& well_state,
const double tolerance_wells,
const double tolerance_pressure_ms_wells,
const Scalar tolerance_wells,
const Scalar tolerance_pressure_ms_wells,
DeferredLogger& deferred_logger) const
{
double control_tolerance = 0.;
Scalar control_tolerance = 0.;
const int well_index = baseif_.indexOfWell();
const auto& ws = well_state.well(well_index);
@ -538,19 +539,19 @@ getControlTolerance(const WellState<Scalar>& well_state,
}
template<typename FluidSystem, typename Indices>
double
typename MultisegmentWellEval<FluidSystem,Indices>::Scalar
MultisegmentWellEval<FluidSystem,Indices>::
getResidualMeasureValue(const WellState<Scalar>& well_state,
const std::vector<double>& residuals,
const double tolerance_wells,
const double tolerance_pressure_ms_wells,
const std::vector<Scalar>& residuals,
const Scalar tolerance_wells,
const Scalar tolerance_pressure_ms_wells,
DeferredLogger& deferred_logger) const
{
assert(int(residuals.size()) == numWellEq + 1);
const double rate_tolerance = tolerance_wells;
const Scalar rate_tolerance = tolerance_wells;
[[maybe_unused]] int count = 0;
double sum = 0;
Scalar sum = 0;
for (int eq_idx = 0; eq_idx < numWellEq - 1; ++eq_idx) {
if (residuals[eq_idx] > rate_tolerance) {
sum += residuals[eq_idx] / rate_tolerance;
@ -558,13 +559,13 @@ getResidualMeasureValue(const WellState<Scalar>& well_state,
}
}
const double pressure_tolerance = tolerance_pressure_ms_wells;
const Scalar pressure_tolerance = tolerance_pressure_ms_wells;
if (residuals[SPres] > pressure_tolerance) {
sum += residuals[SPres] / pressure_tolerance;
++count;
}
const double control_tolerance = getControlTolerance(well_state,
const Scalar control_tolerance = getControlTolerance(well_state,
tolerance_wells,
tolerance_pressure_ms_wells,
deferred_logger);

33
opm/simulators/wells/MultisegmentWellEval.hpp
View File

@ -32,8 +32,7 @@
#include <utility>
#include <vector>
namespace Opm
{
namespace Opm {
class ConvergenceReport;
class Schedule;
@ -100,13 +99,13 @@ protected:
/// check whether the well equations get converged for this well
ConvergenceReport getWellConvergence(const WellState<Scalar>& well_state,
const std::vector<double>& B_avg,
const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger,
const double max_residual_allowed,
const double tolerance_wells,
const double relaxed_inner_tolerance_flow_ms_well,
const double tolerance_pressure_ms_wells,
const double relaxed_inner_tolerance_pressure_ms_well,
const Scalar max_residual_allowed,
const Scalar tolerance_wells,
const Scalar relaxed_inner_tolerance_flow_ms_well,
const Scalar tolerance_pressure_ms_wells,
const Scalar relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance,
const bool well_is_stopped) const;
@ -114,15 +113,15 @@ protected:
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const;
double getControlTolerance(const WellState<Scalar>& well_state,
const double tolerance_wells,
const double tolerance_pressure_ms_wells,
Scalar getControlTolerance(const WellState<Scalar>& well_state,
const Scalar tolerance_wells,
const Scalar tolerance_pressure_ms_wells,
DeferredLogger& deferred_logger) const;
double getResidualMeasureValue(const WellState<Scalar>& well_state,
const std::vector<double>& residuals,
const double tolerance_wells,
const double tolerance_pressure_ms_wells,
Scalar getResidualMeasureValue(const WellState<Scalar>& well_state,
const std::vector<Scalar>& residuals,
const Scalar tolerance_wells,
const Scalar tolerance_pressure_ms_wells,
DeferredLogger& deferred_logger) const;
void assembleAccelerationPressureLoss(const int seg,
@ -141,10 +140,10 @@ protected:
MSWSegments segments_; //!< Segment properties
// depth difference between perforations and the perforated grid cells
std::vector<double> cell_perforation_depth_diffs_;
std::vector<Scalar> cell_perforation_depth_diffs_;
// pressure correction due to the different depth of the perforation and
// center depth of the grid block
std::vector<double> cell_perforation_pressure_diffs_;
std::vector<Scalar> cell_perforation_pressure_diffs_;
};
}

42
opm/simulators/wells/MultisegmentWellGeneric.cpp
View File

@ -58,8 +58,8 @@ scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inle
auto& segments = ws.segments;
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 = ws.surface_rates[phase];
const Scalar unscaled_top_seg_rate = segment_rates[phase];
const Scalar well_phase_rate = ws.surface_rates[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;
@ -67,20 +67,20 @@ scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inle
} else {
// Due to various reasons, the well/top segment rate can be zero for this phase.
// We can not scale this rate directly. The following approach is used to initialize the segment rates.
double sumTw = 0;
Scalar sumTw = 0;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
sumTw += baseif_.wellIndex()[perf];
}
// only handling this specific phase
constexpr double num_single_phase = 1;
std::vector<double> perforation_rates(num_single_phase * baseif_.numPerfs(), 0.0);
const double perf_phaserate_scaled = ws.surface_rates[phase] / sumTw;
constexpr Scalar num_single_phase = 1;
std::vector<Scalar> perforation_rates(num_single_phase * baseif_.numPerfs(), 0.0);
const Scalar perf_phaserate_scaled = ws.surface_rates[phase] / sumTw;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
perforation_rates[perf] = baseif_.wellIndex()[perf] * perf_phaserate_scaled;
}
std::vector<double> rates;
std::vector<Scalar> rates;
WellState<Scalar>::calculateSegmentRates(segment_inlets,
segment_perforations,
perforation_rates,
@ -137,7 +137,7 @@ segmentNumberToIndex(const int segment_number) const
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
detectOscillations(const std::vector<double>& measure_history, bool& oscillate, bool& stagnate) const
detectOscillations(const std::vector<Scalar>& measure_history, bool& oscillate, bool& stagnate) const
{
const auto it = measure_history.size() - 1;
if ( it < 2 ) {
@ -147,16 +147,16 @@ detectOscillations(const std::vector<double>& measure_history, bool& oscillate,
}
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 Scalar F0 = measure_history[it];
const Scalar F1 = measure_history[it - 1];
const Scalar F2 = measure_history[it - 2];
const Scalar d1 = std::abs((F0 - F2) / F0);
const Scalar d2 = std::abs((F0 - F1) / F0);
const double oscillaton_rel_tol = 0.2;
const Scalar oscillaton_rel_tol = 0.2;
oscillate = (d1 < oscillaton_rel_tol) && (oscillaton_rel_tol < d2);
const double stagnation_rel_tol = 1.e-2;
const Scalar stagnation_rel_tol = 1.e-2;
stagnate = std::abs((F1 - F2) / F2) <= stagnation_rel_tol;
}
@ -179,15 +179,15 @@ accelerationalPressureLossConsidered() const
template<class Scalar>
double
Scalar
MultisegmentWellGeneric<Scalar>::getSegmentDp(const int seg,
const double density,
const std::vector<double>& seg_dp) const
const Scalar density,
const std::vector<Scalar>& seg_dp) const
{
const double segment_depth = this->segmentSet()[seg].depth();
const Scalar segment_depth = this->segmentSet()[seg].depth();
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const double segment_depth_outlet = seg == 0 ? baseif_.refDepth() : this->segmentSet()[outlet_segment_index].depth();
double dp = wellhelpers::computeHydrostaticCorrection(segment_depth_outlet, segment_depth,
const Scalar segment_depth_outlet = seg == 0 ? baseif_.refDepth() : this->segmentSet()[outlet_segment_index].depth();
Scalar dp = wellhelpers::computeHydrostaticCorrection(segment_depth_outlet, segment_depth,
density, baseif_.gravity());
// we add the hydrostatic correction from the outlet segment
// in order to get the correction all the way to the bhp ref depth.

8
opm/simulators/wells/MultisegmentWellGeneric.hpp
View File

@ -64,16 +64,16 @@ protected:
WellSegmentCompPressureDrop compPressureDrop() const;
/// Detect oscillation or stagnation based on the residual measure history
void detectOscillations(const std::vector<double>& measure_history,
void detectOscillations(const std::vector<Scalar>& measure_history,
bool& oscillate,
bool& stagnate) const;
bool accelerationalPressureLossConsidered() const;
bool frictionalPressureLossConsidered() const;
double getSegmentDp(const int seg,
const double density,
const std::vector<double>& seg_dp) const;
Scalar getSegmentDp(const int seg,
const Scalar density,
const std::vector<Scalar>& seg_dp) const;
const WellInterfaceGeneric<Scalar>& baseif_;
};

62
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
View File

@ -91,7 +91,7 @@ update(const WellState<Scalar>& well_state,
for (std::size_t seg = 0; seg < value_.size(); ++seg) {
// calculate the total rate for each segment
double total_seg_rate = 0.0;
Scalar total_seg_rate = 0.0;
// the segment pressure
value_[seg][SPres] = segment_pressure[seg];
// TODO: under what kind of circustances, the following will be wrong?
@ -102,9 +102,9 @@ update(const WellState<Scalar>& well_state,
if (seg == 0) {
if (well_.isInjector()) {
total_seg_rate = std::max(total_seg_rate, 0.);
total_seg_rate = std::max(total_seg_rate, Scalar{0.});
} else {
total_seg_rate = std::min(total_seg_rate, 0.);
total_seg_rate = std::min(total_seg_rate, Scalar{0.});
}
}
value_[seg][WQTotal] = total_seg_rate;
@ -157,23 +157,23 @@ update(const WellState<Scalar>& well_state,
template<class FluidSystem, class Indices>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
updateNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double dFLimit,
const Scalar relaxation_factor,
const Scalar dFLimit,
const bool stop_or_zero_rate_target,
const double max_pressure_change)
const Scalar max_pressure_change)
{
const std::vector<std::array<double, numWellEq>> old_primary_variables = value_;
const std::vector<std::array<Scalar, numWellEq>> old_primary_variables = value_;
for (std::size_t seg = 0; seg < value_.size(); ++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);
const Scalar dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
value_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
}
if (has_gfrac_variable) {
const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
const Scalar dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
value_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
}
@ -183,10 +183,10 @@ updateNewton(const BVectorWell& dwells,
// update the segment pressure
{
const int sign = dwells[seg][SPres] > 0.? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
const Scalar dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
// some cases might have defaulted bhp constraint of 1 bar, we use a slightly smaller value as the bhp lower limit for Newton update
// so that bhp constaint can be an active control when needed.
const double lower_limit = (seg == 0) ? bhp_lower_limit : seg_pres_lower_limit;
const Scalar lower_limit = (seg == 0) ? bhp_lower_limit : seg_pres_lower_limit;
value_[seg][SPres] = std::max(old_primary_variables[seg][SPres] - dx_limited, lower_limit);
}
@ -197,9 +197,9 @@ updateNewton(const BVectorWell& dwells,
// make sure that no injector produce and no producer inject
if (seg == 0) {
if (well_.isInjector()) {
value_[seg][WQTotal] = std::max(value_[seg][WQTotal], 0.0);
value_[seg][WQTotal] = std::max(value_[seg][WQTotal], Scalar{0.0});
} else {
value_[seg][WQTotal] = std::min(value_[seg][WQTotal], 0.0);
value_[seg][WQTotal] = std::min(value_[seg][WQTotal], Scalar{0.0});
}
}
}
@ -213,7 +213,7 @@ updateNewton(const BVectorWell& dwells,
template<class FluidSystem, class Indices>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices>::
copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const double rho,
const Scalar rho,
const bool stop_or_zero_rate_target,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
@ -236,7 +236,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
auto& vapoil = segments.vaporized_oil_rate;
auto& segment_pressure = segments.pressure;
for (std::size_t seg = 0; seg < value_.size(); ++seg) {
std::vector<double> fractions(well_.numPhases(), 0.0);
std::vector<Scalar> fractions(well_.numPhases(), 0.0);
fractions[oil_pos] = 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@ -253,7 +253,7 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
// convert the fractions to be Q_p / G_total to calculate the phase rates
for (int p = 0; p < well_.numPhases(); ++p) {
const double scale = well_.scalingFactor(p);
const Scalar scale = well_.scalingFactor(p);
// for injection wells, there should only one non-zero scaling factor
if (scale > 0.) {
fractions[p] /= scale;
@ -264,9 +264,9 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
}
// calculate the phase rates based on the primary variables
const double g_total = value_[seg][WQTotal];
const Scalar g_total = value_[seg][WQTotal];
for (int p = 0; p < well_.numPhases(); ++p) {
const double phase_rate = g_total * fractions[p];
const Scalar phase_rate = g_total * fractions[p];
segment_rates[seg * well_.numPhases() + p] = phase_rate;
if (seg == 0) { // top segment
ws.surface_rates[p] = phase_rate;
@ -282,12 +282,12 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
// Calculate other per-phase dynamic quantities.
const auto temperature = 0.0; // Ignore thermal effects
const auto saltConc = 0.0; // Ignore salt precipitation
const auto Rvw = 0.0; // Ignore vaporised water.
const Scalar temperature = 0.0; // Ignore thermal effects
const Scalar saltConc = 0.0; // Ignore salt precipitation
const Scalar Rvw = 0.0; // Ignore vaporised water.
auto rsMax = 0.0;
auto rvMax = 0.0;
Scalar rsMax = 0.0;
Scalar rvMax = 0.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
// Both oil and gas active.
rsMax = FluidSystem::oilPvt()
@ -322,10 +322,10 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
well_.rateConverter().calcReservoirVoidageRates
(pvtReg, segment_pressure[seg],
std::max(0.0, Rs),
std::max(0.0, Rv),
0.0, // Rsw
0.0, // Rvw
std::max(Scalar{0.0}, Rs),
std::max(Scalar{0.0}, Rv),
Scalar{0.0}, // Rsw
Scalar{0.0}, // Rvw
temperature, saltConc, surf_rates, resv_rates);
}
@ -361,7 +361,9 @@ copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
segments.phase_viscosity[seg * well_.numPhases() + pu.phase_pos[Water]] =
FluidSystem::waterPvt().viscosity(pvtReg, temperature, segment_pressure[seg], 0.0 /*Rsw*/, saltConc);
FluidSystem::waterPvt().viscosity(pvtReg, temperature,
segment_pressure[seg],
Scalar{0.0} /*Rsw*/, saltConc);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
@ -421,7 +423,7 @@ processFractions(const int seg)
const PhaseUsage& pu = well_.phaseUsage();
std::vector<double> fractions(well_.numPhases(), 0.0);
std::vector<Scalar> fractions(well_.numPhases(), 0.0);
assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
const int oil_pos = pu.phase_pos[Oil];
@ -524,7 +526,7 @@ volumeFractionScaled(const int seg,
// For reservoir rate control, the distr in well control is used for the
// rate conversion coefficients. For the injection well, only the distr of the injection
// phase is not zero.
const double scale = well_.scalingFactor(well_.modelCompIdxToFlowCompIdx(comp_idx));
const Scalar scale = well_.scalingFactor(well_.modelCompIdxToFlowCompIdx(comp_idx));
if (scale > 0.) {
return this->volumeFraction(seg, comp_idx) / scale;
}

12
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
View File

@ -76,7 +76,7 @@ public:
static constexpr int numWellEq = Indices::numPhases + 1;
using Scalar = typename FluidSystem::Scalar;
using EvalWell = DenseAd::Evaluation<double, /*size=*/Indices::numEq + numWellEq>;
using EvalWell = DenseAd::Evaluation<Scalar, /*size=*/Indices::numEq + numWellEq>;
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
using BVectorWell = typename Equations::BVectorWell;
@ -97,14 +97,14 @@ public:
//! \brief Update values from newton update vector.
void updateNewton(const BVectorWell& dwells,
const double relaxation_factor,
const double DFLimit,
const Scalar relaxation_factor,
const Scalar DFLimit,
const bool stop_or_zero_rate_target,
const double max_pressure_change);
const Scalar max_pressure_change);
//! \brief Copy values to well state.
void copyToWellState(const MultisegmentWellGeneric<Scalar>& mswell,
const double rho,
const Scalar rho,
const bool stop_or_zero_rate_target,
WellState<Scalar>& well_state,
const SummaryState& summary_state,
@ -166,7 +166,7 @@ private:
//! \brief The values for the primary variables
//! \details Based on different solution strategies, the wells can have different primary variables
std::vector<std::array<double, numWellEq>> value_;
std::vector<std::array<Scalar, numWellEq>> value_;
//! \brief The Evaluation for the well primary variables.
//! \details Contains derivatives and are used in AD calculation

64
opm/simulators/wells/MultisegmentWellSegments.cpp
View File

@ -100,7 +100,7 @@ MultisegmentWellSegments(const int numSegments,
}
perforations_[segment_index].push_back(i_perf_wells);
well.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segment_set[segment_index].depth();
const Scalar segment_depth = segment_set[segment_index].depth();
perforation_depth_diffs_[i_perf_wells] = well_.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
@ -120,10 +120,10 @@ MultisegmentWellSegments(const int numSegments,
// 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 < numSegments; ++seg) {
const double segment_depth = segment_set[seg].depth();
const Scalar segment_depth = segment_set[seg].depth();
const int outlet_segment_number = segment_set[seg].outletSegment();
const Segment& outlet_segment = segment_set[segment_set.segmentNumberToIndex(outlet_segment_number)];
const double outlet_depth = outlet_segment.depth();
const Scalar outlet_depth = outlet_segment.depth();
depth_diffs_[seg] = segment_depth - outlet_depth;
}
}
@ -136,7 +136,7 @@ computeFluidProperties(const EvalWell& temperature,
int pvt_region_index,
DeferredLogger& deferred_logger)
{
std::vector<double> surf_dens(well_.numComponents());
std::vector<Scalar> surf_dens(well_.numComponents());
// Surface density.
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@ -490,7 +490,7 @@ getSurfaceVolume(const EvalWell& temperature,
}
// We increase the segment volume with a factor 10 to stabilize the system.
const double volume = well_.wellEcl().getSegments()[seg_idx].volume();
const Scalar volume = well_.wellEcl().getSegments()[seg_idx].volume();
return volume / vol_ratio;
}
@ -536,13 +536,13 @@ getFrictionPressureLoss(const int seg,
const auto& segment_set = well_.wellEcl().getSegments();
const int outlet_segment_index = segment_set.segmentNumberToIndex(segment_set[seg].outletSegment());
const double length = segment_set[seg].totalLength() - segment_set[outlet_segment_index].totalLength();
const Scalar length = segment_set[seg].totalLength() - segment_set[outlet_segment_index].totalLength();
assert(length > 0.);
const double roughness = segment_set[seg].roughness();
const double area = segment_set[seg].crossArea();
const double diameter = segment_set[seg].internalDiameter();
const Scalar roughness = segment_set[seg].roughness();
const Scalar area = segment_set[seg].crossArea();
const Scalar diameter = segment_set[seg].internalDiameter();
const double sign = mass_rate < 0. ? 1.0 : - 1.0;
const Scalar sign = mass_rate < 0. ? 1.0 : - 1.0;
return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
}
@ -632,11 +632,11 @@ pressureDropSpiralICD(const int seg,
const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
const double viscosity_cali = sicd.viscosityCalibration();
const Scalar viscosity_cali = sicd.viscosityCalibration();
using MathTool = MathToolbox<EvalWell>;
const double density_cali = sicd.densityCalibration();
const Scalar density_cali = sicd.densityCalibration();
// make sure we don't pass negative base to powers
const EvalWell temp_value1 = density > 0.0 ? MathTool::pow(density / density_cali, 0.75) : 0.0;
const EvalWell temp_value2 = mixture_viscosity > 0.0 ? MathTool::pow(mixture_viscosity / viscosity_cali, 0.25) : 0.0;
@ -644,9 +644,9 @@ pressureDropSpiralICD(const int seg,
// formulation before 2016, base_strength is used
// const double base_strength = sicd.strength() / density_cali;
// formulation since 2016, strength is used instead
const double strength = sicd.strength();
const Scalar strength = sicd.strength();
const double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0;
const Scalar sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0;
return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd;
}
@ -737,7 +737,7 @@ pressureDropAutoICD(const int seg,
using MathTool = MathToolbox<EvalWell>;
// make sure we don't pass negative base to powers
auto safe_pow = [](const auto& a, const double b) {
auto safe_pow = [](const auto& a, const Scalar b) {
return a > 0.0 ? MathTool::pow(a,b) : 0.0;
};
@ -749,13 +749,13 @@ pressureDropAutoICD(const int seg,
+ safe_pow(oil_fraction, aicd.oilDensityExponent()) * oil_density
+ safe_pow(gas_fraction, aicd.gasDensityExponent()) * gas_density;
const double rho_reference = aicd.densityCalibration();
const double visc_reference = aicd.viscosityCalibration();
const Scalar rho_reference = aicd.densityCalibration();
const Scalar visc_reference = aicd.viscosityCalibration();
const auto volume_rate_icd = mass_rate * aicd.scalingFactor() / mixture_density;
const double sign = volume_rate_icd <= 0. ? 1.0 : -1.0;
const Scalar sign = volume_rate_icd <= 0. ? 1.0 : -1.0;
// convert 1 unit volume rate
using M = UnitSystem::measure;
const double unit_volume_rate = unit_system.to_si(M::geometric_volume_rate, 1.);
const Scalar unit_volume_rate = unit_system.to_si(M::geometric_volume_rate, 1.);
// TODO: we did not consider the maximum allowed rate here
const auto result = sign / rho_reference * mixture_density * mixture_density
@ -807,22 +807,22 @@ pressureDropValve(const int seg,
}
}
const double additional_length = valve.pipeAdditionalLength();
const double roughness = valve.pipeRoughness();
const double diameter = valve.pipeDiameter();
const double area = valve.pipeCrossArea();
const Scalar additional_length = valve.pipeAdditionalLength();
const Scalar roughness = valve.pipeRoughness();
const Scalar diameter = valve.pipeDiameter();
const Scalar area = valve.pipeCrossArea();
const EvalWell friction_pressure_loss =
mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc);
const ValveUDAEval uda_eval {summary_state, this->well_.name(), static_cast<std::size_t>(segment_set[seg].segmentNumber())};
const double area_con = valve.conCrossArea(uda_eval);
const double cv = valve.conFlowCoefficient();
const Scalar area_con = valve.conCrossArea(uda_eval);
const Scalar cv = valve.conFlowCoefficient();
const EvalWell constriction_pressure_loss =
mswellhelpers::valveContrictionPressureLoss(mass_rate, density, area_con, cv);
const double sign = mass_rate <= 0. ? 1.0 : -1.0;
const Scalar sign = mass_rate <= 0. ? 1.0 : -1.0;
return sign * (friction_pressure_loss + constriction_pressure_loss);
}
@ -830,7 +830,7 @@ template<class FluidSystem, class Indices>
typename MultisegmentWellSegments<FluidSystem,Indices>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices>::
accelerationPressureLossContribution(const int seg,
const double area,
const Scalar area,
const bool extra_reverse_flow_derivatives /*false*/) const
{
// Compute the *signed* velocity head for given segment (sign is positive for flow towards surface, i.e., negative rate)
@ -856,7 +856,7 @@ accelerationPressureLossContribution(const int seg,
mass_rate.clearDerivatives();
}
}
const double sign = mass_rate > 0 ? -1.0 : 1.0;
const Scalar sign = mass_rate > 0 ? -1.0 : 1.0;
return sign*mswellhelpers::velocityHead(area, mass_rate, density);
}
@ -900,7 +900,7 @@ void
MultisegmentWellSegments<FluidSystem,Indices>::
copyPhaseDensities(const unsigned phaseIdx,
const std::size_t stride,
double* dens) const
Scalar* dens) const
{
const auto compIdx = Indices::canonicalToActiveComponentIndex
(FluidSystem::solventComponentIndex(phaseIdx));
@ -912,7 +912,7 @@ copyPhaseDensities(const unsigned phaseIdx,
}
template <class FluidSystem, class Indices>
double
typename MultisegmentWellSegments<FluidSystem,Indices>::Scalar
MultisegmentWellSegments<FluidSystem,Indices>::
mixtureDensity(const int seg) const
{
@ -941,7 +941,7 @@ mixtureDensity(const int seg) const
}
template <class FluidSystem, class Indices>
double
typename MultisegmentWellSegments<FluidSystem,Indices>::Scalar
MultisegmentWellSegments<FluidSystem,Indices>::
mixtureDensityWithExponents(const int seg) const
{
@ -957,7 +957,7 @@ mixtureDensityWithExponents(const int seg) const
}
template <class FluidSystem, class Indices>
double
typename MultisegmentWellSegments<FluidSystem,Indices>::Scalar
MultisegmentWellSegments<FluidSystem,Indices>::
mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const
{

10
opm/simulators/wells/MultisegmentWellSegments.hpp
View File

@ -92,7 +92,7 @@ public:
// pressure loss contribution due to acceleration
EvalWell accelerationPressureLossContribution(const int seg,
const double area,
const Scalar area,
const bool extra_reverse_flow_derivatives = false) const;
const std::vector<std::vector<int>>& inlets() const
@ -176,11 +176,11 @@ private:
void copyPhaseDensities(const unsigned phaseIdx,
const std::size_t stride,
double* dens) const;
Scalar* dens) const;
double mixtureDensity(const int seg) const;
double mixtureDensityWithExponents(const int seg) const;
double mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const;
Scalar mixtureDensity(const int seg) const;
Scalar mixtureDensityWithExponents(const int seg) const;
Scalar mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const;
};
} // namespace Opm

209
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

@ -69,7 +69,7 @@ namespace Opm
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
, MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices>&>(*this))
, regularize_(false)
, segment_fluid_initial_(this->numberOfSegments(), std::vector<double>(this->num_components_, 0.0))
, segment_fluid_initial_(this->numberOfSegments(), std::vector<Scalar>(this->num_components_, 0.0))
{
// not handling solvent or polymer for now with multisegment well
if constexpr (has_solvent) {
@ -116,8 +116,8 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
init(const PhaseUsage* phase_usage_arg,
const std::vector<double>& depth_arg,
const double gravity_arg,
const std::vector<Scalar>& depth_arg,
const Scalar gravity_arg,
const int num_cells,
const std::vector< Scalar >& B_avg,
const bool changed_to_open_this_step)
@ -201,7 +201,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
getWellConvergence(const SummaryState& /* summary_state */,
const WellState<Scalar>& well_state,
const std::vector<double>& B_avg,
const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger,
const bool relax_tolerance) const
{
@ -283,7 +283,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
computeWellPotentials(const Simulator& simulator,
const WellState<Scalar>& well_state,
std::vector<double>& well_potentials,
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger)
{
const auto [compute_potential, bhp_controlled_well] =
@ -327,7 +327,7 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
computeWellRatesAtBhpLimit(const Simulator& simulator,
std::vector<double>& well_flux,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const
{
if (this->well_ecl_.isInjector()) {
@ -343,11 +343,10 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
computeWellRatesWithBhp(const Simulator& simulator,
const double& bhp,
std::vector<double>& well_flux,
const Scalar& bhp,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const
{
const int np = this->number_of_phases_;
well_flux.resize(np, 0.0);
@ -365,7 +364,7 @@ namespace Opm
// flux for each perforation
std::vector<Scalar> mob(this->num_components_, 0.);
getMobility(simulator, perf, mob, deferred_logger);
const double trans_mult = simulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(intQuants, cell_idx);
const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
const Scalar seg_pressure = segment_pressure[seg];
@ -389,7 +388,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
computeWellRatesWithBhpIterations(const Simulator& simulator,
const Scalar& bhp,
std::vector<double>& well_flux,
std::vector<Scalar>& well_flux,
DeferredLogger& deferred_logger) const
{
// creating a copy of the well itself, to avoid messing up the explicit information
@ -429,7 +428,7 @@ namespace Opm
trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
}
if (!trivial) {
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase) {
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
@ -457,21 +456,22 @@ namespace Opm
template<typename TypeTag>
std::vector<double>
std::vector<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
computeWellPotentialWithTHP(const WellState<Scalar>& well_state,
const Simulator& simulator,
DeferredLogger& deferred_logger) const
{
std::vector<double> potentials(this->number_of_phases_, 0.0);
std::vector<Scalar> potentials(this->number_of_phases_, 0.0);
const auto& summary_state = simulator.vanguard().summaryState();
const auto& well = this->well_ecl_;
if (well.isInjector()){
if (well.isInjector()) {
auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = well.injectionControls(summary_state);
const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
const Scalar bhp = std::min(*bhp_at_thp_limit,
static_cast<Scalar>(controls.bhp_limit));
computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
deferred_logger.debug("Converged thp based potential calculation for well "
+ this->name() + ", at bhp = " + std::to_string(bhp));
@ -480,7 +480,7 @@ namespace Opm
"Failed in getting converged thp based potential calculation for well "
+ this->name() + ". Instead the bhp based value is used");
const auto& controls = well.injectionControls(summary_state);
const double bhp = controls.bhp_limit;
const Scalar bhp = controls.bhp_limit;
computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
}
} else {
@ -488,7 +488,8 @@ namespace Opm
well_state, simulator, summary_state, deferred_logger);
if (bhp_at_thp_limit) {
const auto& controls = well.productionControls(summary_state);
const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
const Scalar bhp = std::max(*bhp_at_thp_limit,
static_cast<Scalar>(controls.bhp_limit));
computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
deferred_logger.debug("Converged thp based potential calculation for well "
+ this->name() + ", at bhp = " + std::to_string(bhp));
@ -497,7 +498,7 @@ namespace Opm
"Failed in getting converged thp based potential calculation for well "
+ this->name() + ". Instead the bhp based value is used");
const auto& controls = well.productionControls(summary_state);
const double bhp = controls.bhp_limit;
const Scalar bhp = controls.bhp_limit;
computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger);
}
}
@ -509,7 +510,7 @@ namespace Opm
bool
MultisegmentWell<TypeTag>::
computeWellPotentialsImplicit(const Simulator& simulator,
std::vector<double>& well_potentials,
std::vector<Scalar>& well_potentials,
DeferredLogger& deferred_logger) const
{
// Create a copy of the well.
@ -544,7 +545,7 @@ namespace Opm
trivial = trivial && (ws.well_potentials[phase] == 0.0) ;
}
if (!trivial) {
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
const Scalar sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase) {
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
@ -611,14 +612,14 @@ namespace Opm
{
for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
std::vector<double> kr(this->number_of_phases_, 0.0);
std::vector<double> density(this->number_of_phases_, 0.0);
std::vector<Scalar> kr(this->number_of_phases_, 0.0);
std::vector<Scalar> density(this->number_of_phases_, 0.0);
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
const auto& fs = intQuants.fluidState();
double sum_kr = 0.;
Scalar sum_kr = 0.;
const PhaseUsage& pu = this->phaseUsage();
if (pu.phase_used[Water]) {
@ -645,7 +646,7 @@ namespace Opm
assert(sum_kr != 0.);
// calculate the average density
double average_density = 0.;
Scalar average_density = 0.;
for (int p = 0; p < this->number_of_phases_; ++p) {
average_density += kr[p] * density[p];
}
@ -666,7 +667,7 @@ namespace Opm
{
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// TODO: trying to reduce the times for the surfaceVolumeFraction calculation
const double surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
}
@ -684,12 +685,12 @@ namespace Opm
const BVectorWell& dwells,
WellState<Scalar>& well_state,
DeferredLogger& deferred_logger,
const double relaxation_factor)
const Scalar relaxation_factor)
{
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
const double dFLimit = this->param_.dwell_fraction_max_;
const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
const Scalar dFLimit = this->param_.dwell_fraction_max_;
const Scalar max_pressure_change = this->param_.max_pressure_change_ms_wells_;
const bool stop_or_zero_rate_target = this->stopppedOrZeroRateTarget(summary_state, well_state);
this->primary_variables_.updateNewton(dwells,
relaxation_factor,
@ -746,12 +747,12 @@ namespace Opm
};
const int np = this->number_of_phases_;
auto setToZero = [np](double* x) -> void
auto setToZero = [np](Scalar* x) -> void
{
std::fill_n(x, np, 0.0);
};
auto addVector = [np](const double* src, double* dest) -> void
auto addVector = [np](const Scalar* src, Scalar* dest) -> void
{
std::transform(src, src + np, dest, dest, std::plus<>{});
};
@ -769,7 +770,7 @@ namespace Opm
for ( const auto& perf : *this->perf_data_){
auto allPerfID = perf.ecl_index;
auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
auto connPICalc = [&wellPICalc, allPerfID](const Scalar mobility) -> Scalar
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
@ -803,7 +804,7 @@ namespace Opm
template<typename TypeTag>
double
typename MultisegmentWell<TypeTag>::Scalar
MultisegmentWell<TypeTag>::
connectionDensity(const int globalConnIdx,
[[maybe_unused]] const int openConnIdx) const
@ -877,7 +878,7 @@ namespace Opm
const Value perf_seg_press_diff = this->gravity() * segment_density *
this->segments_.perforation_depth_diff(perf);
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
// perforation pressure is the wellbore pressure corrected to perforation depth
// (positive sign due to convention in segments_.perforation_depth_diff() )
@ -982,7 +983,7 @@ namespace Opm
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs);
const Scalar d = 1.0 - getValue(rv) * getValue(rs);
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_rates.vap_oil = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
@ -1132,11 +1133,11 @@ namespace Opm
// from the reference results, it looks like MSW uses segment pressure instead of BHP here
// Note: this is against the documented definition.
// we can change this depending on what we want
const double segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
const double perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
const Scalar segment_pres = this->primary_variables_.getSegmentPressure(seg).value();
const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value()
* this->segments_.perforation_depth_diff(perf);
const double perf_press = segment_pres + perf_seg_press_diff;
const double multiplier = this->getInjMult(perf, segment_pres, perf_press);
const Scalar perf_press = segment_pres + perf_seg_press_diff;
const Scalar multiplier = this->getInjMult(perf, segment_pres, perf_press);
for (std::size_t i = 0; i < mob.size(); ++i) {
mob[i] *= multiplier;
}
@ -1146,7 +1147,7 @@ namespace Opm
template<typename TypeTag>
double
typename MultisegmentWell<TypeTag>::Scalar
MultisegmentWell<TypeTag>::
getRefDensity() const
{
@ -1161,14 +1162,14 @@ namespace Opm
DeferredLogger& deferred_logger)
{
const auto& summaryState = simulator.vanguard().summaryState();
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
// Crude but works: default is one atmosphere.
// TODO: a better way to detect whether the BHP is defaulted or not
const bool bhp_limit_not_defaulted = bhp_limit > 1.5 * unit::barsa;
if ( bhp_limit_not_defaulted || !this->wellHasTHPConstraints(summaryState) ) {
// if the BHP limit is not defaulted or the well does not have a THP limit
// we need to check the BHP limit
double total_ipr_mass_rate = 0.0;
Scalar total_ipr_mass_rate = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
{
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@ -1176,9 +1177,9 @@ namespace Opm
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
const double ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
const Scalar ipr_rate = this->ipr_a_[compIdx] - this->ipr_b_[compIdx] * bhp_limit;
const double rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
const Scalar rho = FluidSystem::referenceDensity( phaseIdx, Base::pvtRegionIdx() );
total_ipr_mass_rate += ipr_rate * rho;
}
if ( (this->isProducer() && total_ipr_mass_rate < 0.) || (this->isInjector() && total_ipr_mass_rate > 0.) ) {
@ -1190,11 +1191,11 @@ namespace Opm
// option 1: calculate well rates based on the BHP limit.
// option 2: stick with the above IPR curve
// we use IPR here
std::vector<double> well_rates_bhp_limit;
std::vector<Scalar> well_rates_bhp_limit;
computeWellRatesWithBhp(simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
const double thp_limit = this->getTHPConstraint(summaryState);
const double thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
const Scalar thp_limit = this->getTHPConstraint(summaryState);
const Scalar thp = WellBhpThpCalculator(*this).calculateThpFromBhp(well_rates_bhp_limit,
bhp_limit,
this->getRefDensity(),
this->wellEcl().alq_value(summaryState),
@ -1231,10 +1232,10 @@ namespace Opm
std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
const int nseg = this->numberOfSegments();
std::vector<double> seg_dp(nseg, 0.0);
std::vector<Scalar> seg_dp(nseg, 0.0);
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const double dp = this->getSegmentDp(seg,
const Scalar dp = this->getSegmentDp(seg,
this->segments_.density(seg).value(),
seg_dp);
seg_dp[seg] = dp;
@ -1248,13 +1249,13 @@ namespace Opm
const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
const auto& fs = int_quantities.fluidState();
// pressure difference between the segment and the perforation
const double perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
const Scalar perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = this->getPerfCellPressure(fs).value();
const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
std::vector<Scalar> b_perf(this->num_components_);
for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
@ -1264,8 +1265,8 @@ namespace Opm
}
// the pressure difference between the connection and BHP
const double h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
const double pressure_diff = pressure_cell - h_perf;
const Scalar h_perf = cell_perf_press_diff + perf_seg_press_diff + dp;
const Scalar pressure_diff = pressure_cell - h_perf;
// do not take into consideration the crossflow here.
if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
@ -1274,13 +1275,13 @@ namespace Opm
}
// the well index associated with the connection
const double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quantities, cell_idx);
const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> ipr_b_perf(this->ipr_b_.size());
std::vector<Scalar> ipr_a_perf(this->ipr_a_.size());
std::vector<Scalar> ipr_b_perf(this->ipr_b_.size());
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
const Scalar tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
@ -1289,17 +1290,17 @@ namespace Opm
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned oil_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const double rs = (fs.Rs()).value();
const double rv = (fs.Rv()).value();
const Scalar rs = (fs.Rs()).value();
const Scalar rv = (fs.Rv()).value();
const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
const Scalar dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const Scalar vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
ipr_a_perf[gas_comp_idx] += dis_gas_a;
ipr_a_perf[oil_comp_idx] += vap_oil_a;
const double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const double vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
const Scalar dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const Scalar vap_oil_b = rv * ipr_b_perf[gas_comp_idx];
ipr_b_perf[gas_comp_idx] += dis_gas_b;
ipr_b_perf[oil_comp_idx] += vap_oil_b;
@ -1398,10 +1399,10 @@ namespace Opm
if (obtain_bhp) {
this->operability_status_.can_obtain_bhp_with_thp_limit = true;
const double bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summaryState);
this->operability_status_.obey_bhp_limit_with_thp_limit = (*obtain_bhp >= bhp_limit);
const double thp_limit = this->getTHPConstraint(summaryState);
const Scalar thp_limit = this->getTHPConstraint(summaryState);
if (this->isProducer() && *obtain_bhp < thp_limit) {
const std::string msg = " obtained bhp " + std::to_string(unit::convert::to(*obtain_bhp, unit::barsa))
+ " bars is SMALLER than thp limit "
@ -1422,7 +1423,7 @@ namespace Opm
this->operability_status_.can_obtain_bhp_with_thp_limit = false;
this->operability_status_.obey_bhp_limit_with_thp_limit = false;
if (!this->wellIsStopped()) {
const double thp_limit = this->getTHPConstraint(summaryState);
const Scalar thp_limit = this->getTHPConstraint(summaryState);
deferred_logger.debug(" could not find bhp value at thp limit "
+ std::to_string(unit::convert::to(thp_limit, unit::barsa))
+ " bar for well " + this->name() + ", the well might need to be closed ");
@ -1457,11 +1458,11 @@ namespace Opm
}
std::vector<std::vector<Scalar> > residual_history;
std::vector<double> measure_history;
std::vector<Scalar> measure_history;
int it = 0;
// relaxation factor
double relaxation_factor = 1.;
const double min_relaxation_factor = 0.6;
Scalar relaxation_factor = 1.;
const Scalar min_relaxation_factor = 0.6;
bool converged = false;
int stagnate_count = 0;
bool relax_convergence = false;
@ -1536,7 +1537,7 @@ namespace Opm
}
// a factor value to reduce the relaxation_factor
const double reduction_mutliplier = 0.9;
const Scalar reduction_mutliplier = 0.9;
relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
// debug output
@ -1609,11 +1610,11 @@ namespace Opm
}
std::vector<std::vector<Scalar> > residual_history;
std::vector<double> measure_history;
std::vector<Scalar> measure_history;
int it = 0;
// relaxation factor
double relaxation_factor = 1.;
const double min_relaxation_factor = 0.6;
Scalar relaxation_factor = 1.;
const Scalar min_relaxation_factor = 0.6;
bool converged = false;
[[maybe_unused]] int stagnate_count = 0;
bool relax_convergence = false;
@ -1645,8 +1646,10 @@ namespace Opm
for (; it < max_iter_number; ++it, ++debug_cost_counter_) {
its_since_last_switch++;
if (allow_switching && its_since_last_switch >= min_its_after_switch){
const double wqTotal = this->primary_variables_.getWQTotal().value();
changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state, inj_controls, prod_controls, wqTotal, deferred_logger, fixed_control, fixed_status);
const Scalar wqTotal = this->primary_variables_.getWQTotal().value();
changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state,
inj_controls, prod_controls, wqTotal,
deferred_logger, fixed_control, fixed_status);
if (changed){
its_since_last_switch = 0;
switch_count++;
@ -1726,7 +1729,7 @@ namespace Opm
}
// a factor value to reduce the relaxation_factor
constexpr double reduction_mutliplier = 0.9;
constexpr Scalar reduction_mutliplier = 0.9;
relaxation_factor = std::max(relaxation_factor * reduction_mutliplier, min_relaxation_factor);
// debug output
@ -1869,7 +1872,7 @@ namespace Opm
const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<EvalWell> mob(this->num_components_, 0.0);
getMobility(simulator, perf, mob, deferred_logger);
const double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
std::vector<EvalWell> cq_s(this->num_components_, 0.0);
@ -1959,10 +1962,10 @@ namespace Opm
// pressure difference between the segment and the perforation
const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf);
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = this->getPerfCellPressure(fs).value();
const double perf_press = pressure_cell - cell_perf_press_diff;
const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
const Scalar perf_press = pressure_cell - cell_perf_press_diff;
// Pressure drawdown (also used to determine direction of flow)
// TODO: not 100% sure about the sign of the seg_perf_press_diff
const EvalWell drawdown = perf_press - (segment_pressure + perf_seg_press_diff);
@ -2023,7 +2026,7 @@ namespace Opm
template<typename TypeTag>
std::optional<double>
std::optional<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitProd(const WellState<Scalar>& well_state,
const Simulator& simulator,
@ -2040,21 +2043,21 @@ namespace Opm
template<typename TypeTag>
std::optional<double>
std::optional<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitProdWithAlq(const Simulator& simulator,
const SummaryState& summary_state,
const double alq_value,
const Scalar alq_value,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &simulator, &deferred_logger](const double bhp) {
auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
// Not solving the well equations here, which means we are
// calculating at the current Fg/Fw values of the
// well. This does not matter unless the well is
// crossflowing, and then it is likely still a good
// approximation.
std::vector<double> rates(3);
std::vector<Scalar> rates(3);
computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
return rates;
};
@ -2071,11 +2074,11 @@ namespace Opm
if (bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &simulator, &deferred_logger](const double bhp) {
auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
// Solver the well iterations to see if we are
// able to get a solution with an update
// solution
std::vector<double> rates(3);
std::vector<Scalar> rates(3);
computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
return rates;
};
@ -2091,20 +2094,20 @@ namespace Opm
}
template<typename TypeTag>
std::optional<double>
std::optional<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
computeBhpAtThpLimitInj(const Simulator& simulator,
const SummaryState& summary_state,
DeferredLogger& deferred_logger) const
{
// Make the frates() function.
auto frates = [this, &simulator, &deferred_logger](const double bhp) {
auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) {
// Not solving the well equations here, which means we are
// calculating at the current Fg/Fw values of the
// well. This does not matter unless the well is
// crossflowing, and then it is likely still a good
// approximation.
std::vector<double> rates(3);
std::vector<Scalar> rates(3);
computeWellRatesWithBhp(simulator, bhp, rates, deferred_logger);
return rates;
};
@ -2121,11 +2124,11 @@ namespace Opm
if (bhpAtLimit)
return bhpAtLimit;
auto fratesIter = [this, &simulator, &deferred_logger](const double bhp) {
auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) {
// Solver the well iterations to see if we are
// able to get a solution with an update
// solution
std::vector<double> rates(3);
std::vector<Scalar> rates(3);
computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger);
return rates;
};
@ -2145,18 +2148,18 @@ namespace Opm
template<typename TypeTag>
double
typename MultisegmentWell<TypeTag>::Scalar
MultisegmentWell<TypeTag>::
maxPerfPress(const Simulator& simulator) const
{
double max_pressure = 0.0;
Scalar max_pressure = 0.0;
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
const auto& fs = int_quants.fluidState();
double pressure_cell = this->getPerfCellPressure(fs).value();
Scalar pressure_cell = this->getPerfCellPressure(fs).value();
max_pressure = std::max(max_pressure, pressure_cell);
}
}
@ -2168,7 +2171,7 @@ namespace Opm
template<typename TypeTag>
std::vector<double>
std::vector<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
computeCurrentWellRates(const Simulator& simulator,
DeferredLogger& deferred_logger) const
@ -2185,7 +2188,7 @@ namespace Opm
const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<Scalar> mob(this->num_components_, 0.0);
getMobility(simulator, perf, mob, deferred_logger);
const double trans_mult = simulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
const Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quants, cell_idx);
const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
const std::vector<Scalar> Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
std::vector<Scalar> cq_s(this->num_components_, 0.0);
@ -2203,13 +2206,13 @@ namespace Opm
template <typename TypeTag>
std::vector<double>
std::vector<typename MultisegmentWell<TypeTag>::Scalar>
MultisegmentWell<TypeTag>::
getPrimaryVars() const
{
const int num_seg = this->numberOfSegments();
constexpr int num_eq = MSWEval::numWellEq;
std::vector<double> retval(num_seg * num_eq);
std::vector<Scalar> retval(num_seg * num_eq);
for (int ii = 0; ii < num_seg; ++ii) {
const auto& pv = this->primary_variables_.value(ii);
std::copy(pv.begin(), pv.end(), retval.begin() + ii * num_eq);
@ -2223,11 +2226,11 @@ namespace Opm
template <typename TypeTag>
int
MultisegmentWell<TypeTag>::
setPrimaryVars(std::vector<double>::const_iterator it)
setPrimaryVars(typename std::vector<Scalar>::const_iterator it)
{
const int num_seg = this->numberOfSegments();
constexpr int num_eq = MSWEval::numWellEq;
std::array<double, num_eq> tmp;
std::array<Scalar, num_eq> tmp;
for (int ii = 0; ii < num_seg; ++ii) {
const auto start = it + num_seg * num_eq;
std::copy(start, start + num_eq, tmp.begin());