All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
SimulatorFullyImplicitBlackoilOutput.hpp
1 /*
2  Copyright (c) 2014 SINTEF ICT, Applied Mathematics.
3  Copyright (c) 2015 IRIS AS
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 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
21 #define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
22 #include <opm/core/grid.h>
23 #include <opm/simulators/timestepping/SimulatorTimerInterface.hpp>
24 #include <opm/core/simulator/WellState.hpp>
25 #include <opm/autodiff/Compat.hpp>
26 #include <opm/core/utility/DataMap.hpp>
27 #include <opm/common/ErrorMacros.hpp>
28 #include <opm/common/OpmLog/OpmLog.hpp>
29 #include <opm/core/utility/miscUtilities.hpp>
30 #include <opm/core/utility/parameters/ParameterGroup.hpp>
31 #include <opm/core/wells/DynamicListEconLimited.hpp>
32 #include <opm/core/simulator/BlackoilState.hpp>
33 #include <opm/core/simulator/SimulatorReport.hpp>
34 
35 #include <opm/output/data/Cells.hpp>
36 #include <opm/output/data/Solution.hpp>
37 #include <opm/output/eclipse/EclipseIO.hpp>
38 
39 #include <opm/autodiff/GridHelpers.hpp>
40 #include <opm/autodiff/ParallelDebugOutput.hpp>
41 
42 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
43 #include <opm/autodiff/ThreadHandle.hpp>
44 #include <opm/autodiff/AutoDiffBlock.hpp>
45 
46 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
47 #include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
48 #include <opm/simulators/ensureDirectoryExists.hpp>
49 
50 #include <string>
51 #include <sstream>
52 #include <iomanip>
53 #include <fstream>
54 #include <thread>
55 #include <map>
56 
57 #include <boost/filesystem.hpp>
58 
59 #ifdef HAVE_OPM_GRID
60 #include <dune/grid/CpGrid.hpp>
61 #endif
62 namespace Opm
63 {
64 
65  class SimulationDataContainer;
66  class BlackoilState;
67 
68  void outputStateVtk(const UnstructuredGrid& grid,
69  const Opm::SimulationDataContainer& state,
70  const int step,
71  const std::string& output_dir);
72 
73  void outputWellStateMatlab(const Opm::WellState& well_state,
74  const int step,
75  const std::string& output_dir);
76 #ifdef HAVE_OPM_GRID
77  void outputStateVtk(const Dune::CpGrid& grid,
78  const Opm::SimulationDataContainer& state,
79  const int step,
80  const std::string& output_dir);
81 #endif
82 
83  template<class Grid>
84  void outputStateMatlab(const Grid& grid,
85  const Opm::SimulationDataContainer& state,
86  const int step,
87  const std::string& output_dir)
88  {
89  Opm::DataMap dm;
90  dm["saturation"] = &state.saturation();
91  dm["pressure"] = &state.pressure();
92  for (const auto& pair : state.cellData())
93  {
94  const std::string& name = pair.first;
95  std::string key;
96  if( name == "SURFACEVOL" ) {
97  key = "surfvolume";
98  }
99  else if( name == "RV" ) {
100  key = "rv";
101  }
102  else if( name == "GASOILRATIO" ) {
103  key = "rs";
104  }
105  else { // otherwise skip entry
106  continue;
107  }
108  // set data to datmap
109  dm[ key ] = &pair.second;
110  }
111 
112  std::vector<double> cell_velocity;
113  Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
114  AutoDiffGrid::numFaces(grid),
115  AutoDiffGrid::beginFaceCentroids(grid),
116  UgGridHelpers::faceCells(grid),
117  AutoDiffGrid::beginCellCentroids(grid),
118  AutoDiffGrid::beginCellVolumes(grid),
119  AutoDiffGrid::dimensions(grid),
120  state.faceflux(), cell_velocity);
121  dm["velocity"] = &cell_velocity;
122 
123  // Write data (not grid) in Matlab format
124  for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
125  std::ostringstream fname;
126  fname << output_dir << "/" << it->first;
127  ensureDirectoryExists(fname.str());
128  fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
129  std::ofstream file(fname.str().c_str());
130  if (!file) {
131  OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
132  }
133  file.precision(15);
134  const std::vector<double>& d = *(it->second);
135  std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
136  }
137  }
138 
140  public:
141  BlackoilSubWriter( const std::string& outputDir )
142  : outputDir_( outputDir )
143  {}
144 
145  virtual void writeTimeStep(const SimulatorTimerInterface& timer,
146  const SimulationDataContainer& state,
148  bool /*substep*/ = false) = 0;
149  protected:
150  const std::string outputDir_;
151  };
152 
153  template< class Grid >
155  public:
156  BlackoilVTKWriter( const Grid& grid,
157  const std::string& outputDir )
158  : BlackoilSubWriter( outputDir )
159  , grid_( grid )
160  {}
161 
162  void writeTimeStep(const SimulatorTimerInterface& timer,
163  const SimulationDataContainer& state,
165  bool /*substep*/ = false) override
166  {
167  outputStateVtk(grid_, state, timer.currentStepNum(), outputDir_);
168  }
169 
170  protected:
171  const Grid& grid_;
172  };
173 
174  template< typename Grid >
176  {
177  public:
178  BlackoilMatlabWriter( const Grid& grid,
179  const std::string& outputDir )
180  : BlackoilSubWriter( outputDir )
181  , grid_( grid )
182  {}
183 
184  void writeTimeStep(const SimulatorTimerInterface& timer,
185  const SimulationDataContainer& reservoirState,
186  const WellStateFullyImplicitBlackoil& wellState,
187  bool /*substep*/ = false) override
188  {
189  outputStateMatlab(grid_, reservoirState, timer.currentStepNum(), outputDir_);
190  outputWellStateMatlab(wellState, timer.currentStepNum(), outputDir_);
191  }
192 
193  protected:
194  const Grid& grid_;
195  };
196 
197 
199  struct ExtraData
200  {
201  double suggested_step = -1.0;
202  };
203 
204 
207  {
208 
209  public:
210  // constructor creating different sub writers
211  template <class Grid>
212  BlackoilOutputWriter(const Grid& grid,
213  const ParameterGroup& param,
214  const Opm::EclipseState& eclipseState,
215  std::unique_ptr<EclipseIO>&& eclIO,
216  const Opm::PhaseUsage &phaseUsage);
217 
219  void writeInit(const data::Solution& simProps, const NNC& nnc);
220 
227  template<class Model>
228  void writeTimeStep(const SimulatorTimerInterface& timer,
229  const SimulationDataContainer& reservoirState,
230  const Opm::WellStateFullyImplicitBlackoil& wellState,
231  const Model& physicalModel,
232  const bool substep = false,
233  const double nextstep = -1.0,
234  const SimulatorReport& simulatorReport = SimulatorReport());
235 
236 
243  const SimulatorTimerInterface& timer,
244  const SimulationDataContainer& reservoirState,
245  const data::Solution& cellData,
246  const Opm::WellStateFullyImplicitBlackoil& wellState,
247  const std::map<std::string, double>& miscSummaryData,
248  const std::map<std::string, std::vector<double>>& extraRestartData,
249  bool substep = false);
250 
257  const SimulatorTimerInterface& timer,
258  const SimulationDataContainer& reservoirState,
259  const Opm::WellStateFullyImplicitBlackoil& wellState,
260  const std::map<std::string, double>& miscSummaryData,
261  const std::map<std::string, std::vector<double>>& extraRestartData,
262  bool substep = false);
263 
270  const SimulationDataContainer& reservoirState,
271  const Opm::WellStateFullyImplicitBlackoil& wellState,
272  const data::Solution& simProps,
273  const std::map<std::string, double>& miscSummaryData,
274  const std::map<std::string, std::vector<double>>& extraRestartData,
275  bool substep );
276 
278  const std::string& outputDirectory() const { return outputDir_; }
279 
281  bool output () const { return output_; }
282 
284  bool isIORank () const
285  {
286  return parallelOutput_->isIORank();
287  }
288 
289  void restore(SimulatorTimerInterface& timer,
290  BlackoilState& state,
292  const std::string& filename,
293  const int desiredReportStep);
294 
295 
296  template <class Grid, class WellState>
297  void initFromRestartFile(const PhaseUsage& phaseUsage,
298  const Grid& grid,
299  SimulationDataContainer& simulatorstate,
300  WellState& wellstate,
301  ExtraData& extra);
302 
303  bool isRestart() const;
304 
305  bool requireFIPNUM() const;
306 
307  protected:
308  const bool output_;
309  std::unique_ptr< ParallelDebugOutputInterface > parallelOutput_;
310 
311  // Parameters for output.
312  const std::string outputDir_;
313  const bool restart_double_si_;
314 
315  int lastBackupReportStep_;
316 
317  std::ofstream backupfile_;
318  Opm::PhaseUsage phaseUsage_;
319  std::unique_ptr< BlackoilSubWriter > vtkWriter_;
320  std::unique_ptr< BlackoilSubWriter > matlabWriter_;
321  std::unique_ptr< EclipseIO > eclIO_;
322  const EclipseState& eclipseState_;
323 
324  std::unique_ptr< ThreadHandle > asyncOutput_;
325  };
326 
327 
329  //
330  // Implementation
331  //
333  template <class Grid>
334  inline
335  BlackoilOutputWriter::
336  BlackoilOutputWriter(const Grid& grid,
337  const ParameterGroup& param,
338  const Opm::EclipseState& eclipseState,
339  std::unique_ptr<EclipseIO>&& eclIO,
340  const Opm::PhaseUsage &phaseUsage)
341  : output_( [ &param ] () -> bool {
342  // If output parameter is true or all, then we do output
343  const std::string outputString = param.getDefault("output", std::string("all"));
344  return ( outputString == "all" || outputString == "true" );
345  }()
346  ),
347  parallelOutput_( output_ ? new ParallelDebugOutput< Grid >( grid, eclipseState, phaseUsage.num_phases, phaseUsage ) : 0 ),
348  outputDir_( eclipseState.getIOConfig().getOutputDir() ),
349  restart_double_si_( output_ ? param.getDefault("restart_double_si", false) : false ),
350  lastBackupReportStep_( -1 ),
351  phaseUsage_( phaseUsage ),
352  eclipseState_(eclipseState),
353  asyncOutput_()
354  {
355  // For output.
356  if ( output_ )
357  {
358  if ( param.getDefault("output_vtk",false) )
359  {
360  vtkWriter_
361  .reset(new BlackoilVTKWriter< Grid >( grid, outputDir_ ));
362  }
363 
364  auto output_matlab = param.getDefault("output_matlab", false );
365 
366  if ( parallelOutput_->isParallel() && output_matlab )
367  {
368  Opm::OpmLog::warning("Parallel Output Config",
369  "Velocity output for matlab is broken in parallel.");
370  }
371 
372  if( parallelOutput_->isIORank() ) {
373 
374  if ( output_matlab )
375  {
376  matlabWriter_
377  .reset(new BlackoilMatlabWriter< Grid >( grid, outputDir_ ));
378  }
379 
380  eclIO_ = std::move(eclIO);
381 
382  // Ensure that output dir exists
383  ensureDirectoryExists(outputDir_);
384 
385  std::string backupfilename = param.getDefault("backupfile", std::string("") );
386  if( ! backupfilename.empty() )
387  {
388  backupfile_.open( backupfilename.c_str() );
389  }
390  }
391 
392  // create output thread if enabled and rank is I/O rank
393  // async output is enabled by default if pthread are enabled
394 #if HAVE_PTHREAD
395  const bool asyncOutputDefault = true;
396 #else
397  const bool asyncOutputDefault = false;
398 #endif
399  if( param.getDefault("async_output", asyncOutputDefault ) )
400  {
401  const bool isIORank = parallelOutput_ ? parallelOutput_->isIORank() : true;
402 #if HAVE_PTHREAD
403  asyncOutput_.reset( new ThreadHandle( isIORank ) );
404 #else
405  OPM_THROW(std::runtime_error,"Pthreads were not found, cannot enable async_output");
406 #endif
407  }
408  }
409  }
410 
411 
412  template <class Grid, class WellState>
413  inline void
414  BlackoilOutputWriter::
415  initFromRestartFile( const PhaseUsage& phaseUsage,
416  const Grid& grid,
417  SimulationDataContainer& simulatorstate,
418  WellState& wellstate,
419  ExtraData& extra )
420  {
421  std::map<std::string, RestartKey> solution_keys {{"PRESSURE" , RestartKey(UnitSystem::measure::pressure)},
422  {"SWAT" , RestartKey(UnitSystem::measure::identity)},
423  {"SGAS" , RestartKey(UnitSystem::measure::identity)},
424  {"TEMP" , RestartKey(UnitSystem::measure::temperature)},
425  {"RS" , RestartKey(UnitSystem::measure::gas_oil_ratio)},
426  {"RV" , RestartKey(UnitSystem::measure::oil_gas_ratio)},
427  {"SOMAX", {UnitSystem::measure::identity, false}},
428  {"PCSWM_OW", {UnitSystem::measure::identity, false}},
429  {"KRNSW_OW", {UnitSystem::measure::identity, false}},
430  {"PCSWM_GO", {UnitSystem::measure::identity, false}},
431  {"KRNSW_GO", {UnitSystem::measure::identity, false}}};
432 
433  std::map<std::string, bool> extra_keys {
434  {"OPMEXTRA" , false}
435  };
436 
437  if (restart_double_si_) {
438  // Avoid any unit conversions, treat restart input as SI units.
439  for (auto& elem : solution_keys) {
440  elem.second = RestartKey(UnitSystem::measure::identity);
441  }
442  }
443 
444  // gives a dummy dynamic_list_econ_limited
445  DynamicListEconLimited dummy_list_econ_limited;
446  WellsManager wellsmanager(eclipseState_,
447  eclipseState_.getInitConfig().getRestartStep(),
448  Opm::UgGridHelpers::numCells(grid),
449  Opm::UgGridHelpers::globalCell(grid),
450  Opm::UgGridHelpers::cartDims(grid),
451  Opm::UgGridHelpers::dimensions(grid),
452  Opm::UgGridHelpers::cell2Faces(grid),
453  Opm::UgGridHelpers::beginFaceCentroids(grid),
454  dummy_list_econ_limited
455  // We need to pass the optionaly arguments
456  // as we get the following error otherwise
457  // with c++ (Debian 4.9.2-10) 4.9.2 and -std=c++11
458  // converting to ‘const std::unordered_set<std::basic_string<char> >’ from initializer list would use explicit constructo
459  , false,
460  std::unordered_set<std::string>());
461 
462  const Wells* wells = wellsmanager.c_wells();
463  wellstate.resize(wells, simulatorstate, phaseUsage ); //Resize for restart step
464  auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
465 
466  solutionToSim( restart_values.solution, restart_values.extra, phaseUsage, simulatorstate );
467  wellsToState( restart_values.wells, phaseUsage, wellstate );
468 
469  const auto opmextra_iter = restart_values.extra.find("OPMEXTRA");
470  if (opmextra_iter != restart_values.extra.end()) {
471  std::vector<double> opmextra = opmextra_iter->second;
472  assert(opmextra.size() == 1);
473  extra.suggested_step = opmextra[0];
474  } else {
475  OpmLog::warning("Restart data is missing OPMEXTRA field, restart run may deviate from original run.");
476  extra.suggested_step = -1.0;
477  }
478  }
479 
480 
481 
482 
483 
484  namespace detail {
485 
486 
487  template <class V>
488  void addToSimData( SimulationDataContainer& simData,
489  const std::string& name,
490  const V& vec )
491  {
492  if (vec.size() == 0) {
493  return;
494  }
495 
496  typedef std::vector< double > OutputVectorType;
497 
498  // get data map
499  auto& dataMap = simData.cellData();
500 
501  // insert name,vector into data map
502  dataMap.insert( std::make_pair( name, OutputVectorType( vec.data(), vec.data() + vec.size() ) ) );
503  }
504 
505  template <class Scalar>
506  void addToSimData( SimulationDataContainer& simData,
507  const std::string& name,
508  const AutoDiffBlock<Scalar>& adb )
509  {
510  // forward value of ADB to output
511  addToSimData( simData, name, adb.value() );
512  }
513 
514 
515  // this method basically converts all Eigen vectors to std::vectors
516  // stored in a SimulationDataContainer
517  template <class SimulatorData>
518  SimulationDataContainer
519  convertToSimulationDataContainer( const SimulatorData& sd,
520  const SimulationDataContainer& localState,
521  const Opm::PhaseUsage& phaseUsage )
522  {
523  // copy local state and then add missing data
524  SimulationDataContainer simData( localState );
525 
526  //Get shorthands for water, oil, gas
527  const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
528  const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
529  const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
530 
531  const int aqua_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Aqua];
532  const int liquid_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Liquid];
533  const int vapour_idx = phaseUsage.phase_pos[Opm::PhaseUsage::Vapour];
534 
535  // WATER
536  if( aqua_active ) {
537  addToSimData( simData, "1OVERBW", sd.rq[aqua_idx].b );
538  addToSimData( simData, "WAT_DEN", sd.rq[aqua_idx].rho );
539  addToSimData( simData, "WAT_VISC", sd.rq[aqua_idx].mu );
540  addToSimData( simData, "WATKR", sd.rq[aqua_idx].kr );
541  }
542 
543  // OIL
544  if( liquid_active ) {
545  addToSimData( simData, "1OVERBO", sd.rq[liquid_idx].b );
546  addToSimData( simData, "OIL_DEN", sd.rq[liquid_idx].rho );
547  addToSimData( simData, "OIL_VISC", sd.rq[liquid_idx].mu );
548  addToSimData( simData, "OILKR", sd.rq[liquid_idx].kr );
549  }
550 
551  // GAS
552  if( vapour_active ) {
553  addToSimData( simData, "1OVERBG", sd.rq[vapour_idx].b );
554  addToSimData( simData, "GAS_DEN", sd.rq[vapour_idx].rho );
555  addToSimData( simData, "GAS_VISC", sd.rq[vapour_idx].mu );
556  addToSimData( simData, "GASKR", sd.rq[vapour_idx].kr );
557  }
558 
559  // RS and RV
560  addToSimData( simData, "RSSAT", sd.rsSat );
561  addToSimData( simData, "RVSAT", sd.rvSat );
562 
563  addToSimData( simData, "SOMAX", sd.soMax );
564  addToSimData( simData, "PBUB", sd.Pb );
565  addToSimData( simData, "PDEW", sd.Pd );
566  addToSimData( simData, "PCSWMDC_OW", sd.pcswmdc_ow );
567  addToSimData( simData, "KRNSWMDC_OW", sd.krnswdc_ow );
568  addToSimData( simData, "PCSWMDC_GO", sd.pcswmdc_go );
569  addToSimData( simData, "KRNSWMDC_GO", sd.krnswdc_go );
570 
571  return simData;
572  }
573 
574  // in case the data is already in a SimulationDataContainer no
575  // conversion is needed
576  inline
577  SimulationDataContainer&&
578  convertToSimulationDataContainer( SimulationDataContainer&& sd,
579  const SimulationDataContainer& ,
580  const Opm::PhaseUsage& )
581  {
582  return std::move( sd );
583  }
584 
590  template<class Model>
591  void getRestartData(data::Solution& output,
592  SimulationDataContainer&& sd,
593  const Opm::PhaseUsage& /* phaseUsage */,
594  const Model& /* physicalModel */,
595  const RestartConfig& restartConfig,
596  const int reportStepNum,
597  const bool log)
598  {
599  //Get the value of each of the keys for the restart keywords
600  std::map<std::string, int> rstKeywords = restartConfig.getRestartKeywords(reportStepNum);
601  for (auto& keyValue : rstKeywords) {
602  keyValue.second = restartConfig.getKeyword(keyValue.first, reportStepNum);
603  }
604 
605  const bool aqua_active = sd.hasCellData("1OVERBW");
606  const bool liquid_active = sd.hasCellData("1OVERBO");
607  const bool vapour_active = sd.hasCellData("1OVERBG");
608 
609  assert( aqua_active == (sd.hasCellData("WAT_DEN") &&
610  sd.hasCellData("WAT_VISC") &&
611  sd.hasCellData("WATKR")
612  )
613  );
614  assert( liquid_active == (sd.hasCellData("OIL_DEN") &&
615  sd.hasCellData("OIL_VISC") &&
616  sd.hasCellData("OILKR")
617  )
618  );
619  assert( vapour_active == (sd.hasCellData("GAS_DEN") &&
620  sd.hasCellData("GAS_VISC") &&
621  sd.hasCellData("GASKR")
622  )
623  );
624 
628  if (aqua_active && rstKeywords["BW"] > 0) {
629  rstKeywords["BW"] = 0;
630  output.insert("1OVERBW",
631  Opm::UnitSystem::measure::water_inverse_formation_volume_factor,
632  std::move( sd.getCellData("1OVERBW") ),
633  data::TargetType::RESTART_AUXILIARY);
634  }
635  if (liquid_active && rstKeywords["BO"] > 0) {
636  rstKeywords["BO"] = 0;
637  output.insert("1OVERBO",
638  Opm::UnitSystem::measure::oil_inverse_formation_volume_factor,
639  std::move( sd.getCellData("1OVERBO") ),
640  data::TargetType::RESTART_AUXILIARY);
641  }
642  if (vapour_active && rstKeywords["BG"] > 0) {
643  rstKeywords["BG"] = 0;
644  output.insert("1OVERBG",
645  Opm::UnitSystem::measure::gas_inverse_formation_volume_factor,
646  std::move( sd.getCellData("1OVERBG") ),
647  data::TargetType::RESTART_AUXILIARY);
648  }
649 
653  if (rstKeywords["DEN"] > 0) {
654  rstKeywords["DEN"] = 0;
655  if (aqua_active) {
656  output.insert("WAT_DEN",
657  Opm::UnitSystem::measure::density,
658  std::move( sd.getCellData("WAT_DEN") ),
659  data::TargetType::RESTART_AUXILIARY);
660  }
661  if (liquid_active) {
662  output.insert("OIL_DEN",
663  Opm::UnitSystem::measure::density,
664  std::move( sd.getCellData("OIL_DEN") ),
665  data::TargetType::RESTART_AUXILIARY);
666  }
667  if (vapour_active) {
668  output.insert("GAS_DEN",
669  Opm::UnitSystem::measure::density,
670  std::move( sd.getCellData("GAS_DEN") ),
671  data::TargetType::RESTART_AUXILIARY);
672  }
673  }
674 
678  {
679  const bool has_vwat = (rstKeywords["VISC"] > 0) || (rstKeywords["VWAT"] > 0);
680  const bool has_voil = (rstKeywords["VISC"] > 0) || (rstKeywords["VOIL"] > 0);
681  const bool has_vgas = (rstKeywords["VISC"] > 0) || (rstKeywords["VGAS"] > 0);
682  rstKeywords["VISC"] = 0;
683  if (aqua_active && has_vwat) {
684  output.insert("WAT_VISC",
685  Opm::UnitSystem::measure::viscosity,
686  std::move( sd.getCellData("WAT_VISC") ),
687  data::TargetType::RESTART_AUXILIARY);
688  rstKeywords["VWAT"] = 0;
689  }
690  if (liquid_active && has_voil) {
691  output.insert("OIL_VISC",
692  Opm::UnitSystem::measure::viscosity,
693  std::move( sd.getCellData("OIL_VISC") ),
694  data::TargetType::RESTART_AUXILIARY);
695  rstKeywords["VOIL"] = 0;
696  }
697  if (vapour_active && has_vgas) {
698  output.insert("GAS_VISC",
699  Opm::UnitSystem::measure::viscosity,
700  std::move( sd.getCellData("GAS_VISC") ),
701  data::TargetType::RESTART_AUXILIARY);
702  rstKeywords["VGAS"] = 0;
703  }
704  }
705 
709  if (aqua_active && rstKeywords["KRW"] > 0) {
710  auto& krWater = sd.getCellData("WATKR");
711  if (krWater.size() > 0) {
712  rstKeywords["KRW"] = 0;
713  output.insert("WATKR", // WAT_KR ???
714  Opm::UnitSystem::measure::identity,
715  std::move( krWater ),
716  data::TargetType::RESTART_AUXILIARY);
717  }
718  else {
719  if ( log )
720  {
721  Opm::OpmLog::warning("Empty:WATKR",
722  "Not emitting empty Water Rel-Perm");
723  }
724  }
725  }
726  if (liquid_active && rstKeywords["KRO"] > 0) {
727  auto& krOil = sd.getCellData("OILKR");
728  if (krOil.size() > 0) {
729  rstKeywords["KRO"] = 0;
730  output.insert("OILKR",
731  Opm::UnitSystem::measure::identity,
732  std::move( krOil ),
733  data::TargetType::RESTART_AUXILIARY);
734  }
735  else {
736  if ( log )
737  {
738  Opm::OpmLog::warning("Empty:OILKR",
739  "Not emitting empty Oil Rel-Perm");
740  }
741  }
742  }
743  if (vapour_active && rstKeywords["KRG"] > 0) {
744  auto& krGas = sd.getCellData("GASKR");
745  if (krGas.size() > 0) {
746  rstKeywords["KRG"] = 0;
747  output.insert("GASKR",
748  Opm::UnitSystem::measure::identity,
749  std::move( krGas ),
750  data::TargetType::RESTART_AUXILIARY);
751  }
752  else {
753  if ( log )
754  {
755  Opm::OpmLog::warning("Empty:GASKR",
756  "Not emitting empty Gas Rel-Perm");
757  }
758  }
759  }
760 
764  if (vapour_active && liquid_active && rstKeywords["RSSAT"] > 0) {
765  rstKeywords["RSSAT"] = 0;
766  output.insert("RSSAT",
767  Opm::UnitSystem::measure::gas_oil_ratio,
768  std::move( sd.getCellData("RSSAT") ),
769  data::TargetType::RESTART_AUXILIARY);
770  }
771  if (vapour_active && liquid_active && rstKeywords["RVSAT"] > 0) {
772  rstKeywords["RVSAT"] = 0;
773  output.insert("RVSAT",
774  Opm::UnitSystem::measure::oil_gas_ratio,
775  std::move( sd.getCellData("RVSAT") ),
776  data::TargetType::RESTART_AUXILIARY);
777  }
778 
779 
783  if (vapour_active && liquid_active && rstKeywords["PBPD"] > 0) {
784  rstKeywords["PBPD"] = 0;
785  if (sd.hasCellData("PBUB")) {
786  output.insert("PBUB",
787  Opm::UnitSystem::measure::pressure,
788  std::move( sd.getCellData("PBUB") ),
789  data::TargetType::RESTART_AUXILIARY);
790  }
791  else if (log) {
792  Opm::OpmLog::warning("Bubble point pressure unavailable", "Output of bubble point pressure requested but not available in this simulator. Ignoring.");
793  }
794 
795  if (sd.hasCellData("PDEW")) {
796  output.insert("PDEW",
797  Opm::UnitSystem::measure::pressure,
798  std::move( sd.getCellData("PDEW") ),
799  data::TargetType::RESTART_AUXILIARY);
800  }
801  else if (log) {
802  Opm::OpmLog::warning("Dew point pressure unavailable", "Output of dew point pressure requested but not available in this simulator. Ignoring.");
803  }
804 
805  }
806 
807  if (sd.hasCellData("SOMAX")) {
808  output.insert("SOMAX",
809  Opm::UnitSystem::measure::identity,
810  std::move( sd.getCellData("SOMAX") ),
811  data::TargetType::RESTART_AUXILIARY);
812  }
813 
814  if (sd.hasCellData("PCSWMDC_OW")) {
815  output.insert("PCSWM_OW", //FIXME: Eight-long variable name
816  Opm::UnitSystem::measure::identity,
817  std::move( sd.getCellData("PCSWMDC_OW") ),
818  data::TargetType::RESTART_AUXILIARY);
819  }
820  if (sd.hasCellData("KRNSWMDC_OW")) {
821  output.insert("KRNSW_OW",
822  Opm::UnitSystem::measure::identity,
823  std::move( sd.getCellData("KRNSWMDC_OW") ),
824  data::TargetType::RESTART_AUXILIARY);
825  }
826 
827  if (sd.hasCellData("PCSWMDC_GO")) {
828  output.insert("PCSWM_GO", //FIXME: Eight-long variable name
829  Opm::UnitSystem::measure::identity,
830  std::move( sd.getCellData("PCSWMDC_GO") ),
831  data::TargetType::RESTART_AUXILIARY);
832  }
833  if (sd.hasCellData("KRNSWMDC_GO")) {
834  output.insert("KRNSW_GO",
835  Opm::UnitSystem::measure::identity,
836  std::move( sd.getCellData("KRNSWMDC_GO") ),
837  data::TargetType::RESTART_AUXILIARY);
838  }
839 
840  //Warn for any unhandled keyword
841  if (log) {
842  for (auto& keyValue : rstKeywords) {
843  if (keyValue.second > 0) {
844  std::string logstring = "Keyword '";
845  logstring.append(keyValue.first);
846  logstring.append("' is unhandled for output to file.");
847  Opm::OpmLog::warning("Unhandled output keyword", logstring);
848  }
849  }
850  }
851  }
852 
853 
854 
855 
859  inline bool hasFRBKeyword(const SummaryConfig& summaryConfig, const std::string keyword) {
860  std::string field_kw = "F" + keyword;
861  std::string region_kw = "R" + keyword;
862  std::string block_kw = "B" + keyword;
863  return summaryConfig.hasKeyword(field_kw)
864  || summaryConfig.hasKeyword(region_kw)
865  || summaryConfig.hasKeyword(block_kw);
866  }
867 
868 
872  template<class Model>
873  void getSummaryData(data::Solution& output,
874  const Opm::PhaseUsage& phaseUsage,
875  const Model& physicalModel,
876  const SummaryConfig& summaryConfig) {
877 
878  typedef typename Model::FIPDataType FIPDataType;
879  typedef typename FIPDataType::VectorType VectorType;
880 
881  FIPDataType fd = physicalModel.getFIPData();
882 
883  //Get shorthands for water, oil, gas
884  const int aqua_active = phaseUsage.phase_used[Opm::PhaseUsage::Aqua];
885  const int liquid_active = phaseUsage.phase_used[Opm::PhaseUsage::Liquid];
886  const int vapour_active = phaseUsage.phase_used[Opm::PhaseUsage::Vapour];
887 
891  // Water in place
892  if (aqua_active && hasFRBKeyword(summaryConfig, "WIP")) {
893  output.insert("WIP",
894  Opm::UnitSystem::measure::volume,
895  std::move( fd.fip[ FIPDataType::FIP_AQUA ] ),
896  data::TargetType::SUMMARY );
897  }
898  if (liquid_active) {
899  const VectorType& oipl = fd.fip[FIPDataType::FIP_LIQUID];
900  VectorType oip ( oipl );
901  const size_t size = oip.size();
902 
903  const VectorType& oipg = vapour_active ? fd.fip[FIPDataType::FIP_VAPORIZED_OIL] : VectorType(size, 0.0);
904  if( vapour_active )
905  {
906  // oip = oipl + oipg
907  for( size_t i=0; i<size; ++ i ) {
908  oip[ i ] += oipg[ i ];
909  }
910  }
911 
912  //Oil in place (liquid phase only)
913  if (hasFRBKeyword(summaryConfig, "OIPL")) {
914  output.insert("OIPL",
915  Opm::UnitSystem::measure::volume,
916  std::move( oipl ),
917  data::TargetType::SUMMARY );
918  }
919  //Oil in place (gas phase only)
920  if (hasFRBKeyword(summaryConfig, "OIPG")) {
921  output.insert("OIPG",
922  Opm::UnitSystem::measure::volume,
923  std::move( oipg ),
924  data::TargetType::SUMMARY );
925  }
926  // Oil in place (in liquid and gas phases)
927  if (hasFRBKeyword(summaryConfig, "OIP") || hasFRBKeyword(summaryConfig, "OE")) {
928  output.insert("OIP",
929  Opm::UnitSystem::measure::volume,
930  std::move( oip ),
931  data::TargetType::SUMMARY );
932  }
933  }
934  if (vapour_active) {
935  const VectorType& gipg = fd.fip[ FIPDataType::FIP_VAPOUR];
936  VectorType gip( gipg );
937  const size_t size = gip.size();
938 
939  const VectorType& gipl = liquid_active ? fd.fip[ FIPDataType::FIP_DISSOLVED_GAS ] : VectorType(size,0.0);
940  if( liquid_active )
941  {
942  // gip = gipg + gipl
943  for( size_t i=0; i<size; ++ i ) {
944  gip[ i ] += gipl[ i ];
945  }
946  }
947 
948  // Gas in place (gas phase only)
949  if (hasFRBKeyword(summaryConfig, "GIPG")) {
950  output.insert("GIPG",
951  Opm::UnitSystem::measure::volume,
952  std::move( gipg ),
953  data::TargetType::SUMMARY );
954  }
955 
956  // Gas in place (liquid phase only)
957  if (hasFRBKeyword(summaryConfig, "GIPL")) {
958  output.insert("GIPL",
959  Opm::UnitSystem::measure::volume,
960  std::move( gipl ),
961  data::TargetType::SUMMARY );
962  }
963  // Gas in place (in both liquid and gas phases)
964  if (hasFRBKeyword(summaryConfig, "GIP")) {
965  output.insert("GIP",
966  Opm::UnitSystem::measure::volume,
967  std::move( gip ),
968  data::TargetType::SUMMARY );
969  }
970  }
971  // Cell pore volume in reservoir conditions
972  if (hasFRBKeyword(summaryConfig, "RPV")) {
973  output.insert("RPV",
974  Opm::UnitSystem::measure::volume,
975  std::move( fd.fip[FIPDataType::FIP_PV]),
976  data::TargetType::SUMMARY );
977  }
978  // Pressure averaged value (hydrocarbon pore volume weighted)
979  if (summaryConfig.hasKeyword("FPRH") || summaryConfig.hasKeyword("RPRH")) {
980  output.insert("PRH",
981  Opm::UnitSystem::measure::pressure,
982  std::move(fd.fip[FIPDataType::FIP_WEIGHTED_PRESSURE]),
983  data::TargetType::SUMMARY );
984  }
985  }
986 
987  }
988 
989 
990 
991 
992  template<class Model>
993  inline void
996  const SimulationDataContainer& localState,
997  const WellStateFullyImplicitBlackoil& localWellState,
998  const Model& physicalModel,
999  const bool substep,
1000  const double nextstep,
1001  const SimulatorReport& simulatorReport)
1002  {
1003  data::Solution localCellData{};
1004  const RestartConfig& restartConfig = eclipseState_.getRestartConfig();
1005  const SummaryConfig& summaryConfig = eclipseState_.getSummaryConfig();
1006  const int reportStepNum = timer.reportStepNum();
1007  bool logMessages = output_ && parallelOutput_->isIORank();
1008  std::map<std::string, std::vector<double>> extraRestartData;
1009  std::map<std::string, double> miscSummaryData;
1010 
1011  if( output_ )
1012  {
1013  {
1014  // get all data that need to be included in output from the model
1015  // for flow_legacy and polymer this is a struct holding the data
1016  // while for flow_ebos a SimulationDataContainer is returned
1017  // this is addressed in the above specialized methods
1018  SimulationDataContainer sd =
1019  detail::convertToSimulationDataContainer( physicalModel.getSimulatorData(localState), localState, phaseUsage_ );
1020 
1021  localCellData = simToSolution( sd, restart_double_si_, phaseUsage_); // Get "normal" data (SWAT, PRESSURE, ...);
1022 
1023  detail::getRestartData( localCellData, std::move(sd), phaseUsage_, physicalModel,
1024  restartConfig, reportStepNum, logMessages );
1025  // sd will be invalid after getRestartData has been called
1026  }
1027  detail::getSummaryData( localCellData, phaseUsage_, physicalModel, summaryConfig );
1028  assert(!localCellData.empty());
1029 
1030  // Add suggested next timestep to extra data.
1031  extraRestartData["OPMEXTRA"] = std::vector<double>(1, nextstep);
1032 
1033  // Add TCPU if simulatorReport is not defaulted.
1034  const double totalSolverTime = simulatorReport.solver_time;
1035  if (totalSolverTime != 0.0) {
1036  miscSummaryData["TCPU"] = totalSolverTime;
1037  }
1038  }
1039 
1040  writeTimeStepWithCellProperties(timer, localState, localCellData, localWellState, miscSummaryData, extraRestartData, substep);
1041  }
1042 }
1043 #endif
Extra data to read/write for OPM restarting.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:199
void writeTimeStepWithCellProperties(const SimulatorTimerInterface &timer, const SimulationDataContainer &reservoirState, const data::Solution &cellData, const Opm::WellStateFullyImplicitBlackoil &wellState, const std::map< std::string, double > &miscSummaryData, const std::map< std::string, std::vector< double >> &extraRestartData, bool substep=false)
Write a blackoil reservoir state to disk for later inspection with visualization tools like ResInsigh...
Definition: SimulatorFullyImplicitBlackoilOutput.cpp:244
bool isIORank() const
Whether this process does write to disk.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:284
void writeTimeStepSerial(const SimulatorTimerInterface &timer, const SimulationDataContainer &reservoirState, const Opm::WellStateFullyImplicitBlackoil &wellState, const data::Solution &simProps, const std::map< std::string, double > &miscSummaryData, const std::map< std::string, std::vector< double >> &extraRestartData, bool substep)
Write a blackoil reservoir state to disk for later inspection with visualization tools like ResInsigh...
Definition: SimulatorFullyImplicitBlackoilOutput.cpp:317
void wellsToState(const data::Wells &wells, PhaseUsage phases, WellStateFullyImplicitBlackoil &state)
Copies the following fields from wells into state.
Definition: Compat.cpp:221
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:175
data::Solution simToSolution(const SimulationDataContainer &reservoir, const bool use_si_units, PhaseUsage phases)
Returns Solution with the following fields: PRESSURE, TEMP (unconditionally) SWAT, SGAS, RS, RV, SSOL (if appropriate fields present in input) If use_si_units is true, the fields will have the measure UnitSystem::measure::identity, and therefore not be converted to customary units (depending on family) upon output.
Definition: Compat.cpp:77
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
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:139
bool isIORank(const boost::any &parallel_info)
Return true if this is a serial run, or rank zero on an MPI run.
Definition: NewtonIterationUtilities.cpp:294
void writeTimeStepWithoutCellProperties(const SimulatorTimerInterface &timer, const SimulationDataContainer &reservoirState, const Opm::WellStateFullyImplicitBlackoil &wellState, const std::map< std::string, double > &miscSummaryData, const std::map< std::string, std::vector< double >> &extraRestartData, bool substep=false)
Write a blackoil reservoir state to disk for later inspection with visualization tools like ResInsigh...
Definition: SimulatorFullyImplicitBlackoilOutput.cpp:221
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
void writeInit(const data::Solution &simProps, const NNC &nnc)
bool output() const
return true if output is enabled
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:281
Interface class for SimulatorTimer objects, to be improved.
Definition: SimulatorTimerInterface.hpp:35
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:154
virtual int currentStepNum() const =0
Current step number.
Wrapper class for VTK, Matlab, and ECL output.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:206
const std::string & outputDirectory() const
return output directory
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:278
void solutionToSim(const data::Solution &sol, const std::map< std::string, std::vector< double > > &, PhaseUsage phases, SimulationDataContainer &state)
Copies the following fields from sol into state (all conditionally): PRESSURE, TEMP, SWAT, SGAS, RS, RV, SSOL Also handles extra data such as hysteresis parameters, SOMAX, etc.
Definition: Compat.cpp:134
virtual int reportStepNum() const
Current report step number. This might differ from currentStepNum in case of sub stepping.
Definition: SimulatorTimerInterface.hpp:52