From 7672e14890f9d3743eb68f8e41869d7f04a2e25f Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 08:34:02 +0100 Subject: [PATCH 1/6] RateConvert: introduce compile unit move global rate reductions into it --- CMakeLists_files.cmake | 1 + opm/simulators/wells/RateConverter.cpp | 65 ++++++++++++++++++++++++++ opm/simulators/wells/RateConverter.hpp | 29 ++++-------- 3 files changed, 75 insertions(+), 20 deletions(-) create mode 100644 opm/simulators/wells/RateConverter.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index e936a91fa..60b5c1603 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -105,6 +105,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/MultisegmentWellSegments.cpp opm/simulators/wells/ParallelWellInfo.cpp opm/simulators/wells/PerfData.cpp + opm/simulators/wells/RateConverter.cpp opm/simulators/wells/SegmentState.cpp opm/simulators/wells/SingleWellState.cpp opm/simulators/wells/StandardWellAssemble.cpp diff --git a/opm/simulators/wells/RateConverter.cpp b/opm/simulators/wells/RateConverter.cpp new file mode 100644 index 000000000..ad5902624 --- /dev/null +++ b/opm/simulators/wells/RateConverter.cpp @@ -0,0 +1,65 @@ +/* + Copyright 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Statoil ASA. + Copyright 2017, IRIS + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ +#include +#include + +#include +#include + +namespace Opm { +namespace RateConverter { + +template +void SurfaceToReservoirVoidage:: +sumRates(std::unordered_map& attributes_hpv, + std::unordered_map& attributes_pv, + Parallel::Communication comm) +{ + for (const auto& reg : rmap_.activeRegions()) { + // Calculating first using the hydrocarbon pore volumes + auto& ra = attr_.attributes(reg); + const auto& attri_hpv = attributes_hpv[reg]; + ra.data = attri_hpv.data; + comm.sum(ra.data.data(), ra.data.size()); + // TODO: should we have some epsilon here instead of zero? + // No hydrocarbon pore volume, redo but this time using full pore volumes. + if (ra.pv <= 0.) { + // using the pore volume to do the averaging + const auto& attri_pv = attributes_pv[reg]; + ra.data = attri_pv.data; + comm.sum(ra.data.data(), ra.data.size()); + assert(ra.pv > 0.); + } + const double pv_sum = ra.pv; + for (double& d : ra.data) + d /= pv_sum; + ra.pv = pv_sum; + } +} + +using FS = BlackOilFluidSystem; +template void SurfaceToReservoirVoidage>:: + sumRates(std::unordered_map&, + std::unordered_map&, + Parallel::Communication); + +} // namespace RateConverter +} // namespace Opm diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 562df3877..0b004b2a2 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -27,6 +27,8 @@ #include #include +#include + #include #include @@ -204,26 +206,9 @@ namespace Opm { OPM_END_PARALLEL_TRY_CATCH("SurfaceToReservoirVoidage::defineState() failed: ", simulator.vanguard().grid().comm()); - for (const auto& reg : rmap_.activeRegions()) { - // Calculating first using the hydrocarbon pore volumes - auto& ra = attr_.attributes(reg); - const auto& attri_hpv = attributes_hpv[reg]; - ra.data = attri_hpv.data; - comm.sum(ra.data.data(), ra.data.size()); - // TODO: should we have some epsilon here instead of zero? - // No hydrocarbon pore volume, redo but this time using full pore volumes. - if (ra.pv <= 0.) { - // using the pore volume to do the averaging - const auto& attri_pv = attributes_pv[reg]; - ra.data = attri_pv.data; - comm.sum(ra.data.data(), ra.data.size()); - assert(ra.pv > 0.); - } - const double pv_sum = ra.pv; - for (double& d : ra.data) - d /= pv_sum; - ra.pv = pv_sum; - } + this->sumRates(attributes_hpv, + attributes_pv, + comm); } /** @@ -625,6 +610,10 @@ namespace Opm { double& saltConcentration; }; + void sumRates(std::unordered_map& attributes_hpv, + std::unordered_map& attributes_pv, + Parallel::Communication comm); + RegionAttributeHelpers::RegionAttributes attr_; template From 11602bbdda3441fab68dcbbf37d7acb9a0628c16 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 08:49:25 +0100 Subject: [PATCH 2/6] RateConverter: move calcInjCoeff to compile unit --- opm/simulators/wells/RateConverter.cpp | 41 ++++++++++++++++++++++++++ opm/simulators/wells/RateConverter.hpp | 34 +-------------------- 2 files changed, 42 insertions(+), 33 deletions(-) diff --git a/opm/simulators/wells/RateConverter.cpp b/opm/simulators/wells/RateConverter.cpp index ad5902624..19f8a8da3 100644 --- a/opm/simulators/wells/RateConverter.cpp +++ b/opm/simulators/wells/RateConverter.cpp @@ -24,6 +24,8 @@ #include #include +#include + namespace Opm { namespace RateConverter { @@ -55,11 +57,50 @@ sumRates(std::unordered_map& attributes_hpv, } } +template +template +void SurfaceToReservoirVoidage:: +calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const +{ + const auto& pu = phaseUsage_; + const auto& ra = attr_.attributes(r); + const double p = ra.pressure; + const double T = ra.temperature; + const double saltConcentration = ra.saltConcentration; + + const int iw = RegionAttributeHelpers::PhasePos::water(pu); + const int io = RegionAttributeHelpers::PhasePos::oil (pu); + const int ig = RegionAttributeHelpers::PhasePos::gas (pu); + + std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); + + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + // q[w]_r = q[w]_s / bw + + const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration); + + coeff[iw] = 1.0 / bw; + } + + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0); + coeff[io] += 1.0 / bo; + } + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0); + coeff[ig] += 1.0 / bg; + } +} + using FS = BlackOilFluidSystem; template void SurfaceToReservoirVoidage>:: sumRates(std::unordered_map&, std::unordered_map&, Parallel::Communication); +template void SurfaceToReservoirVoidage>:: + calcInjCoeff>(const int, const int, std::vector&) const; + } // namespace RateConverter } // namespace Opm diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 0b004b2a2..0bd22a5a0 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -110,7 +110,6 @@ namespace Opm { ra.saltConcentration = 0.0; ra.rsw = 0.0; ra.rvw = 0.0; - } // quantities for pore volume average @@ -329,38 +328,7 @@ namespace Opm { template void - calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const - { - const auto& pu = phaseUsage_; - const auto& ra = attr_.attributes(r); - const double p = ra.pressure; - const double T = ra.temperature; - const double saltConcentration = ra.saltConcentration; - - const int iw = RegionAttributeHelpers::PhasePos::water(pu); - const int io = RegionAttributeHelpers::PhasePos::oil (pu); - const int ig = RegionAttributeHelpers::PhasePos::gas (pu); - - std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); - - if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - // q[w]_r = q[w]_s / bw - - const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, saltConcentration); - - coeff[iw] = 1.0 / bw; - } - - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0); - coeff[io] += 1.0 / bo; - } - - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0, 0.0); - coeff[ig] += 1.0 / bg; - } - } + calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const; /** * Convert surface volume flow rates to reservoir voidage flow From ae9ed506e807640d3bb32b7179ec7a362e6dd71e Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 08:49:25 +0100 Subject: [PATCH 3/6] RateConverter: move calcCoeff to compile unit --- opm/simulators/wells/RateConverter.cpp | 85 ++++++++++++++++++++++++++ opm/simulators/wells/RateConverter.hpp | 80 +----------------------- 2 files changed, 86 insertions(+), 79 deletions(-) diff --git a/opm/simulators/wells/RateConverter.cpp b/opm/simulators/wells/RateConverter.cpp index 19f8a8da3..e4ef72120 100644 --- a/opm/simulators/wells/RateConverter.cpp +++ b/opm/simulators/wells/RateConverter.cpp @@ -93,6 +93,89 @@ calcInjCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const } } +template +template +void SurfaceToReservoirVoidage:: +calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const +{ + const auto& pu = phaseUsage_; + const auto& ra = attr_.attributes(r); + const double p = ra.pressure; + const double T = ra.temperature; + const double saltConcentration = ra.saltConcentration; + + const int iw = RegionAttributeHelpers::PhasePos::water(pu); + const int io = RegionAttributeHelpers::PhasePos::oil (pu); + const int ig = RegionAttributeHelpers::PhasePos::gas (pu); + + std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); + + // Actual Rsw and Rvw: + double Rsw = ra.rsw; + double Rvw = ra.rvw; + // Determinant of 'R' matrix + const double detRw = 1.0 - (Rsw * Rvw); + + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) + + const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration); + + const double den = bw * detRw; + + coeff[iw] += 1.0 / den; + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + coeff[ig] -= ra.rvw / den; + } + } + + // Actual Rs and Rv: + double Rs = ra.rs; + double Rv = ra.rv; + + // Determinant of 'R' matrix + const double detR = 1.0 - (Rs * Rv); + + // Currently we only support either gas in water or gas in oil + // not both + if (detR != 1 && detRw != 1) { + std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); + throw std::range_error(msg); + } + + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) + + const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + const double den = bo * detR; + + coeff[io] += 1.0 / den; + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + coeff[ig] -= ra.rv / den; + } + } + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) + const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, Rvw); + if (FluidSystem::enableDissolvedGasInWater()) { + const double denw = bg * detRw; + coeff[ig] += 1.0 / denw; + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + coeff[iw] -= ra.rsw / denw; + } + } else { + const double den = bg * detR; + coeff[ig] += 1.0 / den; + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + coeff[io] -= ra.rs / den; + } + } + } +} + using FS = BlackOilFluidSystem; template void SurfaceToReservoirVoidage>:: sumRates(std::unordered_map&, @@ -101,6 +184,8 @@ template void SurfaceToReservoirVoidage>:: template void SurfaceToReservoirVoidage>:: calcInjCoeff>(const int, const int, std::vector&) const; +template void SurfaceToReservoirVoidage>:: + calcCoeff>(const int, const int, std::vector&) const; } // namespace RateConverter } // namespace Opm diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 0bd22a5a0..4d72d4f1a 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -246,85 +246,7 @@ namespace Opm { */ template void - calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const - { - const auto& pu = phaseUsage_; - const auto& ra = attr_.attributes(r); - const double p = ra.pressure; - const double T = ra.temperature; - const double saltConcentration = ra.saltConcentration; - - const int iw = RegionAttributeHelpers::PhasePos::water(pu); - const int io = RegionAttributeHelpers::PhasePos::oil (pu); - const int ig = RegionAttributeHelpers::PhasePos::gas (pu); - - std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); - - // Actual Rsw and Rvw: - double Rsw = ra.rsw; - double Rvw = ra.rvw; - // Determinant of 'R' matrix - const double detRw = 1.0 - (Rsw * Rvw); - - if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) - - const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rsw, saltConcentration); - - const double den = bw * detRw; - - coeff[iw] += 1.0 / den; - - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - coeff[ig] -= ra.rvw / den; - } - } - - // Actual Rs and Rv: - double Rs = ra.rs; - double Rv = ra.rv; - - // Determinant of 'R' matrix - const double detR = 1.0 - (Rs * Rv); - - // Currently we only support either gas in water or gas in oil - // not both - if (detR != 1 && detRw != 1) { - std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); - throw std::range_error(msg); - } - - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) - - const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); - const double den = bo * detR; - - coeff[io] += 1.0 / den; - - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - coeff[ig] -= ra.rv / den; - } - } - - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) - const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv, Rvw); - if (FluidSystem::enableDissolvedGasInWater()) { - const double denw = bg * detRw; - coeff[ig] += 1.0 / denw; - if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - coeff[iw] -= ra.rsw / denw; - } - } else { - const double den = bg * detR; - coeff[ig] += 1.0 / den; - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - coeff[io] -= ra.rs / den; - } - } - } - } + calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const; template void From 8c233ffba21acf2ca0562f4e41dba66df4c72192 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 08:49:25 +0100 Subject: [PATCH 4/6] RateConverter: move calcReservoirVoidageRates to compile unit --- opm/simulators/wells/RateConverter.cpp | 139 +++++++++++++++++++++++++ opm/simulators/wells/RateConverter.hpp | 93 +---------------- 2 files changed, 141 insertions(+), 91 deletions(-) diff --git a/opm/simulators/wells/RateConverter.cpp b/opm/simulators/wells/RateConverter.cpp index e4ef72120..5b7dc18ef 100644 --- a/opm/simulators/wells/RateConverter.cpp +++ b/opm/simulators/wells/RateConverter.cpp @@ -176,6 +176,117 @@ calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff& coeff) const } } +template +template +void SurfaceToReservoirVoidage:: +calcReservoirVoidageRates(const int pvtRegionIdx, + const double p, + const double rs, + const double rv, + const double rsw, + const double rvw, + const double T, + const double saltConcentration, + const SurfaceRates& surface_rates, + VoidageRates& voidage_rates) const +{ + const auto& pu = this->phaseUsage_; + const auto iw = RegionAttributeHelpers::PhasePos::water(pu); + const auto io = RegionAttributeHelpers::PhasePos::oil (pu); + const auto ig = RegionAttributeHelpers::PhasePos::gas (pu); + + const auto [Rs, Rv] = this-> + dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates); + + const auto [Rsw, Rvw] = this-> + dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates); + + + std::fill_n(&voidage_rates[0], pu.num_phases, 0.0); + + + // Determinant of 'R' matrix + const auto detRw = 1.0 - (Rsw * Rvw); + + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) + voidage_rates[iw] = surface_rates[iw]; + + const auto bw = FluidSystem::waterPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, + Rsw, + saltConcentration); + + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + voidage_rates[iw] -= Rvw * surface_rates[ig]; + } + voidage_rates[iw] /= bw * detRw; + } + + // Determinant of 'R' matrix + const auto detR = 1.0 - (Rs * Rv); + + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) + voidage_rates[io] = surface_rates[io]; + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + voidage_rates[io] -= Rv * surface_rates[ig]; + } + + const auto bo = FluidSystem::oilPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + + voidage_rates[io] /= bo * detR; + } + + // we only support either gas in water + // or gas in oil + if (detR != 1 && detRw != 1) { + std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); + throw std::range_error(msg); + } + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) + voidage_rates[ig] = surface_rates[ig]; + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + voidage_rates[ig] -= Rs * surface_rates[io]; + } + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + voidage_rates[ig] -= Rsw * surface_rates[iw]; + } + + const auto bg = FluidSystem::gasPvt() + .inverseFormationVolumeFactor(pvtRegionIdx, T, p, + Rv, Rvw); + + // we only support either gas in water or gas in oil + if (detRw == 1) { + voidage_rates[ig] /= bg * detR; + } else { + voidage_rates[ig] /= bg * detRw; + } + } +} + +template +template +void SurfaceToReservoirVoidage:: +calcReservoirVoidageRates(const RegionId r, + const int pvtRegionIdx, + const Rates& surface_rates, + Rates& voidage_rates) const +{ + const auto& ra = this->attr_.attributes(r); + + this->calcReservoirVoidageRates(pvtRegionIdx, + ra.pressure, ra.rs, ra.rv, + ra.rsw, ra.rvw, + ra.temperature, + ra.saltConcentration, + surface_rates, + voidage_rates); +} + using FS = BlackOilFluidSystem; template void SurfaceToReservoirVoidage>:: sumRates(std::unordered_map&, @@ -186,6 +297,34 @@ template void SurfaceToReservoirVoidage>:: calcInjCoeff>(const int, const int, std::vector&) const; template void SurfaceToReservoirVoidage>:: calcCoeff>(const int, const int, std::vector&) const; +template void SurfaceToReservoirVoidage>:: + calcReservoirVoidageRates,std::vector>(const int, + const double, + const double, + const double, + const double, + const double, + const double, + const double, + const std::vector&, + std::vector&) const; +template void SurfaceToReservoirVoidage>:: + calcReservoirVoidageRates(const int, + const double, + const double, + const double, + const double, + const double, + const double, + const double, + double const* const&, + double*&) const; + +template void SurfaceToReservoirVoidage>:: + calcReservoirVoidageRates>(const int, + const int, + const std::vector&, + std::vector&) const; } // namespace RateConverter } // namespace Opm diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 4d72d4f1a..f5a7e4c3d 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -278,18 +278,7 @@ namespace Opm { void calcReservoirVoidageRates(const RegionId r, const int pvtRegionIdx, const Rates& surface_rates, - Rates& voidage_rates) const - { - const auto& ra = this->attr_.attributes(r); - - this->calcReservoirVoidageRates(pvtRegionIdx, - ra.pressure, ra.rs, ra.rv, - ra.rsw, ra.rvw, - ra.temperature, - ra.saltConcentration, - surface_rates, - voidage_rates); - } + Rates& voidage_rates) const; /** * Convert surface volume flow rates to reservoir voidage flow @@ -335,85 +324,7 @@ namespace Opm { const double T, const double saltConcentration, const SurfaceRates& surface_rates, - VoidageRates& voidage_rates) const - { - const auto& pu = this->phaseUsage_; - const auto iw = RegionAttributeHelpers::PhasePos::water(pu); - const auto io = RegionAttributeHelpers::PhasePos::oil (pu); - const auto ig = RegionAttributeHelpers::PhasePos::gas (pu); - - const auto [Rs, Rv] = this-> - dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates); - - const auto [Rsw, Rvw] = this-> - dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates); - - - std::fill_n(&voidage_rates[0], pu.num_phases, 0.0); - - - // Determinant of 'R' matrix - const auto detRw = 1.0 - (Rsw * Rvw); - - if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - // q[w]_r = 1/(bw * (1 - rsw*rvw)) * (q[w]_s - rvw*q[g]_s) - voidage_rates[iw] = surface_rates[iw]; - - const auto bw = FluidSystem::waterPvt() - .inverseFormationVolumeFactor(pvtRegionIdx, T, p, - Rsw, - saltConcentration); - - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - voidage_rates[iw] -= Rvw * surface_rates[ig]; - } - voidage_rates[iw] /= bw * detRw; - } - - // Determinant of 'R' matrix - const auto detR = 1.0 - (Rs * Rv); - - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) - voidage_rates[io] = surface_rates[io]; - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - voidage_rates[io] -= Rv * surface_rates[ig]; - } - - const auto bo = FluidSystem::oilPvt() - .inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); - - voidage_rates[io] /= bo * detR; - } - - // we only support either gas in water - // or gas in oil - if (detR != 1 && detRw != 1) { - std::string msg = "only support " + std::to_string(detR) + " " + std::to_string(detR); - throw std::range_error(msg); - } - if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { - // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) - voidage_rates[ig] = surface_rates[ig]; - if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { - voidage_rates[ig] -= Rs * surface_rates[io]; - } - if (RegionAttributeHelpers::PhaseUsed::water(pu)) { - voidage_rates[ig] -= Rsw * surface_rates[iw]; - } - - const auto bg = FluidSystem::gasPvt() - .inverseFormationVolumeFactor(pvtRegionIdx, T, p, - Rv, Rvw); - - // we only support either gas in water or gas in oil - if (detRw == 1) { - voidage_rates[ig] /= bg * detR; - } else { - voidage_rates[ig] /= bg * detRw; - } - } - } + VoidageRates& voidage_rates) const; template std::pair From 173e3900d497aaa37b254f754d342cb40f452676 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 08:49:25 +0100 Subject: [PATCH 5/6] RateConverter: move inferDissolvedVaporisedRatio to compile unit --- opm/simulators/wells/RateConverter.cpp | 53 +++++++++++++++++++++++++- opm/simulators/wells/RateConverter.hpp | 35 +---------------- 2 files changed, 52 insertions(+), 36 deletions(-) diff --git a/opm/simulators/wells/RateConverter.cpp b/opm/simulators/wells/RateConverter.cpp index 5b7dc18ef..f36a48005 100644 --- a/opm/simulators/wells/RateConverter.cpp +++ b/opm/simulators/wells/RateConverter.cpp @@ -25,6 +25,36 @@ #include #include +#include +#include +#include + +namespace { + +template +std::pair +dissolvedVaporisedRatio(const int io, + const int ig, + const double rs, + const double rv, + const Rates& surface_rates) +{ + if ((io < 0) || (ig < 0)) { + return { rs, rv }; + } + auto eps = std::copysign(1.0e-15, surface_rates[io]); + const auto Rs = surface_rates[ig] / (surface_rates[io] + eps); + + eps = std::copysign(1.0e-15, surface_rates[ig]); + const auto Rv = surface_rates[io] / (surface_rates[ig] + eps); + + return { + std::clamp(static_cast(Rs), 0.0, rs), + std::clamp(static_cast(Rv), 0.0, rv) + }; +} + +} namespace Opm { namespace RateConverter { @@ -195,10 +225,10 @@ calcReservoirVoidageRates(const int pvtRegionIdx, const auto io = RegionAttributeHelpers::PhasePos::oil (pu); const auto ig = RegionAttributeHelpers::PhasePos::gas (pu); - const auto [Rs, Rv] = this-> + const auto [Rs, Rv] = dissolvedVaporisedRatio(io, ig, rs, rv, surface_rates); - const auto [Rsw, Rvw] = this-> + const auto [Rsw, Rvw] = dissolvedVaporisedRatio(iw, ig, rsw, rvw, surface_rates); @@ -287,6 +317,19 @@ calcReservoirVoidageRates(const RegionId r, voidage_rates); } +template +template +std::pair +SurfaceToReservoirVoidage:: +inferDissolvedVaporisedRatio(const double rsMax, + const double rvMax, + const Rates& surface_rates) const +{ + const auto io = RegionAttributeHelpers::PhasePos::oil(this->phaseUsage_); + const auto ig = RegionAttributeHelpers::PhasePos::gas(this->phaseUsage_); + return dissolvedVaporisedRatio(io, ig, rsMax, rvMax, surface_rates); +} + using FS = BlackOilFluidSystem; template void SurfaceToReservoirVoidage>:: sumRates(std::unordered_map&, @@ -326,5 +369,11 @@ template void SurfaceToReservoirVoidage>:: const std::vector&, std::vector&) const; +template std::pair +SurfaceToReservoirVoidage>:: + inferDissolvedVaporisedRatio::iterator>(const double, + const double, + const std::vector::iterator&) const; + } // namespace RateConverter } // namespace Opm diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index f5a7e4c3d..dbcba62ca 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -32,13 +32,8 @@ #include #include -#include #include #include -#include -#include -#include -#include #include #include #include @@ -330,12 +325,7 @@ namespace Opm { std::pair inferDissolvedVaporisedRatio(const double rsMax, const double rvMax, - const Rates& surface_rates) const - { - const auto io = RegionAttributeHelpers::PhasePos::oil(this->phaseUsage_); - const auto ig = RegionAttributeHelpers::PhasePos::gas(this->phaseUsage_); - return this->dissolvedVaporisedRatio(io, ig, rsMax, rvMax, surface_rates); - } + const Rates& surface_rates) const; /** * Compute coefficients for surface-to-reservoir voidage @@ -416,29 +406,6 @@ namespace Opm { Parallel::Communication comm); RegionAttributeHelpers::RegionAttributes attr_; - - template - std::pair - dissolvedVaporisedRatio(const int io, - const int ig, - const double rs, - const double rv, - const Rates& surface_rates) const - { - if ((io < 0) || (ig < 0)) { - return { rs, rv }; - } - auto eps = std::copysign(1.0e-15, surface_rates[io]); - const auto Rs = surface_rates[ig] / (surface_rates[io] + eps); - - eps = std::copysign(1.0e-15, surface_rates[ig]); - const auto Rv = surface_rates[io] / (surface_rates[ig] + eps); - - return { - std::clamp(static_cast(Rs), 0.0, rs), - std::clamp(static_cast(Rv), 0.0, rv) - }; - } }; } // namespace RateConverter From 3967eb001ba6792cd129b0e996be9c9647f90650 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 6 Mar 2023 09:34:23 +0100 Subject: [PATCH 6/6] RateConverter: prefer using --- opm/simulators/wells/RateConverter.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index dbcba62ca..3ddbba9ae 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -210,7 +210,7 @@ namespace Opm { * * Integral type. */ - typedef typename RegionMapping::RegionId RegionId; + using RegionId = typename RegionMapping::RegionId; /** * Compute coefficients for surface-to-reservoir voidage