Fix incorrect indices for 1 and 2-phase cases with MSW.

This commit is contained in:
Atgeirr Flø Rasmussen 2021-05-18 15:31:27 +02:00
parent 842e0a53a4
commit 32d0854f14
2 changed files with 35 additions and 20 deletions

View File

@ -66,13 +66,28 @@ namespace Opm
// TODO: we need to have order for the primary variables and also the order for the well equations.
// sometimes, they are similar, while sometimes, they can have very different forms.
static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
// Table showing the primary variable indices, depending on what phases are present:
//
// WOG OG WG WO W/O/G (single phase)
// GTotal 0 0 0 0 0
// WFrac 1 -1000 1 1 -1000
// GFrac 2 1 -1000 -1000 -1000
// Spres 3 2 2 2 1
static constexpr bool has_water = (Indices::waterSaturationIdx >= 0);
static constexpr bool has_gas = (Indices::compositionSwitchIdx >= 0);
static constexpr bool has_oil = (numPhases - has_gas - has_water) > 0;
// In the implementation, one should use has_wfrac_variable
// rather than has_water to check if you should do something
// with the variable at the WFrac location, similar for GFrac.
static constexpr bool has_wfrac_variable = has_water && numPhases > 1;
static constexpr bool has_gfrac_variable = has_gas && has_oil;
static constexpr int GTotal = 0;
static constexpr int WFrac = has_water ? 1: -1000;
static constexpr int GFrac = has_gas ? has_water + 1 : -1000;
static constexpr int SPres = has_gas + has_water + 1;
static constexpr int WFrac = has_wfrac_variable ? 1 : -1000;
static constexpr int GFrac = has_gfrac_variable ? has_wfrac_variable + 1 : -1000;
static constexpr int SPres = has_wfrac_variable + has_gfrac_variable + 1;
// the number of well equations TODO: it should have a more general strategy for it
static const int numWellEq = numPhases + 1;

View File

@ -741,11 +741,11 @@ namespace Opm
primary_variables_[seg][GTotal] = total_seg_rate;
if (std::abs(total_seg_rate) > 0.) {
if (has_water) {
if (has_wfrac_variable) {
const int water_pos = pu.phase_pos[Water];
primary_variables_[seg][WFrac] = scalingFactor(water_pos) * segment_rates[number_of_phases_ * seg_index + water_pos] / total_seg_rate;
}
if (has_gas) {
if (has_gfrac_variable) {
const int gas_pos = pu.phase_pos[Gas];
primary_variables_[seg][GFrac] = scalingFactor(gas_pos) * segment_rates[number_of_phases_ * seg_index + gas_pos] / total_seg_rate;
}
@ -754,7 +754,7 @@ namespace Opm
// only single phase injection handled
auto phase = well.getInjectionProperties().injectorType;
if (has_water) {
if (has_wfrac_variable) {
if (phase == InjectorType::WATER) {
primary_variables_[seg][WFrac] = 1.0;
} else {
@ -762,7 +762,7 @@ namespace Opm
}
}
if (has_gas) {
if (has_gfrac_variable) {
if (phase == InjectorType::GAS) {
primary_variables_[seg][GFrac] = 1.0;
} else {
@ -771,11 +771,11 @@ namespace Opm
}
} else if (this->isProducer()) { // producers
if (has_water) {
if (has_wfrac_variable) {
primary_variables_[seg][WFrac] = 1.0 / number_of_phases_;
}
if (has_gas) {
if (has_gfrac_variable) {
primary_variables_[seg][GFrac] = 1.0 / number_of_phases_;
}
}
@ -911,13 +911,13 @@ namespace Opm
const std::vector<std::array<double, numWellEq> > old_primary_variables = primary_variables_;
for (int seg = 0; seg < numberOfSegments(); ++seg) {
if (has_water) {
if (has_wfrac_variable) {
const int sign = dwells[seg][WFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][WFrac]) * relaxation_factor, dFLimit);
primary_variables_[seg][WFrac] = old_primary_variables[seg][WFrac] - dx_limited;
}
if (has_gas) {
if (has_gfrac_variable) {
const int sign = dwells[seg][GFrac] > 0. ? 1 : -1;
const double dx_limited = sign * std::min(std::abs(dwells[seg][GFrac]) * relaxation_factor, dFLimit);
primary_variables_[seg][GFrac] = old_primary_variables[seg][GFrac] - dx_limited;
@ -1148,21 +1148,21 @@ namespace Opm
volumeFraction(const int seg, const unsigned compIdx) const
{
if (has_water && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
if (has_wfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx)) {
return primary_variables_evaluation_[seg][WFrac];
}
if (has_gas && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
if (has_gfrac_variable && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) {
return primary_variables_evaluation_[seg][GFrac];
}
// Oil fraction
EvalWell oil_fraction = 1.0;
if (has_water) {
if (has_wfrac_variable) {
oil_fraction -= primary_variables_evaluation_[seg][WFrac];
}
if (has_gas) {
if (has_gfrac_variable) {
oil_fraction -= primary_variables_evaluation_[seg][GFrac];
}
/* if (has_solvent) {
@ -1941,10 +1941,10 @@ namespace Opm
const int seg_upwind = upwinding_segments_[seg];
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
if (has_water) {
if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
}
if (has_gas) {
if (has_gfrac_variable) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
}
@ -2052,10 +2052,10 @@ namespace Opm
resWell_[seg][SPres] -= accelerationPressureLoss.value();
duneD_[seg][seg][SPres][SPres] -= accelerationPressureLoss.derivative(SPres + numEq);
duneD_[seg][seg][SPres][GTotal] -= accelerationPressureLoss.derivative(GTotal + numEq);
if (has_water) {
if (has_wfrac_variable) {
duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq);
}
if (has_gas) {
if (has_gfrac_variable) {
duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq);
}
}