00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
00022 #define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
00023
00024 #include <opm/autodiff/BlackoilModelEnums.hpp>
00025 #include <opm/core/wells.h>
00026 #include <opm/core/well_controls.h>
00027 #include <opm/core/simulator/WellState.hpp>
00028 #include <opm/core/props/BlackoilPhases.hpp>
00029 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
00030 #include <opm/common/ErrorMacros.hpp>
00031 #include <vector>
00032 #include <cassert>
00033 #include <string>
00034 #include <utility>
00035 #include <map>
00036 #include <algorithm>
00037 #include <array>
00038
00039 namespace Opm
00040 {
00041
00044 class WellStateFullyImplicitBlackoil
00045 : public WellState
00046 {
00047 typedef WellState BaseType;
00048 public:
00049 typedef BaseType :: WellMapType WellMapType;
00050
00051 using BaseType :: wellRates;
00052 using BaseType :: bhp;
00053 using BaseType :: perfPress;
00054 using BaseType :: wellMap;
00055 using BaseType :: numWells;
00056 using BaseType :: numPhases;
00057
00058 template <class State, class PrevWellState>
00059 void init(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
00060 {
00061 init(wells, state.pressure(), prevState, pu);
00062 }
00063
00067 template <class PrevWellState>
00068 void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
00069 {
00070
00071 BaseType :: init(wells, cellPressures);
00072
00073
00074 if (wells == 0) {
00075 return;
00076 }
00077
00078 const int nw = wells->number_of_wells;
00079 if( nw == 0 ) return ;
00080
00081
00082 const int np = wells->number_of_phases;
00083 const int nperf = wells->well_connpos[nw];
00084
00085 perfphaserates_.clear();
00086 perfphaserates_.resize(nperf * np, 0.0);
00087 for (int w = 0; w < nw; ++w) {
00088 assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
00089 const WellControls* ctrl = wells->ctrls[w];
00090
00091 if (well_controls_well_is_stopped(ctrl)) {
00092
00093 } else {
00094 const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
00095
00096
00097 for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
00098 for (int p = 0; p < np; ++p) {
00099 perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
00100 }
00101 perfPress()[perf] = cellPressures[wells->well_cells[perf]];
00102 }
00103 }
00104 }
00105
00106
00107
00108
00109 current_controls_.resize(nw);
00110 for (int w = 0; w < nw; ++w) {
00111 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
00112 }
00113
00114 is_new_well_.resize(nw, true);
00115
00116 perfRateSolvent_.clear();
00117 perfRateSolvent_.resize(nperf, 0.0);
00118
00119
00120
00121 if( ! prevState.wellMap().empty() )
00122 {
00123 typedef typename WellMapType :: const_iterator const_iterator;
00124 const_iterator end = prevState.wellMap().end();
00125 for (int w = 0; w < nw; ++w) {
00126 std::string name( wells->name[ w ] );
00127 const_iterator it = prevState.wellMap().find( name );
00128 if( it != end )
00129 {
00130
00131 is_new_well_[w] = false;
00132
00133 const int oldIndex = (*it).second[ 0 ];
00134 const int newIndex = w;
00135
00136
00137 bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
00138
00139
00140 thp()[ newIndex ] = prevState.thp()[ oldIndex ];
00141
00142
00143 for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
00144 {
00145 wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
00146 }
00147
00148
00149 const int oldPerf_idx_beg = (*it).second[ 1 ];
00150 const int num_perf_old_well = (*it).second[ 2 ];
00151 const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
00152
00153
00154 if( num_perf_old_well == num_perf_this_well )
00155 {
00156 int old_perf_phase_idx = oldPerf_idx_beg *np;
00157 for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
00158 perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
00159 {
00160 perfPhaseRates()[ perf_phase_idx ] = prevState.perfPhaseRates()[ old_perf_phase_idx ];
00161 }
00162 } else {
00163 for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
00164 for (int p = 0; p < np; ++p) {
00165 perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
00166 }
00167 }
00168 }
00169
00170 if( num_perf_old_well == num_perf_this_well )
00171 {
00172 int oldPerf_idx = oldPerf_idx_beg;
00173 for (int perf = wells->well_connpos[ newIndex ];
00174 perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
00175 {
00176 perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
00177 }
00178 }
00179
00180 if (pu.has_solvent) {
00181 if( num_perf_old_well == num_perf_this_well )
00182 {
00183 int oldPerf_idx = oldPerf_idx_beg;
00184 for (int perf = wells->well_connpos[ newIndex ];
00185 perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
00186 {
00187 perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
00188 }
00189 }
00190 }
00191 }
00192
00193
00194
00195
00196 const WellControls* ctrl = wells->ctrls[w];
00197 const int nwc = well_controls_get_num(ctrl);
00198 int ctrl_index = 0;
00199 for (; ctrl_index < nwc; ++ctrl_index) {
00200 if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
00201 break;
00202 }
00203 }
00204
00205 if (ctrl_index == nwc) {
00206 thp()[w] = 0.;
00207 }
00208 }
00209 }
00210
00211 {
00212
00213
00214 top_segment_index_.reserve(nw);
00215 for (int w = 0; w < nw; ++w) {
00216 top_segment_index_.push_back(w);
00217 }
00218 segpress_ = bhp();
00219 segrates_ = wellRates();
00220 }
00221 }
00222
00223 template <class State>
00224 void resize(const Wells* wells, const State& state, const PhaseUsage& pu) {
00225 const WellStateFullyImplicitBlackoil dummy_state{};
00226 init(wells, state, dummy_state, pu) ;
00227 }
00228
00230 std::vector<double>& perfPhaseRates() { return perfphaserates_; }
00231 const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
00232
00234 std::vector<int>& currentControls() { return current_controls_; }
00235 const std::vector<int>& currentControls() const { return current_controls_; }
00236
00237 data::Wells report(const PhaseUsage &pu) const override {
00238 data::Wells res = WellState::report(pu);
00239
00240 const int nw = this->numWells();
00241 if( nw == 0 ) return res;
00242 const int np = pu.num_phases;
00243
00244
00245 using rt = data::Rates::opt;
00246 std::vector< rt > phs( np );
00247 if( pu.phase_used[BlackoilPhases::Aqua] ) {
00248 phs.at( pu.phase_pos[BlackoilPhases::Aqua] ) = rt::wat;
00249 }
00250
00251 if( pu.phase_used[BlackoilPhases::Liquid] ) {
00252 phs.at( pu.phase_pos[BlackoilPhases::Liquid] ) = rt::oil;
00253 }
00254
00255 if( pu.phase_used[BlackoilPhases::Vapour] ) {
00256 phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
00257 }
00258
00259 if (pu.has_solvent) {
00260
00261 for( int w = 0; w < nw; ++w ) {
00262 res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
00263 }
00264 }
00265
00266
00267
00268
00269
00270
00271
00272
00273 for( const auto& wt : this->wellMap() ) {
00274 const auto w = wt.second[ 0 ];
00275 auto& well = res.at( wt.first );
00276 well.control = this->currentControls()[ w ];
00277
00278 int local_comp_index = 0;
00279 for( auto& comp : well.completions ) {
00280 const auto rates = this->perfPhaseRates().begin()
00281 + (np * wt.second[ 1 ])
00282 + (np * local_comp_index);
00283 ++local_comp_index;
00284
00285 for( int i = 0; i < np; ++i ) {
00286 comp.rates.set( phs[ i ], *(rates + i) );
00287 }
00288 }
00289 assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
00290 }
00291
00292 return res;
00293 }
00294
00295
00297 template <typename PrevWellState>
00298 void initWellStateMSWell(const Wells* wells, const std::vector<const Well*>& wells_ecl,
00299 const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state)
00300 {
00301
00302 const int nw = wells->number_of_wells;
00303 if (nw == 0) {
00304 return;
00305 }
00306
00307 top_segment_index_.clear();
00308 top_segment_index_.reserve(nw);
00309 segpress_.clear();
00310 segpress_.reserve(nw);
00311 segrates_.clear();
00312 segrates_.reserve(nw * numPhases());
00313
00314 nseg_ = 0;
00315
00316
00317 for (int w = 0; w < nw; ++w) {
00318 const int nw_wells_ecl = wells_ecl.size();
00319 int index_well_ecl = 0;
00320 const std::string well_name(wells->name[w]);
00321 for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
00322 if (well_name == wells_ecl[index_well_ecl]->name()) {
00323 break;
00324 }
00325 }
00326
00327
00328 if (index_well_ecl == nw_wells_ecl) {
00329 OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
00330 }
00331
00332 const Well* well_ecl = wells_ecl[index_well_ecl];
00333 top_segment_index_.push_back(nseg_);
00334 if ( !well_ecl->isMultiSegment(time_step) ) {
00335 nseg_ += 1;
00336 segpress_.push_back(bhp()[w]);
00337 const int np = numPhases();
00338 for (int p = 0; p < np; ++p) {
00339 segrates_.push_back(wellRates()[np * w + p]);
00340 }
00341 } else {
00342 const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step);
00343
00344 const CompletionSet& completion_set = well_ecl->getCompletions(time_step);
00345
00346 const int well_nseg = segment_set.numberSegment();
00347 const int nperf = completion_set.size();
00348 nseg_ += well_nseg;
00349
00350
00351 std::vector<std::vector<int>> segment_perforations(well_nseg);
00352 for (int perf = 0; perf < nperf; ++perf) {
00353 const Completion& completion = completion_set.get(perf);
00354 const int segment_number = completion.getSegmentNumber();
00355 const int segment_index = segment_set.segmentNumberToIndex(segment_number);
00356 segment_perforations[segment_index].push_back(perf);
00357 }
00358
00359 std::vector<std::vector<int>> segment_inlets(well_nseg);
00360 for (int seg = 0; seg < well_nseg; ++seg) {
00361 const Segment& segment = segment_set[seg];
00362 const int segment_number = segment.segmentNumber();
00363 const int outlet_segment_number = segment.outletSegment();
00364 if (outlet_segment_number > 0) {
00365 const int segment_index = segment_set.segmentNumberToIndex(segment_number);
00366 const int outlet_segment_index = segment_set.segmentNumberToIndex(outlet_segment_number);
00367 segment_inlets[outlet_segment_index].push_back(segment_index);
00368 }
00369 }
00370
00371
00372
00373 {
00374 const int np = numPhases();
00375 const int start_perf = wells->well_connpos[w];
00376 const int start_perf_next_well = wells->well_connpos[w + 1];
00377 assert(nperf == (start_perf_next_well - start_perf));
00378 if (pu.phase_used[Gas]) {
00379 const int gaspos = pu.phase_pos[Gas];
00380
00381
00382
00383
00384
00385 for (int perf = 0; perf < nperf; perf++) {
00386 const int perf_pos = start_perf + perf;
00387 perfPhaseRates()[np * perf_pos + gaspos] *= 100.;
00388 }
00389 }
00390
00391 const std::vector<double> perforation_rates(perfPhaseRates().begin() + np * start_perf,
00392 perfPhaseRates().begin() + np * start_perf_next_well);
00393 std::vector<double> segment_rates;
00394 calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 , segment_rates);
00395 std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
00396 }
00397
00398
00399
00400
00401
00402
00403 {
00404
00405 segpress_.push_back(bhp()[w]);
00406 const int top_segment = top_segment_index_[w];
00407 const int start_perf = wells->well_connpos[w];
00408 for (int seg = 1; seg < well_nseg; ++seg) {
00409 if ( !segment_perforations[seg].empty() ) {
00410 const int first_perf = segment_perforations[seg][0];
00411 segpress_.push_back(perfPress()[start_perf + first_perf]);
00412 } else {
00413
00414
00415 const int outlet_seg = segment_set[seg].outletSegment();
00416 segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
00417 }
00418 }
00419 }
00420 }
00421 }
00422 assert(int(segpress_.size()) == nseg_);
00423 assert(int(segrates_.size()) == nseg_ * numPhases() );
00424
00425 if (!prev_well_state.wellMap().empty()) {
00426
00427 const auto& end = prev_well_state.wellMap().end();
00428 const int np = numPhases();
00429 for (int w = 0; w < nw; ++w) {
00430 const std::string name( wells->name[w] );
00431 const auto& it = prev_well_state.wellMap().find( name );
00432
00433 if (it != end) {
00434
00435
00436
00437 const int old_index_well = (*it).second[0];
00438 const int new_index_well = w;
00439 const int old_top_segment_index = prev_well_state.topSegmentIndex(old_index_well);
00440 const int new_top_segmnet_index = topSegmentIndex(new_index_well);
00441 int number_of_segment = 0;
00442
00443 if (new_index_well == int(top_segment_index_.size()) - 1) {
00444 number_of_segment = nseg_ - new_top_segmnet_index;
00445 } else {
00446 number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index;
00447 }
00448
00449 for (int i = 0; i < number_of_segment * np; ++i) {
00450 segrates_[new_top_segmnet_index * np + i] = prev_well_state.segRates()[old_top_segment_index * np + i];
00451 }
00452
00453 for (int i = 0; i < number_of_segment; ++i) {
00454 segpress_[new_top_segmnet_index + i] = prev_well_state.segPress()[old_top_segment_index + i];
00455 }
00456 }
00457 }
00458 }
00459 }
00460
00461
00462 static void calculateSegmentRates(const std::vector<std::vector<int>>& segment_inlets, const std::vector<std::vector<int>>&segment_perforations,
00463 const std::vector<double>& perforation_rates, const int np, const int segment, std::vector<double>& segment_rates)
00464 {
00465
00466
00467 assert(segment_inlets.size() == segment_perforations.size());
00468 const int well_nseg = segment_inlets.size();
00469 if (segment == 0) {
00470 segment_rates.resize(np * well_nseg, 0.0);
00471 }
00472
00473 for (const int& perf : segment_perforations[segment]) {
00474 for (int p = 0; p < np; ++p) {
00475 segment_rates[np * segment + p] += perforation_rates[np * perf + p];
00476 }
00477 }
00478 for (const int& inlet_seg : segment_inlets[segment]) {
00479 calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, inlet_seg, segment_rates);
00480 for (int p = 0; p < np; ++p) {
00481 segment_rates[np * segment + p] += segment_rates[np * inlet_seg + p];
00482 }
00483 }
00484 }
00485
00486
00487 bool isNewWell(const int w) const {
00488 return is_new_well_[w];
00489 }
00490
00491
00492 void setNewWell(const int w, const bool is_new_well) {
00493 is_new_well_[w] = is_new_well;
00494 }
00495
00496
00498 std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
00499 const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
00500
00502 double solventWellRate(const int w) const {
00503 double solvent_well_rate = 0.0;
00504 for (int perf = wells_->well_connpos[w]; perf < wells_->well_connpos[w+1]; ++perf ) {
00505 solvent_well_rate += perfRateSolvent_[perf];
00506 }
00507 return solvent_well_rate;
00508 }
00509
00510 const std::vector<double>& segRates() const
00511 {
00512 return segrates_;
00513 }
00514
00515 std::vector<double>& segRates()
00516 {
00517 return segrates_;
00518 }
00519
00520 const std::vector<double>& segPress() const
00521 {
00522 return segpress_;
00523 }
00524
00525 std::vector<double>& segPress()
00526 {
00527 return segpress_;
00528 }
00529
00530 int numSegment() const
00531 {
00532 return nseg_;
00533 }
00534
00535 int topSegmentIndex(const int w) const
00536 {
00537 assert(w < int(top_segment_index_.size()) );
00538
00539 return top_segment_index_[w];
00540 }
00541
00542 private:
00543 std::vector<double> perfphaserates_;
00544 std::vector<int> current_controls_;
00545 std::vector<double> perfRateSolvent_;
00546
00547
00548
00549
00550
00551 std::vector<bool> is_new_well_;
00552
00553
00554
00555 std::vector<double> segrates_;
00556 std::vector<double> segpress_;
00557
00558
00559 std::vector<int> top_segment_index_;
00560 int nseg_;
00561
00562 };
00563
00564 }
00565
00566
00567 #endif // OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED