From 2e91d2f353a20d36c1393b880621b3cf28f432a9 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 30 Mar 2022 08:23:41 +0200 Subject: [PATCH 1/4] Do not throw when d = 0 instead continue with rs/rv = 0 --- opm/simulators/wells/MultisegmentWellEval.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 3f1267f2b..d56878359 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -1194,9 +1194,10 @@ 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 without 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? From 5b53fcd8a6738d4c73b8b0c94b216a41b7e0cc4e Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 30 Mar 2022 08:25:35 +0200 Subject: [PATCH 2/4] guard against d = 0 --- opm/simulators/wells/MultisegmentWellEval.cpp | 4 ++-- opm/simulators/wells/StandardWellEval.cpp | 9 +++++---- opm/simulators/wells/StandardWell_impl.hpp | 10 ++++++---- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index d56878359..352c3fff3 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -758,10 +758,10 @@ computeSegmentFluidProperties(const EvalWell& temperature, const EvalWell d = 1.0 - rs * rv; - if (rs != 0.0) { // rs > 0.0? + if (rs > 0.0 && d > 0.0) { mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; } - if (rv != 0.0) { // rv > 0.0? + if (rv > 0.0 && d > 0.0) { mix[oilCompIdx] = (mix_s[oilCompIdx] - mix_s[gasCompIdx] * rv) / d; } } diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 1ee1b54be..8e191e09b 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -950,13 +950,14 @@ computeConnectionDensities(const std::vector& perfComponentRates, if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); } - if (rs != 0.0) { + double d = 1.0 - rs*rv; + if (rs > 0.0 && d > 0.0) { // Subtract gas in oil from gas mixture - x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/(1.0 - rs*rv); + x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; } - if (rv != 0.0) { + if (rv > 0.0 && d > 0.0) { // Subtract oil in gas from oil mixture - x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/(1.0 - rs*rv);; + x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d; } } double volrat = 0.0; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 5e247b9ca..636fdf078 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -356,10 +356,14 @@ namespace Opm 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; + if (d > 0.0) { + 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) { + perf_dis_gas_rate = getValue(rs) * (getValue(cq_s[oilCompIdx]) - getValue(rv) * getValue(cq_s[gasCompIdx])) / d; + } } } } @@ -580,8 +584,6 @@ 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 - // 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)) ); From 69ffed06dee08388d83fe92909caeaff5271f07e Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 1 Apr 2022 08:35:31 +0200 Subject: [PATCH 3/4] Address comments from review 1) Add debug messages 2) Fix bug of missing else in the code --- opm/simulators/wells/MultisegmentWell.hpp | 3 +- opm/simulators/wells/MultisegmentWellEval.cpp | 41 +++++++----- opm/simulators/wells/MultisegmentWellEval.hpp | 3 +- .../wells/MultisegmentWell_impl.hpp | 7 ++- opm/simulators/wells/StandardWell.hpp | 6 +- opm/simulators/wells/StandardWellEval.cpp | 15 ++++- opm/simulators/wells/StandardWellEval.hpp | 3 +- opm/simulators/wells/StandardWell_impl.hpp | 63 ++++++++++++++----- 8 files changed, 103 insertions(+), 38 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 656e1f02d..097cf8115 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -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, diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 352c3fff3..72469cfbb 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -639,7 +639,8 @@ void MultisegmentWellEval:: computeSegmentFluidProperties(const EvalWell& temperature, const EvalWell& saltConcentration, - int pvt_region_index) + int pvt_region_index, + DeferredLogger& deferred_logger) { std::vector 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 && d > 0.0) { - mix[gasCompIdx] = (mix_s[gasCompIdx] - mix_s[oilCompIdx] * rs) / d; - } - if (rv > 0.0 && d > 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; + } } } @@ -1195,16 +1206,16 @@ getSegmentSurfaceVolume(const EvalWell& temperature, << " during conversion to surface volume with rs " << rs << ", rv " << rv << " and pressure " << seg_pressure << " 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."; OpmLog::debug(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; + } } } diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index 0940a7e4e..0654e1e77 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -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; diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 6603427bb..112f44e24 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -1000,7 +1000,7 @@ namespace Opm template void MultisegmentWell:: - 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; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index add189b0d..9cefbfdca 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -264,10 +264,12 @@ namespace Opm const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens_perf); + const std::vector& 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& mob, diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 8e191e09b..1986cbc68 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -831,7 +831,8 @@ computeConnectionDensities(const std::vector& perfComponentRates, const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens_perf) + const std::vector& surf_dens_perf, + DeferredLogger& deferred_logger) { // Verify that we have consistent input. const int nperf = baseif_.numPerfs(); @@ -951,6 +952,18 @@ computeConnectionDensities(const std::vector& perfComponentRates, rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); } 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) { // Subtract gas in oil from gas mixture x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; diff --git a/opm/simulators/wells/StandardWellEval.hpp b/opm/simulators/wells/StandardWellEval.hpp index ac114765f..27fdeeb43 100644 --- a/opm/simulators/wells/StandardWellEval.hpp +++ b/opm/simulators/wells/StandardWellEval.hpp @@ -144,7 +144,8 @@ protected: const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens_perf); + const std::vector& surf_dens_perf, + DeferredLogger& deferred_logger); ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector& B_avg, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 636fdf078..e161d84d2 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -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,6 +359,17 @@ namespace Opm // q_gr = 1 / (b_g * d) * (q_gs - rs * q_os) 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 // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / d if (d > 0.0) { @@ -585,11 +601,28 @@ namespace Opm // 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) - 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)) ); // 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)) ); + } 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 @@ -1444,7 +1477,8 @@ namespace Opm const std::vector& b_perf, const std::vector& rsmax_perf, const std::vector& rvmax_perf, - const std::vector& surf_dens_perf) + const std::vector& surf_dens_perf, + DeferredLogger& deferred_logger) { // Compute densities 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(); } @@ -1512,7 +1546,8 @@ namespace Opm void StandardWell:: 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 @@ -1522,7 +1557,7 @@ namespace Opm std::vector rvmax_perf; std::vector 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); initPrimaryVariablesEvaluation(); - computeWellConnectionPressures(ebosSimulator, well_state); + computeWellConnectionPressures(ebosSimulator, well_state, deferred_logger); this->computeAccumWell(); } @@ -1727,7 +1762,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); } From fba904620127037a2164f5776e04767bfc9df397 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 4 Apr 2022 11:37:31 +0200 Subject: [PATCH 4/4] Some more clean-up based on review --- opm/simulators/wells/StandardWellEval.cpp | 21 +++++----- opm/simulators/wells/StandardWell_impl.hpp | 45 ++++++++++------------ 2 files changed, 30 insertions(+), 36 deletions(-) diff --git a/opm/simulators/wells/StandardWellEval.cpp b/opm/simulators/wells/StandardWellEval.cpp index 1986cbc68..c57c6b83f 100644 --- a/opm/simulators/wells/StandardWellEval.cpp +++ b/opm/simulators/wells/StandardWellEval.cpp @@ -951,8 +951,7 @@ computeConnectionDensities(const std::vector& perfComponentRates, if (!rvmax_perf.empty() && mix[gaspos] > 1e-12) { rv = std::min(mix[oilpos]/mix[gaspos], rvmax_perf[perf]); } - double d = 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() @@ -962,15 +961,15 @@ computeConnectionDensities(const std::vector& perfComponentRates, << " 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) { - // Subtract gas in oil from gas mixture - x[gaspos] = (mix[gaspos] - mix[oilpos]*rs)/d; - } - if (rv > 0.0 && d > 0.0) { - // Subtract oil in gas from oil mixture - x[oilpos] = (mix[oilpos] - mix[gaspos]*rv)/d; + } 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; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index e161d84d2..362d13813 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -361,23 +361,20 @@ namespace Opm 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 - // rv * q_gr * b_g = rv * (q_gs - rs * q_os) / 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 - if (d > 0.0) { + // 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; } } @@ -612,16 +609,14 @@ namespace Opm << " 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)) ); - // q_or = 1 / (b_o * d) * (q_os - rv * q_gs) - } 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)) ); - } 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)); + } 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)) ); + } } }