diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 5612a34f66d..c0f0cd99298 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -837,32 +837,14 @@ initWellStateMSWell(const std::vector& wells_ecl, ws.parallel_info.get().communication().sum(perforation_pressures.data(), perforation_pressures.size()); } + // NOTE: with default settings, producer segment rates computed here are later overwritten by + // initializeProducerWellState / scaleSegmentRatesWithWellRates calculateSegmentRates(ws.parallel_info, segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates); - // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment - // if there is no perforation associated with this segment, it uses the pressure from the outlet segment - // which requres the ordering is successful - // Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be - // improved during the solveWellEq process - { - // top segment is always the first one, and its pressure is the well bhp - auto& segment_pressure = ws.segments.pressure; - segment_pressure[0] = ws.bhp; - // The segment_indices contain the indices of the segments, that are only available on one process. - std::vector segment_indices; - for (int seg = 1; seg < well_nseg; ++seg) { - if (!segment_perforations[seg].empty()) { - const int first_perf_global_index = segment_perforations[seg][0]; - segment_pressure[seg] = perforation_pressures[first_perf_global_index]; - segment_indices.push_back(seg); - } else { - // seg_press_.push_back(bhp); // may not be a good decision - // using the outlet segment pressure // it needs the ordering is correct - const int outlet_seg = segment_set[seg].outletSegment(); - segment_pressure[seg] = segment_pressure[segment_set.segmentNumberToIndex(outlet_seg)]; - } - } - } + // Segment pressure initialization + setSegmentPressuresFromPerforations(segment_inlets, segment_perforations, perforation_pressures, ws.bhp, ws.segments.pressure); + // finally, set the well bhp equal to the top segment pressure + ws.bhp = ws.segments.pressure[0]; } } @@ -917,6 +899,74 @@ initWellStateMSWell(const std::vector& wells_ecl, } } +template +void WellState:: +setSegmentPressuresFromPerforations(const std::vector>& segment_inlets, + const std::vector>& segment_perforations, + const std::vector& perforation_pressures, + const Scalar& well_bhp, + std::vector& segment_pressures) +{ + // 1. If a segment has perforations, it gets the pressure from the first perforation + // 2. If a segment has no perforations, its pressure is set equal to the pressure of + // its first inlet segment recursively + // 3. If there still are unset pressures after step 2, the pressure of such segments + // are set equal to the pressure of its outlet segment recursively + const int well_nseg = segment_inlets.size(); + for (int seg = 0; seg < well_nseg; ++seg) { + if (!segment_perforations[seg].empty()) { + const int first_perf_global_index = segment_perforations[seg][0]; + segment_pressures[seg] = perforation_pressures[first_perf_global_index]; + } else { + // set to negative value, fetch from inlets (or outlets) below + segment_pressures[seg] = -1.0; + } + } + // Populate unset pressures upward from inlet segments + setSegmentPressuresFromInlets(segment_pressures, segment_inlets, 0 /* segment */); + // If there still are unset segment pressures, it means that the well has "dead ends", + // i.e., some segments have no perforations in inlet direction. In this case, just set these + // pressures equal to the closest set pressure in the outlet direction. + if (std::any_of(segment_pressures.begin(), segment_pressures.end(), + [](const Scalar p){ return p < 0.0; })) { + // If top segment pressure is not set at this point, it means the well has no perforations. + // This shouldn't happen, but in case it does, just set the top segment pressure to the well bhp. + if (segment_pressures[0] < 0.0) { + segment_pressures[0] = well_bhp; + } + setSegmentPressuresFromOutlets(segment_pressures, segment_inlets, 0 /* segment */, 0 /* outlet segment */); + } +} + +template +void WellState:: +setSegmentPressuresFromInlets(std::vector& segment_pressures, + const std::vector>& segment_inlets, + const int segment) +{ + for (const int& inlet_seg : segment_inlets[segment]) { + setSegmentPressuresFromInlets(segment_pressures, segment_inlets, inlet_seg); + if (segment_pressures[segment] < 0.0) { + segment_pressures[segment] = segment_pressures[inlet_seg]; + } + } +} + +template +void WellState:: +setSegmentPressuresFromOutlets(std::vector& segment_pressures, + const std::vector>& segment_inlets, + const int segment, + const int outlet_segment) +{ + if (segment_pressures[segment] < 0.0) { + segment_pressures[segment] = segment_pressures[outlet_segment]; + } + for (const int& inlet_seg : segment_inlets[segment]) { + setSegmentPressuresFromOutlets(segment_pressures, segment_inlets, inlet_seg, segment); + } +} + template void WellState:: calculateSegmentRatesBeforeSum(const ParallelWellInfo& pw_info, diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 3b309ec08a6..0f85767ceee 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -168,6 +168,12 @@ class WellState void initWellStateMSWell(const std::vector& wells_ecl, const WellState* prev_well_state); + void setSegmentPressuresFromPerforations(const std::vector>& segment_inlets, + const std::vector>& segment_perforations, + const std::vector& perforation_pressures, + const Scalar& well_bhp, + std::vector& segment_pressures); + static void calculateSegmentRates(const ParallelWellInfo& pw_info, const std::vector>& segment_inlets, const std::vector>& segment_perforations, @@ -176,7 +182,6 @@ class WellState const int segment, std::vector& segment_rates); - void communicateGroupRates(const Parallel::Communication& comm); void updateGlobalIsGrup(const Parallel::Communication& comm, @@ -409,6 +414,15 @@ class WellState const int segment, std::vector& segment_rates); + void setSegmentPressuresFromInlets(std::vector& segment_pressures, + const std::vector>& segment_inlets, + const int segment); + + void setSegmentPressuresFromOutlets(std::vector& segment_pressures, + const std::vector>& segment_outlets, + const int segment, + const int outlet_segment); + void reportConnectionFactors(const std::size_t well_index, std::vector& connections) const;