00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
00022 #define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
00023
00024 #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
00025 #include <opm/autodiff/IterationReport.hpp>
00026 #include <opm/autodiff/NonlinearSolver.hpp>
00027 #include <opm/autodiff/BlackoilModelEbos.hpp>
00028 #include <opm/autodiff/BlackoilModelParameters.hpp>
00029 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
00030 #include <opm/autodiff/BlackoilWellModel.hpp>
00031 #include <opm/autodiff/RateConverter.hpp>
00032 #include <opm/autodiff/SimFIBODetails.hpp>
00033 #include <opm/autodiff/moduleVersion.hpp>
00034 #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
00035 #include <opm/core/utility/initHydroCarbonState.hpp>
00036 #include <opm/core/utility/StopWatch.hpp>
00037
00038 #include <opm/common/Exceptions.hpp>
00039 #include <opm/common/ErrorMacros.hpp>
00040
00041 #include <dune/common/unused.hh>
00042
00043 namespace Opm {
00044
00045
00047 template<class TypeTag>
00048 class SimulatorFullyImplicitBlackoilEbos
00049 {
00050 public:
00051 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
00052 typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
00053 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
00054 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
00055 typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
00056 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
00057 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
00058 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
00059 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
00060
00061 typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
00062
00063 typedef WellStateFullyImplicitBlackoil WellState;
00064 typedef BlackoilState ReservoirState;
00065 typedef BlackoilOutputWriter OutputWriter;
00066 typedef BlackoilModelEbos<TypeTag> Model;
00067 typedef BlackoilModelParameters ModelParameters;
00068 typedef NonlinearSolver<Model> Solver;
00069 typedef BlackoilWellModel<TypeTag> WellModel;
00070 typedef RateConverter::SurfaceToReservoirVoidage<FluidSystem, std::vector<int> > RateConverterType;
00071
00072
00096 SimulatorFullyImplicitBlackoilEbos(Simulator& ebosSimulator,
00097 const ParameterGroup& param,
00098 NewtonIterationBlackoilInterface& linsolver,
00099 const bool has_disgas,
00100 const bool has_vapoil,
00101 const EclipseState& ,
00102 OutputWriter& output_writer,
00103 const std::unordered_set<std::string>& defunct_well_names)
00104 : ebosSimulator_(ebosSimulator),
00105 param_(param),
00106 model_param_(param),
00107 solver_param_(param),
00108 solver_(linsolver),
00109 phaseUsage_(phaseUsageFromDeck(eclState())),
00110 has_disgas_(has_disgas),
00111 has_vapoil_(has_vapoil),
00112 terminal_output_(param.getDefault("output_terminal", true)),
00113 output_writer_(output_writer),
00114 rateConverter_(createRateConverter_()),
00115 defunct_well_names_( defunct_well_names ),
00116 is_parallel_run_( false )
00117 {
00118 #if HAVE_MPI
00119 if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00120 {
00121 const ParallelISTLInformation& info =
00122 boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
00123
00124 terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
00125 is_parallel_run_ = ( info.communicator().size() > 1 );
00126 }
00127 #endif
00128 }
00129
00136 SimulatorReport run(SimulatorTimer& timer,
00137 ReservoirState& state)
00138 {
00139 WellState prev_well_state;
00140
00141 ExtraData extra;
00142
00143 failureReport_ = SimulatorReport();
00144 extractLegacyDepth_();
00145
00146
00147 if (timer.initialStep()) {
00148 convertInput(0, state, ebosSimulator_ );
00149 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
00150 }
00151
00152 if (output_writer_.isRestart()) {
00153
00154 output_writer_.initFromRestartFile(phaseUsage_, grid(), state, prev_well_state, extra);
00155 initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
00156 initHysteresisParams(state);
00157
00158 convertInput(0, state, ebosSimulator_ );
00159 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
00160 }
00161
00162
00163
00164
00165 ebosSimulator_.model().syncOverlap();
00166
00167
00168 Opm::time::StopWatch solver_timer;
00169 Opm::time::StopWatch step_timer;
00170 Opm::time::StopWatch total_timer;
00171 total_timer.start();
00172 std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt";
00173 std::ofstream tstep_os;
00174
00175 if ( output_writer_.output() && output_writer_.isIORank() )
00176 {
00177 tstep_os.open(tstep_filename.c_str());
00178 }
00179
00180 const auto& schedule = eclState().getSchedule();
00181
00182
00183 const auto& events = schedule.getEvents();
00184 std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
00185 if( param_.getDefault("timestep.adaptive", true ) )
00186 {
00187
00188 if (param_.getDefault("use_TUNING", false)) {
00189 adaptiveTimeStepping.reset( new AdaptiveTimeStepping( schedule.getTuning(), timer.currentStepNum(), param_, terminal_output_ ) );
00190 } else {
00191 adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_, terminal_output_ ) );
00192 }
00193
00194 if (output_writer_.isRestart()) {
00195 if (extra.suggested_step > 0.0) {
00196 adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
00197 }
00198 }
00199 }
00200
00201 std::string restorefilename = param_.getDefault("restorefile", std::string("") );
00202 if( ! restorefilename.empty() )
00203 {
00204
00205 const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) );
00206
00207 output_writer_.restore( timer,
00208 state,
00209 prev_well_state,
00210 restorefilename,
00211 desiredRestoreStep );
00212 initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
00213 initHysteresisParams(state);
00214
00215 convertInput(0, state, ebosSimulator_);
00216 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
00217 }
00218
00219 DynamicListEconLimited dynamic_list_econ_limited;
00220 SimulatorReport report;
00221 SimulatorReport stepReport;
00222
00223 std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty("FIPNUM").getData();
00224
00225 std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
00226 if (fipnum_global.empty()) {
00227 std::fill(fipnum.begin(), fipnum.end(), 0);
00228 } else {
00229 for (size_t c = 0; c < fipnum.size(); ++c) {
00230 fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
00231 }
00232 }
00233 std::vector<std::vector<double>> originalFluidInPlace;
00234 std::vector<double> originalFluidInPlaceTotals;
00235
00236
00237 while (!timer.done()) {
00238
00239 step_timer.start();
00240 if ( terminal_output_ )
00241 {
00242 std::ostringstream ss;
00243 timer.report(ss);
00244 OpmLog::note(ss.str());
00245 }
00246
00247
00248 WellsManager wells_manager(eclState(),
00249 timer.currentStepNum(),
00250 Opm::UgGridHelpers::numCells(grid()),
00251 Opm::UgGridHelpers::globalCell(grid()),
00252 Opm::UgGridHelpers::cartDims(grid()),
00253 Opm::UgGridHelpers::dimensions(grid()),
00254 Opm::UgGridHelpers::cell2Faces(grid()),
00255 Opm::UgGridHelpers::beginFaceCentroids(grid()),
00256 dynamic_list_econ_limited,
00257 is_parallel_run_,
00258 defunct_well_names_ );
00259 const Wells* wells = wells_manager.c_wells();
00260 WellState well_state;
00261
00262
00263
00264 size_t nc = Opm::UgGridHelpers::numCells(grid());
00265 std::vector<double> cellPressures(nc, 0.0);
00266 const auto& gridView = ebosSimulator_.gridManager().gridView();
00267 ElementContext elemCtx(ebosSimulator_);
00268 const auto& elemEndIt = gridView.template end<0>();
00269 for (auto elemIt = gridView.template begin</*codim=*/0>();
00270 elemIt != elemEndIt;
00271 ++elemIt)
00272 {
00273 const auto& elem = *elemIt;
00274 if (elem.partitionType() != Dune::InteriorEntity) {
00275 continue;
00276 }
00277
00278 elemCtx.updatePrimaryStencil(elem);
00279 elemCtx.updatePrimaryIntensiveQuantities(0);
00280
00281 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
00282 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
00283 const auto& fs = intQuants.fluidState();
00284
00285 const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
00286 cellPressures[cellIdx] = p;
00287 }
00288 well_state.init(wells, cellPressures, prev_well_state, phaseUsage_);
00289
00290
00291 handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
00292
00293
00294 computeRESV(timer.currentStepNum(), wells, well_state);
00295
00296
00297 solver_timer.start();
00298
00299 const auto& wells_ecl = eclState().getSchedule().getWells(timer.currentStepNum());
00300 extractLegacyCellPvtRegionIndex_();
00301 WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_,
00302 rateConverter_, terminal_output_, timer.currentStepNum(), legacyCellPvtRegionIdx_);
00303
00304
00305 if (model_param_.use_multisegment_well_) {
00306 for (const auto& well : wells_ecl) {
00307
00308
00309
00310 if (well->isMultiSegment(timer.currentStepNum()) ) {
00311 well_state.initWellStateMSWell(wells, wells_ecl, timer.currentStepNum(), phaseUsage_, prev_well_state);
00312 break;
00313 }
00314 }
00315 }
00316
00317
00318 auto solver = createSolver(well_model);
00319
00320 std::vector<std::vector<double>> currentFluidInPlace;
00321 std::vector<double> currentFluidInPlaceTotals;
00322
00323
00324 if (originalFluidInPlace.empty()) {
00325 originalFluidInPlace = solver->computeFluidInPlace(fipnum);
00326 originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace);
00327 FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
00328 FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
00329
00330 currentFluidInPlace = originalFluidInPlace;
00331 currentFluidInPlaceTotals = originalFluidInPlaceTotals;
00332 }
00333
00334
00335 if (timer.initialStep()) {
00336 Dune::Timer perfTimer;
00337 perfTimer.start();
00338
00339 if (terminal_output_) {
00340 outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
00341 for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
00342 outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
00343 }
00344 }
00345
00346
00347
00348 output_writer_.writeTimeStep( timer, state, well_state, solver->model() );
00349
00350 report.output_write_time += perfTimer.stop();
00351 }
00352
00353 if( terminal_output_ )
00354 {
00355 std::ostringstream step_msg;
00356 boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d-%b-%Y");
00357 step_msg.imbue(std::locale(std::locale::classic(), facet));
00358 step_msg << "\nTime step " << std::setw(4) <<timer.currentStepNum()
00359 << " at day " << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day)
00360 << "/" << (double)unit::convert::to(timer.totalTime(), unit::day)
00361 << ", date = " << timer.currentDateTime();
00362 OpmLog::info(step_msg.str());
00363 }
00364
00365 solver->model().beginReportStep();
00366
00367
00368
00369
00370
00371
00372 if( adaptiveTimeStepping ) {
00373 bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.currentStepNum()) ||
00374 events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.currentStepNum()) ||
00375 events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.currentStepNum()) ||
00376 events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.currentStepNum());
00377 stepReport = adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
00378 output_writer_.requireFIPNUM() ? &fipnum : nullptr );
00379 report += stepReport;
00380 failureReport_ += adaptiveTimeStepping->failureReport();
00381 }
00382 else {
00383
00384 stepReport = solver->step(timer, state, well_state);
00385 report += stepReport;
00386 failureReport_ += solver->failureReport();
00387
00388 if( terminal_output_ )
00389 {
00390
00391 std::ostringstream iter_msg;
00392 iter_msg << "Stepsize " << (double)unit::convert::to(timer.currentStepLength(), unit::day);
00393 if (solver->wellIterations() != 0) {
00394 iter_msg << " days well iterations = " << solver->wellIterations() << ", ";
00395 }
00396 iter_msg << "non-linear iterations = " << solver->nonlinearIterations()
00397 << ", total linear iterations = " << solver->linearIterations()
00398 << "\n";
00399 OpmLog::info(iter_msg.str());
00400 }
00401 }
00402
00403 solver->model().endReportStep();
00404
00405
00406 solver_timer.stop();
00407
00408
00409 report.solver_time += solver_timer.secsSinceStart();
00410
00411 if ( output_writer_.output() && output_writer_.isIORank() )
00412 {
00413 stepReport.reportParam(tstep_os);
00414 }
00415
00416
00417
00418 if (timer.initialStep()) {
00419 ReservoirState stateTrivial(0,0,0);
00420 state = stateTrivial;
00421 }
00422
00423
00424 ++timer;
00425
00426
00427 currentFluidInPlace = solver->computeFluidInPlace(fipnum);
00428 currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace);
00429
00430 const std::string version = moduleVersionName();
00431
00432 FIPUnitConvert(eclState().getUnits(), currentFluidInPlace);
00433 FIPUnitConvert(eclState().getUnits(), currentFluidInPlaceTotals);
00434
00435 if (terminal_output_ )
00436 {
00437 outputTimestampFIP(timer, version);
00438 outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
00439 for (size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
00440 outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
00441 }
00442
00443 std::string msg;
00444 msg =
00445 "Time step took " + std::to_string(solver_timer.secsSinceStart()) + " seconds; "
00446 "total solver time " + std::to_string(report.solver_time) + " seconds.";
00447 OpmLog::note(msg);
00448 }
00449
00450
00451 Dune::Timer perfTimer;
00452 perfTimer.start();
00453 const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
00454 output_writer_.writeTimeStep( timer, state, well_state, solver->model(), false, nextstep, report);
00455 report.output_write_time += perfTimer.stop();
00456
00457 prev_well_state = well_state;
00458
00459 updateListEconLimited(solver, eclState().getSchedule(), timer.currentStepNum(), wells,
00460 well_state, dynamic_list_econ_limited);
00461 }
00462
00463
00464 total_timer.stop();
00465 report.total_time = total_timer.secsSinceStart();
00466 report.converged = true;
00467 return report;
00468 }
00469
00472 const SimulatorReport& failureReport() const { return failureReport_; };
00473
00474 const Grid& grid() const
00475 { return ebosSimulator_.gridManager().grid(); }
00476
00477 protected:
00478 void handleAdditionalWellInflow(SimulatorTimer& ,
00479 WellsManager& ,
00480 WellState& ,
00481 const Wells* )
00482 {
00483 }
00484
00485 std::unique_ptr<Solver> createSolver(WellModel& well_model)
00486 {
00487 const auto& gridView = ebosSimulator_.gridView();
00488 const PhaseUsage& phaseUsage = phaseUsage_;
00489 const std::vector<bool> activePhases = detail::activePhases(phaseUsage);
00490 const double gravity = ebosSimulator_.problem().gravity()[2];
00491
00492
00493
00494
00495 int globalNumCells = gridView.size(0);
00496 globalNumCells = gridView.comm().sum(globalNumCells);
00497
00498 well_model.init(phaseUsage,
00499 activePhases,
00500 gravity,
00501 legacyDepth_,
00502 globalNumCells,
00503 grid());
00504 auto model = std::unique_ptr<Model>(new Model(ebosSimulator_,
00505 model_param_,
00506 well_model,
00507 rateConverter_,
00508 solver_,
00509 terminal_output_));
00510
00511 return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
00512 }
00513
00514 void computeRESV(const std::size_t step,
00515 const Wells* wells,
00516 WellState& xw)
00517 {
00518 typedef SimFIBODetails::WellMap WellMap;
00519
00520 const auto w_ecl = eclState().getSchedule().getWells(step);
00521 const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
00522
00523 const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
00524
00525 const std::size_t number_resv_wells = resv_wells.size();
00526 std::size_t global_number_resv_wells = number_resv_wells;
00527 #if HAVE_MPI
00528 if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) )
00529 {
00530 const auto& info =
00531 boost::any_cast<const ParallelISTLInformation&>(solver_.parallelInformation());
00532 global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
00533 if ( global_number_resv_wells )
00534 {
00535
00536
00537
00538
00539 rateConverter_.template defineState<ElementContext>(ebosSimulator_);
00540 }
00541 }
00542 else
00543 #endif
00544 {
00545 if ( global_number_resv_wells )
00546 {
00547 rateConverter_.template defineState<ElementContext>(ebosSimulator_);
00548 }
00549 }
00550
00551 if (! resv_wells.empty()) {
00552 const PhaseUsage& pu = phaseUsage_;
00553 const std::vector<double>::size_type np = phaseUsage_.num_phases;
00554
00555 std::vector<double> distr (np);
00556 std::vector<double> hrates(np);
00557 std::vector<double> prates(np);
00558
00559 for (std::vector<int>::const_iterator
00560 rp = resv_wells.begin(), e = resv_wells.end();
00561 rp != e; ++rp)
00562 {
00563 WellControls* ctrl = wells->ctrls[*rp];
00564 const bool is_producer = wells->type[*rp] == PRODUCER;
00565 const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
00566 const auto& eclProblem = ebosSimulator_.problem();
00567 const int pvtreg = eclProblem.pvtRegionIndex(well_cell_top);
00568
00569
00570 {
00571 const int rctrl = SimFIBODetails::resv_control(ctrl);
00572
00573 if (0 <= rctrl) {
00574 const std::vector<double>::size_type off = (*rp) * np;
00575
00576 if (is_producer) {
00577
00578
00579 std::transform(xw.wellRates().begin() + (off + 0*np),
00580 xw.wellRates().begin() + (off + 1*np),
00581 prates.begin(), std::negate<double>());
00582 } else {
00583 std::copy(xw.wellRates().begin() + (off + 0*np),
00584 xw.wellRates().begin() + (off + 1*np),
00585 prates.begin());
00586 }
00587
00588 const int fipreg = 0;
00589 rateConverter_.calcCoeff(fipreg, pvtreg, distr);
00590
00591 well_controls_iset_distr(ctrl, rctrl, & distr[0]);
00592 }
00593 }
00594
00595
00596
00597 if (is_producer && wells->name[*rp] != 0) {
00598 WellMap::const_iterator i = wmap.find(wells->name[*rp]);
00599
00600 if (i != wmap.end()) {
00601 const auto* wp = i->second;
00602
00603 const WellProductionProperties& p =
00604 wp->getProductionProperties(step);
00605
00606 if (! p.predictionMode) {
00607
00608 SimFIBODetails::historyRates(pu, p, hrates);
00609
00610 const int fipreg = 0;
00611 rateConverter_.calcCoeff(fipreg, pvtreg, distr);
00612
00613
00614
00615
00616
00617 const double target =
00618 - std::inner_product(distr.begin(), distr.end(),
00619 hrates.begin(), 0.0);
00620
00621 well_controls_clear(ctrl);
00622 well_controls_assert_number_of_phases(ctrl, int(np));
00623
00624 static const double invalid_alq = -std::numeric_limits<double>::max();
00625 static const int invalid_vfp = -std::numeric_limits<int>::max();
00626
00627 const int ok_resv =
00628 well_controls_add_new(RESERVOIR_RATE, target,
00629 invalid_alq, invalid_vfp,
00630 & distr[0], ctrl);
00631
00632
00633
00634 double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
00635 const int ok_bhp =
00636 well_controls_add_new(BHP, bhp_limit,
00637 invalid_alq, invalid_vfp,
00638 NULL, ctrl);
00639
00640 if (ok_resv != 0 && ok_bhp != 0) {
00641 xw.currentControls()[*rp] = 0;
00642 well_controls_set_current(ctrl, 0);
00643 }
00644 }
00645 }
00646 }
00647 }
00648 }
00649
00650 if( wells )
00651 {
00652 for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
00653 WellControls* ctrl = wells->ctrls[w];
00654 const bool is_producer = wells->type[w] == PRODUCER;
00655 if (!is_producer && wells->name[w] != 0) {
00656 WellMap::const_iterator i = wmap.find(wells->name[w]);
00657 if (i != wmap.end()) {
00658 const auto* wp = i->second;
00659 const WellInjectionProperties& injector = wp->getInjectionProperties(step);
00660 if (!injector.predictionMode) {
00661
00662 static const double invalid_alq = -std::numeric_limits<double>::max();
00663 static const int invalid_vfp = -std::numeric_limits<int>::max();
00664
00665
00666 double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
00667 const int ok_bhp =
00668 well_controls_add_new(BHP, bhp_limit,
00669 invalid_alq, invalid_vfp,
00670 NULL, ctrl);
00671 if (!ok_bhp) {
00672 OPM_THROW(std::runtime_error, "Failed to add well control.");
00673 }
00674 }
00675 }
00676 }
00677 }
00678 }
00679 }
00680
00681
00682
00683 void updateListEconLimited(const std::unique_ptr<Solver>& solver,
00684 const Schedule& schedule,
00685 const int current_step,
00686 const Wells* wells,
00687 const WellState& well_state,
00688 DynamicListEconLimited& list_econ_limited) const
00689 {
00690 solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
00691 well_state, list_econ_limited);
00692 }
00693
00694 void FIPUnitConvert(const UnitSystem& units,
00695 std::vector<std::vector<double>>& fip)
00696 {
00697 for (size_t i = 0; i < fip.size(); ++i) {
00698 FIPUnitConvert(units, fip[i]);
00699 }
00700 }
00701
00702
00703 void FIPUnitConvert(const UnitSystem& units,
00704 std::vector<double>& fip)
00705 {
00706 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
00707 fip[0] = unit::convert::to(fip[0], unit::stb);
00708 fip[1] = unit::convert::to(fip[1], unit::stb);
00709 fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
00710 fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
00711 fip[4] = unit::convert::to(fip[4], unit::stb);
00712 fip[5] = unit::convert::to(fip[5], unit::stb);
00713 fip[6] = unit::convert::to(fip[6], unit::psia);
00714 }
00715 else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
00716 fip[6] = unit::convert::to(fip[6], unit::barsa);
00717 }
00718 else {
00719 OPM_THROW(std::runtime_error, "Unsupported unit type for fluid in place output.");
00720 }
00721 }
00722
00723
00724 std::vector<double> FIPTotals(const std::vector<std::vector<double>>& fip)
00725 {
00726 std::vector<double> totals(7,0.0);
00727 for (int i = 0; i < 5; ++i) {
00728 for (size_t reg = 0; reg < fip.size(); ++reg) {
00729 totals[i] += fip[reg][i];
00730 }
00731 }
00732
00733 const auto& gridView = ebosSimulator_.gridManager().gridView();
00734 const auto& comm = gridView.comm();
00735 double pv_hydrocarbon_sum = 0.0;
00736 double p_pv_hydrocarbon_sum = 0.0;
00737
00738 ElementContext elemCtx(ebosSimulator_);
00739 const auto& elemEndIt = gridView.template end<0>();
00740 for (auto elemIt = gridView.template begin</*codim=*/0>();
00741 elemIt != elemEndIt;
00742 ++elemIt)
00743 {
00744 const auto& elem = *elemIt;
00745 if (elem.partitionType() != Dune::InteriorEntity) {
00746 continue;
00747 }
00748
00749 elemCtx.updatePrimaryStencil(elem);
00750 elemCtx.updatePrimaryIntensiveQuantities(0);
00751
00752 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
00753 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
00754 const auto& fs = intQuants.fluidState();
00755
00756 const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
00757 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
00758
00759
00760
00761
00762
00763
00764
00765
00766 const double pv =
00767 ebosSimulator_.model().dofTotalVolume(cellIdx)
00768 * intQuants.porosity().value();
00769
00770 totals[5] += pv;
00771 pv_hydrocarbon_sum += pv*hydrocarbon;
00772 p_pv_hydrocarbon_sum += p*pv*hydrocarbon;
00773 }
00774
00775 pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum);
00776 p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum);
00777 totals[5] = comm.sum(totals[5]);
00778 totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum);
00779
00780 return totals;
00781 }
00782
00783
00784 void outputTimestampFIP(SimulatorTimer& timer, const std::string version)
00785 {
00786 std::ostringstream ss;
00787 boost::posix_time::time_facet* facet = new boost::posix_time::time_facet("%d %b %Y");
00788 ss.imbue(std::locale(std::locale::classic(), facet));
00789 ss << "\n **************************************************************************\n"
00790 << " Balance at" << std::setw(10) << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day) << " Days"
00791 << " *" << std::setw(30) << eclState().getTitle() << " *\n"
00792 << " Report " << std::setw(4) << timer.reportStepNum() << " " << timer.currentDateTime()
00793 << " * Flow version " << std::setw(11) << version << " *\n"
00794 << " **************************************************************************\n";
00795 OpmLog::note(ss.str());
00796 }
00797
00798
00799 void outputFluidInPlace(const std::vector<double>& oip, const std::vector<double>& cip, const UnitSystem& units, const int reg)
00800 {
00801 std::ostringstream ss;
00802 if (!reg) {
00803 ss << " ===================================================\n"
00804 << " : Field Totals :\n";
00805 } else {
00806 ss << " ===================================================\n"
00807 << " : FIPNUM report region "
00808 << std::setw(2) << reg << " :\n";
00809 }
00810 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
00811 ss << " : PAV =" << std::setw(14) << cip[6] << " BARSA :\n"
00812 << std::fixed << std::setprecision(0)
00813 << " : PORV =" << std::setw(14) << cip[5] << " RM3 :\n";
00814 if (!reg) {
00815 ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
00816 << " : Porv volumes are taken at reference conditions :\n";
00817 }
00818 ss << " :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
00819 }
00820 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
00821 ss << " : PAV =" << std::setw(14) << cip[6] << " PSIA :\n"
00822 << std::fixed << std::setprecision(0)
00823 << " : PORV =" << std::setw(14) << cip[5] << " RB :\n";
00824 if (!reg) {
00825 ss << " : Pressure is weighted by hydrocarbon pore volume :\n"
00826 << " : Pore volumes are taken at reference conditions :\n";
00827 }
00828 ss << " :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
00829 }
00830 ss << " : Liquid Vapour Total : Total : Free Dissolved Total :" << "\n"
00831 << ":------------------------:------------------------------------------:----------------:------------------------------------------:" << "\n"
00832 << ":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) << ":"
00833 << std::setw(13) << cip[0] << " :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) << ":\n"
00834 << ":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
00835 << ":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) << ":"
00836 << std::setw(13) << oip[0] << " :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) << ":\n"
00837 << ":========================:==========================================:================:==========================================:\n";
00838 OpmLog::note(ss.str());
00839 }
00840
00841
00842 const EclipseState& eclState() const
00843 { return ebosSimulator_.gridManager().eclState(); }
00844
00845 void extractLegacyCellPvtRegionIndex_()
00846 {
00847 const auto& grid = ebosSimulator_.gridManager().grid();
00848 const auto& eclProblem = ebosSimulator_.problem();
00849 const unsigned numCells = grid.size(0);
00850
00851 legacyCellPvtRegionIdx_.resize(numCells);
00852 for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00853 legacyCellPvtRegionIdx_[cellIdx] =
00854 eclProblem.pvtRegionIndex(cellIdx);
00855 }
00856 }
00857
00858 void initHysteresisParams(ReservoirState& state) {
00859 const int num_cells = Opm::UgGridHelpers::numCells(grid());
00860
00861 typedef std::vector<double> VectorType;
00862
00863 const VectorType& somax = state.getCellData( "SOMAX" );
00864
00865 for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
00866 ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
00867 }
00868
00869 if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) {
00870 auto matLawManager = ebosSimulator_.problem().materialLawManager();
00871
00872 VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
00873 VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
00874
00875 VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
00876 VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
00877
00878 for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
00879 matLawManager->setOilWaterHysteresisParams(
00880 pcSwMdc_ow[cellIdx],
00881 krnSwMdc_ow[cellIdx],
00882 cellIdx);
00883 matLawManager->setGasOilHysteresisParams(
00884 pcSwMdc_go[cellIdx],
00885 krnSwMdc_go[cellIdx],
00886 cellIdx);
00887 }
00888 }
00889 }
00890
00891 void extractLegacyDepth_()
00892 {
00893 const auto& grid = ebosSimulator_.gridManager().grid();
00894 const unsigned numCells = grid.size(0);
00895
00896 legacyDepth_.resize(numCells);
00897 for (unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00898 legacyDepth_[cellIdx] =
00899 grid.cellCenterDepth(cellIdx);
00900 }
00901 }
00902
00903
00904 void convertInput( const int iterationIdx,
00905 const ReservoirState& reservoirState,
00906 Simulator& simulator ) const
00907 {
00908 SolutionVector& solution = simulator.model().solution( 0 );
00909 const Opm::PhaseUsage pu = phaseUsage_;
00910
00911 const std::vector<bool> active = detail::activePhases(pu);
00912 bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
00913 bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
00914
00915 const int numCells = reservoirState.numCells();
00916 const int numPhases = phaseUsage_.num_phases;
00917 const auto& oilPressure = reservoirState.pressure();
00918 const auto& saturations = reservoirState.saturation();
00919 const auto& rs = reservoirState.gasoilratio();
00920 const auto& rv = reservoirState.rv();
00921 for( int cellIdx = 0; cellIdx<numCells; ++cellIdx )
00922 {
00923
00924 PrimaryVariables& cellPv = solution[ cellIdx ];
00925
00926 if ( active[Water] ) {
00927 cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
00928 }
00929
00930 if (has_solvent) {
00931 cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
00932 }
00933
00934 if (has_polymer) {
00935 cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
00936 }
00937
00938
00939
00940 if ( active[Gas] ) {
00941 if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
00942 {
00943 cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
00944 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
00945 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
00946 }
00947 else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
00948 {
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959 typedef Opm::SimpleModularFluidState<double,
00960 3,
00961 3,
00962 FluidSystem,
00963 false,
00964 false,
00965 false,
00966 false,
00967 true,
00968 false,
00969 false,
00970 false> SatOnlyFluidState;
00971 SatOnlyFluidState fluidState;
00972 if ( active[Water] ) {
00973 fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
00974 }
00975 else {
00976 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
00977 }
00978 fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
00979 fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
00980
00981 double pC[3] = { 0.0, 0.0, 0.0 };
00982 const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
00983 MaterialLaw::capillaryPressures(pC, matParams, fluidState);
00984 double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
00985
00986 cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
00987 cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
00988 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
00989 }
00990 else
00991 {
00992 assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
00993 cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
00994 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
00995 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
00996 }
00997 } else {
00998
00999 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
01000 }
01001 }
01002
01003
01004 if( iterationIdx == 0 )
01005 {
01006 simulator.model().solution( 1 ) = solution;
01007 }
01008 }
01009
01010 RateConverterType createRateConverter_() {
01011 RateConverterType rate_converter(phaseUsage_,
01012 std::vector<int>(AutoDiffGrid::numCells(grid()), 0));
01013 return rate_converter;
01014 }
01015
01016
01017
01018 Simulator& ebosSimulator_;
01019
01020 std::vector<int> legacyCellPvtRegionIdx_;
01021 std::vector<double> legacyDepth_;
01022 typedef typename Solver::SolverParameters SolverParameters;
01023
01024 SimulatorReport failureReport_;
01025
01026 const ParameterGroup param_;
01027 ModelParameters model_param_;
01028 SolverParameters solver_param_;
01029
01030
01031 NewtonIterationBlackoilInterface& solver_;
01032 PhaseUsage phaseUsage_;
01033
01034 const bool has_disgas_;
01035 const bool has_vapoil_;
01036 bool terminal_output_;
01037
01038 OutputWriter& output_writer_;
01039 RateConverterType rateConverter_;
01040
01041
01042
01043 std::unordered_set<std::string> defunct_well_names_;
01044
01045
01046 bool is_parallel_run_;
01047
01048 };
01049
01050 }
01051
01052 #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED