OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #4630 from akva2/unify_mobility

Browse Source 
changed: unify getMobility(Eval|Scalar)
This commit is contained in:
Kai Bao 2023-05-16 00:07:35 +02:00 committed by GitHub
commit 809f6d19ef
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 132 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -235,14 +235,11 @@ namespace Opm
DeferredLogger& deferred_logger);
// get the mobility for specific perforation
void getMobilityEval(const Simulator& ebosSimulator,
template<class Value>
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob) const;
// get the mobility for specific perforation
void getMobilityScalar(const Simulator& ebosSimulator,
const int perf,
std::vector<Scalar>& mob) const;
std::vector<Value>& mob,
DeferredLogger& deferred_logger) const;
void computeWellRatesAtBhpLimit(const Simulator& ebosSimulator,
std::vector<double>& well_flux,

118
opm/simulators/wells/MultisegmentWell_impl.hpp
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.);
getMobilityScalar(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);
getMobilityScalar(ebosSimulator, static_cast<int>(subsetPerfID), mob);
getMobility(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
const auto& fs = fluidState(subsetPerfID);
setToZero(connPI);
@ -1066,113 +1066,25 @@ namespace Opm
deferred_logger);
}
template <typename TypeTag>
template<class Value>
void
MultisegmentWell<TypeTag>::
getMobilityEval(const Simulator& ebosSimulator,
getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob) const
std::vector<Value>& 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<Eval,3> 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));
}
}
}
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
getMobilityScalar(const Simulator& ebosSimulator,
const int perf,
std::vector<Scalar>& mob) const
auto obtain = [this](const Eval& value)
{
// 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());
// }
if constexpr (std::is_same_v<Value, Scalar>) {
return getValue(value);
} else {
const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
std::array<Scalar,3> 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;
return this->extendEval(value);
}
const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
mob[activeCompIdx] = relativePerms[phaseIdx] / getValue(intQuants.fluidState().viscosity(phaseIdx));
};
WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob, obtain, deferred_logger);
}
}
}
@ -1270,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
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<EvalWell> mob(this->num_components_, 0.0);
getMobilityEval(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);
@ -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<Scalar> mob(this->num_components_, 0.0);
getMobilityScalar(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);

12
opm/simulators/wells/StandardWell.hpp
View File

@ -333,18 +333,12 @@ namespace Opm
virtual double getRefDensity() const override;
// get the mobility for specific perforation
void getMobilityEval(const Simulator& ebosSimulator,
template<class Value>
void getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob,
std::vector<Value>& mob,
DeferredLogger& deferred_logger) const;
// get the mobility for specific perforation
void getMobilityScalar(const Simulator& ebosSimulator,
const int perf,
std::vector<Scalar>& mob,
DeferredLogger& deferred_logger) const;
void updateWaterMobilityWithPolymer(const Simulator& ebos_simulator,
const int perf,
std::vector<EvalWell>& mob_water,

141
opm/simulators/wells/StandardWell_impl.hpp
View File

@ -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<EvalWell> 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<double>(intQuants, cell_idx);
@ -820,130 +820,25 @@ namespace Opm
template<typename TypeTag>
template<class Value>
void
StandardWell<TypeTag>::
getMobilityEval(const Simulator& ebosSimulator,
getMobility(const Simulator& ebosSimulator,
const int perf,
std::vector<EvalWell>& mob,
std::vector<Value>& 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());
}
auto obtain = [this](const Eval& value)
{
if constexpr (std::is_same_v<Value, Scalar>) {
return getValue(value);
} else {
const auto& paramsCell = materialLawManager->connectionMaterialLawParams(satid, cell_idx);
std::array<Eval,3> 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);
}
}
// 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) {
updateWaterMobilityWithPolymer(ebosSimulator, perf, mob, deferred_logger);
}
}
}
template<typename TypeTag>
void
StandardWell<TypeTag>::
getMobilityScalar(const Simulator& ebosSimulator,
const int perf,
std::vector<Scalar>& 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<Eval,3> 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);
}
return this->extendEval(value);
}
};
WellInterface<TypeTag>::getMobility(ebosSimulator, perf, mob,
obtain, deferred_logger);
// modify the water mobility if polymer is present
if constexpr (has_polymer) {
@ -954,6 +849,7 @@ 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) {
if constexpr (std::is_same_v<Value, Scalar>) {
std::vector<EvalWell> 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]);
@ -961,6 +857,9 @@ namespace Opm
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<Scalar> 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<EvalWell> mob(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.0});
getMobilityEval(ebosSimulator, static_cast<int>(subsetPerfID), mob, deferred_logger);
getMobility(ebosSimulator, static_cast<int>(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<Scalar> mob(this->num_components_, 0.);
getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
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;
@ -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<Scalar> mob(this->num_components_, 0.);
getMobilityScalar(ebosSimulator, perf, mob, deferred_logger);
getMobility(ebosSimulator, perf, mob, deferred_logger);
std::vector<Scalar> cq_s(this->num_components_, 0.);
double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
const double Tw = this->well_index_[perf] * trans_mult;

13
opm/simulators/wells/WellInterface.hpp
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;
};
}

64
opm/simulators/wells/WellInterface_impl.hpp
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