21 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
22 #define OPM_SIMULATORFULLYIMPLICITBLACKOILEBOS_HEADER_INCLUDED
24 #include <opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp>
25 #include <opm/autodiff/IterationReport.hpp>
26 #include <opm/autodiff/NonlinearSolver.hpp>
27 #include <opm/autodiff/BlackoilModelEbos.hpp>
28 #include <opm/autodiff/BlackoilModelParameters.hpp>
29 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
30 #include <opm/autodiff/BlackoilWellModel.hpp>
32 #include <opm/autodiff/SimFIBODetails.hpp>
33 #include <opm/autodiff/moduleVersion.hpp>
34 #include <opm/simulators/timestepping/AdaptiveTimeStepping.hpp>
35 #include <opm/core/utility/initHydroCarbonState.hpp>
36 #include <opm/core/utility/StopWatch.hpp>
38 #include <opm/common/Exceptions.hpp>
39 #include <opm/common/ErrorMacros.hpp>
41 #include <dune/common/unused.hh>
47 template<
class TypeTag>
51 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
52 typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
53 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
54 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
55 typedef typename GET_PROP_TYPE(TypeTag, Indices) BlackoilIndices;
56 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
57 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
58 typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector ;
59 typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
61 typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
64 typedef BlackoilState ReservoirState;
97 const ParameterGroup& param,
99 const bool has_disgas,
100 const bool has_vapoil,
101 const EclipseState& ,
103 const std::unordered_set<std::string>& defunct_well_names)
104 : ebosSimulator_(ebosSimulator),
107 solver_param_(param),
109 phaseUsage_(phaseUsageFromDeck(eclState())),
110 has_disgas_(has_disgas),
111 has_vapoil_(has_vapoil),
112 terminal_output_(param.getDefault(
"output_terminal", true)),
113 output_writer_(output_writer),
114 rateConverter_(createRateConverter_()),
115 defunct_well_names_( defunct_well_names ),
116 is_parallel_run_( false )
121 const ParallelISTLInformation& info =
124 terminal_output_ = terminal_output_ && ( info.communicator().rank() == 0 );
125 is_parallel_run_ = ( info.communicator().size() > 1 );
137 ReservoirState& state)
143 failureReport_ = SimulatorReport();
144 extractLegacyDepth_();
148 convertInput(0, state, ebosSimulator_ );
149 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
152 if (output_writer_.isRestart()) {
154 output_writer_.initFromRestartFile(phaseUsage_, grid(), state, prev_well_state, extra);
155 initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
156 initHysteresisParams(state);
158 convertInput(0, state, ebosSimulator_ );
159 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
165 ebosSimulator_.model().syncOverlap();
168 Opm::time::StopWatch solver_timer;
169 Opm::time::StopWatch step_timer;
170 Opm::time::StopWatch total_timer;
172 std::string tstep_filename = output_writer_.
outputDirectory() +
"/step_timing.txt";
173 std::ofstream tstep_os;
177 tstep_os.open(tstep_filename.c_str());
180 const auto& schedule = eclState().getSchedule();
183 const auto& events = schedule.getEvents();
184 std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping;
185 if( param_.getDefault(
"timestep.adaptive",
true ) )
188 if (param_.getDefault(
"use_TUNING",
false)) {
194 if (output_writer_.isRestart()) {
195 if (extra.suggested_step > 0.0) {
196 adaptiveTimeStepping->setSuggestedNextStep(extra.suggested_step);
201 std::string restorefilename = param_.getDefault(
"restorefile", std::string(
"") );
202 if( ! restorefilename.empty() )
205 const int desiredRestoreStep = param_.getDefault(
"restorestep",
int(-1) );
207 output_writer_.restore( timer,
211 desiredRestoreStep );
212 initHydroCarbonState(state, phaseUsage_, Opm::UgGridHelpers::numCells(grid()), has_disgas_, has_vapoil_);
213 initHysteresisParams(state);
215 convertInput(0, state, ebosSimulator_);
216 ebosSimulator_.model().invalidateIntensiveQuantitiesCache(0);
219 DynamicListEconLimited dynamic_list_econ_limited;
220 SimulatorReport report;
221 SimulatorReport stepReport;
223 std::vector<int> fipnum_global = eclState().get3DProperties().getIntGridProperty(
"FIPNUM").getData();
225 std::vector<int> fipnum(Opm::UgGridHelpers::numCells(grid()));
226 if (fipnum_global.empty()) {
227 std::fill(fipnum.begin(), fipnum.end(), 0);
229 for (
size_t c = 0; c < fipnum.size(); ++c) {
230 fipnum[c] = fipnum_global[Opm::UgGridHelpers::globalCell(grid())[c]];
233 std::vector<std::vector<double>> originalFluidInPlace;
234 std::vector<double> originalFluidInPlaceTotals;
237 while (!timer.
done()) {
240 if ( terminal_output_ )
242 std::ostringstream ss;
244 OpmLog::note(ss.str());
248 WellsManager wells_manager(eclState(),
250 Opm::UgGridHelpers::numCells(grid()),
251 Opm::UgGridHelpers::globalCell(grid()),
252 Opm::UgGridHelpers::cartDims(grid()),
253 Opm::UgGridHelpers::dimensions(grid()),
254 Opm::UgGridHelpers::cell2Faces(grid()),
255 Opm::UgGridHelpers::beginFaceCentroids(grid()),
256 dynamic_list_econ_limited,
258 defunct_well_names_ );
259 const Wells* wells = wells_manager.c_wells();
264 size_t nc = Opm::UgGridHelpers::numCells(grid());
265 std::vector<double> cellPressures(nc, 0.0);
266 const auto& gridView = ebosSimulator_.gridManager().gridView();
267 ElementContext elemCtx(ebosSimulator_);
268 const auto& elemEndIt = gridView.template end<0>();
269 for (
auto elemIt = gridView.template begin</*codim=*/0>();
273 const auto& elem = *elemIt;
274 if (elem.partitionType() != Dune::InteriorEntity) {
278 elemCtx.updatePrimaryStencil(elem);
279 elemCtx.updatePrimaryIntensiveQuantities(0);
281 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
282 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
283 const auto& fs = intQuants.fluidState();
285 const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
286 cellPressures[cellIdx] = p;
288 well_state.init(wells, cellPressures, prev_well_state, phaseUsage_);
291 handleAdditionalWellInflow(timer, wells_manager, well_state, wells);
297 solver_timer.start();
299 const auto& wells_ecl = eclState().getSchedule().getWells(timer.
currentStepNum());
300 extractLegacyCellPvtRegionIndex_();
301 WellModel well_model(wells, &(wells_manager.wellCollection()), wells_ecl, model_param_,
302 rateConverter_, terminal_output_, timer.
currentStepNum(), legacyCellPvtRegionIdx_);
306 for (
const auto& well : wells_ecl) {
318 auto solver = createSolver(well_model);
320 std::vector<std::vector<double>> currentFluidInPlace;
321 std::vector<double> currentFluidInPlaceTotals;
324 if (originalFluidInPlace.empty()) {
325 originalFluidInPlace = solver->computeFluidInPlace(fipnum);
326 originalFluidInPlaceTotals = FIPTotals(originalFluidInPlace);
327 FIPUnitConvert(eclState().getUnits(), originalFluidInPlace);
328 FIPUnitConvert(eclState().getUnits(), originalFluidInPlaceTotals);
330 currentFluidInPlace = originalFluidInPlace;
331 currentFluidInPlaceTotals = originalFluidInPlaceTotals;
336 Dune::Timer perfTimer;
339 if (terminal_output_) {
340 outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
341 for (
size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
342 outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
348 output_writer_.
writeTimeStep( timer, state, well_state, solver->model() );
350 report.output_write_time += perfTimer.stop();
353 if( terminal_output_ )
355 std::ostringstream step_msg;
356 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
357 step_msg.imbue(std::locale(std::locale::classic(), facet));
358 step_msg <<
"\nTime step " << std::setw(4) <<timer.
currentStepNum()
360 <<
"/" << (
double)unit::convert::to(timer.
totalTime(), unit::day)
362 OpmLog::info(step_msg.str());
365 solver->model().beginReportStep();
372 if( adaptiveTimeStepping ) {
373 bool event = events.hasEvent(ScheduleEvents::NEW_WELL, timer.
currentStepNum()) ||
374 events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE, timer.
currentStepNum()) ||
375 events.hasEvent(ScheduleEvents::INJECTION_UPDATE, timer.
currentStepNum()) ||
376 events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE, timer.
currentStepNum());
377 stepReport = adaptiveTimeStepping->step( timer, *solver, state, well_state, event, output_writer_,
378 output_writer_.requireFIPNUM() ? &fipnum : nullptr );
379 report += stepReport;
380 failureReport_ += adaptiveTimeStepping->failureReport();
384 stepReport = solver->step(timer, state, well_state);
385 report += stepReport;
386 failureReport_ += solver->failureReport();
388 if( terminal_output_ )
391 std::ostringstream iter_msg;
392 iter_msg <<
"Stepsize " << (double)unit::convert::to(timer.
currentStepLength(), unit::day);
393 if (solver->wellIterations() != 0) {
394 iter_msg <<
" days well iterations = " << solver->wellIterations() <<
", ";
396 iter_msg <<
"non-linear iterations = " << solver->nonlinearIterations()
397 <<
", total linear iterations = " << solver->linearIterations()
399 OpmLog::info(iter_msg.str());
403 solver->model().endReportStep();
409 report.solver_time += solver_timer.secsSinceStart();
413 stepReport.reportParam(tstep_os);
419 ReservoirState stateTrivial(0,0,0);
420 state = stateTrivial;
427 currentFluidInPlace = solver->computeFluidInPlace(fipnum);
428 currentFluidInPlaceTotals = FIPTotals(currentFluidInPlace);
432 FIPUnitConvert(eclState().getUnits(), currentFluidInPlace);
433 FIPUnitConvert(eclState().getUnits(), currentFluidInPlaceTotals);
435 if (terminal_output_ )
437 outputTimestampFIP(timer, version);
438 outputFluidInPlace(originalFluidInPlaceTotals, currentFluidInPlaceTotals,eclState().getUnits(), 0);
439 for (
size_t reg = 0; reg < originalFluidInPlace.size(); ++reg) {
440 outputFluidInPlace(originalFluidInPlace[reg], currentFluidInPlace[reg], eclState().getUnits(), reg+1);
445 "Time step took " + std::to_string(solver_timer.secsSinceStart()) +
" seconds; "
446 "total solver time " + std::to_string(report.solver_time) +
" seconds.";
451 Dune::Timer perfTimer;
453 const double nextstep = adaptiveTimeStepping ? adaptiveTimeStepping->suggestedNextStep() : -1.0;
454 output_writer_.
writeTimeStep( timer, state, well_state, solver->model(),
false, nextstep, report);
455 report.output_write_time += perfTimer.stop();
457 prev_well_state = well_state;
459 updateListEconLimited(solver, eclState().getSchedule(), timer.
currentStepNum(), wells,
460 well_state, dynamic_list_econ_limited);
465 report.total_time = total_timer.secsSinceStart();
466 report.converged =
true;
474 const Grid& grid()
const
475 {
return ebosSimulator_.gridManager().grid(); }
478 void handleAdditionalWellInflow(SimulatorTimer& ,
485 std::unique_ptr<Solver> createSolver(WellModel& well_model)
487 const auto& gridView = ebosSimulator_.gridView();
488 const PhaseUsage& phaseUsage = phaseUsage_;
489 const std::vector<bool> activePhases = detail::activePhases(phaseUsage);
490 const double gravity = ebosSimulator_.problem().gravity()[2];
495 int globalNumCells = gridView.size(0);
496 globalNumCells = gridView.comm().sum(globalNumCells);
498 well_model.init(phaseUsage,
504 auto model = std::unique_ptr<Model>(
new Model(ebosSimulator_,
511 return std::unique_ptr<Solver>(
new Solver(solver_param_, std::move(model)));
514 void computeRESV(
const std::size_t step,
518 typedef SimFIBODetails::WellMap WellMap;
520 const auto w_ecl = eclState().getSchedule().getWells(step);
521 const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
523 const std::vector<int>& resv_wells = SimFIBODetails::resvWells(wells, step, wmap);
525 const std::size_t number_resv_wells = resv_wells.size();
526 std::size_t global_number_resv_wells = number_resv_wells;
532 global_number_resv_wells = info.communicator().sum(global_number_resv_wells);
533 if ( global_number_resv_wells )
539 rateConverter_.template defineState<ElementContext>(ebosSimulator_);
545 if ( global_number_resv_wells )
547 rateConverter_.template defineState<ElementContext>(ebosSimulator_);
551 if (! resv_wells.empty()) {
552 const PhaseUsage& pu = phaseUsage_;
553 const std::vector<double>::size_type np = phaseUsage_.num_phases;
555 std::vector<double> distr (np);
556 std::vector<double> hrates(np);
557 std::vector<double> prates(np);
559 for (std::vector<int>::const_iterator
560 rp = resv_wells.begin(), e = resv_wells.end();
563 WellControls* ctrl = wells->ctrls[*rp];
564 const bool is_producer = wells->type[*rp] == PRODUCER;
565 const int well_cell_top = wells->well_cells[wells->well_connpos[*rp]];
566 const auto& eclProblem = ebosSimulator_.problem();
567 const int pvtreg = eclProblem.pvtRegionIndex(well_cell_top);
571 const int rctrl = SimFIBODetails::resv_control(ctrl);
574 const std::vector<double>::size_type off = (*rp) * np;
579 std::transform(xw.wellRates().begin() + (off + 0*np),
580 xw.wellRates().begin() + (off + 1*np),
581 prates.begin(), std::negate<double>());
583 std::copy(xw.wellRates().begin() + (off + 0*np),
584 xw.wellRates().begin() + (off + 1*np),
588 const int fipreg = 0;
589 rateConverter_.
calcCoeff(fipreg, pvtreg, distr);
591 well_controls_iset_distr(ctrl, rctrl, & distr[0]);
597 if (is_producer && wells->name[*rp] != 0) {
598 WellMap::const_iterator i = wmap.find(wells->name[*rp]);
600 if (i != wmap.end()) {
601 const auto* wp = i->second;
603 const WellProductionProperties& p =
604 wp->getProductionProperties(step);
606 if (! p.predictionMode) {
608 SimFIBODetails::historyRates(pu, p, hrates);
610 const int fipreg = 0;
611 rateConverter_.
calcCoeff(fipreg, pvtreg, distr);
617 const double target =
618 - std::inner_product(distr.begin(), distr.end(),
619 hrates.begin(), 0.0);
621 well_controls_clear(ctrl);
622 well_controls_assert_number_of_phases(ctrl,
int(np));
624 static const double invalid_alq = -std::numeric_limits<double>::max();
625 static const int invalid_vfp = -std::numeric_limits<int>::max();
628 well_controls_add_new(RESERVOIR_RATE, target,
629 invalid_alq, invalid_vfp,
634 double bhp_limit = (p.BHPLimit > 0) ? p.BHPLimit : unit::convert::from(1.0, unit::atm);
636 well_controls_add_new(BHP, bhp_limit,
637 invalid_alq, invalid_vfp,
640 if (ok_resv != 0 && ok_bhp != 0) {
641 xw.currentControls()[*rp] = 0;
642 well_controls_set_current(ctrl, 0);
652 for (
int w = 0, nw = wells->number_of_wells; w < nw; ++w) {
653 WellControls* ctrl = wells->ctrls[w];
654 const bool is_producer = wells->type[w] == PRODUCER;
655 if (!is_producer && wells->name[w] != 0) {
656 WellMap::const_iterator i = wmap.find(wells->name[w]);
657 if (i != wmap.end()) {
658 const auto* wp = i->second;
659 const WellInjectionProperties& injector = wp->getInjectionProperties(step);
660 if (!injector.predictionMode) {
662 static const double invalid_alq = -std::numeric_limits<double>::max();
663 static const int invalid_vfp = -std::numeric_limits<int>::max();
666 double bhp_limit = (injector.BHPLimit > 0) ? injector.BHPLimit : std::numeric_limits<double>::max();
668 well_controls_add_new(BHP, bhp_limit,
669 invalid_alq, invalid_vfp,
672 OPM_THROW(std::runtime_error,
"Failed to add well control.");
683 void updateListEconLimited(
const std::unique_ptr<Solver>& solver,
684 const Schedule& schedule,
685 const int current_step,
687 const WellState& well_state,
688 DynamicListEconLimited& list_econ_limited)
const
690 solver->model().wellModel().updateListEconLimited(schedule, current_step, wells,
691 well_state, list_econ_limited);
694 void FIPUnitConvert(
const UnitSystem& units,
695 std::vector<std::vector<double>>& fip)
697 for (
size_t i = 0; i < fip.size(); ++i) {
698 FIPUnitConvert(units, fip[i]);
703 void FIPUnitConvert(
const UnitSystem& units,
704 std::vector<double>& fip)
706 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
707 fip[0] = unit::convert::to(fip[0], unit::stb);
708 fip[1] = unit::convert::to(fip[1], unit::stb);
709 fip[2] = unit::convert::to(fip[2], 1000*unit::cubic(unit::feet));
710 fip[3] = unit::convert::to(fip[3], 1000*unit::cubic(unit::feet));
711 fip[4] = unit::convert::to(fip[4], unit::stb);
712 fip[5] = unit::convert::to(fip[5], unit::stb);
713 fip[6] = unit::convert::to(fip[6], unit::psia);
715 else if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
716 fip[6] = unit::convert::to(fip[6], unit::barsa);
719 OPM_THROW(std::runtime_error,
"Unsupported unit type for fluid in place output.");
724 std::vector<double> FIPTotals(
const std::vector<std::vector<double>>& fip)
726 std::vector<double> totals(7,0.0);
727 for (
int i = 0; i < 5; ++i) {
728 for (
size_t reg = 0; reg < fip.size(); ++reg) {
729 totals[i] += fip[reg][i];
733 const auto& gridView = ebosSimulator_.gridManager().gridView();
734 const auto& comm = gridView.comm();
735 double pv_hydrocarbon_sum = 0.0;
736 double p_pv_hydrocarbon_sum = 0.0;
738 ElementContext elemCtx(ebosSimulator_);
739 const auto& elemEndIt = gridView.template end<0>();
740 for (
auto elemIt = gridView.template begin</*codim=*/0>();
744 const auto& elem = *elemIt;
745 if (elem.partitionType() != Dune::InteriorEntity) {
749 elemCtx.updatePrimaryStencil(elem);
750 elemCtx.updatePrimaryIntensiveQuantities(0);
752 const unsigned cellIdx = elemCtx.globalSpaceIndex(0, 0);
753 const auto& intQuants = elemCtx.intensiveQuantities(0, 0);
754 const auto& fs = intQuants.fluidState();
756 const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
757 const double hydrocarbon = fs.saturation(FluidSystem::oilPhaseIdx).value() + fs.saturation(FluidSystem::gasPhaseIdx).value();
767 ebosSimulator_.model().dofTotalVolume(cellIdx)
768 * intQuants.porosity().value();
771 pv_hydrocarbon_sum += pv*hydrocarbon;
772 p_pv_hydrocarbon_sum += p*pv*hydrocarbon;
775 pv_hydrocarbon_sum = comm.sum(pv_hydrocarbon_sum);
776 p_pv_hydrocarbon_sum = comm.sum(p_pv_hydrocarbon_sum);
777 totals[5] = comm.sum(totals[5]);
778 totals[6] = (p_pv_hydrocarbon_sum / pv_hydrocarbon_sum);
784 void outputTimestampFIP(SimulatorTimer& timer,
const std::string version)
786 std::ostringstream ss;
787 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d %b %Y");
788 ss.imbue(std::locale(std::locale::classic(), facet));
789 ss <<
"\n **************************************************************************\n"
790 <<
" Balance at" << std::setw(10) << (double)unit::convert::to(timer.simulationTimeElapsed(), unit::day) <<
" Days"
791 <<
" *" << std::setw(30) << eclState().getTitle() <<
" *\n"
792 <<
" Report " << std::setw(4) << timer.reportStepNum() <<
" " << timer.currentDateTime()
793 <<
" * Flow version " << std::setw(11) << version <<
" *\n"
794 <<
" **************************************************************************\n";
795 OpmLog::note(ss.str());
799 void outputFluidInPlace(
const std::vector<double>& oip,
const std::vector<double>& cip,
const UnitSystem& units,
const int reg)
801 std::ostringstream ss;
803 ss <<
" ===================================================\n"
804 <<
" : Field Totals :\n";
806 ss <<
" ===================================================\n"
807 <<
" : FIPNUM report region "
808 << std::setw(2) << reg <<
" :\n";
810 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_METRIC) {
811 ss <<
" : PAV =" << std::setw(14) << cip[6] <<
" BARSA :\n"
812 << std::fixed << std::setprecision(0)
813 <<
" : PORV =" << std::setw(14) << cip[5] <<
" RM3 :\n";
815 ss <<
" : Pressure is weighted by hydrocarbon pore volume :\n"
816 <<
" : Porv volumes are taken at reference conditions :\n";
818 ss <<
" :--------------- Oil SM3 ---------------:-- Wat SM3 --:--------------- Gas SM3 ---------------:\n";
820 if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) {
821 ss <<
" : PAV =" << std::setw(14) << cip[6] <<
" PSIA :\n"
822 << std::fixed << std::setprecision(0)
823 <<
" : PORV =" << std::setw(14) << cip[5] <<
" RB :\n";
825 ss <<
" : Pressure is weighted by hydrocarbon pore volume :\n"
826 <<
" : Pore volumes are taken at reference conditions :\n";
828 ss <<
" :--------------- Oil STB ---------------:-- Wat STB --:--------------- Gas MSCF ---------------:\n";
830 ss <<
" : Liquid Vapour Total : Total : Free Dissolved Total :" <<
"\n"
831 <<
":------------------------:------------------------------------------:----------------:------------------------------------------:" <<
"\n"
832 <<
":Currently in place :" << std::setw(14) << cip[1] << std::setw(14) << cip[4] << std::setw(14) << (cip[1]+cip[4]) <<
":"
833 << std::setw(13) << cip[0] <<
" :" << std::setw(14) << (cip[2]) << std::setw(14) << cip[3] << std::setw(14) << (cip[2] + cip[3]) <<
":\n"
834 <<
":------------------------:------------------------------------------:----------------:------------------------------------------:\n"
835 <<
":Originally in place :" << std::setw(14) << oip[1] << std::setw(14) << oip[4] << std::setw(14) << (oip[1]+oip[4]) <<
":"
836 << std::setw(13) << oip[0] <<
" :" << std::setw(14) << oip[2] << std::setw(14) << oip[3] << std::setw(14) << (oip[2] + oip[3]) <<
":\n"
837 <<
":========================:==========================================:================:==========================================:\n";
838 OpmLog::note(ss.str());
842 const EclipseState& eclState()
const
843 {
return ebosSimulator_.gridManager().eclState(); }
845 void extractLegacyCellPvtRegionIndex_()
847 const auto& grid = ebosSimulator_.gridManager().grid();
848 const auto& eclProblem = ebosSimulator_.problem();
849 const unsigned numCells = grid.size(0);
851 legacyCellPvtRegionIdx_.resize(numCells);
852 for (
unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
853 legacyCellPvtRegionIdx_[cellIdx] =
854 eclProblem.pvtRegionIndex(cellIdx);
858 void initHysteresisParams(ReservoirState& state) {
859 const int num_cells = Opm::UgGridHelpers::numCells(grid());
861 typedef std::vector<double> VectorType;
863 const VectorType& somax = state.getCellData(
"SOMAX" );
865 for (
int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
866 ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
869 if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) {
870 auto matLawManager = ebosSimulator_.problem().materialLawManager();
872 VectorType& pcSwMdc_ow = state.getCellData(
"PCSWMDC_OW" );
873 VectorType& krnSwMdc_ow = state.getCellData(
"KRNSWMDC_OW" );
875 VectorType& pcSwMdc_go = state.getCellData(
"PCSWMDC_GO" );
876 VectorType& krnSwMdc_go = state.getCellData(
"KRNSWMDC_GO" );
878 for (
int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
879 matLawManager->setOilWaterHysteresisParams(
881 krnSwMdc_ow[cellIdx],
883 matLawManager->setGasOilHysteresisParams(
885 krnSwMdc_go[cellIdx],
891 void extractLegacyDepth_()
893 const auto& grid = ebosSimulator_.gridManager().grid();
894 const unsigned numCells = grid.size(0);
896 legacyDepth_.resize(numCells);
897 for (
unsigned cellIdx = 0; cellIdx < numCells; ++cellIdx) {
898 legacyDepth_[cellIdx] =
899 grid.cellCenterDepth(cellIdx);
904 void convertInput(
const int iterationIdx,
905 const ReservoirState& reservoirState,
906 Simulator& simulator )
const
908 SolutionVector& solution = simulator.model().solution( 0 );
909 const Opm::PhaseUsage pu = phaseUsage_;
911 const std::vector<bool> active = detail::activePhases(pu);
912 bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
913 bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);
915 const int numCells = reservoirState.numCells();
916 const int numPhases = phaseUsage_.num_phases;
917 const auto& oilPressure = reservoirState.pressure();
918 const auto& saturations = reservoirState.saturation();
919 const auto& rs = reservoirState.gasoilratio();
920 const auto& rv = reservoirState.rv();
921 for(
int cellIdx = 0; cellIdx<numCells; ++cellIdx )
924 PrimaryVariables& cellPv = solution[ cellIdx ];
926 if ( active[Water] ) {
927 cellPv[BlackoilIndices::waterSaturationIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Water]];
931 cellPv[BlackoilIndices::solventSaturationIdx] = reservoirState.getCellData( reservoirState.SSOL )[cellIdx];
935 cellPv[BlackoilIndices::polymerConcentrationIdx] = reservoirState.getCellData( reservoirState.POLYMER )[cellIdx];
941 if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::OilOnly && has_disgas_ )
943 cellPv[BlackoilIndices::compositionSwitchIdx] = rs[cellIdx];
944 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
945 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Rs );
947 else if( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasOnly && has_vapoil_ )
959 typedef Opm::SimpleModularFluidState<double,
970 false> SatOnlyFluidState;
971 SatOnlyFluidState fluidState;
972 if ( active[Water] ) {
973 fluidState.setSaturation(FluidSystem::waterPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Water]]);
976 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
978 fluidState.setSaturation(FluidSystem::oilPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Oil]]);
979 fluidState.setSaturation(FluidSystem::gasPhaseIdx, saturations[cellIdx*numPhases + pu.phase_pos[Gas]]);
981 double pC[3] = { 0.0, 0.0, 0.0 };
982 const MaterialLawParams& matParams = simulator.problem().materialLawParams(cellIdx);
983 MaterialLaw::capillaryPressures(pC, matParams, fluidState);
984 double pg = oilPressure[cellIdx] + (pC[FluidSystem::gasPhaseIdx] - pC[FluidSystem::oilPhaseIdx]);
986 cellPv[BlackoilIndices::compositionSwitchIdx] = rv[cellIdx];
987 cellPv[BlackoilIndices::pressureSwitchIdx] = pg;
988 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_pg_Rv );
992 assert( reservoirState.hydroCarbonState()[cellIdx] == HydroCarbonState::GasAndOil);
993 cellPv[BlackoilIndices::compositionSwitchIdx] = saturations[cellIdx*numPhases + pu.phase_pos[Gas]];
994 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[ cellIdx ];
995 cellPv.setPrimaryVarsMeaning( PrimaryVariables::Sw_po_Sg );
999 cellPv[BlackoilIndices::pressureSwitchIdx] = oilPressure[cellIdx];
1004 if( iterationIdx == 0 )
1006 simulator.model().solution( 1 ) = solution;
1010 RateConverterType createRateConverter_() {
1011 RateConverterType rate_converter(phaseUsage_,
1012 std::vector<int>(AutoDiffGrid::numCells(grid()), 0));
1013 return rate_converter;
1018 Simulator& ebosSimulator_;
1020 std::vector<int> legacyCellPvtRegionIdx_;
1021 std::vector<double> legacyDepth_;
1022 typedef typename Solver::SolverParameters SolverParameters;
1024 SimulatorReport failureReport_;
1026 const ParameterGroup param_;
1027 ModelParameters model_param_;
1028 SolverParameters solver_param_;
1031 NewtonIterationBlackoilInterface& solver_;
1032 PhaseUsage phaseUsage_;
1034 const bool has_disgas_;
1035 const bool has_vapoil_;
1036 bool terminal_output_;
1038 OutputWriter& output_writer_;
1039 RateConverterType rateConverter_;
1043 std::unordered_set<std::string> defunct_well_names_;
1046 bool is_parallel_run_;
1052 #endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
A nonlinear solver class suitable for general fully-implicit models, as well as pressure, transport and sequential models.
Definition: NonlinearSolver.hpp:37
int currentStepNum() const
Current step number.
Definition: SimulatorTimer.cpp:77
bool isIORank() const
Whether this process does write to disk.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:284
Convert component rates at surface conditions to phase (voidage) rates at reservoir conditions...
Definition: RateConverter.hpp:407
Class for handling the blackoil well model.
Definition: BlackoilWellModel.hpp:66
SimulatorFullyImplicitBlackoilEbos(Simulator &ebosSimulator, const ParameterGroup ¶m, NewtonIterationBlackoilInterface &linsolver, const bool has_disgas, const bool has_vapoil, const EclipseState &, OutputWriter &output_writer, const std::unordered_set< std::string > &defunct_well_names)
Initialise from parameters and objects to observe.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:96
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016...
Definition: moduleVersion.cpp:28
virtual const boost::any & parallelInformation() const =0
Get the information about the parallelization of the grid.
double totalTime() const
Total time.
Definition: SimulatorTimer.cpp:121
A model implementation for three-phase black oil.
Definition: BlackoilModelEbos.hpp:104
bool done() const
Return true if op++() has been called numSteps() times.
Definition: SimulatorTimer.cpp:155
Definition: AdaptiveTimeStepping.hpp:38
double simulationTimeElapsed() const
Time elapsed since the start of the simulation until the beginning of the current time step [s]...
Definition: SimulatorTimer.cpp:104
Facility for converting component rates at surface conditions to phase (voidage) rates at reservoir c...
Solver parameters for the BlackoilModel.
Definition: BlackoilModelParameters.hpp:31
const SimulatorReport & failureReport() const
Returns the simulator report for the failed substeps of the simulation.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:472
a simulator for the blackoil model
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:48
double currentStepLength() const
Current step length.
Definition: SimulatorTimer.cpp:91
boost::posix_time::ptime currentDateTime() const
Return current date.
Definition: SimulatorTimer.cpp:115
bool initialStep() const
Whether the current step is the first step.
Definition: SimulatorTimer.cpp:65
void report(std::ostream &os) const
Print a report with current and total time etc.
Definition: SimulatorTimer.cpp:136
void initWellStateMSWell(const Wells *wells, const std::vector< const Well * > &wells_ecl, const int time_step, const PhaseUsage &pu, const PrevWellState &prev_well_state)
init the MS well related.
Definition: WellStateFullyImplicitBlackoil.hpp:298
void writeTimeStep(const SimulatorTimerInterface &timer, const SimulationDataContainer &reservoirState, const Opm::WellStateFullyImplicitBlackoil &wellState, const Model &physicalModel, const bool substep=false, const double nextstep=-1.0, const SimulatorReport &simulatorReport=SimulatorReport())
Write a blackoil reservoir state to disk for later inspection with visualization tools like ResInsigh...
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:995
SimulatorReport run(SimulatorTimer &timer, ReservoirState &state)
Run the simulation.
Definition: SimulatorFullyImplicitBlackoilEbos.hpp:136
bool output() const
return true if output is enabled
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:281
The state of a set of wells, tailored for use by the fully implicit blackoil simulator.
Definition: WellStateFullyImplicitBlackoil.hpp:44
Wrapper class for VTK, Matlab, and ECL output.
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:206
const std::string & outputDirectory() const
return output directory
Definition: SimulatorFullyImplicitBlackoilOutput.hpp:278
Interface class for (linear) solvers for the fully implicit black-oil system.
Definition: NewtonIterationBlackoilInterface.hpp:31
void calcCoeff(const RegionId r, const int pvtRegionIdx, Coeff &coeff) const
Compute coefficients for surface-to-reservoir voidage conversion.
Definition: RateConverter.hpp:589
Definition: SimulatorTimer.hpp:34
bool use_multisegment_well_
Whether to use MultisegmentWell to handle multisegment wells it is something temporary before the mul...
Definition: BlackoilModelParameters.hpp:89