Address comments from review

1) Add debug messages
2) Fix bug of missing else in the code
This commit is contained in:
Tor Harald Sandve
2022-04-01 08:35:31 +02:00
parent 5b53fcd8a6
commit 69ffed06de
8 changed files with 103 additions and 38 deletions

View File

@@ -228,7 +228,8 @@ namespace Opm
// compute the fluid properties, such as densities, viscosities, and so on, in the segments // 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 // 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 // get the mobility for specific perforation
void getMobilityEval(const Simulator& ebosSimulator, void getMobilityEval(const Simulator& ebosSimulator,

View File

@@ -639,7 +639,8 @@ void
MultisegmentWellEval<FluidSystem,Indices,Scalar>:: MultisegmentWellEval<FluidSystem,Indices,Scalar>::
computeSegmentFluidProperties(const EvalWell& temperature, computeSegmentFluidProperties(const EvalWell& temperature,
const EvalWell& saltConcentration, const EvalWell& saltConcentration,
int pvt_region_index) int pvt_region_index,
DeferredLogger& deferred_logger)
{ {
std::vector<double> surf_dens(baseif_.numComponents()); std::vector<double> surf_dens(baseif_.numComponents());
// Surface density. // Surface density.
@@ -757,12 +758,22 @@ computeSegmentFluidProperties(const EvalWell& temperature,
const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
const EvalWell d = 1.0 - rs * rv; const EvalWell d = 1.0 - rs * rv;
if (d <= 0.0) {
if (rs > 0.0 && d > 0.0) { std::ostringstream sstr;
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; sstr << "Problematic d value " << d << " obtained for well " << baseif_.name()
} << " during segment density calculations with rs " << rs
if (rv > 0.0 && d > 0.0) { << ", rv " << rv << " and pressure " << seg_pressure
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; << " 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;
}
} }
} }
@@ -1195,16 +1206,16 @@ getSegmentSurfaceVolume(const EvalWell& temperature,
<< " during conversion to surface volume with rs " << rs << " during conversion to surface volume with rs " << rs
<< ", rv " << rv << " and pressure " << seg_pressure << ", rv " << rv << " and pressure " << seg_pressure
<< " obtaining d " << d << " obtaining d " << d
<< " Continue without dissolution (rs = 0) and vaporization (rv = 0) " << " Continue as if no dissolution (rs = 0) and vaporization (rv = 0) "
<< " for this connection."; << " for this connection.";
OpmLog::debug(sstr.str()); OpmLog::debug(sstr.str());
} } else {
if (rs > 0.0) {
if (rs > 0.0) { // rs > 0.0? mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d;
mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; }
} if (rv > 0.0) {
if (rv > 0.0) { // rv > 0.0? mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d;
mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; }
} }
} }

View File

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

View File

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

View File

@@ -264,10 +264,12 @@ namespace Opm
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_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, void computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state); const WellState& well_state,
DeferredLogger& deferred_logger);
void computePerfRateEval(const IntensiveQuantities& intQuants, void computePerfRateEval(const IntensiveQuantities& intQuants,
const std::vector<EvalWell>& mob, 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>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_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. // Verify that we have consistent input.
const int nperf = baseif_.numPerfs(); const int nperf = baseif_.numPerfs();
@@ -951,6 +952,18 @@ computeConnectionDensities(const std::vector<double>& perfComponentRates,
rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]);
} }
double d = 1.0 - rs*rv; 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());
}
if (rs > 0.0 && d > 0.0) { if (rs > 0.0 && d > 0.0) {
// Subtract gas in oil from gas mixture // Subtract gas in oil from gas mixture
x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d;

View File

@@ -144,7 +144,8 @@ protected:
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_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, ConvergenceReport getWellConvergence(const WellState& well_state,
const std::vector<double>& B_avg, const std::vector<double>& B_avg,

View File

@@ -312,15 +312,20 @@ namespace Opm
// Incorporate RS/RV factors if both oil and gas active // Incorporate RS/RV factors if both oil and gas active
const Value d = 1.0 - rv * rs; const Value d = 1.0 - rv * rs;
if (getValue(d) == 0.0) { if (d <= 0.0) {
OPM_DEFLOG_THROW(NumericalIssue, "Zero d value obtained for well " << this->name() << " during flux calcuation" std::ostringstream sstr;
<< " with rs " << rs << " and rv " << rv, deferred_logger); 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 = d > 0.0? (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d : cmix_s[oilCompIdx];
const Value tmp_oil = (cmix_s[oilCompIdx] - rv * cmix_s[gasCompIdx]) / d;
volumeRatio += tmp_oil / b_perfcells_dense[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]; volumeRatio += tmp_gas / b_perfcells_dense[gasCompIdx];
} }
else { else {
@@ -354,6 +359,17 @@ namespace Opm
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
const double d = 1.0 - getValue(rv) * getValue(rs); const double d = 1.0 - getValue(rv) * getValue(rs);
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());
}
// vaporized oil into gas // vaporized oil into gas
// rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d
if (d > 0.0) { if (d > 0.0) {
@@ -585,11 +601,28 @@ namespace Opm
// q_os = q_or * b_o + rv * q_gr * b_g // q_os = q_or * b_o + rv * q_gr * b_g
// q_gs = q_gr * g_g + rs * q_or * b_o // q_gs = q_gr * g_g + rs * q_or * b_o
// q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os)
if(FluidSystem::gasPhaseIdx == phaseIdx) // d = 1.0 - rs * rv
const EvalWell d = this->extendEval(1.0 - fs.Rv() * fs.Rs());
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());
}
if(FluidSystem::gasPhaseIdx == phaseIdx && d > 0.0) {
cq_r_thermal = (cq_s[gasCompIdx] - this->extendEval(fs.Rs()) * cq_s[oilCompIdx]) / (d * this->extendEval(fs.invB(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) // q_or = 1 / (b_o * d) * (q_os - rv * q_gs)
if(FluidSystem::oilPhaseIdx == phaseIdx) } else if(FluidSystem::oilPhaseIdx == phaseIdx && d > 0.0) {
cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) ); cq_r_thermal = (cq_s[oilCompIdx] - this->extendEval(fs.Rv()) * cq_s[gasCompIdx]) / (d * this->extendEval(fs.invB(phaseIdx)) );
} else {
assert(d <= 0.0);
// for d <= 0.0 we continue as if rs = 0 and rv = 0
cq_r_thermal = cq_s[activeCompIdx] / this->extendEval(fs.invB(phaseIdx));
}
} }
// change temperature for injecting fluids // change temperature for injecting fluids
@@ -1444,7 +1477,8 @@ namespace Opm
const std::vector<double>& b_perf, const std::vector<double>& b_perf,
const std::vector<double>& rsmax_perf, const std::vector<double>& rsmax_perf,
const std::vector<double>& rvmax_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 // Compute densities
const int nperf = this->number_of_perforations_; const int nperf = this->number_of_perforations_;
@@ -1499,7 +1533,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(); this->computeConnectionPressureDelta();
} }
@@ -1512,7 +1546,8 @@ namespace Opm
void void
StandardWell<TypeTag>:: StandardWell<TypeTag>::
computeWellConnectionPressures(const Simulator& ebosSimulator, computeWellConnectionPressures(const Simulator& ebosSimulator,
const WellState& well_state) const WellState& well_state,
DeferredLogger& deferred_logger)
{ {
// 1. Compute properties required by computeConnectionPressureDelta(). // 1. Compute properties required by computeConnectionPressureDelta().
// Note that some of the complexity of this part is due to the function // Note that some of the complexity of this part is due to the function
@@ -1522,7 +1557,7 @@ namespace Opm
std::vector<double> rvmax_perf; std::vector<double> rvmax_perf;
std::vector<double> surf_dens_perf; std::vector<double> surf_dens_perf;
computePropertiesForWellConnectionPressures(ebosSimulator, well_state, b_perf, rsmax_perf, rvmax_perf, 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);
} }
@@ -1558,7 +1593,7 @@ namespace Opm
{ {
updatePrimaryVariables(well_state, deferred_logger); updatePrimaryVariables(well_state, deferred_logger);
initPrimaryVariablesEvaluation(); initPrimaryVariablesEvaluation();
computeWellConnectionPressures(ebosSimulator, well_state); computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger);
this->computeAccumWell(); this->computeAccumWell();
} }
@@ -1727,7 +1762,7 @@ namespace Opm
deferred_logger.debug(msg); deferred_logger.debug(msg);
} }
well.updatePrimaryVariables(well_state_copy, deferred_logger); well.updatePrimaryVariables(well_state_copy, deferred_logger);
well.computeWellConnectionPressures(ebosSimulator, well_state_copy); well.computeWellConnectionPressures(ebosSimulator, well_state_copy, deferred_logger);
well.initPrimaryVariablesEvaluation(); well.initPrimaryVariablesEvaluation();
well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger); well.computeWellRatesWithBhp(ebosSimulator, bhp, well_flux, deferred_logger);
} }