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;