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 #4744 from bska/smry-phase-densities

Browse Source 
Calculate Segment Phase and Mixture Densities for Summary Output
This commit is contained in:
Bård Skaflestad 2023-07-05 13:05:57 +02:00 committed by GitHub
commit ac9c5fb13c
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 303 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

134
ebos/ecloutputblackoilmodule.hh
View File

@ -731,8 +731,8 @@ public:
const auto& problem = elemCtx.simulator().problem();
for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
// Adding block data
unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
// Adding block data
const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
@ -756,14 +756,16 @@ public:
val.second = getValue(fs.pressure(gasPhaseIdx));
else if (FluidSystem::phaseIsActive(waterPhaseIdx))
val.second = getValue(fs.pressure(waterPhaseIdx));
} else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) {
}
else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) {
if (FluidSystem::phaseIsActive(oilPhaseIdx))
val.second = getValue(fs.temperature(oilPhaseIdx));
else if (FluidSystem::phaseIsActive(gasPhaseIdx))
val.second = getValue(fs.temperature(gasPhaseIdx));
else if (FluidSystem::phaseIsActive(waterPhaseIdx))
val.second = getValue(fs.temperature(waterPhaseIdx));
} else if (key.first == "BWKR" || key.first == "BKRW")
}
else if (key.first == "BWKR" || key.first == "BKRW")
val.second = getValue(intQuants.relativePermeability(waterPhaseIdx));
else if (key.first == "BGKR" || key.first == "BKRG")
val.second = getValue(intQuants.relativePermeability(gasPhaseIdx));
@ -774,12 +776,14 @@ public:
const auto krog
= MaterialLaw::template relpermOilInOilGasSystem<Evaluation>(materialParams, fs);
val.second = getValue(krog);
} else if (key.first == "BKROW") {
}
else if (key.first == "BKROW") {
const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0);
const auto krow
= MaterialLaw::template relpermOilInOilWaterSystem<Evaluation>(materialParams, fs);
val.second = getValue(krow);
} else if (key.first == "BWPC")
}
else if (key.first == "BWPC")
val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
else if (key.first == "BGPC")
val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
@ -793,28 +797,42 @@ public:
val.second = getValue(fs.viscosity(gasPhaseIdx));
else if (key.first == "BVOIL" || key.first == "BOVIS")
val.second = getValue(fs.viscosity(oilPhaseIdx));
else if ((key.first == "BRPV") || (key.first == "BOPV") || (key.first == "BWPV")
|| (key.first == "BGPV")) {
else if ((key.first == "BODEN") || (key.first == "BDENO"))
val.second = getValue(fs.density(oilPhaseIdx));
else if ((key.first == "BGDEN") || (key.first == "BDENG"))
val.second = getValue(fs.density(gasPhaseIdx));
else if ((key.first == "BWDEN") || (key.first == "BDENW"))
val.second = getValue(fs.density(waterPhaseIdx));
else if ((key.first == "BRPV") ||
(key.first == "BOPV") ||
(key.first == "BWPV") ||
(key.first == "BGPV"))
{
if (key.first == "BRPV") {
val.second = 1.0;
} else if (key.first == "BOPV") {
}
else if (key.first == "BOPV") {
val.second = getValue(fs.saturation(oilPhaseIdx));
} else if (key.first == "BWPV") {
}
else if (key.first == "BWPV") {
val.second = getValue(fs.saturation(waterPhaseIdx));
} else {
}
else {
val.second = getValue(fs.saturation(gasPhaseIdx));
}
// Include active pore-volume.
val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
* getValue(intQuants.porosity());
} else if (key.first == "BRS")
}
else if (key.first == "BRS")
val.second = getValue(fs.Rs());
else if (key.first == "BRV")
val.second = getValue(fs.Rv());
else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG")
|| (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG")
|| (key.first == "BWIP")) {
else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG") ||
(key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG") ||
(key.first == "BWIP"))
{
if ((key.first == "BOIP") || (key.first == "BOIPL")) {
val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx));
@ -822,27 +840,32 @@ public:
val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
* getValue(fs.saturation(gasPhaseIdx));
}
} else if (key.first == "BOIPG") {
}
else if (key.first == "BOIPG") {
val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
* getValue(fs.saturation(gasPhaseIdx));
} else if ((key.first == "BGIP") || (key.first == "BGIPG")) {
}
else if ((key.first == "BGIP") || (key.first == "BGIPG")) {
val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx));
if (key.first == "BGIP") {
val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
* getValue(fs.saturation(oilPhaseIdx));
}
} else if (key.first == "BGIPL") {
}
else if (key.first == "BGIPL") {
val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
* getValue(fs.saturation(oilPhaseIdx));
} else { // BWIP
}
else { // BWIP
val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx));
}
// Include active pore-volume.
val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
* getValue(intQuants.porosity());
} else {
}
else {
std::string logstring = "Keyword '";
logstring.append(key.first);
logstring.append("' is unhandled for output to file.");
@ -882,14 +905,21 @@ public:
* to globally unique Cartesian cell/element index.
*/
template <class ActiveIndex, class CartesianIndex>
void processFluxes(const ElementContext& elemCtx, ActiveIndex&& activeIndex, CartesianIndex&& cartesianIndex)
void processFluxes(const ElementContext& elemCtx,
ActiveIndex&& activeIndex,
CartesianIndex&& cartesianIndex)
{
OPM_TIMEBLOCK_LOCAL(processFluxes);
const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem) -> EclInterRegFlowMap::Cell {
const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
-> EclInterRegFlowMap::Cell
{
const auto cellIndex = activeIndex(elem);
return {
static_cast<int>(cellIndex), cartesianIndex(cellIndex), elem.partitionType() == Dune::InteriorEntity};
static_cast<int>(cellIndex),
cartesianIndex(cellIndex),
elem.partitionType() == Dune::InteriorEntity
};
};
const auto timeIdx = 0u;
@ -898,10 +928,11 @@ public:
for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
const auto& face = stencil.interiorFace(scvfIdx);
const auto left = identifyCell(stencil.element(face.interiorIndex()));
const auto left = identifyCell(stencil.element(face.interiorIndex()));
const auto right = identifyCell(stencil.element(face.exteriorIndex()));
const auto rates = this->getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
const auto rates = this->
getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
this->interRegionFlows_.addConnection(left, right, rates);
}
@ -1004,10 +1035,12 @@ public:
updateFluidInPlace_(elemCtx, dofIdx);
}
}
void updateFluidInPlace(unsigned globalDofIdx,const IntensiveQuantities& intQuants, double totVolume)
{
this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
}
private:
bool isDefunctParallelWell(std::string wname) const override
{
@ -1231,10 +1264,11 @@ private:
*
* \return Surface level component flow rates.
*/
data::InterRegFlowMap::FlowRates getComponentSurfaceRates(const ElementContext& elemCtx,
const Scalar faceArea,
const std::size_t scvfIdx,
const std::size_t timeIdx) const
data::InterRegFlowMap::FlowRates
getComponentSurfaceRates(const ElementContext& elemCtx,
const Scalar faceArea,
const std::size_t scvfIdx,
const std::size_t timeIdx) const
{
using Component = data::InterRegFlowMap::Component;
@ -1245,57 +1279,71 @@ private:
const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
const auto& up = elemCtx
.intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
using FluidState = std::remove_cv_t<std::remove_reference_t<
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex();
const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), oilPhaseIdx, pvtReg));
const auto bO = getValue(getInvB_<FluidSystem, FluidState, Scalar>
(up.fluidState(), oilPhaseIdx, pvtReg));
const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
rates[Component::Oil] += qO;
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
const auto Rs = getValue(BlackOil::getRs_<FluidSystem, FluidState, Scalar>(up.fluidState(), pvtReg));
const auto Rs = getValue(
BlackOil::getRs_<FluidSystem, FluidState, Scalar>
(up.fluidState(), pvtReg));
rates[Component::Gas] += qO * Rs;
rates[Component::Gas] += qO * Rs;
rates[Component::Disgas] += qO * Rs;
}
}
if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
const auto& up = elemCtx
.intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
using FluidState = std::remove_cv_t<std::remove_reference_t<
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex();
const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), gasPhaseIdx, pvtReg));
const auto bG = getValue(getInvB_<FluidSystem, FluidState, Scalar>
(up.fluidState(), gasPhaseIdx, pvtReg));
const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
rates[Component::Gas] += qG;
if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
const auto Rv = getValue(BlackOil::getRv_<FluidSystem, FluidState, Scalar>(up.fluidState(), pvtReg));
const auto Rv = getValue(
BlackOil::getRv_<FluidSystem, FluidState, Scalar>
(up.fluidState(), pvtReg));
rates[Component::Oil] += qG * Rv;
rates[Component::Oil] += qG * Rv;
rates[Component::Vapoil] += qG * Rv;
}
}
if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
const auto& up = elemCtx.intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
const auto& up = elemCtx
.intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(up.fluidState())>>;
using FluidState = std::remove_cv_t<std::remove_reference_t<
decltype(up.fluidState())>>;
const auto pvtReg = up.pvtRegionIndex();
const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), waterPhaseIdx, pvtReg));
const auto bW = getValue(getInvB_<FluidSystem, FluidState, Scalar>
(up.fluidState(), waterPhaseIdx, pvtReg));
rates[Component::Water] += alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
rates[Component::Water] +=
alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
}
return rates;

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

@ -23,6 +23,7 @@
#include <opm/common/ErrorMacros.hpp>
#include <opm/input/eclipse/Schedule/MSW/AICD.hpp>
#include <opm/input/eclipse/Schedule/MSW/WellSegments.hpp>
#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
@ -37,10 +38,23 @@
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/wells/MSWellHelpers.hpp>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <fmt/format.h>
#include <algorithm>
#include <array>
#include <cmath>
#include <cstddef>
#include <functional>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
namespace Opm
{
@ -317,17 +331,15 @@ updateUpwindingSegments(const PrimaryVariables& primary_variables)
}
}
template<class FluidSystem, class Indices, class Scalar>
typename MultisegmentWellSegments<FluidSystem,Indices,Scalar>::EvalWell
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getHydroPressureLoss(const int seg,
const int seg_density) const
{
const int seg_density) const
{
return densities_[seg_density] * well_.gravity() * depth_diffs_[seg];
}
template<class FluidSystem, class Indices, class Scalar>
Scalar MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
getPressureDiffSegPerf(const int seg,
@ -779,6 +791,133 @@ accelerationPressureLoss(const int seg) const
return accelerationPressureLoss;
}
template <class FluidSystem, class Indices, class Scalar>
void
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
copyPhaseDensities(const PhaseUsage& pu, SegmentState& segSol) const
{
auto* rho = segSol.phase_density.data();
const auto phaseMap = std::vector {
std::pair { BlackoilPhases::Liquid, FluidSystem::oilPhaseIdx },
std::pair { BlackoilPhases::Vapour, FluidSystem::gasPhaseIdx },
std::pair { BlackoilPhases::Aqua , FluidSystem::waterPhaseIdx },
};
// Densities stored in 'rho' as
// [{ p0, p1, ..., (np - 1), mixture, mixture_with_exponents },
// { p0, p1, ..., (np - 1), mixture, mixture_with_exponents },
// ...
// { p0, p1, ..., (np - 1), mixture, mixture_with_exponents }]
// Stride is np + 2.
for (const auto& [boPhase, fsPhaseIdx] : phaseMap) {
if (pu.phase_used[boPhase]) {
this->copyPhaseDensities(fsPhaseIdx, pu.num_phases + 2,
rho + pu.phase_pos[boPhase]);
}
}
// Mixture densities.
for (auto seg = 0*this->densities_.size(); seg < this->densities_.size(); ++seg) {
const auto mixOffset = seg*(pu.num_phases + 2) + pu.num_phases;
rho[mixOffset + 0] = this->mixtureDensity(seg);
rho[mixOffset + 1] = this->mixtureDensityWithExponents(seg);
}
}
template <class FluidSystem, class Indices, class Scalar>
void
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
copyPhaseDensities(const unsigned phaseIdx,
const std::size_t stride,
double* dens) const
{
const auto compIdx = Indices::canonicalToActiveComponentIndex
(FluidSystem::solventComponentIndex(phaseIdx));
for (const auto& phase_density : this->phase_densities_) {
*dens = phase_density[compIdx].value();
dens += stride;
}
}
template <class FluidSystem, class Indices, class Scalar>
double
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
mixtureDensity(const int seg) const
{
auto mixDens = 0.0;
const auto& rho = this->phase_densities_[seg];
const auto& q = this->phase_fractions_[seg];
for (const auto& phIdx : {
FluidSystem::oilPhaseIdx,
FluidSystem::gasPhaseIdx,
FluidSystem::waterPhaseIdx
})
{
if (! FluidSystem::phaseIsActive(phIdx)) {
continue;
}
const auto compIdx = Indices::
canonicalToActiveComponentIndex(phIdx);
mixDens += q[compIdx].value() * rho[compIdx].value();
}
return mixDens;
}
template <class FluidSystem, class Indices, class Scalar>
double
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
mixtureDensityWithExponents(const int seg) const
{
if (const auto& segment = this->well_.wellEcl().getSegments()[seg];
segment.isAICD())
{
return this->mixtureDensityWithExponents(segment.autoICD(), seg);
}
// No other segment type includes exponents of flowing fractions when
// calculating the mixture/emulsion density.
return this->mixtureDensity(seg);
}
template <class FluidSystem, class Indices, class Scalar>
double
MultisegmentWellSegments<FluidSystem,Indices,Scalar>::
mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const
{
auto mixDens = 0.0;
const auto& rho = this->phase_densities_[seg];
const auto& q = this->phase_fractions_[seg];
constexpr auto densityExponents = std::array {
std::pair { FluidSystem::oilPhaseIdx , &AutoICD::oilDensityExponent },
std::pair { FluidSystem::gasPhaseIdx , &AutoICD::gasDensityExponent },
std::pair { FluidSystem::waterPhaseIdx, &AutoICD::waterDensityExponent },
};
for (const auto& [fsPhaseIdx, densityExponent] : densityExponents) {
if (FluidSystem::phaseIsActive(fsPhaseIdx)) {
const auto compIdx = Indices::
canonicalToActiveComponentIndex(fsPhaseIdx);
// exp = (aicd.*densityExponent)() in native syntax.
const auto exp = std::invoke(densityExponent, aicd);
mixDens += std::pow(q[compIdx].value(), exp) * rho[compIdx].value();
}
}
return mixDens;
}
#define INSTANCE(...) \
template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

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

@ -24,13 +24,20 @@
#include <opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp>
#include <cstddef>
#include <vector>
namespace Opm
{
namespace Opm {
class UnitSystem;
class WellInterfaceGeneric;
class AutoICD;
struct PhaseUsage;
class SegmentState;
class UnitSystem;
class WellInterfaceGeneric;
} // namespace Opm
namespace Opm {
template<typename FluidSystem, typename Indices, typename Scalar>
class MultisegmentWellSegments
@ -52,7 +59,7 @@ public:
void updateUpwindingSegments(const PrimaryVariables& primary_variables);
EvalWell getHydroPressureLoss(const int seg,
const int seg_side) const;
const int seg_side) const;
//! Pressure difference between segment and perforation.
Scalar getPressureDiffSegPerf(const int seg,
@ -109,6 +116,9 @@ public:
return perforation_depth_diffs_[perf];
}
void copyPhaseDensities(const PhaseUsage& pu,
SegmentState& segSol) const;
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.
@ -149,8 +159,16 @@ private:
std::vector<std::vector<EvalWell>> phase_viscosities_;
const WellInterfaceGeneric& well_;
void copyPhaseDensities(const unsigned phaseIdx,
const std::size_t stride,
double* dens) const;
double mixtureDensity(const int seg) const;
double mixtureDensityWithExponents(const int seg) const;
double mixtureDensityWithExponents(const AutoICD& aicd, const int seg) const;
};
}
} // namespace Opm
#endif // OPM_MULTISEGMENTWELL_SEGMENTS_HEADER_INCLUDED

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

@ -600,6 +600,12 @@ namespace Opm
this->primary_variables_.copyToWellState(*this, getRefDensity(), stop_or_zero_rate_target,
well_state, summary_state, deferred_logger);
{
auto& ws = well_state.well(this->index_of_well_);
this->segments_.copyPhaseDensities(ws.pu, ws.segments);
}
Base::calculateReservoirRates(well_state.well(this->index_of_well_));
}

15
opm/simulators/wells/SegmentState.cpp
View File

@ -61,6 +61,7 @@ SegmentState::SegmentState(int num_phases, const WellSegments& segments)
, phase_velocity (segments.size() * num_phases)
, phase_holdup (segments.size() * num_phases)
, phase_viscosity (segments.size() * num_phases)
, phase_density (segments.size() * (num_phases + 2)) // +2 for mixture with and without exponents
, pressure (segments.size())
, pressure_drop_friction (segments.size())
, pressure_drop_hydrostatic(segments.size())
@ -77,12 +78,13 @@ SegmentState SegmentState::serializationTestObject()
result.phase_resv_rates = {7.0, 8.0};
result.phase_velocity = {9.0};
result.phase_holdup = {10.0, 11.0};
result.phase_viscosity = {12.0};
result.pressure = {13.0, 14.0};
result.pressure_drop_friction = {15.0};
result.pressure_drop_hydrostatic = {16.0, 17.0};
result.pressure_drop_accel = {18.0};
result.m_segment_number = {19, 20};
result.phase_viscosity = {12.0, 12.5};
result.phase_density = {13.0, 13.5, 13.6, 13.75};
result.pressure = {14.0, 15.0};
result.pressure_drop_friction = {16.0};
result.pressure_drop_hydrostatic = {17.0, 18.0};
result.pressure_drop_accel = {19.0};
result.m_segment_number = {20, 21};
return result;
}
@ -124,6 +126,7 @@ bool SegmentState::operator==(const SegmentState& rhs) const
this->phase_velocity == rhs.phase_velocity &&
this->phase_holdup == rhs.phase_holdup &&
this->phase_viscosity == rhs.phase_viscosity &&
this->phase_density == rhs.phase_density &&
this->pressure == rhs.pressure &&
this->pressure_drop_friction == rhs.pressure_drop_friction &&
this->pressure_drop_hydrostatic == rhs.pressure_drop_hydrostatic &&

19
opm/simulators/wells/SegmentState.hpp
View File

@ -56,6 +56,7 @@ public:
serializer(phase_velocity);
serializer(phase_holdup);
serializer(phase_viscosity);
serializer(phase_density);
serializer(pressure);
serializer(pressure_drop_friction);
serializer(pressure_drop_hydrostatic);
@ -81,6 +82,24 @@ public:
/// Segment condition phase viscosities.
std::vector<double> phase_viscosity;
/// Segment condition phase densities.
///
/// \code num_phases + 2 \endcode elements per segment, with the two
/// additional elements being fluid mixture densities. Fluid mixture
/// densities are calculated with and without ICD-related flowing phase
/// fraction exponents, respectively. Mixture density *without* flowing
/// phase fraction exponents is stored at offset \c num_phases, and
/// mixture density *with* flowing phase fraction exponents is stored at
/// offset \code num_phases + 1 \endcode.
///
/// In other words, the phase and mixture densities are stored as
///
/// [{ p0, p1, ..., (np - 1), mixture, mixture_with_exponents },
/// { p0, p1, ..., (np - 1), mixture, mixture_with_exponents },
/// ...
/// { p0, p1, ..., (np - 1), mixture, mixture_with_exponents }]
std::vector<double> phase_density;
std::vector<double> pressure;
std::vector<double> pressure_drop_friction;
std::vector<double> pressure_drop_hydrostatic;

11
opm/simulators/wells/WellState.cpp
View File

@ -887,6 +887,7 @@ WellState::reportSegmentResults(const int well_id,
const int seg_no) const
{
using PhaseQuant = data::SegmentPhaseQuantity::Item;
using PhaseDensity = data::SegmentPhaseDensity::Item;
const auto& segments = this->well(well_id).segments;
if (segments.empty()) {
@ -911,6 +912,7 @@ WellState::reportSegmentResults(const int well_id,
const auto* velocity = &segments.phase_velocity[seg_ix * pu.num_phases];
const auto* holdup = &segments.phase_holdup[seg_ix * pu.num_phases];
const auto* viscosity = &segments.phase_viscosity[seg_ix * pu.num_phases];
const auto* density = &segments.phase_density[seg_ix * (pu.num_phases + 2)]; // +2 for mixture densities
if (pu.phase_used[Water]) {
const auto iw = pu.phase_pos[Water];
@ -920,6 +922,7 @@ WellState::reportSegmentResults(const int well_id,
seg_res.velocity.set(PhaseQuant::Water, velocity[iw]);
seg_res.holdup.set(PhaseQuant::Water, holdup[iw]);
seg_res.viscosity.set(PhaseQuant::Water, viscosity[iw]);
seg_res.density.set(PhaseDensity::Water, density[iw]);
}
if (pu.phase_used[Oil]) {
@ -931,6 +934,7 @@ WellState::reportSegmentResults(const int well_id,
seg_res.velocity.set(PhaseQuant::Oil, velocity[io]);
seg_res.holdup.set(PhaseQuant::Oil, holdup[io]);
seg_res.viscosity.set(PhaseQuant::Oil, viscosity[io]);
seg_res.density.set(PhaseDensity::Oil, density[io]);
}
if (pu.phase_used[Gas]) {
@ -942,10 +946,17 @@ WellState::reportSegmentResults(const int well_id,
seg_res.velocity.set(PhaseQuant::Gas, velocity[ig]);
seg_res.holdup.set(PhaseQuant::Gas, holdup[ig]);
seg_res.viscosity.set(PhaseQuant::Gas, viscosity[ig]);
seg_res.density.set(PhaseDensity::Gas, density[ig]);
}
seg_res.segNumber = seg_no;
// Recall: The mixture density *without* exponents is stored at offset
// 'num_phases' and the mixture density *with* exponents is stored at
// offset 'num_phases + 1'.
seg_res.density.set(PhaseDensity::Mixture, density[pu.num_phases]);
seg_res.density.set(PhaseDensity::MixtureWithExponents, density[pu.num_phases + 1]);
return seg_res;
}