27 template <
class Implementation>
32 const RockCompressibility* rock_comp_props,
34 const double* gravity,
35 const bool has_disgas,
36 const bool has_vapoil,
37 std::shared_ptr<EclipseState> eclipse_state,
39 const std::vector<double>& threshold_pressures_by_face)
45 rock_comp_props_(rock_comp_props),
49 has_disgas_(has_disgas),
50 has_vapoil_(has_vapoil),
51 terminal_output_(param.getDefault(
"output_terminal", true)),
52 eclipse_state_(eclipse_state),
53 output_writer_(output_writer),
54 rateConverter_(props_,
std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
55 threshold_pressures_by_face_(threshold_pressures_by_face),
56 is_parallel_run_( false )
59 const int num_cells = AutoDiffGrid::numCells(grid);
61 for (
int cell = 0; cell < num_cells; ++cell) {
67 const ParallelISTLInformation& info =
76 template <
class Implementation>
83 Opm::time::StopWatch solver_timer;
85 Opm::time::StopWatch step_timer;
86 Opm::time::StopWatch total_timer;
88 std::string tstep_filename = output_writer_.outputDirectory() +
"/step_timing.txt";
89 std::ofstream tstep_os(tstep_filename.c_str());
92 std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
93 if( param_.getDefault(
"timestep.adaptive",
true ) )
95 adaptiveTimeStepping.reset(
new AdaptiveTimeStepping( param_, solver_.parallelInformation(), terminal_output_ ) );
99 output_writer_.writeInit( timer );
101 std::string restorefilename = param_.getDefault(
"restorefile", std::string(
"") );
102 if( ! restorefilename.empty() )
105 const int desiredRestoreStep = param_.getDefault(
"restorestep",
int(-1) );
106 output_writer_.restore( timer, state, prev_well_state, restorefilename, desiredRestoreStep );
109 unsigned int totalNewtonIterations = 0;
110 unsigned int totalLinearIterations = 0;
113 while (!timer.done()) {
116 if ( terminal_output_ )
118 timer.report(std::cout);
122 WellsManager wells_manager(eclipse_state_,
123 timer.currentStepNum(),
124 Opm::UgGridHelpers::numCells(grid_),
125 Opm::UgGridHelpers::globalCell(grid_),
126 Opm::UgGridHelpers::cartDims(grid_),
127 Opm::UgGridHelpers::dimensions(grid_),
128 Opm::UgGridHelpers::cell2Faces(grid_),
129 Opm::UgGridHelpers::beginFaceCentroids(grid_),
130 props_.permeability(),
132 const Wells* wells = wells_manager.c_wells();
134 well_state.init(wells, state, prev_well_state);
137 asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
140 output_writer_.writeTimeStep( timer, state, well_state );
143 props_.updateSatOilMax(state.saturation());
144 props_.updateSatHyst(state.saturation(), allcells_);
147 asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
150 solver_timer.start();
152 auto solver = asImpl().createSolver(wells);
159 if( adaptiveTimeStepping ) {
160 adaptiveTimeStepping->step( timer, *solver, state, well_state, output_writer_ );
164 solver->step(timer.currentStepLength(), state, well_state);
171 totalNewtonIterations += solver->newtonIterations();
172 totalLinearIterations += solver->linearIterations();
175 const double st = solver_timer.secsSinceStart();
177 if ( terminal_output_ )
179 std::cout <<
"Fully implicit solver took: " << st <<
" seconds." << std::endl;
183 if ( output_writer_.output() ) {
184 SimulatorReport step_report;
185 step_report.pressure_time = st;
186 step_report.total_time = step_timer.secsSinceStart();
187 step_report.reportParam(tstep_os);
192 prev_well_state = well_state;
196 output_writer_.writeTimeStep( timer, state, prev_well_state );
200 SimulatorReport report;
201 report.pressure_time = stime;
202 report.transport_time = 0.0;
203 report.total_time = total_timer.secsSinceStart();
204 report.total_newton_iterations = totalNewtonIterations;
205 report.total_linear_iterations = totalLinearIterations;
209 namespace SimFIBODetails {
210 typedef std::unordered_map<std::string, WellConstPtr>
WellMap;
217 for (std::vector<WellConstPtr>::const_iterator
218 w = wells.begin(), e = wells.end();
221 wmap.insert(std::make_pair((*w)->name(), *w));
230 int i, n = well_controls_get_num(ctrl);
233 for (i = 0; (! match) && (i < n); ++i) {
234 match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
237 if (! match) { i = 0; }
251 const std::string& name,
252 const std::size_t step)
256 WellMap::const_iterator i = wmap.find(name);
258 if (i != wmap.end()) {
259 WellConstPtr wp = i->second;
261 match = (wp->isProducer(step) &&
262 wp->getProductionProperties(step)
263 .hasProductionControl(WellProducer::RESV))
264 || (wp->isInjector(step) &&
265 wp->getInjectionProperties(step)
266 .hasInjectionControl(WellInjector::RESV));
272 inline std::vector<int>
274 const std::size_t step,
277 std::vector<int> resv_wells;
280 for (
int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
282 ((wells->name[w] != 0) &&
283 is_resv(wmap, wells->name[w], step)))
285 resv_wells.push_back(w);
295 const WellProductionProperties& p,
296 std::vector<double>& rates)
298 assert (! p.predictionMode);
299 assert (rates.size() ==
300 std::vector<double>::size_type(pu.num_phases));
302 if (pu.phase_used[ BlackoilPhases::Aqua ]) {
303 const std::vector<double>::size_type
304 i = pu.phase_pos[ BlackoilPhases::Aqua ];
306 rates[i] = p.WaterRate;
309 if (pu.phase_used[ BlackoilPhases::Liquid ]) {
310 const std::vector<double>::size_type
311 i = pu.phase_pos[ BlackoilPhases::Liquid ];
313 rates[i] = p.OilRate;
316 if (pu.phase_used[ BlackoilPhases::Vapour ]) {
317 const std::vector<double>::size_type
318 i = pu.phase_pos[ BlackoilPhases::Vapour ];
320 rates[i] = p.GasRate;
325 template <
class Implementation>
332 template <
class Implementation>
334 -> std::unique_ptr<Solver>
336 auto model = std::unique_ptr<Model>(
new Model(model_param_,
348 if (!threshold_pressures_by_face_.empty()) {
349 model->setThresholdPressures(threshold_pressures_by_face_);
352 return std::unique_ptr<Solver>(
new Solver(solver_param_, std::move(model)));
355 template <
class Implementation>
358 const BlackoilState& x,
363 const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
368 if (! resv_wells.empty()) {
369 const PhaseUsage& pu = props_.phaseUsage();
370 const std::vector<double>::size_type np = props_.numPhases();
372 rateConverter_.defineState(x);
374 std::vector<double> distr (np);
375 std::vector<double> hrates(np);
376 std::vector<double> prates(np);
378 for (std::vector<int>::const_iterator
379 rp = resv_wells.begin(), e = resv_wells.end();
382 WellControls* ctrl = wells->ctrls[*rp];
383 const bool is_producer = wells->type[*rp] == PRODUCER;
390 const std::vector<double>::size_type off = (*rp) * np;
395 std::transform(xw.wellRates().begin() + (off + 0*np),
396 xw.wellRates().begin() + (off + 1*np),
397 prates.begin(), std::negate<double>());
399 std::copy(xw.wellRates().begin() + (off + 0*np),
400 xw.wellRates().begin() + (off + 1*np),
404 const int fipreg = 0;
405 rateConverter_.calcCoeff(prates, fipreg, distr);
407 well_controls_iset_distr(ctrl, rctrl, & distr[0]);
413 if (is_producer && wells->name[*rp] != 0) {
414 WellMap::const_iterator i = wmap.find(wells->name[*rp]);
416 if (i != wmap.end()) {
417 WellConstPtr wp = i->second;
419 const WellProductionProperties& p =
420 wp->getProductionProperties(step);
422 if (! p.predictionMode) {
426 const int fipreg = 0;
427 rateConverter_.calcCoeff(hrates, fipreg, distr);
433 const double target =
434 - std::inner_product(distr.begin(), distr.end(),
435 hrates.begin(), 0.0);
437 well_controls_clear(ctrl);
438 well_controls_assert_number_of_phases(ctrl,
int(np));
440 static const double invalid_alq = -std::numeric_limits<double>::max();
441 static const int invalid_vfp = -std::numeric_limits<int>::max();
444 well_controls_add_new(RESERVOIR_RATE, target,
445 invalid_alq, invalid_vfp,
450 double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
452 well_controls_add_new(BHP, bhp_limit,
453 invalid_alq, invalid_vfp,
456 if (ok_resv != 0 && ok_bhp != 0) {
457 xw.currentControls()[*rp] = 0;
458 well_controls_set_current(ctrl, 0);
468 for (
int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
469 WellControls* ctrl = wells->ctrls[w];
470 const bool is_producer = wells->type[w] == PRODUCER;
471 if (!is_producer && wells->name[w] != 0) {
472 WellMap::const_iterator i = wmap.find(wells->name[w]);
473 if (i != wmap.end()) {
474 WellConstPtr wp = i->second;
475 const WellInjectionProperties& injector = wp->getInjectionProperties(step);
476 if (!injector.predictionMode) {
478 static const double invalid_alq = -std::numeric_limits<double>::max();
479 static const int invalid_vfp = -std::numeric_limits<int>::max();
482 double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
484 well_controls_add_new(BHP, bhp_limit,
485 invalid_alq, invalid_vfp,
488 OPM_THROW(std::runtime_error,
"Failed to add well control.");
Definition: GeoProps.hpp:53
Definition: BlackoilPropsAdInterface.hpp:38
std::vector< int > resvWells(const Wells *wells, const std::size_t step, const WellMap &wmap)
Definition: SimulatorBase_impl.hpp:273
Definition: AdditionalObjectDeleter.hpp:22
Traits::Grid Grid
Definition: SimulatorBase.hpp:92
bool terminal_output_
Definition: SimulatorBase.hpp:186
SimulatorBase(const parameter::ParameterGroup ¶m, const Grid &grid, const DerivedGeology &geo, BlackoilPropsAdInterface &props, const RockCompressibility *rock_comp_props, NewtonIterationBlackoilInterface &linsolver, const double *gravity, const bool disgas, const bool vapoil, std::shared_ptr< EclipseState > eclipse_state, OutputWriter &output_writer, const std::vector< double > &threshold_pressures_by_face)
Definition: SimulatorBase_impl.hpp:28
bool is_resv(const Wells &wells, const int w)
Definition: SimulatorBase_impl.hpp:243
void historyRates(const PhaseUsage &pu, const WellProductionProperties &p, std::vector< double > &rates)
Definition: SimulatorBase_impl.hpp:294
std::unique_ptr< Solver > createSolver(const Wells *wells)
Definition: SimulatorBase_impl.hpp:333
Traits::Solver Solver
Definition: SimulatorBase.hpp:93
bool is_parallel_run_
Definition: SimulatorBase.hpp:195
void handleAdditionalWellInflow(SimulatorTimer &timer, WellsManager &wells_manager, WellState &well_state, const Wells *wells)
Definition: SimulatorBase_impl.hpp:326
virtual const boost::any & parallelInformation() const =0
Get the information about the parallelization of the grid.
Traits::OutputWriter OutputWriter
Definition: SimulatorBase.hpp:91
NewtonIterationBlackoilInterface & solver_
Definition: SimulatorBase.hpp:181
Traits::ReservoirState ReservoirState
Definition: SimulatorBase.hpp:89
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
std::vector< int > allcells_
Definition: SimulatorBase.hpp:183
void computeRESV(const std::size_t step, const Wells *wells, const BlackoilState &x, WellState &xw)
Definition: SimulatorBase_impl.hpp:356
Traits::Model Model
Definition: SimulatorBase.hpp:166
std::unordered_map< std::string, WellConstPtr > WellMap
Definition: SimulatorBase_impl.hpp:210
int resv_control(const WellControls *ctrl)
Definition: SimulatorBase_impl.hpp:228
SimulatorReport run(SimulatorTimer &timer, ReservoirState &state)
Definition: SimulatorBase_impl.hpp:77
WellMap mapWells(const std::vector< WellConstPtr > &wells)
Definition: SimulatorBase_impl.hpp:213
Traits::WellState WellState
Definition: SimulatorBase.hpp:90