21 #ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED 22 #define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED 24 #include <opm/autodiff/BlackoilModelEnums.hpp> 25 #include <opm/core/wells.h> 26 #include <opm/core/well_controls.h> 27 #include <opm/core/simulator/WellState.hpp> 28 #include <opm/core/props/BlackoilPhases.hpp> 29 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp> 30 #include <opm/common/ErrorMacros.hpp> 47 typedef WellState BaseType;
49 typedef BaseType :: WellMapType WellMapType;
51 using BaseType :: wellRates;
52 using BaseType :: bhp;
53 using BaseType :: perfPress;
54 using BaseType :: wellMap;
55 using BaseType :: numWells;
56 using BaseType :: numPhases;
58 template <
class State,
class PrevWellState>
59 void init(
const Wells* wells,
const State& state,
const PrevWellState& prevState,
const PhaseUsage& pu)
61 init(wells, state.pressure(), prevState, pu);
67 template <
class PrevWellState>
68 void init(
const Wells* wells,
const std::vector<double>& cellPressures ,
const PrevWellState& prevState,
const PhaseUsage& pu)
71 BaseType :: init(wells, cellPressures);
78 const int nw = wells->number_of_wells;
79 if( nw == 0 ) return ;
82 const int np = wells->number_of_phases;
83 const int nperf = wells->well_connpos[nw];
85 perfphaserates_.clear();
86 perfphaserates_.resize(nperf * np, 0.0);
87 for (
int w = 0; w < nw; ++w) {
88 assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
89 const WellControls* ctrl = wells->ctrls[w];
91 if (well_controls_well_is_stopped(ctrl)) {
94 const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
97 for (
int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
98 for (
int p = 0; p < np; ++p) {
99 perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
101 perfPress()[perf] = cellPressures[wells->well_cells[perf]];
109 current_controls_.resize(nw);
110 for (
int w = 0; w < nw; ++w) {
111 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
114 is_new_well_.resize(nw,
true);
116 perfRateSolvent_.clear();
117 perfRateSolvent_.resize(nperf, 0.0);
121 if( ! prevState.wellMap().empty() )
123 typedef typename WellMapType :: const_iterator const_iterator;
124 const_iterator end = prevState.wellMap().end();
125 for (
int w = 0; w < nw; ++w) {
126 std::string name( wells->name[ w ] );
127 const_iterator it = prevState.wellMap().find( name );
131 is_new_well_[w] =
false;
133 const int oldIndex = (*it).second[ 0 ];
134 const int newIndex = w;
137 bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
140 thp()[ newIndex ] = prevState.thp()[ oldIndex ];
143 for(
int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
145 wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
149 const int oldPerf_idx_beg = (*it).second[ 1 ];
150 const int num_perf_old_well = (*it).second[ 2 ];
151 const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
154 if( num_perf_old_well == num_perf_this_well )
156 int old_perf_phase_idx = oldPerf_idx_beg *np;
157 for (
int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
158 perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
160 perfPhaseRates()[ perf_phase_idx ] = prevState.perfPhaseRates()[ old_perf_phase_idx ];
163 for (
int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
164 for (
int p = 0; p < np; ++p) {
165 perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
170 if( num_perf_old_well == num_perf_this_well )
172 int oldPerf_idx = oldPerf_idx_beg;
173 for (
int perf = wells->well_connpos[ newIndex ];
174 perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
176 perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
180 if (pu.has_solvent) {
181 if( num_perf_old_well == num_perf_this_well )
183 int oldPerf_idx = oldPerf_idx_beg;
184 for (
int perf = wells->well_connpos[ newIndex ];
185 perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
187 perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
196 const WellControls* ctrl = wells->ctrls[w];
197 const int nwc = well_controls_get_num(ctrl);
199 for (; ctrl_index < nwc; ++ctrl_index) {
200 if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
205 if (ctrl_index == nwc) {
214 top_segment_index_.reserve(nw);
215 for (
int w = 0; w < nw; ++w) {
216 top_segment_index_.push_back(w);
219 segrates_ = wellRates();
223 template <
class State>
224 void resize(
const Wells* wells,
const State& state,
const PhaseUsage& pu) {
226 init(wells, state, dummy_state, pu) ;
231 const std::vector<double>&
perfPhaseRates()
const {
return perfphaserates_; }
235 const std::vector<int>&
currentControls()
const {
return current_controls_; }
237 data::Wells report(
const PhaseUsage &pu)
const override {
238 data::Wells res = WellState::report(pu);
240 const int nw = this->numWells();
241 if( nw == 0 )
return res;
242 const int np = pu.num_phases;
245 using rt = data::Rates::opt;
246 std::vector< rt > phs( np );
247 if( pu.phase_used[BlackoilPhases::Aqua] ) {
248 phs.at( pu.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
251 if( pu.phase_used[BlackoilPhases::Liquid] ) {
252 phs.at( pu.phase_pos[BlackoilPhases::Liquid] ) =
rt::oil;
255 if( pu.phase_used[BlackoilPhases::Vapour] ) {
256 phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) =
rt::gas;
259 if (pu.has_solvent) {
261 for(
int w = 0; w < nw; ++w ) {
262 res.at( wells_->name[ w ]).rates.set( rt::solvent,
solventWellRate(w) );
273 for(
const auto& wt : this->wellMap() ) {
274 const auto w = wt.second[ 0 ];
275 auto& well = res.at( wt.first );
278 int local_comp_index = 0;
279 for(
auto& comp : well.completions ) {
281 + (np * wt.second[ 1 ])
282 + (np * local_comp_index);
285 for(
int i = 0; i < np; ++i ) {
286 comp.rates.set( phs[ i ], *(rates + i) );
289 assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
297 template <
typename PrevWellState>
299 const int time_step,
const PhaseUsage& pu,
const PrevWellState& prev_well_state)
302 const int nw = wells->number_of_wells;
307 top_segment_index_.clear();
308 top_segment_index_.reserve(nw);
310 segpress_.reserve(nw);
312 segrates_.reserve(nw * numPhases());
317 for (
int w = 0; w < nw; ++w) {
318 const int nw_wells_ecl = wells_ecl.size();
319 int index_well_ecl = 0;
320 const std::string well_name(wells->name[w]);
321 for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
322 if (well_name == wells_ecl[index_well_ecl]->name()) {
328 if (index_well_ecl == nw_wells_ecl) {
329 OPM_THROW(std::logic_error,
"Could not find well " << well_name <<
" in wells_ecl ");
332 const Well* well_ecl = wells_ecl[index_well_ecl];
333 top_segment_index_.push_back(nseg_);
334 if ( !well_ecl->isMultiSegment(time_step) ) {
336 segpress_.push_back(bhp()[w]);
337 const int np = numPhases();
338 for (
int p = 0; p < np; ++p) {
339 segrates_.push_back(wellRates()[np * w + p]);
342 const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step);
344 const CompletionSet& completion_set = well_ecl->getCompletions(time_step);
346 const int well_nseg = segment_set.numberSegment();
347 const int nperf = completion_set.size();
351 std::vector<std::vector<int>> segment_perforations(well_nseg);
352 for (
int perf = 0; perf < nperf; ++perf) {
353 const Completion& completion = completion_set.get(perf);
354 const int segment_number = completion.getSegmentNumber();
355 const int segment_index = segment_set.segmentNumberToIndex(segment_number);
356 segment_perforations[segment_index].push_back(perf);
359 std::vector<std::vector<int>> segment_inlets(well_nseg);
360 for (
int seg = 0; seg < well_nseg; ++seg) {
361 const Segment& segment = segment_set[seg];
362 const int segment_number = segment.segmentNumber();
363 const int outlet_segment_number = segment.outletSegment();
364 if (outlet_segment_number > 0) {
365 const int segment_index = segment_set.segmentNumberToIndex(segment_number);
366 const int outlet_segment_index = segment_set.segmentNumberToIndex(outlet_segment_number);
367 segment_inlets[outlet_segment_index].push_back(segment_index);
374 const int np = numPhases();
375 const int start_perf = wells->well_connpos[w];
376 const int start_perf_next_well = wells->well_connpos[w + 1];
377 assert(nperf == (start_perf_next_well - start_perf));
378 if (pu.phase_used[Gas]) {
379 const int gaspos = pu.phase_pos[Gas];
385 for (
int perf = 0; perf < nperf; perf++) {
386 const int perf_pos = start_perf + perf;
391 const std::vector<double> perforation_rates(
perfPhaseRates().begin() + np * start_perf,
393 std::vector<double> segment_rates;
394 calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 , segment_rates);
395 std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
405 segpress_.push_back(bhp()[w]);
406 const int top_segment = top_segment_index_[w];
407 const int start_perf = wells->well_connpos[w];
408 for (
int seg = 1; seg < well_nseg; ++seg) {
409 if ( !segment_perforations[seg].empty() ) {
410 const int first_perf = segment_perforations[seg][0];
411 segpress_.push_back(perfPress()[start_perf + first_perf]);
415 const int outlet_seg = segment_set[seg].outletSegment();
416 segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
422 assert(
int(segpress_.size()) == nseg_);
423 assert(
int(segrates_.size()) == nseg_ * numPhases() );
425 if (!prev_well_state.wellMap().empty()) {
427 const auto& end = prev_well_state.wellMap().end();
428 const int np = numPhases();
429 for (
int w = 0; w < nw; ++w) {
430 const std::string name( wells->name[w] );
431 const auto& it = prev_well_state.wellMap().find( name );
437 const int old_index_well = (*it).second[0];
438 const int new_index_well = w;
439 const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well);
440 const int new_top_segmnet_index = topSegmentIndex(new_index_well);
441 int number_of_segment = 0;
443 if (new_index_well ==
int(top_segment_index_.size()) - 1) {
444 number_of_segment = nseg_ - new_top_segmnet_index;
446 number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index;
449 for (
int i = 0; i < number_of_segment * np; ++i) {
450 segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i];
453 for (
int i = 0; i < number_of_segment; ++i) {
454 segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i];
462 static void calculateSegmentRates(
const std::vector<std::vector<int>>& segment_inlets,
const std::vector<std::vector<int>>&segment_perforations,
463 const std::vector<double>& perforation_rates,
const int np,
const int segment, std::vector<double>& segment_rates)
467 assert(segment_inlets.size() == segment_perforations.size());
468 const int well_nseg = segment_inlets.size();
470 segment_rates.resize(np * well_nseg, 0.0);
473 for (
const int& perf : segment_perforations[segment]) {
474 for (
int p = 0; p < np; ++p) {
475 segment_rates[np * segment + p] += perforation_rates[np * perf + p];
478 for (
const int& inlet_seg : segment_inlets[segment]) {
479 calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates);
480 for (
int p = 0; p < np; ++p) {
481 segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p];
487 bool isNewWell(
const int w)
const {
488 return is_new_well_[w];
492 void setNewWell(
const int w,
const bool is_new_well) {
493 is_new_well_[w] = is_new_well;
499 const std::vector<double>&
perfRateSolvent()
const {
return perfRateSolvent_; }
503 double solvent_well_rate = 0.0;
504 for (
int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
505 solvent_well_rate += perfRateSolvent_[perf];
507 return solvent_well_rate;
510 const std::vector<double>& segRates()
const 515 std::vector<double>& segRates()
520 const std::vector<double>& segPress()
const 525 std::vector<double>& segPress()
530 int numSegment()
const 535 int topSegmentIndex(
const int w)
const 537 assert(w <
int(top_segment_index_.size()) );
539 return top_segment_index_[w];
543 std::vector<double> perfphaserates_;
544 std::vector<int> current_controls_;
545 std::vector<double> perfRateSolvent_;
551 std::vector<bool> is_new_well_;
555 std::vector<double> segrates_;
556 std::vector<double> segpress_;
559 std::vector<int> top_segment_index_;
567 #endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED double solventWellRate(const int w) const
One rate pr well.
Definition: WellStateFullyImplicitBlackoil.hpp:502
bool gas(const PhaseUsage &pu)
Active gas predicate.
Definition: RateConverter.hpp:318
void init(const Wells *wells, const std::vector< double > &cellPressures, const PrevWellState &prevState, const PhaseUsage &pu)
Allocate and initialize if wells is non-null.
Definition: WellStateFullyImplicitBlackoil.hpp:68
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
std::vector< double > & perfRateSolvent()
One rate pr well connection.
Definition: WellStateFullyImplicitBlackoil.hpp:498
bool oil(const PhaseUsage &pu)
Active oil predicate.
Definition: RateConverter.hpp:305
void initWellStateMSWell(const Wells *wells, const std::vector< const Well *> &wells_ecl, const int time_step, const PhaseUsage &pu, const PrevWellState &prev_well_state)
init the MS well related.
Definition: WellStateFullyImplicitBlackoil.hpp:298
std::vector< double > & perfPhaseRates()
One rate per phase and well connection.
Definition: WellStateFullyImplicitBlackoil.hpp:230
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
std::vector< int > & currentControls()
One current control per well.
Definition: WellStateFullyImplicitBlackoil.hpp:234