diff --git a/opm/simulators/wells/MultisegmentWellGeneric.cpp b/opm/simulators/wells/MultisegmentWellGeneric.cpp index 3d12bbfa91..525304cfeb 100644 --- a/opm/simulators/wells/MultisegmentWellGeneric.cpp +++ b/opm/simulators/wells/MultisegmentWellGeneric.cpp @@ -77,13 +77,17 @@ scaleSegmentRatesWithWellRates(const std::vector>& segment_inle for (int seg = 0; seg < numberOfSegments(); ++seg) { segment_rates[baseif_.numPhases() * seg + phase] *= well_phase_rate / unscaled_top_seg_rate; } - } else { // In this case, we calculated sumTw + } else if (baseif_.numPerfs() > 0) { // Only go here if there are actually active perforations on this process + // In this case, we calculated sumTw // only handling this specific phase constexpr Scalar num_single_phase = 1; - std::vector perforation_rates(num_single_phase * baseif_.numPerfs(), 0.0); + const int global_num_perf_this_well = ws.parallel_info.get().communication().sum(baseif_.numPerfs()); + std::vector perforation_rates(num_single_phase * global_num_perf_this_well, 0.0); + std::cout << "perforation_rates.size() = " << perforation_rates.size() << std::endl; const Scalar perf_phaserate_scaled = ws.surface_rates[phase] / sumTw; - for (int perf = 0; perf < baseif_.numPerfs(); ++perf) { + for (int perf = 0; perf < global_num_perf_this_well; ++perf) { //strangely this loop does not go over the whole array... perforation_rates[perf] = baseif_.wellIndex()[perf] * perf_phaserate_scaled; + std::cout << "perforation_rates[" << perf << "] = " << perforation_rates[perf] << std::endl; } std::vector rates; diff --git a/opm/simulators/wells/MultisegmentWellSegments.cpp b/opm/simulators/wells/MultisegmentWellSegments.cpp index 4babbfa717..b5e8d2e5a7 100644 --- a/opm/simulators/wells/MultisegmentWellSegments.cpp +++ b/opm/simulators/wells/MultisegmentWellSegments.cpp @@ -80,6 +80,9 @@ MultisegmentWellSegments(const int numSegments, // initialize the segment_perforations_ and update perforation_segment_depth_diffs_ const WellConnections& completion_set = well_.wellEcl().getConnections(); + std::cout << "well_.numPerfs() " << well_.numPerfs() << std::endl; + std::cout << "well_.wellEcl().getConnections().size() " << well_.wellEcl().getConnections().size() << std::endl; + std::cout << "parallel_well_info.numLocalPerfs() " << parallel_well_info.numLocalPerfs() << std::endl; // index of the perforation within wells struct // there might be some perforations not active, which causes the number of the perforations in // well_ecl_ and wells struct different @@ -106,10 +109,14 @@ MultisegmentWellSegments(const int numSegments, const Scalar segment_depth = segment_set[segment_index].depth(); int local_perf_index = parallel_well_info.globalToLocal(i_perf_wells); if (local_perf_index > -1) // If local_perf_index == -1, then the perforation is not on this process - local_perforation_depth_diffs_[local_perf_index] = well_.perfDepth()[i_perf_wells] - segment_depth; + if (local_perforation_depth_diffs_.size() > static_cast(local_perf_index)) { + // If local_perforation_depth_diffs_.size() is not large enough, this is a SHUT connection + local_perforation_depth_diffs_[local_perf_index] = well_.perfDepth()[i_perf_wells] - segment_depth; + } i_perf_wells++; } } + std::cout << "First loop of MultisegmentWellSegments constructor done, local_perforation_depth_diffs_.size() = " << local_perforation_depth_diffs_.size() << std::endl; // initialize the segment_inlets_ for (const Segment& segment : segment_set) { @@ -121,6 +128,7 @@ MultisegmentWellSegments(const int numSegments, inlets_[outlet_segment_index].push_back(segment_index); } } + std::cout << "Second loop of MultisegmentWellSegments constructor done, segment_set.size() = " << segment_set.size() << std::endl; // calculating the depth difference between the segment and its oulet_segments // for the top segment, we will make its zero unless we find other purpose to use this value @@ -131,6 +139,7 @@ MultisegmentWellSegments(const int numSegments, const Scalar outlet_depth = outlet_segment.depth(); depth_diffs_[seg] = segment_depth - outlet_depth; } + std::cout << "Third loop of MultisegmentWellSegments constructor done, numSegments = " << numSegments << std::endl; } template diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index 3a5e9a4a9f..c7854b1082 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -350,6 +350,11 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip computeWellRatesWithBhp, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute computeWellRatesWithBhp, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; const int cell_idx = this->well_cells_[local_perf_index]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); // flux for each perforation @@ -887,6 +892,11 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process return; + if (this->cell_perforation_pressure_diffs_.size() <= static_cast(local_perf_index)) { // then the perforation belongs to a shut connection + std::cout << "skip computePerfRate, local_perf_index = " << local_perf_index << ", << cell_perforation_depth_diffs_.size() = " << this->cell_perforation_pressure_diffs_.size() << std::endl; + return; + } + std::cout << "execute computePerfRate, local_perf_index = " << local_perf_index << ", << cell_perforation_depth_diffs_.size() = " << this->cell_perforation_pressure_diffs_.size() << std::endl; // pressure difference between the segment and the perforation const Value perf_seg_press_diff = this->gravity() * segment_density * @@ -1102,25 +1112,37 @@ namespace Opm // perforated cell // TODO: later to investigate how to handle the pvt region int pvt_region_index; - { - // using the first perforated cell, so we look for global index 0 + + Scalar fsTemperature; + using SaltConcType = typename std::decay().saltConcentration())>::type; + SaltConcType fsSaltConcentration; + std::cout << "computeSegmentFluidProperties, this->well_cells_.size() = " << this->well_cells_.size() << std::endl; + if (this->well_cells_.size() > 0) { + // this->well_cells_ is empty if this process does not contain active perforations + // using the pvt region of first perforated cell, so we look for global index 0 + // TODO: it should be a member of the WellInterface, initialized properly const int cell_idx = this->well_cells_[0]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - // The following broadcast calls are neccessary to ensure that processes that do *not* contain - // the first perforation get the correct temperature, saltConcentration and pvt_region_index - auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); - fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature); - temperature.setValue(fsTemperature); - - auto fsSaltConcentration = fs.saltConcentration(); - fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration); - saltConcentration = this->extendEval(fsSaltConcentration); - + fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); + fsSaltConcentration = fs.saltConcentration(); pvt_region_index = fs.pvtRegionIndex(); - pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index); } + std::cout << "will broadcastFirstPerforationValue now" << std::endl; + + // The following broadcast calls are neccessary to ensure that processes that do *not* contain + // the first perforation get the correct temperature, saltConcentration and pvt_region_index + fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature); + temperature.setValue(fsTemperature); + std::cout << "-" << std::endl; + + fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration); + saltConcentration = this->extendEval(fsSaltConcentration); + std::cout << "-" << std::endl; + + pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index); + std::cout << "done" << std::endl; this->segments_.computeFluidProperties(temperature, saltConcentration, @@ -1267,6 +1289,12 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip updateIPR, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute updateIPR, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + 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 @@ -1850,6 +1878,12 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip assembleWellEqWithoutIteration, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute assembleWellEqWithoutIteration, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + 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); @@ -2000,6 +2034,11 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip allDrawDownWrongDirection, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute allDrawDownWrongDirection, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; const int cell_idx = this->well_cells_[local_perf_index]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); @@ -2058,26 +2097,37 @@ namespace Opm EvalWell temperature; EvalWell saltConcentration; int pvt_region_index; - { + + Scalar fsTemperature; + using SaltConcType = typename std::decay().saltConcentration())>::type; + SaltConcType fsSaltConcentration; + + std::cout << "getSegmentSurfaceVolume, this->well_cells_.size() = " << this->well_cells_.size() << std::endl; + if (this->well_cells_.size() > 0) { + // this->well_cells_ is empty if this process does not contain active perforations // using the pvt region of first perforated cell, so we look for global index 0 // TODO: it should be a member of the WellInterface, initialized properly const int cell_idx = this->well_cells_[0]; const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/0); const auto& fs = intQuants.fluidState(); - // The following broadcast calls are neccessary to ensure that processes that do *not* contain - // the first perforation get the correct temperature, saltConcentration and pvt_region_index - auto fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); - fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature); - temperature.setValue(fsTemperature); - - auto fsSaltConcentration = fs.saltConcentration(); - fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration); - saltConcentration = this->extendEval(fsSaltConcentration); - + fsTemperature = fs.temperature(FluidSystem::oilPhaseIdx).value(); + fsSaltConcentration = fs.saltConcentration(); pvt_region_index = fs.pvtRegionIndex(); - pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index); } + std::cout << "will broadcastFirstPerforationValue now" << std::endl; + // The following broadcast calls are neccessary to ensure that processes that do *not* contain + // the first perforation get the correct temperature, saltConcentration and pvt_region_index + fsTemperature = this->pw_info_.broadcastFirstPerforationValue(fsTemperature); + temperature.setValue(fsTemperature); + std::cout << "-" << std::endl; + + fsSaltConcentration = this->pw_info_.broadcastFirstPerforationValue(fsSaltConcentration); + saltConcentration = this->extendEval(fsSaltConcentration); + std::cout << "-" << std::endl; + + pvt_region_index = this->pw_info_.broadcastFirstPerforationValue(pvt_region_index); + std::cout << "done" << std::endl; return this->segments_.getSurfaceVolume(temperature, saltConcentration, @@ -2227,6 +2277,11 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip maxPerfPress, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute maxPerfPress, for seg " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() " << this->well_cells_.size() << std::endl; const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); @@ -2259,6 +2314,11 @@ namespace Opm const int local_perf_index = this->pw_info_.globalToLocal(perf); if (local_perf_index < 0) // then the perforation is not on this process continue; + if (this->well_cells_.size() <= static_cast(local_perf_index)) {// then the perforation belongs to a shut connection + std::cout << "skip computeCurrentWellRates for segment " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() = " << this->well_cells_.size() << std::endl; + continue; + } + std::cout << "execute computeCurrentWellRates for segment " << seg << ", local_perf_index = " << local_perf_index << ", this->well_cells_.size() = " << this->well_cells_.size() << std::endl; const int cell_idx = this->well_cells_[local_perf_index]; const auto& int_quants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0); diff --git a/opm/simulators/wells/ParallelWellInfo.cpp b/opm/simulators/wells/ParallelWellInfo.cpp index 80c8fcd218..a45253d921 100644 --- a/opm/simulators/wells/ParallelWellInfo.cpp +++ b/opm/simulators/wells/ParallelWellInfo.cpp @@ -74,7 +74,7 @@ GlobalPerfContainerFactory:: GlobalPerfContainerFactory(const IndexSet& local_indices, const Parallel::Communication comm, const int num_local_perfs) - : local_indices_(local_indices), comm_(comm) + : local_indices_(local_indices), comm_(comm), num_local_perfs_(num_local_perfs) { if ( comm_.size() > 1 ) { @@ -165,6 +165,12 @@ int GlobalPerfContainerFactory::globalToLocal(const int globalIndex) con return it->second; } +template +int GlobalPerfContainerFactory::numLocalPerfs() const { + return num_local_perfs_; +} + + template std::vector GlobalPerfContainerFactory:: createGlobal(const std::vector& local_perf_container, @@ -539,6 +545,14 @@ int ParallelWellInfo::globalToLocal(const int globalIndex) const { return globalIndex; } +template +int ParallelWellInfo::numLocalPerfs() const { + if(globalPerfCont_) + return globalPerfCont_->numLocalPerfs(); + else // If globalPerfCont_ is not set up, then this is a sequential run and local and the number of indices are the same. + return 0; +} + template void ParallelWellInfo::communicateFirstPerforation(bool hasFirst) { diff --git a/opm/simulators/wells/ParallelWellInfo.hpp b/opm/simulators/wells/ParallelWellInfo.hpp index 04bdeacda2..65b907e817 100644 --- a/opm/simulators/wells/ParallelWellInfo.hpp +++ b/opm/simulators/wells/ParallelWellInfo.hpp @@ -164,6 +164,7 @@ class GlobalPerfContainerFactory int numGlobalPerfs() const; int globalToLocal(const int globalIndex) const; int localToGlobal(std::size_t localIndex) const; + int numLocalPerfs() const; private: void buildLocalToGlobalMap() const; @@ -174,6 +175,7 @@ class GlobalPerfContainerFactory mutable bool g2l_map_built_ = false; const IndexSet& local_indices_; Parallel::Communication comm_; + int num_local_perfs_; int num_global_perfs_; /// \brief sizes for allgatherv std::vector sizes_; @@ -219,6 +221,7 @@ class ParallelWellInfo int globalToLocal(const int globalIndex) const; int localToGlobal(std::size_t localIndex) const; + int numLocalPerfs() const; /// If the well does not have any open connections the member rankWithFirstPerf /// is not initialized, and no broadcast is performed. In this case the argument diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index c4f6f6c0d9..e008438513 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -92,6 +92,8 @@ WellInterfaceGeneric(const Well& well, ref_depth_ = well.getRefDepth(); + std::cout << "In WellInterfaceGeneric: this->well_ecl_.getConnections().size() = " << this->well_ecl_.getConnections().size() << std::endl; + std::cout << "In WellInterfaceGeneric: perf_data.size() = " << perf_data.size() << std::endl; // We do not want to count SHUT perforations here, so // it would be wrong to use wells.getConnections().size(). // This is the number_of_local_perforations_ on this process only! diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index f3d94f7fd1..95e4d3eb2d 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -679,13 +679,21 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, ws.segments = SegmentState{np, segment_set}; const int well_nseg = segment_set.size(); int n_activeperf = 0; + int n_activeperf_local = 0; // we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment // that is why I think we should use a well model to initialize the WellState here std::vector> segment_perforations(well_nseg); + std::cout << "Will now loop over perforations, from 0 to " << completion_set.size() << std::endl; + std::unordered_map active_perf_index_local_to_global = {}; for (std::size_t perf = 0; perf < completion_set.size(); ++perf) { + if (ws.parallel_info.get().globalToLocal(perf) == 0) { + // Reset the counter for the active perforations on this process, because the local perforation numbering is also reset + n_activeperf_local = 0; + } const Connection& connection = completion_set.get(perf); if (connection.state() == Connection::State::OPEN) { + std::cout << "activeperf_global = " << n_activeperf << ", activeperf_local = " << n_activeperf_local << ", perf_global = " << perf << ", perf_local = " << ws.parallel_info.get().globalToLocal(perf) << std::endl; const int segment_index = segment_set.segmentNumberToIndex(connection.segment()); if (segment_index == -1) { OPM_THROW(std::logic_error, @@ -696,13 +704,20 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, } segment_perforations[segment_index].push_back(n_activeperf); + if (ws.parallel_info.get().globalToLocal(perf) > -1) { + active_perf_index_local_to_global.insert({n_activeperf_local, n_activeperf}); + n_activeperf_local++; + } n_activeperf++; + } else { + std::cout << "no activeperf, perf = " << perf << ", local = " << ws.parallel_info.get().globalToLocal(perf) << std::endl; } } if (!this->enableDistributedWells_ && static_cast(ws.perf_data.size()) != n_activeperf) throw std::logic_error("Distributed multi-segment wells cannot be initialized properly yet."); + std::cout << "(should be the same?) n_activeperf: " << n_activeperf << std::endl; std::vector> segment_inlets(well_nseg); for (int seg = 0; seg < well_nseg; ++seg) { @@ -715,10 +730,20 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, segment_inlets[outlet_segment_index].push_back(segment_index); } } + std::cout << "Now: segment_inlets (of size " << well_nseg << "): " << std::endl; + for (size_t i = 0; i < segment_inlets.size(); ++i) { + std::cout << "segment_inlets[" << i << "]: "; + for (int inlet : segment_inlets[i]) { + std::cout << inlet << ", "; + } + std::cout << std::endl; + } auto& perf_data = ws.perf_data; + std::cout << "in the loop for gas" << std::endl; // for the seg_rates_, now it becomes a recursive solution procedure. if (pu.phase_used[Gas]) { + std::cout << "in the loop for gas" << std::endl; auto& perf_rates = perf_data.phase_rates; const int gaspos = pu.phase_pos[Gas]; // scale the phase rates for Gas to avoid too bad initial guess for gas fraction @@ -730,8 +755,62 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, perf_rates[perf*np + gaspos] *= 100; } + std::cout << "perf_data.size() = " << perf_data.size() << std::endl; + std::cout << "np = " << np << std::endl; const auto& perf_rates = perf_data.phase_rates; - std::vector perforation_rates(perf_rates.begin(), perf_rates.end()); + const auto& perf_press = perf_data.pressure; + // The function calculateSegmentRates as well as the loop filling the segment_pressure work + // with *global* containers. Now we create global vectors containing the phase_rates and + // pressures of all processes. + std::cout << "perf_rates.size() = " << perf_rates.size() << std::endl; + std::cout << "perf_press.size() = " << perf_press.size() << std::endl; + size_t number_of_global_perfs = 0; + + if (ws.parallel_info.get().communication().size() > 1) { + number_of_global_perfs = ws.parallel_info.get().communication().sum(perf_data.size()); + } else { + number_of_global_perfs = perf_data.size(); + } + + std::vector perforation_rates(number_of_global_perfs * np, 0.0); + std::vector perforation_pressures(number_of_global_perfs, 0.0); + std::cout << "After communication: perforation_rates.size() = " << perforation_rates.size() << std::endl; + + assert(perf_data.size() == perf_press.size()); + assert(perf_data.size() * np == perf_rates.size()); + for (size_t perf = 0; perf < perf_data.size(); ++perf) { + std::cout << "looping over perf_data.size(), perf = " << perf; + if (active_perf_index_local_to_global.count(perf) > 0) { + const int global_active_perf_index = active_perf_index_local_to_global.at(perf); + std::cout << ", global_active_perf_index = " << global_active_perf_index << std::endl; + perforation_pressures[global_active_perf_index] = perf_press[perf]; + for (int i = 0; i < np; i++) { + perforation_rates[global_active_perf_index * np + i] = perf_rates[perf * np + i]; + } + } else { + std::cout << std::endl; + OPM_THROW(std::logic_error,fmt::format("Error when initializing MS Well state, there is no active perforation index for the local index {}", perf)); + } + } + + std::cout << "Before communication: calculateSegmentRates: perforation_rates.size() " << perforation_rates.size() << ", and the content:" << std::endl; + std::for_each(perforation_rates.begin(), perforation_rates.end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; + std::cout << "Before communication: calculateSegmentRates: perforation_pressures.size() " << perforation_pressures.size() << ", and the content:" << std::endl; + std::for_each(perforation_pressures.begin(), perforation_pressures.end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; + + if (ws.parallel_info.get().communication().size() > 1) { + ws.parallel_info.get().communication().sum(perforation_rates.data(), perforation_rates.size()); + ws.parallel_info.get().communication().sum(perforation_pressures.data(), perforation_pressures.size()); + } + + std::cout << "After communication: calculateSegmentRates: perforation_rates.size() " << perforation_rates.size() << ", and the content:" << std::endl; + std::for_each(perforation_rates.begin(), perforation_rates.end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; + std::cout << "After communication: calculateSegmentRates: perforation_pressures.size() " << perforation_pressures.size() << ", and the content:" << std::endl; + std::for_each(perforation_pressures.begin(), perforation_pressures.end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; calculateSegmentRates(ws.parallel_info, segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates); @@ -744,38 +823,20 @@ void WellState::initWellStateMSWell(const std::vector& wells_ecl, // 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; - const auto& perf_press = perf_data.pressure; + std::cout << "segment_pressure.size() " << segment_pressure.size() << std::endl; + std::cout << "perf_press.size() " << perforation_pressures.size() << ", perforation_pressures:" << std::endl; + std::for_each(perforation_pressures.begin(), perforation_pressures.end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; + // The segment_indices contain the indices of the segments, that are only available on one process. std::vector segment_indices; for (int seg = 1; seg < well_nseg; ++seg) { if (!segment_perforations[seg].empty()) { - 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]; - } else { - segment_pressure[seg] = 0.0; // setting this to 0 here, this will later be filled by the communication below - } + const int first_perf_global_index = segment_perforations[seg][0]; + std::cout << "first perf of segment " << seg << " seg: " << first_perf_global_index << std::endl; + segment_pressure[seg] = perforation_pressures[first_perf_global_index]; segment_indices.push_back(seg); - } - } - if (ws.parallel_info.get().communication().size() > 1) { - // Communicate the segment_pressure values - std::vector 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]; - } - } - // Before addressing the segments with !segment_perforations[seg].empty(), we need to communicate such that the - // vector segment_pressure contains info from all processes - for (int seg = 1; seg < well_nseg; ++ seg) { - if (segment_perforations[seg].empty()) { + } 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(); @@ -823,6 +884,7 @@ calculateSegmentRatesBeforeSum(const ParallelWellInfo& pw_info, const int np, const int segment, std::vector& segment_rates) { + std::cout << "calculateSegmentRatesBeforeSum for the segment " << segment << std::endl; // 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. assert(segment_inlets.size() == segment_perforations.size()); @@ -831,19 +893,30 @@ calculateSegmentRatesBeforeSum(const ParallelWellInfo& pw_info, segment_rates.resize(np * well_nseg, 0.0); } // contributions from the perforations belong to this segment + std::cout << "Will now look at the perforations of segment " << segment << std::endl; for (const int& perf : segment_perforations[segment]) { auto local_perf = pw_info.globalToLocal(perf); + std::cout << "perf = " << perf << ", local_perf = " << local_perf << std::endl; // 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) { + std::cout << "segment_rates.size() = " << segment_rates.size() << ", perforation_rates.size() = " << perforation_rates.size() << std::endl; for (int p = 0; p < np; ++p) { + std::cout << "will access segment_rates[" << np * segment + p << "] and perforation_rates[" << np * local_perf + p << "]" << std::endl; segment_rates[np * segment + p] += perforation_rates[np * local_perf + p]; } } } + std::cout << "Will now loop over the inlet_segs of segment " << segment << std::endl; + std::cout << "The inlet segs are: "; + std::for_each(segment_inlets[segment].begin(), segment_inlets[segment].end(), [](const auto& entry) {std::cout << entry << ", ";}); + std::cout << std::endl; for (const int& inlet_seg : segment_inlets[segment]) { + std::cout << "Now call calculateSegmentRatesBeforeSum for the inlet_seg " << inlet_seg << std::endl; calculateSegmentRatesBeforeSum(pw_info, segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates); + std::cout << "segment_rates.size() = " << segment_rates.size() << std::endl; for (int p = 0; p < np; ++p) { + std::cout << "will access segment_rates[" << np * segment + p << "] and segment_rates[" << np * inlet_seg + p << "]" << std::endl; segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p]; } } @@ -857,6 +930,9 @@ calculateSegmentRates(const ParallelWellInfo& pw_info, const int np, const int segment, std::vector& segment_rates) { + std::cout << "calculateSegmentRates called, segment = " << segment << std::endl; + std::cout << "segment_perforations.size() = " << segment_perforations.size() << std::endl; + std::cout << "(should be the same across all proceses:) segment_rates.size() = " << segment_rates.size() << std::endl; calculateSegmentRatesBeforeSum(pw_info, segment_inlets, segment_perforations, perforation_rates, np, segment, segment_rates); pw_info.communication().sum(segment_rates.data(), segment_rates.size()); }