Calculate the segment rates of the well state correctly, i.e. sum when the well is distributed over multiple processes

This commit is contained in:
Lisa Julia Nebel 2024-11-21 10:53:45 +01:00
parent 9f4fdd09fa
commit 655f4011c7
3 changed files with 65 additions and 15 deletions

View File

@ -81,7 +81,8 @@ scaleSegmentRatesWithWellRates(const std::vector<std::vector<int>>& segment_inle
}
std::vector<Scalar> rates;
WellState<Scalar>::calculateSegmentRates(segment_inlets,
WellState<Scalar>::calculateSegmentRates(ws.parallel_info,
segment_inlets,
segment_perforations,
perforation_rates,
num_single_phase, 0, rates);

View File

@ -747,14 +747,15 @@ void WellState<Scalar>::initWellStateMSWell(const std::vector<Well>& wells_ecl,
// TODO: to see if this strategy can benefit StandardWell too
// TODO: it might cause big problem for gas rate control or if there is a gas rate limit
// maybe the best way is to initialize the fractions first then get the rates
for (int perf = 0; perf < n_activeperf; perf++)
for (std::size_t perf = 0; perf < perf_data.size(); perf++)
perf_rates[perf*np + gaspos] *= 100;
}
const auto& perf_rates = perf_data.phase_rates;
std::vector<Scalar> perforation_rates(perf_rates.begin(), perf_rates.end());
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates);
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
@ -765,10 +766,15 @@ void WellState<Scalar>::initWellStateMSWell(const std::vector<Well>& wells_ecl,
auto& segment_pressure = ws.segments.pressure;
segment_pressure[0] = ws.bhp;
const auto& perf_press = perf_data.pressure;
// 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 = segment_perforations[seg][0];
segment_pressure[seg] = perf_press[first_perf];
const int first_perf = ws.parallel_info.get().globalToLocal(segment_perforations[seg][0]);
if (first_perf > -1) { //-1 indicates that the global id is not on this process
segment_pressure[seg] = perf_press[first_perf];
}
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
@ -776,6 +782,21 @@ void WellState<Scalar>::initWellStateMSWell(const std::vector<Well>& wells_ecl,
segment_pressure[seg] = segment_pressure[segment_set.segmentNumberToIndex(outlet_seg)];
}
}
if (ws.parallel_info.get().communication().size() > 1) {
// Communicate the segment_pressure values
std::vector<Scalar> values_to_combine(segment_indices.size(), 0.0);
for (size_t i = 0; i < segment_indices.size(); ++i) {
values_to_combine[i] = segment_pressure[segment_indices[i]];
}
ws.parallel_info.get().communication().sum(values_to_combine.data(), values_to_combine.size());
// Now make segment_pressure equal across all processes
for (size_t i = 0; i < segment_indices.size(); ++i) {
segment_pressure[segment_indices[i]] = values_to_combine[i];
}
}
}
}
}
@ -810,11 +831,12 @@ void WellState<Scalar>::initWellStateMSWell(const std::vector<Well>& wells_ecl,
template<class Scalar>
void WellState<Scalar>::
calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_rates,
const int np, const int segment,
std::vector<Scalar>& segment_rates)
calculateSegmentRatesBeforeSum(const ParallelWellInfo<Scalar>& pw_info,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_rates,
const int np, const int segment,
std::vector<Scalar>& segment_rates)
{
// the rate of the segment equals to the sum of the contribution from the perforations and inlet segment rates.
// the first segment is always the top segment, its rates should be equal to the well rates.
@ -825,17 +847,34 @@ calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets,
}
// contributions from the perforations belong to this segment
for (const int& perf : segment_perforations[segment]) {
for (int p = 0; p < np; ++p) {
segment_rates[np * segment + p] += perforation_rates[np * perf + p];
auto local_perf = pw_info.globalToLocal(perf);
// If local_perf == -1, then the perforation is not on this process.
// The perforation of the other processes are added in calculateSegmentRates.
if (local_perf > -1) {
for (int p = 0; p < np; ++p) {
segment_rates[np * segment + p] += perforation_rates[np * local_perf + p];
}
}
}
for (const int& inlet_seg : segment_inlets[segment]) {
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates);
calculateSegmentRatesBeforeSum(pw_info, segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates);
for (int p = 0; p < np; ++p) {
segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p];
}
}
}
template<class Scalar>
void WellState<Scalar>::
calculateSegmentRates(const ParallelWellInfo<Scalar>& pw_info,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_rates,
const int np, const int segment,
std::vector<Scalar>& segment_rates)
{
calculateSegmentRatesBeforeSum(pw_info, segment_inlets, segment_perforations, perforation_rates, np, segment, segment_rates);
pw_info.communication().sum(segment_rates.data(), segment_rates.size());
}
template<class Scalar>
void WellState<Scalar>::stopWell(int well_index)

View File

@ -151,8 +151,9 @@ public:
void initWellStateMSWell(const std::vector<Well>& wells_ecl,
const WellState* prev_well_state);
static void calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets, const
std::vector<std::vector<int>>& segment_perforations,
static void calculateSegmentRates(const ParallelWellInfo<Scalar>& pw_info,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_rates,
const int np,
const int segment,
@ -409,6 +410,15 @@ private:
Scalar temperature_first_connection,
const std::vector<PerforationData<Scalar>>& well_perf_data,
const SummaryState& summary_state);
static void calculateSegmentRatesBeforeSum(const ParallelWellInfo<Scalar>& pw_info,
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>& segment_perforations,
const std::vector<Scalar>& perforation_rates,
const int np,
const int segment,
std::vector<Scalar>& segment_rates);
};
} // namespace Opm