guards against non-exsisting phases

This commit is contained in:
Tor Harald Sandve 2020-09-10 08:57:27 +02:00
parent a323487b3b
commit 1a7b617074
2 changed files with 37 additions and 14 deletions

View File

@ -60,11 +60,10 @@ namespace Opm
// sometimes, they are similar, while sometimes, they can have very different forms.
// TODO: the following system looks not rather flexible. Looking into all kinds of possibilities
// TODO: gas is always there? how about oil water case?
// Is it gas oil two phase case?
static const bool gasoil = numEq == 2 && (Indices::compositionSwitchIdx >= 0);
static const int GTotal = 0;
static const int WFrac = gasoil? -1000: 1;
//For oilwater system the GFrac is set to 2 (the same as SPres) but it is never used
static const int GFrac = gasoil? 1 : 2;
static const int SPres = numEq;

View File

@ -2069,8 +2069,12 @@ 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);
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
}
// contribution from the outlet segment
const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
@ -2176,8 +2180,12 @@ 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);
duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq);
duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] -= accelerationPressureLoss.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][GFrac] -= accelerationPressureLoss.derivative(GFrac + numEq);
}
}
@ -2555,8 +2563,12 @@ namespace Opm
// and WFrac and GFrac in seg_upwind
resWell_[seg][comp_idx] -= segment_rate.value();
duneD_[seg][seg][comp_idx][GTotal] -= segment_rate.derivative(GTotal + numEq);
duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq);
duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][comp_idx][WFrac] -= segment_rate.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][comp_idx][GFrac] -= segment_rate.derivative(GFrac + numEq);
}
// pressure derivative should be zero
}
}
@ -2572,8 +2584,12 @@ namespace Opm
// and WFrac and GFrac in inlet_upwind
resWell_[seg][comp_idx] += inlet_rate.value();
duneD_[seg][inlet][comp_idx][GTotal] += inlet_rate.derivative(GTotal + numEq);
duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq);
duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][inlet_upwind][comp_idx][WFrac] += inlet_rate.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][inlet_upwind][comp_idx][GFrac] += inlet_rate.derivative(GFrac + numEq);
}
// pressure derivative should be zero
}
}
@ -3194,8 +3210,12 @@ namespace Opm
resWell_[seg][SPres] = pressure_equation.value();
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
}
// contribution from the outlet segment
const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());
@ -3236,8 +3256,12 @@ namespace Opm
resWell_[seg][SPres] = pressure_equation.value();
duneD_[seg][seg][SPres][SPres] += pressure_equation.derivative(SPres + numEq);
duneD_[seg][seg][SPres][GTotal] += pressure_equation.derivative(GTotal + numEq);
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][WFrac] += pressure_equation.derivative(WFrac + numEq);
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
duneD_[seg][seg_upwind][SPres][GFrac] += pressure_equation.derivative(GFrac + numEq);
}
// contribution from the outlet segment
const int outlet_segment_index = segmentNumberToIndex(segmentSet()[seg].outletSegment());