All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
SimulatorFullyImplicitBlackoilEbos.hpp
1 /*
2  Copyright 2013, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 Andreas Lauser
4 
5  This file is part of the Open Porous Media project (OPM).
6 
7  OPM is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  OPM is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with OPM. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
22 #define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
23 
24 #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
25 #include <opm/autodiff/IterationReport.hpp>
26 #include <opm/autodiff/NonlinearSolver.hpp>
27 #include <opm/autodiff/BlackoilModelEbos.hpp>
28 #include <opm/autodiff/BlackoilModelParameters.hpp>
29 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
30 #include <opm/autodiff/BlackoilWellModel.hpp>
32 #include <opm/autodiff/SimFIBODetails.hpp>
33 #include <opm/autodiff/moduleVersion.hpp>
34 #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
35 #include <opm/core/utility/initHydroCarbonState.hpp>
36 #include <opm/core/utility/StopWatch.hpp>
37 
38 #include <opm/common/Exceptions.hpp>
39 #include <opm/common/ErrorMacros.hpp>
40 
41 #include <dune/common/unused.hh>
42 
43 namespace Opm {
44 
45 
47 template<class TypeTag>
49 {
50 public:
51  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
52  typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
53  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
54  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
55  typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
56  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
57  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
58  typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
59  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
60 
61  typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
62 
64  typedef BlackoilState ReservoirState;
71 
72 
96  SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
97  const ParameterGroup& param,
99  const bool has_disgas,
100  const bool has_vapoil,
101  const EclipseState& /* eclState */,
102  OutputWriter& output_writer,
103  const std::unordered_set<std::string>& defunct_well_names)
104  : ebosSimulator_(ebosSimulator),
105  param_(param),
106  model_param_(param),
107  solver_param_(param),
108  solver_(linsolver),
109  phaseUsage_(phaseUsageFromDeck(eclState())),
110  has_disgas_(has_disgas),
111  has_vapoil_(has_vapoil),
112  terminal_output_(param.getDefault("output_terminal", true)),
113  output_writer_(output_writer),
114  rateConverter_(createRateConverter_()),
115  defunct_well_names_( defunct_well_names ),
116  is_parallel_run_( false )
117  {
118 #if HAVE_MPI
119  if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
120  {
121  const ParallelISTLInformation& info =
122  boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
123  // Only rank 0 does print to std::cout
124  terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
125  is_parallel_run_ = ( info.communicator().size() > 1 );
126  }
127 #endif
128  }
129 
136  SimulatorReport run(SimulatorTimer& timer,
137  ReservoirState& state)
138  {
139  WellState prev_well_state;
140 
141  ExtraData extra;
142 
143  failureReport_ = SimulatorReport();
144  extractLegacyDepth_();
145 
146  // communicate the initial solution to ebos
147  if (timer.initialStep()) {
148  convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
149  ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
150  }
151 
152  if (output_writer_.isRestart()) {
153  // This is a restart, populate WellState and ReservoirState state objects from restart file
154  output_writer_.initFromRestartFile(phaseUsage_, grid(), state, prev_well_state, extra);
155  initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
156  initHysteresisParams(state);
157  // communicate the restart solution to ebos
158  convertInput(/*iterationIdx=*/0, state, ebosSimulator_ );
159  ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
160  }
161 
162  // Sync the overlap region of the inital solution. It was generated
163  // from the ReservoirState which has wrong values in the ghost region
164  // for some models (SPE9, Norne, Model 2)
165  ebosSimulator_.model().syncOverlap();
166 
167  // Create timers and file for writing timing info.
168  Opm::time::StopWatch solver_timer;
169  Opm::time::StopWatch step_timer;
170  Opm::time::StopWatch total_timer;
171  total_timer.start();
172  std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
173  std::ofstream tstep_os;
174 
175  if ( output_writer_.output() && output_writer_.isIORank() )
176  {
177  tstep_os.open(tstep_filename.c_str());
178  }
179 
180  const auto& schedule = eclState().getSchedule();
181 
182  // adaptive time stepping
183  const auto& events = schedule.getEvents();
184  std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
185  if( param_.getDefault("timestep.adaptive", true ) )
186  {
187 
188  if (param_.getDefault("use_TUNING", false)) {
189  adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule.getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
190  } else {
191  adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
192  }
193 
194  if (output_writer_.isRestart()) {
195  if (extra.suggested_step > 0.0) {
196  adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
197  }
198  }
199  }
200 
201  std::string restorefilename = param_.getDefault("restorefile", std::string("") );
202  if( ! restorefilename.empty() )
203  {
204  // -1 means that we'll take the last report step that was written
205  const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
206 
207  output_writer_.restore( timer,
208  state,
209  prev_well_state,
210  restorefilename,
211  desiredRestoreStep );
212  initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
213  initHysteresisParams(state);
214  // communicate the restart solution to ebos
215  convertInput(0, state, ebosSimulator_);
216  ebosSimulator_.model().invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
217  }
218 
219  DynamicListEconLimited dynamic_list_econ_limited;
220  SimulatorReport report;
221  SimulatorReport stepReport;
222 
223  std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
224  //Get compressed cell fipnum.
225  std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
226  if (fipnum_global.empty()) {
227  std::fill(fipnum.begin(), fipnum.end(), 0);
228  } else {
229  for (size_t c = 0; c < fipnum.size(); ++c) {
230  fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
231  }
232  }
233  std::vector<std::vector<double>> originalFluidInPlace;
234  std::vector<double> originalFluidInPlaceTotals;
235 
236  // Main simulation loop.
237  while (!timer.done()) {
238  // Report timestep.
239  step_timer.start();
240  if ( terminal_output_ )
241  {
242  std::ostringstream ss;
243  timer.report(ss);
244  OpmLog::note(ss.str());
245  }
246 
247  // Create wells and well state.
248  WellsManager wells_manager(eclState(),
249  timer.currentStepNum(),
250  Opm::UgGridHelpers::numCells(grid()),
251  Opm::UgGridHelpers::globalCell(grid()),
252  Opm::UgGridHelpers::cartDims(grid()),
253  Opm::UgGridHelpers::dimensions(grid()),
254  Opm::UgGridHelpers::cell2Faces(grid()),
255  Opm::UgGridHelpers::beginFaceCentroids(grid()),
256  dynamic_list_econ_limited,
257  is_parallel_run_,
258  defunct_well_names_ );
259  const Wells* wells = wells_manager.c_wells();
260  WellState well_state;
261 
262  // The well state initialize bhp with the cell pressure in the top cell.
263  // We must therefore provide it with updated cell pressures
264  size_t nc = Opm::UgGridHelpers::numCells(grid());
265  std::vector<double> cellPressures(nc, 0.0);
266  const auto& gridView = ebosSimulator_.gridManager().gridView();
267  ElementContext elemCtx(ebosSimulator_);
268  const auto& elemEndIt = gridView.template end</*codim=*/0>();
269  for (auto elemIt = gridView.template begin</*codim=*/0>();
270  elemIt != elemEndIt;
271  ++elemIt)
272  {
273  const auto& elem = *elemIt;
274  if (elem.partitionType() != Dune::InteriorEntity) {
275  continue;
276  }
277 
278  elemCtx.updatePrimaryStencil(elem);
279  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
280 
281  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
282  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
283  const auto& fs = intQuants.fluidState();
284 
285  const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
286  cellPressures[cellIdx] = p;
287  }
288  well_state.init(wells, cellPressures, prev_well_state, phaseUsage_);
289 
290  // give the polymer and surfactant simulators the chance to do their stuff
291  handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
292 
293  // Compute reservoir volumes for RESV controls.
294  computeRESV(timer.currentStepNum(), wells, well_state);
295 
296  // Run a multiple steps of the solver depending on the time step control.
297  solver_timer.start();
298 
299  const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum());
300  extractLegacyCellPvtRegionIndex_();
301  WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_,
302  rateConverter_, terminal_output_, timer.currentStepNum(), legacyCellPvtRegionIdx_);
303 
304  // handling MS well related
305  if (model_param_.use_multisegment_well_) { // if we use MultisegmentWell model
306  for (const auto& well : wells_ecl) {
307  // TODO: this is acutally not very accurate, because sometimes a deck just claims a MS well
308  // while keep the well shut. More accurately, we should check if the well exisits in the Wells
309  // structure here
310  if (well->isMultiSegment(timer.currentStepNum()) ) { // there is one well is MS well
311  well_state.initWellStateMSWell(wells, wells_ecl, timer.currentStepNum(), phaseUsage_, prev_well_state);
312  break;
313  }
314  }
315  }
316 
317 
318  auto solver = createSolver(well_model);
319 
320  std::vector<std::vector<double>> currentFluidInPlace;
321  std::vector<double> currentFluidInPlaceTotals;
322 
323  // Compute orignal fluid in place if this has not been done yet
324  if (originalFluidInPlace.empty()) {
325  originalFluidInPlace = solver->computeFluidInPlace(fipnum);
326  originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace);
327  FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
328  FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
329 
330  currentFluidInPlace = originalFluidInPlace;
331  currentFluidInPlaceTotals = originalFluidInPlaceTotals;
332  }
333 
334  // write the inital state at the report stage
335  if (timer.initialStep()) {
336  Dune::Timer perfTimer;
337  perfTimer.start();
338 
339  if (terminal_output_) {
340  outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
341  for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
342  outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
343  }
344  }
345 
346  // No per cell data is written for initial step, but will be
347  // for subsequent steps, when we have started simulating
348  output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
349 
350  report.output_write_time += perfTimer.stop();
351  }
352 
353  if( terminal_output_ )
354  {
355  std::ostringstream step_msg;
356  boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
357  step_msg.imbue(std::locale(std::locale::classic(), facet));
358  step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
359  << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
360  << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
361  << ", date = " << timer.currentDateTime();
362  OpmLog::info(step_msg.str());
363  }
364 
365  solver->model().beginReportStep();
366 
367  // If sub stepping is enabled allow the solver to sub cycle
368  // in case the report steps are too large for the solver to converge
369  //
370  // \Note: The report steps are met in any case
371  // \Note: The sub stepping will require a copy of the state variables
372  if( adaptiveTimeStepping ) {
373  bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
374  events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
375  events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
376  events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
377  stepReport = adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
378  output_writer_.requireFIPNUM() ? &fipnum : nullptr );
379  report += stepReport;
380  failureReport_ += adaptiveTimeStepping->failureReport();
381  }
382  else {
383  // solve for complete report step
384  stepReport = solver->step(timer, state, well_state);
385  report += stepReport;
386  failureReport_ += solver->failureReport();
387 
388  if( terminal_output_ )
389  {
390  //stepReport.briefReport();
391  std::ostringstream iter_msg;
392  iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
393  if (solver->wellIterations() != 0) {
394  iter_msg << " days well iterations = " << solver->wellIterations() << ", ";
395  }
396  iter_msg << "non-linear iterations = " << solver->nonlinearIterations()
397  << ", total linear iterations = " << solver->linearIterations()
398  << "\n";
399  OpmLog::info(iter_msg.str());
400  }
401  }
402 
403  solver->model().endReportStep();
404 
405  // take time that was used to solve system for this reportStep
406  solver_timer.stop();
407 
408  // update timing.
409  report.solver_time += solver_timer.secsSinceStart();
410 
411  if ( output_writer_.output() && output_writer_.isIORank() )
412  {
413  stepReport.reportParam(tstep_os);
414  }
415 
416  // We don't need the reservoir state anymore. It is just passed around to avoid
417  // code duplication. Pass empty state instead.
418  if (timer.initialStep()) {
419  ReservoirState stateTrivial(0,0,0);
420  state = stateTrivial;
421  }
422 
423  // Increment timer, remember well state.
424  ++timer;
425 
426  // Compute current fluid in place.
427  currentFluidInPlace = solver->computeFluidInPlace(fipnum);
428  currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace);
429 
430  const std::string version = moduleVersionName();
431 
432  FIPUnitConvert(eclState().getUnits(), currentFluidInPlace);
433  FIPUnitConvert(eclState().getUnits(), currentFluidInPlaceTotals);
434 
435  if (terminal_output_ )
436  {
437  outputTimestampFIP(timer, version);
438  outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
439  for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
440  outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
441  }
442 
443  std::string msg;
444  msg =
445  "Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; "
446  "total solver time " + std::to_string(report.solver_time) + " seconds.";
447  OpmLog::note(msg);
448  }
449 
450  // write simulation state at the report stage
451  Dune::Timer perfTimer;
452  perfTimer.start();
453  const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
454  output_writer_.writeTimeStep( timer, state, well_state, solver->model(), false, nextstep, report);
455  report.output_write_time += perfTimer.stop();
456 
457  prev_well_state = well_state;
458 
459  updateListEconLimited(solver, eclState().getSchedule(), timer.currentStepNum(), wells,
460  well_state, dynamic_list_econ_limited);
461  }
462 
463  // Stop timer and create timing report
464  total_timer.stop();
465  report.total_time = total_timer.secsSinceStart();
466  report.converged = true;
467  return report;
468  }
469 
472  const SimulatorReport& failureReport() const { return failureReport_; };
473 
474  const Grid& grid() const
475  { return ebosSimulator_.gridManager().grid(); }
476 
477 protected:
478  void handleAdditionalWellInflow(SimulatorTimer& /*timer*/,
479  WellsManager& /* wells_manager */,
480  WellState& /* well_state */,
481  const Wells* /* wells */)
482  {
483  }
484 
485  std::unique_ptr<Solver> createSolver(WellModel& well_model)
486  {
487  const auto& gridView = ebosSimulator_.gridView();
488  const PhaseUsage& phaseUsage = phaseUsage_;
489  const std::vector<bool> activePhases = detail::activePhases(phaseUsage);
490  const double gravity = ebosSimulator_.problem().gravity()[2];
491 
492  // calculate the number of elements of the compressed sequential grid. this needs
493  // to be done in two steps because the dune communicator expects a reference as
494  // argument for sum()
495  int globalNumCells = gridView.size(/*codim=*/0);
496  globalNumCells = gridView.comm().sum(globalNumCells);
497 
498  well_model.init(phaseUsage,
499  activePhases,
500  gravity,
501  legacyDepth_,
502  globalNumCells,
503  grid());
504  auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
505  model_param_,
506  well_model,
507  rateConverter_,
508  solver_,
509  terminal_output_));
510 
511  return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
512  }
513 
514  void computeRESV(const std::size_t step,
515  const Wells* wells,
516  WellState& xw)
517  {
518  typedef SimFIBODetails::WellMap WellMap;
519 
520  const auto w_ecl = eclState().getSchedule().getWells(step);
521  const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
522 
523  const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
524 
525  const std::size_t number_resv_wells = resv_wells.size();
526  std::size_t global_number_resv_wells = number_resv_wells;
527 #if HAVE_MPI
528  if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
529  {
530  const auto& info =
531  boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
532  global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
533  if ( global_number_resv_wells )
534  {
535  // At least one process has resv wells. Therefore rate converter needs
536  // to calculate averages over regions that might cross process
537  // borders. This needs to be done by all processes and therefore
538  // outside of the next if statement.
539  rateConverter_.template defineState<ElementContext>(ebosSimulator_);
540  }
541  }
542  else
543 #endif
544  {
545  if ( global_number_resv_wells )
546  {
547  rateConverter_.template defineState<ElementContext>(ebosSimulator_);
548  }
549  }
550 
551  if (! resv_wells.empty()) {
552  const PhaseUsage& pu = phaseUsage_;
553  const std::vector<double>::size_type np = phaseUsage_.num_phases;
554 
555  std::vector<double> distr (np);
556  std::vector<double> hrates(np);
557  std::vector<double> prates(np);
558 
559  for (std::vector<int>::const_iterator
560  rp = resv_wells.begin(), e = resv_wells.end();
561  rp != e; ++rp)
562  {
563  WellControls* ctrl = wells->ctrls[*rp];
564  const bool is_producer = wells->type[*rp] == PRODUCER;
565  const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
566  const auto& eclProblem = ebosSimulator_.problem();
567  const int pvtreg = eclProblem.pvtRegionIndex(well_cell_top);
568 
569  // RESV control mode, all wells
570  {
571  const int rctrl = SimFIBODetails::resv_control(ctrl);
572 
573  if (0 <= rctrl) {
574  const std::vector<double>::size_type off = (*rp) * np;
575 
576  if (is_producer) {
577  // Convert to positive rates to avoid issues
578  // in coefficient calculations.
579  std::transform(xw.wellRates().begin() + (off + 0*np),
580  xw.wellRates().begin() + (off + 1*np),
581  prates.begin(), std::negate<double>());
582  } else {
583  std::copy(xw.wellRates().begin() + (off + 0*np),
584  xw.wellRates().begin() + (off + 1*np),
585  prates.begin());
586  }
587 
588  const int fipreg = 0; // Hack. Ignore FIP regions.
589  rateConverter_.calcCoeff(fipreg, pvtreg, distr);
590 
591  well_controls_iset_distr(ctrl, rctrl, & distr[0]);
592  }
593  }
594 
595  // RESV control, WCONHIST wells. A bit of duplicate
596  // work, regrettably.
597  if (is_producer && wells->name[*rp] != 0) {
598  WellMap::const_iterator i = wmap.find(wells->name[*rp]);
599 
600  if (i != wmap.end()) {
601  const auto* wp = i->second;
602 
603  const WellProductionProperties& p =
604  wp->getProductionProperties(step);
605 
606  if (! p.predictionMode) {
607  // History matching (WCONHIST/RESV)
608  SimFIBODetails::historyRates(pu, p, hrates);
609 
610  const int fipreg = 0; // Hack. Ignore FIP regions.
611  rateConverter_.calcCoeff(fipreg, pvtreg, distr);
612 
613  // WCONHIST/RESV target is sum of all
614  // observed phase rates translated to
615  // reservoir conditions. Recall sign
616  // convention: Negative for producers.
617  const double target =
618  - std::inner_product(distr.begin(), distr.end(),
619  hrates.begin(), 0.0);
620 
621  well_controls_clear(ctrl);
622  well_controls_assert_number_of_phases(ctrl, int(np));
623 
624  static const double invalid_alq = -std::numeric_limits<double>::max();
625  static const int invalid_vfp = -std::numeric_limits<int>::max();
626 
627  const int ok_resv =
628  well_controls_add_new(RESERVOIR_RATE, target,
629  invalid_alq, invalid_vfp,
630  & distr[0], ctrl);
631 
632  // For WCONHIST the BHP limit is set to 1 atm.
633  // or a value specified using WELTARG
634  double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
635  const int ok_bhp =
636  well_controls_add_new(BHP, bhp_limit,
637  invalid_alq, invalid_vfp,
638  NULL, ctrl);
639 
640  if (ok_resv != 0 && ok_bhp != 0) {
641  xw.currentControls()[*rp] = 0;
642  well_controls_set_current(ctrl, 0);
643  }
644  }
645  }
646  }
647  }
648  }
649 
650  if( wells )
651  {
652  for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
653  WellControls* ctrl = wells->ctrls[w];
654  const bool is_producer = wells->type[w] == PRODUCER;
655  if (!is_producer && wells->name[w] != 0) {
656  WellMap::const_iterator i = wmap.find(wells->name[w]);
657  if (i != wmap.end()) {
658  const auto* wp = i->second;
659  const WellInjectionProperties& injector = wp->getInjectionProperties(step);
660  if (!injector.predictionMode) {
661  //History matching WCONINJEH
662  static const double invalid_alq = -std::numeric_limits<double>::max();
663  static const int invalid_vfp = -std::numeric_limits<int>::max();
664  // For WCONINJEH the BHP limit is set to a large number
665  // or a value specified using WELTARG
666  double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
667  const int ok_bhp =
668  well_controls_add_new(BHP, bhp_limit,
669  invalid_alq, invalid_vfp,
670  NULL, ctrl);
671  if (!ok_bhp) {
672  OPM_THROW(std::runtime_error, "Failed to add well control.");
673  }
674  }
675  }
676  }
677  }
678  }
679  }
680 
681 
682 
683  void updateListEconLimited(const std::unique_ptr<Solver>& solver,
684  const Schedule& schedule,
685  const int current_step,
686  const Wells* wells,
687  const WellState& well_state,
688  DynamicListEconLimited& list_econ_limited) const
689  {
690  solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
691  well_state, list_econ_limited);
692  }
693 
694  void FIPUnitConvert(const UnitSystem& units,
695  std::vector<std::vector<double>>& fip)
696  {
697  for (size_t i = 0; i < fip.size(); ++i) {
698  FIPUnitConvert(units, fip[i]);
699  }
700  }
701 
702 
703  void FIPUnitConvert(const UnitSystem& units,
704  std::vector<double>& fip)
705  {
706  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
707  fip[0] = unit::convert::to(fip[0], unit::stb);
708  fip[1] = unit::convert::to(fip[1], unit::stb);
709  fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
710  fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
711  fip[4] = unit::convert::to(fip[4], unit::stb);
712  fip[5] = unit::convert::to(fip[5], unit::stb);
713  fip[6] = unit::convert::to(fip[6], unit::psia);
714  }
715  else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
716  fip[6] = unit::convert::to(fip[6], unit::barsa);
717  }
718  else {
719  OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
720  }
721  }
722 
723 
724  std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip)
725  {
726  std::vector<double> totals(7,0.0);
727  for (int i = 0; i < 5; ++i) {
728  for (size_t reg = 0; reg < fip.size(); ++reg) {
729  totals[i] += fip[reg][i];
730  }
731  }
732 
733  const auto& gridView = ebosSimulator_.gridManager().gridView();
734  const auto& comm = gridView.comm();
735  double pv_hydrocarbon_sum = 0.0;
736  double p_pv_hydrocarbon_sum = 0.0;
737 
738  ElementContext elemCtx(ebosSimulator_);
739  const auto& elemEndIt = gridView.template end</*codim=*/0>();
740  for (auto elemIt = gridView.template begin</*codim=*/0>();
741  elemIt != elemEndIt;
742  ++elemIt)
743  {
744  const auto& elem = *elemIt;
745  if (elem.partitionType() != Dune::InteriorEntity) {
746  continue;
747  }
748 
749  elemCtx.updatePrimaryStencil(elem);
750  elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
751 
752  const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
753  const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
754  const auto& fs = intQuants.fluidState();
755 
756  const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
757  const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
758 
759  // calculate the pore volume of the current cell. Note that the
760  // porosity returned by the intensive quantities is defined as the
761  // ratio of pore space to total cell volume and includes all pressure
762  // dependent (-> rock compressibility) and static modifiers (MULTPV,
763  // MULTREGP, NTG, PORV, MINPV and friends). Also note that because of
764  // this, the porosity returned by the intensive quantities can be
765  // outside of the physical range [0, 1] in pathetic cases.
766  const double pv =
767  ebosSimulator_.model().dofTotalVolume(cellIdx)
768  * intQuants.porosity().value();
769 
770  totals[5] += pv;
771  pv_hydrocarbon_sum += pv*hydrocarbon;
772  p_pv_hydrocarbon_sum += p*pv*hydrocarbon;
773  }
774 
775  pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum);
776  p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum);
777  totals[5] = comm.sum(totals[5]);
778  totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum);
779 
780  return totals;
781  }
782 
783 
784  void outputTimestampFIP(SimulatorTimer& timer, const std::string version)
785  {
786  std::ostringstream ss;
787  boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d %b %Y");
788  ss.imbue(std::locale(std::locale::classic(), facet));
789  ss << "\n **************************************************************************\n"
790  << " Balance at" << std::setw(10) << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day) << " Days"
791  << " *" << std::setw(30) << eclState().getTitle() << " *\n"
792  << " Report " << std::setw(4) << timer.reportStepNum() << " " << timer.currentDateTime()
793  << " * Flow version " << std::setw(11) << version << " *\n"
794  << " **************************************************************************\n";
795  OpmLog::note(ss.str());
796  }
797 
798 
799  void outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
800  {
801  std::ostringstream ss;
802  if (!reg) {
803  ss << " ===================================================\n"
804  << " : Field Totals :\n";
805  } else {
806  ss << " ===================================================\n"
807  << " : FIPNUM report region "
808  << std::setw(2) << reg << " :\n";
809  }
810  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
811  ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
812  << std::fixed << std::setprecision(0)
813  << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
814  if (!reg) {
815  ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
816  << " : Porv volumes are taken at reference conditions :\n";
817  }
818  ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
819  }
820  if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
821  ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
822  << std::fixed << std::setprecision(0)
823  << " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
824  if (!reg) {
825  ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
826  << " : Pore volumes are taken at reference conditions :\n";
827  }
828  ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
829  }
830  ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
831  << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
832  << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
833  << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
834  << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
835  << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
836  << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
837  << ":========================:==========================================:================:==========================================:\n";
838  OpmLog::note(ss.str());
839  }
840 
841 
842  const EclipseState& eclState() const
843  { return ebosSimulator_.gridManager().eclState(); }
844 
845  void extractLegacyCellPvtRegionIndex_()
846  {
847  const auto& grid = ebosSimulator_.gridManager().grid();
848  const auto& eclProblem = ebosSimulator_.problem();
849  const unsigned numCells = grid.size(/*codim=*/0);
850 
851  legacyCellPvtRegionIdx_.resize(numCells);
852  for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
853  legacyCellPvtRegionIdx_[cellIdx] =
854  eclProblem.pvtRegionIndex(cellIdx);
855  }
856  }
857 
858  void initHysteresisParams(ReservoirState& state) {
859  const int num_cells = Opm::UgGridHelpers::numCells(grid());
860 
861  typedef std::vector<double> VectorType;
862 
863  const VectorType& somax = state.getCellData( "SOMAX" );
864 
865  for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
866  ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
867  }
868 
869  if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) {
870  auto matLawManager = ebosSimulator_.problem().materialLawManager();
871 
872  VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
873  VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
874 
875  VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
876  VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
877 
878  for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
879  matLawManager->setOilWaterHysteresisParams(
880  pcSwMdc_ow[cellIdx],
881  krnSwMdc_ow[cellIdx],
882  cellIdx);
883  matLawManager->setGasOilHysteresisParams(
884  pcSwMdc_go[cellIdx],
885  krnSwMdc_go[cellIdx],
886  cellIdx);
887  }
888  }
889  }
890 
891  void extractLegacyDepth_()
892  {
893  const auto& grid = ebosSimulator_.gridManager().grid();
894  const unsigned numCells = grid.size(/*codim=*/0);
895 
896  legacyDepth_.resize(numCells);
897  for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
898  legacyDepth_[cellIdx] =
899  grid.cellCenterDepth(cellIdx);
900  }
901  }
902 
903  // Used to convert initial Reservoirstate to primary variables in the SolutionVector
904  void convertInput( const int iterationIdx,
905  const ReservoirState& reservoirState,
906  Simulator& simulator ) const
907  {
908  SolutionVector& solution = simulator.model().solution( 0 /* timeIdx */ );
909  const Opm::PhaseUsage pu = phaseUsage_;
910 
911  const std::vector<bool> active = detail::activePhases(pu);
912  bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
913  bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
914 
915  const int numCells = reservoirState.numCells();
916  const int numPhases = phaseUsage_.num_phases;
917  const auto& oilPressure = reservoirState.pressure();
918  const auto& saturations = reservoirState.saturation();
919  const auto& rs = reservoirState.gasoilratio();
920  const auto& rv = reservoirState.rv();
921  for( int cellIdx = 0; cellIdx<numCells; ++cellIdx )
922  {
923  // set non-switching primary variables
924  PrimaryVariables& cellPv = solution[ cellIdx ];
925  // set water saturation
926  if ( active[Water] ) {
927  cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
928  }
929 
930  if (has_solvent) {
931  cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
932  }
933 
934  if (has_polymer) {
935  cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
936  }
937 
938 
939  // set switching variable and interpretation
940  if ( active[Gas] ) {
941  if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
942  {
943  cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
944  cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
945  cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
946  }
947  else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
948  {
949  // this case (-> gas only with vaporized oil in the gas) is
950  // relatively expensive as it requires to compute the capillary
951  // pressure in order to get the gas phase pressure. (the reason why
952  // ebos uses the gas pressure here is that it makes the common case
953  // of the primary variable switching code fast because to determine
954  // whether the oil phase appears one needs to compute the Rv value
955  // for the saturated gas phase and if this is not available as a
956  // primary variable, it needs to be computed.) luckily for here, the
957  // gas-only case is not too common, so the performance impact of this
958  // is limited.
959  typedef Opm::SimpleModularFluidState<double,
960  /*numPhases=*/3,
961  /*numComponents=*/3,
962  FluidSystem,
963  /*storePressure=*/false,
964  /*storeTemperature=*/false,
965  /*storeComposition=*/false,
966  /*storeFugacity=*/false,
967  /*storeSaturation=*/true,
968  /*storeDensity=*/false,
969  /*storeViscosity=*/false,
970  /*storeEnthalpy=*/false> SatOnlyFluidState;
971  SatOnlyFluidState fluidState;
972  if ( active[Water] ) {
973  fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
974  }
975  else {
976  fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
977  }
978  fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
979  fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
980 
981  double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 };
982  const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
983  MaterialLaw::capillaryPressures(pC, matParams, fluidState);
984  double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
985 
986  cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
987  cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
988  cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
989  }
990  else
991  {
992  assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
993  cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
994  cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
995  cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
996  }
997  } else {
998  // for oil-water case oil pressure should be used as primary variable
999  cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
1000  }
1001  }
1002 
1003  // store the solution at the beginning of the time step
1004  if( iterationIdx == 0 )
1005  {
1006  simulator.model().solution( 1 /* timeIdx */ ) = solution;
1007  }
1008  }
1009 
1010  RateConverterType createRateConverter_() {
1011  RateConverterType rate_converter(phaseUsage_,
1012  std::vector<int>(AutoDiffGrid::numCells(grid()), 0)); // FIP = 0
1013  return rate_converter;
1014  }
1015 
1016 
1017  // Data.
1018  Simulator& ebosSimulator_;
1019 
1020  std::vector<int> legacyCellPvtRegionIdx_;
1021  std::vector<double> legacyDepth_;
1022  typedef typename Solver::SolverParameters SolverParameters;
1023 
1024  SimulatorReport failureReport_;
1025 
1026  const ParameterGroup param_;
1027  ModelParameters model_param_;
1028  SolverParameters solver_param_;
1029 
1030  // Observed objects.
1031  NewtonIterationBlackoilInterface& solver_;
1032  PhaseUsage phaseUsage_;
1033  // Misc. data
1034  const bool has_disgas_;
1035  const bool has_vapoil_;
1036  bool terminal_output_;
1037  // output_writer
1038  OutputWriter& output_writer_;
1039  RateConverterType rateConverter_;
1040  // The names of wells that should be defunct
1041  // (e.g. in a parallel run when they are handeled by
1042  // a different process)
1043  std::unordered_set<std::string> defunct_well_names_;
1044 
1045  // Whether this a parallel simulation or not
1046  bool is_parallel_run_;
1047 
1048 };
1049 
1050 } // namespace Opm
1051 
1052 #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
A nonlinear solver class suitable for general fully-implicit models, as well as pressure, transport and sequential models.
Definition: NonlinearSolver.hpp:37
Extra data to read/write for OPM restarting.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:199
int currentStepNum() const
Current step number.
Definition: SimulatorTimer.cpp:77
bool isIORank() const
Whether this process does write to disk.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:284
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions...
Definition: RateConverter.hpp:407
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:66
SimulatorFullyImplicitBlackoilEbos(Simulator &ebosSimulator, const ParameterGroup &param, NewtonIterationBlackoilInterface &linsolver, const bool has_disgas, const bool has_vapoil, const EclipseState &, OutputWriter &output_writer, const std::unordered_set< std::string > &defunct_well_names)
Initialise from parameters and objects to observe.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:96
std::string moduleVersionName()
Return the version name of the module, for example &quot;2015.10&quot; (for a release branch) or &quot;2016...
Definition: moduleVersion.cpp:28
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
A model implementation for three-phase black oil.
Definition: BlackoilModelEbos.hpp:104
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
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
const SimulatorReport & failureReport() const
Returns the simulator report for the failed substeps of the simulation.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:472
a simulator for the blackoil model
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:48
double currentStepLength() const
Current step length.
Definition: SimulatorTimer.cpp:91
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
void initWellStateMSWell(const Wells *wells, const std::vector< const Well * > &wells_ecl, const int time_step, const PhaseUsage &pu, const PrevWellState &prev_well_state)
init the MS well related.
Definition: WellStateFullyImplicitBlackoil.hpp:298
void writeTimeStep(const SimulatorTimerInterface &timer, const SimulationDataContainer &reservoirState, const Opm::WellStateFullyImplicitBlackoil &wellState, const Model &physicalModel, const bool substep=false, const double nextstep=-1.0, const SimulatorReport &simulatorReport=SimulatorReport())
Write a blackoil reservoir state to disk for later inspection with visualization tools like ResInsigh...
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:995
SimulatorReport run(SimulatorTimer &timer, ReservoirState &state)
Run the simulation.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:136
bool output() const
return true if output is enabled
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:281
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Wrapper class for VTK, Matlab, and ECL output.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:206
const std::string & outputDirectory() const
return output directory
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:278
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:589
Definition: SimulatorTimer.hpp:34
bool use_multisegment_well_
Whether to use MultisegmentWell to handle multisegment wells it is something temporary before the mul...
Definition: BlackoilModelParameters.hpp:89