WellStateFullyImplicitBlackoil.hpp
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3  Copyright 2017 IRIS AS
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
22 #define OPM_WELLSTATEFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
23 
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>
31 #include <vector>
32 #include <cassert>
33 #include <string>
34 #include <utility>
35 #include <map>
36 #include <algorithm>
37 #include <array>
38 
39 namespace Opm
40 {
41 
45  : public WellState
46  {
47  typedef WellState BaseType;
48  public:
49  typedef BaseType :: WellMapType WellMapType;
50 
51  using BaseType :: wellRates;
52  using BaseType :: bhp;
53  using BaseType :: perfPress;
54  using BaseType :: wellMap;
55  using BaseType :: numWells;
56  using BaseType :: numPhases;
57 
58  template <class State, class PrevWellState>
59  void init(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
60  {
61  init(wells, state.pressure(), prevState, pu);
62  }
63 
67  template <class PrevWellState>
68  void init(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
69  {
70  // call init on base class
71  BaseType :: init(wells, cellPressures);
72 
73  // if there are no well, do nothing in init
74  if (wells == 0) {
75  return;
76  }
77 
78  const int nw = wells->number_of_wells;
79  if( nw == 0 ) return ;
80 
81  // Initialize perfphaserates_, which must be done here.
82  const int np = wells->number_of_phases;
83  const int nperf = wells->well_connpos[nw];
84  // Ensure that we start out with zero rates by default.
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];
90 
91  if (well_controls_well_is_stopped(ctrl)) {
92  // Shut well: perfphaserates_ are all zero.
93  } else {
94  const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
95  // Open well: Initialize perfphaserates_ to well
96  // rates divided by the number of perforations.
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);
100  }
101  perfPress()[perf] = cellPressures[wells->well_cells[perf]];
102  }
103  }
104  }
105 
106  // Initialize current_controls_.
107  // The controls set in the Wells object are treated as defaults,
108  // and also used for initial values.
109  current_controls_.resize(nw);
110  for (int w = 0; w < nw; ++w) {
111  current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
112  }
113 
114  is_new_well_.resize(nw, true);
115 
116  perfRateSolvent_.clear();
117  perfRateSolvent_.resize(nperf, 0.0);
118 
119  // intialize wells that have been there before
120  // order may change so the mapping is based on the well name
121  if( ! prevState.wellMap().empty() )
122  {
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 );
128  if( it != end )
129  {
130  // this is not a new added well
131  is_new_well_[w] = false;
132 
133  const int oldIndex = (*it).second[ 0 ];
134  const int newIndex = w;
135 
136  // bhp
137  bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
138 
139  // thp
140  thp()[ newIndex ] = prevState.thp()[ oldIndex ];
141 
142  // wellrates
143  for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
144  {
145  wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
146  }
147 
148  // perfPhaseRates
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];
152  // copy perforation rates when the number of perforations is equal,
153  // otherwise initialize perfphaserates to well rates divided by the number of perforations.
154  if( num_perf_old_well == num_perf_this_well )
155  {
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 )
159  {
160  perfPhaseRates()[ perf_phase_idx ] = prevState.perfPhaseRates()[ old_perf_phase_idx ];
161  }
162  } else {
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);
166  }
167  }
168  }
169  // perfPressures
170  if( num_perf_old_well == num_perf_this_well )
171  {
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 )
175  {
176  perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
177  }
178  }
179  // perfSolventRates
180  if (pu.has_solvent) {
181  if( num_perf_old_well == num_perf_this_well )
182  {
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 )
186  {
187  perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
188  }
189  }
190  }
191  }
192 
193 
194  // If in the new step, there is no THP related target/limit anymore, its thp value should be
195  // set to zero.
196  const WellControls* ctrl = wells->ctrls[w];
197  const int nwc = well_controls_get_num(ctrl);
198  int ctrl_index = 0;
199  for (; ctrl_index < nwc; ++ctrl_index) {
200  if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
201  break;
202  }
203  }
204  // not finding any thp related control/limits
205  if (ctrl_index == nwc) {
206  thp()[w] = 0.;
207  }
208  }
209  }
210 
211  {
212  // we need to create a trival segment related values to avoid there will be some
213  // multi-segment wells added later.
214  top_segment_index_.reserve(nw);
215  for (int w = 0; w < nw; ++w) {
216  top_segment_index_.push_back(w);
217  }
218  segpress_ = bhp();
219  segrates_ = wellRates();
220  }
221  }
222 
223  template <class State>
224  void resize(const Wells* wells, const State& state, const PhaseUsage& pu) {
225  const WellStateFullyImplicitBlackoil dummy_state{}; // Init with an empty previous state only resizes
226  init(wells, state, dummy_state, pu) ;
227  }
228 
230  std::vector<double>& perfPhaseRates() { return perfphaserates_; }
231  const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
232 
234  std::vector<int>& currentControls() { return current_controls_; }
235  const std::vector<int>& currentControls() const { return current_controls_; }
236 
237  data::Wells report(const PhaseUsage &pu) const override {
238  data::Wells res = WellState::report(pu);
239 
240  const int nw = this->numWells();
241  if( nw == 0 ) return res;
242  const int np = pu.num_phases;
243 
244 
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;
249  }
250 
251  if( pu.phase_used[BlackoilPhases::Liquid] ) {
252  phs.at( pu.phase_pos[BlackoilPhases::Liquid] ) = rt::oil;
253  }
254 
255  if( pu.phase_used[BlackoilPhases::Vapour] ) {
256  phs.at( pu.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
257  }
258 
259  if (pu.has_solvent) {
260  // add solvent component
261  for( int w = 0; w < nw; ++w ) {
262  res.at( wells_->name[ w ]).rates.set( rt::solvent, solventWellRate(w) );
263  }
264  }
265 
266  /* this is a reference or example on **how** to convert from
267  * WellState to something understood by opm-output. it is intended
268  * to be properly implemented and maintained as a part of
269  * simulators, as it relies on simulator internals, details and
270  * representations.
271  */
272 
273  for( const auto& wt : this->wellMap() ) {
274  const auto w = wt.second[ 0 ];
275  auto& well = res.at( wt.first );
276  well.control = this->currentControls()[ w ];
277 
278  int local_comp_index = 0;
279  for( auto& comp : well.completions ) {
280  const auto rates = this->perfPhaseRates().begin()
281  + (np * wt.second[ 1 ])
282  + (np * local_comp_index);
283  ++local_comp_index;
284 
285  for( int i = 0; i < np; ++i ) {
286  comp.rates.set( phs[ i ], *(rates + i) );
287  }
288  }
289  assert(local_comp_index == this->wells_->well_connpos[ w + 1 ] - this->wells_->well_connpos[ w ]);
290  }
291 
292  return res;
293  }
294 
295 
297  template <typename PrevWellState>
298  void initWellStateMSWell(const Wells* wells, const std::vector<const Well*>& wells_ecl,
299  const int time_step, const PhaseUsage& pu, const PrevWellState& prev_well_state)
300  {
301  // still using the order in wells
302  const int nw = wells->number_of_wells;
303  if (nw == 0) {
304  return;
305  }
306 
307  top_segment_index_.clear();
308  top_segment_index_.reserve(nw);
309  segpress_.clear();
310  segpress_.reserve(nw);
311  segrates_.clear();
312  segrates_.reserve(nw * numPhases());
313 
314  nseg_ = 0;
315  // in the init function, the well rates and perforation rates have been initialized or copied from prevState
316  // what we do here, is to set the segment rates and perforation rates
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()) {
323  break;
324  }
325  }
326 
327  // It should be able to find in wells_ecl.
328  if (index_well_ecl == nw_wells_ecl) {
329  OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
330  }
331 
332  const Well* well_ecl = wells_ecl[index_well_ecl];
333  top_segment_index_.push_back(nseg_);
334  if ( !well_ecl->isMultiSegment(time_step) ) { // not multi-segment well
335  nseg_ += 1;
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]);
340  }
341  } else { // it is a multi-segment well
342  const SegmentSet& segment_set = well_ecl->getSegmentSet(time_step);
343  // assuming the order of the perforations in well_ecl is the same with Wells
344  const CompletionSet& completion_set = well_ecl->getCompletions(time_step);
345  // number of segment for this single well
346  const int well_nseg = segment_set.numberSegment();
347  const int nperf = completion_set.size();
348  nseg_ += well_nseg;
349  // we need to know for each segment, how many perforation it has and how many segments using it as outlet_segment
350  // that is why I think we should use a well model to initialize the WellState here
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);
357  }
358 
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);
368  }
369  }
370 
371 
372  // for the segrates_, now it becomes a recursive solution procedure.
373  {
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)); // make sure the information from wells_ecl consistent with wells
378  if (pu.phase_used[Gas]) {
379  const int gaspos = pu.phase_pos[Gas];
380  // scale the phase rates for Gas to avoid too bad initial guess for gas fraction
381  // it will probably benefit the standard well too, while it needs to be justified
382  // TODO: to see if this strategy can benefit StandardWell too
383  // TODO: it might cause big problem for gas rate control or if there is a gas rate limit
384  // maybe the best way is to initialize the fractions first then get the rates
385  for (int perf = 0; perf < nperf; perf++) {
386  const int perf_pos = start_perf + perf;
387  perfPhaseRates()[np * perf_pos + gaspos] *= 100.;
388  }
389  }
390 
391  const std::vector<double> perforation_rates(perfPhaseRates().begin() + np * start_perf,
392  perfPhaseRates().begin() + np * start_perf_next_well); // the perforation rates for this well
393  std::vector<double> segment_rates;
394  calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, segment_rates);
395  std::copy(segment_rates.begin(), segment_rates.end(), std::back_inserter(segrates_));
396  }
397 
398  // for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment
399  // if there is no perforation associated with this segment, it uses the pressure from the outlet segment
400  // which requres the ordering is successful
401  // Not sure what is the best way to handle the initialization, hopefully, the bad initialization can be
402  // improved during the solveWellEq process
403  {
404  // top segment is always the first one, and its pressure is the well bhp
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]);
412  } else {
413  // segpress_.push_back(bhp); // may not be a good decision
414  // using the outlet segment pressure // it needs the ordering is correct
415  const int outlet_seg = segment_set[seg].outletSegment();
416  segpress_.push_back(segpress_[top_segment + segment_set.segmentNumberToIndex(outlet_seg)]);
417  }
418  }
419  }
420  }
421  }
422  assert(int(segpress_.size()) == nseg_);
423  assert(int(segrates_.size()) == nseg_ * numPhases() );
424 
425  if (!prev_well_state.wellMap().empty()) {
426  // copying MS well related
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 );
432 
433  if (it != end) { // the well is found in the prev_well_state
434  // TODO: the well with same name can change a lot, like they might not have same number of segments
435  // we need to handle that later.
436  // for now, we just copy them.
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;
442  // if it is the last well in list
443  if (new_index_well == int(top_segment_index_.size()) - 1) {
444  number_of_segment = nseg_ - new_top_segmnet_index;
445  } else {
446  number_of_segment = topSegmentIndex(new_index_well + 1) - new_top_segmnet_index;
447  }
448 
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];
451  }
452 
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];
455  }
456  }
457  }
458  }
459  }
460 
461 
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)
464  {
465  // the rate of the segment equals to the sum of the contribution from the perforations and inlet segment rates.
466  // the first segment is always the top segment, its rates should be equal to the well rates.
467  assert(segment_inlets.size() == segment_perforations.size());
468  const int well_nseg = segment_inlets.size();
469  if (segment == 0) { // beginning the calculation
470  segment_rates.resize(np * well_nseg, 0.0);
471  }
472  // contributions from the perforations belong to this segment
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];
476  }
477  }
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];
482  }
483  }
484  }
485 
486 
487  bool isNewWell(const int w) const {
488  return is_new_well_[w];
489  }
490 
491 
492  void setNewWell(const int w, const bool is_new_well) {
493  is_new_well_[w] = is_new_well;
494  }
495 
496 
498  std::vector<double>& perfRateSolvent() { return perfRateSolvent_; }
499  const std::vector<double>& perfRateSolvent() const { return perfRateSolvent_; }
500 
502  double solventWellRate(const int w) const {
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];
506  }
507  return solvent_well_rate;
508  }
509 
510  const std::vector<double>& segRates() const
511  {
512  return segrates_;
513  }
514 
515  std::vector<double>& segRates()
516  {
517  return segrates_;
518  }
519 
520  const std::vector<double>& segPress() const
521  {
522  return segpress_;
523  }
524 
525  std::vector<double>& segPress()
526  {
527  return segpress_;
528  }
529 
530  int numSegment() const
531  {
532  return nseg_;
533  }
534 
535  int topSegmentIndex(const int w) const
536  {
537  assert(w < int(top_segment_index_.size()) );
538 
539  return top_segment_index_[w];
540  }
541 
542  private:
543  std::vector<double> perfphaserates_;
544  std::vector<int> current_controls_;
545  std::vector<double> perfRateSolvent_;
546 
547  // marking whether the well is just added
548  // for newly added well, the current initialized rates from WellState
549  // will have very wrong compositions for production wells, will mostly cause
550  // problem with VFP interpolation
551  std::vector<bool> is_new_well_;
552 
553  // MS well related
554  // for StandardWell, the number of segments will be one
555  std::vector<double> segrates_;
556  std::vector<double> segpress_;
557  // the index of the top segments, which is used to locate the
558  // multisegment well related information in WellState
559  std::vector<int> top_segment_index_;
560  int nseg_; // total number of the segments
561 
562  };
563 
564 } // namespace Opm
565 
566 
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