PI/II: Switch to Using Values Only

We don't need to do the calculations in terms of EvalWell when we're
going to reduce this to the .value() before calling the PI/II
calculation routine.  We can also get by with a simpler approach to
computing the II by assuming we always inject pure phases and no
cross flow in injectors.

Suggested by: [at]atgeirr
This commit is contained in:
Bård Skaflestad 2020-11-06 21:13:10 +01:00
parent 92589a697b
commit bd79d4b9d5
4 changed files with 80 additions and 179 deletions

View File

@ -22,9 +22,10 @@
#ifndef OPM_MULTISEGMENTWELL_HEADER_INCLUDED
#define OPM_MULTISEGMENTWELL_HEADER_INCLUDED
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
namespace Opm
{
@ -185,14 +186,14 @@ namespace Opm
DeferredLogger& deferred_logger) const override;
void computeConnLevelProdInd(const FluidState& fs,
const std::function<double(const EvalWell&)>& connPICalc,
const std::function<double(const double)>& connPICalc,
const std::vector<EvalWell>& mobility,
double* connPI) const;
void computeConnLevelInjInd(const FluidState& fs,
const std::function<double(const EvalWell&)>& connIICalc,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::vector<EvalWell>& mobility,
const int segIx,
double* connII,
DeferredLogger& deferred_logger) const;
protected:

View File

@ -23,6 +23,8 @@
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
#include <algorithm>
namespace Opm
{
@ -1191,18 +1193,6 @@ namespace Opm
.cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
};
auto isInjecting = [this](const typename StandardWell<TypeTag>::FluidState& fs,
const int seg,
const int perf) -> bool
{
const auto pressure_cell = this->extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
const auto perf_seg_press_diff = this->gravity_ * this->segment_densities_[seg] * this->perforation_segment_depth_diffs_[perf];
const auto cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf];
const auto drawdown = pressure_cell - cell_perf_press_diff - (this->getSegmentPressure(seg) + perf_seg_press_diff);
return drawdown.value() < 0.0;
};
const int np = this->number_of_phases_;
auto setToZero = [np](double* x) -> void
{
@ -1214,14 +1204,13 @@ namespace Opm
std::transform(src, src + np, dest, dest, std::plus<>{});
};
const auto allow_xflow = this->getAllowCrossFlow()
|| this->openCrossFlowAvoidSingularity(ebosSimulator);
auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
setToZero(wellPI);
const auto preferred_phase = this->well_ecl_.getPreferredPhase();
const auto& allConn = this->well_ecl_.getConnections();
const auto nPerf = allConn.size();
auto subsetPerfID = 0*nPerf;
@ -1230,29 +1219,23 @@ namespace Opm
continue;
}
auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<EvalWell> mob(this->num_components_, 0.0);
getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob);
const auto& fs = fluidState(subsetPerfID);
auto connPICalc = [&wellPICalc, allPerfID](const EvalWell& mobility) -> double
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility.value());
};
setToZero(connPI);
const auto segIx = this->segmentNumberToIndex(allConn[allPerfID].segment());
if (isInjecting(fs, segIx, subsetPerfID)) { // Connection injects into reservoir
if (this->isInjector() || allow_xflow) {
this->computeConnLevelInjInd(fs, connPICalc, mob, segIx,
connPI, deferred_logger);
}
if (this->isInjector()) {
this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
mob, connPI, deferred_logger);
}
else { // Production or zero flow rate
if (this->isProducer() || allow_xflow) {
this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
}
this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
}
addVector(connPI, wellPI);
@ -3992,7 +3975,7 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
computeConnLevelProdInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
const std::function<double(const EvalWell&)>& connPICalc,
const std::function<double(const double)>& connPICalc,
const std::vector<EvalWell>& mobility,
double* connPI) const
{
@ -4001,8 +3984,9 @@ namespace Opm
for (int p = 0; p < np; ++p) {
// Note: E100's notion of PI value phase mobility includes
// the reciprocal FVF.
const auto connMob = mobility[ flowPhaseToEbosCompIdx(p) ]
* this->extendEval(fs.invB(p));
const auto connMob =
mobility[ flowPhaseToEbosCompIdx(p) ].value()
* fs.invB(p).value();
connPI[p] = connPICalc(connMob);
}
@ -4029,72 +4013,36 @@ namespace Opm
void
MultisegmentWell<TypeTag>::
computeConnLevelInjInd(const typename MultisegmentWell<TypeTag>::FluidState& fs,
const std::function<double(const EvalWell&)>& connIICalc,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::vector<EvalWell>& mobility,
const int segIx,
double* connII,
DeferredLogger& deferred_logger) const
{
// Assumes single phase injection
const auto& pu = this->phaseUsage();
const int np = this->number_of_phases_;
const auto zero = EvalWell { 0.0 };
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
auto cmix_s = std::vector<EvalWell>(this->num_components_, zero);
for (auto componentIdx = 0*this->num_components_;
componentIdx < this->num_components_; ++componentIdx)
{
cmix_s[componentIdx] = this->surfaceVolumeFraction(segIx, componentIdx);
auto phase_pos = 0;
if (preferred_phase == Phase::GAS) {
phase_pos = pu.phase_pos[Gas];
}
auto volumeRatio = zero;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const auto waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volumeRatio += cmix_s[waterCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Water]));
else if (preferred_phase == Phase::OIL) {
phase_pos = pu.phase_pos[Oil];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
{
const auto oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const auto gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const auto rs = this->extendEval(fs.Rs());
const auto rv = this->extendEval(fs.Rv());
// Incorporate RS/RV factors if both oil and gas active
const auto d = 1.0 - rv*rs;
if (d.value() == 0.0) {
OPM_DEFLOG_THROW(Opm::NumericalIssue,
"Zero d value obtained for well " << this->name()
<< " during connection I.I. calculation"
<< " with rs " << rs << " and rv " << rv,
deferred_logger);
}
const auto tmp_oil = (cmix_s[oilCompIdx] - rv*cmix_s[gasCompIdx]) / d;
volumeRatio += tmp_oil / this->extendEval(fs.invB(pu.phase_pos[Oil]));
const auto tmp_gas = (cmix_s[gasCompIdx] - rs*cmix_s[oilCompIdx]) / d;
volumeRatio += tmp_gas / this->extendEval(fs.invB(pu.phase_pos[Gas]));
else if (preferred_phase == Phase::WATER) {
phase_pos = pu.phase_pos[Water];
}
else {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const auto oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
volumeRatio += cmix_s[oilCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Oil]));
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const auto gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
volumeRatio += cmix_s[gasCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Gas]));
}
OPM_DEFLOG_THROW(Opm::NotImplemented,
"Unsupported Injector Type ("
<< static_cast<int>(preferred_phase)
<< ") for well " << this->name()
<< " during connection I.I. calculation",
deferred_logger);
}
const auto cqt_is = mt / volumeRatio;
for (int p = 0; p < np; ++p) {
connII[p] = connIICalc(cmix_s[ flowPhaseToEbosCompIdx(p) ] * cqt_is);
}
const auto zero = EvalWell { 0.0 };
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
connII[phase_pos] = connIICalc(mt.value() * fs.invB(phase_pos).value());
}
} // namespace Opm

View File

@ -39,6 +39,7 @@
#include <opm/models/blackoil/blackoilbrinemodules.hh>
#include <opm/material/densead/DynamicEvaluation.hpp>
#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleTypes.hpp>
#include <dune/common/dynvector.hh>
@ -311,12 +312,13 @@ namespace Opm
DeferredLogger& deferred_logger) const override;
void computeConnLevelProdInd(const FluidState& fs,
const std::function<double(const EvalWell&)>& connPICalc,
const std::function<double(const double)>& connPICalc,
const std::vector<EvalWell>& mobility,
double* connPI) const;
void computeConnLevelInjInd(const FluidState& fs,
const std::function<double(const EvalWell&)>& connIICalc,
void computeConnLevelInjInd(const typename StandardWell<TypeTag>::FluidState& fs,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::vector<EvalWell>& mobility,
double* connII,
DeferredLogger& deferred_logger) const;

View File

@ -2299,16 +2299,6 @@ namespace Opm
.cachedIntensiveQuantities(cell_idx, /*timeIdx=*/ 0)->fluidState();
};
auto isInjecting = [this](const typename StandardWell<TypeTag>::FluidState& fs,
const int perf) -> bool
{
const auto cell_press = this->extendEval(this->getPerfCellPressure(fs));
const auto conn_press = this->getBhp() + this->perf_pressure_diffs_[perf];
const auto drawdown = cell_press - conn_press;
return drawdown.value() < 0.0;
};
const int np = this->number_of_phases_;
auto setToZero = [np](double* x) -> void
{
@ -2320,14 +2310,13 @@ namespace Opm
std::transform(src, src + np, dest, dest, std::plus<>{});
};
const auto allow_xflow = this->getAllowCrossFlow()
|| this->openCrossFlowAvoidSingularity(ebosSimulator);
auto* wellPI = &well_state.productivityIndex()[this->index_of_well_*np + 0];
auto* connPI = &well_state.connectionProductivityIndex()[this->first_perf_*np + 0];
setToZero(wellPI);
const auto preferred_phase = this->well_ecl_.getPreferredPhase();
const auto& allConn = this->well_ecl_.getConnections();
const auto nPerf = allConn.size();
auto subsetPerfID = 0*nPerf;
@ -2336,27 +2325,23 @@ namespace Opm
continue;
}
auto connPICalc = [&wellPICalc, allPerfID](const double mobility) -> double
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility);
};
std::vector<EvalWell> mob(num_components_, {numWellEq_ + numEq, 0.0});
getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID);
auto connPICalc = [&wellPICalc, allPerfID](const EvalWell& mobility) -> double
{
return wellPICalc.connectionProdIndStandard(allPerfID, mobility.value());
};
setToZero(connPI);
if (isInjecting(fs, subsetPerfID)) { // Connection injects into reservoir
if (this->isInjector() || allow_xflow) {
this->computeConnLevelInjInd(fs, connPICalc, mob,
connPI, deferred_logger);
}
if (this->isInjector()) {
this->computeConnLevelInjInd(fs, preferred_phase, connPICalc,
mob, connPI, deferred_logger);
}
else { // Production or zero flow rate
if (this->isProducer() || allow_xflow) {
this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
}
this->computeConnLevelProdInd(fs, connPICalc, mob, connPI);
}
addVector(connPI, wellPI);
@ -4283,7 +4268,7 @@ namespace Opm
void
StandardWell<TypeTag>::
computeConnLevelProdInd(const typename StandardWell<TypeTag>::FluidState& fs,
const std::function<double(const EvalWell&)>& connPICalc,
const std::function<double(const double)>& connPICalc,
const std::vector<EvalWell>& mobility,
double* connPI) const
{
@ -4292,8 +4277,8 @@ namespace Opm
for (int p = 0; p < np; ++p) {
// Note: E100's notion of PI value phase mobility includes
// the reciprocal FVF.
const auto connMob = mobility[ flowPhaseToEbosCompIdx(p) ]
* this->extendEval(fs.invB(p));
const auto connMob =
mobility[ flowPhaseToEbosCompIdx(p) ].value() * fs.invB(p).value();
connPI[p] = connPICalc(connMob);
}
@ -4320,71 +4305,36 @@ namespace Opm
void
StandardWell<TypeTag>::
computeConnLevelInjInd(const typename StandardWell<TypeTag>::FluidState& fs,
const std::function<double(const EvalWell&)>& connIICalc,
const Phase preferred_phase,
const std::function<double(const double)>& connIICalc,
const std::vector<EvalWell>& mobility,
double* connII,
DeferredLogger& deferred_logger) const
{
// Assumes single phase injection
const auto& pu = this->phaseUsage();
const int np = this->number_of_phases_;
const auto zero = EvalWell { this->numWellEq_ + this->numEq, 0.0 };
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
auto cmix_s = std::vector<EvalWell>(this->num_components_, zero);
for (auto componentIdx = 0*this->num_components_;
componentIdx < this->num_components_; ++componentIdx)
{
cmix_s[componentIdx] = this->wellSurfaceVolumeFraction(componentIdx);
auto phase_pos = 0;
if (preferred_phase == Phase::GAS) {
phase_pos = pu.phase_pos[Gas];
}
auto volumeRatio = zero;
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
const auto waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
volumeRatio += cmix_s[waterCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Water]));
else if (preferred_phase == Phase::OIL) {
phase_pos = pu.phase_pos[Oil];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))
{
const auto oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const auto gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
const auto rs = this->extendEval(fs.Rs());
const auto rv = this->extendEval(fs.Rv());
// Incorporate RS/RV factors if both oil and gas active
const auto d = EvalWell { this->numWellEq_ + this->numEq, 1.0 } - rv*rs;
if (d.value() == 0.0) {
OPM_DEFLOG_THROW(Opm::NumericalIssue,
"Zero d value obtained for well " << this->name()
<< " during connection I.I. calculation"
<< " with rs " << rs << " and rv " << rv,
deferred_logger);
}
const auto tmp_oil = (cmix_s[oilCompIdx] - rv*cmix_s[gasCompIdx]) / d;
volumeRatio += tmp_oil / this->extendEval(fs.invB(pu.phase_pos[Oil]));
const auto tmp_gas = (cmix_s[gasCompIdx] - rs*cmix_s[oilCompIdx]) / d;
volumeRatio += tmp_gas / this->extendEval(fs.invB(pu.phase_pos[Gas]));
else if (preferred_phase == Phase::WATER) {
phase_pos = pu.phase_pos[Water];
}
else {
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
const auto oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
volumeRatio += cmix_s[oilCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Oil]));
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
const auto gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
volumeRatio += cmix_s[gasCompIdx] / this->extendEval(fs.invB(pu.phase_pos[Gas]));
}
OPM_DEFLOG_THROW(Opm::NotImplemented,
"Unsupported Injector Type ("
<< static_cast<int>(preferred_phase)
<< ") for well " << this->name()
<< " during connection I.I. calculation",
deferred_logger);
}
const auto cqt_is = mt / volumeRatio;
for (int p = 0; p < np; ++p) {
connII[p] = connIICalc(cmix_s[ flowPhaseToEbosCompIdx(p) ] * cqt_is);
}
const auto zero = EvalWell { this->numWellEq_ + this->numEq, 0.0 };
const auto mt = std::accumulate(mobility.begin(), mobility.end(), zero);
connII[phase_pos] = connIICalc(mt.value() * fs.invB(phase_pos).value());
}
} // namespace Opm