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

fixing the perf rates initialization with zero total_mobility

Browse Source
This commit is contained in:
Kai Bao 2023-06-01 22:39:51 +02:00
parent e7abbe28cb
commit 1d7b33aa55
1 changed files with 17 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
opm/simulators/wells/StandardWellConnections.cpp
View File

@ -35,6 +35,7 @@
#include <opm/simulators/wells/WellInterfaceIndices.hpp> #include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <opm/simulators/wells/WellState.hpp> #include <opm/simulators/wells/WellState.hpp>
#include <limits>
#include <numeric> #include <numeric>
#include <sstream> #include <sstream>
@ -475,21 +476,32 @@ computeProperties(const WellState& well_state,
const int cell_idx = well_.cells()[perf]; const int cell_idx = well_.cells()[perf];
const double well_tw_fraction = well_.wellIndex()[perf] / total_tw; const double well_tw_fraction = well_.wellIndex()[perf] / total_tw;
double total_mobility = 0.0; double total_mobility = 0.0;
double total_invB = 0.;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
int ebosPhaseIdx = well_.flowPhaseToEbosPhaseIdx(p); int ebosPhaseIdx = well_.flowPhaseToEbosPhaseIdx(p);
total_mobility += invB(cell_idx, ebosPhaseIdx) * mobility(cell_idx, ebosPhaseIdx); total_mobility += invB(cell_idx, ebosPhaseIdx) * mobility(cell_idx, ebosPhaseIdx);
total_invB += invB(cell_idx, ebosPhaseIdx);
} }
if constexpr (Indices::enableSolvent) { if constexpr (Indices::enableSolvent) {
total_mobility += solventInverseFormationVolumeFactor(cell_idx) * solventMobility(cell_idx); total_mobility += solventInverseFormationVolumeFactor(cell_idx) * solventMobility(cell_idx);
total_invB += solventInverseFormationVolumeFactor(cell_idx);
} }
const bool non_zero_total_mobility = total_mobility > std::numeric_limits<double>::min();
assert(total_invB > std::numeric_limits<double>::min());
// for the perforation having zero mobility for all the phases, we use a small value to generate a small
// perforation rates for those perforations, at the same time, we can use the rates to recover the mixing
// ratios for those perforations.
constexpr double small_value = 1.e-10;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
int ebosPhaseIdx = well_.flowPhaseToEbosPhaseIdx(p); const int ebosPhaseIdx = well_.flowPhaseToEbosPhaseIdx(p);
perfRates[perf * well_.numComponents() + p] = well_tw_fraction * const double factor = non_zero_total_mobility ?
invB(cell_idx, ebosPhaseIdx) * mobility(cell_idx, ebosPhaseIdx) / total_mobility; invB(cell_idx, ebosPhaseIdx) * mobility(cell_idx, ebosPhaseIdx) / total_mobility : invB(cell_idx, ebosPhaseIdx) * small_value / total_invB;
perfRates[perf * well_.numComponents() + p] = well_tw_fraction * factor;
} }
if constexpr (Indices::enableSolvent) { if constexpr (Indices::enableSolvent) {
perfRates[perf * well_.numComponents() + Indices::contiSolventEqIdx] = well_tw_fraction * const double factor = non_zero_total_mobility ?
solventInverseFormationVolumeFactor(cell_idx) * solventMobility(cell_idx) / total_mobility; solventInverseFormationVolumeFactor(cell_idx) * solventMobility(cell_idx) / total_mobility : solventInverseFormationVolumeFactor(cell_idx) * small_value / total_invB;
perfRates[perf * well_.numComponents() + Indices::contiSolventEqIdx] = well_tw_fraction * factor;
} }
} }
} }