diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp index 22fa6bfa8..acec5c1a8 100644 --- a/opm/autodiff/StandardWell.hpp +++ b/opm/autodiff/StandardWell.hpp @@ -281,7 +281,8 @@ namespace Opm void computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob_perfcells_dense, const double Tw, const EvalWell& bhp, const double& cdp, - const bool& allow_cf, std::vector& cq_s) const; + const bool& allow_cf, std::vector& cq_s, + double& perf_dis_gas_rate, double& perf_vap_oil_rate) const; // TODO: maybe we should provide a light version of computePerfRate, which does not include the // calculation of the derivatives diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp index db37cc05a..26a73b646 100644 --- a/opm/autodiff/StandardWell_impl.hpp +++ b/opm/autodiff/StandardWell_impl.hpp @@ -395,7 +395,8 @@ namespace Opm computePerfRate(const IntensiveQuantities& intQuants, const std::vector& mob_perfcells_dense, const double Tw, const EvalWell& bhp, const double& cdp, - const bool& allow_cf, std::vector& cq_s) const + const bool& allow_cf, std::vector& cq_s, + double& perf_dis_gas_rate, double& perf_vap_oil_rate) const { std::vector cmix_s(num_components_,0.0); for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { @@ -440,8 +441,17 @@ namespace Opm const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); const EvalWell cq_sOil = cq_s[oilCompIdx]; const EvalWell cq_sGas = cq_s[gasCompIdx]; - cq_s[gasCompIdx] += rs * cq_sOil; - cq_s[oilCompIdx] += rv * cq_sGas; + const EvalWell dis_gas = rs * cq_sOil; + const EvalWell vap_oil = rv * cq_sGas; + + cq_s[gasCompIdx] += dis_gas; + cq_s[oilCompIdx] += vap_oil; + + // recording the perforation solution gas rate and solution oil rates + if (well_type_ == PRODUCER) { + perf_dis_gas_rate = dis_gas.value(); + perf_vap_oil_rate = vap_oil.value(); + } } } else { @@ -506,6 +516,18 @@ namespace Opm for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is; // * b_perfcells_dense[phase]; } + + // calculating the perforation solution gas rate and solution oil rates + if (well_type_ == PRODUCER) { + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + const unsigned oilCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx); + const unsigned gasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx); + // TODO: not sure if this is the correct way + const double d = 1.0 - rv.value() * rs.value(); + perf_vap_oil_rate = (rv.value() * cq_s[gasCompIdx].value()) / d; + perf_dis_gas_rate = (rs.value() * cq_s[oilCompIdx].value()) / d; + } + } } } @@ -541,6 +563,10 @@ namespace Opm const EvalWell& bhp = getBhp(); + // the solution gas rate and solution oil rate needs to be updated for well_state. + std::fill(well_state.wellVaporizedOilRates().begin(), well_state.wellVaporizedOilRates().end(), 0.0); + std::fill(well_state.wellDissolvedGasRates().begin(), well_state.wellDissolvedGasRates().end(), 0.0); + for (int perf = 0; perf < number_of_perforations_; ++perf) { const int cell_idx = well_cells_[perf]; @@ -548,7 +574,16 @@ namespace Opm std::vector cq_s(num_components_,0.0); std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); - computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, + cq_s, perf_dis_gas_rate, perf_vap_oil_rate); + + // updating the solution gas rate and solution oil rate + if (well_type_ == PRODUCER) { + well_state.wellDissolvedGasRates()[index_of_well_] += perf_dis_gas_rate; + well_state.wellVaporizedOilRates()[index_of_well_] += perf_vap_oil_rate; + } for (int componentIdx = 0; componentIdx < num_components_; ++componentIdx) { // the cq_s entering mass balance equations need to consider the efficiency factors. @@ -1609,7 +1644,10 @@ namespace Opm std::vector cq_s(num_components_, 0.0); std::vector mob(num_components_, 0.0); getMobility(ebosSimulator, perf, mob); - computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + computePerfRate(intQuants, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, + cq_s, perf_dis_gas_rate, perf_vap_oil_rate); for(int p = 0; p < np; ++p) { well_flux[ebosCompIdxToFlowCompIdx(p)] += cq_s[p].value(); @@ -1995,7 +2033,10 @@ namespace Opm const bool allow_cf = crossFlowAllowed(ebos_simulator); const EvalWell& bhp = getBhp(); std::vector cq_s(num_components_,0.0); - computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, cq_s); + double perf_dis_gas_rate = 0.; + double perf_vap_oil_rate = 0.; + computePerfRate(int_quant, mob, well_index_[perf], bhp, perf_pressure_diffs_[perf], allow_cf, + cq_s, perf_dis_gas_rate, perf_vap_oil_rate); // TODO: make area a member const double area = 2 * M_PI * perf_rep_radius_[perf] * perf_length_[perf]; const auto& material_law_manager = ebos_simulator.problem().materialLawManager(); diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index d5ec6b868..2faba5e95 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -534,31 +534,16 @@ namespace Opm return solvent_well_rate; } - /* const std::vector& wellReservoirRates() const - { - return well_reservoir_rates_; - } */ - std::vector& wellReservoirRates() { return well_reservoir_rates_; } - /* const std::vector& wellDissolvedGasRates() const - { - return well_dissolved_gas_rates_; - } */ - std::vector& wellDissolvedGasRates() { return well_dissolved_gas_rates_; } - /* const std::vector& wellVaporizedOilRates() const - { - return well_vaporized_oil_rates_; - } */ - std::vector& wellVaporizedOilRates() { return well_vaporized_oil_rates_;