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 #4334 from akva2/msw_segments

Browse Source 
added: MultisegmentWellSegments
This commit is contained in:
Bård Skaflestad 2022-12-19 17:48:43 +01:00 committed by GitHub
commit 74811e8780
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
12 changed files with 1067 additions and 865 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -96,6 +96,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/MultisegmentWellEval.cpp
opm/simulators/wells/MultisegmentWellGeneric.cpp
opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp
opm/simulators/wells/MultisegmentWellSegments.cpp
opm/simulators/wells/ParallelWellInfo.cpp
opm/simulators/wells/PerfData.cpp
opm/simulators/wells/SegmentState.cpp
@ -382,6 +383,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/MultisegmentWellEval.hpp
opm/simulators/wells/MultisegmentWellGeneric.hpp
opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp
opm/simulators/wells/MultisegmentWellSegments.hpp
opm/simulators/wells/ParallelWellInfo.hpp
opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp

12
opm/simulators/wells/MultisegmentWellEquations.cpp
View File

@ -50,7 +50,9 @@ template<class Scalar, int numWellEq, int numEq>
void MultisegmentWellEquations<Scalar,numWellEq,numEq>::
init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells)
const std::vector<int>& cells,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& perforations)
{
duneB_.setBuildMode(OffDiagMatWell::row_wise);
duneC_.setBuildMode(OffDiagMatWell::row_wise);
@ -64,7 +66,7 @@ init(const int num_cells,
// NNZ = number_of_segments + 2 * (number_of_inlets / number_of_outlets)
{
int nnz_d = well_.numberOfSegments();
for (const std::vector<int>& inlets : well_.segmentInlets()) {
for (const std::vector<int>& inlets : segment_inlets) {
nnz_d += 2 * inlets.size();
}
duneD_.setSize(well_.numberOfSegments(), well_.numberOfSegments(), nnz_d);
@ -89,7 +91,7 @@ init(const int num_cells,
row.insert(seg);
// insert the item related to its inlets
for (const int& inlet : well_.segmentInlets()[seg]) {
for (const int& inlet : segment_inlets[seg]) {
row.insert(inlet);
}
}
@ -98,7 +100,7 @@ init(const int num_cells,
for (auto row = duneC_.createbegin(),
end = duneC_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}
@ -108,7 +110,7 @@ init(const int num_cells,
for (auto row = duneB_.createbegin(),
end = duneB_.createend(); row != end; ++row) {
// the number of the row corresponds to the segment number now.
for (const int& perf : well_.segmentPerforations()[row.index()]) {
for (const int& perf : perforations[row.index()]) {
const int cell_idx = cells[perf];
row.insert(cell_idx);
}

6
opm/simulators/wells/MultisegmentWellEquations.hpp
View File

@ -71,9 +71,13 @@ public:
//! \param num_cells Total number of cells
//! \param numPerfs Number of perforations
//! \param cells Cell indices for perforations
//! \param segment_inlets Cell indices for segment inlets
//! \param segment_perforations Cell indices for segment perforations
void init(const int num_cells,
const int numPerfs,
const std::vector<int>& cells);
const std::vector<int>& cells,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations);
//! \brief Set all coefficients to 0.
void clear();

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

@ -57,13 +57,7 @@ MultisegmentWellEval(WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif)
, baseif_(baseif)
, linSys_(*this)
, primary_variables_(baseif)
, upwinding_segments_(this->numberOfSegments(), 0)
, segment_densities_(this->numberOfSegments(), 0.0)
, segment_mass_rates_(this->numberOfSegments(), 0.0)
, segment_viscosities_(this->numberOfSegments(), 0.0)
, segment_phase_densities_(this->numberOfSegments(), std::vector<EvalWell>(baseif_.numComponents(), 0.0)) // number of phase here?
, segment_phase_fractions_(this->numberOfSegments(), std::vector<EvalWell>(baseif_.numComponents(), 0.0)) // number of phase here?
, segment_phase_viscosities_(this->numberOfSegments(), std::vector<EvalWell>(baseif_.numComponents(), 0.0)) // number of phase here?
, segments_(this->numberOfSegments(), baseif)
, cell_perforation_depth_diffs_(baseif_.numPerfs(), 0.0)
, cell_perforation_pressure_diffs_(baseif_.numPerfs(), 0.0)
{
@ -74,7 +68,9 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
initMatrixAndVectors(const int num_cells)
{
linSys_.init(num_cells, baseif_.numPerfs(), baseif_.cells());
linSys_.init(num_cells, baseif_.numPerfs(),
baseif_.cells(), segments_.inlets(),
segments_.perforations());
primary_variables_.resize(this->numberOfSegments());
}
@ -183,607 +179,21 @@ extendEval(const Eval& in) const
return out;
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
int pvt_region_index,
DeferredLogger& deferred_logger)
{
std::vector<double> surf_dens(baseif_.numComponents());
// Surface density.
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index);
}
for (int seg = 0; seg < this->numberOfSegments(); ++seg) {
// the compostion of the components inside wellbore under surface condition
std::vector<EvalWell> mix_s(baseif_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
mix_s[comp_idx] = primary_variables_.surfaceVolumeFraction(seg, comp_idx);
}
std::vector<EvalWell> b(baseif_.numComponents(), 0.0);
std::vector<EvalWell> visc(baseif_.numComponents(), 0.0);
std::vector<EvalWell>& phase_densities = segment_phase_densities_[seg];
const EvalWell seg_pressure = primary_variables_.getSegmentPressure(seg);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[waterCompIdx] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
visc[waterCompIdx] =
FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration);
// TODO: double check here
// TODO: should not we use phaseIndex here?
phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx];
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[oilCompIdx] > 0.0) {
if (mix_s[gasCompIdx] > 0.0) {
rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
}
if (rv > rvmax) {
rv = rvmax;
}
b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv, rvw);
visc[gasCompIdx] =
FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv, rvw);
phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx]
+ rv * b[gasCompIdx] * surf_dens[oilCompIdx];
} else { // no oil exists
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx];
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
}
}
EvalWell rs(0.0);
// oil phase
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[gasCompIdx] > 0.0) {
if (mix_s[oilCompIdx] > 0.0) {
rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
}
if (rs > rsmax) {
rs = rsmax;
}
b[oilCompIdx] =
FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
visc[oilCompIdx] =
FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx]
+ rs * b[oilCompIdx] * surf_dens[gasCompIdx];
} else { // no oil exists
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx];
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
}
}
segment_phase_viscosities_[seg] = visc;
std::vector<EvalWell> mix(mix_s);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << baseif_.name()
<< " during segment density calculations with rs " << rs
<< ", rv " << rv << " and pressure " << seg_pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
} else {
if (rs > 0.0) {
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv > 0.0) {
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
}
}
}
EvalWell volrat(0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
volrat += mix[comp_idx] / b[comp_idx];
}
this->segment_viscosities_[seg] = 0.;
// calculate the average viscosity
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++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
this->segment_phase_fractions_[seg][comp_idx] = fraction; // >= 0.0 ? fraction : 0.0;
this->segment_viscosities_[seg] += visc[comp_idx] * this->segment_phase_fractions_[seg][comp_idx];
}
EvalWell density(0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
density += surf_dens[comp_idx] * mix_s[comp_idx];
}
this->segment_densities_[seg] = density / volrat;
// calculate the mass rates
segment_mass_rates_[seg] = 0.;
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
const int upwind_seg = upwinding_segments_[seg];
const EvalWell rate = primary_variables_.getSegmentRateUpwinding(seg,
upwind_seg,
comp_idx);
this->segment_mass_rates_[seg] += rate * surf_dens[comp_idx];
}
}
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getHydroPressureLoss(const int seg) const
{
return segment_densities_[seg] * baseif_.gravity() * this->segment_depth_diffs_[seg];
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getFrictionPressureLoss(const int seg) const
{
const EvalWell mass_rate = segment_mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = segment_densities_[seg_upwind];
EvalWell visc = segment_viscosities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
density.clearDerivatives();
visc.clearDerivatives();
}
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const double length = this->segmentSet()[seg].totalLength() - this->segmentSet()[outlet_segment_index].totalLength();
assert(length > 0.);
const double roughness = this->segmentSet()[seg].roughness();
const double area = this->segmentSet()[seg].crossArea();
const double diameter = this->segmentSet()[seg].internalDiameter();
const double sign = mass_rate < 0. ? 1.0 : - 1.0;
return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
pressureDropSpiralICD(const int seg) const
{
const SICD& sicd = this->segmentSet()[seg].spiralICD();
const int seg_upwind = upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = segment_phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = segment_phase_viscosities_[seg_upwind];
EvalWell water_fraction = 0.;
EvalWell water_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_fraction = phase_fractions[water_pos];
water_viscosity = phase_viscosities[water_pos];
}
EvalWell oil_fraction = 0.;
EvalWell oil_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
oil_fraction = phase_fractions[oil_pos];
oil_viscosity = phase_viscosities[oil_pos];
}
EvalWell gas_fraction = 0.;
EvalWell gas_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
}
EvalWell density = segment_densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
density.clearDerivatives();
}
const EvalWell liquid_fraction = water_fraction + oil_fraction;
// viscosity contribution from the liquid
const EvalWell liquid_viscosity_fraction = liquid_fraction < 1.e-30 ? oil_fraction * oil_viscosity + water_fraction * water_viscosity :
liquid_fraction * mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, oil_fraction, oil_viscosity, sicd);
const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity;
const EvalWell reservoir_rate = segment_mass_rates_[seg] / density;
const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
const double viscosity_cali = sicd.viscosityCalibration();
using MathTool = MathToolbox<EvalWell>;
const double density_cali = sicd.densityCalibration();
const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75);
const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25);
// 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 double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0;
return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const
{
const AutoICD& aicd = this->segmentSet()[seg].autoICD();
const int seg_upwind = this->upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = this->segment_phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = this->segment_phase_viscosities_[seg_upwind];
const std::vector<EvalWell>& phase_densities = this->segment_phase_densities_[seg_upwind];
EvalWell water_fraction = 0.;
EvalWell water_viscosity = 0.;
EvalWell water_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_fraction = phase_fractions[water_pos];
water_viscosity = phase_viscosities[water_pos];
water_density = phase_densities[water_pos];
}
EvalWell oil_fraction = 0.;
EvalWell oil_viscosity = 0.;
EvalWell oil_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
oil_fraction = phase_fractions[oil_pos];
oil_viscosity = phase_viscosities[oil_pos];
oil_density = phase_densities[oil_pos];
}
EvalWell gas_fraction = 0.;
EvalWell gas_viscosity = 0.;
EvalWell gas_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
gas_density = phase_densities[gas_pos];
}
EvalWell density = segment_densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
water_density.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
oil_density.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
gas_density.clearDerivatives();
density.clearDerivatives();
}
using MathTool = MathToolbox<EvalWell>;
const EvalWell mixture_viscosity = MathTool::pow(water_fraction, aicd.waterViscExponent()) * water_viscosity
+ MathTool::pow(oil_fraction, aicd.oilViscExponent()) * oil_viscosity
+ MathTool::pow(gas_fraction, aicd.gasViscExponent()) * gas_viscosity;
const EvalWell mixture_density = MathTool::pow(water_fraction, aicd.waterDensityExponent()) * water_density
+ MathTool::pow(oil_fraction, aicd.oilDensityExponent()) * oil_density
+ MathTool::pow(gas_fraction, aicd.gasDensityExponent()) * gas_density;
const double rho_reference = aicd.densityCalibration();
const double visc_reference = aicd.viscosityCalibration();
const auto volume_rate_icd = this->segment_mass_rates_[seg] * aicd.scalingFactor() / mixture_density;
const double 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.);
// TODO: we did not consider the maximum allowed rate here
const auto result = sign / rho_reference * mixture_density * mixture_density
* MathTool::pow(visc_reference/mixture_viscosity, aicd.viscExponent())
* aicd.strength() * MathTool::pow( -sign * volume_rate_icd, aicd.flowRateExponent())
* std::pow(unit_volume_rate, (2. - aicd.flowRateExponent())) ;
return result;
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
pressureDropValve(const int seg) const
{
const Valve& valve = this->segmentSet()[seg].valve();
const EvalWell& mass_rate = segment_mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell visc = segment_viscosities_[seg_upwind];
EvalWell density = segment_densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
visc.clearDerivatives();
density.clearDerivatives();
}
const double additional_length = valve.pipeAdditionalLength();
const double roughness = valve.pipeRoughness();
const double diameter = valve.pipeDiameter();
const double area = valve.pipeCrossArea();
const EvalWell friction_pressure_loss =
mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc);
const double area_con = valve.conCrossArea();
const double 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;
return sign * (friction_pressure_loss + constriction_pressure_loss);
}
template<typename FluidSystem, typename Indices, typename Scalar>
typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
getSegmentSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const int pvt_region_index,
const int seg_idx) const
{
const EvalWell seg_pressure = primary_variables_.getSegmentPressure(seg_idx);
std::vector<EvalWell> mix_s(baseif_.numComponents(), 0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
mix_s[comp_idx] = primary_variables_.surfaceVolumeFraction(seg_idx, comp_idx);
}
std::vector<EvalWell> b(baseif_.numComponents(), 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[waterCompIdx] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
saltConcentration);
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index,
temperature,
seg_pressure);
if (rvmax < 0.0) { // negative rvmax can happen if the seg_pressure is outside the range of the table
rvmax = 0.0;
}
if (mix_s[oilCompIdx] > 0.0) {
if (mix_s[gasCompIdx] > 0.0) {
rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
}
if (rv > rvmax) {
rv = rvmax;
}
b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
rv,
rvw);
} else { // no oil exists
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
}
EvalWell rs(0.0);
// oil phase
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index,
temperature,
seg_pressure);
if (rsmax < 0.0) { // negative rsmax can happen if the seg_pressure is outside the range of the table
rsmax = 0.0;
}
if (mix_s[gasCompIdx] > 0.0) {
if (mix_s[oilCompIdx] > 0.0) {
rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
}
// std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl;
if (rs > rsmax) {
rs = rsmax;
}
b[oilCompIdx] =
FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
rs);
} else { // no oil exists
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
} else { // no gas phase
// it is the same with zero mix_s[Gas]
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
}
std::vector<EvalWell> mix(mix_s);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv;
if (d <= 0.0 || d > 1.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << baseif_.name()
<< " during conversion to surface volume with rs " << rs
<< ", rv " << rv << " and pressure " << seg_pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
OpmLog::debug(sstr.str());
} else {
if (rs > 0.0) {
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv > 0.0) {
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
}
}
}
EvalWell vol_ratio(0.0);
for (int comp_idx = 0; comp_idx < baseif_.numComponents(); ++comp_idx) {
vol_ratio += mix[comp_idx] / b[comp_idx];
}
// We increase the segment volume with a factor 10 to stabilize the system.
const double volume = this->segmentSet()[seg_idx].volume();
return volume / vol_ratio;
}
template<typename FluidSystem, typename Indices, typename Scalar>
void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
handleAccelerationPressureLoss(const int seg,
WellState& well_state)
{
const double area = this->segmentSet()[seg].crossArea();
const EvalWell mass_rate = segment_mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = segment_densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
density.clearDerivatives();
}
EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
// handling the velocity head of intlet segments
for (const int inlet : this->segment_inlets_[seg]) {
const int seg_upwind_inlet = upwinding_segments_[inlet];
const double inlet_area = this->segmentSet()[inlet].crossArea();
EvalWell inlet_density = this->segment_densities_[seg_upwind_inlet];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (inlet != seg_upwind_inlet) {
inlet_density.clearDerivatives();
}
const EvalWell inlet_mass_rate = segment_mass_rates_[inlet];
accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
}
// We change the sign of the accelerationPressureLoss for injectors.
// Is this correct? Testing indicates that this is what the reference simulator does
const double sign = mass_rate < 0. ? 1.0 : - 1.0;
accelerationPressureLoss *= sign;
const EvalWell accelerationPressureLoss = segments_.accelerationPressureLoss(seg);
auto& segments = well_state.well(baseif_.indexOfWell()).segments;
segments.pressure_drop_accel[seg] = accelerationPressureLoss.value();
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureLoss(seg, seg_upwind, accelerationPressureLoss, linSys_);
assemblePressureLoss(seg,
segments_.upwinding_segment(seg),
accelerationPressureLoss, linSys_);
}
template<typename FluidSystem, typename Indices, typename Scalar>
@ -801,14 +211,14 @@ assembleDefaultPressureEq(const int seg,
// we only consider the hydrostatic pressure loss first
// TODO: we might be able to add member variables to store these values, then we update well state
// after converged
const auto hydro_pressure_drop = getHydroPressureLoss(seg);
const auto hydro_pressure_drop = segments_.getHydroPressureLoss(seg);
auto& ws = well_state.well(baseif_.indexOfWell());
auto& segments = ws.segments;
segments.pressure_drop_hydrostatic[seg] = hydro_pressure_drop.value();
pressure_equation -= hydro_pressure_drop;
if (this->frictionalPressureLossConsidered()) {
const auto friction_pressure_drop = getFrictionPressureLoss(seg);
const auto friction_pressure_drop = segments_.getFrictionPressureLoss(seg);
pressure_equation -= friction_pressure_drop;
segments.pressure_drop_friction[seg] = friction_pressure_drop.value();
}
@ -817,7 +227,7 @@ assembleDefaultPressureEq(const int seg,
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index);
const int seg_upwind = upwinding_segments_[seg];
const int seg_upwind = segments_.upwinding_segment(seg);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEq(seg, seg_upwind, outlet_segment_index,
pressure_equation, outlet_pressure, linSys_);
@ -859,13 +269,13 @@ assembleICDPressureEq(const int seg,
EvalWell icd_pressure_drop;
switch(this->segmentSet()[seg].segmentType()) {
case Segment::SegmentType::SICD :
icd_pressure_drop = pressureDropSpiralICD(seg);
icd_pressure_drop = segments_.pressureDropSpiralICD(seg);
break;
case Segment::SegmentType::AICD :
icd_pressure_drop = pressureDropAutoICD(seg, unit_system);
icd_pressure_drop = segments_.pressureDropAutoICD(seg, unit_system);
break;
case Segment::SegmentType::VALVE :
icd_pressure_drop = pressureDropValve(seg);
icd_pressure_drop = segments_.pressureDropValve(seg);
break;
default: {
OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(this->segmentSet()[seg].segmentNumber())
@ -880,7 +290,7 @@ assembleICDPressureEq(const int seg,
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const EvalWell outlet_pressure = primary_variables_.getSegmentPressure(outlet_segment_index);
const int seg_upwind = upwinding_segments_[seg];
const int seg_upwind = segments_.upwinding_segment(seg);
MultisegmentWellAssemble<FluidSystem,Indices,Scalar>(baseif_).
assemblePressureEq(seg, seg_upwind, outlet_segment_index,
pressure_equation, outlet_pressure,

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

@ -25,6 +25,7 @@
#include <opm/simulators/wells/MultisegmentWellEquations.hpp>
#include <opm/simulators/wells/MultisegmentWellGeneric.hpp>
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <opm/simulators/wells/MultisegmentWellSegments.hpp>
#include <opm/material/densead/Evaluation.hpp>
@ -35,10 +36,6 @@
#include <utility>
#include <vector>
namespace Dune {
template<class Matrix> class UMFPack;
}
namespace Opm
{
@ -59,6 +56,7 @@ protected:
static constexpr int WQTotal = PrimaryVariables::WQTotal;
using Equations = MultisegmentWellEquations<Scalar,numWellEq,Indices::numEq>;
using MSWSegments = MultisegmentWellSegments<FluidSystem,Indices,Scalar>;
using BVector = typename Equations::BVector;
using BVectorWell = typename Equations::BVectorWell;
@ -105,18 +103,6 @@ protected:
const double relaxed_inner_tolerance_pressure_ms_well,
const bool relax_tolerance) const;
void computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
int pvt_region_index,
DeferredLogger& deferred_logger);
EvalWell getFrictionPressureLoss(const int seg) const;
EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getSegmentSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const int pvt_region_index,
const int seg_idx) const;
std::pair<bool, std::vector<Scalar> >
getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
DeferredLogger& deferred_logger) const;
@ -135,41 +121,17 @@ protected:
void handleAccelerationPressureLoss(const int seg,
WellState& well_state);
// pressure drop for Autonomous ICD segment (WSEGAICD)
EvalWell pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
// convert a Eval from reservoir to contain the derivative related to wells
EvalWell extendEval(const Eval& in) const;
const WellInterfaceIndices<FluidSystem,Indices,Scalar>& baseif_;
Equations linSys_; //!< The equation system
PrimaryVariables primary_variables_; //!< The primary variables
// the upwinding segment for each segment based on the flow direction
std::vector<int> upwinding_segments_;
// the densities of segment fluids
// we should not have this member variable
std::vector<EvalWell> segment_densities_;
// the mass rate of the segments
std::vector<EvalWell> segment_mass_rates_;
// the viscosity of the segments
std::vector<EvalWell> segment_viscosities_;
std::vector<std::vector<EvalWell>> segment_phase_densities_;
std::vector<std::vector<EvalWell>> segment_phase_fractions_;
std::vector<std::vector<EvalWell>> segment_phase_viscosities_;
MSWSegments segments_; //!< Segment properties
// depth difference between perforations and the perforated grid cells
std::vector<double> cell_perforation_depth_diffs_;

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

@ -46,69 +46,15 @@ template<typename Scalar>
MultisegmentWellGeneric<Scalar>::
MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
: baseif_(baseif)
, segment_perforations_(numberOfSegments())
, segment_inlets_(numberOfSegments())
, segment_depth_diffs_(numberOfSegments(), 0.0)
, perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
// for other facilities needed but not available from parser, we need to process them here
// initialize the segment_perforations_ and update perforation_segment_depth_diffs_
const WellConnections& completion_set = baseif_.wellEcl().getConnections();
// index of the perforation within wells struct
// there might be some perforations not active, which causes the number of the perforations in
// well_ecl_ and wells struct different
// the current implementation is a temporary solution for now, it should be corrected from the parser
// side
int i_perf_wells = 0;
baseif.perfDepth().resize(baseif_.numPerfs(), 0.);
for (size_t perf = 0; perf < completion_set.size(); ++perf) {
const Connection& connection = completion_set.get(perf);
if (connection.state() == Connection::State::OPEN) {
const int segment_index = segmentNumberToIndex(connection.segment());
if ( segment_index == -1) {
OPM_THROW(std::logic_error,
fmt::format("COMPSEGS: Well {} has connection in cell {}, {}, {} "
"without associated segment.", baseif_.wellEcl().name(),
connection.getI() + 1, connection.getJ() + 1,
connection.getK() + 1));
}
segment_perforations_[segment_index].push_back(i_perf_wells);
baseif.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segmentSet()[segment_index].depth();
perforation_segment_depth_diffs_[i_perf_wells] = baseif.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
}
// initialize the segment_inlets_
for (int seg = 0; seg < numberOfSegments(); ++seg) {
const Segment& segment = segmentSet()[seg];
const int segment_number = segment.segmentNumber();
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) {
const int segment_index = segmentNumberToIndex(segment_number);
const int outlet_segment_index = segmentNumberToIndex(outlet_segment_number);
segment_inlets_[outlet_segment_index].push_back(segment_index);
}
}
// calculating the depth difference between the segment and its oulet_segments
// for the top segment, we will make its zero unless we find other purpose to use this value
for (int seg = 1; seg < numberOfSegments(); ++seg) {
const double segment_depth = segmentSet()[seg].depth();
const int outlet_segment_number = segmentSet()[seg].outletSegment();
const Segment& outlet_segment = segmentSet()[segmentNumberToIndex(outlet_segment_number)];
const double outlet_depth = outlet_segment.depth();
segment_depth_diffs_[seg] = segment_depth - outlet_depth;
}
}
template<typename Scalar>
void
MultisegmentWellGeneric<Scalar>::
scaleSegmentRatesWithWellRates(WellState& well_state) const
scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
WellState& well_state) const
{
auto& ws = well_state.well(baseif_.indexOfWell());
auto& segments = ws.segments;
@ -137,7 +83,10 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
}
std::vector<double> rates;
WellState::calculateSegmentRates(segment_inlets_, segment_perforations_, perforation_rates, num_single_phase, 0, rates);
WellState::calculateSegmentRates(segment_inlets,
segment_perforations,
perforation_rates,
num_single_phase, 0, rates);
for (int seg = 0; seg < numberOfSegments(); ++seg) {
segment_rates[baseif_.numPhases() * seg + phase] = rates[seg];
}
@ -179,22 +128,6 @@ compPressureDrop() const
return segmentSet().compPressureDrop();
}
template<typename Scalar>
const std::vector<std::vector<int>>&
MultisegmentWellGeneric<Scalar>::
segmentInlets() const
{
return segment_inlets_;
}
template<typename Scalar>
const std::vector<std::vector<int>>&
MultisegmentWellGeneric<Scalar>::
segmentPerforations() const
{
return segment_perforations_;
}
template<typename Scalar>
int
MultisegmentWellGeneric<Scalar>::

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

@ -41,12 +41,6 @@ template <typename Scalar>
class MultisegmentWellGeneric
{
public:
//! \brief Returns the inlet segments for each segment.
const std::vector<std::vector<int>>& segmentInlets() const;
//! \brief Returns the perforation index for each segment.
const std::vector<std::vector<int>>& segmentPerforations() const;
// get the WellSegments from the well_ecl_
const WellSegments& segmentSet() const;
@ -61,7 +55,9 @@ protected:
MultisegmentWellGeneric(WellInterfaceGeneric& baseif);
// scale the segment rates and pressure based on well rates and bhp
void scaleSegmentRatesWithWellRates(WellState& well_state) const;
void scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
WellState& well_state) const;
void scaleSegmentPressuresWithBhp(WellState& well_state) const;
// components of the pressure drop to be included
@ -77,27 +73,6 @@ protected:
bool frictionalPressureLossConsidered() const;
const WellInterfaceGeneric& baseif_;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int>> segment_perforations_;
// the inlet segments for each segment. It is for convenience and efficiency reason
std::vector<std::vector<int>> segment_inlets_;
std::vector<double> segment_depth_diffs_;
// depth difference between the segment and the perforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<double> perforation_segment_depth_diffs_;
};
}

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

@ -600,31 +600,6 @@ getWQTotal() const
return evaluation_[0][WQTotal];
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>::
updateUpwindingSegments(const MultisegmentWellGeneric<Scalar>& mswell,
std::vector<int>& upwinding_segments) const
{
for (size_t seg = 0; seg < value_.size(); ++seg) {
// special treatment is needed for segment 0
if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors
assert(!(well_.isProducer() && evaluation_[seg][WQTotal] > 0.));
assert(!(well_.isInjector() && evaluation_[seg][WQTotal] < 0.));
upwinding_segments[seg] = seg;
continue;
}
// for other normal segments
if (evaluation_[seg][WQTotal] <= 0.) {
upwinding_segments[seg] = seg;
} else {
const int outlet_segment_index = mswell.segmentNumberToIndex(mswell.segmentSet()[seg].outletSegment());
upwinding_segments[seg] = outlet_segment_index;
}
}
}
#define INSTANCE(...) \
template class MultisegmentWellPrimaryVariables<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

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

@ -133,10 +133,6 @@ public:
//! \brief Get WQTotal.
EvalWell getWQTotal() const;
//! \brief Update upwinding segments.
void updateUpwindingSegments(const MultisegmentWellGeneric<Scalar>& mswell,
std::vector<int>& upwinding_segments) const;
//! \brief Returns a const ref to an evaluation.
const std::array<EvalWell,numWellEq>& eval(const int idx) const
{ return evaluation_[idx]; }

781
opm/simulators/wells/MultisegmentWellSegments.cpp Normal file
View File

@ -0,0 +1,781 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/MultisegmentWellSegments.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/models/blackoil/blackoilindices.hh>
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <fmt/format.h>
#include <sstream>
namespace Opm
{
template<class FluidSystem, class Indices, class Scalar>
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
MultisegmentWellSegments(const int numSegments,
WellInterfaceGeneric& well)
: perforations_(numSegments)
, perforation_depth_diffs_(well.numPerfs(), 0.0)
, inlets_(well.wellEcl().getSegments().size())
, depth_diffs_(numSegments, 0.0)
, densities_(numSegments, 0.0)
, 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?
, well_(well)
{
// since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
// for other facilities needed but not available from parser, we need to process them here
// initialize the segment_perforations_ and update perforation_segment_depth_diffs_
const WellConnections& completion_set = well_.wellEcl().getConnections();
// index of the perforation within wells struct
// there might be some perforations not active, which causes the number of the perforations in
// well_ecl_ and wells struct different
// the current implementation is a temporary solution for now, it should be corrected from the parser
// side
int i_perf_wells = 0;
well.perfDepth().resize(well_.numPerfs(), 0.);
const auto& segment_set = well_.wellEcl().getSegments();
for (size_t perf = 0; perf < completion_set.size(); ++perf) {
const Connection& connection = completion_set.get(perf);
if (connection.state() == Connection::State::OPEN) {
const int segment_index = segment_set.segmentNumberToIndex(connection.segment());
if (segment_index == -1) {
OPM_THROW(std::logic_error,
fmt::format("COMPSEGS: Well {} has connection in cell {}, {}, {} "
"without associated segment.", well_.wellEcl().name(),
connection.getI() + 1, connection.getJ() + 1,
connection.getK() + 1));
}
perforations_[segment_index].push_back(i_perf_wells);
well.perfDepth()[i_perf_wells] = connection.depth();
const double segment_depth = segment_set[segment_index].depth();
perforation_depth_diffs_[i_perf_wells] = well_.perfDepth()[i_perf_wells] - segment_depth;
i_perf_wells++;
}
}
// initialize the segment_inlets_
for (const Segment& segment : segment_set) {
const int segment_number = segment.segmentNumber();
const int outlet_segment_number = segment.outletSegment();
if (outlet_segment_number > 0) {
const int segment_index = segment_set.segmentNumberToIndex(segment_number);
const int outlet_segment_index = segment_set.segmentNumberToIndex(outlet_segment_number);
inlets_[outlet_segment_index].push_back(segment_index);
}
}
// calculating the depth difference between the segment and its oulet_segments
// for the top segment, we will make its zero unless we find other purpose to use this value
for (int seg = 1; seg < numSegments; ++seg) {
const double 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();
depth_diffs_[seg] = segment_depth - outlet_depth;
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
computeFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
const PrimaryVariables& primary_variables,
int pvt_region_index,
DeferredLogger& deferred_logger)
{
std::vector<double> surf_dens(well_.numComponents());
// Surface density.
for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
if (!FluidSystem::phaseIsActive(phaseIdx)) {
continue;
}
const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
surf_dens[compIdx] = FluidSystem::referenceDensity( phaseIdx, pvt_region_index);
}
for (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) {
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>& phase_densities = phase_densities_[seg];
const EvalWell seg_pressure = primary_variables.getSegmentPressure(seg);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[waterCompIdx] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, saltConcentration);
visc[waterCompIdx] =
FluidSystem::waterPvt().viscosity(pvt_region_index, temperature, seg_pressure, saltConcentration);
// TODO: double check here
// TODO: should not we use phaseIndex here?
phase_densities[waterCompIdx] = b[waterCompIdx] * surf_dens[waterCompIdx];
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[oilCompIdx] > 0.0) {
if (mix_s[gasCompIdx] > 0.0) {
rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
}
if (rv > rvmax) {
rv = rvmax;
}
b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rv, rvw);
visc[gasCompIdx] =
FluidSystem::gasPvt().viscosity(pvt_region_index, temperature, seg_pressure, rv, rvw);
phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx]
+ rv * b[gasCompIdx] * surf_dens[oilCompIdx];
} else { // no oil exists
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
phase_densities[gasCompIdx] = b[gasCompIdx] * surf_dens[gasCompIdx];
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[gasCompIdx] =
FluidSystem::gasPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
}
}
EvalWell rs(0.0);
// oil phase
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index, temperature, seg_pressure);
if (mix_s[gasCompIdx] > 0.0) {
if (mix_s[oilCompIdx] > 0.0) {
rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
}
if (rs > rsmax) {
rs = rsmax;
}
b[oilCompIdx] =
FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure, rs);
visc[oilCompIdx] =
FluidSystem::oilPvt().viscosity(pvt_region_index, temperature, seg_pressure, rs);
phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx]
+ rs * b[oilCompIdx] * surf_dens[gasCompIdx];
} else { // no oil exists
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
phase_densities[oilCompIdx] = b[oilCompIdx] * surf_dens[oilCompIdx];
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index, temperature, seg_pressure);
visc[oilCompIdx] =
FluidSystem::oilPvt().saturatedViscosity(pvt_region_index, temperature, seg_pressure);
}
}
phase_viscosities_[seg] = visc;
std::vector<EvalWell> mix(mix_s);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << well_.name()
<< " during segment density calculations with rs " << rs
<< ", rv " << rv << " and pressure " << seg_pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
} else {
if (rs > 0.0) {
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv > 0.0) {
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
}
}
}
EvalWell volrat(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++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) {
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;
viscosities_[seg] += visc[comp_idx] * phase_fractions_[seg][comp_idx];
}
EvalWell density(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++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) {
const int upwind_seg = upwinding_segments_[seg];
const EvalWell rate = primary_variables.getSegmentRateUpwinding(seg,
upwind_seg,
comp_idx);
mass_rates_[seg] += rate * surf_dens[comp_idx];
}
}
}
template<class FluidSystem, class Indices, class Scalar>
void MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
updateUpwindingSegments(const PrimaryVariables& primary_variables)
{
for (size_t seg = 0; seg < perforations_.size(); ++seg) {
// special treatment is needed for segment 0
if (seg == 0) {
// we are not supposed to have injecting producers and producing injectors
assert(!(well_.isProducer() && primary_variables.eval(seg)[primary_variables.WQTotal] > 0.));
assert(!(well_.isInjector() && primary_variables.eval(seg)[primary_variables.WQTotal] < 0.));
upwinding_segments_[seg] = seg;
continue;
}
// for other normal segments
if (primary_variables.eval(seg)[primary_variables.WQTotal] <= 0.) {
upwinding_segments_[seg] = seg;
} else {
const auto& segment_set = well_.wellEcl().getSegments();
const int outlet_segment_index = segment_set.segmentNumberToIndex(segment_set[seg].outletSegment());
upwinding_segments_[seg] = outlet_segment_index;
}
}
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getHydroPressureLoss(const int seg) const
{
return densities_[seg] * well_.gravity() * depth_diffs_[seg];
}
template<class FluidSystem, class Indices, class Scalar>
Scalar MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getPressureDiffSegPerf(const int seg,
const int perf) const
{
return well_.gravity() * densities_[seg].value() * perforation_depth_diffs_[perf];
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const PrimaryVariables& primary_variables,
const int pvt_region_index,
const int seg_idx) const
{
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) {
mix_s[comp_idx] = primary_variables.surfaceVolumeFraction(seg_idx, comp_idx);
}
std::vector<EvalWell> b(well_.numComponents(), 0.);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
b[waterCompIdx] =
FluidSystem::waterPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
saltConcentration);
}
EvalWell rv(0.0);
EvalWell rvw(0.0);
// gas phase
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
EvalWell rvmax = FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvt_region_index,
temperature,
seg_pressure);
if (rvmax < 0.0) { // negative rvmax can happen if the seg_pressure is outside the range of the table
rvmax = 0.0;
}
if (mix_s[oilCompIdx] > 0.0) {
if (mix_s[gasCompIdx] > 0.0) {
rv = mix_s[oilCompIdx] / mix_s[gasCompIdx];
}
if (rv > rvmax) {
rv = rvmax;
}
b[gasCompIdx] =
FluidSystem::gasPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
rv,
rvw);
} else { // no oil exists
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
} else { // no Liquid phase
// it is the same with zero mix_s[Oil]
b[gasCompIdx] =
FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
}
EvalWell rs(0.0);
// oil phase
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
EvalWell rsmax = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvt_region_index,
temperature,
seg_pressure);
if (rsmax < 0.0) { // negative rsmax can happen if the seg_pressure is outside the range of the table
rsmax = 0.0;
}
if (mix_s[gasCompIdx] > 0.0) {
if (mix_s[oilCompIdx] > 0.0) {
rs = mix_s[gasCompIdx] / mix_s[oilCompIdx];
}
// std::cout << " rs " << rs.value() << " rsmax " << rsmax.value() << std::endl;
if (rs > rsmax) {
rs = rsmax;
}
b[oilCompIdx] =
FluidSystem::oilPvt().inverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure,
rs);
} else { // no oil exists
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
} else { // no gas phase
// it is the same with zero mix_s[Gas]
b[oilCompIdx] =
FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvt_region_index,
temperature,
seg_pressure);
}
}
std::vector<EvalWell> mix(mix_s);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv;
if (d <= 0.0 || d > 1.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << well_.name()
<< " during conversion to surface volume with rs " << rs
<< ", rv " << rv << " and pressure " << seg_pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
OpmLog::debug(sstr.str());
} else {
if (rs > 0.0) {
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv > 0.0) {
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
}
}
}
EvalWell vol_ratio(0.0);
for (int comp_idx = 0; comp_idx < well_.numComponents(); ++comp_idx) {
vol_ratio += mix[comp_idx] / b[comp_idx];
}
// We increase the segment volume with a factor 10 to stabilize the system.
const double volume = well_.wellEcl().getSegments()[seg_idx].volume();
return volume / vol_ratio;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getFrictionPressureLoss(const int seg) const
{
const EvalWell mass_rate = mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = densities_[seg_upwind];
EvalWell visc = viscosities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
density.clearDerivatives();
visc.clearDerivatives();
}
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();
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 double sign = mass_rate < 0. ? 1.0 : - 1.0;
return sign * mswellhelpers::frictionPressureLoss(length, diameter, area, roughness, density, mass_rate, visc);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
pressureDropSpiralICD(const int seg) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const SICD& sicd = segment_set[seg].spiralICD();
const int seg_upwind = upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = phase_viscosities_[seg_upwind];
EvalWell water_fraction = 0.;
EvalWell water_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_fraction = phase_fractions[water_pos];
water_viscosity = phase_viscosities[water_pos];
}
EvalWell oil_fraction = 0.;
EvalWell oil_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
oil_fraction = phase_fractions[oil_pos];
oil_viscosity = phase_viscosities[oil_pos];
}
EvalWell gas_fraction = 0.;
EvalWell gas_viscosity = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
}
EvalWell density = densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
density.clearDerivatives();
}
const EvalWell liquid_fraction = water_fraction + oil_fraction;
// viscosity contribution from the liquid
const EvalWell liquid_viscosity_fraction = liquid_fraction < 1.e-30 ? oil_fraction * oil_viscosity + water_fraction * water_viscosity :
liquid_fraction * mswellhelpers::emulsionViscosity(water_fraction, water_viscosity, oil_fraction, oil_viscosity, sicd);
const EvalWell mixture_viscosity = liquid_viscosity_fraction + gas_fraction * gas_viscosity;
const EvalWell reservoir_rate = mass_rates_[seg] / density;
const EvalWell reservoir_rate_icd = reservoir_rate * sicd.scalingFactor();
const double viscosity_cali = sicd.viscosityCalibration();
using MathTool = MathToolbox<EvalWell>;
const double density_cali = sicd.densityCalibration();
const EvalWell temp_value1 = MathTool::pow(density / density_cali, 0.75);
const EvalWell temp_value2 = MathTool::pow(mixture_viscosity / viscosity_cali, 0.25);
// 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 double sign = reservoir_rate_icd <= 0. ? 1.0 : -1.0;
return sign * temp_value1 * temp_value2 * strength * reservoir_rate_icd * reservoir_rate_icd;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const AutoICD& aicd = segment_set[seg].autoICD();
const int seg_upwind = upwinding_segments_[seg];
const std::vector<EvalWell>& phase_fractions = phase_fractions_[seg_upwind];
const std::vector<EvalWell>& phase_viscosities = phase_viscosities_[seg_upwind];
const std::vector<EvalWell>& phase_densities = phase_densities_[seg_upwind];
EvalWell water_fraction = 0.;
EvalWell water_viscosity = 0.;
EvalWell water_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const int water_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
water_fraction = phase_fractions[water_pos];
water_viscosity = phase_viscosities[water_pos];
water_density = phase_densities[water_pos];
}
EvalWell oil_fraction = 0.;
EvalWell oil_viscosity = 0.;
EvalWell oil_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const int oil_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
oil_fraction = phase_fractions[oil_pos];
oil_viscosity = phase_viscosities[oil_pos];
oil_density = phase_densities[oil_pos];
}
EvalWell gas_fraction = 0.;
EvalWell gas_viscosity = 0.;
EvalWell gas_density = 0.;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const int gas_pos = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
gas_fraction = phase_fractions[gas_pos];
gas_viscosity = phase_viscosities[gas_pos];
gas_density = phase_densities[gas_pos];
}
EvalWell density = densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
water_fraction.clearDerivatives();
water_viscosity.clearDerivatives();
water_density.clearDerivatives();
oil_fraction.clearDerivatives();
oil_viscosity.clearDerivatives();
oil_density.clearDerivatives();
gas_fraction.clearDerivatives();
gas_viscosity.clearDerivatives();
gas_density.clearDerivatives();
density.clearDerivatives();
}
using MathTool = MathToolbox<EvalWell>;
const EvalWell mixture_viscosity = MathTool::pow(water_fraction, aicd.waterViscExponent()) * water_viscosity
+ MathTool::pow(oil_fraction, aicd.oilViscExponent()) * oil_viscosity
+ MathTool::pow(gas_fraction, aicd.gasViscExponent()) * gas_viscosity;
const EvalWell mixture_density = MathTool::pow(water_fraction, aicd.waterDensityExponent()) * water_density
+ MathTool::pow(oil_fraction, aicd.oilDensityExponent()) * oil_density
+ MathTool::pow(gas_fraction, aicd.gasDensityExponent()) * gas_density;
const double rho_reference = aicd.densityCalibration();
const double visc_reference = aicd.viscosityCalibration();
const auto volume_rate_icd = mass_rates_[seg] * aicd.scalingFactor() / mixture_density;
const double 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.);
// TODO: we did not consider the maximum allowed rate here
const auto result = sign / rho_reference * mixture_density * mixture_density
* MathTool::pow(visc_reference/mixture_viscosity, aicd.viscExponent())
* aicd.strength() * MathTool::pow( -sign * volume_rate_icd, aicd.flowRateExponent())
* std::pow(unit_volume_rate, (2. - aicd.flowRateExponent())) ;
return result;
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
pressureDropValve(const int seg) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const Valve& valve = segment_set[seg].valve();
const EvalWell& mass_rate = mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell visc = viscosities_[seg_upwind];
EvalWell density = densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
visc.clearDerivatives();
density.clearDerivatives();
}
const double additional_length = valve.pipeAdditionalLength();
const double roughness = valve.pipeRoughness();
const double diameter = valve.pipeDiameter();
const double area = valve.pipeCrossArea();
const EvalWell friction_pressure_loss =
mswellhelpers::frictionPressureLoss(additional_length, diameter, area, roughness, density, mass_rate, visc);
const double area_con = valve.conCrossArea();
const double 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;
return sign * (friction_pressure_loss + constriction_pressure_loss);
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
accelerationPressureLoss(const int seg) const
{
const auto& segment_set = well_.wellEcl().getSegments();
const double area = segment_set[seg].crossArea();
const EvalWell mass_rate = mass_rates_[seg];
const int seg_upwind = upwinding_segments_[seg];
EvalWell density = densities_[seg_upwind];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (seg != seg_upwind) {
density.clearDerivatives();
}
EvalWell accelerationPressureLoss = mswellhelpers::velocityHead(area, mass_rate, density);
// handling the velocity head of intlet segments
for (const int inlet : inlets_[seg]) {
const int seg_upwind_inlet = upwinding_segments_[inlet];
const double inlet_area = segment_set[inlet].crossArea();
EvalWell inlet_density = densities_[seg_upwind_inlet];
// WARNING
// We disregard the derivatives from the upwind density to make sure derivatives
// wrt. to different segments dont get mixed.
if (inlet != seg_upwind_inlet) {
inlet_density.clearDerivatives();
}
const EvalWell inlet_mass_rate = mass_rates_[inlet];
accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
}
// We change the sign of the accelerationPressureLoss for injectors.
// Is this correct? Testing indicates that this is what the reference simulator does
const double sign = mass_rate < 0. ? 1.0 : - 1.0;
accelerationPressureLoss *= sign;
return accelerationPressureLoss;
}
#define INSTANCE(...) \
template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
// One phase
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
INSTANCE(BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>)
// Two phase
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>)
INSTANCE(BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)
// Blackoil
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,true,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,true,2u,0u>)
INSTANCE(BlackOilIndices<1u,0u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,1u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,1u,0u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,false,0u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,0u,false,false,1u,0u>)
INSTANCE(BlackOilIndices<0u,0u,0u,1u,false,true,0u,0u>)
} // namespace Opm

155
opm/simulators/wells/MultisegmentWellSegments.hpp Normal file
View File

@ -0,0 +1,155 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_MULTISEGMENTWELL_SEGMENTS_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_SEGMENTS_HEADER_INCLUDED
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <vector>
namespace Opm
{
class UnitSystem;
class WellInterfaceGeneric;
template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellSegments
{
using PrimaryVariables = MultisegmentWellPrimaryVariables<FluidSystem,Indices,Scalar>;
using EvalWell = typename PrimaryVariables::EvalWell;
public:
MultisegmentWellSegments(const int numSegments,
WellInterfaceGeneric& well);
void computeFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
const PrimaryVariables& primary_variables,
int pvt_region_index,
DeferredLogger& deferred_logger);
//! \brief Update upwinding segments.
void updateUpwindingSegments(const PrimaryVariables& primary_variables);
EvalWell getHydroPressureLoss(const int seg) const;
//! Pressure difference between segment and perforation.
Scalar getPressureDiffSegPerf(const int seg,
const int perf) const;
EvalWell getSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const PrimaryVariables& primary_variables,
const int pvt_region_index,
const int seg_idx) const;
EvalWell getFrictionPressureLoss(const int seg) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;
// pressure drop for Autonomous ICD segment (WSEGAICD)
EvalWell pressureDropAutoICD(const int seg,
const UnitSystem& unit_system) const;
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
// pressure loss due to acceleration
EvalWell accelerationPressureLoss(const int seg) const;
const std::vector<std::vector<int>>& inlets() const
{
return inlets_;
}
const std::vector<std::vector<int>>& perforations() const
{
return perforations_;
}
int upwinding_segment(const int seg) const
{
return upwinding_segments_[seg];
}
Scalar getRefDensity() const
{
return densities_[0].value();
}
const EvalWell& density(const int seg) const
{
return densities_[seg];
}
Scalar perforation_depth_diff(const int perf) const
{
return perforation_depth_diffs_[perf];
}
private:
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
// the completions that is related to each segment
// the completions's ids are their index in the vector well_index_, well_cell_
// This is also assuming the order of the completions in Well is the same with
// the order of the completions in wells.
// it is for convinience reason. we can just calcuate the inforation for segment once then using it for all the perofrations
// belonging to this segment
std::vector<std::vector<int>> perforations_;
// depth difference between the segment and the perforation
// or in another way, the depth difference between the perforation and
// the segment the perforation belongs to
std::vector<Scalar> perforation_depth_diffs_;
// the inlet segments for each segment. It is for convenience and efficiency reason
std::vector<std::vector<int>> inlets_;
std::vector<Scalar> depth_diffs_;
// the densities of segment fluids
// we should not have this member variable
std::vector<EvalWell> densities_;
// the mass rate of the segments
std::vector<EvalWell> mass_rates_;
// the viscosity of the segments
std::vector<EvalWell> viscosities_;
// the upwinding segment for each segment based on the flow direction
std::vector<int> upwinding_segments_;
std::vector<std::vector<EvalWell>> phase_densities_;
std::vector<std::vector<EvalWell>> phase_fractions_;
std::vector<std::vector<EvalWell>> phase_viscosities_;
const WellInterfaceGeneric& well_;
};
}
#endif // OPM_MULTISEGMENTWELL_SEGMENTS_HEADER_INCLUDED

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

@ -161,7 +161,9 @@ namespace Opm
Base::updateWellStateWithTarget(ebos_simulator, group_state, well_state, deferred_logger);
// scale segment rates based on the wellRates
// and segment pressure based on bhp
this->scaleSegmentRatesWithWellRates(well_state);
this->scaleSegmentRatesWithWellRates(this->segments_.inlets(),
this->segments_.perforations(),
well_state);
this->scaleSegmentPressuresWithBhp(well_state);
}
@ -368,7 +370,7 @@ namespace Opm
segments_copy.scale_pressure(bhp);
const auto& segment_pressure = segments_copy.pressure;
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
// flux for each perforation
@ -441,7 +443,9 @@ namespace Opm
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
well_copy.scaleSegmentRatesWithWellRates(this->segments_.inlets(),
this->segments_.perforations(),
well_state_copy);
well_copy.calculateExplicitQuantities(ebosSimulator, well_state_copy, deferred_logger);
const double dt = ebosSimulator.timeStepSize();
@ -770,7 +774,8 @@ namespace Opm
DeferredLogger& deferred_logger) const
{
// pressure difference between the segment and the perforation
const Value perf_seg_press_diff = this->gravity() * segment_density * this->perforation_segment_depth_diffs_[perf];
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];
@ -931,7 +936,7 @@ namespace Opm
Tw,
perf,
segment_pressure,
this->segment_densities_[seg],
this->segments_.density(seg),
allow_cf,
cmix_s,
cq_s,
@ -992,7 +997,7 @@ namespace Opm
Tw,
perf,
segment_pressure,
getValue(this->segment_densities_[seg]),
getValue(this->segments_.density(seg)),
allow_cf,
cmix_s,
cq_s,
@ -1032,10 +1037,11 @@ namespace Opm
pvt_region_index = fs.pvtRegionIndex();
}
this->MSWEval::computeSegmentFluidProperties(temperature,
saltConcentration,
pvt_region_index,
deferred_logger);
this->segments_.computeFluidProperties(temperature,
saltConcentration,
this->primary_variables_,
pvt_region_index,
deferred_logger);
}
@ -1153,7 +1159,7 @@ namespace Opm
MultisegmentWell<TypeTag>::
getRefDensity() const
{
return this->segment_densities_[0].value();
return this->segments_.getRefDensity();
}
template<typename TypeTag>
@ -1238,83 +1244,84 @@ namespace Opm
const double segment_depth = this->segmentSet()[seg].depth();
const int outlet_segment_index = this->segmentNumberToIndex(this->segmentSet()[seg].outletSegment());
const double segment_depth_outlet = seg == 0? ref_depth : this->segmentSet()[outlet_segment_index].depth();
double dp = wellhelpers::computeHydrostaticCorrection(segment_depth_outlet, segment_depth, this->segment_densities_[seg].value(), this->gravity_);
double dp = wellhelpers::computeHydrostaticCorrection(segment_depth_outlet, segment_depth,
this->segments_.density(seg).value(), this->gravity_);
// we add the hydrostatic correction from the outlet segment
// in order to get the correction all the way to the bhp ref depth.
if (seg > 0) {
dp += seg_dp[outlet_segment_index];
}
seg_dp[seg] = dp;
for (const int perf : this->segment_perforations_[seg]) {
std::vector<Scalar> mob(this->num_components_, 0.0);
for (const int perf : this->segments_.perforations()[seg]) {
std::vector<Scalar> mob(this->num_components_, 0.0);
// TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityScalar(ebos_simulator, perf, mob);
// TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration
getMobilityScalar(ebos_simulator, perf, mob);
const int cell_idx = this->well_cells_[perf];
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(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->gravity_ * this->segment_densities_[seg].value() * this->perforation_segment_depth_diffs_[perf];
// pressure difference between the perforation and the grid cell
const double cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const double pressure_cell = this->getPerfCellPressure(fs).value();
const int cell_idx = this->well_cells_[perf];
const auto& int_quantities = *(ebos_simulator.model().cachedIntensiveQuantities(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);
// 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();
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
// calculating the b for the connection
std::vector<double> b_perf(this->num_components_);
for (size_t phase = 0; phase < FluidSystem::numPhases; ++phase) {
if (!FluidSystem::phaseIsActive(phase)) {
continue;
}
const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
b_perf[comp_idx] = fs.invB(phase).value();
}
const unsigned comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phase));
b_perf[comp_idx] = fs.invB(phase).value();
}
// 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;
// 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;
// do not take into consideration the crossflow here.
if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
deferred_logger.debug("CROSSFLOW_IPR",
"cross flow found when updateIPR for well " + this->name());
}
// do not take into consideration the crossflow here.
if ( (this->isProducer() && pressure_diff < 0.) || (this->isInjector() && pressure_diff > 0.) ) {
deferred_logger.debug("CROSSFLOW_IPR",
"cross flow found when updateIPR for well " + this->name());
}
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
// the well index associated with the connection
const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> 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 * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
std::vector<double> ipr_a_perf(this->ipr_a_.size());
std::vector<double> 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 * mob[comp_idx] * b_perf[comp_idx];
ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
ipr_b_perf[comp_idx] += tw_mob;
}
// we need to handle the rs and rv when both oil and gas are present
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();
// we need to handle the rs and rv when both oil and gas are present
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 double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double vap_oil_a = rv * ipr_a_perf[gas_comp_idx];
const double dis_gas_a = rs * ipr_a_perf[oil_comp_idx];
const double 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;
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 double dis_gas_b = rs * ipr_b_perf[oil_comp_idx];
const double 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;
}
ipr_b_perf[gas_comp_idx] += dis_gas_b;
ipr_b_perf[oil_comp_idx] += vap_oil_b;
}
for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
}
for (size_t comp_idx = 0; comp_idx < ipr_a_perf.size(); ++comp_idx) {
this->ipr_a_[comp_idx] += ipr_a_perf[comp_idx];
this->ipr_b_[comp_idx] += ipr_b_perf[comp_idx];
}
}
}
}
@ -1529,8 +1536,7 @@ namespace Opm
if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
// update the upwinding segments
this->primary_variables_.updateUpwindingSegments(*this,
this->upwinding_segments_);
this->segments_.updateUpwindingSegments(this->primary_variables_);
// calculate the fluid properties needed.
computeSegmentFluidProperties(ebosSimulator, deferred_logger);
@ -1572,7 +1578,7 @@ namespace Opm
}
// considering the contributions due to flowing out from the segment
{
const int seg_upwind = this->upwinding_segments_[seg];
const int seg_upwind = this->segments_.upwinding_segment(seg);
for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
const EvalWell segment_rate =
this->primary_variables_.getSegmentRateUpwinding(seg,
@ -1586,8 +1592,8 @@ namespace Opm
// considering the contributions from the inlet segments
{
for (const int inlet : this->segment_inlets_[seg]) {
const int inlet_upwind = this->upwinding_segments_[inlet];
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) {
const EvalWell inlet_rate =
this->primary_variables_.getSegmentRateUpwinding(inlet,
@ -1605,7 +1611,7 @@ namespace Opm
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
auto& perf_press_state = perf_data.pressure;
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<EvalWell> mob(this->num_components_, 0.0);
@ -1687,14 +1693,14 @@ namespace Opm
for (int seg = 0; seg < nseg; ++seg) {
const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg);
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& intQuants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = intQuants.fluidState();
// pressure difference between the segment and the perforation
const EvalWell perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
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];
@ -1751,10 +1757,11 @@ namespace Opm
pvt_region_index = fs.pvtRegionIndex();
}
return this->MSWEval::getSegmentSurfaceVolume(temperature,
saltConcentration,
pvt_region_index,
seg_idx);
return this->segments_.getSurfaceVolume(temperature,
saltConcentration,
this->primary_variables_,
pvt_region_index,
seg_idx);
}
@ -1888,7 +1895,7 @@ namespace Opm
double max_pressure = 0.0;
const int nseg = this->numberOfSegments();
for (int seg = 0; seg < nseg; ++seg) {
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebos_simulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
const auto& fs = int_quants.fluidState();
@ -1916,7 +1923,7 @@ namespace Opm
for (int seg = 0; seg < nseg; ++seg) {
// calculating the perforation rate for each perforation that belongs to this segment
const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg));
for (const int perf : this->segment_perforations_[seg]) {
for (const int perf : this->segments_.perforations()[seg]) {
const int cell_idx = this->well_cells_[perf];
const auto& int_quants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0));
std::vector<Scalar> mob(this->num_components_, 0.0);