Skip to content

Rework initialization of segment pressures #6587

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
steink wants to merge 4 commits into
base: master
Choose a base branch
from
Draft

Rework initialization of segment pressures #6587

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
ed4e182
Rework initialization of segment pressures
steink Nov 3, 2025
cd5a082
Updates according to PR comments
steink Nov 6, 2025
91650e6
Merge branch 'OPM:master' into Rework-initialization-of-segement-pres…
steink Nov 6, 2025
88f7cfc
Merge branch 'OPM:master' into Rework-initialization-of-segement-pres…
steink Nov 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
98 changes: 74 additions & 24 deletions opm/simulators/wells/WellState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -837,32 +837,14 @@ initWellStateMSWell(const std::vector<Well>& 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<int> 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];
}
}

Expand Down Expand Up @@ -917,6 +899,74 @@ initWellStateMSWell(const std::vector<Well>& wells_ecl,
}
}

template<typename Scalar, typename IndexTraits>
void WellState<Scalar, IndexTraits>::
setSegmentPressuresFromPerforations(const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_pressures,
const Scalar& well_bhp,
std::vector<Scalar>& 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<typename Scalar, typename IndexTraits>
void WellState<Scalar, IndexTraits>::
setSegmentPressuresFromInlets(std::vector<Scalar>& segment_pressures,
const std::vector<std::vector<int>>& 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<typename Scalar, typename IndexTraits>
void WellState<Scalar, IndexTraits>::
setSegmentPressuresFromOutlets(std::vector<Scalar>& segment_pressures,
const std::vector<std::vector<int>>& 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<typename Scalar, typename IndexTraits>
void WellState<Scalar, IndexTraits>::
calculateSegmentRatesBeforeSum(const ParallelWellInfo<Scalar>& pw_info,
Expand Down
16 changes: 15 additions & 1 deletion opm/simulators/wells/WellState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,12 @@ class WellState
void initWellStateMSWell(const std::vector<Well>& wells_ecl,
const WellState* prev_well_state);

void setSegmentPressuresFromPerforations(const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_pressures,
const Scalar& well_bhp,
std::vector<Scalar>& segment_pressures);

static void calculateSegmentRates(const ParallelWellInfo<Scalar>& pw_info,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
Expand All @@ -176,7 +182,6 @@ class WellState
const int segment,
std::vector<Scalar>& segment_rates);


void communicateGroupRates(const Parallel::Communication& comm);

void updateGlobalIsGrup(const Parallel::Communication& comm,
Expand Down Expand Up @@ -409,6 +414,15 @@ class WellState
const int segment,
std::vector<Scalar>& segment_rates);

void setSegmentPressuresFromInlets(std::vector<Scalar>& segment_pressures,
const std::vector<std::vector<int>>& segment_inlets,
const int segment);

void setSegmentPressuresFromOutlets(std::vector<Scalar>& segment_pressures,
const std::vector<std::vector<int>>& segment_outlets,
const int segment,
const int outlet_segment);

void reportConnectionFactors(const std::size_t well_index,
std::vector<data::Connection>& connections) const;

Expand Down