20#ifndef OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
21#define OPM_BLACKOILSIMULATOR_HEADER_INCLUDED
26#include <dune/grid/io/file/vtk/vtkwriter.hh>
27#include <opm/core/utility/Units.hpp>
28#include <opm/core/utility/parameters/ParameterGroup.hpp>
32#include <opm/parser/eclipse/Parser/Parser.hpp>
33#include <opm/parser/eclipse/Parser/ParseContext.hpp>
34#include <opm/parser/eclipse/Deck/Deck.hpp>
35#include <boost/filesystem/convenience.hpp>
36#include <boost/lexical_cast.hpp>
46 template<
class Gr
idT,
class Rock,
class Flu
idT,
class Wells,
class FlowSolver,
class TransportSolver>
50 void init(
const Opm::parameter::ParameterGroup& param);
74 FlowSolver flow_solver_;
75 TransportSolver transport_solver_;
78 typename Grid::Vector gravity_;
79 std::vector<double> src_;
87 double initial_stepsize_;
88 bool increase_stepsize_;
89 double stepsize_increase_factor_;
90 double minimum_stepsize_;
91 double maximum_stepsize_;
92 std::vector<double> report_times_;
94 bool ignore_impes_stability_;
95 std::string output_dir_;
98 static void output(
const Grid& grid,
100 const State& simstate,
101 const std::vector<double>& face_flux,
103 const std::string& filebase);
112template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
115init(
const Opm::parameter::ParameterGroup& param)
118 std::string fileformat = param.getDefault<std::string>(
"fileformat",
"cartesian");
119 if (fileformat ==
"eclipse") {
120 Opm::ParseContext parseContext;
121 Opm::ParserPtr parser(
new Opm::Parser());
122 Opm::DeckConstPtr deck = parser->parseFile(param.get<std::string>(
"filename") , parseContext);
123 double z_tolerance = param.getDefault<
double>(
"z_tolerance", 0.0);
124 bool periodic_extension = param.getDefault<
bool>(
"periodic_extension",
false);
125 bool turn_normals = param.getDefault<
bool>(
"turn_normals",
false);
126 grid_.processEclipseFormat(deck, z_tolerance, periodic_extension, turn_normals);
127 double perm_threshold_md = param.getDefault(
"perm_threshold_md", 0.0);
128 double perm_threshold = Opm::unit::convert::from(perm_threshold_md, Opm::prefix::milli*Opm::unit::darcy);
129 rock_.init(deck, grid_.globalCell(), perm_threshold);
131 wells_.init(deck, grid_, rock_);
132 }
else if (fileformat ==
"cartesian") {
133 std::array<int, 3> dims = {{ param.getDefault<
int>(
"nx", 1),
134 param.getDefault<
int>(
"ny", 1),
135 param.getDefault<
int>(
"nz", 1) }};
136 std::array<double, 3> cellsz = {{ param.getDefault<
double>(
"dx", 1.0),
137 param.getDefault<
double>(
"dy", 1.0),
138 param.getDefault<
double>(
"dz", 1.0) }};
139 grid_.createCartesian(dims, cellsz);
140 double default_poro = param.getDefault(
"default_poro", 1.0);
141 double default_perm_md = param.getDefault(
"default_perm_md", 100.0);
142 double default_perm = unit::convert::from(default_perm_md, prefix::milli*unit::darcy);
143 OPM_MESSAGE(
"Warning: For generated cartesian grids, we use uniform rock properties.");
144 rock_.init(grid_.size(0), default_poro, default_perm);
146 Opm::ParseContext parseContext;
147 Opm::ParserPtr parser(
new Opm::Parser());
148 Opm::DeckConstPtr deck = parser->parseFile(param.get<std::string>(
"filename") , parseContext);
150 wells_.init(deck, grid_, rock_);
152 OPM_THROW(std::runtime_error,
"Unknown file format string: " << fileformat);
154 flow_solver_.init(param);
155 transport_solver_.init(param);
156 if (param.has(
"timestep_file")) {
157 std::ifstream is(param.get<std::string>(
"timestep_file").c_str());
158 std::istream_iterator<double> beg(is);
159 std::istream_iterator<double> end;
160 report_times_.assign(beg, end);
162 std::partial_sum(report_times_.begin(), report_times_.end(), report_times_.begin());
163 assert(!report_times_.empty());
164 total_time_ = report_times_.back();
165 initial_stepsize_ = report_times_.front();
167 total_time_ = param.getDefault(
"total_time", 30*unit::day);
168 initial_stepsize_ = param.getDefault(
"initial_stepsize", 1.0*unit::day);
169 increase_stepsize_ = param.getDefault(
"increase_stepsize",
false);
170 if (increase_stepsize_) {
171 stepsize_increase_factor_ = param.getDefault(
"stepsize_increase_factor", 1.5);
172 maximum_stepsize_ = param.getDefault(
"maximum_stepsize", 1.0*unit::day);
174 stepsize_increase_factor_ = 1.0;
175 maximum_stepsize_ = 1e100;
178 minimum_stepsize_ = param.getDefault(
"minimum_stepsize", 0.0);
179 do_impes_ = param.getDefault(
"do_impes",
false);
181 ignore_impes_stability_ = param.getDefault(
"ignore_impes_stability",
false);
183 output_dir_ = param.getDefault<std::string>(
"output_dir",
"output");
184 output_interval_ = param.getDefault(
"output_interval", 1);
189 bool bdy_dirichlet = param.getDefault(
"bdy_dirichlet",
false);
191 flow_bc_.flowCond(1) = BC(BC::Dirichlet, param.get<
double>(
"bdy_pressure_left"));
192 flow_bc_.flowCond(2) = BC(BC::Dirichlet, param.get<
double>(
"bdy_pressure_right"));
193 }
else if (param.getDefault(
"lateral_dirichlet",
false)) {
194 flow_bc_.flowCond(1) = BC(BC::Dirichlet, -17.0);
195 flow_bc_.flowCond(2) = BC(BC::Dirichlet, -17.0);
196 flow_bc_.flowCond(3) = BC(BC::Dirichlet, -17.0);
197 flow_bc_.flowCond(4) = BC(BC::Dirichlet, -17.0);
202 if (param.has(
"gravity")) {
203 std::string g = param.get<std::string>(
"gravity");
204 if (g ==
"standard") {
205 gravity_[2] = Opm::unit::gravity;
207 gravity_[2] = boost::lexical_cast<double>(g);
212 if (param.getDefault(
"spe9_init",
false)) {
214 initializer.
init(param, grid_, fluid_, gravity_, state_);
217 initializer.
init(param, grid_, fluid_, gravity_, state_);
243 bdy_z_ = flow_solver_.inflowMixture();
244 bdy_pressure_ = 300.0*Opm::unit::barsa;
248 for (
int cell = 0; cell < grid_.numCells(); ++cell) {
249 typename Fluid::FluidState state = fluid_.computeState(state_.cell_pressure_[cell], state_.cell_z_[cell]);
250 double fluid_vol_dens = state.total_phase_volume_density_;
251 state_.cell_z_[cell] *= 1.0/fluid_vol_dens;
253 int num_faces = grid_.numFaces();
254 state_.face_pressure_.resize(num_faces);
255 for (
int face = 0; face < num_faces; ++face) {
256 int bid = grid_.boundaryId(face);
257 if (flow_bc_.flowCond(bid).isDirichlet() && flow_bc_.flowCond(bid).pressure() >= 0.0) {
258 state_.face_pressure_[face] = flow_bc_.flowCond(bid).pressure();
260 int c[2] = { grid_.faceCell(face, 0), grid_.faceCell(face, 1) };
261 state_.face_pressure_[face] = 0.0;
263 for (
int j = 0; j < 2; ++j) {
265 state_.face_pressure_[face] += state_.cell_pressure_[c[j]];
269 state_.face_pressure_[face] /= double(num);
274 flow_solver_.setup(grid_, rock_, fluid_, wells_, gravity_, flow_bc_, &(state_.face_pressure_));
277 transport_solver_.setup(grid_, rock_, fluid_, wells_, flow_solver_.faceTransmissibilities(), gravity_);
280 src_.resize(grid_.numCells(), 0.0);
286 state_.well_perf_pressure_.clear();
287 for (
int well = 0; well < wells_.numWells(); ++well) {
288 int num_perf = wells_.numPerforations(well);
289 for (
int perf = 0; perf < num_perf; ++perf) {
290 int cell = wells_.wellCell(well, perf);
291 state_.well_perf_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
294 state_.well_bhp_pressure_.push_back(wells_.target(well));
296 int cell = wells_.wellCell(well, 0);
297 state_.well_bhp_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
300 state_.well_perf_flux_.clear();
301 state_.well_perf_flux_.resize(state_.well_perf_pressure_.size(), 0.0);
302 wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
305 if (param.anyUnused()) {
306 std::cout <<
"***** WARNING: Unused parameters: *****\n";
307 param.displayUsage();
311 std::string paramfilename = output_dir_ +
"/simulator-parameters.param";
312 boost::filesystem::path fpath(paramfilename);
313 if (fpath.has_branch_path()) {
314 create_directories(fpath.branch_path());
316 param.writeParam(paramfilename);
325template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
330 double voldisclimit = flow_solver_.volumeDiscrepancyLimit();
331 double stepsize = initial_stepsize_;
332 double current_time = 0.0;
334 std::vector<double> face_flux;
336 std::string output_name = output_dir_ +
"/" +
"blackoil-output";
337 while (current_time < total_time_) {
338 start_state = state_;
341 if (current_time + stepsize > total_time_) {
342 stepsize = total_time_ - current_time;
344 std::cout <<
"\n\n================ Simulation step number " << step
345 <<
" ==============="
346 <<
"\n Current time (days) " << Opm::unit::convert::to(current_time, Opm::unit::day)
347 <<
"\n Current stepsize (days) " << Opm::unit::convert::to(stepsize, Opm::unit::day)
348 <<
"\n Total time (days) " << Opm::unit::convert::to(total_time_, Opm::unit::day)
349 <<
"\n" << std::endl;
352 enum FlowSolver::ReturnCode result
353 = flow_solver_.solve(state_.cell_pressure_, state_.face_pressure_, state_.cell_z_, face_flux,
354 state_.well_bhp_pressure_,
355 state_.well_perf_pressure_, state_.well_perf_flux_, src_, stepsize);
358 if (result == FlowSolver::VolumeDiscrepancyTooLarge) {
359 OPM_THROW(std::runtime_error,
"Flow solver refused to run due to too large volume discrepancy.");
360 }
else if (result == FlowSolver::FailedToConverge) {
361 std::cout <<
"********* Nonlinear convergence failure: Shortening (pressure) stepsize, redoing step number " << step <<
" **********" << std::endl;
363 state_ = start_state;
367 assert(result == FlowSolver::SolveOk);
370 wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
373 bool voldisc_ok =
true;
375 double actual_computed_time
376 = transport_solver_.transport(bdy_pressure_, bdy_z_,
377 face_flux, state_.cell_pressure_, state_.face_pressure_,
378 stepsize, voldisclimit, state_.cell_z_);
379 voldisc_ok = (actual_computed_time == stepsize);
382 flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
383 state_.well_perf_pressure_, state_.cell_z_, stepsize);
387 double max_dt = ignore_impes_stability_ ? 1e100 : flow_solver_.stableStepIMPES();
388 if (ignore_impes_stability_) {
389 std::cout <<
"Timestep was " << stepsize <<
" and max stepsize was not computed." << std::endl;
391 std::cout <<
"Timestep was " << stepsize <<
" and max stepsize was " << max_dt << std::endl;
393 if (stepsize < max_dt || stepsize <= minimum_stepsize_) {
394 flow_solver_.doStepIMPES(state_.cell_z_, stepsize);
395 voldisc_ok = flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
396 state_.well_perf_pressure_, state_.cell_z_, stepsize);
399 stepsize = max_dt/1.5;
400 std::cout <<
"Restarting pressure step with new timestep " << stepsize << std::endl;
401 state_ = start_state;
409 std::cout <<
"********* Too large volume discrepancy: Shortening (pressure) stepsize, redoing step number " << step <<
" **********" << std::endl;
411 state_ = start_state;
417 current_time += stepsize;
418 if (voldisc_ok && increase_stepsize_ && stepsize < maximum_stepsize_) {
419 stepsize *= stepsize_increase_factor_;
420 stepsize = std::min(maximum_stepsize_, stepsize);
423 if (!report_times_.empty()) {
424 if (current_time >= report_times_[step]) {
425 bool output_now = ((step + 1) % output_interval_ == 0);
427 output(grid_, fluid_, state_, face_flux, step, output_name);
430 if (step ==
int(report_times_.size())) {
434 stepsize = report_times_[step] - current_time;
436 bool output_now = ((step + 1) % output_interval_ == 0);
438 output(grid_, fluid_, state_, face_flux, step, output_name);
443 if (step % output_interval_ != 0) {
445 output(grid_, fluid_, state_, face_flux, step - 1, output_name);
454template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
459 const State& simstate,
460 const std::vector<double>& face_flux,
462 const std::string& filebase)
465 int num_cells = grid.numCells();
466 std::vector<typename Fluid::PhaseVec> sat(num_cells);
467 std::vector<typename Fluid::PhaseVec> mass_frac(num_cells);
468 std::vector<double> totflvol_dens(num_cells);
469 for (
int cell = 0; cell < num_cells; ++cell) {
470 typename Fluid::FluidState fstate = fluid.computeState(simstate.cell_pressure_[cell], simstate.cell_z_[cell]);
471 sat[cell] = fstate.saturation_;
472 totflvol_dens[cell] = fstate.total_phase_volume_density_;
473 double totMass_dens = simstate.cell_z_[cell]*fluid.surfaceDensities();
474 mass_frac[cell][Fluid::Water] = simstate.cell_z_[cell][Fluid::Water]*fluid.surfaceDensities()[Fluid::Water]/totMass_dens;
475 mass_frac[cell][Fluid::Oil] = simstate.cell_z_[cell][Fluid::Oil]*fluid.surfaceDensities()[Fluid::Oil]/totMass_dens;
476 mass_frac[cell][Fluid::Gas] = simstate.cell_z_[cell][Fluid::Gas]*fluid.surfaceDensities()[Fluid::Gas]/totMass_dens;
480 boost::filesystem::path fpath(filebase);
481 if (fpath.has_branch_path()) {
482 create_directories(fpath.branch_path());
486 std::vector<typename Grid::Vector> cell_velocity;
489 std::vector<double> cell_pressure_flat(&*simstate.cell_pressure_.front().begin(),
490 &*simstate.cell_pressure_.back().end());
491 std::vector<double> cell_velocity_flat(&*cell_velocity.front().begin(),
492 &*cell_velocity.back().end());
493 std::vector<double> z_flat(&*simstate.cell_z_.front().begin(),
494 &*simstate.cell_z_.back().end());
495 std::vector<double> sat_flat(&*sat.front().begin(),
497 std::vector<double> mass_frac_flat(&*mass_frac.front().begin(),
498 &*mass_frac.back().end());
499#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
500 Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafGridView());
502 Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
504 vtkwriter.addCellData(cell_pressure_flat,
"pressure", Fluid::numPhases);
505 vtkwriter.addCellData(cell_velocity_flat,
"velocity", Grid::dimension);
506 vtkwriter.addCellData(z_flat,
"z", Fluid::numComponents);
507 vtkwriter.addCellData(sat_flat,
"sat", Fluid::numPhases);
508 vtkwriter.addCellData(mass_frac_flat,
"massFrac", Fluid::numComponents);
509 vtkwriter.addCellData(totflvol_dens,
"total fl. vol.");
510 vtkwriter.write(filebase +
'-' + boost::lexical_cast<std::string>(step),
514 std::vector<double> zv[Fluid::numComponents];
515 for (
int comp = 0; comp < Fluid::numComponents; ++comp) {
516 zv[comp].resize(grid.numCells());
517 for (
int cell = 0; cell < grid.numCells(); ++cell) {
518 zv[comp][cell] = simstate.cell_z_[cell][comp];
521 std::vector<double> sv[Fluid::numPhases];
522 for (
int phase = 0; phase < Fluid::numPhases; ++phase) {
523 sv[phase].resize(grid.numCells());
524 for (
int cell = 0; cell < grid.numCells(); ++cell) {
525 sv[phase][cell] = sat[cell][phase];
528 std::string matlabdumpname(filebase +
"-");
529 matlabdumpname += boost::lexical_cast<std::string>(step);
530 matlabdumpname +=
".dat";
531 std::ofstream dump(matlabdumpname.c_str());
533 std::vector<double> liq_press(num_cells);
534 for (
int cell = 0; cell < num_cells; ++cell) {
535 liq_press[cell] = simstate.cell_pressure_[cell][Fluid::Liquid];
538 std::copy(liq_press.begin(), liq_press.end(),
539 std::ostream_iterator<double>(dump,
" "));
542 for (
int comp = 0; comp < Fluid::numComponents; ++comp) {
543 std::copy(zv[comp].begin(), zv[comp].end(),
544 std::ostream_iterator<double>(dump,
" "));
548 for (
int phase = 0; phase < Fluid::numPhases; ++phase) {
549 std::copy(sv[phase].begin(), sv[phase].end(),
550 std::ostream_iterator<double>(dump,
" "));
554 std::copy(totflvol_dens.begin(), totflvol_dens.end(),
555 std::ostream_iterator<double>(dump,
" "));
558 const double seconds_pr_day = 3600.*24.;
559 for (
unsigned int perf=0; perf<Wells::WellReport::report()->perfPressure.size(); ++perf) {
560 dump << std::setw(8) << Wells::WellReport::report()->cellId[perf] <<
" "
561 << std::setw(22) << Wells::WellReport::report()->perfPressure[perf] <<
" "
562 << std::setw(22) << Wells::WellReport::report()->cellPressure[perf] <<
" "
563 << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Water] <<
" "
564 << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Oil] <<
" "
565 << std::setw(22) << seconds_pr_day*Wells::WellReport::report()->massRate[perf][Fluid::Gas] <<
'\n';
Initialize basic cases.
Definition: BlackoilInitialization.hpp:53
virtual void init(const Opm::parameter::ParameterGroup ¶m, const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, State &simstate)
Definition: BlackoilInitialization.hpp:59
Definition: BlackoilSimulator.hpp:48
Fluid::CompVec CompVec
Definition: BlackoilSimulator.hpp:56
GridT Grid
Definition: BlackoilSimulator.hpp:53
void simulate()
Definition: BlackoilSimulator.hpp:328
FluidT Fluid
Definition: BlackoilSimulator.hpp:54
Fluid::PhaseVec PhaseVec
Definition: BlackoilSimulator.hpp:57
void init(const Opm::parameter::ParameterGroup ¶m)
Definition: BlackoilSimulator.hpp:115
A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:122
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:80
Initialize SPE9 type case.
Definition: BlackoilInitialization.hpp:194
virtual void init(const Opm::parameter::ParameterGroup ¶m, const Grid &grid, const Fluid &fluid, typename Grid::Vector gravity, State &simstate)
Definition: BlackoilInitialization.hpp:200
@ Pressure
Definition: Wells.hpp:46
Definition: BlackoilFluid.hpp:32
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:90
Definition: BlackoilSimulator.hpp:60
std::vector< double > well_bhp_pressure_
Definition: BlackoilSimulator.hpp:63
std::vector< CompVec > cell_z_
Definition: BlackoilSimulator.hpp:66
std::vector< PhaseVec > cell_pressure_
Definition: BlackoilSimulator.hpp:61
std::vector< PhaseVec > face_pressure_
Definition: BlackoilSimulator.hpp:62
std::vector< double > well_perf_pressure_
Definition: BlackoilSimulator.hpp:64
std::vector< double > well_perf_flux_
Definition: BlackoilSimulator.hpp:65