Merge pull request #3853 from totto82/fix_d_issues

Fix d issues
This commit is contained in:
Tor Harald Sandve
2022-04-04 12:50:56 +02:00
committed by GitHub
8 changed files with 119 additions and 56 deletions

View File

@@ -228,7 +228,8 @@ namespace Opm
// compute the fluid properties, such as densities, viscosities, and so on, in the segments
// They will be treated implicitly, so they need to be of Evaluation type
void computeSegmentFluidProperties(const Simulator& ebosSimulator);
void computeSegmentFluidProperties(const Simulator& ebosSimulator,
DeferredLogger& deferred_logger);
// get the mobility for specific perforation
void getMobilityEval(const Simulator& ebosSimulator,

View File

@@ -639,7 +639,8 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>::
computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
int pvt_region_index)
int pvt_region_index,
DeferredLogger& deferred_logger)
{
std::vector<double> surf_dens(baseif_.numComponents());
// Surface density.
@@ -757,12 +758,22 @@ computeSegmentFluidProperties(const EvalWell& temperature,
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv;
if (rs != 0.0) { // rs > 0.0?
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv != 0.0) { // rv > 0.0?
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
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;
}
}
}
@@ -1194,16 +1205,17 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
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;
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
OpmLog::debug(sstr.str());
OPM_THROW_NOLOG(NumericalIssue, sstr.str());
}
if (rs > 0.0) { // rs > 0.0?
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
}
if (rv > 0.0) { // rv > 0.0?
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
} 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;
}
}
}

View File

@@ -180,7 +180,8 @@ protected:
void computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration,
int pvt_region_index);
int pvt_region_index,
DeferredLogger& deferred_logger);
EvalWell getBhp() const;
EvalWell getFrictionPressureLoss(const int seg) const;

View File

@@ -1000,7 +1000,7 @@ namespace Opm
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
computeSegmentFluidProperties(const Simulator& ebosSimulator)
computeSegmentFluidProperties(const Simulator& ebosSimulator, DeferredLogger& deferred_logger)
{
// TODO: the concept of phases and components are rather confusing in this function.
// needs to be addressed sooner or later.
@@ -1029,7 +1029,8 @@ namespace Opm
this->MSWEval::computeSegmentFluidProperties(temperature,
saltConcentration,
pvt_region_index);
pvt_region_index,
deferred_logger);
}
@@ -1494,7 +1495,7 @@ namespace Opm
this->updateUpwindingSegments();
// calculate the fluid properties needed.
computeSegmentFluidProperties(ebosSimulator);
computeSegmentFluidProperties(ebosSimulator, deferred_logger);
// clear all entries
this->duneB_ = 0.0;

View File

@@ -264,10 +264,12 @@ namespace Opm
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger);
void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state);
const WellState& well_state,
DeferredLogger& deferred_logger);
void computePerfRateEval(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob,

View File

@@ -831,7 +831,8 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf)
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger)
{
// Verify that we have consistent input.
const int nperf = baseif_.numPerfs();
@@ -950,13 +951,25 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) {
rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
}
if (rs != 0.0) {
// Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv);
}
if (rv != 0.0) {
// Subtract oil in gas from oil mixture
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);;
const double d = 1.0 - rs*rv;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << baseif_.name()
<< " during ccomputeConnectionDensities with rs " << rs
<< ", rv " << rv
<< " 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) {
// Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d;
}
if (rv > 0.0) {
// Subtract oil in gas from oil mixture
x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d;
}
}
}
double volrat = 0.0;

View File

@@ -144,7 +144,8 @@ protected:
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf);
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger);
ConvergenceReport getWellConvergence(const WellState& well_state,
const std::vector<double>& B_avg,

View File

@@ -312,15 +312,20 @@ namespace Opm
// Incorporate RS/RV factors if both oil and gas active
const Value d = 1.0 - rv * rs;
if (getValue(d) == 0.0) {
OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << this->name() << " during flux calcuation"
<< " with rs " << rs << " and rv " << rv, deferred_logger);
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during computePerfRate calculations with rs " << rs
<< ", rv " << rv << " and pressure " << pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
}
const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
const Value tmp_oil = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
volumeRatio += tmp_oil / b_perfcells_dense[oilCompIdx];
const Value tmp_gas = (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d;
const Value tmp_gas = d > 0.0? (cmix_s[gasCompIdx] - rs * cmix_s[oilCompIdx]) / d : cmix_s[gasCompIdx];
volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
}
else {
@@ -354,12 +359,24 @@ namespace Opm
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs);
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during computePerfRate calculations with rs " << rs
<< ", rv " << rv << " and pressure " << pressure
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
} else {
// vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
perf_vap_oil_rate = getValue(rv) * (getValue(cq_s[gasCompIdx]) - getValue(rs) * getValue(cq_s[oilCompIdx])) / d;
// dissolved of gas in oil
// rs * q_or * b_o = rs * (q_os - rv * q_gs) / d
perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d;
}
}
}
}
@@ -580,14 +597,27 @@ namespace Opm
const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
// d = 1.0 - rs * rv
const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
if(FluidSystem::gasPhaseIdx == phaseIdx)
cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
if(FluidSystem::oilPhaseIdx == phaseIdx)
cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
if (d <= 0.0) {
std::ostringstream sstr;
sstr << "Problematic d value " << d << " obtained for well " << this->name()
<< " during calculateSinglePerf with rs " << fs.Rs()
<< ", rv " << fs.Rv()
<< " obtaining d " << d
<< " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection.";
deferred_logger.debug(sstr.str());
cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
} else {
if(FluidSystem::gasPhaseIdx == phaseIdx) {
cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
} else if(FluidSystem::oilPhaseIdx == phaseIdx) {
// q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
}
}
}
// change temperature for injecting fluids
@@ -1442,7 +1472,8 @@ namespace Opm
const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_perf,
const std::vector<double>& surf_dens_perf)
const std::vector<double>& surf_dens_perf,
DeferredLogger& deferred_logger)
{
// Compute densities
const int nperf = this->number_of_perforations_;
@@ -1497,7 +1528,7 @@ namespace Opm
}
}
this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
this->computeConnectionDensities(perfRates, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
this->computeConnectionPressureDelta();
}
@@ -1510,7 +1541,8 @@ namespace Opm
void
StandardWell<TypeTag>::
computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state)
const WellState& well_state,
DeferredLogger& deferred_logger)
{
// 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function
@@ -1520,7 +1552,7 @@ namespace Opm
std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
computeWellConnectionDensitesPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, deferred_logger);
}
@@ -1556,7 +1588,7 @@ namespace Opm
{
updatePrimaryVariables(well_state, deferred_logger);
initPrimaryVariablesEvaluation();
computeWellConnectionPressures(ebosSimulator, well_state);
computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
this->computeAccumWell();
}
@@ -1725,7 +1757,7 @@ namespace Opm
deferred_logger.debug(msg);
}
well.updatePrimaryVariables(well_state_copy, deferred_logger);
well.computeWellConnectionPressures(ebosSimulator, well_state_copy);
well.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
well.initPrimaryVariablesEvaluation();
well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
}