From dadfe3a6348a9fac92710e166d60c163ce0d5dc3 Mon Sep 17 00:00:00 2001 From: Lisa Julia Nebel Date: Tue, 15 Oct 2024 14:11:34 +0200 Subject: [PATCH] Use local_perf_index instead of perf index where applicable The vectors that contain info about the perforations of a multisegment well are *local* vectors, yet the way we access them by looping over the perforations is global, so when accessing a value we need to get the local perforation index (and with this also check if the perforation actually resides on the own process) --- .../wells/MultisegmentWell_impl.hpp | 91 ++++++++++++------- opm/simulators/wells/WellInterfaceGeneric.cpp | 1 + opm/simulators/wells/WellInterfaceGeneric.hpp | 2 +- opm/simulators/wells/WellInterface_impl.hpp | 7 +- 4 files changed, 68 insertions(+), 33 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index ffb1f5223..8f53efce5 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -136,9 +136,10 @@ namespace Opm this->initMatrixAndVectors(); // calculate the depth difference between the perforations and the perforated grid block - for (int perf = 0; perf < this->number_of_perforations_; ++perf) { - const int cell_idx = this->well_cells_[perf]; - this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf]; + for (int local_perf_index = 0; local_perf_index < this->number_of_perforations_; ++local_perf_index) { + // This variable loops over the number_of_perforations_ of *this* process, hence it is *local*. + const int cell_idx = this->well_cells_[local_perf_index]; + this->cell_perforation_depth_diffs_[local_perf_index] = depth_arg[cell_idx] - this->perf_depth_[local_perf_index]; } } @@ -358,14 +359,17 @@ namespace Opm const auto& segment_pressure = segments_copy.pressure; for (int seg = 0; seg < nseg; ++seg) { for (const int perf : this->segments_.perforations()[seg]) { - const int cell_idx = this->well_cells_[perf]; + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; + const int cell_idx = this->well_cells_[local_perf_index]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); // flux for each perforation std::vector mob(this->num_components_, 0.); - getMobility(simulator, perf, mob, deferred_logger); + getMobility(simulator, local_perf_index, mob, deferred_logger); const Scalar trans_mult = simulator.problem().template wellTransMultiplier(intQuants, cell_idx); const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const std::vector Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(local_perf_index, intQuants, trans_mult, wellstate_nupcol); const Scalar seg_pressure = segment_pressure[seg]; std::vector cq_s(this->num_components_, 0.); Scalar perf_press = 0.0; @@ -780,6 +784,10 @@ namespace Opm }; std::vector mob(this->num_components_, 0.0); + // The subsetPerfID loops over 0 .. this->perf_data_->size(). + // *(this->perf_data_) contains info about the local processes only, + // hence subsetPerfID is a local perf id and we can call getMobility + // directly with that. getMobility(simulator, static_cast(subsetPerfID), mob, deferred_logger); const auto& fs = fluidState(subsetPerfID); @@ -878,11 +886,15 @@ namespace Opm PerforationRates& perf_rates, DeferredLogger& deferred_logger) const { + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + return; + // pressure difference between the segment and the perforation const Value perf_seg_press_diff = this->gravity() * segment_density * this->segments_.perforation_depth_diff(perf); // pressure difference between the perforation and the grid cell - const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; + const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index]; // perforation pressure is the wellbore pressure corrected to perforation depth // (positive sign due to convention in segments_.perforation_depth_diff() ) @@ -1114,7 +1126,7 @@ namespace Opm void MultisegmentWell:: getMobility(const Simulator& simulator, - const int perf, + const int local_perf_index, std::vector& mob, DeferredLogger& deferred_logger) const { @@ -1128,10 +1140,10 @@ namespace Opm } }; - WellInterface::getMobility(simulator, perf, mob, obtain, deferred_logger); + WellInterface::getMobility(simulator, local_perf_index, mob, obtain, deferred_logger); if (this->isInjector() && this->well_ecl_.getInjMultMode() != Well::InjMultMode::NONE) { - const auto perf_ecl_index = this->perforationData()[perf].ecl_index; + const auto perf_ecl_index = this->perforationData()[local_perf_index].ecl_index; const Connection& con = this->well_ecl_.getConnections()[perf_ecl_index]; const int seg = this->segmentNumberToIndex(con.segment()); // from the reference results, it looks like MSW uses segment pressure instead of BHP here @@ -1139,9 +1151,9 @@ namespace Opm // we can change this depending on what we want const Scalar segment_pres = this->primary_variables_.getSegmentPressure(seg).value(); const Scalar perf_seg_press_diff = this->gravity() * this->segments_.density(seg).value() - * this->segments_.perforation_depth_diff(perf); + * this->segments_.perforation_depth_diff(local_perf_index); const Scalar perf_press = segment_pres + perf_seg_press_diff; - const Scalar multiplier = this->getInjMult(perf, segment_pres, perf_press, deferred_logger); + const Scalar multiplier = this->getInjMult(local_perf_index, segment_pres, perf_press, deferred_logger); for (std::size_t i = 0; i < mob.size(); ++i) { mob[i] *= multiplier; } @@ -1244,18 +1256,21 @@ namespace Opm seg_dp); seg_dp[seg] = dp; for (const int perf : this->segments_.perforations()[seg]) { + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; std::vector mob(this->num_components_, 0.0); // TODO: maybe we should store the mobility somewhere, so that we only need to calculate it one per iteration - getMobility(simulator, perf, mob, deferred_logger); + getMobility(simulator, local_perf_index, mob, deferred_logger); - const int cell_idx = this->well_cells_[perf]; + const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quantities = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); const auto& fs = int_quantities.fluidState(); // pressure difference between the segment and the perforation const Scalar perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf); // pressure difference between the perforation and the grid cell - const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; + const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index]; const Scalar pressure_cell = this->getPerfCellPressure(fs).value(); // calculating the b for the connection @@ -1281,7 +1296,7 @@ namespace Opm // the well index associated with the connection const Scalar trans_mult = simulator.problem().template wellTransMultiplier(int_quantities, cell_idx); const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const std::vector tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol); + const std::vector tw_perf = this->wellIndex(local_perf_index, int_quantities, trans_mult, wellstate_nupcol); std::vector ipr_a_perf(this->ipr_a_.size()); std::vector ipr_b_perf(this->ipr_b_.size()); for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { @@ -1871,13 +1886,16 @@ namespace Opm auto& perf_rates = perf_data.phase_rates; auto& perf_press_state = perf_data.pressure; for (const int perf : this->segments_.perforations()[seg]) { - const int cell_idx = this->well_cells_[perf]; + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; + const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobility(simulator, perf, mob, deferred_logger); + getMobility(simulator, local_perf_index, mob, deferred_logger); const Scalar trans_mult = simulator.problem().template wellTransMultiplier(int_quants, cell_idx); const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const std::vector Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol); std::vector cq_s(this->num_components_, 0.0); EvalWell perf_press; PerforationRates perfRates; @@ -1888,24 +1906,24 @@ namespace Opm if (this->isProducer()) { ws.phase_mixing_rates[ws.dissolved_gas] += perfRates.dis_gas; ws.phase_mixing_rates[ws.vaporized_oil] += perfRates.vap_oil; - perf_data.phase_mixing_rates[perf][ws.dissolved_gas] = perfRates.dis_gas; - perf_data.phase_mixing_rates[perf][ws.vaporized_oil] = perfRates.vap_oil; + perf_data.phase_mixing_rates[local_perf_index][ws.dissolved_gas] = perfRates.dis_gas; + perf_data.phase_mixing_rates[local_perf_index][ws.vaporized_oil] = perfRates.vap_oil; } // store the perf pressure and rates for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { - perf_rates[perf*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value(); + perf_rates[local_perf_index*this->number_of_phases_ + this->modelCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value(); } - perf_press_state[perf] = perf_press.value(); + perf_press_state[local_perf_index] = perf_press.value(); for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) { // the cq_s entering mass balance equations need to consider the efficiency factors. const EvalWell cq_s_effective = cq_s[comp_idx] * this->well_efficiency_factor_; - this->connectionRates_[perf][comp_idx] = Base::restrictEval(cq_s_effective); + this->connectionRates_[local_perf_index][comp_idx] = Base::restrictEval(cq_s_effective); MultisegmentWellAssemble(*this). - assemblePerforationEq(seg, perf, comp_idx, cq_s_effective, this->linSys_); + assemblePerforationEq(seg, local_perf_index, comp_idx, cq_s_effective, this->linSys_); } } @@ -1959,15 +1977,18 @@ namespace Opm for (int seg = 0; seg < nseg; ++seg) { const EvalWell segment_pressure = this->primary_variables_.getSegmentPressure(seg); for (const int perf : this->segments_.perforations()[seg]) { + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; - const int cell_idx = this->well_cells_[perf]; + const int cell_idx = this->well_cells_[local_perf_index]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); const auto& fs = intQuants.fluidState(); // pressure difference between the segment and the perforation const EvalWell perf_seg_press_diff = this->segments_.getPressureDiffSegPerf(seg, perf); // pressure difference between the perforation and the grid cell - const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[perf]; + const Scalar cell_perf_press_diff = this->cell_perforation_pressure_diffs_[local_perf_index]; const Scalar pressure_cell = this->getPerfCellPressure(fs).value(); const Scalar perf_press = pressure_cell - cell_perf_press_diff; @@ -2166,7 +2187,11 @@ namespace Opm const int nseg = this->numberOfSegments(); for (int seg = 0; seg < nseg; ++seg) { for (const int perf : this->segments_.perforations()[seg]) { - const int cell_idx = this->well_cells_[perf]; + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; + + const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); const auto& fs = int_quants.fluidState(); Scalar pressure_cell = this->getPerfCellPressure(fs).value(); @@ -2194,13 +2219,17 @@ namespace Opm // calculating the perforation rate for each perforation that belongs to this segment const Scalar seg_pressure = getValue(this->primary_variables_.getSegmentPressure(seg)); for (const int perf : this->segments_.perforations()[seg]) { - const int cell_idx = this->well_cells_[perf]; + const int local_perf_index = this->pw_info_.globalToLocal(perf); + if (local_perf_index < 0) // then the perforation is not on this process + continue; + + const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); std::vector mob(this->num_components_, 0.0); - getMobility(simulator, perf, mob, deferred_logger); + getMobility(simulator, local_perf_index, mob, deferred_logger); const Scalar trans_mult = simulator.problem().template wellTransMultiplier(int_quants, cell_idx); const auto& wellstate_nupcol = simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_); - const std::vector Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol); + const std::vector Tw = this->wellIndex(local_perf_index, int_quants, trans_mult, wellstate_nupcol); std::vector cq_s(this->num_components_, 0.0); Scalar perf_press = 0.0; PerforationRates perf_rates; diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index b5b2fbbc6..8a9748768 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -88,6 +88,7 @@ WellInterfaceGeneric(const Well& well, // We do not want to count SHUT perforations here, so // it would be wrong to use wells.getConnections().size(). + // This is the number_of_perforations_ on this process only! number_of_perforations_ = perf_data.size(); // perforations related diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 632965cce..b1268847b 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -311,7 +311,7 @@ protected: // well index for each perforation std::vector well_index_; - // number of the perforations for this well + // number of the perforations for this well on this process int number_of_perforations_; // depth for each perforation diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index a63c61c97..03d27af72 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1627,7 +1627,9 @@ namespace Opm { // Add a Forchheimer term to the gas phase CTF if the run uses // either of the WDFAC or the WDFACCOR keywords. - + if (static_cast(perf) >= this->well_cells_.size()) { + OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly wellIndex was called with a global instead of a local perforation index!"); + } auto wi = std::vector (this->num_components_, this->well_index_[perf] * trans_mult); @@ -1818,6 +1820,9 @@ namespace Opm return std::array{}; } }; + if (static_cast(perf) >= this->well_cells_.size()) { + OPM_THROW(std::invalid_argument,"The perforation index exceeds the size of the local containers - possibly getMobility was called with a global instead of a local perforation index!"); + } const int cell_idx = this->well_cells_[perf]; assert (int(mob.size()) == this->num_components_); const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0);