All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
SimulatorBase_impl.hpp
1 /*
2  Copyright 2013 SINTEF ICT, Applied Mathematics.
3  Copyright 2014-2016 IRIS AS
4  Copyright 2015 Andreas Lauser
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #include <utility>
23 #include <functional>
24 #include <algorithm>
25 #include <locale>
26 #include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
27 #include <opm/core/utility/initHydroCarbonState.hpp>
28 #include <opm/core/well_controls.h>
29 #include <opm/core/wells/DynamicListEconLimited.hpp>
30 #include <opm/autodiff/BlackoilModel.hpp>
31 
32 namespace Opm
33 {
34 
35  template <class Implementation>
36  SimulatorBase<Implementation>::SimulatorBase(const ParameterGroup& param,
37  const Grid& grid,
38  DerivedGeology& geo,
40  const RockCompressibility* rock_comp_props,
42  const double* gravity,
43  const bool has_disgas,
44  const bool has_vapoil,
45  std::shared_ptr<EclipseState> eclipse_state,
46  OutputWriter& output_writer,
47  const std::vector<double>& threshold_pressures_by_face,
48  const std::unordered_set<std::string>& defunct_well_names)
49  : param_(param),
50  model_param_(param),
51  solver_param_(param),
52  grid_(grid),
53  props_(props),
54  rock_comp_props_(rock_comp_props),
55  gravity_(gravity),
56  geo_(geo),
57  solver_(linsolver),
58  has_disgas_(has_disgas),
59  has_vapoil_(has_vapoil),
60  terminal_output_(param.getDefault("output_terminal", true)),
61  eclipse_state_(eclipse_state),
62  output_writer_(output_writer),
63  rateConverter_(props_.phaseUsage(), std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
64  threshold_pressures_by_face_(threshold_pressures_by_face),
65  is_parallel_run_( false ),
66  defunct_well_names_(defunct_well_names)
67  {
68  // Misc init.
69  const int num_cells = AutoDiffGrid::numCells(grid);
70  allcells_.resize(num_cells);
71  for (int cell = 0; cell < num_cells; ++cell) {
72  allcells_[cell] = cell;
73  }
74 #if HAVE_MPI
75  if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
76  {
77  const ParallelISTLInformation& info =
78  boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
79  // Only rank 0 does print to std::cout
80  terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
81  is_parallel_run_ = ( info.communicator().size() > 1 );
82  }
83 #endif
84  }
85 
86  template <class Implementation>
88  ReservoirState& state)
89  {
90  WellState prev_well_state;
91 
92  ExtraData extra;
93  if (output_writer_.isRestart()) {
94  // This is a restart, populate WellState and ReservoirState state objects from restart file
95  output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state, extra);
96  initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid_), has_disgas_, has_vapoil_);
97  initHysteresisParams(state);
98  }
99 
100  // Create timers and file for writing timing info.
101  Opm::time::StopWatch solver_timer;
102  Opm::time::StopWatch step_timer;
103  Opm::time::StopWatch total_timer;
104  total_timer.start();
105  std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
106  std::ofstream tstep_os;
107 
108  if ( output_writer_.output() ) {
109  if ( output_writer_.isIORank() )
110  tstep_os.open(tstep_filename.c_str());
111  }
112 
113  const auto& schedule = eclipse_state_->getSchedule();
114 
115  // adaptive time stepping
116  const auto& events = schedule.getEvents();
117  std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
118  if( param_.getDefault("timestep.adaptive", true ) )
119  {
120 
121  if (param_.getDefault("use_TUNING", false)) {
122  adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule.getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
123  } else {
124  adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
125  }
126  if (output_writer_.isRestart()) {
127  if (extra.suggested_step > 0.0) {
128  adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
129  }
130  }
131  }
132 
133  std::string restorefilename = param_.getDefault("restorefile", std::string("") );
134  if( ! restorefilename.empty() )
135  {
136  // -1 means that we'll take the last report step that was written
137  const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
138 
139  output_writer_.restore( timer,
140  state,
141  prev_well_state,
142  restorefilename,
143  desiredRestoreStep );
144  }
145 
146  DynamicListEconLimited dynamic_list_econ_limited;
147  SimulatorReport report;
148  SimulatorReport stepReport;
149 
150  bool ooip_computed = false;
151  std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
152  //Get compressed cell fipnum.
153  std::vector<int> fipnum(AutoDiffGrid::numCells(grid_));
154  if (fipnum_global.empty()) {
155  std::fill(fipnum.begin(), fipnum.end(), 0);
156  } else {
157  for (size_t c = 0; c < fipnum.size(); ++c) {
158  fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
159  }
160  }
161  std::vector<std::vector<double> > OOIP;
162  // Main simulation loop.
163  while (!timer.done()) {
164  // Report timestep.
165  step_timer.start();
166  if ( terminal_output_ )
167  {
168  std::ostringstream ss;
169  timer.report(ss);
170  OpmLog::note(ss.str());
171  }
172 
173  // Create wells and well state.
174  WellsManager wells_manager(*eclipse_state_,
175  timer.currentStepNum(),
176  Opm::UgGridHelpers::numCells(grid_),
177  Opm::UgGridHelpers::globalCell(grid_),
178  Opm::UgGridHelpers::cartDims(grid_),
179  Opm::UgGridHelpers::dimensions(grid_),
180  Opm::UgGridHelpers::cell2Faces(grid_),
181  Opm::UgGridHelpers::beginFaceCentroids(grid_),
182  dynamic_list_econ_limited,
183  is_parallel_run_,
184  defunct_well_names_);
185  const Wells* wells = wells_manager.c_wells();
186  WellState well_state;
187  well_state.init(wells, state, prev_well_state, props_.phaseUsage());
188 
189  // give the polymer and surfactant simulators the chance to do their stuff
190  asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
191 
192  // write the inital state at the report stage
193  if (timer.initialStep()) {
194  Dune::Timer perfTimer;
195  perfTimer.start();
196 
197  // No per cell data is written for initial step, but will be
198  // for subsequent steps, when we have started simulating
199  output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state, {}, {} );
200 
201  report.output_write_time += perfTimer.stop();
202  }
203 
204  // Max oil saturation (for VPPARS), hysteresis update.
205  props_.updateSatOilMax(state.saturation());
206  props_.updateSatHyst(state.saturation(), allcells_);
207 
208  // Compute reservoir volumes for RESV controls.
209  asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
210 
211  // Run a multiple steps of the solver depending on the time step control.
212  solver_timer.start();
213 
214  const WellModel well_model(wells, &(wells_manager.wellCollection()));
215 
216  std::unique_ptr<Solver> solver = asImpl().createSolver(well_model);
217 
218  // Compute orignal FIP;
219  if (!ooip_computed) {
220  OOIP = solver->computeFluidInPlace(state, fipnum);
221  FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
222  ooip_computed = true;
223  }
224 
225  if( terminal_output_ )
226  {
227  std::ostringstream step_msg;
228  boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
229  step_msg.imbue(std::locale(std::locale::classic(), facet));
230  step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
231  << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
232  << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
233  << ", date = " << timer.currentDateTime();
234  OpmLog::info(step_msg.str());
235  }
236 
237  // If sub stepping is enabled allow the solver to sub cycle
238  // in case the report steps are too large for the solver to converge
239  //
240  // \Note: The report steps are met in any case
241  // \Note: The sub stepping will require a copy of the state variables
242  if( adaptiveTimeStepping ) {
243  bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
244  events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
245  events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
246  events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
247  report += adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
248  output_writer_.requireFIPNUM() ? &fipnum : nullptr );
249  }
250  else {
251  // solve for complete report step
252  stepReport = solver->step(timer, state, well_state);
253  report += stepReport;
254 
255  if( terminal_output_ )
256  {
257  std::ostringstream iter_msg;
258  iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
259  if (solver->wellIterations() != 0) {
260  iter_msg << " days well iterations = " << solver->wellIterations() << ", ";
261  }
262  iter_msg << "non-linear iterations = " << solver->nonlinearIterations()
263  << ", total linear iterations = " << solver->linearIterations()
264  << "\n";
265  OpmLog::info(iter_msg.str());
266  }
267  }
268 
269  // update the derived geology (transmissibilities, pore volumes, etc) if the
270  // has geology changed for the next report step
271  const int nextTimeStepIdx = timer.currentStepNum() + 1;
272  if (nextTimeStepIdx < timer.numSteps()
273  && events.hasEvent(ScheduleEvents::GEO_MODIFIER, nextTimeStepIdx)) {
274  // bring the contents of the keywords to the current state of the SCHEDULE
275  // section
276  //
277  // TODO (?): handle the parallel case (maybe this works out of the box)
278  const auto& miniDeck = schedule.getModifierDeck(nextTimeStepIdx);
279  eclipse_state_->applyModifierDeck(miniDeck);
280  geo_.update(grid_, props_, *eclipse_state_, gravity_);
281  }
282 
283  // take time that was used to solve system for this reportStep
284  solver_timer.stop();
285 
286  // update timing.
287  report.solver_time += solver_timer.secsSinceStart();
288 
289  // Compute current FIP.
290  std::vector<std::vector<double> > COIP;
291  COIP = solver->computeFluidInPlace(state, fipnum);
292  std::vector<double> OOIP_totals = FIPTotals(OOIP, state);
293  std::vector<double> COIP_totals = FIPTotals(COIP, state);
294 
295  //Convert to correct units
296  FIPUnitConvert(eclipse_state_->getUnits(), COIP);
297  FIPUnitConvert(eclipse_state_->getUnits(), OOIP_totals);
298  FIPUnitConvert(eclipse_state_->getUnits(), COIP_totals);
299 
300  if ( terminal_output_ )
301  {
302  outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
303  for (size_t reg = 0; reg < OOIP.size(); ++reg) {
304  outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
305  }
306  }
307 
308  if ( terminal_output_ )
309  {
310  std::string msg;
311  msg = "Fully implicit solver took: " + std::to_string(stepReport.solver_time) + " seconds. Total solver time taken: " + std::to_string(report.solver_time) + " seconds.";
312  OpmLog::note(msg);
313  }
314 
315  if ( tstep_os.is_open() ) {
316  stepReport.reportParam(tstep_os);
317  }
318 
319  // Increment timer, remember well state.
320  ++timer;
321 
322  // write simulation state at the report stage
323  Dune::Timer perfTimer;
324  perfTimer.start();
325  const auto& physicalModel = solver->model();
326  output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
327  report.output_write_time += perfTimer.stop();
328 
329  prev_well_state = well_state;
330 
331  asImpl().updateListEconLimited(solver, eclipse_state_->getSchedule(), timer.currentStepNum(), wells,
332  well_state, dynamic_list_econ_limited);
333  }
334 
335  // Stop timer and create timing report
336  total_timer.stop();
337  report.total_time = total_timer.secsSinceStart();
338  report.converged = true;
339  return report;
340  }
341 
342  namespace SimFIBODetails {
343  typedef std::unordered_map<std::string, const Well* > WellMap;
344 
345  inline WellMap
346  mapWells(const std::vector< const Well* >& wells)
347  {
348  WellMap wmap;
349 
350  for (std::vector< const Well* >::const_iterator
351  w = wells.begin(), e = wells.end();
352  w != e; ++w)
353  {
354  wmap.insert(std::make_pair((*w)->name(), *w));
355  }
356 
357  return wmap;
358  }
359 
360  inline int
361  resv_control(const WellControls* ctrl)
362  {
363  int i, n = well_controls_get_num(ctrl);
364 
365  bool match = false;
366  for (i = 0; (! match) && (i < n); ++i) {
367  match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
368  }
369 
370  if (! match) { i = 0; }
371 
372  return i - 1; // -1 if no match, undo final "++" otherwise
373  }
374 
375  inline bool
376  is_resv(const Wells& wells,
377  const int w)
378  {
379  return (0 <= resv_control(wells.ctrls[w]));
380  }
381 
382  inline bool
383  is_resv(const WellMap& wmap,
384  const std::string& name,
385  const std::size_t step)
386  {
387  bool match = false;
388 
389  WellMap::const_iterator i = wmap.find(name);
390 
391  if (i != wmap.end()) {
392  const Well* wp = i->second;
393 
394  match = (wp->isProducer(step) &&
395  wp->getProductionProperties(step)
396  .hasProductionControl(WellProducer::RESV))
397  || (wp->isInjector(step) &&
398  wp->getInjectionProperties(step)
399  .hasInjectionControl(WellInjector::RESV));
400  }
401 
402  return match;
403  }
404 
405  inline std::vector<int>
406  resvWells(const Wells* wells,
407  const std::size_t step,
408  const WellMap& wmap)
409  {
410  std::vector<int> resv_wells;
411  if( wells )
412  {
413  for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
414  if (is_resv(*wells, w) ||
415  ((wells->name[w] != 0) &&
416  is_resv(wmap, wells->name[w], step)))
417  {
418  resv_wells.push_back(w);
419  }
420  }
421  }
422 
423  return resv_wells;
424  }
425 
426  inline void
427  historyRates(const PhaseUsage& pu,
428  const WellProductionProperties& p,
429  std::vector<double>& rates)
430  {
431  assert (! p.predictionMode);
432  assert (rates.size() ==
433  std::vector<double>::size_type(pu.num_phases));
434 
435  if (pu.phase_used[ BlackoilPhases::Aqua ]) {
436  const std::vector<double>::size_type
437  i = pu.phase_pos[ BlackoilPhases::Aqua ];
438 
439  rates[i] = p.WaterRate;
440  }
441 
442  if (pu.phase_used[ BlackoilPhases::Liquid ]) {
443  const std::vector<double>::size_type
444  i = pu.phase_pos[ BlackoilPhases::Liquid ];
445 
446  rates[i] = p.OilRate;
447  }
448 
449  if (pu.phase_used[ BlackoilPhases::Vapour ]) {
450  const std::vector<double>::size_type
451  i = pu.phase_pos[ BlackoilPhases::Vapour ];
452 
453  rates[i] = p.GasRate;
454  }
455  }
456  } // namespace SimFIBODetails
457 
458  template <class Implementation>
459  void SimulatorBase<Implementation>::handleAdditionalWellInflow(SimulatorTimer& /* timer */,
460  WellsManager& /* wells_manager */,
461  WellState& /* well_state */,
462  const Wells* /* wells */)
463  { }
464 
465  template <class Implementation>
466  auto SimulatorBase<Implementation>::createSolver(const WellModel& well_model)
467  -> std::unique_ptr<Solver>
468  {
469  auto model = std::unique_ptr<Model>(new Model(model_param_,
470  grid_,
471  props_,
472  geo_,
473  rock_comp_props_,
474  well_model,
475  solver_,
476  eclipse_state_,
477  has_disgas_,
478  has_vapoil_,
479  terminal_output_));
480 
481  if (!threshold_pressures_by_face_.empty()) {
482  model->setThresholdPressures(threshold_pressures_by_face_);
483  }
484 
485  return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
486  }
487 
488 
489  template <class Implementation>
490  void SimulatorBase<Implementation>::computeRESV(const std::size_t step,
491  const Wells* wells,
492  const BlackoilState& x,
493  WellState& xw)
494  {
495  typedef SimFIBODetails::WellMap WellMap;
496 
497  const auto w_ecl = eclipse_state_->getSchedule().getWells(step);
498  const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
499 
500  const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
501 
502  const std::size_t number_resv_wells = resv_wells.size();
503  std::size_t global_number_resv_wells = number_resv_wells;
504 #if HAVE_MPI
505  if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
506  {
507  const auto& info =
508  boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
509  global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
510  if ( global_number_resv_wells )
511  {
512  // At least one process has resv wells. Therefore rate converter needs
513  // to calculate averages over regions that might cross process
514  // borders. This needs to be done by all processes and therefore
515  // outside of the next if statement.
516  rateConverter_.defineState(x, boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()));
517  }
518  }
519  else
520 #endif
521  {
522  if ( global_number_resv_wells )
523  {
524  rateConverter_.defineState(x);
525  }
526  }
527 
528  if (! resv_wells.empty()) {
529  const PhaseUsage& pu = props_.phaseUsage();
530  const std::vector<double>::size_type np = props_.numPhases();
531 
532  std::vector<double> distr (np);
533  std::vector<double> hrates(np);
534  std::vector<double> prates(np);
535 
536  for (std::vector<int>::const_iterator
537  rp = resv_wells.begin(), e = resv_wells.end();
538  rp != e; ++rp)
539  {
540  WellControls* ctrl = wells->ctrls[*rp];
541  const bool is_producer = wells->type[*rp] == PRODUCER;
542 
543  // RESV control mode, all wells
544  {
545  const int rctrl = SimFIBODetails::resv_control(ctrl);
546 
547  if (0 <= rctrl) {
548  const std::vector<double>::size_type off = (*rp) * np;
549 
550  if (is_producer) {
551  // Convert to positive rates to avoid issues
552  // in coefficient calculations.
553  std::transform(xw.wellRates().begin() + (off + 0*np),
554  xw.wellRates().begin() + (off + 1*np),
555  prates.begin(), std::negate<double>());
556  } else {
557  std::copy(xw.wellRates().begin() + (off + 0*np),
558  xw.wellRates().begin() + (off + 1*np),
559  prates.begin());
560  }
561 
562  const int fipreg = 0; // Hack. Ignore FIP regions.
563  const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
564  const int pvtreg = props_.cellPvtRegionIndex()[well_cell_top];
565  rateConverter_.calcCoeff(fipreg, pvtreg, distr);
566 
567  well_controls_iset_distr(ctrl, rctrl, & distr[0]);
568  }
569  }
570 
571  // RESV control, WCONHIST wells. A bit of duplicate
572  // work, regrettably.
573  if (is_producer && wells->name[*rp] != 0) {
574  WellMap::const_iterator i = wmap.find(wells->name[*rp]);
575 
576  if (i != wmap.end()) {
577  const auto* wp = i->second;
578 
579  const WellProductionProperties& p =
580  wp->getProductionProperties(step);
581 
582  if (! p.predictionMode) {
583  // History matching (WCONHIST/RESV)
584  SimFIBODetails::historyRates(pu, p, hrates);
585 
586  const int fipreg = 0; // Hack. Ignore FIP regions.
587  const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
588  const int pvtreg = props_.cellPvtRegionIndex()[well_cell_top];
589  rateConverter_.calcCoeff(fipreg, pvtreg, distr);
590 
591  // WCONHIST/RESV target is sum of all
592  // observed phase rates translated to
593  // reservoir conditions. Recall sign
594  // convention: Negative for producers.
595  const double target =
596  - std::inner_product(distr.begin(), distr.end(),
597  hrates.begin(), 0.0);
598 
599  well_controls_clear(ctrl);
600  well_controls_assert_number_of_phases(ctrl, int(np));
601 
602  static const double invalid_alq = -std::numeric_limits<double>::max();
603  static const int invalid_vfp = -std::numeric_limits<int>::max();
604 
605  const int ok_resv =
606  well_controls_add_new(RESERVOIR_RATE, target,
607  invalid_alq, invalid_vfp,
608  & distr[0], ctrl);
609 
610  // For WCONHIST the BHP limit is set to 1 atm.
611  // or a value specified using WELTARG
612  double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
613  const int ok_bhp =
614  well_controls_add_new(BHP, bhp_limit,
615  invalid_alq, invalid_vfp,
616  NULL, ctrl);
617 
618  if (ok_resv != 0 && ok_bhp != 0) {
619  xw.currentControls()[*rp] = 0;
620  well_controls_set_current(ctrl, 0);
621  }
622  }
623  }
624  }
625  }
626  }
627 
628  if( wells )
629  {
630  for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
631  WellControls* ctrl = wells->ctrls[w];
632  const bool is_producer = wells->type[w] == PRODUCER;
633  if (!is_producer && wells->name[w] != 0) {
634  WellMap::const_iterator i = wmap.find(wells->name[w]);
635  if (i != wmap.end()) {
636  const auto* wp = i->second;
637  const WellInjectionProperties& injector = wp->getInjectionProperties(step);
638  if (!injector.predictionMode) {
639  //History matching WCONINJEH
640  static const double invalid_alq = -std::numeric_limits<double>::max();
641  static const int invalid_vfp = -std::numeric_limits<int>::max();
642  // For WCONINJEH the BHP limit is set to a large number
643  // or a value specified using WELTARG
644  double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
645  const int ok_bhp =
646  well_controls_add_new(BHP, bhp_limit,
647  invalid_alq, invalid_vfp,
648  NULL, ctrl);
649  if (!ok_bhp) {
650  OPM_THROW(std::runtime_error, "Failed to add well control.");
651  }
652  }
653  }
654  }
655  }
656  }
657  }
658 
659  template <class Implementation>
660  void
661  SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units,
662  std::vector<std::vector<double> >& fip)
663  {
664  for (size_t i = 0; i < fip.size(); ++i) {
665  FIPUnitConvert(units, fip[i]);
666  }
667  }
668 
669  template <class Implementation>
670  void
671  SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units, std::vector<double>& fip)
672  {
673  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
674  fip[0] = unit::convert::to(fip[0], unit::stb);
675  fip[1] = unit::convert::to(fip[1], unit::stb);
676  fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
677  fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
678  fip[4] = unit::convert::to(fip[4], unit::stb);
679  fip[5] = unit::convert::to(fip[5], unit::stb);
680  fip[6] = unit::convert::to(fip[6], unit::psia);
681  }
682  else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
683  fip[6] = unit::convert::to(fip[6], unit::barsa);
684  }
685  else {
686  OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
687  }
688  }
689 
690 
691  template <class Implementation>
692  std::vector<double>
693  SimulatorBase<Implementation>::FIPTotals(const std::vector<std::vector<double> >& fip, const ReservoirState& state)
694  {
695  std::vector<double> totals(7, 0.0);
696  for (int i = 0; i < 5; ++i) {
697  for (size_t reg = 0; reg < fip.size(); ++reg) {
698  totals[i] += fip[reg][i];
699  }
700  }
701  const int nc = Opm::AutoDiffGrid::numCells(grid_);
702  const int np = state.numPhases();
703  const PhaseUsage& pu = props_.phaseUsage();
704  const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
705  std::vector<double> so(nc);
706  std::vector<double> sg(nc);
707  std::vector<double> hydrocarbon(nc);
708  // Using dummy indices if phase not used, the columns will not be accessed below if unused.
709  const int oilpos = pu.phase_used[BlackoilPhases::Liquid] ? pu.phase_pos[BlackoilPhases::Liquid] : 0;
710  const int gaspos = pu.phase_used[BlackoilPhases::Vapour] ? pu.phase_pos[BlackoilPhases::Vapour] : 0;
711  const auto& soCol = s.col(oilpos);
712  const auto& sgCol = s.col(gaspos);
713  for (unsigned c = 0; c < so.size(); ++ c) {
714  double mySo = 0.0;
715  if (pu.phase_used[BlackoilPhases::Liquid]) {
716  mySo = soCol[c];
717  }
718  double mySg = 0.0;
719  if (pu.phase_used[BlackoilPhases::Vapour]) {
720  mySg = sgCol[c];
721  }
722  so[c] = mySo;
723  sg[c] = mySg;
724  hydrocarbon[c] = mySo + mySg;
725  }
726  const std::vector<double> p = state.pressure();
727  if ( ! is_parallel_run_ )
728  {
729  double tmp = 0.0;
730  double tmp2 = 0.0;
731  for (unsigned i = 0; i < p.size(); ++i) {
732  tmp += p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
733  tmp2 += geo_.poreVolume()[i] * hydrocarbon[i];
734  }
735  totals[5] = geo_.poreVolume().sum();
736  totals[6] = tmp/tmp2;
737  }
738  else
739  {
740 #if HAVE_MPI
741  const auto & pinfo =
742  boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
743  auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
744  Opm::Reduction::makeGlobalSumFunctor<double>(),
745  Opm::Reduction::makeGlobalSumFunctor<double>());
746  std::vector<double> pav_nom(p.size());
747  std::vector<double> pav_denom(pav_nom.size());
748  for (unsigned i = 0; i < p.size(); ++i) {
749  pav_nom[i] = p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
750  pav_denom[i] = geo_.poreVolume()[i] * hydrocarbon[i];
751  }
752 
753  // using ref cref to prevent copying
754  auto inputs = std::make_tuple(std::cref(geo_.poreVolume()),
755  std::cref(pav_nom), std::cref(pav_denom));
756  std::tuple<double, double, double> results(0.0, 0.0, 0.0);
757 
758  pinfo.computeReduction(inputs, operators, results);
759  using std::get;
760  totals[5] = get<0>(results);
761  totals[6] = get<1>(results)/get<2>(results);
762 
763 #else
764  // This should never happen!
765  OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel");
766 #endif
767  }
768 
769  return totals;
770  }
771 
772 
773 
774  template <class Implementation>
775  void
776  SimulatorBase<Implementation>::outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
777  {
778  std::ostringstream ss;
779  if (!reg) {
780  ss << " ===================================================\n"
781  << " : Field Totals :\n";
782  } else {
783  ss << " ===================================================\n"
784  << " : FIPNUM report region "
785  << std::setw(2) << reg << " :\n";
786  }
787  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
788  ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
789  << std::fixed << std::setprecision(0)
790  << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
791  if (!reg) {
792  ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
793  << " : Porv volumes are taken at reference conditions :\n";
794  }
795  ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
796  }
797  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
798  ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
799  << std::fixed << std::setprecision(0)
800  << " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
801  if (!reg) {
802  ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
803  << " : Pore volumes are taken at reference conditions :\n";
804  }
805  ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
806  }
807  ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
808  << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
809  << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
810  << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
811  << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
812  << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
813  << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
814  << ":========================:==========================================:================:==========================================:\n";
815  OpmLog::note(ss.str());
816  }
817 
818 
819  template <class Implementation>
820  void
821  SimulatorBase<Implementation>::
822  updateListEconLimited(const std::unique_ptr<Solver>& solver,
823  const Schedule& schedule,
824  const int current_step,
825  const Wells* wells,
826  const WellState& well_state,
827  DynamicListEconLimited& list_econ_limited) const
828  {
829 
830  solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
831  well_state, list_econ_limited);
832  }
833 
834  template <class Implementation>
835  void
836  SimulatorBase<Implementation>::
837  initHysteresisParams(ReservoirState& state)
838  {
839  typedef std::vector<double> VectorType;
840 
841  const VectorType& somax = state.getCellData( "SOMAX" );
842 
843  VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
844  VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
845 
846  VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
847  VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
848 
849  props_.setSatOilMax(somax);
850  props_.setOilWaterHystParams(pcSwMdc_ow, krnSwMdc_ow, allcells_);
851  props_.setGasOilHystParams(pcSwMdc_go, krnSwMdc_go, allcells_);
852  }
853 
854 
855 } // namespace Opm
Extra data to read/write for OPM restarting.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:199
int currentStepNum() const
Current step number.
Definition: SimulatorTimer.cpp:77
virtual const boost::any & parallelInformation() const =0
Get the information about the parallelization of the grid.
double totalTime() const
Total time.
Definition: SimulatorTimer.cpp:121
bool done() const
Return true if op++() has been called numSteps() times.
Definition: SimulatorTimer.cpp:155
Definition: AdaptiveTimeStepping.hpp:38
double simulationTimeElapsed() const
Time elapsed since the start of the simulation until the beginning of the current time step [s]...
Definition: SimulatorTimer.cpp:104
SimulatorBase(const ParameterGroup &param, const Grid &grid, DerivedGeology &geo, BlackoilPropsAdFromDeck &props, const RockCompressibility *rock_comp_props, NewtonIterationBlackoilInterface &linsolver, const double *gravity, const bool disgas, const bool vapoil, std::shared_ptr< EclipseState > eclipse_state, OutputWriter &output_writer, const std::vector< double > &threshold_pressures_by_face, const std::unordered_set< std::string > &defunct_well_names)
Initialise from parameters and objects to observe.
Definition: SimulatorBase_impl.hpp:36
double currentStepLength() const
Current step length.
Definition: SimulatorTimer.cpp:91
int numSteps() const
Total number of steps.
Definition: SimulatorTimer.cpp:71
boost::posix_time::ptime currentDateTime() const
Return current date.
Definition: SimulatorTimer.cpp:115
bool initialStep() const
Whether the current step is the first step.
Definition: SimulatorTimer.cpp:65
void report(std::ostream &os) const
Print a report with current and total time etc.
Definition: SimulatorTimer.cpp:136
Class containing static geological properties that are derived from grid and petrophysical properties...
Definition: GeoProps.hpp:59
SimulatorReport run(SimulatorTimer &timer, ReservoirState &state)
Run the simulation.
Definition: SimulatorBase_impl.hpp:87
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
This class implements the AD-adapted fluid interface for three-phase black-oil.
Definition: BlackoilPropsAdFromDeck.hpp:61
Definition: SimulatorTimer.hpp:34