Merge pull request #3008 from GitPaean/cleaningup_after_aicd_pr

refactoring the pressure assemble equations for ICD
This commit is contained in:
Atgeirr Flø Rasmussen 2021-01-22 08:26:16 +01:00 committed by GitHub
commit 18a8d78f02
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
2 changed files with 55 additions and 117 deletions

View File

@ -424,7 +424,10 @@ namespace Opm
const Well::ProductionControls& prod_controls,
Opm::DeferredLogger& deferred_logger);
void assemblePressureEq(const int seg, WellState& well_state) const;
void assemblePressureEq(const int seg, const UnitSystem& unit_system,
WellState& well_state, DeferredLogger& deferred_logger) const;
void assembleDefaultPressureEq(const int seg, WellState& well_state) const;
// hytrostatic pressure loss
EvalWell getHydroPressureLoss(const int seg) const;
@ -503,17 +506,19 @@ namespace Opm
double maxPerfPress(const Simulator& ebos_simulator) const;
void assembleSICDPressureEq(const int seg, WellState& well_state) const;
// pressure drop for Spiral ICD segment (WSEGSICD)
EvalWell pressureDropSpiralICD(const int seg) const;
void assembleAICDPressureEq(const int seg, WellState& well_state, const UnitSystem& unit_system) const;
// pressure drop for Autonomous ICD segment (WSEGAICD)
EvalWell pressureDropAutoICD(const int seg, const UnitSystem& unit_system) const;
// assemble the pressure equation for sub-critical valve (WSEGVALV)
void assembleValvePressureEq(const int seg, WellState& well_state) const;
// pressure drop for sub-critical valve (WSEGVALV)
EvalWell pressureDropValve(const int seg) const;
// assemble pressure equation for ICD segments
void assembleICDPressureEq(const int seg, const UnitSystem& unit_system,
WellState& well_state, DeferredLogger& deferred_logger) const;
// check whether the well is operable under BHP limit with current reservoir condition
virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger) override;

View File

@ -2213,7 +2213,7 @@ namespace Opm
template <typename TypeTag>
void
MultisegmentWell<TypeTag>::
assemblePressureEq(const int seg, WellState& well_state) const
assembleDefaultPressureEq(const int seg, WellState& well_state) const
{
assert(seg != 0); // not top segment
@ -2996,22 +2996,8 @@ namespace Opm
const Opm::Schedule& schedule = ebosSimulator.vanguard().schedule();
assembleControlEq(well_state, schedule, summaryState, inj_controls, prod_controls, deferred_logger);
} else {
// TODO: maybe the following should go to the function assemblePressureEq()
switch(segmentSet()[seg].segmentType()) {
case Segment::SegmentType::SICD :
assembleSICDPressureEq(seg, well_state);
break;
case Segment::SegmentType::AICD : {
const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
assembleAICDPressureEq(seg, well_state, unit_system);
break;
}
case Segment::SegmentType::VALVE :
assembleValvePressureEq(seg, well_state);
break;
default :
assemblePressureEq(seg, well_state);
}
const UnitSystem& unit_system = ebosSimulator.vanguard().eclState().getDeckUnitSystem();
assemblePressureEq(seg, unit_system, well_state, deferred_logger);
}
well_state.segPressDrop()[seg] = well_state.segPressDropHydroStatic()[seg] +
@ -3022,6 +3008,27 @@ namespace Opm
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
assemblePressureEq(const int seg, const UnitSystem& unit_system,
WellState& well_state, DeferredLogger& deferred_logger) const
{
switch(segmentSet()[seg].segmentType()) {
case Segment::SegmentType::SICD :
case Segment::SegmentType::AICD :
case Segment::SegmentType::VALVE : {
assembleICDPressureEq(seg, unit_system, well_state,deferred_logger);
break;
}
default :
assembleDefaultPressureEq(seg, well_state);
}
}
template<typename TypeTag>
bool
MultisegmentWell<TypeTag>::
@ -3512,13 +3519,11 @@ namespace Opm
// TODO: the three ICD functions probably can be refactored a little bit to reduce some code duplication
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
assembleSICDPressureEq(const int seg, WellState& well_state) const
assembleICDPressureEq(const int seg, const UnitSystem& unit_system,
WellState& well_state, DeferredLogger& deferred_logger) const
{
// TODO: upwinding needs to be taken care of
// top segment can not be a spiral ICD device
@ -3530,9 +3535,24 @@ namespace Opm
EvalWell pressure_equation = getSegmentPressure(seg);
const auto sicd_pressure_drop = pressureDropSpiralICD(seg);
pressure_equation = pressure_equation - sicd_pressure_drop;
well_state.segPressDropFriction()[seg] = sicd_pressure_drop.value();
EvalWell icd_pressure_drop;
switch(segmentSet()[seg].segmentType()) {
case Segment::SegmentType::SICD :
icd_pressure_drop = pressureDropSpiralICD(seg);
break;
case Segment::SegmentType::AICD :
icd_pressure_drop = pressureDropAutoICD(seg, unit_system);
break;
case Segment::SegmentType::VALVE :
icd_pressure_drop = pressureDropValve(seg);
break;
default: {
OPM_DEFLOG_THROW(std::runtime_error, "Segment " + std::to_string(segmentSet()[seg].segmentNumber())
+ " for well " + name() + " is not of ICD type", deferred_logger);
}
}
pressure_equation = pressure_equation - icd_pressure_drop;
well_state.segPressDropFriction()[seg] = icd_pressure_drop.value();
const int seg_upwind = upwinding_segments_[seg];
resWell_[seg][SPres] = pressure_equation.value();
@ -3558,93 +3578,6 @@ namespace Opm
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
assembleAICDPressureEq(const int seg, WellState& well_state, const UnitSystem& unit_system) const
{
// top segment can not be an AICD device
assert(seg != 0);
// the pressure equation is something like
// p_seg - deltaP - p_outlet = 0.
// the major part is how to calculate the deltaP
EvalWell pressure_equation = getSegmentPressure(seg);
const auto aicd_pressure_drop = pressureDropAutoICD(seg, unit_system);
pressure_equation = pressure_equation - aicd_pressure_drop;
well_state.segPressDropFriction()[seg] = aicd_pressure_drop.value();
const int seg_upwind = upwinding_segments_[seg];
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);
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());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
resWell_[seg][SPres] -= outlet_pressure.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
}
}
template<typename TypeTag>
void
MultisegmentWell<TypeTag>::
assembleValvePressureEq(const int seg, WellState& well_state) const
{
// top segment can not be a spiral ICD device
assert(seg != 0);
// const Valve& valve = *segmentSet()[seg].Valve();
// the pressure equation is something like
// p_seg - deltaP - p_outlet = 0.
// the major part is how to calculate the deltaP
EvalWell pressure_equation = getSegmentPressure(seg);
const int seg_upwind = upwinding_segments_[seg];
const auto valve_pressure_drop = pressureDropValve(seg);
pressure_equation = pressure_equation - valve_pressure_drop;
well_state.segPressDropFriction()[seg] = valve_pressure_drop.value();
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);
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());
const EvalWell outlet_pressure = getSegmentPressure(outlet_segment_index);
resWell_[seg][SPres] -= outlet_pressure.value();
for (int pv_idx = 0; pv_idx < numWellEq; ++pv_idx) {
duneD_[seg][outlet_segment_index][SPres][pv_idx] = -outlet_pressure.derivative(pv_idx + numEq);
}
}
template<typename TypeTag>
std::optional<double>