00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <utility>
00023 #include <functional>
00024 #include <algorithm>
00025 #include <locale>
00026 #include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
00027 #include <opm/core/utility/initHydroCarbonState.hpp>
00028 #include <opm/core/well_controls.h>
00029 #include <opm/core/wells/DynamicListEconLimited.hpp>
00030 #include <opm/autodiff/BlackoilModel.hpp>
00031
00032 namespace Opm
00033 {
00034
00035 template <class Implementation>
00036 SimulatorBase<Implementation>::SimulatorBase(const ParameterGroup& param,
00037 const Grid& grid,
00038 DerivedGeology& geo,
00039 BlackoilPropsAdFromDeck& props,
00040 const RockCompressibility* rock_comp_props,
00041 NewtonIterationBlackoilInterface& linsolver,
00042 const double* gravity,
00043 const bool has_disgas,
00044 const bool has_vapoil,
00045 std::shared_ptr<EclipseState> eclipse_state,
00046 OutputWriter& output_writer,
00047 const std::vector<double>& threshold_pressures_by_face,
00048 const std::unordered_set<std::string>& defunct_well_names)
00049 : param_(param),
00050 model_param_(param),
00051 solver_param_(param),
00052 grid_(grid),
00053 props_(props),
00054 rock_comp_props_(rock_comp_props),
00055 gravity_(gravity),
00056 geo_(geo),
00057 solver_(linsolver),
00058 has_disgas_(has_disgas),
00059 has_vapoil_(has_vapoil),
00060 terminal_output_(param.getDefault("output_terminal", true)),
00061 eclipse_state_(eclipse_state),
00062 output_writer_(output_writer),
00063 rateConverter_(props_.phaseUsage(), std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
00064 threshold_pressures_by_face_(threshold_pressures_by_face),
00065 is_parallel_run_( false ),
00066 defunct_well_names_(defunct_well_names)
00067 {
00068
00069 const int num_cells = AutoDiffGrid::numCells(grid);
00070 allcells_.resize(num_cells);
00071 for (int cell = 0; cell < num_cells; ++cell) {
00072 allcells_[cell] = cell;
00073 }
00074 #if HAVE_MPI
00075 if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00076 {
00077 const ParallelISTLInformation& info =
00078 boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
00079
00080 terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
00081 is_parallel_run_ = ( info.communicator().size() > 1 );
00082 }
00083 #endif
00084 }
00085
00086 template <class Implementation>
00087 SimulatorReport SimulatorBase<Implementation>::run(SimulatorTimer& timer,
00088 ReservoirState& state)
00089 {
00090 WellState prev_well_state;
00091
00092 ExtraData extra;
00093 if (output_writer_.isRestart()) {
00094
00095 output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state, extra);
00096 initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid_), has_disgas_, has_vapoil_);
00097 initHysteresisParams(state);
00098 }
00099
00100
00101 Opm::time::StopWatch solver_timer;
00102 Opm::time::StopWatch step_timer;
00103 Opm::time::StopWatch total_timer;
00104 total_timer.start();
00105 std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
00106 std::ofstream tstep_os;
00107
00108 if ( output_writer_.output() ) {
00109 if ( output_writer_.isIORank() )
00110 tstep_os.open(tstep_filename.c_str());
00111 }
00112
00113 const auto& schedule = eclipse_state_->getSchedule();
00114
00115
00116 const auto& events = schedule.getEvents();
00117 std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
00118 if( param_.getDefault("timestep.adaptive", true ) )
00119 {
00120
00121 if (param_.getDefault("use_TUNING", false)) {
00122 adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule.getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
00123 } else {
00124 adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
00125 }
00126 if (output_writer_.isRestart()) {
00127 if (extra.suggested_step > 0.0) {
00128 adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
00129 }
00130 }
00131 }
00132
00133 std::string restorefilename = param_.getDefault("restorefile", std::string("") );
00134 if( ! restorefilename.empty() )
00135 {
00136
00137 const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
00138
00139 output_writer_.restore( timer,
00140 state,
00141 prev_well_state,
00142 restorefilename,
00143 desiredRestoreStep );
00144 }
00145
00146 DynamicListEconLimited dynamic_list_econ_limited;
00147 SimulatorReport report;
00148 SimulatorReport stepReport;
00149
00150 bool ooip_computed = false;
00151 std::vector<int> fipnum_global = eclipse_state_->get3DProperties().getIntGridProperty("FIPNUM").getData();
00152
00153 std::vector<int> fipnum(AutoDiffGrid::numCells(grid_));
00154 if (fipnum_global.empty()) {
00155 std::fill(fipnum.begin(), fipnum.end(), 0);
00156 } else {
00157 for (size_t c = 0; c < fipnum.size(); ++c) {
00158 fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]];
00159 }
00160 }
00161 std::vector<std::vector<double> > OOIP;
00162
00163 while (!timer.done()) {
00164
00165 step_timer.start();
00166 if ( terminal_output_ )
00167 {
00168 std::ostringstream ss;
00169 timer.report(ss);
00170 OpmLog::note(ss.str());
00171 }
00172
00173
00174 WellsManager wells_manager(*eclipse_state_,
00175 timer.currentStepNum(),
00176 Opm::UgGridHelpers::numCells(grid_),
00177 Opm::UgGridHelpers::globalCell(grid_),
00178 Opm::UgGridHelpers::cartDims(grid_),
00179 Opm::UgGridHelpers::dimensions(grid_),
00180 Opm::UgGridHelpers::cell2Faces(grid_),
00181 Opm::UgGridHelpers::beginFaceCentroids(grid_),
00182 dynamic_list_econ_limited,
00183 is_parallel_run_,
00184 defunct_well_names_);
00185 const Wells* wells = wells_manager.c_wells();
00186 WellState well_state;
00187 well_state.init(wells, state, prev_well_state, props_.phaseUsage());
00188
00189
00190 asImpl().handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
00191
00192
00193 if (timer.initialStep()) {
00194 Dune::Timer perfTimer;
00195 perfTimer.start();
00196
00197
00198
00199 output_writer_.writeTimeStepWithoutCellProperties( timer, state, well_state, {}, {} );
00200
00201 report.output_write_time += perfTimer.stop();
00202 }
00203
00204
00205 props_.updateSatOilMax(state.saturation());
00206 props_.updateSatHyst(state.saturation(), allcells_);
00207
00208
00209 asImpl().computeRESV(timer.currentStepNum(), wells, state, well_state);
00210
00211
00212 solver_timer.start();
00213
00214 const WellModel well_model(wells, &(wells_manager.wellCollection()));
00215
00216 std::unique_ptr<Solver> solver = asImpl().createSolver(well_model);
00217
00218
00219 if (!ooip_computed) {
00220 OOIP = solver->computeFluidInPlace(state, fipnum);
00221 FIPUnitConvert(eclipse_state_->getUnits(), OOIP);
00222 ooip_computed = true;
00223 }
00224
00225 if( terminal_output_ )
00226 {
00227 std::ostringstream step_msg;
00228 boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
00229 step_msg.imbue(std::locale(std::locale::classic(), facet));
00230 step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
00231 << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
00232 << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
00233 << ", date = " << timer.currentDateTime();
00234 OpmLog::info(step_msg.str());
00235 }
00236
00237
00238
00239
00240
00241
00242 if( adaptiveTimeStepping ) {
00243 bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
00244 events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
00245 events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
00246 events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
00247 report += adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
00248 output_writer_.requireFIPNUM() ? &fipnum : nullptr );
00249 }
00250 else {
00251
00252 stepReport = solver->step(timer, state, well_state);
00253 report += stepReport;
00254
00255 if( terminal_output_ )
00256 {
00257 std::ostringstream iter_msg;
00258 iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
00259 if (solver->wellIterations() != 0) {
00260 iter_msg << " days well iterations = " << solver->wellIterations() << ", ";
00261 }
00262 iter_msg << "non-linear iterations = " << solver->nonlinearIterations()
00263 << ", total linear iterations = " << solver->linearIterations()
00264 << "\n";
00265 OpmLog::info(iter_msg.str());
00266 }
00267 }
00268
00269
00270
00271 const int nextTimeStepIdx = timer.currentStepNum() + 1;
00272 if (nextTimeStepIdx < timer.numSteps()
00273 && events.hasEvent(ScheduleEvents::GEO_MODIFIER, nextTimeStepIdx)) {
00274
00275
00276
00277
00278 const auto& miniDeck = schedule.getModifierDeck(nextTimeStepIdx);
00279 eclipse_state_->applyModifierDeck(miniDeck);
00280 geo_.update(grid_, props_, *eclipse_state_, gravity_);
00281 }
00282
00283
00284 solver_timer.stop();
00285
00286
00287 report.solver_time += solver_timer.secsSinceStart();
00288
00289
00290 std::vector<std::vector<double> > COIP;
00291 COIP = solver->computeFluidInPlace(state, fipnum);
00292 std::vector<double> OOIP_totals = FIPTotals(OOIP, state);
00293 std::vector<double> COIP_totals = FIPTotals(COIP, state);
00294
00295
00296 FIPUnitConvert(eclipse_state_->getUnits(), COIP);
00297 FIPUnitConvert(eclipse_state_->getUnits(), OOIP_totals);
00298 FIPUnitConvert(eclipse_state_->getUnits(), COIP_totals);
00299
00300 if ( terminal_output_ )
00301 {
00302 outputFluidInPlace(OOIP_totals, COIP_totals,eclipse_state_->getUnits(), 0);
00303 for (size_t reg = 0; reg < OOIP.size(); ++reg) {
00304 outputFluidInPlace(OOIP[reg], COIP[reg], eclipse_state_->getUnits(), reg+1);
00305 }
00306 }
00307
00308 if ( terminal_output_ )
00309 {
00310 std::string msg;
00311 msg = "Fully implicit solver took: " + std::to_string(stepReport.solver_time) + " seconds. Total solver time taken: " + std::to_string(report.solver_time) + " seconds.";
00312 OpmLog::note(msg);
00313 }
00314
00315 if ( tstep_os.is_open() ) {
00316 stepReport.reportParam(tstep_os);
00317 }
00318
00319
00320 ++timer;
00321
00322
00323 Dune::Timer perfTimer;
00324 perfTimer.start();
00325 const auto& physicalModel = solver->model();
00326 output_writer_.writeTimeStep( timer, state, well_state, physicalModel );
00327 report.output_write_time += perfTimer.stop();
00328
00329 prev_well_state = well_state;
00330
00331 asImpl().updateListEconLimited(solver, eclipse_state_->getSchedule(), timer.currentStepNum(), wells,
00332 well_state, dynamic_list_econ_limited);
00333 }
00334
00335
00336 total_timer.stop();
00337 report.total_time = total_timer.secsSinceStart();
00338 report.converged = true;
00339 return report;
00340 }
00341
00342 namespace SimFIBODetails {
00343 typedef std::unordered_map<std::string, const Well* > WellMap;
00344
00345 inline WellMap
00346 mapWells(const std::vector< const Well* >& wells)
00347 {
00348 WellMap wmap;
00349
00350 for (std::vector< const Well* >::const_iterator
00351 w = wells.begin(), e = wells.end();
00352 w != e; ++w)
00353 {
00354 wmap.insert(std::make_pair((*w)->name(), *w));
00355 }
00356
00357 return wmap;
00358 }
00359
00360 inline int
00361 resv_control(const WellControls* ctrl)
00362 {
00363 int i, n = well_controls_get_num(ctrl);
00364
00365 bool match = false;
00366 for (i = 0; (! match) && (i < n); ++i) {
00367 match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
00368 }
00369
00370 if (! match) { i = 0; }
00371
00372 return i - 1;
00373 }
00374
00375 inline bool
00376 is_resv(const Wells& wells,
00377 const int w)
00378 {
00379 return (0 <= resv_control(wells.ctrls[w]));
00380 }
00381
00382 inline bool
00383 is_resv(const WellMap& wmap,
00384 const std::string& name,
00385 const std::size_t step)
00386 {
00387 bool match = false;
00388
00389 WellMap::const_iterator i = wmap.find(name);
00390
00391 if (i != wmap.end()) {
00392 const Well* wp = i->second;
00393
00394 match = (wp->isProducer(step) &&
00395 wp->getProductionProperties(step)
00396 .hasProductionControl(WellProducer::RESV))
00397 || (wp->isInjector(step) &&
00398 wp->getInjectionProperties(step)
00399 .hasInjectionControl(WellInjector::RESV));
00400 }
00401
00402 return match;
00403 }
00404
00405 inline std::vector<int>
00406 resvWells(const Wells* wells,
00407 const std::size_t step,
00408 const WellMap& wmap)
00409 {
00410 std::vector<int> resv_wells;
00411 if( wells )
00412 {
00413 for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
00414 if (is_resv(*wells, w) ||
00415 ((wells->name[w] != 0) &&
00416 is_resv(wmap, wells->name[w], step)))
00417 {
00418 resv_wells.push_back(w);
00419 }
00420 }
00421 }
00422
00423 return resv_wells;
00424 }
00425
00426 inline void
00427 historyRates(const PhaseUsage& pu,
00428 const WellProductionProperties& p,
00429 std::vector<double>& rates)
00430 {
00431 assert (! p.predictionMode);
00432 assert (rates.size() ==
00433 std::vector<double>::size_type(pu.num_phases));
00434
00435 if (pu.phase_used[ BlackoilPhases::Aqua ]) {
00436 const std::vector<double>::size_type
00437 i = pu.phase_pos[ BlackoilPhases::Aqua ];
00438
00439 rates[i] = p.WaterRate;
00440 }
00441
00442 if (pu.phase_used[ BlackoilPhases::Liquid ]) {
00443 const std::vector<double>::size_type
00444 i = pu.phase_pos[ BlackoilPhases::Liquid ];
00445
00446 rates[i] = p.OilRate;
00447 }
00448
00449 if (pu.phase_used[ BlackoilPhases::Vapour ]) {
00450 const std::vector<double>::size_type
00451 i = pu.phase_pos[ BlackoilPhases::Vapour ];
00452
00453 rates[i] = p.GasRate;
00454 }
00455 }
00456 }
00457
00458 template <class Implementation>
00459 void SimulatorBase<Implementation>::handleAdditionalWellInflow(SimulatorTimer& ,
00460 WellsManager& ,
00461 WellState& ,
00462 const Wells* )
00463 { }
00464
00465 template <class Implementation>
00466 auto SimulatorBase<Implementation>::createSolver(const WellModel& well_model)
00467 -> std::unique_ptr<Solver>
00468 {
00469 auto model = std::unique_ptr<Model>(new Model(model_param_,
00470 grid_,
00471 props_,
00472 geo_,
00473 rock_comp_props_,
00474 well_model,
00475 solver_,
00476 eclipse_state_,
00477 has_disgas_,
00478 has_vapoil_,
00479 terminal_output_));
00480
00481 if (!threshold_pressures_by_face_.empty()) {
00482 model->setThresholdPressures(threshold_pressures_by_face_);
00483 }
00484
00485 return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
00486 }
00487
00488
00489 template <class Implementation>
00490 void SimulatorBase<Implementation>::computeRESV(const std::size_t step,
00491 const Wells* wells,
00492 const BlackoilState& x,
00493 WellState& xw)
00494 {
00495 typedef SimFIBODetails::WellMap WellMap;
00496
00497 const auto w_ecl = eclipse_state_->getSchedule().getWells(step);
00498 const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
00499
00500 const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
00501
00502 const std::size_t number_resv_wells = resv_wells.size();
00503 std::size_t global_number_resv_wells = number_resv_wells;
00504 #if HAVE_MPI
00505 if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00506 {
00507 const auto& info =
00508 boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
00509 global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
00510 if ( global_number_resv_wells )
00511 {
00512
00513
00514
00515
00516 rateConverter_.defineState(x, boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation()));
00517 }
00518 }
00519 else
00520 #endif
00521 {
00522 if ( global_number_resv_wells )
00523 {
00524 rateConverter_.defineState(x);
00525 }
00526 }
00527
00528 if (! resv_wells.empty()) {
00529 const PhaseUsage& pu = props_.phaseUsage();
00530 const std::vector<double>::size_type np = props_.numPhases();
00531
00532 std::vector<double> distr (np);
00533 std::vector<double> hrates(np);
00534 std::vector<double> prates(np);
00535
00536 for (std::vector<int>::const_iterator
00537 rp = resv_wells.begin(), e = resv_wells.end();
00538 rp != e; ++rp)
00539 {
00540 WellControls* ctrl = wells->ctrls[*rp];
00541 const bool is_producer = wells->type[*rp] == PRODUCER;
00542
00543
00544 {
00545 const int rctrl = SimFIBODetails::resv_control(ctrl);
00546
00547 if (0 <= rctrl) {
00548 const std::vector<double>::size_type off = (*rp) * np;
00549
00550 if (is_producer) {
00551
00552
00553 std::transform(xw.wellRates().begin() + (off + 0*np),
00554 xw.wellRates().begin() + (off + 1*np),
00555 prates.begin(), std::negate<double>());
00556 } else {
00557 std::copy(xw.wellRates().begin() + (off + 0*np),
00558 xw.wellRates().begin() + (off + 1*np),
00559 prates.begin());
00560 }
00561
00562 const int fipreg = 0;
00563 const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
00564 const int pvtreg = props_.cellPvtRegionIndex()[well_cell_top];
00565 rateConverter_.calcCoeff(fipreg, pvtreg, distr);
00566
00567 well_controls_iset_distr(ctrl, rctrl, & distr[0]);
00568 }
00569 }
00570
00571
00572
00573 if (is_producer && wells->name[*rp] != 0) {
00574 WellMap::const_iterator i = wmap.find(wells->name[*rp]);
00575
00576 if (i != wmap.end()) {
00577 const auto* wp = i->second;
00578
00579 const WellProductionProperties& p =
00580 wp->getProductionProperties(step);
00581
00582 if (! p.predictionMode) {
00583
00584 SimFIBODetails::historyRates(pu, p, hrates);
00585
00586 const int fipreg = 0;
00587 const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
00588 const int pvtreg = props_.cellPvtRegionIndex()[well_cell_top];
00589 rateConverter_.calcCoeff(fipreg, pvtreg, distr);
00590
00591
00592
00593
00594
00595 const double target =
00596 - std::inner_product(distr.begin(), distr.end(),
00597 hrates.begin(), 0.0);
00598
00599 well_controls_clear(ctrl);
00600 well_controls_assert_number_of_phases(ctrl, int(np));
00601
00602 static const double invalid_alq = -std::numeric_limits<double>::max();
00603 static const int invalid_vfp = -std::numeric_limits<int>::max();
00604
00605 const int ok_resv =
00606 well_controls_add_new(RESERVOIR_RATE, target,
00607 invalid_alq, invalid_vfp,
00608 & distr[0], ctrl);
00609
00610
00611
00612 double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
00613 const int ok_bhp =
00614 well_controls_add_new(BHP, bhp_limit,
00615 invalid_alq, invalid_vfp,
00616 NULL, ctrl);
00617
00618 if (ok_resv != 0 && ok_bhp != 0) {
00619 xw.currentControls()[*rp] = 0;
00620 well_controls_set_current(ctrl, 0);
00621 }
00622 }
00623 }
00624 }
00625 }
00626 }
00627
00628 if( wells )
00629 {
00630 for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
00631 WellControls* ctrl = wells->ctrls[w];
00632 const bool is_producer = wells->type[w] == PRODUCER;
00633 if (!is_producer && wells->name[w] != 0) {
00634 WellMap::const_iterator i = wmap.find(wells->name[w]);
00635 if (i != wmap.end()) {
00636 const auto* wp = i->second;
00637 const WellInjectionProperties& injector = wp->getInjectionProperties(step);
00638 if (!injector.predictionMode) {
00639
00640 static const double invalid_alq = -std::numeric_limits<double>::max();
00641 static const int invalid_vfp = -std::numeric_limits<int>::max();
00642
00643
00644 double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
00645 const int ok_bhp =
00646 well_controls_add_new(BHP, bhp_limit,
00647 invalid_alq, invalid_vfp,
00648 NULL, ctrl);
00649 if (!ok_bhp) {
00650 OPM_THROW(std::runtime_error, "Failed to add well control.");
00651 }
00652 }
00653 }
00654 }
00655 }
00656 }
00657 }
00658
00659 template <class Implementation>
00660 void
00661 SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units,
00662 std::vector<std::vector<double> >& fip)
00663 {
00664 for (size_t i = 0; i < fip.size(); ++i) {
00665 FIPUnitConvert(units, fip[i]);
00666 }
00667 }
00668
00669 template <class Implementation>
00670 void
00671 SimulatorBase<Implementation>::FIPUnitConvert(const UnitSystem& units, std::vector<double>& fip)
00672 {
00673 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
00674 fip[0] = unit::convert::to(fip[0], unit::stb);
00675 fip[1] = unit::convert::to(fip[1], unit::stb);
00676 fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
00677 fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
00678 fip[4] = unit::convert::to(fip[4], unit::stb);
00679 fip[5] = unit::convert::to(fip[5], unit::stb);
00680 fip[6] = unit::convert::to(fip[6], unit::psia);
00681 }
00682 else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
00683 fip[6] = unit::convert::to(fip[6], unit::barsa);
00684 }
00685 else {
00686 OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
00687 }
00688 }
00689
00690
00691 template <class Implementation>
00692 std::vector<double>
00693 SimulatorBase<Implementation>::FIPTotals(const std::vector<std::vector<double> >& fip, const ReservoirState& state)
00694 {
00695 std::vector<double> totals(7, 0.0);
00696 for (int i = 0; i < 5; ++i) {
00697 for (size_t reg = 0; reg < fip.size(); ++reg) {
00698 totals[i] += fip[reg][i];
00699 }
00700 }
00701 const int nc = Opm::AutoDiffGrid::numCells(grid_);
00702 const int np = state.numPhases();
00703 const PhaseUsage& pu = props_.phaseUsage();
00704 const DataBlock s = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
00705 std::vector<double> so(nc);
00706 std::vector<double> sg(nc);
00707 std::vector<double> hydrocarbon(nc);
00708
00709 const int oilpos = pu.phase_used[BlackoilPhases::Liquid] ? pu.phase_pos[BlackoilPhases::Liquid] : 0;
00710 const int gaspos = pu.phase_used[BlackoilPhases::Vapour] ? pu.phase_pos[BlackoilPhases::Vapour] : 0;
00711 const auto& soCol = s.col(oilpos);
00712 const auto& sgCol = s.col(gaspos);
00713 for (unsigned c = 0; c < so.size(); ++ c) {
00714 double mySo = 0.0;
00715 if (pu.phase_used[BlackoilPhases::Liquid]) {
00716 mySo = soCol[c];
00717 }
00718 double mySg = 0.0;
00719 if (pu.phase_used[BlackoilPhases::Vapour]) {
00720 mySg = sgCol[c];
00721 }
00722 so[c] = mySo;
00723 sg[c] = mySg;
00724 hydrocarbon[c] = mySo + mySg;
00725 }
00726 const std::vector<double> p = state.pressure();
00727 if ( ! is_parallel_run_ )
00728 {
00729 double tmp = 0.0;
00730 double tmp2 = 0.0;
00731 for (unsigned i = 0; i < p.size(); ++i) {
00732 tmp += p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
00733 tmp2 += geo_.poreVolume()[i] * hydrocarbon[i];
00734 }
00735 totals[5] = geo_.poreVolume().sum();
00736 totals[6] = tmp/tmp2;
00737 }
00738 else
00739 {
00740 #if HAVE_MPI
00741 const auto & pinfo =
00742 boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
00743 auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor<double>(),
00744 Opm::Reduction::makeGlobalSumFunctor<double>(),
00745 Opm::Reduction::makeGlobalSumFunctor<double>());
00746 std::vector<double> pav_nom(p.size());
00747 std::vector<double> pav_denom(pav_nom.size());
00748 for (unsigned i = 0; i < p.size(); ++i) {
00749 pav_nom[i] = p[i] * geo_.poreVolume()[i] * hydrocarbon[i];
00750 pav_denom[i] = geo_.poreVolume()[i] * hydrocarbon[i];
00751 }
00752
00753
00754 auto inputs = std::make_tuple(std::cref(geo_.poreVolume()),
00755 std::cref(pav_nom), std::cref(pav_denom));
00756 std::tuple<double, double, double> results(0.0, 0.0, 0.0);
00757
00758 pinfo.computeReduction(inputs, operators, results);
00759 using std::get;
00760 totals[5] = get<0>(results);
00761 totals[6] = get<1>(results)/get<2>(results);
00762
00763 #else
00764
00765 OPM_THROW(std::logic_error, "HAVE_MPI should be defined if we are running in parallel");
00766 #endif
00767 }
00768
00769 return totals;
00770 }
00771
00772
00773
00774 template <class Implementation>
00775 void
00776 SimulatorBase<Implementation>::outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
00777 {
00778 std::ostringstream ss;
00779 if (!reg) {
00780 ss << " ===================================================\n"
00781 << " : Field Totals :\n";
00782 } else {
00783 ss << " ===================================================\n"
00784 << " : FIPNUM report region "
00785 << std::setw(2) << reg << " :\n";
00786 }
00787 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
00788 ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
00789 << std::fixed << std::setprecision(0)
00790 << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
00791 if (!reg) {
00792 ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
00793 << " : Porv volumes are taken at reference conditions :\n";
00794 }
00795 ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
00796 }
00797 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
00798 ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
00799 << std::fixed << std::setprecision(0)
00800 << " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
00801 if (!reg) {
00802 ss << " : Pressure is weighted by hydrocarbon pore voulme :\n"
00803 << " : Pore volumes are taken at reference conditions :\n";
00804 }
00805 ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
00806 }
00807 ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
00808 << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
00809 << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
00810 << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
00811 << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
00812 << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
00813 << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
00814 << ":========================:==========================================:================:==========================================:\n";
00815 OpmLog::note(ss.str());
00816 }
00817
00818
00819 template <class Implementation>
00820 void
00821 SimulatorBase<Implementation>::
00822 updateListEconLimited(const std::unique_ptr<Solver>& solver,
00823 const Schedule& schedule,
00824 const int current_step,
00825 const Wells* wells,
00826 const WellState& well_state,
00827 DynamicListEconLimited& list_econ_limited) const
00828 {
00829
00830 solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
00831 well_state, list_econ_limited);
00832 }
00833
00834 template <class Implementation>
00835 void
00836 SimulatorBase<Implementation>::
00837 initHysteresisParams(ReservoirState& state)
00838 {
00839 typedef std::vector<double> VectorType;
00840
00841 const VectorType& somax = state.getCellData( "SOMAX" );
00842
00843 VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
00844 VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
00845
00846 VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
00847 VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
00848
00849 props_.setSatOilMax(somax);
00850 props_.setOilWaterHystParams(pcSwMdc_ow, krnSwMdc_ow, allcells_);
00851 props_.setGasOilHystParams(pcSwMdc_go, krnSwMdc_go, allcells_);
00852 }
00853
00854
00855 }