All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
FlowMainEbos.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_FLOW_MAIN_EBOS_HEADER_INCLUDED
24 #define OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
25 
26 
27 #include <sys/utsname.h>
28 
29 
30 #include <opm/simulators/ParallelFileMerger.hpp>
31 #include <opm/simulators/ensureDirectoryExists.hpp>
32 
33 #include <opm/autodiff/BlackoilModelEbos.hpp>
34 #include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
35 #include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
36 #include <opm/autodiff/NewtonIterationBlackoilInterleaved.hpp>
37 #include <opm/autodiff/MissingFeatures.hpp>
38 #include <opm/autodiff/moduleVersion.hpp>
39 #include <opm/autodiff/ExtractParallelGridInformationToISTL.hpp>
40 #include <opm/autodiff/RedistributeDataHandles.hpp>
41 #include <opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp>
42 
43 #include <opm/core/props/satfunc/RelpermDiagnostics.hpp>
44 
45 #include <opm/common/OpmLog/OpmLog.hpp>
46 #include <opm/common/OpmLog/EclipsePRTLog.hpp>
47 #include <opm/common/OpmLog/LogUtil.hpp>
48 
49 #include <opm/parser/eclipse/Deck/Deck.hpp>
50 #include <opm/parser/eclipse/Parser/Parser.hpp>
51 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
52 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
53 #include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
54 #include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
55 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
56 
57 #if HAVE_DUNE_FEM
58 #include <dune/fem/misc/mpimanager.hh>
59 #else
60 #include <dune/common/parallel/mpihelper.hh>
61 #endif
62 
63 namespace Opm
64 {
65  // The FlowMain class is the ebos based black-oil simulator.
66  template <class TypeTag>
68  {
69  enum FileOutputValue{
71  OUTPUT_NONE = 0,
73  OUTPUT_LOG_ONLY = 1,
75  OUTPUT_ALL = 3
76  };
77 
78  public:
79  typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager MaterialLawManager;
80  typedef typename GET_PROP_TYPE(TypeTag, Simulator) EbosSimulator;
81  typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
82  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
83  typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
84  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
85  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
86 
88  typedef typename Simulator::ReservoirState ReservoirState;
89  typedef typename Simulator::OutputWriter OutputWriter;
90 
96  int execute(int argc, char** argv)
97  {
98  try {
99  setupParallelism();
100  printStartupMessage();
101  const bool ok = setupParameters(argc, argv);
102  if (!ok) {
103  return EXIT_FAILURE;
104  }
105 
106  setupEbosSimulator();
107  setupOutput();
108  setupLogging();
109  printPRTHeader();
110  extractMessages();
111  runDiagnostics();
112  setupState();
113  writeInit();
114  setupOutputWriter();
115  setupLinearSolver();
116  createSimulator();
117 
118  // Run.
119  auto ret = runSimulator();
120 
121  mergeParallelLogFiles();
122 
123  return ret;
124  }
125  catch (const std::exception &e) {
126  std::ostringstream message;
127  message << "Program threw an exception: " << e.what();
128 
129  if( output_cout_ )
130  {
131  // in some cases exceptions are thrown before the logging system is set
132  // up.
133  if (OpmLog::hasBackend("STREAMLOG")) {
134  OpmLog::error(message.str());
135  }
136  else {
137  std::cout << message.str() << "\n";
138  }
139  }
140 
141  return EXIT_FAILURE;
142  }
143  }
144 
145  protected:
146  void setupParallelism()
147  {
148  // determine the rank of the current process and the number of processes
149  // involved in the simulation. MPI must have already been initialized here.
150 #if HAVE_MPI
151  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank_);
152  int mpi_size;
153  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
154 #else
155  mpi_rank_ = 0;
156  const int mpi_size = 1;
157 #endif
158  output_cout_ = ( mpi_rank_ == 0 );
159  must_distribute_ = ( mpi_size > 1 );
160 
161 #ifdef _OPENMP
162  // OpenMP setup.
163  if (!getenv("OMP_NUM_THREADS")) {
164  // Default to at most 4 threads, regardless of
165  // number of cores (unless ENV(OMP_NUM_THREADS) is defined)
166  int num_cores = omp_get_num_procs();
167  int num_threads = std::min(4, num_cores);
168  omp_set_num_threads(num_threads);
169  }
170 #pragma omp parallel
171  if (omp_get_thread_num() == 0) {
172  // omp_get_num_threads() only works as expected within a parallel region.
173  const int num_omp_threads = omp_get_num_threads();
174  if (mpi_size == 1) {
175  std::cout << "OpenMP using " << num_omp_threads << " threads." << std::endl;
176  } else {
177  std::cout << "OpenMP using " << num_omp_threads << " threads on MPI rank " << mpi_rank_ << "." << std::endl;
178  }
179  }
180 #endif
181  }
182 
183  // Print startup message if on output rank.
184  void printStartupMessage()
185  {
186 
187  if (output_cout_) {
188  const int lineLen = 70;
189  const std::string version = moduleVersionName();
190  const std::string banner = "This is flow "+version;
191  const int bannerPreLen = (lineLen - 2 - banner.size())/2;
192  const int bannerPostLen = bannerPreLen + (lineLen - 2 - banner.size())%2;
193  std::cout << "**********************************************************************\n";
194  std::cout << "* *\n";
195  std::cout << "*" << std::string(bannerPreLen, ' ') << banner << std::string(bannerPostLen, ' ') << "*\n";
196  std::cout << "* *\n";
197  std::cout << "* Flow is a simulator for fully implicit three-phase black-oil flow, *\n";
198  std::cout << "* including solvent and polymer capabilities. *\n";
199  std::cout << "* For more information, see http://opm-project.org *\n";
200  std::cout << "* *\n";
201  std::cout << "**********************************************************************\n\n";
202  }
203  }
204 
205  // Read parameters, see if a deck was specified on the command line, and if
206  // it was, insert it into parameters.
207  // Writes to:
208  // param_
209  // Returns true if ok, false if not.
210  bool setupParameters(int argc, char** argv)
211  {
212  param_ = ParameterGroup(argc, argv, false, output_cout_);
213 
214  // See if a deck was specified on the command line.
215  if (!param_.unhandledArguments().empty()) {
216  if (param_.unhandledArguments().size() != 1) {
217  std::cerr << "You can only specify a single input deck on the command line.\n";
218  return false;
219  } else {
220  const auto casename = this->simulationCaseName( param_.unhandledArguments()[ 0 ] );
221  param_.insertParameter("deck_filename", casename.string() );
222  }
223  }
224 
225  // We must have an input deck. Grid and props will be read from that.
226  if (!param_.has("deck_filename")) {
227  std::cerr << "This program must be run with an input deck.\n"
228  "Specify the deck filename either\n"
229  " a) as a command line argument by itself\n"
230  " b) as a command line parameter with the syntax deck_filename=<path to your deck>, or\n"
231  " c) as a parameter in a parameter file (.param or .xml) passed to the program.\n";
232  return false;
233  }
234  return true;
235  }
236 
237  // Set output_to_files_ and set/create output dir. Write parameter file.
238  // Writes to:
239  // output_to_files_
240  // output_dir_
241  // Throws std::runtime_error if failed to create (if requested) output dir.
242  void setupOutput()
243  {
244  const std::string output = param_.getDefault("output", std::string("all"));
245  static std::map<std::string, FileOutputValue> string2OutputEnum =
246  { {"none", OUTPUT_NONE },
247  {"false", OUTPUT_LOG_ONLY },
248  {"log", OUTPUT_LOG_ONLY },
249  {"all" , OUTPUT_ALL },
250  {"true" , OUTPUT_ALL }};
251  auto converted = string2OutputEnum.find(output);
252  if ( converted != string2OutputEnum.end() )
253  {
254  output_ = string2OutputEnum[output];
255  }
256  else
257  {
258  std::cerr << "Value " << output <<
259  " passed to option output was invalid. Using \"all\" instead."
260  << std::endl;
261  }
262 
263  output_to_files_ = output_cout_ && output_ > OUTPUT_NONE;
264 
265  // Setup output directory.
266  auto& ioConfig = eclState().getIOConfig();
267  // Default output directory is the directory where the deck is found.
268  const std::string default_output_dir = ioConfig.getOutputDir();
269  output_dir_ = param_.getDefault("output_dir", default_output_dir);
270  // Override output directory if user specified.
271  ioConfig.setOutputDir(output_dir_);
272 
273  // Write parameters used for later reference. (only if rank is zero)
274  if (output_to_files_) {
275  // Create output directory if needed.
276  ensureDirectoryExists(output_dir_);
277  // Write simulation parameters.
278  param_.writeParam(output_dir_ + "/simulation.param");
279  }
280  }
281 
282  // Setup OpmLog backend with output_dir.
283  void setupLogging()
284  {
285  std::string deck_filename = param_.get<std::string>("deck_filename");
286  // create logFile
287  using boost::filesystem::path;
288  path fpath(deck_filename);
289  std::string baseName;
290  std::ostringstream debugFileStream;
291  std::ostringstream logFileStream;
292 
293  if (boost::to_upper_copy(path(fpath.extension()).string()) == ".DATA") {
294  baseName = path(fpath.stem()).string();
295  } else {
296  baseName = path(fpath.filename()).string();
297  }
298 
299  logFileStream << output_dir_ << "/" << baseName;
300  debugFileStream << output_dir_ << "/" << "." << baseName;
301 
302  if ( must_distribute_ && mpi_rank_ != 0 )
303  {
304  // Added rank to log file for non-zero ranks.
305  // This prevents message loss.
306  debugFileStream << "."<< mpi_rank_;
307  // If the following file appears then there is a bug.
308  logFileStream << "." << mpi_rank_;
309  }
310  logFileStream << ".PRT";
311  debugFileStream << ".DEBUG";
312 
313  logFile_ = logFileStream.str();
314 
315  if( output_ > OUTPUT_NONE)
316  {
317  std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::NoDebugMessageTypes, false, output_cout_);
318  OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
319  prtLog->setMessageLimiter(std::make_shared<MessageLimiter>());
320  prtLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(false));
321  }
322 
323  if( output_ >= OUTPUT_LOG_ONLY && !param_.getDefault("no_debug_log", false) )
324  {
325  std::string debugFile = debugFileStream.str();
326  std::shared_ptr<StreamLog> debugLog = std::make_shared<EclipsePRTLog>(debugFile, Log::DefaultMessageTypes, false, output_cout_);
327  OpmLog::addBackend( "DEBUGLOG" , debugLog);
328  }
329 
330  std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(std::cout, Log::StdoutMessageTypes);
331  OpmLog::addBackend( "STREAMLOG", streamLog);
332  const auto& msgLimits = eclState().getSchedule().getMessageLimits();
333  const std::map<int64_t, int> limits = {{Log::MessageType::Note, msgLimits.getCommentPrintLimit(0)},
334  {Log::MessageType::Info, msgLimits.getMessagePrintLimit(0)},
335  {Log::MessageType::Warning, msgLimits.getWarningPrintLimit(0)},
336  {Log::MessageType::Error, msgLimits.getErrorPrintLimit(0)},
337  {Log::MessageType::Problem, msgLimits.getProblemPrintLimit(0)},
338  {Log::MessageType::Bug, msgLimits.getBugPrintLimit(0)}};
339  streamLog->setMessageLimiter(std::make_shared<MessageLimiter>(10, limits));
340  streamLog->setMessageFormatter(std::make_shared<SimpleMessageFormatter>(true));
341 
342  if ( output_cout_ )
343  {
344  // Read Parameters.
345  OpmLog::debug("\n--------------- Reading parameters ---------------\n");
346  }
347  }
348 
349  void printPRTHeader()
350  {
351  // Print header for PRT file.
352  if ( output_cout_ ) {
353  const std::string version = moduleVersionName();
354  const double megabyte = 1024 * 1024;
355  unsigned num_cpu = std::thread::hardware_concurrency();
356  struct utsname arch;
357  const char* user = getlogin();
358  time_t now = std::time(0);
359  struct tm tstruct;
360  char tmstr[80];
361  tstruct = *localtime(&now);
362  strftime(tmstr, sizeof(tmstr), "%d-%m-%Y at %X", &tstruct);
363  const double mem_size = getTotalSystemMemory() / megabyte;
364  std::ostringstream ss;
365  ss << "\n\n\n ######## # ###### # #\n";
366  ss << " # # # # # # \n";
367  ss << " ##### # # # # # # \n";
368  ss << " # # # # # # # # \n";
369  ss << " # ####### ###### # # \n\n";
370  ss << "Flow is a simulator for fully implicit three-phase black-oil flow,";
371  ss << " and is part of OPM.\nFor more information visit: http://opm-project.org \n\n";
372  ss << "Flow Version = " + version + "\n";
373  if (uname(&arch) == 0) {
374  ss << "System = " << arch.nodename << " (Number of cores: " << num_cpu;
375  ss << ", RAM: " << std::fixed << std::setprecision (2) << mem_size << " MB) \n";
376  ss << "Architecture = " << arch.sysname << " " << arch.machine << " (Release: " << arch.release;
377  ss << ", Version: " << arch.version << " )\n";
378  }
379  if (user) {
380  ss << "User = " << user << std::endl;
381  }
382  ss << "Simulation started on " << tmstr << " hrs\n";
383  OpmLog::note(ss.str());
384  }
385  }
386 
387  void mergeParallelLogFiles()
388  {
389  // force closing of all log files.
390  OpmLog::removeAllBackends();
391 
392  if( mpi_rank_ != 0 || !must_distribute_ || !output_to_files_ )
393  {
394  return;
395  }
396 
397  namespace fs = boost::filesystem;
398  fs::path output_path(".");
399  if ( param_.has("output_dir") )
400  {
401  output_path = fs::path(output_dir_);
402  }
403 
404  fs::path deck_filename(param_.get<std::string>("deck_filename"));
405 
406  std::for_each(fs::directory_iterator(output_path),
407  fs::directory_iterator(),
408  detail::ParallelFileMerger(output_path, deck_filename.stem().string()));
409  }
410 
411  void setupEbosSimulator()
412  {
413  std::string progName("flow_ebos");
414  std::string deckFile("--ecl-deck-file-name=");
415  deckFile += param_.get<std::string>("deck_filename");
416  char* ptr[2];
417  ptr[ 0 ] = const_cast< char * > (progName.c_str());
418  ptr[ 1 ] = const_cast< char * > (deckFile.c_str());
419  EbosSimulator::registerParameters();
420  Ewoms::setupParameters_< TypeTag > ( 2, ptr );
421  ebosSimulator_.reset(new EbosSimulator(/*verbose=*/false));
422  ebosSimulator_->model().applyInitialSolution();
423 
424  // Create a grid with a global view.
425  globalGrid_.reset(new Grid(grid()));
426  globalGrid_->switchToGlobalView();
427 
428  try {
429  if (output_cout_) {
430  MissingFeatures::checkKeywords(deck());
431  }
432 
433  // Possible to force initialization only behavior (NOSIM).
434  if (param_.has("nosim")) {
435  const bool nosim = param_.get<bool>("nosim");
436  auto& ioConfig = eclState().getIOConfig();
437  ioConfig.overrideNOSIM( nosim );
438  }
439  }
440  catch (const std::invalid_argument& e) {
441  std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile_ << std::endl;
442  std::cerr << "Exception caught: " << e.what() << std::endl;
443  throw;
444  }
445 
446  // Possibly override IOConfig setting (from deck) for how often RESTART files should get written to disk (every N report step)
447  if (param_.has("output_interval")) {
448  const int output_interval = param_.get<int>("output_interval");
449  eclState().getRestartConfig().overrideRestartWriteInterval( size_t( output_interval ) );
450  }
451  }
452 
453  const Deck& deck() const
454  { return ebosSimulator_->gridManager().deck(); }
455 
456  Deck& deck()
457  { return ebosSimulator_->gridManager().deck(); }
458 
459  const EclipseState& eclState() const
460  { return ebosSimulator_->gridManager().eclState(); }
461 
462  EclipseState& eclState()
463  { return ebosSimulator_->gridManager().eclState(); }
464 
465  // Initialise the reservoir state. Updated fluid props for SWATINIT.
466  // Writes to:
467  // state_
468  // threshold_pressures_
469  void setupState()
470  {
471  const PhaseUsage pu = Opm::phaseUsageFromDeck(deck());
472  const Grid& grid = this->grid();
473 
474  // Need old-style fluid object for init purposes (only).
475  BlackoilPropertiesFromDeck props(deck(),
476  eclState(),
477  materialLawManager(),
478  grid.size(/*codim=*/0),
479  grid.globalCell().data(),
480  grid.logicalCartesianSize().data(),
481  param_);
482 
483 
484  // Init state variables (saturation and pressure).
485  if (param_.has("init_saturation")) {
486  state_.reset(new ReservoirState(grid.size(/*codim=*/0),
487  grid.numFaces(),
488  props.numPhases()));
489 
490  initStateBasic(grid.size(/*codim=*/0),
491  grid.globalCell().data(),
492  grid.logicalCartesianSize().data(),
493  grid.numFaces(),
494  Opm::UgGridHelpers::faceCells(grid),
495  Opm::UgGridHelpers::beginFaceCentroids(grid),
496  Opm::UgGridHelpers::beginCellCentroids(grid),
497  Grid::dimension,
498  props, param_, gravity(), *state_);
499 
500  initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_);
501 
502  enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
503  if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
504  const int numPhases = props.numPhases();
505  const int numCells = Opm::UgGridHelpers::numCells(grid);
506 
507  // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState.
508  auto& gor = state_->getCellData( BlackoilState::GASOILRATIO );
509  const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL );
510  for (int c = 0; c < numCells; ++c) {
511  // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field.
512  gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]];
513  }
514  }
515  } else if (deck().hasKeyword("EQUIL")) {
516  // Which state class are we really using - what a f... mess?
517  state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
518  Opm::UgGridHelpers::numFaces(grid),
519  props.numPhases()));
520 
521  initStateEquil(grid, props, deck(), eclState(), gravity(), *state_);
522  //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
523  } else {
524  state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
525  Opm::UgGridHelpers::numFaces(grid),
526  props.numPhases()));
527  initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
528  Opm::UgGridHelpers::globalCell(grid),
529  Opm::UgGridHelpers::numFaces(grid),
530  Opm::UgGridHelpers::faceCells(grid),
531  Opm::UgGridHelpers::beginFaceCentroids(grid),
532  Opm::UgGridHelpers::beginCellCentroids(grid),
533  Opm::UgGridHelpers::dimensions(grid),
534  props, deck(), gravity(), *state_);
535  }
536 
537  initHydroCarbonState(*state_, pu, Opm::UgGridHelpers::numCells(grid), deck().hasKeyword("DISGAS"), deck().hasKeyword("VAPOIL"));
538 
539  // Get initial polymer concentration from ebos
540  if (GET_PROP_VALUE(TypeTag, EnablePolymer)) {
541  auto& cpolymer = state_->getCellData( state_->POLYMER );
542  const int numCells = Opm::UgGridHelpers::numCells(grid);
543  for (int c = 0; c < numCells; ++c) {
544  cpolymer[c] = ebosProblem().polymerConcentration(c);
545  }
546  }
547  // Get initial solvent saturation from ebos
548  if (GET_PROP_VALUE(TypeTag, EnableSolvent)) {
549  auto& solvent = state_->getCellData( state_->SSOL );
550  auto& sat = state_->saturation();
551  const int np = props.numPhases();
552  const int numCells = Opm::UgGridHelpers::numCells(grid);
553  for (int c = 0; c < numCells; ++c) {
554  solvent[c] = ebosProblem().solventSaturation(c);
555  sat[c * np + pu.phase_pos[Water]];
556  }
557  }
558 
559  }
560 
561  // Extract messages from parser.
562  // Writes to:
563  // OpmLog singleton.
564  void extractMessages()
565  {
566  if ( !output_cout_ )
567  {
568  return;
569  }
570 
571  auto extractMessage = [this](const Message& msg) {
572  auto log_type = this->convertMessageType(msg.mtype);
573  const auto& location = msg.location;
574  if (location) {
575  OpmLog::addMessage(log_type, Log::fileMessage(location.filename, location.lineno, msg.message));
576  } else {
577  OpmLog::addMessage(log_type, msg.message);
578  }
579  };
580 
581  // Extract messages from Deck.
582  for(const auto& msg : deck().getMessageContainer()) {
583  extractMessage(msg);
584  }
585 
586  // Extract messages from EclipseState.
587  for (const auto& msg : eclState().getMessageContainer()) {
588  extractMessage(msg);
589  }
590  }
591 
592  // Run diagnostics.
593  // Writes to:
594  // OpmLog singleton.
595  void runDiagnostics()
596  {
597  if( ! output_cout_ )
598  {
599  return;
600  }
601 
602  // Run relperm diagnostics
603  RelpermDiagnostics diagnostic;
604  diagnostic.diagnosis(eclState(), deck(), this->grid());
605  }
606 
607  void writeInit()
608  {
609  bool output = ( output_ > OUTPUT_LOG_ONLY );
610  bool output_ecl = param_.getDefault("output_ecl", true);
611  if( output && output_ecl && grid().comm().rank() == 0 )
612  {
613  exportNncStructure_();
614 
615  const EclipseGrid& inputGrid = eclState().getInputGrid();
616  eclIO_.reset(new EclipseIO(eclState(), UgGridHelpers::createEclipseGrid( this->globalGrid() , inputGrid )));
617  eclIO_->writeInitial(computeLegacySimProps_(), nnc_);
618  }
619  }
620 
621  // Setup output writer.
622  // Writes to:
623  // output_writer_
624  void setupOutputWriter()
625  {
626  // create output writer after grid is distributed, otherwise the parallel output
627  // won't work correctly since we need to create a mapping from the distributed to
628  // the global view
629  output_writer_.reset(new OutputWriter(grid(),
630  param_,
631  eclState(),
632  std::move(eclIO_),
633  Opm::phaseUsageFromDeck(deck())) );
634  }
635 
636  // Run the simulator.
637  // Returns EXIT_SUCCESS if it does not throw.
638  int runSimulator()
639  {
640  const auto& schedule = eclState().getSchedule();
641  const auto& timeMap = schedule.getTimeMap();
642  auto& ioConfig = eclState().getIOConfig();
643  SimulatorTimer simtimer;
644 
645  // initialize variables
646  const auto& initConfig = eclState().getInitConfig();
647  simtimer.init(timeMap, (size_t)initConfig.getRestartStep());
648 
649  if (!ioConfig.initOnly()) {
650  if (output_cout_) {
651  std::string msg;
652  msg = "\n\n================ Starting main simulation loop ===============\n";
653  OpmLog::info(msg);
654  }
655 
656  SimulatorReport successReport = simulator_->run(simtimer, *state_);
657  SimulatorReport failureReport = simulator_->failureReport();
658 
659  if (output_cout_) {
660  std::ostringstream ss;
661  ss << "\n\n================ End of simulation ===============\n\n";
662  successReport.reportFullyImplicit(ss, &failureReport);
663  OpmLog::info(ss.str());
664  if (param_.anyUnused()) {
665  // This allows a user to catch typos and misunderstandings in the
666  // use of simulator parameters.
667  std::cout << "-------------------- Unused parameters: --------------------\n";
668  param_.displayUsage();
669  std::cout << "----------------------------------------------------------------" << std::endl;
670  }
671  }
672 
673  if (output_to_files_) {
674  std::string filename = output_dir_ + "/walltime.txt";
675  std::fstream tot_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
676  successReport.reportParam(tot_os);
677  }
678  } else {
679  if (output_cout_) {
680  std::cout << "\n\n================ Simulation turned off ===============\n" << std::flush;
681  }
682 
683  }
684  return EXIT_SUCCESS;
685  }
686 
687  // Setup linear solver.
688  // Writes to:
689  // fis_solver_
690  void setupLinearSolver()
691  {
692  typedef typename BlackoilModelEbos<TypeTag> :: ISTLSolverType ISTLSolverType;
693 
694  extractParallelGridInformationToISTL(grid(), parallel_information_);
695  fis_solver_.reset( new ISTLSolverType( param_, parallel_information_ ) );
696  }
697 
699  // Create simulator instance.
700  // Writes to:
701  // simulator_
703  {
704  // Create the simulator instance.
705  simulator_.reset(new Simulator(*ebosSimulator_,
706  param_,
707  *fis_solver_,
708  FluidSystem::enableDissolvedGas(),
709  FluidSystem::enableVaporizedOil(),
710  eclState(),
711  *output_writer_,
712  defunctWellNames()));
713  }
714 
715  private:
716  boost::filesystem::path simulationCaseName( const std::string& casename ) {
717  namespace fs = boost::filesystem;
718 
719  const auto exists = []( const fs::path& f ) -> bool {
720  if( !fs::exists( f ) ) return false;
721 
722  if( fs::is_regular_file( f ) ) return true;
723 
724  return fs::is_symlink( f )
725  && fs::is_regular_file( fs::read_symlink( f ) );
726  };
727 
728  auto simcase = fs::path( casename );
729 
730  if( exists( simcase ) ) {
731  return simcase;
732  }
733 
734  for( const auto& ext : { std::string("data"), std::string("DATA") } ) {
735  if( exists( simcase.replace_extension( ext ) ) ) {
736  return simcase;
737  }
738  }
739 
740  throw std::invalid_argument( "Cannot find input case " + casename );
741  }
742 
743  unsigned long long getTotalSystemMemory()
744  {
745  long pages = sysconf(_SC_PHYS_PAGES);
746  long page_size = sysconf(_SC_PAGE_SIZE);
747  return pages * page_size;
748  }
749 
750  int64_t convertMessageType(const Message::type& mtype)
751  {
752  switch (mtype) {
753  case Message::type::Debug:
754  return Log::MessageType::Debug;
755  case Message::type::Info:
756  return Log::MessageType::Info;
757  case Message::type::Warning:
758  return Log::MessageType::Warning;
759  case Message::type::Error:
760  return Log::MessageType::Error;
761  case Message::type::Problem:
762  return Log::MessageType::Problem;
763  case Message::type::Bug:
764  return Log::MessageType::Bug;
765  case Message::type::Note:
766  return Log::MessageType::Note;
767  }
768  throw std::logic_error("Invalid messages type!\n");
769  }
770 
771  Grid& grid()
772  { return ebosSimulator_->gridManager().grid(); }
773 
774  const Grid& globalGrid()
775  { return *globalGrid_; }
776 
777  Problem& ebosProblem()
778  { return ebosSimulator_->problem(); }
779 
780  const Problem& ebosProblem() const
781  { return ebosSimulator_->problem(); }
782 
783  std::shared_ptr<MaterialLawManager> materialLawManager()
784  { return ebosProblem().materialLawManager(); }
785 
786  Scalar gravity() const
787  { return ebosProblem().gravity()[2]; }
788 
789  std::unordered_set<std::string> defunctWellNames() const
790  { return ebosSimulator_->gridManager().defunctWellNames(); }
791 
792  data::Solution computeLegacySimProps_()
793  {
794  const int* dims = UgGridHelpers::cartDims(grid());
795  const int globalSize = dims[0]*dims[1]*dims[2];
796 
797  data::CellData tranx = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
798  data::CellData trany = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
799  data::CellData tranz = {UnitSystem::measure::transmissibility, std::vector<double>( globalSize ), data::TargetType::INIT};
800 
801  for (size_t i = 0; i < tranx.data.size(); ++i) {
802  tranx.data[0] = 0.0;
803  trany.data[0] = 0.0;
804  tranz.data[0] = 0.0;
805  }
806 
807  const Grid& globalGrid = this->globalGrid();
808  const auto& globalGridView = globalGrid.leafGridView();
809  typedef typename Grid::LeafGridView GridView;
810  typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
811  ElementMapper globalElemMapper(globalGridView);
812  const auto& cartesianCellIdx = globalGrid.globalCell();
813 
814  const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
815  if (grid().comm().size() < 2) {
816  // in the sequential case we must use the transmissibilites defined by
817  // the problem. (because in the sequential case, the grid manager does
818  // not compute "global" transmissibilities for performance reasons. in
819  // the parallel case, the problem's transmissibilities can't be used
820  // because this object refers to the distributed grid and we need the
821  // sequential version here.)
822  globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
823  }
824 
825  auto elemIt = globalGridView.template begin</*codim=*/0>();
826  const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
827  for (; elemIt != elemEndIt; ++ elemIt) {
828  const auto& elem = *elemIt;
829 
830  auto isIt = globalGridView.ibegin(elem);
831  const auto& isEndIt = globalGridView.iend(elem);
832  for (; isIt != isEndIt; ++ isIt) {
833  const auto& is = *isIt;
834 
835  if (!is.neighbor())
836  {
837  continue; // intersection is on the domain boundary
838  }
839 
840  unsigned c1 = globalElemMapper.index(is.inside());
841  unsigned c2 = globalElemMapper.index(is.outside());
842 
843  if (c1 > c2)
844  {
845  continue; // we only need to handle each connection once, thank you.
846  }
847 
848 
849  int gc1 = std::min(cartesianCellIdx[c1], cartesianCellIdx[c2]);
850  int gc2 = std::max(cartesianCellIdx[c1], cartesianCellIdx[c2]);
851  if (gc2 - gc1 == 1) {
852  tranx.data[gc1] = globalTrans->transmissibility(c1, c2);
853  }
854 
855  if (gc2 - gc1 == dims[0]) {
856  trany.data[gc1] = globalTrans->transmissibility(c1, c2);
857  }
858 
859  if (gc2 - gc1 == dims[0]*dims[1]) {
860  tranz.data[gc1] = globalTrans->transmissibility(c1, c2);
861  }
862  }
863  }
864 
865  return {{"TRANX" , tranx},
866  {"TRANY" , trany} ,
867  {"TRANZ" , tranz}};
868  }
869 
870  void exportNncStructure_()
871  {
872  nnc_ = eclState().getInputNNC();
873  int nx = eclState().getInputGrid().getNX();
874  int ny = eclState().getInputGrid().getNY();
875  //int nz = eclState().getInputGrid().getNZ()
876 
877  const Grid& globalGrid = this->globalGrid();
878  const auto& globalGridView = globalGrid.leafGridView();
879  typedef typename Grid::LeafGridView GridView;
880  typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
881  ElementMapper globalElemMapper(globalGridView);
882 
883  const auto* globalTrans = &(ebosSimulator_->gridManager().globalTransmissibility());
884  if (grid().comm().size() < 2) {
885  // in the sequential case we must use the transmissibilites defined by
886  // the problem. (because in the sequential case, the grid manager does
887  // not compute "global" transmissibilities for performance reasons. in
888  // the parallel case, the problem's transmissibilities can't be used
889  // because this object refers to the distributed grid and we need the
890  // sequential version here.)
891  globalTrans = &ebosSimulator_->problem().eclTransmissibilities();
892  }
893 
894  auto elemIt = globalGridView.template begin</*codim=*/0>();
895  const auto& elemEndIt = globalGridView.template end</*codim=*/0>();
896  for (; elemIt != elemEndIt; ++ elemIt) {
897  const auto& elem = *elemIt;
898 
899  auto isIt = globalGridView.ibegin(elem);
900  const auto& isEndIt = globalGridView.iend(elem);
901  for (; isIt != isEndIt; ++ isIt) {
902  const auto& is = *isIt;
903 
904  if (!is.neighbor())
905  {
906  continue; // intersection is on the domain boundary
907  }
908 
909  unsigned c1 = globalElemMapper.index(is.inside());
910  unsigned c2 = globalElemMapper.index(is.outside());
911 
912  if (c1 > c2)
913  {
914  continue; // we only need to handle each connection once, thank you.
915  }
916 
917  // TODO (?): use the cartesian index mapper to make this code work
918  // with grids other than Dune::CpGrid. The problem is that we need
919  // the a mapper for the sequential grid, not for the distributed one.
920  int cc1 = globalGrid.globalCell()[c1];
921  int cc2 = globalGrid.globalCell()[c2];
922 
923  if (std::abs(cc1 - cc2) != 1 &&
924  std::abs(cc1 - cc2) != nx &&
925  std::abs(cc1 - cc2) != nx*ny)
926  {
927  nnc_.addNNC(cc1, cc2, globalTrans->transmissibility(c1, c2));
928  }
929  }
930  }
931  }
932 
933  std::unique_ptr<EbosSimulator> ebosSimulator_;
934  int mpi_rank_ = 0;
935  bool output_cout_ = false;
936  FileOutputValue output_ = OUTPUT_ALL;
937  bool must_distribute_ = false;
938  ParameterGroup param_;
939  bool output_to_files_ = false;
940  std::string output_dir_ = std::string(".");
941  std::unique_ptr<ReservoirState> state_;
942  NNC nnc_;
943  std::unique_ptr<EclipseIO> eclIO_;
944  std::unique_ptr<OutputWriter> output_writer_;
945  boost::any parallel_information_;
946  std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver_;
947  std::unique_ptr<Simulator> simulator_;
948  std::string logFile_;
949  // Needs to be shared pointer because it gets initialzed before MPI_Init.
950  std::shared_ptr<Grid> globalGrid_;
951  };
952 } // namespace Opm
953 
954 #endif // OPM_FLOW_MAIN_EBOS_HEADER_INCLUDED
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
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
a simulator for the blackoil model
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:48
Definition: FlowMainEbos.hpp:67
void createSimulator()
This is the main function of Flow.
Definition: FlowMainEbos.hpp:702
int execute(int argc, char **argv)
This is the main function of Flow.
Definition: FlowMainEbos.hpp:96
Wrapper class for VTK, Matlab, and ECL output.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:206