OilfieldToolsNet/opm-simulators
4
0 Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-10-25 15:10:09 -05:00
Code Issues Packages Projects Releases Wiki Activity

changing numComponents in wells to be numConservationQuantities

Browse Source 
numComponents has been a very misleading name and something causes
confusion when coming to phase/component indices.

For example, mobility is not way to be associated with compoent.
This commit is contained in:
Kai Bao committed 2025-09-12 09:38:04 +02:00
1 parent 1d75bf47fe
commit 7d5663cd47
22 files changed
+160 -170
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

opm/simulators/wells/BlackoilWellModel.hpp
+2 -2
View File
@@ -511,8 +511,8 @@ template<class Scalar> class WellContributions;
void calculateProductivityIndexValues(const WellInterface<TypeTag>* wellPtr,
DeferredLogger& deferred_logger);
// The number of components in the model.
int numComponents() const;
// The number of conservation quantities.
int numConservationQuantities() const;
int reportStepIndex() const;
opm/simulators/wells/BlackoilWellModel_impl.hpp
+4 -14
View File
@@ -1051,7 +1051,7 @@ namespace Opm {
this->param_,
*this->rateConverter_,
global_pvtreg,
this->numComponents(),
this->numConservationQuantities(),
this->numPhases(),
wellID,
perf_data);
@@ -2170,7 +2170,7 @@ namespace Opm {
BlackoilWellModel<TypeTag>::
updateAverageFormationFactor()
{
std::vector< Scalar > B_avg(numComponents(), Scalar() );
std::vector< Scalar > B_avg(numConservationQuantities(), Scalar() );
const auto& grid = simulator_.vanguard().grid();
const auto& gridView = grid.leafGridView();
ElementContext elemCtx(simulator_);
@@ -2241,19 +2241,9 @@ namespace Opm {
// The number of components in the model.
template<typename TypeTag>
int
BlackoilWellModel<TypeTag>::numComponents() const
BlackoilWellModel<TypeTag>::numConservationQuantities() const
{
// The numComponents here does not reflect the actual number of the components in the system.
// It more or less reflects the number of mass conservation equations for the well equations.
// For example, in the current formulation, we do not have the polymer conservation equation
// in the well equations. As a result, for an oil-water-polymer system, this function will return 2.
// In some way, it makes this function appear to be confusing from its name, and we need
// to revisit/revise this function again when extending the variants of system that flow can simulate.
int numComp = this->numPhases() < 3 ? this->numPhases() : FluidSystem::numComponents;
if constexpr (has_solvent_) {
numComp++;
}
return numComp;
return this->numPhases() + has_solvent_;
}
template<typename TypeTag>
opm/simulators/wells/MultisegmentWell.hpp
+1 -1
View File
@@ -78,7 +78,7 @@ namespace Opm {
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data);
opm/simulators/wells/MultisegmentWellEval.cpp
+5 -5
View File
@@ -90,7 +90,7 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
const bool relax_tolerance,
const bool well_is_stopped) const
{
assert(int(B_avg.size()) == baseif_.numComponents());
assert(int(B_avg.size()) == baseif_.numConservationQuantities());
// checking if any residual is NaN or too large. The two large one is only handled for the well flux
std::vector<std::vector<Scalar>> abs_residual(this->numberOfSegments(),
@@ -107,7 +107,7 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
// 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
if (eq_idx < baseif_.numConservationQuantities()) { // phase or component mass equations
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;
@@ -127,7 +127,7 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& 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
if (eq_idx < baseif_.numConservationQuantities()) { // phase or component mass equations
const Scalar flux_residual = maximum_residual[eq_idx];
// TODO: the report can not handle the segment number yet.
@@ -430,13 +430,13 @@ MultisegmentWellEval<FluidSystem,Indices>::
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const
{
assert(int(B_avg.size() ) == baseif_.numComponents());
assert(int(B_avg.size() ) == baseif_.numConservationQuantities());
std::vector<Scalar> residuals(numWellEq + 1, 0.0);
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
for (int eq_idx = 0; eq_idx < numWellEq; ++eq_idx) {
Scalar residual = 0.;
if (eq_idx < baseif_.numComponents()) {
if (eq_idx < baseif_.numConservationQuantities()) {
residual = std::abs(linSys_.residual()[seg][eq_idx]) * B_avg[eq_idx];
} else {
if (seg > 0) {
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
+1 -1
View File
@@ -591,7 +591,7 @@ surfaceVolumeFraction(const int seg,
const int comp_idx) const
{
EvalWell sum_volume_fraction_scaled = 0.;
for (int idx = 0; idx < well_.numComponents(); ++idx) {
for (int idx = 0; idx < well_.numConservationQuantities(); ++idx) {
sum_volume_fraction_scaled += this->volumeFractionScaled(seg, idx);
}
opm/simulators/wells/MultisegmentWellSegments.cpp
+16 -16
View File
@@ -72,9 +72,9 @@ MultisegmentWellSegments(const int numSegments,
, mass_rates_(numSegments, 0.0)
, viscosities_(numSegments, 0.0)
, upwinding_segments_(numSegments, 0)
, phase_densities_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, phase_fractions_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, phase_viscosities_(numSegments, std::vector<EvalWell>(well.numComponents(), 0.0)) // number of phase here?
, phase_densities_(numSegments, std::vector<EvalWell>(well.numConservationQuantities(), 0.0)) // number of phase here?
, phase_fractions_(numSegments, std::vector<EvalWell>(well.numConservationQuantities(), 0.0)) // number of phase here?
, phase_viscosities_(numSegments, std::vector<EvalWell>(well.numConservationQuantities(), 0.0)) // number of phase here?
, well_(well)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
@@ -143,7 +143,7 @@ computeFluidProperties(const EvalWell& temperature,
int pvt_region_index,
DeferredLogger& deferred_logger)
{
std::vector<Scalar> surf_dens(well_.numComponents());
std::vector<Scalar> surf_dens(well_.numConservationQuantities());
// Surface density.
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@@ -156,13 +156,13 @@ computeFluidProperties(const EvalWell& temperature,
for (std::size_t seg = 0; seg < perforations_.size(); ++seg) {
// the compostion of the components inside wellbore under surface condition
std::vector<EvalWell> mix_s(well_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
std::vector<EvalWell> mix_s(well_.numConservationQuantities(), 0.0);
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
mix_s[comp_idx] = primary_variables.surfaceVolumeFraction(seg, comp_idx);
}
std::vector<EvalWell> b(well_.numComponents(), 0.0);
std::vector<EvalWell> visc(well_.numComponents(), 0.0);
std::vector<EvalWell> b(well_.numConservationQuantities(), 0.0);
std::vector<EvalWell> visc(well_.numConservationQuantities(), 0.0);
std::vector<EvalWell>& phase_densities = phase_densities_[seg];
const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg);
@@ -281,13 +281,13 @@ computeFluidProperties(const EvalWell& temperature,
}
EvalWell volrat(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
volrat += mix[comp_idx] / b[comp_idx];
}
viscosities_[seg] = 0.;
// calculate the average viscosity
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
const EvalWell fraction = mix[comp_idx] / b[comp_idx] / volrat;
// TODO: a little more work needs to be done to handle the negative fractions here
phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0;
@@ -295,14 +295,14 @@ computeFluidProperties(const EvalWell& temperature,
}
EvalWell density(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
density += surf_dens[comp_idx] * mix_s[comp_idx];
}
densities_[seg] = density / volrat;
// calculate the mass rates
mass_rates_[seg] = 0.;
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
const int upwind_seg = upwinding_segments_[seg];
const EvalWell rate = primary_variables.getSegmentRateUpwinding(seg,
upwind_seg,
@@ -366,12 +366,12 @@ getSurfaceVolume(const EvalWell& temperature,
{
const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg_idx);
std::vector<EvalWell> mix_s(well_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
std::vector<EvalWell> mix_s(well_.numConservationQuantities(), 0.0);
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
mix_s[comp_idx] = primary_variables.surfaceVolumeFraction(seg_idx, comp_idx);
}
std::vector<EvalWell> b(well_.numComponents(), 0.);
std::vector<EvalWell> b(well_.numConservationQuantities(), 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
EvalWell rsw(0.0);
@@ -492,7 +492,7 @@ getSurfaceVolume(const EvalWell& temperature,
}
EvalWell vol_ratio(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < well_.numConservationQuantities(); ++comp_idx) {
vol_ratio += mix[comp_idx] / b[comp_idx];
}
opm/simulators/wells/MultisegmentWell_impl.hpp
+31 -31
View File
@@ -65,14 +65,14 @@ namespace Opm
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_conservation_quantities, num_phases, index_of_well, perf_data)
, MSWEval(static_cast<WellInterfaceIndices<FluidSystem,Indices>&>(*this), pw_info)
, regularize_(false)
, segment_fluid_initial_(this->numberOfSegments(), std::vector<Scalar>(this->num_components_, 0.0))
, segment_fluid_initial_(this->numberOfSegments(), std::vector<Scalar>(this->num_conservation_quantities_, 0.0))
{
// not handling solvent or polymer for now with multisegment well
if constexpr (has_solvent) {
@@ -368,13 +368,13 @@ namespace Opm
const int cell_idx = this->well_cells_[local_perf_index];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
// flux for each perforation
std::vector<Scalar> mob(this->num_components_, 0.);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
getMobility(simulator, local_perf_index, mob, deferred_logger);
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(local_perf_index, intQuants, trans_mult, wellstate_nupcol);
const Scalar seg_pressure = segment_pressure[seg];
std::vector<Scalar> cq_s(this->num_components_, 0.);
std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
Scalar perf_press = 0.0;
PerforationRates<Scalar> perf_rates;
computePerfRate(intQuants, mob, Tw, seg, perf, seg_pressure,
@@ -455,7 +455,7 @@ namespace Opm
// compute the potential and store in the flux vector.
well_flux.clear();
well_flux.resize(np, 0.0);
for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
for (int compIdx = 0; compIdx < this->num_conservation_quantities_; ++compIdx) {
const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
well_flux[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
}
@@ -577,7 +577,7 @@ namespace Opm
// fetch potentials (sign is updated on the outside).
well_potentials.clear();
well_potentials.resize(np, 0.0);
for (int compIdx = 0; compIdx < this->num_components_; ++compIdx) {
for (int compIdx = 0; compIdx < this->num_conservation_quantities_; ++compIdx) {
const EvalWell rate = well_copy.primary_variables_.getQs(compIdx);
well_potentials[this->modelCompIdxToFlowCompIdx(compIdx)] = rate.value();
}
@@ -679,7 +679,7 @@ namespace Opm
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// TODO: trying to reduce the times for the surfaceVolumeFraction calculation
const Scalar surface_volume = getSegmentSurfaceVolume(simulator, seg).value();
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
segment_fluid_initial_[seg][comp_idx] = surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx).value();
}
}
@@ -788,7 +788,7 @@ namespace Opm
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<Scalar> mob(this->num_components_, 0.0);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
// The subsetPerfID loops over 0 .. this->perf_data_->size().
// *(this->perf_data_) contains info about the local processes only,
// hence subsetPerfID is a local perf id and we can call getMobility
@@ -933,7 +933,7 @@ namespace Opm
}
// compute component volumetric rates at standard conditions
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->numConservationQuantities(); ++comp_idx) {
const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
}
@@ -954,7 +954,7 @@ namespace Opm
// for injecting perforations, we use total mobility
Value total_mob = mob_perfcells[0];
for (int comp_idx = 1; comp_idx < this->numComponents(); ++comp_idx) {
for (int comp_idx = 1; comp_idx < this->numConservationQuantities(); ++comp_idx) {
total_mob += mob_perfcells[comp_idx];
}
@@ -998,7 +998,7 @@ namespace Opm
}
}
// injecting connections total volumerates at standard conditions
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown);
Value cqt_is = cqt_i / volume_ratio;
cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
@@ -1070,7 +1070,7 @@ namespace Opm
const Value rv = obtain(fs.Rv());
// not using number_of_phases_ because of solvent
std::vector<Value> b_perfcells(this->num_components_, 0.0);
std::vector<Value> b_perfcells(this->num_conservation_quantities_, 0.0);
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
@@ -1081,8 +1081,8 @@ namespace Opm
b_perfcells[compIdx] = obtain(fs.invB(phaseIdx));
}
std::vector<Value> cmix_s(this->numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
std::vector<Value> cmix_s(this->numConservationQuantities(), 0.0);
for (int comp_idx = 0; comp_idx < this->numConservationQuantities(); ++comp_idx) {
cmix_s[comp_idx] = obtainN(this->primary_variables_.surfaceVolumeFraction(seg, comp_idx));
}
@@ -1272,7 +1272,7 @@ namespace Opm
const int local_perf_index = this->pw_info_.activePerfToLocalPerf(perf);
if (local_perf_index < 0) // then the perforation is not on this process
continue;
std::vector<Scalar> mob(this->num_components_, 0.0);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
// TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobility(simulator, local_perf_index, mob, deferred_logger);
@@ -1287,7 +1287,7 @@ namespace Opm
const Scalar pressure_cell = this->getPerfCellPressure(fs).value();
// calculating the b for the connection
std::vector<Scalar> b_perf(this->num_components_);
std::vector<Scalar> b_perf(this->num_conservation_quantities_);
for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
@@ -1312,7 +1312,7 @@ namespace Opm
const std::vector<Scalar> tw_perf = this->wellIndex(local_perf_index, int_quantities, trans_mult, wellstate_nupcol);
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) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++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;
@@ -1374,7 +1374,7 @@ namespace Opm
/*
// could revert to standard approach here:
updateIPR(simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
@@ -1404,7 +1404,7 @@ namespace Opm
const BVectorWell x_well = this->linSys_.solve(rhs);
constexpr int num_eq = MSWEval::numWellEq;
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
const EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
for (size_t pvIdx = 0; pvIdx < num_eq; ++pvIdx) {
@@ -1819,12 +1819,12 @@ namespace Opm
continue;
const int cell_idx = this->well_cells_[local_perf_index];
const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<EvalWell> mob(this->num_components_, 0.0);
std::vector<EvalWell> mob(this->num_conservation_quantities_, 0.0);
getMobility(simulator, local_perf_index, mob, deferred_logger);
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(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
std::vector<EvalWell> cq_s(this->num_components_, 0.0);
std::vector<EvalWell> cq_s(this->num_conservation_quantities_, 0.0);
EvalWell perf_press;
PerforationRates<Scalar> perfRates;
computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
@@ -1839,7 +1839,7 @@ namespace Opm
}
// store the perf pressure and rates
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
perf_rates[local_perf_index*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
}
perf_press_state[local_perf_index] = perf_press.value();
@@ -1849,7 +1849,7 @@ namespace Opm
perf_data.wat_mass_rates[local_perf_index] = cq_s[watCompIdx].value() * rhow;
}
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
// the cq_s entering mass balance equations need to consider the efficiency factors.
const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_;
@@ -1881,7 +1881,7 @@ namespace Opm
// difficult cases
const Scalar regularization_factor = this->regularize_? this->param_.regularization_factor_wells_ : 1.0;
// for each component
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
const EvalWell accumulation_term = regularization_factor * (segment_surface_volume * this->primary_variables_.surfaceVolumeFraction(seg, comp_idx)
- segment_fluid_initial_[seg][comp_idx]) / dt;
MultisegmentWellAssemble(*this).
@@ -1891,7 +1891,7 @@ namespace Opm
// considering the contributions due to flowing out from the segment
{
const int seg_upwind = this->segments_.upwinding_segment(seg);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
const EvalWell segment_rate =
this->primary_variables_.getSegmentRateUpwinding(seg,
seg_upwind,
@@ -1906,7 +1906,7 @@ namespace Opm
{
for (const int inlet : this->segments_.inlets()[seg]) {
const int inlet_upwind = this->segments_.upwinding_segment(inlet);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
const EvalWell inlet_rate =
this->primary_variables_.getSegmentRateUpwinding(inlet,
inlet_upwind,
@@ -2205,7 +2205,7 @@ namespace Opm
DeferredLogger& deferred_logger) const
{
// Calculate the rates that follow from the current primary variables.
std::vector<Scalar> well_q_s(this->num_components_, 0.0);
std::vector<Scalar> well_q_s(this->num_conservation_quantities_, 0.0);
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
@@ -2218,17 +2218,17 @@ namespace Opm
const int cell_idx = this->well_cells_[local_perf_index];
const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<Scalar> mob(this->num_components_, 0.0);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
getMobility(simulator, local_perf_index, mob, deferred_logger);
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(local_perf_index, int_quants, trans_mult, wellstate_nupcol);
std::vector<Scalar> cq_s(this->num_components_, 0.0);
std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.0);
Scalar perf_press = 0.0;
PerforationRates<Scalar> perf_rates;
computePerfRate(int_quants, mob, Tw, seg, perf, seg_pressure,
allow_cf, cq_s, perf_press, perf_rates, deferred_logger);
for (int comp = 0; comp < this->num_components_; ++comp) {
for (int comp = 0; comp < this->num_conservation_quantities_; ++comp) {
well_q_s[comp] += cq_s[comp];
}
}
opm/simulators/wells/ParallelWellInfo.cpp
+13 -13
View File
@@ -172,15 +172,15 @@ int GlobalPerfContainerFactory<Scalar>::globalToLocal(const int globalIndex) con
template<class Scalar>
std::vector<Scalar> GlobalPerfContainerFactory<Scalar>::
createGlobal(const std::vector<Scalar>& local_perf_container,
std::size_t num_components) const
std::size_t num_quantities) const
{
// Could be become templated later.
using Value = Scalar;
if (comm_.size() > 1)
{
std::vector<Value> global_recv(perf_ecl_index_.size() * num_components);
if (num_components == 1)
std::vector<Value> global_recv(perf_ecl_index_.size() * num_quantities);
if (num_quantities == 1)
{
comm_.allgatherv(local_perf_container.data(), local_perf_container.size(),
global_recv.data(), const_cast<int*>(sizes_.data()),
@@ -189,12 +189,12 @@ createGlobal(const std::vector<Scalar>& local_perf_container,
else
{
#if HAVE_MPI
// Create MPI type for sending num_components entries
// Create MPI type for sending num_quantities entries
MPI_Datatype data_type;
MPI_Type_contiguous(num_components, Dune::MPITraits<Value>::getType(), &data_type);
MPI_Type_contiguous(num_quantities, Dune::MPITraits<Value>::getType(), &data_type);
MPI_Type_commit(&data_type);
MPI_Allgatherv(local_perf_container.data(),
local_perf_container.size()/num_components,
local_perf_container.size()/num_quantities,
data_type, global_recv.data(), sizes_.data(),
displ_.data(), data_type, comm_);
MPI_Type_free(&data_type);
@@ -202,13 +202,13 @@ createGlobal(const std::vector<Scalar>& local_perf_container,
}
// reorder by ascending ecl index.
std::vector<Value> global_remapped(perf_ecl_index_.size() * num_components);
std::vector<Value> global_remapped(perf_ecl_index_.size() * num_quantities);
auto global = global_remapped.begin();
for (auto map_entry = map_received_.begin(); map_entry != map_received_.end(); ++map_entry)
{
auto global_index = *map_entry * num_components;
auto global_index = *map_entry * num_quantities;
for(std::size_t i = 0; i < num_components; ++i)
for(std::size_t i = 0; i < num_quantities; ++i)
*(global++) = global_recv[global_index++];
}
assert(global == global_remapped.end());
@@ -224,7 +224,7 @@ template<class Scalar>
void GlobalPerfContainerFactory<Scalar>::
copyGlobalToLocal(const std::vector<Scalar>& global,
std::vector<Scalar>& local,
std::size_t num_components) const
std::size_t num_quantities) const
{
if (global.empty())
{
@@ -241,9 +241,9 @@ copyGlobalToLocal(const std::vector<Scalar>& global,
global_perf = std::lower_bound(global_perf, perf_ecl_index_.end(), pair.global());
assert(global_perf != perf_ecl_index_.end());
assert(*global_perf == pair.global());
auto local_index = pair.local() * num_components;
auto global_index = (global_perf - perf_ecl_index_.begin()) * num_components;
for (std::size_t i = 0; i < num_components; ++i)
auto local_index = pair.local() * num_quantities;
auto global_index = (global_perf - perf_ecl_index_.begin()) * num_quantities;
for (std::size_t i = 0; i < num_quantities; ++i)
local[local_index++] = global[global_index++];
}
}
opm/simulators/wells/ParallelWellInfo.hpp
+4 -4
View File
@@ -151,18 +151,18 @@ public:
/// \brief Creates a container that holds values for all perforations
/// \param local_perf_container Container with values attached to the local perforations.
/// \param num_components the number of components per perforation.
/// \param num_quantities the number of quantities per perforation.
/// \return A container with values attached to all perforations of a well.
/// Values are ordered by the index of the perforation in the ECL schedule.
std::vector<Scalar> createGlobal(const std::vector<Scalar>& local_perf_container,
std::size_t num_components) const;
std::size_t num_quantities) const;
/// \brief Copies the values of the global perforation to the local representation
/// \param global values attached to all peforations of a well (as if the well would live on one process)
/// \param num_components the number of components per perforation.
/// \param num_quantities the number of quantities per perforation.
/// \param[out] local The values attached to the local perforations only.
void copyGlobalToLocal(const std::vector<Scalar>& global, std::vector<Scalar>& local,
std::size_t num_components) const;
std::size_t num_quantities) const;
int numGlobalPerfs() const;
int globalToLocal(const int globalIndex) const;
opm/simulators/wells/StandardWell.hpp
+1 -1
View File
@@ -128,7 +128,7 @@ namespace Opm
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data);
opm/simulators/wells/StandardWellConnections.cpp
+31 -31
View File
@@ -233,7 +233,7 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
: static_cast<Ix>(-1);
const int nperf = this->well_.numLocalPerfs();
const int num_comp = this->well_.numComponents();
const int num_quantities = this->well_.numConservationQuantities();
// 1. Compute the flow (in surface volume units for each component)
// exiting up the wellbore from each perforation, taking into account
@@ -245,11 +245,11 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
// volume ratios (formation factors) for each perforation. Finally
// compute densities for the segments associated with each
// perforation.
auto componentMixture = std::vector<Scalar>(num_comp, Scalar{0});
auto phaseMixture = std::vector<Scalar>(num_comp, Scalar{0});
auto componentMixture = std::vector<Scalar>(num_quantities, Scalar{0});
auto phaseMixture = std::vector<Scalar>(num_quantities, Scalar{0});
for (int perf = 0; perf < nperf; ++perf) {
this->initialiseConnectionMixture(num_comp, perf,
this->initialiseConnectionMixture(num_quantities, perf,
q_out_perf,
phaseMixture,
componentMixture);
@@ -275,13 +275,13 @@ computeDensities(const std::vector<Scalar>& perfComponentRates,
// Compute connection level mixture density as a weighted average of
// phase densities.
const auto* const rho_s = &props.surf_dens_perf[perf*num_comp + 0];
const auto* const b = &props.b_perf [perf*num_comp + 0];
const auto* const rho_s = &props.surf_dens_perf[perf*num_quantities + 0];
const auto* const b = &props.b_perf [perf*num_quantities + 0];
auto& rho = this->perf_densities_[perf];
auto volrat = rho = Scalar{0};
for (auto comp = 0*num_comp; comp < num_comp; ++comp) {
for (auto comp = 0*num_quantities; comp < num_quantities; ++comp) {
rho += componentMixture[comp] * rho_s[comp];
volrat += phaseMixture [comp] / b [comp];
}
@@ -296,9 +296,9 @@ StandardWellConnections<FluidSystem, Indices>::
calculatePerforationOutflow(const std::vector<Scalar>& perfComponentRates) const
{
const int nperf = this->well_.numLocalPerfs();
const int num_comp = this->well_.numComponents();
const int num_quantities = this->well_.numConservationQuantities();
auto q_out_perf = std::vector<Scalar>(nperf * num_comp, Scalar{0});
auto q_out_perf = std::vector<Scalar>(nperf * num_quantities, Scalar{0});
// Component flow rates depend on the order of the perforations. Thus,
// we must use the global view of the well's perforation to get an
@@ -307,17 +307,17 @@ calculatePerforationOutflow(const std::vector<Scalar>& perfComponentRates) const
const auto& factory = this->well_.parallelWellInfo()
.getGlobalPerfContainerFactory();
auto global_q_out_perf = factory.createGlobal(q_out_perf, num_comp);
auto global_q_out_perf = factory.createGlobal(q_out_perf, num_quantities);
const auto global_perf_comp_rates = factory
.createGlobal(perfComponentRates, num_comp);
.createGlobal(perfComponentRates, num_quantities);
// TODO: Investigate whether we should use the following techniques to
// calcuate the composition of flows in the wellbore. Iterate over well
// perforations from bottom to top.
for (int perf = factory.numGlobalPerfs() - 1; perf >= 0; --perf) {
for (int component = 0; component < num_comp; ++component) {
const auto index = perf*num_comp + component;
for (int component = 0; component < num_quantities; ++component) {
const auto index = perf * num_quantities + component;
auto& q_out = global_q_out_perf[index];
// Initialise current perforation's component flow rate to that
@@ -325,7 +325,7 @@ calculatePerforationOutflow(const std::vector<Scalar>& perfComponentRates) const
// component rate.
q_out = (perf == factory.numGlobalPerfs() - 1)
? Scalar{0}
: global_q_out_perf[index + num_comp];
: global_q_out_perf[index + num_quantities];
// Subtract outflow through perforation.
q_out -= global_perf_comp_rates[index];
@@ -333,14 +333,14 @@ calculatePerforationOutflow(const std::vector<Scalar>& perfComponentRates) const
}
// Copy the data back to local view.
factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_comp);
factory.copyGlobalToLocal(global_q_out_perf, q_out_perf, num_quantities);
return q_out_perf;
}
template<typename FluidSystem, typename Indices>
void StandardWellConnections<FluidSystem, Indices>::
initialiseConnectionMixture(const int num_comp,
initialiseConnectionMixture(const int num_quantities,
const int perf,
const std::vector<Scalar>& q_out_perf,
const std::vector<Scalar>& phaseMixture,
@@ -348,18 +348,18 @@ initialiseConnectionMixture(const int num_comp,
{
// Find component mix.
const auto tot_surf_rate =
std::accumulate(q_out_perf.begin() + num_comp*(perf + 0),
q_out_perf.begin() + num_comp*(perf + 1), Scalar{0});
std::accumulate(q_out_perf.begin() + num_quantities*(perf + 0),
q_out_perf.begin() + num_quantities*(perf + 1), Scalar{0});
if (tot_surf_rate != Scalar{0}) {
const auto* const qo = &q_out_perf[perf*num_comp + 0];
const auto* const qo = &q_out_perf[perf*num_quantities + 0];
for (int component = 0; component < num_comp; ++component) {
for (int component = 0; component < num_quantities; ++component) {
componentMixture[component] = std::abs(qo[component] / tot_surf_rate);
}
}
else if (num_comp == 1) {
componentMixture[num_comp - 1] = Scalar{1};
else if (num_quantities == 1) {
componentMixture[num_quantities - 1] = Scalar{1};
}
else {
std::fill(componentMixture.begin(), componentMixture.end(), Scalar{0});
@@ -468,8 +468,8 @@ computePropertiesForPressures(const WellState<Scalar, IndexTraits>& well
const int nperf = well_.numLocalPerfs();
props.b_perf .resize(nperf * this->well_.numComponents());
props.surf_dens_perf.resize(nperf * this->well_.numComponents());
props.b_perf .resize(nperf * this->well_.numConservationQuantities());
props.surf_dens_perf.resize(nperf * this->well_.numConservationQuantities());
const auto& ws = well_state.well(this->well_.indexOfWell());
@@ -521,13 +521,13 @@ computePropertiesForPressures(const WellState<Scalar, IndexTraits>& well
}
}
props.b_perf[waterCompIdx + perf * well_.numComponents()] = FluidSystem::waterPvt()
props.b_perf[waterCompIdx + perf * well_.numConservationQuantities()] = FluidSystem::waterPvt()
.inverseFormationVolumeFactor(region_idx, temperature, p_avg, rsw, saltConcentration);
}
if (gasPresent) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const int gaspos = gasCompIdx + perf * well_.numComponents();
const int gaspos = gasCompIdx + perf * well_.numConservationQuantities();
Scalar rvw = 0.0;
Scalar rv = 0.0;
@@ -574,7 +574,7 @@ computePropertiesForPressures(const WellState<Scalar, IndexTraits>& well
if (oilPresent) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const int oilpos = oilCompIdx + perf * well_.numComponents();
const int oilpos = oilCompIdx + perf * well_.numConservationQuantities();
Scalar rs = 0.0;
if (gasPresent) {
@@ -606,16 +606,16 @@ computePropertiesForPressures(const WellState<Scalar, IndexTraits>& well
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
props.surf_dens_perf[well_.numComponents() * perf + compIdx] =
props.surf_dens_perf[well_.numConservationQuantities() * perf + compIdx] =
FluidSystem::referenceDensity( phaseIdx, region_idx );
}
// We use cell values for solvent injector
if constexpr (Indices::enableSolvent) {
props.b_perf[well_.numComponents() * perf + Indices::contiSolventEqIdx] =
props.b_perf[well_.numConservationQuantities() * perf + Indices::contiSolventEqIdx] =
prop_func.solventInverseFormationVolumeFactor(cell_idx);
props.surf_dens_perf[well_.numComponents() * perf + Indices::contiSolventEqIdx] =
props.surf_dens_perf[well_.numConservationQuantities() * perf + Indices::contiSolventEqIdx] =
prop_func.solventRefDensity(cell_idx);
}
}
@@ -633,7 +633,7 @@ copyInPerforationRates(const Properties& props,
const int nperf = this->well_.numLocalPerfs();
const int np = this->well_.numPhases();
const int nc = this->well_.numComponents();
const int nc = this->well_.numConservationQuantities();
const auto srcIx = [this, np]() {
auto ix = std::vector<int>(np);
opm/simulators/wells/StandardWellEval.cpp
+2 -2
View File
@@ -111,10 +111,10 @@ getWellConvergence(const WellState<Scalar, IndexTraits>& well_state,
res[eq_idx] = std::abs(this->linSys_.residual()[0][eq_idx]);
}
std::vector<Scalar> well_flux_residual(baseif_.numComponents());
std::vector<Scalar> well_flux_residual(baseif_.numConservationQuantities());
// Finish computation
for (int compIdx = 0; compIdx < baseif_.numComponents(); ++compIdx )
for (int compIdx = 0; compIdx < baseif_.numConservationQuantities(); ++compIdx )
{
well_flux_residual[compIdx] = B_avg[compIdx] * res[compIdx];
}
opm/simulators/wells/StandardWellPrimaryVariables.cpp
+2 -2
View File
@@ -519,7 +519,7 @@ StandardWellPrimaryVariables<FluidSystem,Indices>::
surfaceVolumeFraction(const int compIdx) const
{
EvalWell sum_volume_fraction_scaled(numWellEq_ + Indices::numEq, 0.);
for (int idx = 0; idx < well_.numComponents(); ++idx) {
for (int idx = 0; idx < well_.numConservationQuantities(); ++idx) {
sum_volume_fraction_scaled += this->volumeFractionScaled(idx);
}
@@ -534,7 +534,7 @@ StandardWellPrimaryVariables<FluidSystem,Indices>::
getQs(const int comp_idx) const
{
// Note: currently, the WQTotal definition is still depends on Injector/Producer.
assert(comp_idx < well_.numComponents());
assert(comp_idx < well_.numConservationQuantities());
if (well_.isInjector()) { // only single phase injection
Scalar inj_frac = 0.0;
opm/simulators/wells/StandardWell_impl.hpp
+28 -28
View File
@@ -55,15 +55,15 @@ namespace Opm
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data)
: Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_conservation_quantities, num_phases, index_of_well, perf_data)
, StdWellEval(static_cast<const WellInterfaceIndices<FluidSystem,Indices>&>(*this))
, regularize_(false)
{
assert(this->num_components_ == numWellConservationEq);
assert(this->num_conservation_quantities_ == numWellConservationEq);
}
@@ -134,7 +134,7 @@ namespace Opm
const Value rvw = obtain(fs.Rvw());
const Value rsw = obtain(fs.Rsw());
std::vector<Value> b_perfcells_dense(this->numComponents(), zeroElem());
std::vector<Value> b_perfcells_dense(this->numConservationQuantities(), zeroElem());
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
@@ -163,8 +163,8 @@ namespace Opm
}
// surface volume fraction of fluids within wellbore
std::vector<Value> cmix_s(this->numComponents(), zeroElem());
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
std::vector<Value> cmix_s(this->numConservationQuantities(), zeroElem());
for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
cmix_s[componentIdx] = obtainN(this->primary_variables_.surfaceVolumeFraction(componentIdx));
}
@@ -237,7 +237,7 @@ namespace Opm
}
// compute component volumetric rates at standard conditions
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
}
@@ -261,7 +261,7 @@ namespace Opm
// Using total mobilities
Value total_mob_dense = mob[0];
for (int componentIdx = 1; componentIdx < this->numComponents(); ++componentIdx) {
for (int componentIdx = 1; componentIdx < this->numConservationQuantities(); ++componentIdx) {
total_mob_dense += mob[componentIdx];
}
@@ -307,7 +307,7 @@ namespace Opm
}
// injecting connections total volumerates at standard conditions
for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
for (int componentIdx = 0; componentIdx < this->numConservationQuantities(); ++componentIdx) {
const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
Value cqt_is = cqt_i / volumeRatio;
cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
@@ -388,7 +388,7 @@ namespace Opm
auto& perf_rates = perf_data.phase_rates;
for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
// Calculate perforation quantities.
std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
EvalWell water_flux_s{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
EvalWell cq_s_zfrac_effective{this->primary_variables_.numWellEq() + Indices::numEq, 0.0};
calculateSinglePerf(simulator, perf, well_state, connectionRates,
@@ -401,7 +401,7 @@ namespace Opm
water_flux_s, deferred_logger);
}
}
for (int componentIdx = 0; componentIdx < this->num_components_; ++componentIdx) {
for (int componentIdx = 0; componentIdx < this->num_conservation_quantities_; ++componentIdx) {
// the cq_s entering mass balance equations need to consider the efficiency factors.
const EvalWell cq_s_effective = cq_s[componentIdx] * this->well_efficiency_factor_;
@@ -503,7 +503,7 @@ namespace Opm
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
std::vector<EvalWell> mob(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
getMobility(simulator, perf, mob, deferred_logger);
PerforationRates<Scalar> perf_rates;
@@ -683,7 +683,7 @@ namespace Opm
// as a result, the polymer and water share the same viscosity
if constexpr (!Base::has_polymermw) {
if constexpr (std::is_same_v<Value, Scalar>) {
std::vector<EvalWell> mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
std::vector<EvalWell> mob_eval(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
for (std::size_t i = 0; i < mob.size(); ++i) {
mob_eval[i].setValue(mob[i]);
}
@@ -785,7 +785,7 @@ namespace Opm
std::fill(this->ipr_b_.begin(), this->ipr_b_.end(), 0.);
for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
std::vector<Scalar> mob(this->num_components_, 0.0);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
getMobility(simulator, perf, mob, deferred_logger);
const int cell_idx = this->well_cells_[perf];
@@ -795,7 +795,7 @@ namespace Opm
Scalar p_r = this->getPerfCellPressure(fs).value();
// calculating the b for the connection
std::vector<Scalar> b_perf(this->num_components_);
std::vector<Scalar> b_perf(this->num_conservation_quantities_);
for (std::size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
@@ -831,7 +831,7 @@ namespace Opm
wellstate_nupcol);
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) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++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;
@@ -892,7 +892,7 @@ namespace Opm
/*
// could revert to standard approach here:
updateIPR(simulator, deferred_logger);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
ws.implicit_ipr_a[idx] = this->ipr_a_[comp_idx];
ws.implicit_ipr_b[idx] = this->ipr_b_[comp_idx];
@@ -929,7 +929,7 @@ namespace Opm
x_well[0].resize(nEq);
this->linSys_.solve(rhs, x_well);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx){
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx){
EvalWell comp_rate = this->primary_variables_.getQs(comp_idx);
const int idx = this->modelCompIdxToFlowCompIdx(comp_idx);
for (size_t pvIdx = 0; pvIdx < nEq; ++pvIdx) {
@@ -1175,7 +1175,7 @@ namespace Opm
{
// the following implementation assume that the polymer is always after the w-o-g phases
// For the polymer, energy and foam cases, there is one more mass balance equations of reservoir than wells
assert((int(B_avg.size()) == this->num_components_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
assert((int(B_avg.size()) == this->num_conservation_quantities_) || has_polymer || has_energy || has_foam || has_brine || has_zFraction || has_micp);
Scalar tol_wells = this->param_.tolerance_wells_;
// use stricter tolerance for stopped wells and wells under zero rate target control.
@@ -1252,7 +1252,7 @@ namespace Opm
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<Scalar> mob(this->num_components_, 0.0);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.0);
getMobility(simulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID);
@@ -1458,13 +1458,13 @@ namespace Opm
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
// flux for each perforation
std::vector<Scalar> mob(this->num_components_, 0.);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
getMobility(simulator, perf, mob, deferred_logger);
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);
std::vector<Scalar> cq_s(this->num_components_, 0.);
std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
PerforationRates<Scalar> perf_rates;
computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_rates, deferred_logger);
@@ -1646,7 +1646,7 @@ namespace Opm
// fetch potentials (sign is updated on the outside).
well_potentials.clear();
well_potentials.resize(np, 0.0);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
for (int comp_idx = 0; comp_idx < this->num_conservation_quantities_; ++comp_idx) {
if (has_solvent && comp_idx == Indices::contiSolventEqIdx) continue; // we do not store the solvent in the well_potentials
const EvalWell rate = well_copy.primary_variables_.getQs(comp_idx);
well_potentials[this->modelCompIdxToFlowCompIdx(comp_idx)] = rate.value();
@@ -1851,7 +1851,7 @@ namespace Opm
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
std::vector<EvalWell> cq_s(this->num_conservation_quantities_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
PerforationRates<Scalar> perf_rates;
Scalar trans_mult = simulator.problem().template wellTransMultiplier<Scalar>(int_quant, cell_idx);
const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
@@ -2483,22 +2483,22 @@ namespace Opm
DeferredLogger& deferred_logger) const
{
// Calculate the rates that follow from the current primary variables.
std::vector<Scalar> well_q_s(this->num_components_, 0.);
std::vector<Scalar> well_q_s(this->num_conservation_quantities_, 0.);
const EvalWell& bhp = this->primary_variables_.eval(Bhp);
const bool allow_cf = this->getAllowCrossFlow() || openCrossFlowAvoidSingularity(simulator);
for (int perf = 0; perf < this->number_of_local_perforations_; ++perf) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
std::vector<Scalar> mob(this->num_components_, 0.);
std::vector<Scalar> mob(this->num_conservation_quantities_, 0.);
getMobility(simulator, perf, mob, deferred_logger);
std::vector<Scalar> cq_s(this->num_components_, 0.);
std::vector<Scalar> cq_s(this->num_conservation_quantities_, 0.);
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);
PerforationRates<Scalar> perf_rates;
computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
cq_s, perf_rates, deferred_logger);
for (int comp = 0; comp < this->num_components_; ++comp) {
for (int comp = 0; comp < this->num_conservation_quantities_; ++comp) {
well_q_s[comp] += cq_s[comp];
}
}
opm/simulators/wells/WellInterface.hpp
+1 -1
View File
@@ -140,7 +140,7 @@ public:
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data);
opm/simulators/wells/WellInterfaceFluidSystem.cpp
+2 -2
View File
@@ -51,12 +51,12 @@ WellInterfaceFluidSystem(const Well& well,
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data)
: WellInterfaceGeneric<Scalar, IndexTraits>(well, parallel_well_info, time_step, param,
pvtRegionIdx, num_components, num_phases,
pvtRegionIdx, num_conservation_quantities, num_phases,
index_of_well, FluidSystem::phaseUsage(), perf_data)
, rateConverter_(rate_converter)
{
opm/simulators/wells/WellInterfaceFluidSystem.hpp
+1 -1
View File
@@ -74,7 +74,7 @@ protected:
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data);
opm/simulators/wells/WellInterfaceGeneric.cpp
+4 -4
View File
@@ -67,7 +67,7 @@ WellInterfaceGeneric(const Well& well,
const int time_step,
const ModelParameters& param,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const PhaseUsageInfo<IndexTraits>& phase_usage,
@@ -77,13 +77,13 @@ WellInterfaceGeneric(const Well& well,
, current_step_(time_step)
, param_(param)
, pvtRegionIdx_(pvtRegionIdx)
, num_components_(num_components)
, num_conservation_quantities_(num_conservation_quantities)
, number_of_phases_(num_phases)
, index_of_well_(index_of_well)
, phase_usage_(phase_usage)
, perf_data_(&perf_data)
, ipr_a_(num_components)
, ipr_b_(num_components)
, ipr_a_(num_conservation_quantities)
, ipr_b_(num_conservation_quantities)
{
assert(well.name()==pw_info.name());
assert(std::is_sorted(perf_data.begin(), perf_data.end(),
opm/simulators/wells/WellInterfaceGeneric.hpp
+3 -3
View File
@@ -61,7 +61,7 @@ public:
const int time_step,
const ModelParameters& param,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const PhaseUsageInfo<IndexTraits>& phase_usage,
@@ -128,7 +128,7 @@ public:
const GuideRate* guideRate() const { return guide_rate_; }
int numComponents() const { return num_components_; }
int numConservationQuantities() const { return num_conservation_quantities_; }
int numPhases() const { return number_of_phases_; }
@@ -307,7 +307,7 @@ protected:
// We assume a well to not penetrate more than one pvt region.
const int pvtRegionIdx_;
const int num_components_;
const int num_conservation_quantities_;
// number of phases
int number_of_phases_;
opm/simulators/wells/WellInterfaceIndices.cpp
+2 -2
View File
@@ -41,7 +41,7 @@ WellInterfaceIndices(const Well& well,
const ModelParameters& param,
const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data)
@@ -51,7 +51,7 @@ WellInterfaceIndices(const Well& well,
param,
rate_converter,
pvtRegionIdx,
num_components,
num_conservation_quantities,
num_phases,
index_of_well,
perf_data)
opm/simulators/wells/WellInterfaceIndices.hpp
+1 -1
View File
@@ -63,7 +63,7 @@ protected:
const ModelParameters& param,
const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data);
opm/simulators/wells/WellInterface_impl.hpp
+5 -5
View File
@@ -61,7 +61,7 @@ namespace Opm
const ModelParameters& param,
const RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_conservation_quantities,
const int num_phases,
const int index_of_well,
const std::vector<PerforationData<Scalar>>& perf_data)
@@ -71,7 +71,7 @@ namespace Opm
param,
rate_converter,
pvtRegionIdx,
num_components,
num_conservation_quantities,
num_phases,
index_of_well,
perf_data)
@@ -338,7 +338,7 @@ namespace Opm
Scalar inj_limit = inj_controls.bhp_limit;
const bool has_thp = this->wellHasTHPConstraints(summary_state);
if (has_thp){
std::vector<Scalar> rates(this->num_components_);
std::vector<Scalar> rates(this->num_conservation_quantities_);
if (this->isInjector()){
const Scalar bhp_thp = WellBhpThpCalculator(*this).
calculateBhpFromThp(well_state, rates,
@@ -1721,7 +1721,7 @@ namespace Opm
OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly wellIndex was called with a global instead of a local perforation index!");
}
auto wi = std::vector<Scalar>
(this->num_components_, this->well_index_[perf] * trans_mult);
(this->num_conservation_quantities_, this->well_index_[perf] * trans_mult);
if constexpr (! Indices::gasEnabled) {
return wi;
@@ -1914,7 +1914,7 @@ namespace Opm
OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly getMobility was called with a global instead of a local perforation index!");
}
const int cell_idx = this->well_cells_[local_perf_index];
assert (int(mob.size()) == this->num_components_);
assert (int(mob.size()) == this->num_conservation_quantities_);
const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);
const auto& materialLawManager = simulator.problem().materialLawManager();