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); }