From d1c1aecac71bf5f161199fa690a0d2e57a6da735 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 21 Nov 2022 13:20:34 +0100 Subject: [PATCH] move implementation of computePropertiesForWellConnectionPressure to StandardWellConnections --- .../wells/StandardWellConnections.cpp | 172 ++++++++++++++++ .../wells/StandardWellConnections.hpp | 13 ++ opm/simulators/wells/StandardWell_impl.hpp | 185 ++++-------------- 3 files changed, 221 insertions(+), 149 deletions(-) diff --git a/opm/simulators/wells/StandardWellConnections.cpp b/opm/simulators/wells/StandardWellConnections.cpp index 6d75e6185..925a359e3 100644 --- a/opm/simulators/wells/StandardWellConnections.cpp +++ b/opm/simulators/wells/StandardWellConnections.cpp @@ -22,6 +22,8 @@ #include #include +#include + #include #include @@ -31,6 +33,7 @@ #include #include #include +#include #include @@ -273,6 +276,175 @@ computeDensities(const std::vector& perfComponentRates, } } +template +void StandardWellConnections:: +computePropertiesForWellConnectionPressures(const WellState& well_state, + const std::function& getTemperature, + const std::function& getSaltConcentration, + const std::function& pvtRegionIdx, + const std::function& solventInverseFormationVolumeFactor, + const std::function& solventRefDensity, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& rvwmax_perf, + std::vector& surf_dens_perf) const +{ + const int nperf = well_.numPerfs(); + const PhaseUsage& pu = well_.phaseUsage(); + b_perf.resize(nperf * well_.numComponents()); + surf_dens_perf.resize(nperf * well_.numComponents()); + const auto& ws = well_state.well(well_.indexOfWell()); + + static constexpr int Water = BlackoilPhases::Aqua; + static constexpr int Oil = BlackoilPhases::Liquid; + static constexpr int Gas = BlackoilPhases::Vapour; + + const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); + const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + + //rs and rv are only used if both oil and gas is present + if (oilPresent && gasPresent) { + rsmax_perf.resize(nperf); + rvmax_perf.resize(nperf); + } + //rvw is only used if both water and gas is present + if (waterPresent && gasPresent) { + rvwmax_perf.resize(nperf); + } + + // Compute the average pressure in each well block + const auto& perf_press = ws.perf_data.pressure; + auto p_above = well_.parallelWellInfo().communicateAboveValues(ws.bhp, + perf_press.data(), + nperf); + + for (int perf = 0; perf < nperf; ++perf) { + const int cell_idx = well_.cells()[perf]; + + const double p_avg = (perf_press[perf] + p_above[perf])/2; + const double temperature = getTemperature(cell_idx, FluidSystem::oilPhaseIdx); + const double saltConcentration = getSaltConcentration(cell_idx); + const int region_idx = pvtRegionIdx(cell_idx); + + if (waterPresent) { + const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); + b_perf[ waterCompIdx + perf * well_.numComponents()] = + FluidSystem::waterPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, saltConcentration); + } + + if (gasPresent) { + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + const int gaspos = gasCompIdx + perf * well_.numComponents(); + + if (oilPresent && waterPresent) { + const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers + const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(region_idx, temperature, p_avg); + rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(region_idx, temperature, p_avg); + double rv = 0.0; + double rvw = 0.0; + if (oilrate > 0) { + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + } + if (waterrate > 0) { + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + if (gasrate > 0) { + rvw = waterrate / gasrate; + } + rvw = std::min(rvw, rvwmax_perf[perf]); + } + if (rv > 0.0 || rvw > 0.0){ + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rv, rvw); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + } else if (oilPresent) { + //no water + const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers + rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(region_idx, temperature, p_avg); + if (oilrate > 0) { + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + double rv = 0.0; + if (gasrate > 0) { + rv = oilrate / gasrate; + } + rv = std::min(rv, rvmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rv, 0.0 /*Rvw*/); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + } else if (waterPresent) { + //no oil + const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers + rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(region_idx, temperature, p_avg); + if (waterrate > 0) { + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + double rvw = 0.0; + if (gasrate > 0) { + rvw = waterrate / gasrate; + } + rvw = std::min(rvw, rvwmax_perf[perf]); + + b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, 0.0 /*Rv*/, rvw); + } + else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + + } else { + b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + } + + if (oilPresent) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const int oilpos = oilCompIdx + perf * well_.numComponents(); + if (gasPresent) { + rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(region_idx, temperature, p_avg); + const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (Indices::enableSolvent ? ws.sum_solvent_rates() : 0.0); + if (gasrate > 0) { + const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); + double rs = 0.0; + if (oilrate > 0) { + rs = gasrate / oilrate; + } + rs = std::min(rs, rsmax_perf[perf]); + b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(region_idx, temperature, p_avg, rs); + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + } else { + b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(region_idx, temperature, p_avg); + } + } + + // Surface density. + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + surf_dens_perf[well_.numComponents() * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, region_idx ); + } + + // We use cell values for solvent injector + if constexpr (Indices::enableSolvent) { + b_perf[well_.numComponents() * perf + Indices::contiSolventEqIdx] = solventInverseFormationVolumeFactor(cell_idx); + surf_dens_perf[well_.numComponents() * perf + Indices::contiSolventEqIdx] = solventRefDensity(cell_idx); + } + } +} + #define INSTANCE(...) \ template class StandardWellConnections,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/StandardWellConnections.hpp b/opm/simulators/wells/StandardWellConnections.hpp index 2b427577b..73ee9fa67 100644 --- a/opm/simulators/wells/StandardWellConnections.hpp +++ b/opm/simulators/wells/StandardWellConnections.hpp @@ -23,6 +23,7 @@ #ifndef OPM_STANDARDWELL_CONNECTIONS_HEADER_INCLUDED #define OPM_STANDARDWELL_CONNECTIONS_HEADER_INCLUDED +#include #include namespace Opm @@ -49,6 +50,18 @@ public: const std::vector& surf_dens_perf, DeferredLogger& deferred_logger); + void computePropertiesForWellConnectionPressures(const WellState& well_state, + const std::function& getTemperature, + const std::function& getSaltConcentration, + const std::function& pvtRegionIdx, + const std::function& solventInverseFormationVolumeFactor, + const std::function& solventRefDensity, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& rvwmax_perf, + std::vector& surf_dens_perf) const; + Scalar getRho() const { return this->perf_densities_.empty() ? 0.0 : perf_densities_[0]; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c21fadd6c..77db7b5f9 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1288,156 +1288,43 @@ namespace Opm std::vector& rvwmax_perf, std::vector& surf_dens_perf) const { - const int nperf = this->number_of_perforations_; - const PhaseUsage& pu = phaseUsage(); - b_perf.resize(nperf * this->num_components_); - surf_dens_perf.resize(nperf * this->num_components_); - const auto& ws = well_state.well(this->index_of_well_); + std::function getTemperature = + [&ebosSimulator](int cell_idx, int phase_idx) + { + return ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0)->fluidState().temperature(phase_idx).value(); + }; + std::function getSaltConcentration = + [&ebosSimulator](int cell_idx) + { + return ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0)->fluidState().saltConcentration().value(); + }; + std::function getPvtRegionIdx = + [&ebosSimulator](int cell_idx) + { + return ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0)->fluidState().pvtRegionIndex(); + }; + std::function getInvFac = + [&ebosSimulator](int cell_idx) + { + return ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0)->solventInverseFormationVolumeFactor().value(); + }; + std::function getSolventDensity = + [&ebosSimulator](int cell_idx) + { + return ebosSimulator.model().cachedIntensiveQuantities(cell_idx, 0)->solventRefDensity(); + }; - const bool waterPresent = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool oilPresent = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const bool gasPresent = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - - //rs and rv are only used if both oil and gas is present - if (oilPresent && gasPresent) { - rsmax_perf.resize(nperf); - rvmax_perf.resize(nperf); - } - //rvw is only used if both water and gas is present - if (waterPresent && gasPresent) { - rvwmax_perf.resize(nperf); - } - - // Compute the average pressure in each well block - const auto& perf_press = ws.perf_data.pressure; - auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp, - perf_press.data(), - nperf); - - for (int perf = 0; perf < nperf; ++perf) { - const int cell_idx = this->well_cells_[perf]; - const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0)); - const auto& fs = intQuants.fluidState(); - - const double p_avg = (perf_press[perf] + p_above[perf])/2; - const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); - const double saltConcentration = fs.saltConcentration().value(); - - if (waterPresent) { - const unsigned waterCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx); - b_perf[ waterCompIdx + perf * this->num_components_] = - FluidSystem::waterPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, saltConcentration); - } - - if (gasPresent) { - const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); - const int gaspos = gasCompIdx + perf * this->num_components_; - - if (oilPresent && waterPresent) { - const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers - const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers - rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - double rv = 0.0; - double rvw = 0.0; - if (oilrate > 0) { - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); - if (gasrate > 0) { - rv = oilrate / gasrate; - } - rv = std::min(rv, rvmax_perf[perf]); - } - if (waterrate > 0) { - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); - if (gasrate > 0) { - rvw = waterrate / gasrate; - } - rvw = std::min(rvw, rvwmax_perf[perf]); - } - if (rv > 0.0 || rvw > 0.0){ - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, rvw); - } - else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } else if (oilPresent) { - //no water - const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers - rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - if (oilrate > 0) { - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); - double rv = 0.0; - if (gasrate > 0) { - rv = oilrate / gasrate; - } - rv = std::min(rv, rvmax_perf[perf]); - - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rv, 0.0 /*Rvw*/); - } - else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } else if (waterPresent) { - //no oil - const double waterrate = std::abs(ws.surface_rates[pu.phase_pos[Water]]); //in order to handle negative rates in producers - rvwmax_perf[perf] = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg); - if (waterrate > 0) { - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); - double rvw = 0.0; - if (gasrate > 0) { - rvw = waterrate / gasrate; - } - rvw = std::min(rvw, rvwmax_perf[perf]); - - b_perf[gaspos] = FluidSystem::gasPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, 0.0 /*Rv*/, rvw); - } - else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - - } else { - b_perf[gaspos] = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - if (oilPresent) { - const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); - const int oilpos = oilCompIdx + perf * this->num_components_; - if (gasPresent) { - rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg); - const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? ws.sum_solvent_rates() : 0.0); - if (gasrate > 0) { - const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); - double rs = 0.0; - if (oilrate > 0) { - rs = gasrate / oilrate; - } - rs = std::min(rs, rsmax_perf[perf]); - b_perf[oilpos] = FluidSystem::oilPvt().inverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg, rs); - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } else { - b_perf[oilpos] = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(fs.pvtRegionIndex(), temperature, p_avg); - } - } - - // Surface density. - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned compIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - surf_dens_perf[this->num_components_ * perf + compIdx] = FluidSystem::referenceDensity( phaseIdx, fs.pvtRegionIndex() ); - } - - // We use cell values for solvent injector - if constexpr (has_solvent) { - b_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventInverseFormationVolumeFactor().value(); - surf_dens_perf[this->num_components_ * perf + Indices::contiSolventEqIdx] = intQuants.solventRefDensity(); - } - } + this->connections_.computePropertiesForWellConnectionPressures(well_state, + getTemperature, + getSaltConcentration, + getPvtRegionIdx, + getInvFac, + getSolventDensity, + b_perf, + rsmax_perf, + rvmax_perf, + rvwmax_perf, + surf_dens_perf); }