All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
WellStateMultiSegment.hpp
1 /*
2  Copyright 2015 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
21 #define OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDED
22 
23 
24 #include <opm/core/wells.h>
25 #include <opm/core/well_controls.h>
26 #include <opm/common/ErrorMacros.hpp>
27 #include <opm/autodiff/AutoDiffBlock.hpp>
28 // #include <opm/autodiff/WellMultiSegment.hpp>
29 #include <opm/autodiff/MultisegmentWells.hpp>
30 #include <opm/autodiff/WellStateFullyImplicitBlackoil.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 
43  // Since we are avoiding to use the old Wells structure,
44  // it might be a good idea not to relate this State to the WellState much.
47  {
48  public:
49 
51 
52  typedef struct {
53  int well_number;
54  int start_segment;
55  int number_of_segments;
56  int start_perforation;
57  int number_of_perforations;
58  std::vector<int> start_perforation_segment; // the starting position of perforation inside the segment
59  std::vector<int> number_of_perforations_segment; // the numbers for perforations for the segments
61 
62  typedef std::map<std::string, SegmentedMapentryType> SegmentedWellMapType;
63 
67  template <class ReservoirState, class PrevWellState>
68  void init(const MultisegmentWells& ms_wells, const ReservoirState& state, const PrevWellState& prevState, const Wells* legacy_wells_ptr)
69  {
70  // Used by output facilities.
71  this->wells_.reset( clone_wells( legacy_wells_ptr ) );
72 
73  const std::vector<WellMultiSegmentConstPtr>& wells = ms_wells.msWells();
74  const int nw = wells.size();
75  nseg_ = 0;
76  nperf_ = 0;
77  if (nw == 0) {
78  perfPhaseRates().clear();
79  perfPress().clear();
80  segphaserates_.clear();
81  segpress_.clear();
82  return;
83  }
84 
85  const int np = wells[0]->numberOfPhases(); // number of the phases
86 
87  for (int iw = 0; iw < nw; ++iw) {
88  nperf_ += wells[iw]->numberOfPerforations();
89  nseg_ += wells[iw]->numberOfSegments();
90  }
91 
92  bhp().resize(nw);
93  thp().resize(nw);
94  top_segment_loc_.resize(nw);
95  temperature().resize(nw, 273.15 + 20); // standard temperature for now
96 
97  wellRates().resize(nw * np, 0.0);
98 
99  currentControls().resize(nw);
100  for(int iw = 0; iw < nw; ++iw) {
101  currentControls()[iw] = well_controls_get_current(wells[iw]->wellControls());
102  }
103 
104  for (int iw = 0; iw < nw; ++iw) {
105  assert((wells[iw]->wellType() == INJECTOR) || (wells[iw]->wellType() == PRODUCER));
106  }
107 
108  int start_segment = 0;
109  int start_perforation = 0;
110 
111  perfPhaseRates().resize(nperf_ * np, 0.0);
112  perfPress().resize(nperf_, -1.0e100);
113  segphaserates_.resize(nseg_ * np, 0.0);
114  segpress_.resize(nseg_, -1.0e100);
115 
116 
117  for (int w = 0; w < nw; ++w) {
118  assert((wells[w]->wellType() == INJECTOR) || (wells[w]->wellType() == PRODUCER));
119  const WellControls* ctrl = wells[w]->wellControls();
120 
121  std::string well_name(wells[w]->name());
122  // Initialize the wellMap_ here
123  SegmentedMapentryType& wellMapEntry = segmentedWellMap_[well_name];
124  wellMapEntry.well_number = w;
125  wellMapEntry.start_segment = start_segment;
126  wellMapEntry.number_of_segments = wells[w]->numberOfSegments();
127  wellMapEntry.start_perforation = start_perforation;
128  wellMapEntry.number_of_perforations = wells[w]->numberOfPerforations();
129 
130  top_segment_loc_[w] = start_segment;
131 
132  int start_perforation_segment = 0;
133  wellMapEntry.start_perforation_segment.resize(wellMapEntry.number_of_segments);
134  wellMapEntry.number_of_perforations_segment.resize(wellMapEntry.number_of_segments);
135 
136  for (int i = 0; i < wellMapEntry.number_of_segments; ++i) {
137  wellMapEntry.start_perforation_segment[i] = start_perforation_segment;
138  wellMapEntry.number_of_perforations_segment[i] = wells[w]->segmentPerforations()[i].size();
139  start_perforation_segment += wellMapEntry.number_of_perforations_segment[i];
140  }
141  assert(start_perforation_segment == wellMapEntry.number_of_perforations);
142 
143 
144  if (well_controls_well_is_stopped(ctrl)) {
145  // 1. WellRates: 0
146  // 2. Bhp: assign bhp equal to bhp control, if applicable, otherwise
147  // assign equal to first perforation cell pressure.
148  if (well_controls_get_current_type(ctrl) == BHP) {
149  bhp()[w] = well_controls_get_current_target(ctrl);
150  } else {
151  const int first_cell = wells[0]->wellCells()[0];
152  bhp()[w] = state.pressure()[first_cell];
153  }
154  // 3. Thp: assign thp equal to thp control, if applicable,
155  // otherwise assign equal to bhp value.
156  if (well_controls_get_current_type(ctrl) == THP) {
157  thp()[w] = well_controls_get_current_target( ctrl );
158  } else {
159  thp()[w] = bhp()[w];
160  }
161  // 4. Perforation pressures and phase rates
162  // 5. Segment pressures and phase rates
163 
164  } else {
165  // Open Wells
166  // 1. Rates: initialize well rates to match controls if type is SURFACE_RATE. Otherwise, we
167  // cannot set the correct value here, so we aasign a small rate with the correct sign so that any
168  // logic depending on that sign will work as expected.
169  if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
170  const double rate_target = well_controls_get_current_target(ctrl);
171  const double * distr = well_controls_get_current_distr( ctrl );
172  for (int p = 0; p < np; ++p) {
173  wellRates()[np * w + p] = rate_target * distr[p];
174  }
175  } else {
176  const double small_rate = 1e-14;
177  const double sign = (wells[w]->wellType() == INJECTOR) ? 1.0 : -1.0;
178  for (int p = 0; p < np; ++p) {
179  wellRates()[np * w + p] = small_rate * sign;
180  }
181  }
182 
183  // 2. Bhp:
184  if (well_controls_get_current_type(ctrl) == BHP) {
185  bhp()[w] = well_controls_get_current_target(ctrl);
186  } else {
187  const int first_cell = wells[w]->wellCells()[0];
188  const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
189  bhp()[w] = safety_factor* state.pressure()[first_cell];
190  }
191  // 3. Thp:
192  if (well_controls_get_current_type(ctrl) == THP) {
193  thp()[w] = well_controls_get_current_target(ctrl);
194  } else {
195  thp()[w] = bhp()[w];
196  }
197 
198  // 4. Perf rates and pressures
199  int number_of_perforations = wellMapEntry.number_of_perforations;
200  for (int i = 0; i < number_of_perforations; ++i) {
201  for (int p = 0; p < np; ++p) {
202  perfPhaseRates()[np * (i + start_perforation) + p] = wellRates()[np * w + p] / double(number_of_perforations);
203  }
204  if (wells[w]->isMultiSegmented()) {
205  const double safety_factor = (wells[w]->wellType() == INJECTOR) ? 1.01 : 0.99;
206  perfPress()[i + start_perforation] = safety_factor * state.pressure()[wells[w]->wellCells()[i]];
207  } else {
208  perfPress()[i + start_perforation] = state.pressure()[wells[w]->wellCells()[i]];
209  }
210  }
211 
212  // 5. Segment rates and pressures
213  int number_of_segments = wellMapEntry.number_of_segments;
214  // the seg_pressure is the same with the first perf_pressure. For the top segment, it is the same with bhp,
215  // when under bhp control.
216  // the seg_rates will related to the sum of the perforation rates, and also trying to keep consistent with the
217  // well rates. Most importantly, the segment rates of the top segment is the same with the well rates.
218  segpress_[start_segment] = bhp()[w];
219  for (int i = 1; i < number_of_segments; ++i) {
220  int first_perforation_segment = start_perforation + wellMapEntry.start_perforation_segment[i];
221  segpress_[i + start_segment] = perfPress()[first_perforation_segment];
222  // the segmnent pressure of the top segment should be the bhp
223  }
224 
225  for (int p = 0; p < np; ++p) {
226  Eigen::VectorXd v_perf_rates(number_of_perforations);
227  for (int i = 0; i < number_of_perforations; ++i) {
228  v_perf_rates[i] = perfPhaseRates()[np * (i + start_perforation) + p];
229  }
230 
231  Eigen::VectorXd v_segment_rates = wells[w]->wellOps().p2s_gather * v_perf_rates;
232 
233  for (int i = 0; i < number_of_segments; ++i) {
234  segphaserates_[np * (i + start_segment) + p] = v_segment_rates[i];
235  }
236  }
237  }
238 
239  start_segment += wellMapEntry.number_of_segments;
240  start_perforation += wellMapEntry.number_of_perforations;
241 
242  }
243 
244  // Initialize current_controls_.
245  // The controls set in the Wells object are treated as defaults,
246  // and also used for initial values.
247  currentControls().resize(nw);
248  for (int w = 0; w < nw; ++w) {
249  currentControls()[w] = well_controls_get_current(wells[w]->wellControls());
250  }
251 
252  // initialize wells that have been there before
253  // order can change so the mapping is based on the well names
254  if ( !(prevState.segmentedWellMap().empty()) )
255  {
256  typedef typename SegmentedWellMapType::const_iterator const_iterator;
257  const_iterator end_old = prevState.segmentedWellMap().end();
258 
259  for (int w = 0; w < nw; ++w) {
260  std::string well_name(wells[w]->name());
261  const_iterator it_old = prevState.segmentedWellMap().find(well_name);
262  const_iterator it_this = segmentedWellMap().find(well_name);
263 
264  assert(it_this != segmentedWellMap().end()); // the current well must be present in the current well map
265 
266  if (it_old != end_old) {
267  const int oldIndex = (*it_old).second.well_number;
268  const int newIndex = w;
269 
270  // bhp
271  bhp()[newIndex] = prevState.bhp()[oldIndex];
272 
273  // well rates
274  for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
275  {
276  wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
277  }
278 
279  const int num_seg_old_well = (*it_old).second.number_of_segments;
280  const int num_perf_old_well = (*it_old).second.number_of_perforations;
281 
282  const int num_seg_this_well = (*it_this).second.number_of_segments;
283  const int num_perf_this_well = (*it_this).second.number_of_perforations;
284 
285  // determing if the structure of the wells has been changed by comparing the number of segments and perforations
286  // may not be very safe.
287  // The strategy HAS to be changed later with experiments and analysis
288  if ((num_perf_old_well == num_perf_this_well) && (num_seg_old_well == num_seg_this_well)) {
289  const int old_start_perforation = (*it_old).second.start_perforation;
290  const int old_start_segment = (*it_old).second.start_segment;
291 
292  const int this_start_perforation = (*it_this).second.start_perforation;
293  const int this_start_segment = (*it_this).second.start_segment;
294 
295  // this is not good when the the well rates changed dramatically
296  for (int i = 0; i < num_seg_this_well * np; ++i) {
297  segphaserates_[this_start_segment * np + i] = prevState.segPhaseRates()[old_start_segment * np + i];
298  }
299 
300  for (int i = 0; i < num_perf_this_well * np; ++i) {
301  perfPhaseRates()[this_start_perforation * np + i] = prevState.perfPhaseRates()[old_start_perforation * np + i];
302  }
303 
304  // perf_pressures_
305  for (int i = 0; i < num_perf_this_well; ++i) {
306  // p
307  perfPress()[this_start_perforation + i] = prevState.perfPress()[old_start_perforation + i];
308  }
309 
310  // segpress_
311  for (int i = 0; i < num_seg_this_well; ++i) {
312  // p
313  segpress_[this_start_segment + i] = prevState.segPress()[old_start_segment + i];
314  }
315 
316  // current controls
317  const int old_control_index = prevState.currentControls()[ oldIndex ];
318  if (old_control_index < well_controls_get_num(wells[w]->wellControls())) {
319  // If the set of controls have changed, this may not be identical
320  // to the last control, but it must be a valid control.
321  currentControls()[ newIndex ] = old_control_index;
322  }
323  }
324  }
325  }
326  }
327  }
328 
329 
330  std::vector<double>& segPress() { return segpress_; }
331  const std::vector<double>& segPress() const { return segpress_; }
332 
333  std::vector<double>& segPhaseRates() { return segphaserates_; }
334  const std::vector<double>& segPhaseRates() const { return segphaserates_; }
335 
336  const std::vector<int>& topSegmentLoc() const { return top_segment_loc_; };
337 
338  const SegmentedWellMapType& segmentedWellMap() const { return segmentedWellMap_; }
339  SegmentedWellMapType& segmentedWellMap() { return segmentedWellMap_; }
340 
341  int numSegments() const { return nseg_; }
342  int numPerforations() const { return nperf_; }
343 
344  private:
345  // pressure for the segment nodes
346  std::vector<double> segpress_;
347  // phase rates for the segments
348  std::vector<double> segphaserates_;
349 
350  // the location of the top segments within the whole segment list
351  // it is better in the Wells class if we have a class instead of
352  // using a vector for all the wells
353  std::vector<int> top_segment_loc_;
354 
355  SegmentedWellMapType segmentedWellMap_;
356 
357  int nseg_;
358  int nperf_;
359  };
360 
361 } // namespace Opm
362 
363 
364 #endif // OPM_WELLSTATEMULTISEGMENT_HEADER_INCLUDE
void init(const MultisegmentWells &ms_wells, const ReservoirState &state, const PrevWellState &prevState, const Wells *legacy_wells_ptr)
Allocate and initialize if wells is non-null.
Definition: WellStateMultiSegment.hpp:68
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
The state of a set of multi-sgemnet wells.
Definition: WellStateMultiSegment.hpp:45
Definition: WellStateMultiSegment.hpp:52
Class for handling the multi-segment well model.
Definition: MultisegmentWells.hpp:55
Eigen::ArrayXd sign(const Eigen::ArrayXd &x)
Return a vector of (-1.0, 0.0 or 1.0), depending on sign per element.
Definition: AutoDiffHelpers.hpp:719