FlowMain.hpp
1 /*
2  Copyright 2013, 2014, 2015 SINTEF ICT, Applied Mathematics.
3  Copyright 2014 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2015 IRIS AS
5  Copyright 2014 STATOIL ASA.
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 3 of the License, or
12  (at your option) any later version.
13 
14  OPM is distributed in the hope that it will be useful,
15  but WITHOUT ANY WARRANTY; without even the implied warranty of
16  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  GNU General Public License for more details.
18 
19  You should have received a copy of the GNU General Public License
20  along with OPM. If not, see <http://www.gnu.org/licenses/>.
21 */
22 
23 #ifndef OPM_FLOWMAIN_HEADER_INCLUDED
24 #define OPM_FLOWMAIN_HEADER_INCLUDED
25 
26 #include <opm/common/utility/platform_dependent/disable_warnings.h>
27 
28 
29 #include <dune/common/version.hh>
30 #include <dune/common/parallel/mpihelper.hh>
31 
32 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
33 
34 
35 #include <opm/core/grid/GridManager.hpp>
36 #include <opm/autodiff/GridHelpers.hpp>
37 #include <opm/autodiff/createGlobalCellArray.hpp>
38 #include <opm/autodiff/GridInit.hpp>
39 #include <opm/simulators/ParallelFileMerger.hpp>
40 #include <opm/simulators/ensureDirectoryExists.hpp>
41 
42 #include <opm/core/wells.h>
43 #include <opm/core/wells/WellsManager.hpp>
44 #include <opm/common/ErrorMacros.hpp>
45 #include <opm/core/simulator/initState.hpp>
46 #include <opm/core/simulator/initStateEquil.hpp>
47 #include <opm/core/simulator/SimulatorReport.hpp>
48 #include <opm/simulators/timestepping/SimulatorTimer.hpp>
49 #include <opm/core/utility/miscUtilities.hpp>
50 #include <opm/core/utility/parameters/ParameterGroup.hpp>
51 #include <opm/simulators/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available). \TODO: Fix.
52 
53 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
54 #include <opm/core/props/BlackoilPropertiesBasic.hpp>
55 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
56 #include <opm/core/props/rock/RockCompressibility.hpp>
57 #include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
58 #include <opm/core/linalg/LinearSolverFactory.hpp>
59 #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
60 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
61 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
62 
63 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
64 
65 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
66 #include <opm/autodiff/RedistributeDataHandles.hpp>
67 #include <opm/autodiff/moduleVersion.hpp>
68 #include <opm/autodiff/MissingFeatures.hpp>
69 
70 #include <opm/core/utility/share_obj.hpp>
71 #include <opm/core/utility/initHydroCarbonState.hpp>
72 #include <opm/common/OpmLog/OpmLog.hpp>
73 #include <opm/common/OpmLog/EclipsePRTLog.hpp>
74 #include <opm/common/OpmLog/LogUtil.hpp>
75 #include <opm/parser/eclipse/Deck/Deck.hpp>
76 #include <opm/parser/eclipse/Parser/Parser.hpp>
77 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
78 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
79 #include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
80 #include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
81 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
82 #include <opm/common/ResetLocale.hpp>
83 
84 #include <boost/filesystem.hpp>
85 #include <boost/algorithm/string.hpp>
86 
87 
88 #ifdef _OPENMP
89 #include <omp.h>
90 #endif
91 
92 #include <memory>
93 #include <algorithm>
94 #include <cstdlib>
95 #include <iostream>
96 #include <vector>
97 #include <numeric>
98 #include <cstdlib>
99 #include <stdexcept>
100 
101 
102 
103 
104 namespace Opm
105 {
106 
107 
108  namespace detail
109  {
110  boost::filesystem::path simulationCaseName( const std::string& casename );
111  int64_t convertMessageType(const Message::type& mtype);
112  }
113 
114 
115 
116 
117 
118 
121  template <class Implementation, class Grid, class Simulator>
123  {
124  public:
125 
126 
132  int execute(int argc, char** argv)
133  try {
134  // we always want to use the default locale, and thus spare us the trouble
135  // with incorrect locale settings.
136  resetLocale();
137 
138  // Setup.
139  asImpl().setupParallelism(argc, argv);
140  asImpl().printStartupMessage();
141  const bool ok = asImpl().setupParameters(argc, argv);
142  if (!ok) {
143  return EXIT_FAILURE;
144  }
145  asImpl().readDeckInput();
146  asImpl().setupOutput();
147  asImpl().setupLogging();
148  asImpl().extractMessages();
149  asImpl().setupGridAndProps();
150  asImpl().runDiagnostics();
151  asImpl().setupState();
152  asImpl().writeInit();
153  asImpl().distributeData();
154  asImpl().setupOutputWriter();
155  asImpl().setupLinearSolver();
156  asImpl().createSimulator();
157 
158  // Run.
159  auto ret = asImpl().runSimulator();
160 
161  asImpl().mergeParallelLogFiles();
162 
163  return ret;
164  }
165  catch (const std::exception &e) {
166  std::ostringstream message;
167  message << "Program threw an exception: " << e.what();
168 
169  if( output_cout_ )
170  {
171  // in some cases exceptions are thrown before the logging system is set
172  // up.
173  if (OpmLog::hasBackend("STREAMLOG")) {
174  OpmLog::error(message.str());
175  }
176  else {
177  std::cout << message.str() << "\n";
178  }
179  }
180 
181  return EXIT_FAILURE;
182  }
183 
184 
185 
186  protected:
187 
188 
189 
190 
191 
192  // ------------ Types ------------
193 
194 
195  typedef BlackoilPropsAdFromDeck FluidProps;
196  typedef FluidProps::MaterialLawManager MaterialLawManager;
197  typedef typename Simulator::ReservoirState ReservoirState;
198  typedef typename Simulator::OutputWriter OutputWriter;
199 
200 
201  // ------------ Data members ------------
202 
203 
204  // The comments indicate in which method the
205  // members first occur.
206 
207  // setupParallelism()
208  int mpi_rank_ = 0;
209  bool output_cout_ = false;
210  bool must_distribute_ = false;
211  // setupParameters()
212  ParameterGroup param_;
213  // setupOutput()
214  bool output_to_files_ = false;
215  std::string output_dir_ = std::string(".");
216  // readDeckInput()
217  std::shared_ptr<Deck> deck_;
218  std::shared_ptr<EclipseState> eclipse_state_;
219  // setupGridAndProps()
220  std::unique_ptr<GridInit<Grid>> grid_init_;
221  std::shared_ptr<MaterialLawManager> material_law_manager_;
222  std::unique_ptr<FluidProps> fluidprops_;
223  std::unique_ptr<RockCompressibility> rock_comp_;
224  std::array<double, 3> gravity_;
225  bool use_local_perm_ = true;
226  std::unique_ptr<DerivedGeology> geoprops_;
227  // setupState()
228  std::unique_ptr<ReservoirState> state_;
229 
230  std::vector<double> threshold_pressures_;
231  // distributeData()
232  boost::any parallel_information_;
233  // setupOutputWriter()
234  std::unique_ptr<EclipseIO> eclipse_writer_;
235  std::unique_ptr<OutputWriter> output_writer_;
236  // setupLinearSolver
237  std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
238  // createSimulator()
239  std::unique_ptr<Simulator> simulator_;
240  // create log file
241  std::string logFile_;
242  // The names of wells that are artifically defunct in parallel runs.
243  // Those wells are handled on a another process.
244  std::unordered_set<std::string> defunct_well_names_;
245  // ------------ Methods ------------
246 
247 
248  // Set up MPI and OpenMP.
249  // Writes to:
250  // output_cout_
251  // must_distribute_
252  void setupParallelism(int argc, char** argv)
253  {
254  // MPI setup.
255  // Must ensure an instance of the helper is created to initialise MPI.
256  // For a build without MPI the Dune::FakeMPIHelper is used, so rank will
257  // be 0 and size 1.
258  const Dune::MPIHelper& mpi_helper = Dune::MPIHelper::instance(argc, argv);
259  mpi_rank_ = mpi_helper.rank();
260  const int mpi_size = mpi_helper.size();
261  output_cout_ = ( mpi_rank_ == 0 );
262  must_distribute_ = ( mpi_size > 1 );
263 
264 #ifdef _OPENMP
265  // OpenMP setup.
266  if (!getenv("OMP_NUM_THREADS")) {
267  // Default to at most 4 threads, regardless of
268  // number of cores (unless ENV(OMP_NUM_THREADS) is defined)
269  int num_cores = omp_get_num_procs();
270  int num_threads = std::min(4, num_cores);
271  omp_set_num_threads(num_threads);
272  }
273 #pragma omp parallel
274  if (omp_get_thread_num() == 0) {
275  // omp_get_num_threads() only works as expected within a parallel region.
276  const int num_omp_threads = omp_get_num_threads();
277  if (mpi_size == 1) {
278  std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl;
279  } else {
280  std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl;
281  }
282  }
283 #endif
284  }
285 
287  bool cartesianAdjacent(const Grid& grid, int g1, int g2) {
288  // we need cartDims from UgGridHelpers
289  using namespace UgGridHelpers;
290 
291  int diff = std::abs(g1 - g2);
292 
293  const int * dimens = cartDims(grid);
294  if (diff == 1)
295  return true;
296  if (diff == dimens[0])
297  return true;
298  if (diff == dimens[0] * dimens[1])
299  return true;
300 
301  return false;
302  }
303 
304  // Print startup message if on output rank.
305  void printStartupMessage()
306  {
307  if (output_cout_) {
308  const std::string version = moduleVersionName();
309  std::cout << "**********************************************************************\n";
310  std::cout << "* *\n";
311  std::cout << "* This is flow_legacy (version " << version << ")"
312  << std::string(26 - version.size(), ' ') << "*\n";
313  std::cout << "* *\n";
314  std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
315  std::cout << "* and is part of OPM. For more information see: *\n";
316  std::cout << "* http://opm-project.org *\n";
317  std::cout << "* *\n";
318  std::cout << "**********************************************************************\n\n";
319  }
320  }
321 
322  // Read parameters, see if a deck was specified on the command line, and if
323  // it was, insert it into parameters.
324  // Writes to:
325  // param_
326  // Returns true if ok, false if not.
327  bool setupParameters(int argc, char** argv)
328  {
329  param_ = ParameterGroup(argc, argv, false, output_cout_);
330 
331  // See if a deck was specified on the command line.
332  if (!param_.unhandledArguments().empty()) {
333  if (param_.unhandledArguments().size() != 1) {
334  std::cerr << "You can only specify a single input deck on the command line.\n";
335  return false;
336  } else {
337  const auto casename = detail::simulationCaseName( param_.unhandledArguments()[ 0 ] );
338  param_.insertParameter("deck_filename", casename.string() );
339  }
340  }
341 
342  // We must have an input deck. Grid and props will be read from that.
343  if (!param_.has("deck_filename")) {
344  std::cerr << "This program must be run with an input deck.\n"
345  "Specify the deck filename either\n"
346  " a) as a command line argument by itself\n"
347  " b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
348  " c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
349  return false;
350  }
351  return true;
352  }
353 
354 
355 
356 
357 
358  // Set output_to_files_ and set/create output dir. Write parameter file.
359  // Writes to:
360  // output_to_files_
361  // output_dir_
362  // Throws std::runtime_error if failed to create (if requested) output dir.
363  void setupOutput()
364  {
365  output_to_files_ = output_cout_ && param_.getDefault("output", true);
366 
367  // Setup output directory.
368  auto& ioConfig = eclipse_state_->getIOConfig();
369  // Default output directory is the directory where the deck is found.
370  const std::string default_output_dir = ioConfig.getOutputDir();
371  output_dir_ = param_.getDefault("output_dir", default_output_dir);
372  // Override output directory if user specified.
373  ioConfig.setOutputDir(output_dir_);
374 
375  // Write parameters used for later reference. (only if rank is zero)
376  if (output_to_files_) {
377  // Create output directory if needed.
378  ensureDirectoryExists(output_dir_);
379  // Write simulation parameters.
380  param_.writeParam(output_dir_ + "/simulation.param");
381  }
382  }
383 
384 
385 
386 
387 
388  // Setup OpmLog backend with output_dir.
389  void setupLogging()
390  {
391  std::string deck_filename = param_.get<std::string>("deck_filename");
392  // create logFile
393  using boost::filesystem::path;
394  path fpath(deck_filename);
395  std::string baseName;
396  std::ostringstream debugFileStream;
397  std::ostringstream logFileStream;
398 
399  if (boost::to_upper_copy(path(fpath.extension()).string()) == ".DATA") {
400  baseName = path(fpath.stem()).string();
401  } else {
402  baseName = path(fpath.filename()).string();
403  }
404 
405  logFileStream << output_dir_ << "/" << baseName;
406  debugFileStream << output_dir_ << "/" << "." << baseName;
407 
408  if ( must_distribute_ && mpi_rank_ != 0 )
409  {
410  // Added rank to log file for non-zero ranks.
411  // This prevents message loss.
412  debugFileStream << "."<< mpi_rank_;
413  // If the following file appears then there is a bug.
414  logFileStream << "." << mpi_rank_;
415  }
416  logFileStream << ".PRT";
417  debugFileStream << ".DEBUG";
418 
419  std::string debugFile = debugFileStream.str();
420  logFile_ = logFileStream.str();
421 
422  std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::NoDebugMessageTypes, false, output_cout_);
423  const bool all_to_terminal = param_.getDefault("all_messages_to_terminal", false);
424  const auto terminal_msg_types = all_to_terminal ? Log::DefaultMessageTypes : Log::StdoutMessageTypes;
425  std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(std::cout, terminal_msg_types);
426  OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
427  OpmLog::addBackend( "STREAMLOG", streamLog);
428  std::shared_ptr<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
429  OpmLog::addBackend( "DEBUGLOG" , debugLog);
430  const auto& msgLimits = eclipse_state_->getSchedule().getMessageLimits();
431  const std::map<int64_t, int> limits = {{Log::MessageType::Note, msgLimits.getCommentPrintLimit(0)},
432  {Log::MessageType::Info, msgLimits.getMessagePrintLimit(0)},
433  {Log::MessageType::Warning, msgLimits.getWarningPrintLimit(0)},
434  {Log::MessageType::Error, msgLimits.getErrorPrintLimit(0)},
435  {Log::MessageType::Problem, msgLimits.getProblemPrintLimit(0)},
436  {Log::MessageType::Bug, msgLimits.getBugPrintLimit(0)}};
437  prtLog->setMessageLimiter(std::make_shared<MessageLimiter>());
438  prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
439  streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
440  streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true));
441 
442  // Read parameters.
443  if ( output_cout_ )
444  {
445  OpmLog::debug("\n--------------- Reading parameters ---------------\n");
446  }
447  }
448 
449 
450 
451 
452 
453  void mergeParallelLogFiles()
454  {
455  // force closing of all log files.
456  OpmLog::removeAllBackends();
457 
458  if( mpi_rank_ != 0 || !must_distribute_ || !output_to_files_ )
459  {
460  return;
461  }
462 
463  namespace fs = boost::filesystem;
464  fs::path output_path(".");
465  if ( param_.has("output_dir") )
466  {
467  output_path = fs::path(output_dir_);
468  }
469 
470  fs::path deck_filename(param_.get<std::string>("deck_filename"));
471 
472  std::for_each(fs::directory_iterator(output_path),
473  fs::directory_iterator(),
474  detail::ParallelFileMerger(output_path, deck_filename.stem().string()));
475  }
476 
477  // Parser the input and creates the Deck and EclipseState objects.
478  // Writes to:
479  // deck_
480  // eclipse_state_
481  // May throw if errors are encountered, here configured to be somewhat tolerant.
482  void readDeckInput()
483  {
484  std::string deck_filename = param_.get<std::string>("deck_filename");
485 
486  // Create Parser
487  Parser parser;
488 
489  // Create Deck and EclipseState.
490  try {
491  ParseContext parseContext({ { ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE },
492  { ParseContext::PARSE_MISSING_DIMS_KEYWORD, InputError::WARN },
493  { ParseContext::SUMMARY_UNKNOWN_WELL, InputError::WARN },
494  { ParseContext::SUMMARY_UNKNOWN_GROUP, InputError::WARN }});
495  deck_ = std::make_shared< Deck >( parser.parseFile(deck_filename, parseContext) );
496  checkDeck(*deck_, parser);
497 
498  if ( output_cout_)
499  {
500  MissingFeatures::checkKeywords(*deck_);
501  }
502 
503  eclipse_state_.reset(new EclipseState(*deck_, parseContext));
504  }
505  catch (const std::invalid_argument& e) {
506  std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile_ << std::endl;
507  std::cerr << "Exception caught: " << e.what() << std::endl;
508  throw;
509  }
510 
511  // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
512  if (param_.has("output_interval")) {
513  const int output_interval = param_.get<int>("output_interval");
514  eclipse_state_->getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) );
515  }
516 
517  // Possible to force initialization only behavior (NOSIM).
518  if (param_.has("nosim")) {
519  const bool nosim = param_.get<bool>("nosim");
520  auto& ioConfig = eclipse_state_->getIOConfig();
521  ioConfig.overrideNOSIM( nosim );
522  }
523  }
524 
525 
526 
527 
528 
529  // Create grid and property objects.
530  // Writes to:
531  // grid_init_
532  // material_law_manager_
533  // fluidprops_
534  // rock_comp_
535  // gravity_
536  // use_local_perm_
537  // geoprops_
538  void setupGridAndProps()
539  {
540  // Create grid.
541  const std::vector<double>& porv =
542  eclipse_state_->get3DProperties().getDoubleGridProperty("PORV").getData();
543  grid_init_.reset(new GridInit<Grid>(*eclipse_state_, porv));
544  const Grid& grid = grid_init_->grid();
545 
546  // Create material law manager.
547  std::vector<int> compressedToCartesianIdx;
548  Opm::createGlobalCellArray(grid, compressedToCartesianIdx);
549  material_law_manager_.reset(new MaterialLawManager());
550  material_law_manager_->initFromDeck(*deck_, *eclipse_state_, compressedToCartesianIdx);
551 
552  // Rock and fluid properties.
553  fluidprops_.reset(new BlackoilPropsAdFromDeck(*deck_, *eclipse_state_, material_law_manager_, grid));
554 
555  // Rock compressibility.
556  rock_comp_.reset(new RockCompressibility(*eclipse_state_, output_cout_));
557 
558  // Gravity.
559  assert(UgGridHelpers::dimensions(grid) == 3);
560  gravity_.fill(0.0);
561  gravity_[2] = deck_->hasKeyword("NOGRAV")
562  ? param_.getDefault("gravity", 0.0)
563  : param_.getDefault("gravity", unit::gravity);
564 
565  // Geological properties
566  use_local_perm_ = param_.getDefault("use_local_perm", use_local_perm_);
567  geoprops_.reset(new DerivedGeology(grid, *fluidprops_, *eclipse_state_, use_local_perm_, gravity_.data()));
568  }
569 
570 
571 
572 
573 
574  // Initialise the reservoir state. Updated fluid props for SWATINIT.
575  // Writes to:
576  // state_
577  // threshold_pressures_
578  // fluidprops_ (if SWATINIT is used)
579  void setupState()
580  {
581  const PhaseUsage pu = Opm::phaseUsageFromDeck(*deck_);
582  const Grid& grid = grid_init_->grid();
583 
584  // Need old-style fluid object for init purposes (only).
585  BlackoilPropertiesFromDeck props( *deck_, *eclipse_state_, material_law_manager_,
586  Opm::UgGridHelpers::numCells(grid),
587  Opm::UgGridHelpers::globalCell(grid),
588  Opm::UgGridHelpers::cartDims(grid),
589  param_);
590 
591 
592  // Init state variables (saturation and pressure).
593  if (param_.has("init_saturation")) {
594  state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
595  Opm::UgGridHelpers::numFaces(grid),
596  props.numPhases() ));
597 
598  initStateBasic(Opm::UgGridHelpers::numCells(grid),
599  Opm::UgGridHelpers::globalCell(grid),
600  Opm::UgGridHelpers::cartDims(grid),
601  Opm::UgGridHelpers::numFaces(grid),
602  Opm::UgGridHelpers::faceCells(grid),
603  Opm::UgGridHelpers::beginFaceCentroids(grid),
604  Opm::UgGridHelpers::beginCellCentroids(grid),
605  Opm::UgGridHelpers::dimensions(grid),
606  props, param_, gravity_[2], *state_);
607 
608  initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_);
609 
610  enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
611  if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
612  const int numPhases = props.numPhases();
613  const int numCells = Opm::UgGridHelpers::numCells(grid);
614 
615  // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState.
616  auto& gor = state_->getCellData( BlackoilState::GASOILRATIO );
617  const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL );
618  for (int c = 0; c < numCells; ++c) {
619  // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field.
620  gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]];
621  }
622  }
623  } else if (deck_->hasKeyword("EQUIL")) {
624  // Which state class are we really using - what a f... mess?
625  state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
626  Opm::UgGridHelpers::numFaces(grid),
627  props.numPhases()));
628 
629  initStateEquil(grid, props, *deck_, *eclipse_state_, gravity_[2], *state_);
630  //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
631  } else {
632  state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
633  Opm::UgGridHelpers::numFaces(grid),
634  props.numPhases()));
635  initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
636  Opm::UgGridHelpers::globalCell(grid),
637  Opm::UgGridHelpers::numFaces(grid),
638  Opm::UgGridHelpers::faceCells(grid),
639  Opm::UgGridHelpers::beginFaceCentroids(grid),
640  Opm::UgGridHelpers::beginCellCentroids(grid),
641  Opm::UgGridHelpers::dimensions(grid),
642  props, *deck_, gravity_[2], *state_);
643  }
644 
645  // Threshold pressures.
646  std::map<std::pair<int, int>, double> maxDp;
647  computeMaxDp(maxDp, *deck_, *eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]);
648  threshold_pressures_ = thresholdPressures(*deck_, *eclipse_state_, grid, maxDp);
649  std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(*eclipse_state_, geoprops_->nnc(), maxDp);
650  threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
651 
652  // The capillary pressure is scaled in fluidprops_ to match the scaled capillary pressure in props.
653  if (deck_->hasKeyword("SWATINIT")) {
654  const int numCells = Opm::UgGridHelpers::numCells(grid);
655  std::vector<int> cells(numCells);
656  for (int c = 0; c < numCells; ++c) { cells[c] = c; }
657  std::vector<double> pc = state_->saturation();
658  props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr);
659  fluidprops_->setSwatInitScaling(state_->saturation(), pc);
660  }
661  initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck_->hasKeyword("DISGAS"), deck_->hasKeyword("VAPOIL"));
662 
663 
664  }
665 
666 
667 
668 
669 
670  // Distribute the grid, properties and state.
671  // Writes to:
672  // grid_init_->grid()
673  // state_
674  // fluidprops_
675  // geoprops_
676  // material_law_manager_
677  // parallel_information_
678  void distributeData()
679  {
680  // At this point all properties and state variables are correctly initialized
681  // If there are more than one processors involved, we now repartition the grid
682  // and initilialize new properties and states for it.
683  if (must_distribute_) {
684  defunct_well_names_ =
685  distributeGridAndData(grid_init_->grid(), *deck_, *eclipse_state_,
686  *state_, *fluidprops_, *geoprops_,
687  material_law_manager_, threshold_pressures_,
688  parallel_information_, use_local_perm_);
689  }
690  }
691 
692 
693 
694 
695 
696  // Extract messages from parser.
697  // Writes to:
698  // OpmLog singleton.
699  void extractMessages()
700  {
701  if ( !output_cout_ )
702  {
703  return;
704  }
705 
706  auto extractMessage = [](const Message& msg) {
707  auto log_type = detail::convertMessageType(msg.mtype);
708  const auto& location = msg.location;
709  if (location) {
710  OpmLog::addTaggedMessage(log_type, "Parser message", Log::fileMessage(location.filename, location.lineno, msg.message));
711  } else {
712  OpmLog::addTaggedMessage(log_type, "Parser message", msg.message);
713  }
714  };
715 
716  // Extract messages from Deck.
717  for(const auto& msg : deck_->getMessageContainer()) {
718  extractMessage(msg);
719  }
720 
721  // Extract messages from EclipseState.
722  for (const auto& msg : eclipse_state_->getMessageContainer()) {
723  extractMessage(msg);
724  }
725  }
726 
727 
728 
729 
730 
731  // Run diagnostics.
732  // Writes to:
733  // OpmLog singleton.
734  void runDiagnostics()
735  {
736  if( ! output_cout_ )
737  {
738  return;
739  }
740 
741  // Run relperm diagnostics
742  RelpermDiagnostics diagnostic;
743  diagnostic.diagnosis(*eclipse_state_, *deck_, grid_init_->grid());
744  }
745 
746 
747  void writeInit()
748  {
749  bool output = param_.getDefault("output", true);
750  bool output_ecl = param_.getDefault("output_ecl", true);
751  const Grid& grid = grid_init_->grid();
752  if( output && output_ecl && output_cout_)
753  {
754  const EclipseGrid& inputGrid = eclipse_state_->getInputGrid();
755  eclipse_writer_.reset(new EclipseIO(*eclipse_state_, UgGridHelpers::createEclipseGrid( grid , inputGrid )));
756  eclipse_writer_->writeInitial(geoprops_->simProps(grid),
757  geoprops_->nonCartesianConnections());
758  }
759  }
760 
761 
762  // Setup output writer.
763  // Writes to:
764  // output_writer_
765  void setupOutputWriter()
766  {
767  // create output writer after grid is distributed, otherwise the parallel output
768  // won't work correctly since we need to create a mapping from the distributed to
769  // the global view
770  output_writer_.reset(new OutputWriter(grid_init_->grid(),
771  param_,
772  *eclipse_state_,
773  std::move(eclipse_writer_),
774  Opm::phaseUsageFromDeck(*deck_)));
775  }
776 
777 
778 
779 
780 
781  // Setup linear solver.
782  // Writes to:
783  // fis_solver_
784  void setupLinearSolver()
785  {
786  const std::string cprSolver = "cpr";
787  const std::string interleavedSolver = "interleaved";
788  const std::string directSolver = "direct";
789  std::string flowDefaultSolver = interleavedSolver;
790 
791  if (!param_.has("solver_approach")) {
792  if (eclipse_state_->getSimulationConfig().useCPR()) {
793  flowDefaultSolver = cprSolver;
794  }
795  }
796 
797  const std::string solver_approach = param_.getDefault("solver_approach", flowDefaultSolver);
798 
799  if (solver_approach == cprSolver) {
800  fis_solver_.reset(new NewtonIterationBlackoilCPR(param_, parallel_information_));
801  } else if (solver_approach == interleavedSolver) {
802  fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_));
803  } else if (solver_approach == directSolver) {
804  fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_));
805  } else {
806  OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
807  }
808  }
809 
810 
811 
812 
813 
814  // Run the simulator.
815  // Returns EXIT_SUCCESS if it does not throw.
816  int runSimulator()
817  {
818  const auto& schedule = eclipse_state_->getSchedule();
819  const auto& timeMap = schedule.getTimeMap();
820  auto& ioConfig = eclipse_state_->getIOConfig();
821  SimulatorTimer simtimer;
822 
823  // initialize variables
824  const auto& initConfig = eclipse_state_->getInitConfig();
825  simtimer.init(timeMap, (size_t)initConfig.getRestartStep());
826 
827  if (!ioConfig.initOnly()) {
828  if (output_cout_) {
829  std::string msg;
830  msg = "\n\n================ Starting main simulation loop ===============\n";
831  OpmLog::info(msg);
832  }
833 
834  SimulatorReport fullReport = simulator_->run(simtimer, *state_);
835 
836  if (output_cout_) {
837  std::ostringstream ss;
838  ss << "\n\n================ End of simulation ===============\n\n";
839  fullReport.reportFullyImplicit(ss);
840  OpmLog::info(ss.str());
841  if (param_.anyUnused()) {
842  // This allows a user to catch typos and misunderstandings in the
843  // use of simulator parameters.
844  std::cout << "-------------------- Unused parameters: --------------------\n";
845  param_.displayUsage();
846  std::cout << "----------------------------------------------------------------" << std::endl;
847  }
848  }
849 
850  if (output_to_files_) {
851  std::string filename = output_dir_ + "/walltime.txt";
852  std::fstream tot_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
853  fullReport.reportParam(tot_os);
854  }
855  } else {
856  if (output_cout_) {
857  std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
858  }
859 
860  }
861  return EXIT_SUCCESS;
862  }
863 
864 
865 
866 
867 
868  // Access the most-derived class used for
869  // static polymorphism (CRTP).
870  Implementation& asImpl()
871  {
872  return static_cast<Implementation&>(*this);
873  }
874 
875 
876  }; // class FlowMainBase
877 
878 
879 
880 
881 
882 
883  // The FlowMain class is the basic black-oil simulator case.
884  template <class Grid, class Simulator>
885  class FlowMain : public FlowMainBase<FlowMain<Grid, Simulator>, Grid, Simulator>
886  {
887  protected:
888  using Base = FlowMainBase<FlowMain<Grid, Simulator>, Grid, Simulator>;
889  friend Base;
890 
891  // Create simulator instance.
892  // Writes to:
893  // simulator_
894  void createSimulator()
895  {
896  // Create the simulator instance.
897  Base::simulator_.reset(new Simulator(Base::param_,
898  Base::grid_init_->grid(),
899  *Base::geoprops_,
900  *Base::fluidprops_,
901  Base::rock_comp_->isActive() ? Base::rock_comp_.get() : nullptr,
902  *Base::fis_solver_,
903  Base::gravity_.data(),
904  Base::deck_->hasKeyword("DISGAS"),
905  Base::deck_->hasKeyword("VAPOIL"),
906  Base::eclipse_state_,
907  *Base::output_writer_,
908  Base::threshold_pressures_,
909  Base::defunct_well_names_));
910  }
911  };
912 
913 
914 
915 
916 
917 
918  namespace detail
919  {
920 
921  boost::filesystem::path simulationCaseName( const std::string& casename ) {
922  namespace fs = boost::filesystem;
923 
924  const auto exists = []( const fs::path& f ) -> bool {
925  if( !fs::exists( f ) ) return false;
926 
927  if( fs::is_regular_file( f ) ) return true;
928 
929  return fs::is_symlink( f )
930  && fs::is_regular_file( fs::read_symlink( f ) );
931  };
932 
933  auto simcase = fs::path( casename );
934 
935  if( exists( simcase ) ) {
936  return simcase;
937  }
938 
939  for( const auto& ext : { std::string("data"), std::string("DATA") } ) {
940  if( exists( simcase.replace_extension( ext ) ) ) {
941  return simcase;
942  }
943  }
944 
945  throw std::invalid_argument( "Cannot find input case " + casename );
946  }
947 
948 
949 
950 
951 
952  int64_t convertMessageType(const Message::type& mtype)
953  {
954  switch (mtype) {
955  case Message::type::Debug:
956  return Log::MessageType::Debug;
957  case Message::type::Info:
958  return Log::MessageType::Info;
959  case Message::type::Warning:
960  return Log::MessageType::Warning;
961  case Message::type::Error:
962  return Log::MessageType::Error;
963  case Message::type::Problem:
964  return Log::MessageType::Problem;
965  case Message::type::Bug:
966  return Log::MessageType::Bug;
967  case Message::type::Note:
968  return Log::MessageType::Note;
969  }
970  throw std::logic_error("Invalid messages type!\n");
971  }
972 
973 
974  } // namespace detail
975 
976 
977 } // namespace Opm
978 
979 #endif // OPM_FLOWMAIN_HEADER_INCLUDED
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016...
Definition: moduleVersion.cpp:28
This class encapsulates the setup and running of a simulator based on an input deck.
Definition: FlowMain.hpp:122
Definition: FlowMain.hpp:885
bool cartesianAdjacent(const Grid &grid, int g1, int g2)
checks cartesian adjacency of global indices g1 and g2
Definition: FlowMain.hpp:287
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
void ensureDirectoryExists(const boost::filesystem::path &dirpath)
The directory pointed to by &#39;dirpath&#39; will be created if it does not already exist.
Definition: ensureDirectoryExists.cpp:30
int execute(int argc, char **argv)
This is the main function of Flow.
Definition: FlowMain.hpp:132
void computeMaxDp(std::map< std::pair< int, int >, double > &maxDp, const Deck &deck, const EclipseState &eclipseState, const Grid &grid, const BlackoilState &initialState, const BlackoilPropertiesFromDeck &props, const double gravity)
Compute the maximum gravity corrected pressure difference of all equilibration regions given a reserv...
Definition: thresholdPressures.hpp:47
void createGlobalCellArray(const Grid &grid, std::vector< int > &dest)
Create a mapping from a global cell index of a grid to the logically Cartesian index of the ECL deck...
Definition: createGlobalCellArray.hpp:31
std::vector< double > thresholdPressuresNNC(const EclipseState &eclipseState, const NNC &nnc, const std::map< std::pair< int, int >, double > &maxDp)
Get a vector of pressure thresholds from either EclipseState or maxDp (for defaulted values) for all ...
Definition: thresholdPressures.hpp:383
std::vector< double > thresholdPressures(const Deck &, const EclipseState &eclipseState, const Grid &grid, const std::map< std::pair< int, int >, double > &maxDp)
Get a vector of pressure thresholds from EclipseState.
Definition: thresholdPressures.hpp:321