From a5c571e54ec1d5b4bbb9d1cbc8bea62cc1d27ac2 Mon Sep 17 00:00:00 2001 From: Lisa Julia Nebel Date: Mon, 26 Aug 2024 11:03:58 +0200 Subject: [PATCH] First steps --- .../wells/MultisegmentWellEquations.cpp | 11 ++- .../wells/MultisegmentWellEquations.hpp | 15 ++- opm/simulators/wells/MultisegmentWellEval.cpp | 17 +++- opm/simulators/wells/MultisegmentWellEval.hpp | 3 +- .../wells/MultisegmentWellGeneric.cpp | 3 +- .../wells/MultisegmentWellSegments.cpp | 14 +-- .../wells/MultisegmentWell_impl.hpp | 97 ++++++++++++++----- opm/simulators/wells/ParallelWellInfo.cpp | 4 +- opm/simulators/wells/ParallelWellInfo.hpp | 4 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 1 + opm/simulators/wells/WellInterface_impl.hpp | 7 +- opm/simulators/wells/WellState.cpp | 65 ++++++++++--- opm/simulators/wells/WellState.hpp | 15 ++- 13 files changed, 191 insertions(+), 65 deletions(-) diff --git a/opm/simulators/wells/MultisegmentWellEquations.cpp b/opm/simulators/wells/MultisegmentWellEquations.cpp index 596ca127798..f6012836b13 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.cpp +++ b/opm/simulators/wells/MultisegmentWellEquations.cpp @@ -55,9 +55,10 @@ MultisegmentWellEquations(const MultisegmentWellGeneric& well) template void MultisegmentWellEquations:: -init(const int num_cells, - const int numPerfs, - const std::vector& cells, +init(const int num_cells_this_process, //num_cells on this process + [[maybe_unused]] const int num_perfs_this_process, //num_perfs_this_process on this process + const int num_perfs_whole_mswell, + const std::vector& cells, // cells for the perforations of all processes const std::vector>& segment_inlets, const std::vector>& perforations) { @@ -78,8 +79,8 @@ init(const int num_cells, } duneD_.setSize(well_.numberOfSegments(), well_.numberOfSegments(), nnz_d); } - duneB_.setSize(well_.numberOfSegments(), num_cells, numPerfs); - duneC_.setSize(well_.numberOfSegments(), num_cells, numPerfs); + duneB_.setSize(well_.numberOfSegments(), num_cells_this_process, num_perfs_whole_mswell); //means: sgs x num_cells_this_process entries and num_perfs_whole_mswell nonzero entries + duneC_.setSize(well_.numberOfSegments(), num_cells_this_process, num_perfs_whole_mswell); // we need to add the off diagonal ones for (auto row = duneD_.createbegin(), diff --git a/opm/simulators/wells/MultisegmentWellEquations.hpp b/opm/simulators/wells/MultisegmentWellEquations.hpp index f204c069928..eba5eddc9c0 100644 --- a/opm/simulators/wells/MultisegmentWellEquations.hpp +++ b/opm/simulators/wells/MultisegmentWellEquations.hpp @@ -22,6 +22,8 @@ #ifndef OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED #define OPM_MULTISEGMENTWELL_EQUATIONS_HEADER_INCLUDED +#include +#include #include #include #include @@ -70,13 +72,15 @@ class MultisegmentWellEquations MultisegmentWellEquations(const MultisegmentWellGeneric& well); //! \brief Setup sparsity pattern for the matrices. - //! \param num_cells Total number of cells - //! \param numPerfs Number of perforations + //! \param num_cells_this_process Total number of cells on this process + //! \param num_perfs_this_process Number of perforations on this process + //! \param num_perfs_whole_mswell Number of perforations of the whole well //! \param cells Cell indices for perforations //! \param segment_inlets Cell indices for segment inlets //! \param segment_perforations Cell indices for segment perforations - void init(const int num_cells, - const int numPerfs, + void init(const int num_cells_this_process, + const int num_perfs_this_process, + const int num_perfs_whole_mswell, const std::vector& cells, const std::vector>& segment_inlets, const std::vector>& segment_perforations); @@ -122,6 +126,9 @@ class MultisegmentWellEquations const int seg_pressure_var_ind, const WellState& well_state) const; + //! \brief Sum with off-process contribution. + void sumDistributed(Parallel::Communication comm); + //! \brief Returns a const reference to the residual. const BVectorWell& residual() const { diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp index 367a68a3ec9..2e7f9008146 100644 --- a/opm/simulators/wells/MultisegmentWellEval.cpp +++ b/opm/simulators/wells/MultisegmentWellEval.cpp @@ -53,8 +53,9 @@ namespace Opm template MultisegmentWellEval:: -MultisegmentWellEval(WellInterfaceIndices& baseif) +MultisegmentWellEval(WellInterfaceIndices& baseif, const ParallelWellInfo& pw_info) : MultisegmentWellGeneric(baseif) + , pw_info_(pw_info) , baseif_(baseif) , linSys_(*this) , primary_variables_(baseif) @@ -69,8 +70,17 @@ void MultisegmentWellEval:: initMatrixAndVectors(const int num_cells) { - linSys_.init(num_cells, baseif_.numPerfs(), - baseif_.cells(), segments_.inlets(), + int num_perfs_this_process = baseif_.numPerfs(); + int num_perfs_whole_mswell = this->pw_info_.communication().sum(num_perfs_this_process); + + std::vector cells_for_perfs_whole_well(num_perfs_whole_mswell, 0.0); + for (int perf = 0; perf < num_perfs_this_process; perf++) { + cells_for_perfs_whole_well[this->pw_info_.getGlobalPerfContainerFactory().localToGlobal(perf)] = baseif_.cells()[perf]; + } + this->pw_info_.communication().sum(cells_for_perfs_whole_well.data(), num_perfs_whole_mswell); + + linSys_.init(num_cells, num_perfs_this_process, num_perfs_whole_mswell, + cells_for_perfs_whole_well, segments_.inlets(), segments_.perforations()); primary_variables_.resize(this->numberOfSegments()); } @@ -433,6 +443,7 @@ MultisegmentWellEval:: getFiniteWellResiduals(const std::vector& B_avg, DeferredLogger& deferred_logger) const { + std::cout << "B_avg.size() = " << B_avg.size() << " and baseif_.numComponents() = " << baseif_.numComponents() << std::endl; assert(int(B_avg.size() ) == baseif_.numComponents()); std::vector residuals(numWellEq + 1, 0.0); diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp index 8f831e06216..e238eda59e1 100644 --- a/opm/simulators/wells/MultisegmentWellEval.hpp +++ b/opm/simulators/wells/MultisegmentWellEval.hpp @@ -67,9 +67,10 @@ class MultisegmentWellEval : public MultisegmentWellGeneric& pw_info_; protected: - MultisegmentWellEval(WellInterfaceIndices& baseif); + MultisegmentWellEval(WellInterfaceIndices& baseif, const ParallelWellInfo& pw_info); void initMatrixAndVectors(const int num_cells); diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp index 1fdc1523dd9..ea94d1b3899 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.cpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -81,7 +81,8 @@ scaleSegmentRatesWithWellRates(const std::vector>& segment_inle } std::vector rates; - WellState::calculateSegmentRates(segment_inlets, + WellState::calculateSegmentRates(ws, + segment_inlets, segment_perforations, perforation_rates, num_single_phase, 0, rates); diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 5f1ba95fa48..21893bf78ad 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -62,7 +62,7 @@ MultisegmentWellSegments:: MultisegmentWellSegments(const int numSegments, WellInterfaceGeneric& well) : perforations_(numSegments) - , perforation_depth_diffs_(well.numPerfs(), 0.0) + , perforation_depth_diffs_(well.wellEcl().getConnections().size(), 0.0) , inlets_(well.wellEcl().getSegments().size()) , depth_diffs_(numSegments, 0.0) , densities_(numSegments, 0.0) @@ -84,8 +84,8 @@ MultisegmentWellSegments(const int numSegments, // well_ecl_ and wells struct different // the current implementation is a temporary solution for now, it should be corrected from the parser // side - int i_perf_wells = 0; - well.perfDepth().resize(well_.numPerfs(), 0.); + // The perfDepth vector contains all perforations across all processes of this well! + well.perfDepth().resize(completion_set.size(), 0.); const auto& segment_set = well_.wellEcl().getSegments(); for (std::size_t perf = 0; perf < completion_set.size(); ++perf) { const Connection& connection = completion_set.get(perf); @@ -98,11 +98,10 @@ MultisegmentWellSegments(const int numSegments, connection.getI() + 1, connection.getJ() + 1, connection.getK() + 1)); } - perforations_[segment_index].push_back(i_perf_wells); - well.perfDepth()[i_perf_wells] = connection.depth(); + perforations_[segment_index].push_back(perf); + well.perfDepth()[perf] = connection.depth(); const Scalar segment_depth = segment_set[segment_index].depth(); - perforation_depth_diffs_[i_perf_wells] = well_.perfDepth()[i_perf_wells] - segment_depth; - i_perf_wells++; + perforation_depth_diffs_[perf] = well_.perfDepth()[perf] - segment_depth; } } @@ -345,6 +344,7 @@ MultisegmentWellSegments:: getPressureDiffSegPerf(const int seg, const int perf) const { + // Attention, perf here has a global numbering! return well_.gravity() * densities_[seg].value() * perforation_depth_diffs_[perf]; } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 0eb5457ea68..007a2582657 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -67,7 +67,7 @@ namespace Opm const int index_of_well, const std::vector>& perf_data) : Base(well, pw_info, time_step, param, rate_converter, pvtRegionIdx, num_components, num_phases, index_of_well, perf_data) - , MSWEval(static_cast&>(*this)) + , MSWEval(static_cast&>(*this), pw_info) , regularize_(false) , segment_fluid_initial_(this->numberOfSegments(), std::vector(this->num_components_, 0.0)) { @@ -138,6 +138,8 @@ namespace Opm // calculate the depth difference between the perforations and the perforated grid block for (int perf = 0; perf < this->number_of_perforations_; ++perf) { + // The variable perf loops over the number_of_perforations_ of *this* process. + // This is why we do not need to convert perf to a local id, it is already *local*. const int cell_idx = this->well_cells_[perf]; this->cell_perforation_depth_diffs_[perf] = depth_arg[cell_idx] - this->perf_depth_[perf]; } @@ -359,14 +361,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_.getGlobalPerfContainerFactory().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 +785,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); @@ -799,6 +808,12 @@ namespace Opm connPI += np; } + // Sum with communication in case of distributed well. + const auto& comm = this->parallel_well_info_.communication(); + if (comm.size() > 1) { + comm.sum(wellPI, np); + } + assert (static_cast(subsetPerfID) == this->number_of_perforations_ && "Internal logic error in processing connections for PI/II"); } @@ -878,11 +893,15 @@ namespace Opm PerforationRates& perf_rates, DeferredLogger& deferred_logger) const { + const int local_perf_index = this->pw_info_.getGlobalPerfContainerFactory().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 +1133,7 @@ namespace Opm void MultisegmentWell:: getMobility(const Simulator& simulator, - const int perf, + const int perf, // this must be a *local* perf index! std::vector& mob, DeferredLogger& deferred_logger) const { @@ -1244,18 +1263,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_.getGlobalPerfContainerFactory().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 +1303,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) { @@ -1316,6 +1338,8 @@ namespace Opm } } } + this->parallel_well_info_.communication().sum(this->ipr_a_.data(), this->ipr_a_.size()); + this->parallel_well_info_.communication().sum(this->ipr_b_.data(), this->ipr_b_.size()); } template @@ -1871,13 +1895,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_.getGlobalPerfContainerFactory().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,21 +1915,21 @@ 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, cell_idx, comp_idx, cq_s_effective, this->linSys_); @@ -1959,15 +1986,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_.getGlobalPerfContainerFactory().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; @@ -1985,6 +2015,12 @@ namespace Opm } } } + const auto& comm = this->parallel_well_info_.communication(); + if (comm.size() > 1) + { + all_drawdown_wrong_direction = + (comm.min(all_drawdown_wrong_direction ? 1 : 0) == 1); + } return all_drawdown_wrong_direction; } @@ -2161,7 +2197,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_.getGlobalPerfContainerFactory().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(); @@ -2189,13 +2229,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_.getGlobalPerfContainerFactory().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; @@ -2206,6 +2250,11 @@ namespace Opm } } } + const auto& comm = this->parallel_well_info_.communication(); + if (comm.size() > 1) + { + comm.sum(well_q_s.data(), well_q_s.size()); + } return well_q_s; } diff --git a/opm/simulators/wells/ParallelWellInfo.cpp b/opm/simulators/wells/ParallelWellInfo.cpp index 856d51ddb53..b3f48250e92 100644 --- a/opm/simulators/wells/ParallelWellInfo.cpp +++ b/opm/simulators/wells/ParallelWellInfo.cpp @@ -120,7 +120,7 @@ GlobalPerfContainerFactory(const IndexSet& local_indices, } template -int GlobalPerfContainerFactory::localToGlobal(int localIndex) const { +int GlobalPerfContainerFactory::localToGlobal(std::size_t localIndex) const { // todo: check why this does not work: //return local_indices_[localIndex].global(); if (local_indices_.begin() == local_indices_.end()) @@ -132,7 +132,7 @@ int GlobalPerfContainerFactory::localToGlobal(int localIndex) const { } template -int GlobalPerfContainerFactory::globalToLocal(int globalIndex) const { +int GlobalPerfContainerFactory::globalToLocal(const int globalIndex) const { if (local_indices_.begin() == local_indices_.end()) return globalIndex; auto it = std::find_if(local_indices_.begin(), local_indices_.end(), [globalIndex](const auto& index) {return index.global() == globalIndex;}); diff --git a/opm/simulators/wells/ParallelWellInfo.hpp b/opm/simulators/wells/ParallelWellInfo.hpp index 773a80b27b7..86bf691bb17 100644 --- a/opm/simulators/wells/ParallelWellInfo.hpp +++ b/opm/simulators/wells/ParallelWellInfo.hpp @@ -162,8 +162,8 @@ class GlobalPerfContainerFactory std::size_t num_components) const; int numGlobalPerfs() const; - int globalToLocal(int globalIndex) const; - int localToGlobal(int localIndex) const; + int globalToLocal(const int globalIndex) const; + int localToGlobal(std::size_t localIndex) const; private: const IndexSet& local_indices_; diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index c71fb90dd8d..efdb43fb844 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -86,6 +86,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/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index d9e1d7e5ad3..03b4bcd560b 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1589,7 +1589,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 - seemingly wellIndex was called with a gloabl instead of a local perforation index!"); + } auto wi = std::vector (this->num_components_, this->well_index_[perf] * trans_mult); @@ -1780,6 +1782,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 - seemingly getMobility was called with a gloabl 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); diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index 431618d1b84..c15c3430814 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -746,14 +746,15 @@ void WellState::initWellStateMSWell(const std::vector& 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 perforation_rates(perf_rates.begin(), perf_rates.end()); - calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates); + calculateSegmentRates(ws, 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 @@ -764,10 +765,19 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, auto& segment_pressure = ws.segments.pressure; segment_pressure[0] = ws.bhp; const auto& perf_press = perf_data.pressure; + // The indices_to_combine contain the indices that are only availabe on one process, we reduce the segment pressure vector + // later for these indices. + std::vector indices_to_combine; for (int seg = 1; seg < well_nseg; ++seg) { + // todo: this is not the ideal solution yet, but ok + // the problem is: segment_pressure has 25 entries but perf_press only has 5 or 7 (in total 12) + // -> but that is exactly what the array segment_perforations is doing, it gets the perforations for each segment, right? 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().getGlobalPerfContainerFactory().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]; + } + indices_to_combine.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 @@ -775,6 +785,19 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, segment_pressure[seg] = segment_pressure[segment_set.segmentNumberToIndex(outlet_seg)]; } } + + // Commuinicate the segment_pressure values + std::vector values_to_combine(indices_to_combine.size(), 0.0); + + for (size_t i = 0; i < indices_to_combine.size(); ++i) { + values_to_combine[i] = segment_pressure[indices_to_combine[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 < indices_to_combine.size(); ++i) { + segment_pressure[indices_to_combine[i]] = values_to_combine[i]; + } } } } @@ -809,11 +832,12 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, template void WellState:: -calculateSegmentRates(const std::vector>& segment_inlets, - const std::vector>& segment_perforations, - const std::vector& perforation_rates, - const int np, const int segment, - std::vector& segment_rates) +calculateSegmentRatesBeforeSum(const SingleWellState& ws, + const std::vector>& segment_inlets, + const std::vector>& segment_perforations, + const std::vector& perforation_rates, + const int np, const int segment, + std::vector& 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. @@ -824,17 +848,34 @@ calculateSegmentRates(const std::vector>& 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 = ws.parallel_info.get().getGlobalPerfContainerFactory().globalToLocal(perf); + if (local_perf > -1) { + // if local_perf == -1, then the perforation is not on this process, therefore we do not sum that here + // but at the end of the function + 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(ws, 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 +void WellState:: +calculateSegmentRates(const SingleWellState& ws, + const std::vector>& segment_inlets, + const std::vector>& segment_perforations, + const std::vector& perforation_rates, + const int np, const int segment, + std::vector& segment_rates) +{ + calculateSegmentRatesBeforeSum(ws, segment_inlets, segment_perforations, perforation_rates, np, segment, segment_rates); + ws.parallel_info.get().communication().sum(segment_rates.data(), segment_rates.size()); +} template void WellState::stopWell(int well_index) diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 2847ec2e09d..5c719c6c385 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -151,14 +151,14 @@ class WellState void initWellStateMSWell(const std::vector& wells_ecl, const WellState* prev_well_state); - static void calculateSegmentRates(const std::vector>& segment_inlets, const - std::vector>& segment_perforations, + static void calculateSegmentRates(const SingleWellState& ws, + const std::vector>& segment_inlets, + const std::vector>& segment_perforations, const std::vector& perforation_rates, const int np, const int segment, std::vector& segment_rates); - void communicateGroupRates(const Parallel::Communication& comm); void updateGlobalIsGrup(const Parallel::Communication& comm); @@ -404,6 +404,15 @@ class WellState Scalar temperature_first_connection, const std::vector>& well_perf_data, const SummaryState& summary_state); + + static void calculateSegmentRatesBeforeSum(const SingleWellState& ws, + const std::vector>& segment_inlets, + const std::vector>& segment_perforations, + const std::vector& perforation_rates, + const int np, + const int segment, + std::vector& segment_rates); + }; } // namespace Opm