BlackoilWellModel_impl.hpp
1 
2 
3 
4 namespace Opm {
5 
6  template<typename TypeTag>
7  BlackoilWellModel<TypeTag>::
8  BlackoilWellModel(const Wells* wells_arg,
9  WellCollection* well_collection,
10  const std::vector< const Well* >& wells_ecl,
11  const ModelParameters& param,
12  const RateConverterType& rate_converter,
13  const bool terminal_output,
14  const int current_timeIdx,
15  const std::vector<int>& pvt_region_idx)
16  : wells_active_(wells_arg!=nullptr)
17  , wells_(wells_arg)
18  , wells_ecl_(wells_ecl)
19  , number_of_wells_(wells_arg ? (wells_arg->number_of_wells) : 0)
20  , number_of_phases_(wells_arg ? (wells_arg->number_of_phases) : 0) // TODO: not sure if it is proper for this way
21  , param_(param)
22  , well_container_(createWellContainer(wells_arg, wells_ecl, param.use_multisegment_well_, current_timeIdx, param) )
23  , well_collection_(well_collection)
24  , terminal_output_(terminal_output)
25  , has_solvent_(GET_PROP_VALUE(TypeTag, EnableSolvent))
26  , has_polymer_(GET_PROP_VALUE(TypeTag, EnablePolymer))
27  , current_timeIdx_(current_timeIdx)
28  , rate_converter_(rate_converter)
29  , pvt_region_idx_(pvt_region_idx)
30  {
31  }
32 
33 
34 
35 
36 
37  template<typename TypeTag>
38  void
39  BlackoilWellModel<TypeTag>::
40  init(const PhaseUsage phase_usage_arg,
41  const std::vector<bool>& active_arg,
42  const double gravity_arg,
43  const std::vector<double>& depth_arg,
44  long int global_nc,
45  const Grid& grid)
46  {
47  // has to be set always for the convergence check!
48  global_nc_ = global_nc;
49 
50  phase_usage_ = phase_usage_arg;
51  active_ = active_arg;
52 
53  if ( ! localWellsActive() ) {
54  return;
55  }
56 
57  calculateEfficiencyFactors();
58 
59 #ifndef NDEBUG
60  const auto& pu = phase_usage_;
61  const int np = pu.num_phases;
62 
63  // assumes the gas fractions are stored after water fractions
64  // WellVariablePositions needs to be changed for 2p runs
65  assert (np == 3 || (np == 2 && !pu.phase_used[Gas]) );
66 #endif
67 
68  if (has_polymer_)
69  {
70  if (PolymerModule::hasPlyshlog()) {
71  computeRepRadiusPerfLength(grid);
72  }
73  }
74 
75  number_of_cells_ = Opm::UgGridHelpers::numCells(grid);
76  // do the initialization for all the wells
77  // TODO: to see whether we can postpone of the intialization of the well containers to
78  // optimize the usage of the following several member variables
79  for (auto& well : well_container_) {
80  well->init(&phase_usage_, &active_, depth_arg, gravity_arg, number_of_cells_);
81  }
82  }
83 
84 
85 
86 
87 
88  template<typename TypeTag>
89  void
90  BlackoilWellModel<TypeTag>::
91  setVFPProperties(const VFPProperties* vfp_properties_arg)
92  {
93  for (auto& well : well_container_) {
94  well->setVFPProperties(vfp_properties_arg);
95  }
96  }
97 
98 
99 
100 
101 
102 
103  template<typename TypeTag>
104  int
105  BlackoilWellModel<TypeTag>::
106  numWells() const
107  {
108  return number_of_wells_;
109  }
110 
111 
112 
113 
114 
115  template<typename TypeTag>
116  std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
117  BlackoilWellModel<TypeTag>::
118  createWellContainer(const Wells* wells,
119  const std::vector< const Well* >& wells_ecl,
120  const bool use_multisegment_well,
121  const int time_step,
122  const ModelParameters& param)
123  {
124  std::vector<WellInterfacePtr> well_container;
125 
126  const int nw = wells ? (wells->number_of_wells) : 0;
127 
128  if (nw > 0) {
129  well_container.reserve(nw);
130 
131  // With the following way, it will have the same order with wells struct
132  // Hopefully, it can generate the same residual history with master branch
133  for (int w = 0; w < nw; ++w) {
134  const std::string well_name = std::string(wells->name[w]);
135 
136  // finding the location of the well in wells_ecl
137  const int nw_wells_ecl = wells_ecl.size();
138  int index_well = 0;
139  for (; index_well < nw_wells_ecl; ++index_well) {
140  if (well_name == wells_ecl[index_well]->name()) {
141  break;
142  }
143  }
144 
145  // It should be able to find in wells_ecl.
146  if (index_well == nw_wells_ecl) {
147  OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
148  }
149 
150  const Well* well_ecl = wells_ecl[index_well];
151 
152  if ( !well_ecl->isMultiSegment(time_step) || !use_multisegment_well) {
153  well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, time_step, wells, param) );
154  } else {
155  well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, time_step, wells, param) );
156  }
157  }
158  }
159  return well_container;
160  }
161 
162 
163 
164 
165 
166  template<typename TypeTag>
167  SimulatorReport
168  BlackoilWellModel<TypeTag>::
169  assemble(Simulator& ebosSimulator,
170  const int iterationIdx,
171  const double dt,
172  WellState& well_state)
173  {
174  SimulatorReport report;
175  if ( ! wellsActive() ) {
176  return report;
177  }
178 
179  if (iterationIdx == 0) {
180  prepareTimeStep(ebosSimulator, well_state);
181  }
182 
183  updateWellControls(well_state);
184  // Set the well primary variables based on the value of well solutions
185  initPrimaryVariablesEvaluation();
186 
187  if (iterationIdx == 0) {
188  calculateExplicitQuantities(ebosSimulator, well_state);
189  }
190 
191  if (param_.solve_welleq_initially_ && iterationIdx == 0) {
192  // solve the well equations as a pre-processing step
193  report = solveWellEq(ebosSimulator, dt, well_state);
194  }
195  assembleWellEq(ebosSimulator, dt, well_state, false);
196 
197  report.converged = true;
198  return report;
199  }
200 
201 
202 
203 
204 
205  template<typename TypeTag>
206  void
207  BlackoilWellModel<TypeTag>::
208  assembleWellEq(Simulator& ebosSimulator,
209  const double dt,
210  WellState& well_state,
211  bool only_wells) const
212  {
213  for (int w = 0; w < number_of_wells_; ++w) {
214  well_container_[w]->assembleWellEq(ebosSimulator, dt, well_state, only_wells);
215  }
216  }
217 
218 
219 
220 
221 
222  // applying the well residual to reservoir residuals
223  // r = r - duneC_^T * invDuneD_ * resWell_
224  template<typename TypeTag>
225  void
226  BlackoilWellModel<TypeTag>::
227  apply( BVector& r) const
228  {
229  if ( ! localWellsActive() ) {
230  return;
231  }
232 
233  for (auto& well : well_container_) {
234  well->apply(r);
235  }
236  }
237 
238 
239 
240 
241 
242  // Ax = A x - C D^-1 B x
243  template<typename TypeTag>
244  void
245  BlackoilWellModel<TypeTag>::
246  apply(const BVector& x, BVector& Ax) const
247  {
248  // TODO: do we still need localWellsActive()?
249  if ( ! localWellsActive() ) {
250  return;
251  }
252 
253  for (auto& well : well_container_) {
254  well->apply(x, Ax);
255  }
256  }
257 
258 
259 
260 
261 
262  // Ax = Ax - alpha * C D^-1 B x
263  template<typename TypeTag>
264  void
265  BlackoilWellModel<TypeTag>::
266  applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const
267  {
268  if ( ! localWellsActive() ) {
269  return;
270  }
271 
272  if( scaleAddRes_.size() != Ax.size() ) {
273  scaleAddRes_.resize( Ax.size() );
274  }
275 
276  scaleAddRes_ = 0.0;
277  // scaleAddRes_ = - C D^-1 B x
278  apply( x, scaleAddRes_ );
279  // Ax = Ax + alpha * scaleAddRes_
280  Ax.axpy( alpha, scaleAddRes_ );
281  }
282 
283 
284 
285 
286 
287  template<typename TypeTag>
288  void
289  BlackoilWellModel<TypeTag>::
290  recoverWellSolutionAndUpdateWellState(const BVector& x, WellState& well_state) const
291  {
292  for (auto& well : well_container_) {
293  well->recoverWellSolutionAndUpdateWellState(x, well_state);
294  }
295  }
296 
297 
298 
299 
300 
301  template<typename TypeTag>
302  int
303  BlackoilWellModel<TypeTag>::
304  flowPhaseToEbosPhaseIdx( const int phaseIdx ) const
305  {
306  const auto& pu = phase_usage_;
307  if (active_[Water] && pu.phase_pos[Water] == phaseIdx)
308  return FluidSystem::waterPhaseIdx;
309  if (active_[Oil] && pu.phase_pos[Oil] == phaseIdx)
310  return FluidSystem::oilPhaseIdx;
311  if (active_[Gas] && pu.phase_pos[Gas] == phaseIdx)
312  return FluidSystem::gasPhaseIdx;
313 
314  assert(phaseIdx < 3);
315  // for other phases return the index
316  return phaseIdx;
317  }
318 
319 
320 
321 
322 
323  template<typename TypeTag>
324  int
325  BlackoilWellModel<TypeTag>::
326  numPhases() const
327  {
328  return number_of_phases_;
329  }
330 
331 
332 
333 
334 
335  template<typename TypeTag>
336  void
337  BlackoilWellModel<TypeTag>::
338  resetWellControlFromState(const WellState& xw) const
339  {
340  const int nw = wells_->number_of_wells;
341  for (int w = 0; w < nw; ++w) {
342  WellControls* wc = wells_->ctrls[w];
343  well_controls_set_current( wc, xw.currentControls()[w]);
344  }
345  }
346 
347 
348 
349 
350 
351  template<typename TypeTag>
352  bool
355  {
356  return wells_active_;
357  }
358 
359 
360 
361 
362 
363  template<typename TypeTag>
364  void
366  setWellsActive(const bool wells_active)
367  {
368  wells_active_ = wells_active;
369  }
370 
371 
372 
373 
374 
375  template<typename TypeTag>
376  bool
379  {
380  return number_of_wells_ > 0;
381  }
382 
383 
384 
385 
386 
387  template<typename TypeTag>
388  void
391  {
392  for (auto& well : well_container_) {
393  well->initPrimaryVariablesEvaluation();
394  }
395  }
396 
397 
398 
399 
400 
401  template<typename TypeTag>
402  SimulatorReport
403  BlackoilWellModel<TypeTag>::
404  solveWellEq(Simulator& ebosSimulator,
405  const double dt,
406  WellState& well_state) const
407  {
408  const int nw = number_of_wells_;
409  WellState well_state0 = well_state;
410 
411  const int numComp = numComponents();
412  std::vector< Scalar > B_avg( numComp, Scalar() );
413  computeAverageFormationFactor(ebosSimulator, B_avg);
414 
415  const int max_iter = param_.max_welleq_iter_;
416 
417  int it = 0;
418  bool converged;
419  do {
420  assembleWellEq(ebosSimulator, dt, well_state, true);
421 
422  converged = getWellConvergence(ebosSimulator, B_avg);
423 
424  // checking whether the group targets are converged
425  if (wellCollection()->groupControlActive()) {
426  converged = converged && wellCollection()->groupTargetConverged(well_state.wellRates());
427  }
428 
429  if (converged) {
430  break;
431  }
432 
433  ++it;
434  if( localWellsActive() )
435  {
436  for (auto& well : well_container_) {
437  well->solveEqAndUpdateWellState(well_state);
438  }
439  }
440  // updateWellControls uses communication
441  // Therefore the following is executed if there
442  // are active wells anywhere in the global domain.
443  if( wellsActive() )
444  {
445  updateWellControls(well_state);
446  initPrimaryVariablesEvaluation();
447  }
448  } while (it < max_iter);
449 
450  if (converged) {
451  if ( terminal_output_ ) {
452  OpmLog::debug("Well equation solution gets converged with " + std::to_string(it) + " iterations");
453  }
454  } else {
455  if ( terminal_output_ ) {
456  OpmLog::debug("Well equation solution failed in getting converged with " + std::to_string(it) + " iterations");
457  }
458 
459  well_state = well_state0;
460  updatePrimaryVariables(well_state);
461  // also recover the old well controls
462  for (int w = 0; w < nw; ++w) {
463  WellControls* wc = well_container_[w]->wellControls();
464  well_controls_set_current(wc, well_state.currentControls()[w]);
465  }
466  }
467 
468  SimulatorReport report;
469  report.converged = converged;
470  report.total_well_iterations = it;
471  return report;
472  }
473 
474 
475 
476 
477 
478  template<typename TypeTag>
479  bool
480  BlackoilWellModel<TypeTag>::
481  getWellConvergence(const Simulator& ebosSimulator,
482  const std::vector<Scalar>& B_avg) const
483  {
484  ConvergenceReport report;
485 
486  for (const auto& well : well_container_) {
487  report += well->getWellConvergence(B_avg);
488  }
489 
490  // checking NaN residuals
491  {
492  bool nan_residual_found = report.nan_residual_found;
493  const auto& grid = ebosSimulator.gridManager().grid();
494  int value = nan_residual_found ? 1 : 0;
495 
496  nan_residual_found = grid.comm().max(value);
497 
498  if (nan_residual_found) {
499  for (const auto& well : report.nan_residual_wells) {
500  OpmLog::debug("NaN residual found with phase " + well.phase_name + " for well " + well.well_name);
501  }
502  OPM_THROW(Opm::NumericalProblem, "NaN residual found!");
503  }
504  }
505 
506  // checking too large residuals
507  {
508  bool too_large_residual_found = report.too_large_residual_found;
509  const auto& grid = ebosSimulator.gridManager().grid();
510  int value = too_large_residual_found ? 1 : 0;
511 
512  too_large_residual_found = grid.comm().max(value);
513  if (too_large_residual_found) {
514  for (const auto& well : report.too_large_residual_wells) {
515  OpmLog::debug("Too large residual found with phase " + well.phase_name + " fow well " + well.well_name);
516  }
517  OPM_THROW(Opm::NumericalProblem, "Too large residual found!");
518  }
519  }
520 
521  // checking convergence
522  bool converged_well = report.converged;
523  {
524  const auto& grid = ebosSimulator.gridManager().grid();
525  int value = converged_well ? 1 : 0;
526 
527  converged_well = grid.comm().min(value);
528  }
529 
530  return converged_well;
531  }
532 
533 
534 
535 
536 
537  template<typename TypeTag>
538  void
540  calculateExplicitQuantities(const Simulator& ebosSimulator,
541  const WellState& xw) const
542  {
543  for (auto& well : well_container_) {
544  well->calculateExplicitQuantities(ebosSimulator, xw);
545  }
546  }
547 
548 
549 
550 
551 
552  template<typename TypeTag>
553  void
555  updateWellControls(WellState& xw) const
556  {
557  // Even if there no wells active locally, we cannot
558  // return as the Destructor of the WellSwitchingLogger
559  // uses global communication. For no well active globally
560  // we simply return.
561  if( !wellsActive() ) return ;
562 
564 
565  for (const auto& well : well_container_) {
566  well->updateWellControl(xw, logger);
567  }
568 
569  updateGroupControls(xw);
570  }
571 
572 
573 
574 
575 
576  template<typename TypeTag>
577  void
579  updateListEconLimited(const Schedule& schedule,
580  const int current_step,
581  const Wells* wells_struct,
582  const WellState& well_state,
583  DynamicListEconLimited& list_econ_limited) const
584  {
585  for (const auto& well : well_container_) {
586  well->updateListEconLimited(well_state, list_econ_limited);
587  }
588  }
589 
590 
591 
592 
593 
594  template<typename TypeTag>
595  void
597  computeWellPotentials(const Simulator& ebosSimulator,
598  const WellState& well_state,
599  std::vector<double>& well_potentials) const
600  {
601  // number of wells and phases
602  const int nw = number_of_wells_;
603  const int np = number_of_phases_;
604  well_potentials.resize(nw * np, 0.0);
605 
606  for (int w = 0; w < nw; ++w) {
607  std::vector<double> potentials;
608  well_container_[w]->computeWellPotentials(ebosSimulator, well_state, potentials);
609 
610  // putting the sucessfully calculated potentials to the well_potentials
611  for (int p = 0; p < np; ++p) {
612  well_potentials[w * np + p] = std::abs(potentials[p]);
613  }
614  } // end of for (int w = 0; w < nw; ++w)
615  }
616 
617 
618 
619 
620 
621  template<typename TypeTag>
622  void
623  BlackoilWellModel<TypeTag>::
624  prepareTimeStep(const Simulator& ebos_simulator,
625  WellState& well_state) const
626  {
627  // after restarting, the well_controls can be modified while
628  // the well_state still uses the old control index
629  // we need to synchronize these two.
630  // keep in mind that we set the control index of well_state to be the same with
631  // with the wellControls from the deck when we create well_state at the beginning of the report step
632  resetWellControlFromState(well_state);
633 
634  // process group control related
635  prepareGroupControl(ebos_simulator, well_state);
636 
637  // since the controls are all updated, we should update well_state accordingly
638  for (int w = 0; w < number_of_wells_; ++w) {
639  WellControls* wc = well_container_[w]->wellControls();
640  const int control = well_controls_get_current(wc);
641  well_state.currentControls()[w] = control;
642  // TODO: for VFP control, the perf_densities are still zero here, investigate better
643  // way to handle it later.
644  well_container_[w]->updateWellStateWithTarget(control, well_state);
645 
646  // The wells are not considered to be newly added
647  // for next time step
648  if (well_state.isNewWell(w) ) {
649  well_state.setNewWell(w, false);
650  }
651  } // end of for (int w = 0; w < nw; ++w)
652  }
653 
654 
655 
656 
657 
658  template<typename TypeTag>
659  void
660  BlackoilWellModel<TypeTag>::
661  prepareGroupControl(const Simulator& ebos_simulator,
662  WellState& well_state) const
663  {
664  // group control related processing
665  if (well_collection_->groupControlActive()) {
666  for (int w = 0; w < number_of_wells_; ++w) {
667  WellControls* wc = well_container_[w]->wellControls();
668  WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name());
669 
670  // handling the situation that wells do not have a valid control
671  // it happens the well specified with GRUP and restarting due to non-convergencing
672  // putting the well under group control for this situation
673  int ctrl_index = well_controls_get_current(wc);
674 
675  const int group_control_index = well_node.groupControlIndex();
676  if (group_control_index >= 0 && ctrl_index < 0) {
677  // put well under group control
678  well_controls_set_current(wc, group_control_index);
679  well_state.currentControls()[w] = group_control_index;
680  }
681 
682  // Final step, update whehter the well is under group control or individual control
683  // updated ctrl_index from the well control
684  ctrl_index = well_controls_get_current(wc);
685  if (well_node.groupControlIndex() >= 0 && ctrl_index == well_node.groupControlIndex()) {
686  // under group control
687  well_node.setIndividualControl(false);
688  } else {
689  // individual control
690  well_node.setIndividualControl(true);
691  }
692  }
693 
694  if (well_collection_->requireWellPotentials()) {
695 
696  // calculate the well potentials
697  std::vector<double> well_potentials;
698  computeWellPotentials(ebos_simulator, well_state, well_potentials);
699 
700  // update/setup guide rates for each well based on the well_potentials
701  // TODO: this is one of two places that still need Wells struct. In this function, only the well names
702  // well types are used, probably the order of the wells to locate the correct values in well_potentials.
703  well_collection_->setGuideRatesWithPotentials(wells_, phase_usage_, well_potentials);
704  }
705 
706  applyVREPGroupControl(well_state);
707 
708  if (!wellCollection()->groupControlApplied()) {
709  wellCollection()->applyGroupControls();
710  } else {
711  wellCollection()->updateWellTargets(well_state.wellRates());
712  }
713  }
714  }
715 
716 
717 
718 
719 
720  template<typename TypeTag>
721  WellCollection*
722  BlackoilWellModel<TypeTag>::
723  wellCollection() const
724  {
725  return well_collection_;
726  }
727 
728 
729 
730 
731 
732  template<typename TypeTag>
733  void
734  BlackoilWellModel<TypeTag>::
735  calculateEfficiencyFactors()
736  {
737  if ( !localWellsActive() ) {
738  return;
739  }
740 
741  const int nw = number_of_wells_;
742 
743  for (int w = 0; w < nw; ++w) {
744  const std::string well_name = well_container_[w]->name();
745  const WellNode& well_node = wellCollection()->findWellNode(well_name);
746 
747  const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
748 
749  well_container_[w]->setWellEfficiencyFactor(well_efficiency_factor);
750  }
751  }
752 
753 
754 
755 
756 
757  template<typename TypeTag>
758  void
759  BlackoilWellModel<TypeTag>::
760  computeWellVoidageRates(const WellState& well_state,
761  std::vector<double>& well_voidage_rates,
762  std::vector<double>& voidage_conversion_coeffs) const
763  {
764  if ( !localWellsActive() ) {
765  return;
766  }
767  // TODO: for now, we store the voidage rates for all the production wells.
768  // For injection wells, the rates are stored as zero.
769  // We only store the conversion coefficients for all the injection wells.
770  // Later, more delicate model will be implemented here.
771  // And for the moment, group control can only work for serial running.
772  const int nw = numWells();
773  const int np = numPhases();
774 
775  // we calculate the voidage rate for each well, that means the sum of all the phases.
776  well_voidage_rates.resize(nw, 0);
777  // store the conversion coefficients, while only for the use of injection wells.
778  voidage_conversion_coeffs.resize(nw * np, 1.0);
779 
780  std::vector<double> well_rates(np, 0.0);
781  std::vector<double> convert_coeff(np, 1.0);
782 
783  for (int w = 0; w < nw; ++w) {
784  const bool is_producer = well_container_[w]->wellType() == PRODUCER;
785  const int well_cell_top = well_container_[w]->cells()[0];
786  const int pvtRegionIdx = pvt_region_idx_[well_cell_top];
787 
788  // not sure necessary to change all the value to be positive
789  if (is_producer) {
790  std::transform(well_state.wellRates().begin() + np * w,
791  well_state.wellRates().begin() + np * (w + 1),
792  well_rates.begin(), std::negate<double>());
793 
794  // the average hydrocarbon conditions of the whole field will be used
795  const int fipreg = 0; // Not considering FIP for the moment.
796 
797  rate_converter_.calcCoeff(fipreg, pvtRegionIdx, convert_coeff);
798  well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
799  convert_coeff.begin(), 0.0);
800  } else {
801  // TODO: Not sure whether will encounter situation with all zero rates
802  // and whether it will cause problem here.
803  std::copy(well_state.wellRates().begin() + np * w,
804  well_state.wellRates().begin() + np * (w + 1),
805  well_rates.begin());
806  // the average hydrocarbon conditions of the whole field will be used
807  const int fipreg = 0; // Not considering FIP for the moment.
808  rate_converter_.calcCoeff(fipreg, pvtRegionIdx, convert_coeff);
809  std::copy(convert_coeff.begin(), convert_coeff.end(),
810  voidage_conversion_coeffs.begin() + np * w);
811  }
812  }
813  }
814 
815 
816 
817 
818 
819  template<typename TypeTag>
820  void
821  BlackoilWellModel<TypeTag>::
822  applyVREPGroupControl(WellState& well_state) const
823  {
824  if ( wellCollection()->havingVREPGroups() ) {
825  std::vector<double> well_voidage_rates;
826  std::vector<double> voidage_conversion_coeffs;
827  computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs);
828  wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
829 
830  // for the wells under group control, update the control index for the well_state and well_controls
831  for (const WellNode* well_node : wellCollection()->getLeafNodes()) {
832  if (well_node->isInjector() && !well_node->individualControl()) {
833  const int well_index = well_node->selfIndex();
834  well_state.currentControls()[well_index] = well_node->groupControlIndex();
835 
836  WellControls* wc = well_container_[well_index]->wellControls();
837  well_controls_set_current(wc, well_node->groupControlIndex());
838  }
839  }
840  }
841  }
842 
843 
844 
845 
846 
847  template<typename TypeTag>
848  void
849  BlackoilWellModel<TypeTag>::
850  updateGroupControls(WellState& well_state) const
851  {
852 
853  if (wellCollection()->groupControlActive()) {
854  for (int w = 0; w < number_of_wells_; ++w) {
855  // update whether well is under group control
856  // get well node in the well collection
857  WellNode& well_node = well_collection_->findWellNode(well_container_[w]->name());
858 
859  // update whehter the well is under group control or individual control
860  const int current = well_state.currentControls()[w];
861  if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
862  // under group control
863  well_node.setIndividualControl(false);
864  } else {
865  // individual control
866  well_node.setIndividualControl(true);
867  }
868  }
869 
870  applyVREPGroupControl(well_state);
871  // upate the well targets following group controls
872  // it will not change the control mode, only update the targets
873  wellCollection()->updateWellTargets(well_state.wellRates());
874 
875  for (int w = 0; w < number_of_wells_; ++w) {
876  // TODO: check whether we need current argument in updateWellStateWithTarget
877  // maybe there is some circumstances that the current is different from the one
878  // in the WellState.
879  // while probalby, the current argument can be removed
880  const int current = well_state.currentControls()[w];
881  well_container_[w]->updateWellStateWithTarget(current, well_state);
882  }
883  }
884  }
885 
886 
887 
888 
889 
890  template<typename TypeTag>
891  void
892  BlackoilWellModel<TypeTag>::
893  setupCompressedToCartesian(const int* global_cell, int number_of_cells, std::map<int,int>& cartesian_to_compressed ) const
894  {
895  if (global_cell) {
896  for (int i = 0; i < number_of_cells; ++i) {
897  cartesian_to_compressed.insert(std::make_pair(global_cell[i], i));
898  }
899  }
900  else {
901  for (int i = 0; i < number_of_cells; ++i) {
902  cartesian_to_compressed.insert(std::make_pair(i, i));
903  }
904  }
905 
906  }
907 
908  template<typename TypeTag>
909  void
910  BlackoilWellModel<TypeTag>::
911  computeRepRadiusPerfLength(const Grid& grid)
912  {
913  // TODO, the function does not work for parallel running
914  // to be fixed later.
915  const int* global_cell = Opm::UgGridHelpers::globalCell(grid);
916 
917  std::map<int,int> cartesian_to_compressed;
918  setupCompressedToCartesian(global_cell, number_of_cells_,
919  cartesian_to_compressed);
920 
921  for (const auto& well : well_container_) {
922  well->computeRepRadiusPerfLength(grid, cartesian_to_compressed);
923  }
924  }
925 
926 
927 
928 
929 
930  template<typename TypeTag>
931  void
932  BlackoilWellModel<TypeTag>::
933  computeAverageFormationFactor(const Simulator& ebosSimulator,
934  std::vector<double>& B_avg) const
935  {
936  const int np = numPhases();
937 
938  const auto& grid = ebosSimulator.gridManager().grid();
939  const auto& gridView = grid.leafGridView();
940  ElementContext elemCtx(ebosSimulator);
941  const auto& elemEndIt = gridView.template end</*codim=*/0, Dune::Interior_Partition>();
942 
943  for (auto elemIt = gridView.template begin</*codim=*/0, Dune::Interior_Partition>();
944  elemIt != elemEndIt; ++elemIt)
945  {
946  elemCtx.updatePrimaryStencil(*elemIt);
947  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
948 
949  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
950  const auto& fs = intQuants.fluidState();
951 
952  for ( int phaseIdx = 0; phaseIdx < np; ++phaseIdx )
953  {
954  auto& B = B_avg[ phaseIdx ];
955  const int ebosPhaseIdx = flowPhaseToEbosPhaseIdx(phaseIdx);
956 
957  B += 1 / fs.invB(ebosPhaseIdx).value();
958  }
959  if (has_solvent_) {
960  auto& B = B_avg[solventSaturationIdx];
961  B += 1 / intQuants.solventInverseFormationVolumeFactor().value();
962  }
963  }
964 
965  // compute global average
966  grid.comm().sum(B_avg.data(), B_avg.size());
967  for(auto& bval: B_avg)
968  {
969  bval/=global_nc_;
970  }
971  }
972 
973 
974 
975 
976 
977  template<typename TypeTag>
978  void
979  BlackoilWellModel<TypeTag>::
980  updatePrimaryVariables(const WellState& well_state) const
981  {
982  for (const auto& well : well_container_) {
983  well->updatePrimaryVariables(well_state);
984  }
985  }
986 
987 } // namespace Opm
bool localWellsActive() const
return true if wells are available on this process
Definition: BlackoilWellModel_impl.hpp:378
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:66
Utility class to handle the log messages about well switching.
Definition: WellSwitchingLogger.hpp:44
bool wellsActive() const
return true if wells are available in the reservoir
Definition: BlackoilWellModel_impl.hpp:354
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
void calculateExplicitQuantities(const Simulator &ebosSimulator, const WellState &xw) const
Calculating the explict quantities used in the well calculation.
Definition: BlackoilWellModel_impl.hpp:540
void updateListEconLimited(const Schedule &schedule, const int current_step, const Wells *wells_struct, const WellState &well_state, DynamicListEconLimited &list_econ_limited) const
upate the dynamic lists related to economic limits
Definition: BlackoilWellModel_impl.hpp:579