20 #ifndef OPM_BLACKOILSIMULATOR_HEADER_INCLUDED 21 #define OPM_BLACKOILSIMULATOR_HEADER_INCLUDED 25 #include <opm/common/utility/platform_dependent/disable_warnings.h> 27 #include <dune/grid/io/file/vtk/vtkwriter.hh> 28 #include <boost/filesystem/convenience.hpp> 29 #include <boost/lexical_cast.hpp> 31 #include <opm/common/utility/platform_dependent/reenable_warnings.h> 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> 51 template<
class Gr
idT,
class Rock,
class Flu
idT,
class Wells,
class FlowSolver,
class TransportSolver>
55 void init(
const Opm::ParameterGroup& param);
61 typedef typename Fluid::CompVec CompVec;
62 typedef typename Fluid::PhaseVec PhaseVec;
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_;
79 FlowSolver flow_solver_;
80 TransportSolver transport_solver_;
83 typename Grid::Vector gravity_;
84 std::vector<double> src_;
88 PhaseVec bdy_pressure_;
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_;
99 bool ignore_impes_stability_;
100 std::string output_dir_;
101 int output_interval_;
103 static void output(
const Grid& grid,
105 const State& simstate,
106 const std::vector<double>& face_flux,
108 const std::string& filebase);
117 template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
120 init(
const Opm::ParameterGroup& param)
123 std::string fileformat = param.getDefault<std::string>(
"fileformat",
"cartesian");
124 if (fileformat ==
"eclipse") {
125 Opm::ParseContext parseContext;
127 auto deck = parser.parseFile(param.get<std::string>(
"filename") , parseContext);
128 Opm::EclipseGrid inputGrid( deck );
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);
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);
153 Opm::ParseContext parseContext;
155 auto deck = parser.parseFile(param.get<std::string>(
"filename") , parseContext);
157 wells_.init(deck, grid_, rock_);
159 OPM_THROW(std::runtime_error,
"Unknown file format string: " << fileformat);
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);
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();
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);
181 stepsize_increase_factor_ = 1.0;
182 maximum_stepsize_ = 1e100;
185 minimum_stepsize_ = param.getDefault(
"minimum_stepsize", 0.0);
186 do_impes_ = param.getDefault(
"do_impes",
false);
188 ignore_impes_stability_ = param.getDefault(
"ignore_impes_stability",
false);
190 output_dir_ = param.getDefault<std::string>(
"output_dir",
"output");
191 output_interval_ = param.getDefault(
"output_interval", 1);
196 bool bdy_dirichlet = param.getDefault(
"bdy_dirichlet",
false);
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);
202 flow_bc_.flowCond(2) = BC(BC::Dirichlet, -17.0);
203 flow_bc_.flowCond(3) = BC(BC::Dirichlet, -17.0);
204 flow_bc_.flowCond(4) = BC(BC::Dirichlet, -17.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;
214 gravity_[2] = boost::lexical_cast<
double>(g);
219 if (param.getDefault(
"spe9_init",
false)) {
221 initializer.init(param, grid_, fluid_, gravity_, state_);
224 initializer.init(param, grid_, fluid_, gravity_, state_);
250 bdy_z_ = flow_solver_.inflowMixture();
251 bdy_pressure_ = 300.0*Opm::unit::barsa;
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;
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();
267 int c[2] = { grid_.faceCell(face, 0), grid_.faceCell(face, 1) };
268 state_.face_pressure_[face] = 0.0;
270 for (
int j = 0; j < 2; ++j) {
272 state_.face_pressure_[face] += state_.cell_pressure_[c[j]];
276 state_.face_pressure_[face] /= double(num);
281 flow_solver_.setup(grid_, rock_, fluid_, wells_, gravity_, flow_bc_, &(state_.face_pressure_));
284 transport_solver_.setup(grid_, rock_, fluid_, wells_, flow_solver_.faceTransmissibilities(), gravity_);
287 src_.resize(grid_.numCells(), 0.0);
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]);
300 if (wells_.control(well) == Wells::Pressure) {
301 state_.well_bhp_pressure_.push_back(wells_.target(well));
303 int cell = wells_.wellCell(well, 0);
304 state_.well_bhp_pressure_.push_back(state_.cell_pressure_[cell][Fluid::Liquid]);
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_);
312 if (param.anyUnused()) {
313 std::cout <<
"***** WARNING: Unused parameters: *****\n";
314 param.displayUsage();
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());
323 param.writeParam(paramfilename);
332 template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
337 double voldisclimit = flow_solver_.volumeDiscrepancyLimit();
338 double stepsize = initial_stepsize_;
339 double current_time = 0.0;
341 std::vector<double> face_flux;
343 std::string output_name = output_dir_ +
"/" +
"blackoil-output";
344 while (current_time < total_time_) {
345 start_state = state_;
348 if (current_time + stepsize > total_time_) {
349 stepsize = total_time_ - current_time;
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;
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);
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;
370 state_ = start_state;
371 wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
374 assert(result == FlowSolver::SolveOk);
377 wells_.update(grid_.numCells(), state_.well_perf_pressure_, state_.well_perf_flux_);
380 bool voldisc_ok =
true;
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);
389 flow_solver_.volumeDiscrepancyAcceptable(state_.cell_pressure_, state_.face_pressure_,
390 state_.well_perf_pressure_, state_.cell_z_, 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;
398 std::cout <<
"Timestep was " << stepsize <<
" and max stepsize was " << max_dt << std::endl;
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);
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_);
416 std::cout <<
"********* Too large volume discrepancy: Shortening (pressure) stepsize, redoing step number " << step <<
" **********" << std::endl;
418 state_ = start_state;
419 wells_.update(grid_.numCells(), start_state.well_perf_pressure_, start_state.well_perf_flux_);
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);
430 if (!report_times_.empty()) {
431 if (current_time >= report_times_[step]) {
432 bool output_now = ((step + 1) % output_interval_ == 0);
434 output(grid_, fluid_, state_, face_flux, step, output_name);
437 if (step ==
int(report_times_.size())) {
441 stepsize = report_times_[step] - current_time;
443 bool output_now = ((step + 1) % output_interval_ == 0);
445 output(grid_, fluid_, state_, face_flux, step, output_name);
450 if (step % output_interval_ != 0) {
452 output(grid_, fluid_, state_, face_flux, step - 1, output_name);
461 template<
class Gr
id,
class Rock,
class Flu
id,
class Wells,
class FlowSolver,
class TransportSolver>
466 const State& simstate,
467 const std::vector<double>& face_flux,
469 const std::string& filebase)
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;
487 boost::filesystem::path fpath(filebase);
488 if (fpath.has_branch_path()) {
489 create_directories(fpath.branch_path());
493 std::vector<typename Grid::Vector> cell_velocity;
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(),
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());
509 Dune::VTKWriter<typename Grid::LeafGridView> vtkwriter(grid.leafView());
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),
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];
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];
535 std::string matlabdumpname(filebase +
"-");
536 matlabdumpname += boost::lexical_cast<std::string>(step);
537 matlabdumpname +=
".dat";
538 std::ofstream dump(matlabdumpname.c_str());
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];
545 std::copy(liq_press.begin(), liq_press.end(),
546 std::ostream_iterator<double>(dump,
" "));
549 for (
int comp = 0; comp < Fluid::numComponents; ++comp) {
550 std::copy(zv[comp].begin(), zv[comp].end(),
551 std::ostream_iterator<double>(dump,
" "));
555 for (
int phase = 0; phase < Fluid::numPhases; ++phase) {
556 std::copy(sv[phase].begin(), sv[phase].end(),
557 std::ostream_iterator<double>(dump,
" "));
561 std::copy(totflvol_dens.begin(), totflvol_dens.end(),
562 std::ostream_iterator<double>(dump,
" "));
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';
587 #endif // OPM_BLACKOILSIMULATOR_HEADER_INCLUDED A class for representing a flow boundary condition.
Definition: BoundaryConditions.hpp:121
Class for immiscible dead oil and dry gas.
Definition: applier.hpp:18
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