diff --git a/examples/Makefile.am b/examples/Makefile.am index a60d05ea..7a178da6 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -8,7 +8,8 @@ noinst_PROGRAMS = \ refine_wells \ scaneclipsedeck \ sim_wateroil \ -wells_example +wells_example \ +sim_2p_incomp_reorder refine_wells_SOURCES = refine_wells.cpp @@ -18,6 +19,12 @@ $(LDADD) $(LIBS) \ $(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ $(LAPACK_LIBS) $(LIBS) $(LIBS) +sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp +sim_2p_incomp_reorder_LDADD = \ +$(LDADD) $(LIBS) \ +$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB) \ +$(LAPACK_LIBS) $(LIBS) $(LIBS) + wells_example_SOURCES = wells_example.cpp if UMFPACK diff --git a/examples/sim_2p_incomp_reorder.cpp b/examples/sim_2p_incomp_reorder.cpp new file mode 100644 index 00000000..0b1ed7fc --- /dev/null +++ b/examples/sim_2p_incomp_reorder.cpp @@ -0,0 +1,220 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + + +// ----------------- Main program ----------------- +int +main(int argc, char** argv) +{ + using namespace Opm; + + std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n"; + Opm::parameter::ParameterGroup param(argc, argv, false); + std::cout << "--------------- Reading parameters ---------------" << std::endl; + + // If we have a "deck_filename", grid and props will be read from that. + bool use_deck = param.has("deck_filename"); + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr wells; + boost::scoped_ptr rock_comp; + Opm::SimulatorTimer simtimer; + Opm::TwophaseState state; + // bool check_well_controls = false; + // int max_well_control_iterations = 0; + double gravity[3] = { 0.0 }; + if (use_deck) { + std::string deck_filename = param.get("deck_filename"); + Opm::EclipseGridParser deck(deck_filename); + // Grid init + grid.reset(new Opm::GridManager(deck)); + // Rock and fluid init + const int* gc = grid->c_grid()->global_cell; + std::vector global_cell(gc, gc + grid->c_grid()->number_of_cells); + props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); + // Wells init. + wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + // Timer init. + if (deck.hasField("TSTEP")) { + simtimer.init(deck); + } else { + simtimer.init(param); + } + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(deck)); + // Gravity. + gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; + // Init state variables (saturation and pressure). + if (param.has("init_saturation")) { + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } else { + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + } + } else { + // Grid init. + const int nx = param.getDefault("nx", 100); + const int ny = param.getDefault("ny", 100); + const int nz = param.getDefault("nz", 1); + const double dx = param.getDefault("dx", 1.0); + const double dy = param.getDefault("dy", 1.0); + const double dz = param.getDefault("dz", 1.0); + grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); + // Rock and fluid init. + props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + // Wells init. + wells.reset(new Opm::WellsManager()); + // Timer init. + simtimer.init(param); + // Rock compressibility. + rock_comp.reset(new Opm::RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + } + + // Warn if gravity but no density difference. + bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); + if (use_gravity) { + if (props->density()[0] == props->density()[1]) { + std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl; + } + } + + // Source-related variables init. + int num_cells = grid->c_grid()->number_of_cells; + std::vector totmob; + std::vector omega; // Will remain empty if no gravity. + std::vector rc; // Will remain empty if no rock compressibility. + + // Extra rock init. + std::vector porevol; + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } else { + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + } + double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); + + // Initialising src + std::vector src(num_cells, 0.0); + if (wells->c_wells()) { + // Do nothing, wells will be the driving force, not source terms. + // Opm::wellsToSrc(*wells->c_wells(), num_cells, src); + } else { + const double default_injection = use_gravity ? 0.0 : 0.1; + const double flow_per_sec = param.getDefault("injected_porevolumes_per_day", default_injection) + *tot_porevol_init/Opm::unit::day; + src[0] = flow_per_sec; + src[num_cells - 1] = -flow_per_sec; + } + + // Boundary conditions. + Opm::FlowBCManager bcs; + if (param.getDefault("use_pside", false)) { + int pside = param.get("pside"); + double pside_pressure = param.get("pside_pressure"); + bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure); + } + + // Linear solver. + Opm::LinearSolverFactory linsolver(param); + + const double *grav = use_gravity ? &gravity[0] : 0; + + + Opm::SimulatorTwophase simulator(param, + *grid->c_grid(), + *props, + rock_comp->isActive() ? rock_comp.get() : 0, + wells->c_wells(), + src, + bcs.c_bcs(), + linsolver, + grav); + + // Warn if any parameters are unused. + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } + + // Write parameters used for later reference. + // if (output) { + // param.writeParam(output_dir + "/spu_2p.param"); + // } + + WellState well_state; + well_state.init(wells->c_wells(), state); + + simulator.run(simtimer, state, well_state); +} diff --git a/opm/core/simulator/SimulatorTwophase.cpp b/opm/core/simulator/SimulatorTwophase.cpp index a83b57b7..1a617a4b 100644 --- a/opm/core/simulator/SimulatorTwophase.cpp +++ b/opm/core/simulator/SimulatorTwophase.cpp @@ -17,8 +17,40 @@ along with OPM. If not, see . */ +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + #include #include +#include + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + namespace Opm { @@ -26,87 +58,404 @@ namespace Opm class SimulatorTwophase::Impl { public: - Impl() {} - void init(const parameter::ParameterGroup& param); - void run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state); + Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); + + void run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state); + private: -#if 0 + // Data. + + // Parameters for output. bool output_; std::string output_dir_; - int output_interval; + int output_interval_; + // Parameters for pressure solver. + int nl_pressure_maxiter_; + double nl_pressure_tolerance_; + // Parameters for transport solver. + int nl_maxiter_; + double nl_tolerance_; int num_transport_substeps_; - double gravity_[3]; - bool use_gravity_; - bool use_segregation_split = false; // - int nl_pressure_maxiter = 0; - double nl_pressure_tolerance = 0.0; - const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); - const int nl_maxiter = param.getDefault("nl_maxiter", 30); - - // boost::scoped_ptr grid_; - // boost::scoped_ptr props_; - // boost::scoped_ptr wells_; - // boost::scoped_ptr rock_comp_; - // Opm::SimulatorTimer simtimer_; - // Opm::TwophaseState state_; - // std::vector src(num_cells, 0.0); - // Opm::FlowBCManager bcs; - // Opm::LinearSolverFactory linsolver(param); + bool use_segregation_split_; + // Observed objects. + const UnstructuredGrid& grid_; + const IncompPropertiesInterface& props_; + const RockCompressibility* rock_comp_; + const Wells* wells_; + const std::vector& src_; + const FlowBoundaryConditions* bcs_; + const LinearSolverInterface& linsolver_; + const double* gravity_; + // Solvers + IncompTpfa psolver_; + TransportModelTwophase tsolver_; + // Needed by column-based gravity segregation solver. + std::vector< std::vector > columns_; + // Misc. data + std::vector allcells_; + }; + SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + { + pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity)); + } + + + + + + void SimulatorTwophase::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { + pimpl_->run(timer, state, well_state); + } + + + + static void outputState(const UnstructuredGrid& grid, + const Opm::TwophaseState& state, + const int step, + const std::string& output_dir) + { + // Write data in VTK format. + std::ostringstream vtkfilename; + vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; + std::ofstream vtkfile(vtkfilename.str().c_str()); + if (!vtkfile) { + THROW("Failed to open " << vtkfilename.str()); + } + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + std::vector cell_velocity; + Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); + dm["velocity"] = &cell_velocity; + Opm::writeVtkData(grid, dm, vtkfile); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } + } + + + static void outputWaterCut(const Opm::Watercut& watercut, + const std::string& output_dir) + { + // Write water cut curve. + std::string fname = output_dir + "/watercut.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + watercut.write(os); + } + + + static void outputWellReport(const Opm::WellReport& wellreport, + const std::string& output_dir) + { + // Write well report. + std::string fname = output_dir + "/wellreport.txt"; + std::ofstream os(fname.c_str()); + if (!os) { + THROW("Failed to open " << fname); + } + wellreport.write(os); + } + + + + + + SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity) + : grid_(grid), + props_(props), + rock_comp_(rock_comp), + wells_(wells), + src_(src), + bcs_(bcs), + linsolver_(linsolver), + gravity_(gravity), + psolver_(grid, props.permeability(), gravity, linsolver), + tsolver_(grid, props, 1e-9, 30) + { + // For output. + output_ = param.getDefault("output", true); + if (output_) { + output_dir_ = param.getDefault("output_dir", std::string("output")); + // Ensure that output dir exists + boost::filesystem::path fpath(output_dir_); + try { + create_directories(fpath); + } + catch (...) { + THROW("Creating directories failed: " << fpath); + } + output_interval_ = param.getDefault("output_interval", 1); + } + + // For pressure solver + nl_pressure_maxiter_ = param.getDefault("nl_pressure_maxiter", 10); + nl_pressure_tolerance_ = param.getDefault("nl_pressure_tolerance", 1.0); // Pascal + + // For transport solver. + nl_maxiter_ = param.getDefault("nl_maxiter", 30); + nl_tolerance_ = param.getDefault("nl_tolerance", 1e-9); + num_transport_substeps_ = param.getDefault("num_transport_substeps", 1); + use_segregation_split_ = param.getDefault("use_segregation_split", false); + if (gravity != 0 && use_segregation_split_){ + tsolver_.initGravity(gravity); + extractColumn(grid_, columns_); + } + + // Misc init. + const int num_cells = grid.number_of_cells; + allcells_.resize(num_cells); + for (int cell = 0; cell < num_cells; ++cell) { + allcells_[cell] = cell; + } + } + + + + + void SimulatorTwophase::Impl::run(SimulatorTimer& timer, + TwophaseState& state, + WellState& well_state) + { std::vector totmob; std::vector omega; // Will remain empty if no gravity. std::vector rc; // Will remain empty if no rock compressibility. - std::vector porevol; - double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); std::vector transport_src; - Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver); + // Initialisation. + std::vector porevol; + if (rock_comp_ && rock_comp_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + } else { + computePorevolume(grid_, props_.porosity(), porevol); + } + const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); - Opm::TransportModelTwophase tsolver(*grid->c_grid(), *props, nl_tolerance, nl_maxiter); - typedef std::pair, std::vector > > ColMap; - ColMap columns; + // Main simulation loop. + Opm::time::StopWatch pressure_timer; + double ptime = 0.0; + Opm::time::StopWatch transport_timer; + double ttime = 0.0; + Opm::time::StopWatch total_timer; + total_timer.start(); + std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl; + double init_satvol[2] = { 0.0 }; + double satvol[2] = { 0.0 }; + double injected[2] = { 0.0 }; + double produced[2] = { 0.0 }; + double tot_injected[2] = { 0.0 }; + double tot_produced[2] = { 0.0 }; + Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol); + std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init + << " " << init_satvol[1]/tot_porevol_init << std::endl; + Opm::Watercut watercut; + watercut.push(0.0, 0.0, 0.0); + Opm::WellReport wellreport; + std::vector fractional_flows; + std::vector well_resflows_phase; + int num_wells = 0; + if (wells_) { + num_wells = wells_->number_of_wells; + well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); + wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates()); + } + const int num_cells = grid_.number_of_cells; + for (; !timer.done(); ++timer) { + // Report timestep and (optionally) write state to disk. + timer.report(std::cout); + if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + } - std::vector allcells(num_cells); -#endif - }; + // Solve pressure. + if (gravity_) { + computeTotalMobilityOmega(props_, allcells_, state.saturation(), totmob, omega); + } else { + computeTotalMobility(props_, allcells_, state.saturation(), totmob); + } + std::vector wdp; + if (wells_) { + Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(), + gravity_ ? gravity_[2] : 0.0, true, wdp); + } + do { + pressure_timer.start(); + if (rock_comp_ && rock_comp_->isActive()) { + rc.resize(num_cells); + std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(grid_, props_.porosity(), *rock_comp_, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells + num_wells); + std::vector prev_pressure(num_cells + num_wells); + for (int iter = 0; iter < nl_pressure_maxiter_; ++iter) { + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp_->rockComp(state.pressure()[cell]); + } + computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol); + std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); + std::copy(well_state.bhp().begin(), well_state.bhp().end(), prev_pressure.begin() + num_cells); + // prev_pressure = state.pressure(); - SimulatorTwophase::SimulatorTwophase() - { - pimpl_.reset(new Impl); - } + // compute pressure increment + psolver_.solveIncrement(totmob, omega, src_, wdp, bcs_, porevol, rc, + prev_pressure, initial_porevolume, timer.currentStepLength(), + pressure_increment); + + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); + } + for (int well = 0; well < num_wells; ++well) { + well_state.bhp()[well] += pressure_increment[num_cells + well]; + max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well])); + } + + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance_) { + break; + } + } + psolver_.computeFaceFlux(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } else { + psolver_.solve(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(), + well_state.bhp(), well_state.perfRates()); + } + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + } while (false); + + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, + wells_, well_state.perfRates(), transport_src); + + // Solve transport. + transport_timer.start(); + double stepsize = timer.currentStepLength(); + if (num_transport_substeps_ != 1) { + stepsize /= double(num_transport_substeps_); + std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl; + } + for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { + tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0], + stepsize, state.saturation()); + Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced); + if (use_segregation_split_) { + tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation()); + } + } + transport_timer.stop(); + double tt = transport_timer.secsSinceStart(); + std::cout << "Transport solver took: " << tt << " seconds." << std::endl; + ttime += tt; + + // Report volume balances. + Opm::computeSaturatedVol(porevol, state.saturation(), satvol); + tot_injected[0] += injected[0]; + tot_injected[1] += injected[1]; + tot_produced[0] += produced[0]; + tot_produced[1] += produced[1]; + std::cout.precision(5); + const int width = 18; + std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n"; + std::cout << " Saturated volumes: " + << std::setw(width) << satvol[0]/tot_porevol_init + << std::setw(width) << satvol[1]/tot_porevol_init << std::endl; + std::cout << " Injected volumes: " + << std::setw(width) << injected[0]/tot_porevol_init + << std::setw(width) << injected[1]/tot_porevol_init << std::endl; + std::cout << " Produced volumes: " + << std::setw(width) << produced[0]/tot_porevol_init + << std::setw(width) << produced[1]/tot_porevol_init << std::endl; + std::cout << " Total inj volumes: " + << std::setw(width) << tot_injected[0]/tot_porevol_init + << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl; + std::cout << " Total prod volumes: " + << std::setw(width) << tot_produced[0]/tot_porevol_init + << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl; + std::cout << " In-place + prod - inj: " + << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init + << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl; + std::cout << " Init - now - pr + inj: " + << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init + << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init + << std::endl; + std::cout.precision(8); + + watercut.push(timer.currentTime() + timer.currentStepLength(), + produced[0]/(produced[0] + produced[1]), + tot_produced[0]/tot_porevol_init); + if (wells_) { + wellreport.push(props_, *wells_, state.saturation(), + timer.currentTime() + timer.currentStepLength(), + well_state.bhp(), well_state.perfRates()); + } + } + total_timer.stop(); + + std::cout << "\n\n================ End of simulation ===============\n" + << "Total time taken: " << total_timer.secsSinceStart() + << "\n Pressure time: " << ptime + << "\n Transport time: " << ttime << std::endl; + + if (output_) { + outputState(grid_, state, timer.currentStepNum(), output_dir_); + outputWaterCut(watercut, output_dir_); + if (wells_) { + outputWellReport(wellreport, output_dir_); + } + } - void SimulatorTwophase::init(const parameter::ParameterGroup& param) - { - pimpl_->init(param); - } - - - - void SimulatorTwophase::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) - { - pimpl_->run(timer, wells, state, well_state); - } - - - - void SimulatorTwophase::Impl::init(const parameter::ParameterGroup& /*param*/) - { - } - - - - void SimulatorTwophase::Impl::run(const SimulatorTimer& timer, const Wells& wells, TwophaseState& state, WellState& well_state) - { - (void) timer; - (void) wells; - (void) state; - (void) well_state; } } // namespace Opm diff --git a/opm/core/simulator/SimulatorTwophase.hpp b/opm/core/simulator/SimulatorTwophase.hpp index 636d4f94..053a3196 100644 --- a/opm/core/simulator/SimulatorTwophase.hpp +++ b/opm/core/simulator/SimulatorTwophase.hpp @@ -20,14 +20,19 @@ #ifndef OPM_SIMULATORTWOPHASE_HEADER_INCLUDED #define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED -#include +#include +#include struct UnstructuredGrid; struct Wells; +struct FlowBoundaryConditions; namespace Opm { namespace parameter { class ParameterGroup; } + class IncompPropertiesInterface; + class RockCompressibility; + class LinearSolverInterface; class SimulatorTimer; class TwophaseState; class WellState; @@ -36,25 +41,53 @@ namespace Opm class SimulatorTwophase { public: - /// Default constructor. - SimulatorTwophase(); - - /// Initialise from parameters. - void init(const parameter::ParameterGroup& param); + /// Initialise from parameters and objects to observe. + /// \param[in] param parameters, this class accepts the following: + /// parameter (default) effect + /// ----------------------------------------------------------- + /// output (true) write output to files? + /// output_dir ("output") output directoty + /// output_interval (1) output every nth step + /// nl_pressure_maxiter (10) max nonlinear iterations in pressure + /// nl_pressure_tolerance (1.0) pressure solver nonlinear tolerance (in Pascal) + /// nl_maxiter (30) max nonlinear iterations in transport + /// nl_tolerance (1e-9) transport solver absolute residual tolerance + /// num_transport_substeps (1) number of transport steps per pressure step + /// use_segregation_split (false) solve for gravity segregation (if false, + /// segregation is ignored). + /// + /// \param[in] grid grid data structure + /// \param[in] props fluid and rock properties + /// \param[in] rock_comp if non-null, rock compressibility properties + /// \param[in] wells if non-null, wells data structure + /// \param[in] src source terms + /// \param[in] bcs boundary conditions, treat as all noflow if null + /// \param[in] linsolver linear solver + /// \param[in] gravity if non-null, gravity vector + SimulatorTwophase(const parameter::ParameterGroup& param, + const UnstructuredGrid& grid, + const IncompPropertiesInterface& props, + const RockCompressibility* rock_comp, + const Wells* wells, + const std::vector& src, + const FlowBoundaryConditions* bcs, + const LinearSolverInterface& linsolver, + const double* gravity); /// Run the simulation. - /// \param[in] timer governs the requested reporting timesteps - /// \param[in] wells data structure for wells - /// \param[out] state state of reservoir: pressure, fluxes - /// \param[out] well_state state of wells: bhp, perforation rates - void run(const SimulatorTimer& timer, - const Wells& wells, + /// This will run succesive timesteps until timer.done() is true. It will + /// modify the reservoir and well states. + /// \param[in,out] timer governs the requested reporting timesteps + /// \param[in,out] state state of reservoir: pressure, fluxes + /// \param[in,out] well_state state of wells: bhp, perforation rates + void run(SimulatorTimer& timer, TwophaseState& state, WellState& well_state); private: class Impl; - boost::scoped_ptr pimpl_; + // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. + boost::shared_ptr pimpl_; }; } // namespace Opm