WellsManager_impl.hpp
1 #include <opm/parser/eclipse/Units/Units.hpp>
2 #include <opm/core/grid/GridHelpers.hpp>
3 
4 #include <opm/common/ErrorMacros.hpp>
5 #include <opm/common/OpmLog/OpmLog.hpp>
6 #include <opm/core/utility/compressedToCartesian.hpp>
7 #include <opm/core/props/rock/RockFromDeck.hpp>
8 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
9 #include <opm/parser/eclipse/EclipseState/Schedule/CompletionSet.hpp>
10 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
11 
12 #include <algorithm>
13 #include <array>
14 #include <cstddef>
15 #include <exception>
16 #include <iterator>
17 #include <numeric>
18 
19 namespace WellsManagerDetail
20 {
21 
22 
23  namespace ProductionControl
24  {
25  enum Mode { ORAT, WRAT, GRAT,
26  LRAT, CRAT, RESV,
27  BHP , THP , GRUP };
28  /*
29  namespace Details {
30  std::map<std::string, Mode>
31  init_mode_map();
32  } // namespace Details
33  */
34  Mode mode(const std::string& control);
35 
36 
37  Mode mode(Opm::WellProducer::ControlModeEnum controlMode);
38  } // namespace ProductionControl
39 
40 
41  namespace InjectionControl
42  {
43  enum Mode { RATE, RESV, BHP,
44  THP, GRUP };
45  /*
46  namespace Details {
47  std::map<std::string, Mode>
48  init_mode_map();
49  } // namespace Details
50  */
51  Mode mode(const std::string& control);
52 
53  Mode mode(Opm::WellInjector::ControlModeEnum controlMode);
54 
55  } // namespace InjectionControl
56 
57 double computeWellIndex(const double radius,
58  const std::array<double, 3>& cubical,
59  const double* cell_permeability,
60  const double skin_factor,
61  const Opm::WellCompletion::DirectionEnum direction,
62  const double ntg);
63 
64 template <int dim, class C2F, class FC>
65 std::array<double, dim>
66 getCubeDim(const C2F& c2f,
67  FC begin_face_centroids,
68  int cell)
69 {
70  std::array< std::vector<double>, dim > X;
71  {
72  const std::vector<double>::size_type
73  nf = std::distance(c2f[cell].begin(),
74  c2f[cell].end ());
75 
76  for (int d = 0; d < dim; ++d) {
77  X[d].reserve(nf);
78  }
79  }
80 
81  typedef typename C2F::row_type::const_iterator FI;
82 
83  for (FI f = c2f[cell].begin(), e = c2f[cell].end(); f != e; ++f) {
84  using Opm::UgGridHelpers::increment;
85  using Opm::UgGridHelpers::getCoordinate;
86 
87  const FC& fc = increment(begin_face_centroids, *f, dim);
88 
89  for (int d = 0; d < dim; ++d) {
90  X[d].push_back(getCoordinate(fc, d));
91  }
92  }
93 
94  std::array<double, dim> cube;
95  for (int d = 0; d < dim; ++d) {
96  typedef std::vector<double>::iterator VI;
97  typedef std::pair<VI,VI> PVI;
98 
99  const PVI m = std::minmax_element(X[d].begin(), X[d].end());
100 
101  cube[d] = *m.second - *m.first;
102  }
103 
104  return cube;
105 }
106 } // end namespace WellsManagerDetail
107 
108 namespace Opm
109 {
110 template<class C2F, class FC, class NTG>
111 void WellsManager::createWellsFromSpecs(std::vector<const Well*>& wells, size_t timeStep,
112  const C2F& c2f,
113  const int* cart_dims,
114  FC begin_face_centroids,
115  int dimensions,
116  std::vector<double>& dz,
117  std::vector<std::string>& well_names,
118  std::vector<WellData>& well_data,
119  std::map<std::string, int>& well_names_to_index,
120  const PhaseUsage& phaseUsage,
121  const std::map<int,int>& cartesian_to_compressed,
122  const double* permeability,
123  const NTG& ntg,
124  std::vector<int>& wells_on_proc,
125  const std::unordered_set<std::string>& ignored_wells,
126  const DynamicListEconLimited& list_econ_limited)
127 {
128  if (dimensions != 3) {
129  OPM_THROW(std::domain_error,
130  "WellsManager::createWellsFromSpecs() only "
131  "supported in three space dimensions");
132  }
133 
134  std::vector<std::vector<PerfData> > wellperf_data;
135  wellperf_data.resize(wells.size());
136  wells_on_proc.resize(wells.size(), 1);
137 
138  // The well index on the current process.
139  // Note that some wells are deactivated as they live on the interior
140  // domain of another proccess. Therefore this might different from
141  // the index of the well according to the eclipse state
142  int active_well_index = 0;
143  for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
144  const auto* well = (*wellIter);
145 
146  if (well->getStatus(timeStep) == WellCommon::SHUT) {
147  continue;
148  }
149 
150  if ( ignored_wells.find(well->name()) != ignored_wells.end() ) {
151  wells_on_proc[ wellIter - wells.begin() ] = 0;
152  continue;
153  }
154 
155  if (list_econ_limited.wellShutEconLimited(well->name())) {
156  continue;
157  }
158 
159  std::vector<int> cells_connection_closed;
160  if (list_econ_limited.anyConnectionClosedForWell(well->name())) {
161  cells_connection_closed = list_econ_limited.getClosedConnectionsForWell(well->name());
162  }
163 
164  { // COMPDAT handling
165  // shut completions and open ones stored in this process will have 1 others 0.
166 
167  for(const auto& completion : well->getCompletions(timeStep)) {
168  if (completion.getState() == WellCompletion::OPEN) {
169  int i = completion.getI();
170  int j = completion.getJ();
171  int k = completion.getK();
172 
173  const int* cpgdim = cart_dims;
174  int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k);
175  std::map<int, int>::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx);
176  if (cgit == cartesian_to_compressed.end()) {
177  OPM_MESSAGE("****Warning: Cell with i,j,k indices " << i << ' ' << j << ' '
178  << k << " not found in grid. The completion will be igored (well = "
179  << well->name() << ')');
180  }
181  else
182  {
183  int cell = cgit->second;
184  // check if the connection is closed due to economic limits
185  if (!cells_connection_closed.empty()) {
186  const bool connection_found = std::find(cells_connection_closed.begin(),
187  cells_connection_closed.end(), cell)
188  != cells_connection_closed.end();
189  if (connection_found) {
190  continue;
191  }
192  }
193 
194  PerfData pd;
195  pd.cell = cell;
196  {
197  const Value<double>& transmissibilityFactor = completion.getConnectionTransmissibilityFactorAsValueObject();
198  const double wellPi = completion.getWellPi();
199  if (transmissibilityFactor.hasValue()) {
200  pd.well_index = transmissibilityFactor.getValue();
201  } else {
202  double radius = 0.5*completion.getDiameter();
203  if (radius <= 0.0) {
204  radius = 0.5*unit::feet;
205  OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius);
206  }
207 
208  std::array<double, 3> cubical =
209  WellsManagerDetail::getCubeDim<3>(c2f, begin_face_centroids, cell);
210 
211  // overwrite dz values calculated in getCubeDim.
212  if (dz.size() > 0) {
213  cubical[2] = dz[cell];
214  }
215 
216  const double* cell_perm = &permeability[dimensions*dimensions*cell];
217  pd.well_index =
218  WellsManagerDetail::computeWellIndex(radius, cubical, cell_perm,
219  completion.getSkinFactor(),
220  completion.getDirection(),
221  ntg[cell]);
222  }
223  pd.satnumid = completion.getSatTableId();
224  pd.well_index *= wellPi;
225  }
226  wellperf_data[active_well_index].push_back(pd);
227  }
228  } else {
229  if (completion.getState() != WellCompletion::SHUT) {
230  OPM_THROW(std::runtime_error, "Completion state: " << WellCompletion::StateEnum2String( completion.getState() ) << " not handled");
231  }
232  }
233  }
234  }
235  { // WELSPECS handling
236  well_names_to_index[well->name()] = active_well_index;
237  well_names.push_back(well->name());
238  {
239  WellData wd;
240  wd.reference_bhp_depth = well->getRefDepth( timeStep );
241  wd.welspecsline = -1;
242  if (well->isInjector( timeStep ))
243  wd.type = INJECTOR;
244  else
245  wd.type = PRODUCER;
246 
247  wd.allowCrossFlow = well->getAllowCrossFlow();
248  well_data.push_back(wd);
249  }
250  }
251 
252  active_well_index++;
253  }
254  // Set up reference depths that were defaulted. Count perfs.
255 
256  const int num_wells = well_data.size();
257 
258  int num_perfs = 0;
259  assert (dimensions == 3);
260  for (int w = 0; w < num_wells; ++w) {
261  num_perfs += wellperf_data[w].size();
262  }
263 
264  // Create the well data structures.
265  w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs);
266  if (!w_) {
267  OPM_THROW(std::runtime_error, "Failed creating Wells struct.");
268  }
269 
270 
271  // Add wells.
272  for (int w = 0; w < num_wells; ++w) {
273  const int w_num_perf = wellperf_data[w].size();
274  std::vector<int> perf_cells (w_num_perf);
275  std::vector<double> perf_prodind(w_num_perf);
276  std::vector<int> perf_satnumid(w_num_perf);
277 
278  for (int perf = 0; perf < w_num_perf; ++perf) {
279  perf_cells [perf] = wellperf_data[w][perf].cell;
280  perf_prodind[perf] = wellperf_data[w][perf].well_index;
281  perf_satnumid[perf] = wellperf_data[w][perf].satnumid;
282  }
283 
284  const double* comp_frac = NULL;
285 
286  // We initialize all wells with a null component fraction,
287  // and must (for injection wells) overwrite it later.
288  const int ok =
289  add_well(well_data[w].type,
290  well_data[w].reference_bhp_depth,
291  w_num_perf,
292  comp_frac,
293  perf_cells.data(),
294  perf_prodind.data(),
295  perf_satnumid.data(),
296  well_names[w].c_str(),
297  well_data[w].allowCrossFlow,
298  w_);
299 
300  if (!ok) {
301  OPM_THROW(std::runtime_error,
302  "Failed adding well "
303  << well_names[w]
304  << " to Wells data structure.");
305  }
306  }
307 }
308 
309 template <class C2F, class FC>
311 WellsManager(const Opm::EclipseState& eclipseState,
312  const size_t timeStep,
313  int number_of_cells,
314  const int* global_cell,
315  const int* cart_dims,
316  int dimensions,
317  const C2F& cell_to_faces,
318  FC begin_face_centroids,
319  const DynamicListEconLimited& list_econ_limited,
320  bool is_parallel_run,
321  const std::unordered_set<std::string>& deactivated_wells)
322  : w_(0), is_parallel_run_(is_parallel_run)
323 {
324  init(eclipseState, timeStep, number_of_cells, global_cell,
325  cart_dims, dimensions,
326  cell_to_faces, begin_face_centroids, list_econ_limited, deactivated_wells);
327 }
328 
330 template <class C2F, class FC>
331 void
332 WellsManager::init(const Opm::EclipseState& eclipseState,
333  const size_t timeStep,
334  int number_of_cells,
335  const int* global_cell,
336  const int* cart_dims,
337  int dimensions,
338  const C2F& cell_to_faces,
339  FC begin_face_centroids,
340  const DynamicListEconLimited& list_econ_limited,
341  const std::unordered_set<std::string>& deactivated_wells)
342 {
343  if (dimensions != 3) {
344  OPM_THROW(std::runtime_error,
345  "We cannot initialize wells from a deck unless "
346  "the corresponding grid is 3-dimensional.");
347  }
348 
349  if (eclipseState.getSchedule().numWells() == 0) {
350  OPM_MESSAGE("No wells specified in Schedule section, "
351  "initializing no wells");
352  return;
353  }
354 
355  std::map<int,int> cartesian_to_compressed;
356  setupCompressedToCartesian(global_cell, number_of_cells,
357  cartesian_to_compressed);
358 
359  // Obtain phase usage data.
360  PhaseUsage pu = phaseUsageFromDeck(eclipseState);
361 
362  // These data structures will be filled in this constructor,
363  // then used to initialize the Wells struct.
364  std::vector<std::string> well_names;
365  std::vector<WellData> well_data;
366 
367 
368  // For easy lookup:
369  std::map<std::string, int> well_names_to_index;
370 
371  const auto& schedule = eclipseState.getSchedule();
372  auto wells = schedule.getWells(timeStep);
373  std::vector<int> wells_on_proc;
374 
375  well_names.reserve(wells.size());
376  well_data.reserve(wells.size());
377 
378  typedef GridPropertyAccess::ArrayPolicy::ExtractFromDeck<double> DoubleArray;
379  typedef GridPropertyAccess::Compressed<DoubleArray, GridPropertyAccess::Tag::NTG> NTGArray;
380 
381  DoubleArray ntg_glob(eclipseState, "NTG", 1.0);
382  NTGArray ntg(ntg_glob, global_cell);
383 
384  const auto& eclGrid = eclipseState.getInputGrid();
385 
386  // use cell thickness (dz) from eclGrid
387  // dz overwrites values calculated by WellDetails::getCubeDim
388  std::vector<double> dz(number_of_cells);
389  {
390  std::vector<int> gc = compressedToCartesian(number_of_cells, global_cell);
391  for (int cell = 0; cell < number_of_cells; ++cell) {
392  dz[cell] = eclGrid.getCellThicknes(gc[cell]);
393  }
394  }
395 
396 
397  std::vector<double> interleavedPerm;
399  number_of_cells,
400  global_cell,
401  cart_dims,
402  0.0,
403  interleavedPerm);
404 
405  createWellsFromSpecs(wells, timeStep, cell_to_faces,
406  cart_dims,
407  begin_face_centroids,
408  dimensions,
409  dz,
410  well_names, well_data, well_names_to_index,
411  pu, cartesian_to_compressed, interleavedPerm.data(), ntg,
412  wells_on_proc, deactivated_wells, list_econ_limited);
413 
414  setupWellControls(wells, timeStep, well_names, pu, wells_on_proc, list_econ_limited);
415 
416  {
417  const auto& fieldGroup = schedule.getGroup( "FIELD" );
418  well_collection_.addField(fieldGroup, timeStep, pu);
419 
420  const auto& grouptree = schedule.getGroupTree( timeStep );
421  std::vector< std::string > group_stack = { "FIELD" };
422 
423  do {
424  auto parent = group_stack.back();
425  group_stack.pop_back();
426  const auto& children = grouptree.children( parent );
427  group_stack.insert( group_stack.end(), children.begin(), children.end() );
428 
429  for( const auto& child : children ) {
430  well_collection_.addGroup( schedule.getGroup( child ), parent, timeStep, pu );
431  }
432 
433  } while( !group_stack.empty() );
434  }
435 
436  for (size_t i = 0; i < wells_on_proc.size(); ++i) {
437  // wells_on_proc is a vector of flag to indicate whether a well is on the process
438  if (wells_on_proc[i]) {
439  well_collection_.addWell(wells[i], timeStep, pu);
440  }
441  }
442 
443  well_collection_.setWellsPointer(w_);
444 
445  if (well_collection_.groupControlActive()) {
446  // here does not consider the well potentials related guide rate setting
447  setupGuideRates(wells, timeStep, well_data, well_names_to_index);
448  }
449 
450  // Debug output.
451 #define EXTRA_OUTPUT
452 #ifdef EXTRA_OUTPUT
453  /*
454  std::cout << "\t WELL DATA" << std::endl;
455  for(int i = 0; i< num_wells; ++i) {
456  std::cout << i << ": " << well_data[i].type << " "
457  << well_data[i].control << " " << well_data[i].target
458  << std::endl;
459  }
460 
461  std::cout << "\n\t PERF DATA" << std::endl;
462  for(int i=0; i< int(wellperf_data.size()); ++i) {
463  for(int j=0; j< int(wellperf_data[i].size()); ++j) {
464  std::cout << i << ": " << wellperf_data[i][j].cell << " "
465  << wellperf_data[i][j].well_index << std::endl;
466  }
467  }
468  */
469 #endif
470 }
471 
472 } // end namespace Opm
Well constrained by BHP target.
Definition: well_controls.h:35
Well constrained by THP target.
Definition: well_controls.h:36
Definition: WellsManager.cpp:51
static void extractInterleavedPermeability(const Opm::EclipseState &eclState, const int number_of_cells, const int *global_cell, const int *cart_dims, const double perm_threshold, std::vector< double > &permeability)
Convert the permeabilites for the logically Cartesian grid in EclipseState to an array of size number...
Definition: RockFromDeck.cpp:110
Definition: AnisotropicEikonal.cpp:446
bool groupControlActive() const
Whether we have active group control.
Definition: WellCollection.cpp:327
Well is an injector.
Definition: wells.h:42
struct Wells * create_wells(int nphases, int nwells, int nperf)
Construct a Wells object initially capable of managing a given number of wells and total number of we...
Definition: wells.c:262
WellsManager()
Default constructor – no wells.
Definition: WellsManager.cpp:317
Well is a producer.
Definition: wells.h:43
void setWellsPointer(Wells *wells)
Adds the well pointer to each leaf node (does not take ownership).
Definition: WellCollection.cpp:185
int add_well(enum WellType type, double depth_ref, int nperf, const double *comp_frac, const int *cells, const double *WI, const int *sat_table_id, const char *name, int allow_cf, struct Wells *W)
Append a new well to an existing Wells object.
Definition: wells.c:365