diff --git a/opm/autodiff/MultisegmentWell_impl.hpp b/opm/autodiff/MultisegmentWell_impl.hpp index 636ab9ebc..17a44fc2f 100644 --- a/opm/autodiff/MultisegmentWell_impl.hpp +++ b/opm/autodiff/MultisegmentWell_impl.hpp @@ -60,12 +60,25 @@ namespace Opm // since we decide to use the WellSegments from the well parser. we can reuse a lot from it. // for other facilities needed but not available from parser, we need to process them here - // initialize the segment_perforations_ + // initialize the segment_perforations_ and update perforation_segment_depth_diffs_ const WellConnections& completion_set = well_ecl_->getConnections(current_step_); - for (int perf = 0; perf < number_of_perforations_; ++perf) { + // 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 + // the current implementation is a temporary solution for now, it should be corrected from the parser + // side + int i_perf_wells = 0; + perf_depth_.resize(number_of_perforations_, 0.); + for (int perf = 0; perf < completion_set.size(); ++perf) { const Connection& connection = completion_set.get(perf); - const int segment_index = segmentNumberToIndex(connection.segment()); - segment_perforations_[segment_index].push_back(perf); + if (connection.state() == WellCompletion::OPEN) { + const int segment_index = segmentNumberToIndex(connection.segment()); + segment_perforations_[segment_index].push_back(i_perf_wells); + perf_depth_[i_perf_wells] = connection.depth(); + const double segment_depth = segmentSet()[segment_index].depth(); + perforation_segment_depth_diffs_[i_perf_wells] = perf_depth_[i_perf_wells] - segment_depth; + i_perf_wells++; + } } // initialize the segment_inlets_ @@ -80,16 +93,6 @@ namespace Opm } } - // callcuate the depth difference between perforations and their segments - perf_depth_.resize(number_of_perforations_, 0.); - for (int seg = 0; seg < numberOfSegments(); ++seg) { - const double segment_depth = segmentSet()[seg].depth(); - for (const int perf : segment_perforations_[seg]) { - perf_depth_[perf] = completion_set.get(perf).depth(); - perforation_segment_depth_diffs_[perf] = perf_depth_[perf] - segment_depth; - } - } - // 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 for (int seg = 1; seg < numberOfSegments(); ++seg) { diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index ad2d45a86..cb63b7085 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -656,7 +656,8 @@ namespace Opm const WellConnections& completion_set = well_ecl->getConnections(time_step); // number of segment for this single well const int well_nseg = segment_set.size(); - const int nperf = completion_set.size(); + // const int nperf = completion_set.size(); + int n_activeperf = 0; nseg_ += well_nseg; for (auto segID = 0*well_nseg; segID < well_nseg; ++segID) { this->seg_number_.push_back(segment_set[segID].segmentNumber()); @@ -664,10 +665,13 @@ namespace Opm // 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); - for (int perf = 0; perf < nperf; ++perf) { + for (int perf = 0; perf < completion_set.size(); ++perf) { const Connection& connection = completion_set.get(perf); - const int segment_index = segment_set.segmentNumberToIndex(connection.segment()); - segment_perforations[segment_index].push_back(perf); + if (connection.state() == WellCompletion::OPEN) { + const int segment_index = segment_set.segmentNumberToIndex(connection.segment()); + segment_perforations[segment_index].push_back(n_activeperf); + n_activeperf++; + } } std::vector> segment_inlets(well_nseg); @@ -688,7 +692,10 @@ namespace Opm const int np = numPhases(); const int start_perf = wells->well_connpos[w]; const int start_perf_next_well = wells->well_connpos[w + 1]; - assert(nperf == (start_perf_next_well - start_perf)); // make sure the information from wells_ecl consistent with wells + + // make sure the information from wells_ecl consistent with wells + assert(n_activeperf == (start_perf_next_well - start_perf)); + if (pu.phase_used[Gas]) { const int gaspos = pu.phase_pos[Gas]; // scale the phase rates for Gas to avoid too bad initial guess for gas fraction @@ -696,7 +703,7 @@ namespace Opm // 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 < nperf; perf++) { + for (int perf = 0; perf < n_activeperf; perf++) { const int perf_pos = start_perf + perf; perfPhaseRates()[np * perf_pos + gaspos] *= 100.; }