All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
BlackoilSimulator.hpp
1 /*
2  Copyright 2011 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
21 #define OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
22 
23 
24 
25 #include <opm/common/utility/platform_dependent/disable_warnings.h>
26 
27 #include <dune/grid/io/file/vtk/vtkwriter.hh>
28 #include <boost/filesystem/convenience.hpp>
29 #include <boost/lexical_cast.hpp>
30 
31 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
32 
33 
34 #include <opm/parser/eclipse/Units/Units.hpp>
35 #include <opm/core/utility/parameters/ParameterGroup.hpp>
36 #include <opm/porsol/common/BoundaryConditions.hpp>
37 #include <opm/porsol/blackoil/BlackoilInitialization.hpp>
38 #include <opm/porsol/common/SimulatorUtilities.hpp>
39 #include <opm/parser/eclipse/Parser/Parser.hpp>
40 #include <opm/parser/eclipse/Parser/ParseContext.hpp>
41 #include <opm/parser/eclipse/Deck/Deck.hpp>
42 #include <iostream>
43 #include <string>
44 #include <vector>
45 #include <numeric>
46 
47 
48 namespace Opm
49 {
50 
51  template<class GridT, class Rock, class FluidT, class Wells, class FlowSolver, class TransportSolver>
53  {
54  public:
55  void init(const Opm::ParameterGroup& param);
56  void simulate();
57 
58  typedef GridT Grid;
59  typedef FluidT Fluid;
60 
61  typedef typename Fluid::CompVec CompVec;
62  typedef typename Fluid::PhaseVec PhaseVec;
63 
64  struct State
65  {
66  std::vector<PhaseVec> cell_pressure_;
67  std::vector<PhaseVec> face_pressure_;
68  std::vector<double> well_bhp_pressure_;
69  std::vector<double> well_perf_pressure_;
70  std::vector<double> well_perf_flux_;
71  std::vector<CompVec> cell_z_;
72  };
73 
74  private:
75  Grid grid_;
76  Rock rock_;
77  Fluid fluid_;
78  Wells wells_;
79  FlowSolver flow_solver_;
80  TransportSolver transport_solver_;
81 
83  typename Grid::Vector gravity_;
84  std::vector<double> src_;
85 
86  State state_;
87 
88  PhaseVec bdy_pressure_;
89  CompVec bdy_z_;
90 
91  double total_time_;
92  double initial_stepsize_;
93  bool increase_stepsize_;
94  double stepsize_increase_factor_;
95  double minimum_stepsize_;
96  double maximum_stepsize_;
97  std::vector<double> report_times_;
98  bool do_impes_;
99  bool ignore_impes_stability_;
100  std::string output_dir_;
101  int output_interval_;
102 
103  static void output(const Grid& grid,
104  const Fluid& fluid,
105  const State& simstate,
106  const std::vector<double>& face_flux,
107  const int step,
108  const std::string& filebase);
109  };
110 
111 
112 
113 
114  // Method implementations below.
115 
116 
117 template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
118 void
120 init(const Opm::ParameterGroup& param)
121 {
122  using namespace Opm;
123  std::string fileformat = param.getDefault<std::string>("fileformat", "cartesian");
124  if (fileformat == "eclipse") {
125  Opm::ParseContext parseContext;
126  Opm::Parser parser;
127  auto deck = parser.parseFile(param.get<std::string>("filename") , parseContext);
128  Opm::EclipseGrid inputGrid( deck );
129 
130  double z_tolerance = param.getDefault<double>("z_tolerance", 0.0);
131  bool periodic_extension = param.getDefault<bool>("periodic_extension", false);
132  bool turn_normals = param.getDefault<bool>("turn_normals", false);
133  grid_.processEclipseFormat(inputGrid, z_tolerance, periodic_extension, turn_normals);
134  double perm_threshold_md = param.getDefault("perm_threshold_md", 0.0);
135  double perm_threshold = Opm::unit::convert::from(perm_threshold_md, Opm::prefix::milli*Opm::unit::darcy);
136  rock_.init(deck, grid_.globalCell(), perm_threshold);
137  fluid_.init(deck);
138  wells_.init(deck, grid_, rock_);
139  } else if (fileformat == "cartesian") {
140  std::array<int, 3> dims = {{ param.getDefault<int>("nx", 1),
141  param.getDefault<int>("ny", 1),
142  param.getDefault<int>("nz", 1) }};
143  std::array<double, 3> cellsz = {{ param.getDefault<double>("dx", 1.0),
144  param.getDefault<double>("dy", 1.0),
145  param.getDefault<double>("dz", 1.0) }};
146  grid_.createCartesian(dims, cellsz);
147  double default_poro = param.getDefault("default_poro", 1.0);
148  double default_perm_md = param.getDefault("default_perm_md", 100.0);
149  double default_perm = unit::convert::from(default_perm_md, prefix::milli*unit::darcy);
150  OPM_MESSAGE("Warning: For generated cartesian grids, we use uniform rock properties.");
151  rock_.init(grid_.size(0), default_poro, default_perm);
152 
153  Opm::ParseContext parseContext;
154  Opm::Parser parser;
155  auto deck = parser.parseFile(param.get<std::string>("filename") , parseContext);
156  fluid_.init(deck);
157  wells_.init(deck, grid_, rock_);
158  } else {
159  OPM_THROW(std::runtime_error, "Unknown file format string: " << fileformat);
160  }
161  flow_solver_.init(param);
162  transport_solver_.init(param);
163  if (param.has("timestep_file")) {
164  std::ifstream is(param.get<std::string>("timestep_file").c_str());
165  std::istream_iterator<double> beg(is);
166  std::istream_iterator<double> end;
167  report_times_.assign(beg, end);
168  // File contains deltas, we want accumulated times.
169  std::partial_sum(report_times_.begin(), report_times_.end(), report_times_.begin());
170  assert(!report_times_.empty());
171  total_time_ = report_times_.back();
172  initial_stepsize_ = report_times_.front();
173  } else {
174  total_time_ = param.getDefault("total_time", 30*unit::day);
175  initial_stepsize_ = param.getDefault("initial_stepsize", 1.0*unit::day);
176  increase_stepsize_ = param.getDefault("increase_stepsize", false);
177  if (increase_stepsize_) {
178  stepsize_increase_factor_ = param.getDefault("stepsize_increase_factor", 1.5);
179  maximum_stepsize_ = param.getDefault("maximum_stepsize", 1.0*unit::day);
180  } else {
181  stepsize_increase_factor_ = 1.0;
182  maximum_stepsize_ = 1e100;
183  }
184  }
185  minimum_stepsize_ = param.getDefault("minimum_stepsize", 0.0);
186  do_impes_ = param.getDefault("do_impes", false);
187  if (do_impes_) {
188  ignore_impes_stability_ = param.getDefault("ignore_impes_stability", false);
189  }
190  output_dir_ = param.getDefault<std::string>("output_dir", "output");
191  output_interval_ = param.getDefault("output_interval", 1);
192 
193  // Boundary conditions.
194  typedef Opm::FlowBC BC;
195  flow_bc_.resize(7);
196  bool bdy_dirichlet = param.getDefault("bdy_dirichlet", false);
197  if (bdy_dirichlet) {
198  flow_bc_.flowCond(1) = BC(BC::Dirichlet, param.get<double>("bdy_pressure_left"));
199  flow_bc_.flowCond(2) = BC(BC::Dirichlet, param.get<double>("bdy_pressure_right"));
200  } else if (param.getDefault("lateral_dirichlet", false)) {
201  flow_bc_.flowCond(1) = BC(BC::Dirichlet, -17.0); // Use a negative value to instruct flow solver
202  flow_bc_.flowCond(2) = BC(BC::Dirichlet, -17.0); // to use initial face pressures (hydrostatic)
203  flow_bc_.flowCond(3) = BC(BC::Dirichlet, -17.0); // as boundary conditions.
204  flow_bc_.flowCond(4) = BC(BC::Dirichlet, -17.0);
205  }
206 
207  // Gravity.
208  gravity_ = 0.0;
209  if (param.has("gravity")) {
210  std::string g = param.get<std::string>("gravity");
211  if (g == "standard") {
212  gravity_[2] = Opm::unit::gravity;
213  } else {
214  gravity_[2] = boost::lexical_cast<double>(g);
215  }
216  }
217 
218  // Initial state.
219  if (param.getDefault("spe9_init", false)) {
221  initializer.init(param, grid_, fluid_, gravity_, state_);
222  } else {
224  initializer.init(param, grid_, fluid_, gravity_, state_);
225  }
226 
227 
228 
229  // Write initial state to std::cout
230  /*
231  for (int cell = 0; cell < grid_.numCells(); ++cell) {
232  std::cout.precision(2);
233  std::cout << std::fixed << std::showpoint;
234 
235  std::cout << std::setw(5) << cell << std::setw(12) << grid_.cellCentroid(cell)[0]
236  << std::setw(12) << grid_.cellCentroid(cell)[1]
237  << std::setw(12) << grid_.cellCentroid(cell)[2]
238  << std::setw(20) << state_.cell_pressure_[cell][0]
239  << std::setw(15) << state_.cell_z_[cell][0]
240  << std::setw(15) << state_.cell_z_[cell][1]
241  << std::setw(15) << state_.cell_z_[cell][2]
242  << std::endl;
243  if ((cell+1)%nz == 0) {
244  std::cout << "------------------------------------------------------------------------------------------------------------------" << std::endl;
245  }
246 
247  }
248  */
249 
250  bdy_z_ = flow_solver_.inflowMixture();
251  bdy_pressure_ = 300.0*Opm::unit::barsa;
252  // PhaseVec bdy_pressure_(100.0*Opm::unit::barsa); // WELLS
253  // Rescale z values so that pore volume is filled exactly
254  // (to get zero initial volume discrepancy).
255  for (int cell = 0; cell < grid_.numCells(); ++cell) {
256  typename Fluid::FluidState state = fluid_.computeState(state_.cell_pressure_[cell], state_.cell_z_[cell]);
257  double fluid_vol_dens = state.total_phase_volume_density_;
258  state_.cell_z_[cell] *= 1.0/fluid_vol_dens;
259  }
260  int num_faces = grid_.numFaces();
261  state_.face_pressure_.resize(num_faces);
262  for (int face = 0; face < num_faces; ++face) {
263  int bid = grid_.boundaryId(face);
264  if (flow_bc_.flowCond(bid).isDirichlet() && flow_bc_.flowCond(bid).pressure() >= 0.0) {
265  state_.face_pressure_[face] = flow_bc_.flowCond(bid).pressure();
266  } else {
267  int c[2] = { grid_.faceCell(face, 0), grid_.faceCell(face, 1) };
268  state_.face_pressure_[face] = 0.0;
269  int num = 0;
270  for (int j = 0; j < 2; ++j) {
271  if (c[j] >= 0) {
272  state_.face_pressure_[face] += state_.cell_pressure_[c[j]];
273  ++num;
274  }
275  }
276  state_.face_pressure_[face] /= double(num);
277  }
278  }
279 
280  // Flow solver setup.
281  flow_solver_.setup(grid_, rock_, fluid_, wells_, gravity_, flow_bc_, &(state_.face_pressure_));
282 
283  // Transport solver setup.
284  transport_solver_.setup(grid_, rock_, fluid_, wells_, flow_solver_.faceTransmissibilities(), gravity_);
285 
286  // Simple source terms.
287  src_.resize(grid_.numCells(), 0.0);
288 
289  // Set initial well perforation pressures equal to cell pressures,
290  // and perforation fluxes equal to zero.
291  // Set initial well bhp values to the target if bhp well, or to
292  // first perforation pressure if not.
293  state_.well_perf_pressure_.clear();
294  for (int well = 0; well < wells_.numWells(); ++well) {
295  int num_perf = wells_.numPerforations(well);
296  for (int perf = 0; perf < num_perf; ++perf) {
297  int cell = wells_.wellCell(well, perf);
298  state_.well_perf_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
299  }
300  if (wells_.control(well) == Wells::Pressure) {
301  state_.well_bhp_pressure_.push_back(wells_.target(well));
302  } else {
303  int cell = wells_.wellCell(well, 0);
304  state_.well_bhp_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
305  }
306  }
307  state_.well_perf_flux_.clear();
308  state_.well_perf_flux_.resize(state_.well_perf_pressure_.size(), 0.0);
309  wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
310 
311  // Check for unused parameters (potential typos).
312  if (param.anyUnused()) {
313  std::cout << "***** WARNING: Unused parameters: *****\n";
314  param.displayUsage();
315  }
316 
317  // Write parameters used to file, ensuring directory exists.
318  std::string paramfilename = output_dir_ + "/simulator-parameters.param";
319  boost::filesystem::path fpath(paramfilename);
320  if (fpath.has_branch_path()) {
321  create_directories(fpath.branch_path());
322  }
323  param.writeParam(paramfilename);
324 }
325 
326 
327 
328 
329 
330 
331 
332 template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
333 void
335 simulate()
336 {
337  double voldisclimit = flow_solver_.volumeDiscrepancyLimit();
338  double stepsize = initial_stepsize_;
339  double current_time = 0.0;
340  int step = 0;
341  std::vector<double> face_flux;
342  State start_state;
343  std::string output_name = output_dir_ + "/" + "blackoil-output";
344  while (current_time < total_time_) {
345  start_state = state_;
346 
347  // Do not run past total_time_.
348  if (current_time + stepsize > total_time_) {
349  stepsize = total_time_ - current_time;
350  }
351  std::cout << "\n\n================ Simulation step number " << step
352  << " ==============="
353  << "\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
354  << "\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
355  << "\n Total time (days) " << Opm::unit::convert::to(total_time_, Opm::unit::day)
356  << "\n" << std::endl;
357 
358  // Solve flow system.
359  enum FlowSolver::ReturnCode result
360  = flow_solver_.solve(state_.cell_pressure_, state_.face_pressure_, state_.cell_z_, face_flux,
361  state_.well_bhp_pressure_,
362  state_.well_perf_pressure_, state_.well_perf_flux_, src_, stepsize);
363 
364  // Check if the flow solver succeeded.
365  if (result == FlowSolver::VolumeDiscrepancyTooLarge) {
366  OPM_THROW(std::runtime_error, "Flow solver refused to run due to too large volume discrepancy.");
367  } else if (result == FlowSolver::FailedToConverge) {
368  std::cout << "********* Nonlinear convergence failure: Shortening (pressure) stepsize, redoing step number " << step <<" **********" << std::endl;
369  stepsize *= 0.5;
370  state_ = start_state;
371  wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
372  continue;
373  }
374  assert(result == FlowSolver::SolveOk);
375 
376  // Update wells with new perforation pressures and fluxes.
377  wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
378 
379  // Transport and check volume discrepancy.
380  bool voldisc_ok = true;
381  if (!do_impes_) {
382  double actual_computed_time
383  = transport_solver_.transport(bdy_pressure_, bdy_z_,
384  face_flux, state_.cell_pressure_, state_.face_pressure_,
385  stepsize, voldisclimit, state_.cell_z_);
386  voldisc_ok = (actual_computed_time == stepsize);
387  if (voldisc_ok) {
388  // Just for output.
389  flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
390  state_.well_perf_pressure_, state_.cell_z_, stepsize);
391  }
392  } else {
393  // First check IMPES stepsize.
394  double max_dt = ignore_impes_stability_ ? 1e100 : flow_solver_.stableStepIMPES();
395  if (ignore_impes_stability_) {
396  std::cout << "Timestep was " << stepsize << " and max stepsize was not computed." << std::endl;
397  } else {
398  std::cout << "Timestep was " << stepsize << " and max stepsize was " << max_dt << std::endl;
399  }
400  if (stepsize < max_dt || stepsize <= minimum_stepsize_) {
401  flow_solver_.doStepIMPES(state_.cell_z_, stepsize);
402  voldisc_ok = flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
403  state_.well_perf_pressure_, state_.cell_z_, stepsize);
404  } else {
405  // Restarting step.
406  stepsize = max_dt/1.5;
407  std::cout << "Restarting pressure step with new timestep " << stepsize << std::endl;
408  state_ = start_state;
409  wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
410  continue;
411  }
412  }
413 
414  // If discrepancy too large, redo entire pressure step.
415  if (!voldisc_ok) {
416  std::cout << "********* Too large volume discrepancy: Shortening (pressure) stepsize, redoing step number " << step <<" **********" << std::endl;
417  stepsize *= 0.5;
418  state_ = start_state;
419  wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
420  continue;
421  }
422 
423  // Adjust time.
424  current_time += stepsize;
425  if (voldisc_ok && increase_stepsize_ && stepsize < maximum_stepsize_) {
426  stepsize *= stepsize_increase_factor_;
427  stepsize = std::min(maximum_stepsize_, stepsize);
428  }
429  // If using given timesteps, set stepsize to match.
430  if (!report_times_.empty()) {
431  if (current_time >= report_times_[step]) {
432  bool output_now = ((step + 1) % output_interval_ == 0);
433  if (output_now) {
434  output(grid_, fluid_, state_, face_flux, step, output_name);
435  }
436  ++step;
437  if (step == int(report_times_.size())) {
438  break;
439  }
440  }
441  stepsize = report_times_[step] - current_time;
442  } else {
443  bool output_now = ((step + 1) % output_interval_ == 0);
444  if (output_now) {
445  output(grid_, fluid_, state_, face_flux, step, output_name);
446  }
447  ++step;
448  }
449  }
450  if (step % output_interval_ != 0) {
451  // Output was not written at last step, write final output.
452  output(grid_, fluid_, state_, face_flux, step - 1, output_name);
453  }
454 }
455 
456 
457 
458 
459 
460 
461 template<class Grid, class Rock, class Fluid, class Wells, class FlowSolver, class TransportSolver>
462 void
464 output(const Grid& grid,
465  const Fluid& fluid,
466  const State& simstate,
467  const std::vector<double>& face_flux,
468  const int step,
469  const std::string& filebase)
470 {
471  // Compute saturations, total fluid volume density and mass fractions.
472  int num_cells = grid.numCells();
473  std::vector<typename Fluid::PhaseVec> sat(num_cells);
474  std::vector<typename Fluid::PhaseVec> mass_frac(num_cells);
475  std::vector<double> totflvol_dens(num_cells);
476  for (int cell = 0; cell < num_cells; ++cell) {
477  typename Fluid::FluidState fstate = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
478  sat[cell] = fstate.saturation_;
479  totflvol_dens[cell] = fstate.total_phase_volume_density_;
480  double totMass_dens = simstate.cell_z_[cell]*fluid.surfaceDensities();
481  mass_frac[cell][Fluid::Water] = simstate.cell_z_[cell][Fluid::Water]*fluid.surfaceDensities()[Fluid::Water]/totMass_dens;
482  mass_frac[cell][Fluid::Oil] = simstate.cell_z_[cell][Fluid::Oil]*fluid.surfaceDensities()[Fluid::Oil]/totMass_dens;
483  mass_frac[cell][Fluid::Gas] = simstate.cell_z_[cell][Fluid::Gas]*fluid.surfaceDensities()[Fluid::Gas]/totMass_dens;
484  }
485 
486  // Ensure directory exists.
487  boost::filesystem::path fpath(filebase);
488  if (fpath.has_branch_path()) {
489  create_directories(fpath.branch_path());
490  }
491 
492  // Output to VTK.
493  std::vector<typename Grid::Vector> cell_velocity;
494  Opm::estimateCellVelocitySimpleInterface(cell_velocity, grid, face_flux);
495  // Dune's vtk writer wants multi-component data to be flattened.
496  std::vector<double> cell_pressure_flat(&*simstate.cell_pressure_.front().begin(),
497  &*simstate.cell_pressure_.back().end());
498  std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
499  &*cell_velocity.back().end());
500  std::vector<double> z_flat(&*simstate.cell_z_.front().begin(),
501  &*simstate.cell_z_.back().end());
502  std::vector<double> sat_flat(&*sat.front().begin(),
503  &*sat.back().end());
504  std::vector<double> mass_frac_flat(&*mass_frac.front().begin(),
505  &*mass_frac.back().end());
506 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
507  Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafGridView());
508 #else
509  Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
510 #endif
511  vtkwriter.addCellData(cell_pressure_flat, "pressure", Fluid::numPhases);
512  vtkwriter.addCellData(cell_velocity_flat, "velocity", Grid::dimension);
513  vtkwriter.addCellData(z_flat, "z", Fluid::numComponents);
514  vtkwriter.addCellData(sat_flat, "sat", Fluid::numPhases);
515  vtkwriter.addCellData(mass_frac_flat, "massFrac", Fluid::numComponents);
516  vtkwriter.addCellData(totflvol_dens, "total fl. vol.");
517  vtkwriter.write(filebase + '-' + boost::lexical_cast<std::string>(step),
518  Dune::VTK::ascii);
519 
520  // Dump data for Matlab.
521  std::vector<double> zv[Fluid::numComponents];
522  for (int comp = 0; comp < Fluid::numComponents; ++comp) {
523  zv[comp].resize(grid.numCells());
524  for (int cell = 0; cell < grid.numCells(); ++cell) {
525  zv[comp][cell] = simstate.cell_z_[cell][comp];
526  }
527  }
528  std::vector<double> sv[Fluid::numPhases];
529  for (int phase = 0; phase < Fluid::numPhases; ++phase) {
530  sv[phase].resize(grid.numCells());
531  for (int cell = 0; cell < grid.numCells(); ++cell) {
532  sv[phase][cell] = sat[cell][phase];
533  }
534  }
535  std::string matlabdumpname(filebase + "-");
536  matlabdumpname += boost::lexical_cast<std::string>(step);
537  matlabdumpname += ".dat";
538  std::ofstream dump(matlabdumpname.c_str());
539  dump.precision(15);
540  std::vector<double> liq_press(num_cells);
541  for (int cell = 0; cell < num_cells; ++cell) {
542  liq_press[cell] = simstate.cell_pressure_[cell][Fluid::Liquid];
543  }
544  // Liquid phase pressure.
545  std::copy(liq_press.begin(), liq_press.end(),
546  std::ostream_iterator<double>(dump, " "));
547  dump << '\n';
548  // z (3 components)
549  for (int comp = 0; comp < Fluid::numComponents; ++comp) {
550  std::copy(zv[comp].begin(), zv[comp].end(),
551  std::ostream_iterator<double>(dump, " "));
552  dump << '\n';
553  }
554  // s (3 components)
555  for (int phase = 0; phase < Fluid::numPhases; ++phase) {
556  std::copy(sv[phase].begin(), sv[phase].end(),
557  std::ostream_iterator<double>(dump, " "));
558  dump << '\n';
559  }
560  // Total fluid volume
561  std::copy(totflvol_dens.begin(), totflvol_dens.end(),
562  std::ostream_iterator<double>(dump, " "));
563  dump << '\n';
564  // Well report ...
565  const double seconds_pr_day = 3600.*24.;
566  for (unsigned int perf=0; perf<Wells::WellReport::report()->perfPressure.size(); ++perf) {
567  dump << std::setw(8) << Wells::WellReport::report()->cellId[perf] << " "
568  << std::setw(22) << Wells::WellReport::report()->perfPressure[perf] << " "
569  << std::setw(22) << Wells::WellReport::report()->cellPressure[perf] << " "
570  << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Water] << " "
571  << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Oil] << " "
572  << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Gas] << '\n';
573  }
574  dump << '\n';
575 }
576 
577 
578 
579 
580 
581 } // namespace Opm
582 
583 
584 
585 
586 
587 #endif // OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:121
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:82
Initialize basic cases.
Definition: BlackoilInitialization.hpp:52
Initialize SPE9 type case.
Definition: BlackoilInitialization.hpp:193
Definition: BlackoilSimulator.hpp:52
Definition: BlackoilSimulator.hpp:64
A class designed to encapsulate a set of rate- or pressure-controlled wells.
Definition: Wells.hpp:37
void estimateCellVelocitySimpleInterface(std::vector< typename GridInterface::Vector > &cell_velocity, const GridInterface &grid, const std::vector< double > &face_flux)
Estimates a scalar cell velocity from face fluxes.
Definition: SimulatorUtilities.hpp:95