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 1 commit
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
71 changes: 54 additions & 17 deletions opm/simulators/wells/WellState.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -834,29 +834,37 @@ initWellStateMSWell(const std::vector<Well>& wells_ecl,

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
// Segment pressure initialization
// 1. If a segment has perforations, it uses 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
{
// 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) {
auto& segment_pressures = ws.segments.pressure;
for (int seg = 0; 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);
segment_pressures[seg] = perforation_pressures[first_perf_global_index];
} 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)];
// set to negative value, fetch from inlets (or outlets) afterwards
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; })) {
setSegmentPressuresFromOutlets(segment_pressures, segment_inlets, 0 /* segment */, -1 /* outlet segment */);
Copy link
Member

@GitPaean GitPaean Nov 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

to make this code work segment_pressures[0] can not be negative. an safety checking will be needed here.

Copy link
Contributor Author

@steink steink Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, if top segment pressure us unset at this stage, it means the well has no connections which I guess should not happen at this stage, but better to be safe. In this situation arises, I guess it's best just to set it to bhp (which means all segments will get that pressure).

Copy link
Member

@GitPaean GitPaean Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

or more directly, inside the function setSegmentPressuresFromOutlets(), make sure the -1 index should not be used.

}
assert(!std::any_of(segment_pressures.begin(), segment_pressures.end(),
[](const Scalar p){ return p < 0.0; }) &&
"Segment pressures could not be initialized properly and/or reservoir has negative pressures.");
// finally, set the well bhp to be the top segment pressure
ws.bhp = segment_pressures[0];
}
}
}
Expand Down Expand Up @@ -912,6 +920,35 @@ initWellStateMSWell(const std::vector<Well>& wells_ecl,
}
}

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
8 changes: 8 additions & 0 deletions opm/simulators/wells/WellState.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,14 @@ class WellState
const int segment,
std::vector<Scalar>& segment_rates);

void setSegmentPressuresFromInlets(std::vector<Scalar>& segment_pressures,
Copy link
Member

@GitPaean GitPaean Nov 3, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These two functions have to be called together. Calling either one only is not a complete process.

They should be made private, and called through another function (can also be private unless needs to be public) including these two function calls.

Copy link
Contributor Author

@steink steink Nov 4, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will set them private. I'm not 100% sure the second one is needed (that is, if a well is allowed to have toe-segemnts without perforations), but included it just in case.

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 communicateGroupRates(const Parallel::Communication& comm);

Expand Down