From 35c56e4ce48a6c09aa024e7f8c89e89b1b4fa736 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Sun, 7 May 2023 21:56:09 +0200 Subject: [PATCH 1/5] changed: unify StandardWell::getMobility(Eval|Scalar) --- opm/simulators/wells/StandardWell.hpp | 16 +-- opm/simulators/wells/StandardWell_impl.hpp | 123 ++++++--------------- 2 files changed, 37 insertions(+), 102 deletions(-) diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index 95f66a987..d1ad5820b 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -333,17 +333,11 @@ namespace Opm virtual double getRefDensity() const override; // get the mobility for specific perforation - void getMobilityEval(const Simulator& ebosSimulator, - const int perf, - std::vector& mob, - DeferredLogger& deferred_logger) const; - - // get the mobility for specific perforation - void getMobilityScalar(const Simulator& ebosSimulator, - const int perf, - std::vector& mob, - DeferredLogger& deferred_logger) const; - + template + void getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob, + DeferredLogger& deferred_logger) const; void updateWaterMobilityWithPolymer(const Simulator& ebos_simulator, const int perf, diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index a16e07a31..c78ad025a 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -618,7 +618,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.}); - getMobilityEval(ebosSimulator, perf, mob, deferred_logger); + getMobility(ebosSimulator, perf, mob, deferred_logger); PerforationRates perf_rates; double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); @@ -820,15 +820,23 @@ namespace Opm - template + template void StandardWell:: - getMobilityEval(const Simulator& ebosSimulator, - const int perf, - std::vector& mob, - DeferredLogger& deferred_logger) const + getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob, + DeferredLogger& deferred_logger) const { + auto obtain = [this](const Eval& value) + { + if constexpr (std::is_same_v) { + return getValue(value); + } else { + return this->extendEval(value); + } + }; const int cell_idx = this->well_cells_[perf]; assert (int(mob.size()) == this->num_components_); const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); @@ -838,21 +846,19 @@ namespace Opm // based on passing the saturation table index const int satid = this->saturation_table_number_[perf] - 1; const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell - + if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { continue; } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx)); + mob[activeCompIdx] = obtain(intQuants.mobility(phaseIdx)); } if (has_solvent) { - mob[Indices::contiSolventEqIdx] = this->extendEval(intQuants.solventMobility()); + mob[Indices::contiSolventEqIdx] = obtain(intQuants.solventMobility()); } } else { - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); std::array relativePerms = { 0.0, 0.0, 0.0 }; MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); @@ -867,7 +873,7 @@ namespace Opm } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); + mob[activeCompIdx] = obtain(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); } // this may not work if viscosity and relperms has been modified? @@ -885,81 +891,16 @@ namespace Opm // for the cases related to polymer molecular weight, we assume fully mixing // as a result, the polymer and water share the same viscosity if constexpr (!Base::has_polymermw) { - updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger); - } - } - } - - template - void - StandardWell:: - getMobilityScalar(const Simulator& ebosSimulator, - const int perf, - std::vector& mob, - DeferredLogger& deferred_logger) const - { - const int cell_idx = this->well_cells_[perf]; - assert (int(mob.size()) == this->num_components_); - const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - - // either use mobility of the perforation cell or calcualte its own - // based on passing the saturation table index - const int satid = this->saturation_table_number_[perf] - 1; - const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell - - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx)); - } - if (has_solvent) { - mob[Indices::contiSolventEqIdx] = getValue(intQuants.solventMobility()); - } - } else { - - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - std::array relativePerms = { 0.0, 0.0, 0.0 }; - MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); - - // reset the satnumvalue back to original - materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); - - // compute the mobility - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = getValue(relativePerms[phaseIdx]) / getValue(intQuants.fluidState().viscosity(phaseIdx)); - } - - // this may not work if viscosity and relperms has been modified? - if constexpr (has_solvent) { - OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger); - } - } - - // modify the water mobility if polymer is present - if constexpr (has_polymer) { - if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { - OPM_DEFLOG_THROW(std::runtime_error, "Water is required when polymer is active", deferred_logger); - } - - // for the cases related to polymer molecular weight, we assume fully mixing - // as a result, the polymer and water share the same viscosity - if constexpr (!Base::has_polymermw) { - std::vector mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.}); - for (size_t i = 0; i < mob.size(); ++i) - mob_eval[i].setValue(mob[i]); - updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger); - for (size_t i = 0; i < mob.size(); ++i) { - mob[i] = getValue(mob_eval[i]); + if constexpr (std::is_same_v) { + std::vector mob_eval(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.}); + for (size_t i = 0; i < mob.size(); ++i) + mob_eval[i].setValue(mob[i]); + updateWaterMobilityWithPolymer(ebosSimulator, perf, mob_eval, deferred_logger); + for (size_t i = 0; i < mob.size(); ++i) { + mob[i] = getValue(mob_eval[i]); + } + } else { + updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger); } } } @@ -1043,7 +984,7 @@ namespace Opm for (int perf = 0; perf < this->number_of_perforations_; ++perf) { std::vector mob(this->num_components_, 0.0); - getMobilityScalar(ebos_simulator, perf, mob, deferred_logger); + getMobility(ebos_simulator, perf, mob, deferred_logger); const int cell_idx = this->well_cells_[perf]; const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); @@ -1450,7 +1391,7 @@ namespace Opm }; std::vector mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0}); - getMobilityEval(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); const auto& fs = fluidState(subsetPerfID); setToZero(connPI); @@ -1652,7 +1593,7 @@ namespace Opm const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); // flux for each perforation std::vector mob(this->num_components_, 0.); - getMobilityScalar(ebosSimulator, perf, mob, deferred_logger); + getMobility(ebosSimulator, perf, mob, deferred_logger); double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; @@ -2509,7 +2450,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.); - getMobilityScalar(ebosSimulator, perf, mob, deferred_logger); + getMobility(ebosSimulator, perf, mob, deferred_logger); std::vector cq_s(this->num_components_, 0.); double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; From 5126097d7b7223b42e6ecf957a28ad047c136ce3 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Sun, 7 May 2023 21:56:09 +0200 Subject: [PATCH 2/5] changed: unify MultisegmentWell::getMobility(Eval|Scalar) --- opm/simulators/wells/MultisegmentWell.hpp | 12 +-- .../wells/MultisegmentWell_impl.hpp | 100 +++++------------- 2 files changed, 33 insertions(+), 79 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index c22b9a77c..18c74d9e1 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -235,14 +235,10 @@ namespace Opm DeferredLogger& deferred_logger); // get the mobility for specific perforation - void getMobilityEval(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const; - - // get the mobility for specific perforation - void getMobilityScalar(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const; + template + void getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const; void computeWellRatesAtBhpLimit(const Simulator& ebosSimulator, std::vector& well_flux, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 3293c8bc8..2e68fe2d5 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -389,7 +389,7 @@ namespace Opm const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); // flux for each perforation std::vector mob(this->num_components_, 0.); - getMobilityScalar(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob); double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; @@ -716,7 +716,7 @@ namespace Opm }; std::vector mob(this->num_components_, 0.0); - getMobilityScalar(ebosSimulator, static_cast(subsetPerfID), mob); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob); const auto& fs = fluidState(subsetPerfID); setToZero(connPI); @@ -1066,17 +1066,30 @@ namespace Opm deferred_logger); } - - - - template + template void MultisegmentWell:: - getMobilityEval(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const + getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob) const { + auto obtain = [this](const Eval& value) + { + if constexpr (std::is_same_v) { + return getValue(value); + } else { + return this->extendEval(value); + } + }; + auto relpermArray = []() + { + if constexpr (std::is_same_v) { + return std::array{}; + } else { + return std::array{}; + } + }; // TODO: most of this function, if not the whole function, can be moved to the base class const int cell_idx = this->well_cells_[perf]; assert (int(mob.size()) == this->num_components_); @@ -1087,23 +1100,21 @@ namespace Opm // based on passing the saturation table index const int satid = this->saturation_table_number_[perf] - 1; const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell - + if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { if (!FluidSystem::phaseIsActive(phaseIdx)) { continue; } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = this->extendEval(intQuants.mobility(phaseIdx)); + mob[activeCompIdx] = obtain(intQuants.mobility(phaseIdx)); } // if (has_solvent) { // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); // } } else { - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - std::array relativePerms = { 0.0, 0.0, 0.0 }; + auto relativePerms = relpermArray(); MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); // reset the satnumvalue back to original @@ -1116,65 +1127,12 @@ namespace Opm } const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); + mob[activeCompIdx] = obtain(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); } } } - template - void - MultisegmentWell:: - getMobilityScalar(const Simulator& ebosSimulator, - const int perf, - std::vector& mob) const - { - // TODO: most of this function, if not the whole function, can be moved to the base class - const int cell_idx = this->well_cells_[perf]; - assert (int(mob.size()) == this->num_components_); - const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - - // either use mobility of the perforation cell or calcualte its own - // based on passing the saturation table index - const int satid = this->saturation_table_number_[perf] - 1; - const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if( satid == satid_elem ) { // the same saturation number is used. i.e. just use the mobilty from the cell - - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = getValue(intQuants.mobility(phaseIdx)); - } - // if (has_solvent) { - // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); - // } - } else { - - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - std::array relativePerms = { 0.0, 0.0, 0.0 }; - MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); - - // reset the satnumvalue back to original - materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); - - // compute the mobility - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx)); - } - } - } - - - template double @@ -1270,7 +1228,7 @@ namespace Opm std::vector mob(this->num_components_, 0.0); // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration - getMobilityScalar(ebos_simulator, perf, mob); + getMobility(ebos_simulator, perf, mob); const int cell_idx = this->well_cells_[perf]; const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); @@ -1628,7 +1586,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobilityEval(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); @@ -1940,7 +1898,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobilityScalar(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); From c0bc0abc73e8c680ef0170c2aa43ab3f490f21db Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 10 May 2023 14:55:17 +0200 Subject: [PATCH 3/5] StandardWell::getMobility: use Scalar relperms when possible --- opm/simulators/wells/StandardWell_impl.hpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index c78ad025a..e7ba66402 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -837,6 +837,15 @@ namespace Opm return this->extendEval(value); } }; + auto relpermArray = []() + { + if constexpr (std::is_same_v) { + return std::array{}; + } else { + return std::array{}; + } + }; + const int cell_idx = this->well_cells_[perf]; assert (int(mob.size()) == this->num_components_); const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); @@ -860,7 +869,7 @@ namespace Opm } } else { const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - std::array relativePerms = { 0.0, 0.0, 0.0 }; + auto relativePerms = relpermArray(); MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); // reset the satnumvalue back to original From e406d2f0a1ad7e6e803ca11a862399503905067e Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 10 May 2023 14:56:09 +0200 Subject: [PATCH 4/5] StandardWell::getMobility: use if constexpr --- opm/simulators/wells/StandardWell_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index e7ba66402..d9bf534d9 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -864,7 +864,7 @@ namespace Opm const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); mob[activeCompIdx] = obtain(intQuants.mobility(phaseIdx)); } - if (has_solvent) { + if constexpr (has_solvent) { mob[Indices::contiSolventEqIdx] = obtain(intQuants.solventMobility()); } } else { From 0406a033a282f60e75f655c587f4d839ea4a0f5b Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Fri, 12 May 2023 10:02:05 +0200 Subject: [PATCH 5/5] introduce getMobility in WellInterface we now share implementation between StandardWell and MultisegmentWell --- opm/simulators/wells/MultisegmentWell.hpp | 3 +- .../wells/MultisegmentWell_impl.hpp | 64 +++---------------- opm/simulators/wells/StandardWell_impl.hpp | 57 +---------------- opm/simulators/wells/WellInterface.hpp | 13 ++-- opm/simulators/wells/WellInterface_impl.hpp | 64 +++++++++++++++++++ 5 files changed, 87 insertions(+), 114 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 18c74d9e1..50ff0f629 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -238,7 +238,8 @@ namespace Opm template void getMobility(const Simulator& ebosSimulator, const int perf, - std::vector& mob) const; + std::vector& mob, + DeferredLogger& deferred_logger) const; void computeWellRatesAtBhpLimit(const Simulator& ebosSimulator, std::vector& well_flux, diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 2e68fe2d5..c629ead1a 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -389,7 +389,7 @@ namespace Opm const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); // flux for each perforation std::vector mob(this->num_components_, 0.); - getMobility(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob, deferred_logger); double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(intQuants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; @@ -716,7 +716,7 @@ namespace Opm }; std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, static_cast(subsetPerfID), mob); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); const auto& fs = fluidState(subsetPerfID); setToZero(connPI); @@ -1072,64 +1072,18 @@ namespace Opm MultisegmentWell:: getMobility(const Simulator& ebosSimulator, const int perf, - std::vector& mob) const + std::vector& mob, + DeferredLogger& deferred_logger) const { auto obtain = [this](const Eval& value) { if constexpr (std::is_same_v) { return getValue(value); } else { - return this->extendEval(value); + return this->extendEval(value); } }; - auto relpermArray = []() - { - if constexpr (std::is_same_v) { - return std::array{}; - } else { - return std::array{}; - } - }; - // TODO: most of this function, if not the whole function, can be moved to the base class - const int cell_idx = this->well_cells_[perf]; - assert (int(mob.size()) == this->num_components_); - const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - - // either use mobility of the perforation cell or calcualte its own - // based on passing the saturation table index - const int satid = this->saturation_table_number_[perf] - 1; - const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = obtain(intQuants.mobility(phaseIdx)); - } - // if (has_solvent) { - // mob[contiSolventEqIdx] = extendEval(intQuants.solventMobility()); - // } - } else { - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - auto relativePerms = relpermArray(); - MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); - - // reset the satnumvalue back to original - materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); - - // compute the mobility - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = obtain(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); - } - } + WellInterface::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger); } @@ -1228,7 +1182,7 @@ namespace Opm std::vector mob(this->num_components_, 0.0); // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration - getMobility(ebos_simulator, perf, mob); + getMobility(ebos_simulator, perf, mob, deferred_logger); const int cell_idx = this->well_cells_[perf]; const auto& int_quantities = ebos_simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); @@ -1586,7 +1540,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob, deferred_logger); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); @@ -1898,7 +1852,7 @@ namespace Opm const int cell_idx = this->well_cells_[perf]; const auto& int_quants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobility(ebosSimulator, perf, mob); + getMobility(ebosSimulator, perf, mob, deferred_logger); const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier(int_quants, cell_idx); const double Tw = this->well_index_[perf] * trans_mult; std::vector cq_s(this->num_components_, 0.0); diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index d9bf534d9..3acd44f2e 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -834,62 +834,11 @@ namespace Opm if constexpr (std::is_same_v) { return getValue(value); } else { - return this->extendEval(value); + return this->extendEval(value); } }; - auto relpermArray = []() - { - if constexpr (std::is_same_v) { - return std::array{}; - } else { - return std::array{}; - } - }; - - const int cell_idx = this->well_cells_[perf]; - assert (int(mob.size()) == this->num_components_); - const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); - const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); - - // either use mobility of the perforation cell or calcualte its own - // based on passing the saturation table index - const int satid = this->saturation_table_number_[perf] - 1; - const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); - if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = obtain(intQuants.mobility(phaseIdx)); - } - if constexpr (has_solvent) { - mob[Indices::contiSolventEqIdx] = obtain(intQuants.solventMobility()); - } - } else { - const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); - auto relativePerms = relpermArray(); - MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); - - // reset the satnumvalue back to original - materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); - - // compute the mobility - for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - - const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); - mob[activeCompIdx] = obtain(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); - } - - // this may not work if viscosity and relperms has been modified? - if constexpr (has_solvent) { - OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger); - } - } + WellInterface::getMobility(ebosSimulator, perf, mob, + obtain, deferred_logger); // modify the water mobility if polymer is present if constexpr (has_polymer) { diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 5ab9a51cb..d3a2fa336 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -329,14 +329,10 @@ public: } protected: - // simulation parameters const ModelParameters& param_; - std::vector connectionRates_; - std::vector< Scalar > B_avg_; - bool changed_to_stopped_this_step_ = false; double wpolymer() const; @@ -394,6 +390,15 @@ protected: Eval getPerfCellPressure(const FluidState& fs) const; + // get the mobility for specific perforation + template + void getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob, + Callback& extendEval, + [[maybe_unused]] DeferredLogger& deferred_logger) const; + + }; } diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 1489abba9..c9bcdbf4f 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1203,4 +1203,68 @@ namespace Opm } } + template + template + void + WellInterface:: + getMobility(const Simulator& ebosSimulator, + const int perf, + std::vector& mob, + Callback& extendEval, + [[maybe_unused]] DeferredLogger& deferred_logger) const + { + auto relpermArray = []() + { + if constexpr (std::is_same_v) { + return std::array{}; + } else { + return std::array{}; + } + }; + const int cell_idx = this->well_cells_[perf]; + assert (int(mob.size()) == this->num_components_); + const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); + const auto& materialLawManager = ebosSimulator.problem().materialLawManager(); + + // either use mobility of the perforation cell or calculate its own + // based on passing the saturation table index + const int satid = this->saturation_table_number_[perf] - 1; + const int satid_elem = materialLawManager->satnumRegionIdx(cell_idx); + if (satid == satid_elem) { // the same saturation number is used. i.e. just use the mobilty from the cell + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = extendEval(intQuants.mobility(phaseIdx)); + } + if constexpr (has_solvent) { + mob[Indices::contiSolventEqIdx] = extendEval(intQuants.solventMobility()); + } + } else { + const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx); + auto relativePerms = relpermArray(); + MaterialLaw::relativePermeabilities(relativePerms, paramsCell, intQuants.fluidState()); + + // reset the satnumvalue back to original + materialLawManager->connectionMaterialLawParams(satid_elem, cell_idx); + + // compute the mobility + for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + + const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + mob[activeCompIdx] = extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); + } + + // this may not work if viscosity and relperms has been modified? + if constexpr (has_solvent) { + OPM_DEFLOG_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent", deferred_logger); + } + } + } + } // namespace Opm