diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index c22b9a77c..50ff0f629 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -235,14 +235,11 @@ 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, + 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 3293c8bc8..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.); - getMobilityScalar(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); - getMobilityScalar(ebosSimulator, static_cast(subsetPerfID), mob); + getMobility(ebosSimulator, static_cast(subsetPerfID), mob, deferred_logger); const auto& fs = fluidState(subsetPerfID); setToZero(connPI); @@ -1066,115 +1066,27 @@ 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, + DeferredLogger& deferred_logger) 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] = this->extendEval(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] = this->extendEval(relativePerms[phaseIdx] / intQuants.fluidState().viscosity(phaseIdx)); - } - } + auto obtain = [this](const Eval& value) + { + if constexpr (std::is_same_v) { + return getValue(value); + } else { + return this->extendEval(value); + } + }; + WellInterface::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger); } - 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 +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 - getMobilityScalar(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); @@ -1628,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); - getMobilityEval(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); @@ -1940,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); - getMobilityScalar(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.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..3acd44f2e 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,61 +820,25 @@ 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 { - 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] = this->extendEval(intQuants.mobility(phaseIdx)); - } - if (has_solvent) { - mob[Indices::contiSolventEqIdx] = this->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] = this->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); - } - } + auto obtain = [this](const Eval& value) + { + if constexpr (std::is_same_v) { + return getValue(value); + } else { + return this->extendEval(value); + } + }; + WellInterface::getMobility(ebosSimulator, perf, mob, + obtain, deferred_logger); // modify the water mobility if polymer is present if constexpr (has_polymer) { @@ -885,81 +849,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 +942,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 +1349,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 +1551,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 +2408,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; 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