introduce getMobility in WellInterface

we now share implementation between StandardWell and MultisegmentWell
This commit is contained in:
Arne Morten Kvarving 2023-05-12 10:02:05 +02:00
parent e406d2f0a1
commit 0406a033a2
5 changed files with 87 additions and 114 deletions

View File

@ -238,7 +238,8 @@ namespace Opm
template<class Value>
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<Value>& mob) const;
std::vector<Value>& mob,
DeferredLogger& deferred_logger) const;
void computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
std::vector<double>& well_flux,

View File

@ -389,7 +389,7 @@ namespace Opm
const auto& intQuants = ebosSimulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
// flux for each perforation
std::vector<Scalar> mob(this->num_components_, 0.);
getMobility(ebosSimulator, perf, mob);
getMobility(ebosSimulator, perf, mob, deferred_logger);
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
@ -716,7 +716,7 @@ namespace Opm
};
std::vector<Scalar> mob(this->num_components_, 0.0);
getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob);
getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID);
setToZero(connPI);
@ -1072,64 +1072,18 @@ namespace Opm
MultisegmentWell<TypeTag>::
getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<Value>& mob) const
std::vector<Value>& mob,
DeferredLogger& deferred_logger) const
{
auto obtain = [this](const Eval& value)
{
if constexpr (std::is_same_v<Value, Scalar>) {
return getValue(value);
} else {
return this->extendEval(value);
return this->extendEval(value);
}
};
auto relpermArray = []()
{
if constexpr (std::is_same_v<Value, Scalar>) {
return std::array<Scalar,3>{};
} else {
return std::array<Eval,3>{};
}
};
// 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<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
}
@ -1228,7 +1182,7 @@ namespace Opm
std::vector<Scalar> 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<EvalWell> mob(this->num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
getMobility(ebosSimulator, perf, mob, deferred_logger);
const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
std::vector<EvalWell> 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<Scalar> mob(this->num_components_, 0.0);
getMobility(ebosSimulator, perf, mob);
getMobility(ebosSimulator, perf, mob, deferred_logger);
const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;
std::vector<Scalar> cq_s(this->num_components_, 0.0);

View File

@ -834,62 +834,11 @@ namespace Opm
if constexpr (std::is_same_v<Value, Scalar>) {
return getValue(value);
} else {
return this->extendEval(value);
return this->extendEval(value);
}
};
auto relpermArray = []()
{
if constexpr (std::is_same_v<Value, Scalar>) {
return std::array<Scalar,3>{};
} else {
return std::array<Eval,3>{};
}
};
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<TypeTag>::getMobility(ebosSimulator, perf, mob,
obtain, deferred_logger);
// modify the water mobility if polymer is present
if constexpr (has_polymer) {

View File

@ -329,14 +329,10 @@ public:
}
protected:
// simulation parameters
const ModelParameters& param_;
std::vector<RateVector> 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<class Value, class Callback>
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<Value>& mob,
Callback& extendEval,
[[maybe_unused]] DeferredLogger& deferred_logger) const;
};
}

View File

@ -1203,4 +1203,68 @@ namespace Opm
}
}
template <typename TypeTag>
template<class Value, class Callback>
void
WellInterface<TypeTag>::
getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<Value>& mob,
Callback& extendEval,
[[maybe_unused]] DeferredLogger& deferred_logger) const
{
auto relpermArray = []()
{
if constexpr (std::is_same_v<Value, Scalar>) {
return std::array<Scalar,3>{};
} else {
return std::array<Eval,3>{};
}
};
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