move getSegmentSurfaceVolume to MultisegmentWellSegments

This commit is contained in:
Arne Morten Kvarving 2022-12-19 14:20:55 +01:00
parent 9c19120855
commit 4a9cedf452
5 changed files with 157 additions and 152 deletions

View File

@ -585,150 +585,6 @@ pressureDropValve(const int seg) const
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>::

View File

@ -113,10 +113,6 @@ protected:
DeferredLogger& deferred_logger);
EvalWell getFrictionPressureLoss(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,

View File

@ -33,6 +33,7 @@
#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
#include <fmt/format.h>
#include <sstream>
namespace Opm
{
@ -116,6 +117,151 @@ getHydroPressureLoss(const int seg) const
return densities_[seg] * well_.gravity() * depth_diffs_[seg];
}
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;
}
#define INSTANCE(...) \
template class MultisegmentWellSegments<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;

View File

@ -44,6 +44,12 @@ public:
EvalWell getHydroPressureLoss(const int seg) const;
EvalWell getSurfaceVolume(const EvalWell& temperature,
const EvalWell& saltConcentration,
const PrimaryVariables& primary_variables,
const int pvt_region_index,
const int seg_idx) const;
// 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.

View File

@ -1755,8 +1755,9 @@ namespace Opm
pvt_region_index = fs.pvtRegionIndex();
}
return this->MSWEval::getSegmentSurfaceVolume(temperature,
return this->segments_.getSurfaceVolume(temperature,
saltConcentration,
this->primary_variables_,
pvt_region_index,
seg_idx);
}