diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index 914c2ff88..2208636bf 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -58,8 +58,13 @@ #include #include -#include +#include +#include +#include #include +#include +#include +#include #include #include @@ -112,9 +117,48 @@ try double gravity[3] = { 0.0 }; std::string deck_filename = param.get("deck_filename"); - Opm::ParserPtr newParser(new Opm::Parser()); - Opm::DeckConstPtr deck = newParser->parseFile(deck_filename); - std::shared_ptr eclipseState(new EclipseState(deck)); + // Write parameters used for later reference. + bool output = param.getDefault("output", true); + std::string output_dir; + if (output) { + // Create output directory if needed. + output_dir = + param.getDefault("output_dir", std::string("output")); + boost::filesystem::path fpath(output_dir); + try { + create_directories(fpath); + } + catch (...) { + std::cerr << "Creating directories failed: " << fpath << std::endl; + return EXIT_FAILURE; + } + // Write simulation parameters. + param.writeParam(output_dir + "/simulation.param"); + } + + std::string logFile = output_dir + "/LOGFILE.txt"; + Opm::ParserPtr parser(new Opm::Parser()); + { + std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); + std::shared_ptr counterLog = std::make_shared(Opm::Log::DefaultMessageTypes); + + Opm::OpmLog::addBackend( "STREAM" , streamLog ); + Opm::OpmLog::addBackend( "COUNTER" , counterLog ); + } + + Opm::DeckConstPtr deck; + std::shared_ptr eclipseState; + try { + deck = parser->parseFile(deck_filename); + Opm::checkDeck(deck); + eclipseState.reset(new Opm::EclipseState(deck)); + } + catch (const std::invalid_argument& e) { + std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl; + std::cerr << "Exception caught: " << e.what() << std::endl; + return EXIT_FAILURE; + } + // Grid init std::vector porv; if (eclipseState->hasDoubleGridProperty("PORV")) { @@ -123,20 +167,23 @@ try grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); auto &cGrid = *grid->c_grid(); const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); - Opm::EclipseWriter outputWriter(param, - eclipseState, - pu, - cGrid.number_of_cells, - cGrid.global_cell); + Opm::BlackoilOutputWriter outputWriter(cGrid, + param, + eclipseState, + pu ); + // Rock and fluid init props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid())); PolymerProperties polymer_props(deck, eclipseState); PolymerPropsAd polymer_props_ad(polymer_props); + // Rock compressibility. rock_comp.reset(new RockCompressibility(deck, eclipseState)); + // Gravity. gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; + // Init state variables (saturation and pressure). if (param.has("init_saturation")) { initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); @@ -155,22 +202,6 @@ try fis_solver.reset(new NewtonIterationBlackoilSimple(param)); } - // Write parameters used for later reference. - bool output = param.getDefault("output", true); - std::string output_dir; - if (output) { - output_dir = - param.getDefault("output_dir", std::string("output")); - boost::filesystem::path fpath(output_dir); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - param.writeParam(output_dir + "/simulation.param"); - } - Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap()); SimulatorTimer simtimer; simtimer.init(timeMap); @@ -193,17 +224,18 @@ try SimulatorReport fullReport; // Create and run simulator. Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); - SimulatorFullyImplicitCompressiblePolymer simulator(param, - *grid->c_grid(), - geology, - *new_props, - polymer_props_ad, - rock_comp->isActive() ? rock_comp.get() : 0, - eclipseState, - outputWriter, - deck, - *fis_solver, - grav); + SimulatorFullyImplicitCompressiblePolymer + simulator(param, + *grid->c_grid(), + geology, + *new_props, + polymer_props_ad, + rock_comp->isActive() ? rock_comp.get() : 0, + eclipseState, + outputWriter, + deck, + *fis_solver, + grav); fullReport= simulator.run(simtimer, state); std::cout << "\n\n================ End of simulation ===============\n\n"; diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 78caf9c26..a46e51339 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -21,63 +21,75 @@ #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED -#include -#include +#include +#include +#include +#include +#include +#include -struct UnstructuredGrid; -struct Wells; -struct FlowBoundaryConditions; +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +//#include +#include +#include +#include +#include + +#include + +//#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace Opm { - namespace parameter { class ParameterGroup; } - class BlackoilPropsAdInterface; - class RockCompressibility; - class DerivedGeology; - class NewtonIterationBlackoilInterface; - class SimulatorTimer; - class PolymerBlackoilState; - class WellStateFullyImplicitBlackoil; - class EclipseState; - class BlackoilOutputWriter; - class PolymerPropsAd; - class PolymerInflowInterface; - struct SimulatorReport; - /// Class collecting all necessary components for a two-phase simulation. - template + template class SimulatorFullyImplicitBlackoilPolymer + : public SimulatorBase > { + typedef SimulatorFullyImplicitBlackoilPolymer ThisType; + typedef SimulatorBase BaseType; + public: /// \brief The type of the grid that we use. - typedef T Grid; - /// Initialise from parameters and objects to observe. - /// \param[in] param parameters, this class accepts the following: - /// parameter (default) effect - /// ----------------------------------------------------------- - /// output (true) write output to files? - /// output_dir ("output") output directoty - /// output_interval (1) output every nth step - /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) - /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) - /// nl_pressure_maxiter (10) max nonlinear iterations in pressure - /// nl_maxiter (30) max nonlinear iterations in transport - /// nl_tolerance (1e-9) transport solver absolute residual tolerance - /// num_transport_substeps (1) number of transport steps per pressure step - /// use_segregation_split (false) solve for gravity segregation (if false, - /// segregation is ignored). - /// - /// \param[in] grid grid data structure - /// \param[in] geo derived geological properties - /// \param[in] props fluid and rock properties - /// \param[in] rock_comp_props if non-null, rock compressibility properties - /// \param[in] linsolver linear solver - /// \param[in] gravity if non-null, gravity vector - /// \param[in] disgas true for dissolved gas option - /// \param[in] vapoil true for vaporized oil option - /// \param[in] eclipse_state - /// \param[in] output_writer - /// \param[in] threshold_pressures_by_face if nonempty, threshold pressures that inhibit flow + typedef GridT Grid; + SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, const Grid& grid, const DerivedGeology& geo, @@ -94,20 +106,13 @@ namespace Opm Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face); - /// Run the simulation. - /// This will run succesive timesteps until timer.done() is true. It will - /// modify the reservoir and well states. - /// \param[in,out] timer governs the requested reporting timesteps - /// \param[in,out] state state of reservoir: pressure, fluxes - /// \param[in,out] well_state state of wells: bhp, perforation rates - /// \return simulation report, with timing data SimulatorReport run(SimulatorTimer& timer, PolymerBlackoilState& state); private: - class Impl; - // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. - std::shared_ptr pimpl_; + const PolymerPropsAd& polymer_props_; + bool has_polymer_; + DeckConstPtr deck_; }; } // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 0e7abd142..c6801adad 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -19,211 +19,55 @@ along with OPM. If not, see . */ -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -//#include -#include -#include -#include -#include - -#include - -//#include -#include - -#include -#include -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - namespace Opm { template - class SimulatorFullyImplicitBlackoilPolymer::Impl - { - public: - Impl(const parameter::ParameterGroup& param, - const Grid& grid, - const DerivedGeology& geo, - BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity, - bool has_disgas, - bool has_vapoil, - bool has_polymer, - std::shared_ptr eclipse_state, - BlackoilOutputWriter& output_writer, - Opm::DeckConstPtr& deck, - const std::vector& threshold_pressures_by_face); - - SimulatorReport run(SimulatorTimer& timer, - PolymerBlackoilState& state); - - private: - // Data. - typedef RateConverter:: - SurfaceToReservoirVoidage< BlackoilPropsAdInterface, - std::vector > RateConverterType; - - const parameter::ParameterGroup param_; - - // Observed objects. - const Grid& grid_; - BlackoilPropsAdInterface& props_; - const PolymerPropsAd& polymer_props_; - const RockCompressibility* rock_comp_props_; - const double* gravity_; - // Solvers - const DerivedGeology& geo_; - NewtonIterationBlackoilInterface& solver_; - // Misc. data - std::vector allcells_; - const bool has_disgas_; - const bool has_vapoil_; - const bool has_polymer_; - bool terminal_output_; - // eclipse_state - std::shared_ptr eclipse_state_; - // output_writer - BlackoilOutputWriter& output_writer_; - Opm::DeckConstPtr& deck_; - RateConverterType rateConverter_; - // Threshold pressures. - std::vector threshold_pressures_by_face_; - - void - computeRESV(const std::size_t step, - const Wells* wells, - const BlackoilState& x, - WellStateFullyImplicitBlackoilPolymer& xw); - }; - - - - - template - SimulatorFullyImplicitBlackoilPolymer::SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, - const Grid& grid, - const DerivedGeology& geo, - BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity, - const bool has_disgas, - const bool has_vapoil, - const bool has_polymer, - std::shared_ptr eclipse_state, - BlackoilOutputWriter& output_writer, - Opm::DeckConstPtr& deck, - const std::vector& threshold_pressures_by_face) - - { - pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil, has_polymer, - eclipse_state, output_writer, deck, threshold_pressures_by_face)); - } - - - - - - template - SimulatorReport SimulatorFullyImplicitBlackoilPolymer::run(SimulatorTimer& timer, - PolymerBlackoilState& state) - { - return pimpl_->run(timer, state); - } - - - // \TODO: Treat bcs. - template - SimulatorFullyImplicitBlackoilPolymer::Impl::Impl(const parameter::ParameterGroup& param, - const Grid& grid, - const DerivedGeology& geo, - BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity, - const bool has_disgas, - const bool has_vapoil, - const bool has_polymer, - std::shared_ptr eclipse_state, - BlackoilOutputWriter& output_writer, - Opm::DeckConstPtr& deck, - const std::vector& threshold_pressures_by_face) - : param_(param), - grid_(grid), - props_(props), - polymer_props_(polymer_props), - rock_comp_props_(rock_comp_props), - gravity_(gravity), - geo_(geo), - solver_(linsolver), - has_disgas_(has_disgas), - has_vapoil_(has_vapoil), - has_polymer_(has_polymer), - terminal_output_(param.getDefault("output_terminal", true)), - eclipse_state_(eclipse_state), - output_writer_(output_writer), - deck_(deck), - rateConverter_(props_, std::vector(AutoDiffGrid::numCells(grid_), 0)), - threshold_pressures_by_face_(threshold_pressures_by_face) + SimulatorFullyImplicitBlackoilPolymer:: + SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, + const Grid& grid, + const DerivedGeology& geo, + BlackoilPropsAdInterface& props, + const PolymerPropsAd& polymer_props, + const RockCompressibility* rock_comp_props, + NewtonIterationBlackoilInterface& linsolver, + const double* gravity, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + std::shared_ptr eclipse_state, + BlackoilOutputWriter& output_writer, + Opm::DeckConstPtr& deck, + const std::vector& threshold_pressures_by_face) + : BaseType(param, + grid, + geo, + props, + rock_comp_props, + linsolver, + gravity, + has_disgas, + has_vapoil, + eclipse_state, + output_writer, + threshold_pressures_by_face) + , polymer_props_(polymer_props) + , has_polymer_(has_polymer) + , deck_(deck) { // Misc init. const int num_cells = AutoDiffGrid::numCells(grid); - allcells_.resize(num_cells); + BaseType::allcells_.resize(num_cells); for (int cell = 0; cell < num_cells; ++cell) { - allcells_[cell] = cell; + BaseType::allcells_[cell] = cell; } #if HAVE_MPI - if ( terminal_output_ ) { - if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + if ( BaseType::terminal_output_ ) { + if ( BaseType::solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) { const ParallelISTLInformation& info = - boost::any_cast(solver_.parallelInformation()); + boost::any_cast(BaseType::solver_.parallelInformation()); // Only rank 0 does print to std::cout - terminal_output_= (info.communicator().rank()==0); + BaseType::terminal_output_= (info.communicator().rank()==0); } } #endif @@ -233,8 +77,8 @@ namespace Opm template - SimulatorReport SimulatorFullyImplicitBlackoilPolymer::Impl::run(SimulatorTimer& timer, - PolymerBlackoilState& state) + SimulatorReport SimulatorFullyImplicitBlackoilPolymer::run(SimulatorTimer& timer, + PolymerBlackoilState& state) { WellStateFullyImplicitBlackoilPolymer prev_well_state; @@ -244,33 +88,33 @@ namespace Opm Opm::time::StopWatch step_timer; Opm::time::StopWatch total_timer; total_timer.start(); - std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt"; + std::string tstep_filename = BaseType::output_writer_.outputDirectory() + "/step_timing.txt"; std::ofstream tstep_os(tstep_filename.c_str()); typedef T Grid; typedef BlackoilPolymerModel Model; typedef typename Model::ModelParameters ModelParams; - ModelParams modelParams( param_ ); + ModelParams modelParams( BaseType::param_ ); typedef NewtonSolver Solver; typedef typename Solver::SolverParameters SolverParams; - SolverParams solverParams( param_ ); + SolverParams solverParams( BaseType::param_ ); //adaptive time stepping // std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; - // if( param_.getDefault("timestep.adaptive", bool(false) ) ) + // if( BaseType::param_.getDefault("timestep.adaptive", bool(false) ) ) // { - // adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) ); + // adaptiveTimeStepping.reset( new AdaptiveTimeStepping( BaseType::param_ ) ); // } // init output writer - output_writer_.writeInit( timer ); + BaseType::output_writer_.writeInit( timer ); - std::string restorefilename = param_.getDefault("restorefile", std::string("") ); + std::string restorefilename = BaseType::param_.getDefault("restorefile", std::string("") ); if( ! restorefilename.empty() ) { // -1 means that we'll take the last report step that was written - const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); - output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep ); + const int desiredRestoreStep = BaseType::param_.getDefault("restorestep", int(-1) ); + BaseType::output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep ); } unsigned int totalNewtonIterations = 0; @@ -280,21 +124,21 @@ namespace Opm while (!timer.done()) { // Report timestep. step_timer.start(); - if ( terminal_output_ ) + if ( BaseType::terminal_output_ ) { timer.report(std::cout); } // Create wells and well state. - WellsManager wells_manager(eclipse_state_, + WellsManager wells_manager(BaseType::eclipse_state_, timer.currentStepNum(), - Opm::UgGridHelpers::numCells(grid_), - Opm::UgGridHelpers::globalCell(grid_), - Opm::UgGridHelpers::cartDims(grid_), - Opm::UgGridHelpers::dimensions(grid_), - Opm::UgGridHelpers::cell2Faces(grid_), - Opm::UgGridHelpers::beginFaceCentroids(grid_), - props_.permeability()); + Opm::UgGridHelpers::numCells(BaseType::grid_), + Opm::UgGridHelpers::globalCell(BaseType::grid_), + Opm::UgGridHelpers::cartDims(BaseType::grid_), + Opm::UgGridHelpers::dimensions(BaseType::grid_), + Opm::UgGridHelpers::cell2Faces(BaseType::grid_), + Opm::UgGridHelpers::beginFaceCentroids(BaseType::grid_), + BaseType::props_.permeability()); const Wells* wells = wells_manager.c_wells(); WellStateFullyImplicitBlackoilPolymer well_state; well_state.init(wells, state.blackoilState(), prev_well_state); @@ -305,34 +149,34 @@ namespace Opm if (wells_manager.c_wells() == 0) { OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); } - polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, eclipse_state_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum())); + polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, BaseType::eclipse_state_, *wells, Opm::UgGridHelpers::numCells(BaseType::grid_), timer.currentStepNum())); } else { polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day, 1.0*Opm::unit::day, 0.0)); } - std::vector polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_)); + std::vector polymer_inflow_c(Opm::UgGridHelpers::numCells(BaseType::grid_)); polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), timer.simulationTimeElapsed() + timer.currentStepLength(), polymer_inflow_c); well_state.polymerInflow() = polymer_inflow_c; // write simulation state at the report stage - output_writer_.writeTimeStep( timer, state.blackoilState(), well_state ); + BaseType::output_writer_.writeTimeStep( timer, state.blackoilState(), well_state ); // Max oil saturation (for VPPARS), hysteresis update. - props_.updateSatOilMax(state.saturation()); - props_.updateSatHyst(state.saturation(), allcells_); + BaseType::props_.updateSatOilMax(state.saturation()); + BaseType::props_.updateSatHyst(state.saturation(), BaseType::allcells_); // Compute reservoir volumes for RESV controls. - computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state); + BaseType::computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state); // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - Model model(modelParams, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_); - if (!threshold_pressures_by_face_.empty()) { - model.setThresholdPressures(threshold_pressures_by_face_); + Model model(modelParams, BaseType::grid_, BaseType::props_, BaseType::geo_, BaseType::rock_comp_props_, polymer_props_, wells, BaseType::solver_, BaseType::has_disgas_, BaseType::has_vapoil_, has_polymer_, BaseType::terminal_output_); + if (!BaseType::threshold_pressures_by_face_.empty()) { + model.setThresholdPressures(BaseType::threshold_pressures_by_face_); } Solver solver(solverParams, model); @@ -342,7 +186,7 @@ namespace Opm // \Note: The report steps are met in any case // \Note: The sub stepping will require a copy of the state variables // if( adaptiveTimeStepping ) { - // adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ ); + // adaptiveTimeStepping->step( timer, solver, state, well_state, BaseType::output_writer_ ); // } else { // solve for complete report step solver.step(timer.currentStepLength(), state, well_state); @@ -358,13 +202,13 @@ namespace Opm // Report timing. const double st = solver_timer.secsSinceStart(); - if ( terminal_output_ ) + if ( BaseType::terminal_output_ ) { std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; } stime += st; - if ( output_writer_.output() ) { + if ( BaseType::output_writer_.output() ) { SimulatorReport step_report; step_report.pressure_time = st; step_report.total_time = step_timer.secsSinceStart(); @@ -377,7 +221,7 @@ namespace Opm } // Write final simulation state. - output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state ); + BaseType::output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state ); // Stop timer and create timing report total_timer.stop(); @@ -389,229 +233,4 @@ namespace Opm report.total_linear_iterations = totalLinearIterations; return report; } - - namespace SimFIBODetails { - typedef std::unordered_map WellMap; - - inline WellMap - mapWells(const std::vector& wells) - { - WellMap wmap; - - for (std::vector::const_iterator - w = wells.begin(), e = wells.end(); - w != e; ++w) - { - wmap.insert(std::make_pair((*w)->name(), *w)); - } - - return wmap; - } - - inline int - resv_control(const WellControls* ctrl) - { - int i, n = well_controls_get_num(ctrl); - - bool match = false; - for (i = 0; (! match) && (i < n); ++i) { - match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE; - } - - if (! match) { i = 0; } - - return i - 1; // -1 if no match, undo final "++" otherwise - } - - inline bool - is_resv(const Wells& wells, - const int w) - { - return (0 <= resv_control(wells.ctrls[w])); - } - - inline bool - is_resv(const WellMap& wmap, - const std::string& name, - const std::size_t step) - { - bool match = false; - - WellMap::const_iterator i = wmap.find(name); - - if (i != wmap.end()) { - WellConstPtr wp = i->second; - - match = (wp->isProducer(step) && - wp->getProductionProperties(step) - .hasProductionControl(WellProducer::RESV)) - || (wp->isInjector(step) && - wp->getInjectionProperties(step) - .hasInjectionControl(WellInjector::RESV)); - } - - return match; - } - - inline std::vector - resvWells(const Wells* wells, - const std::size_t step, - const WellMap& wmap) - { - std::vector resv_wells; - if( wells ) - { - for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) { - if (is_resv(*wells, w) || - ((wells->name[w] != 0) && - is_resv(wmap, wells->name[w], step))) - { - resv_wells.push_back(w); - } - } - } - - return resv_wells; - } - - inline void - historyRates(const PhaseUsage& pu, - const WellProductionProperties& p, - std::vector& rates) - { - assert (! p.predictionMode); - assert (rates.size() == - std::vector::size_type(pu.num_phases)); - - if (pu.phase_used[ BlackoilPhases::Aqua ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Aqua ]; - - rates[i] = p.WaterRate; - } - - if (pu.phase_used[ BlackoilPhases::Liquid ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Liquid ]; - - rates[i] = p.OilRate; - } - - if (pu.phase_used[ BlackoilPhases::Vapour ]) { - const std::vector::size_type - i = pu.phase_pos[ BlackoilPhases::Vapour ]; - - rates[i] = p.GasRate; - } - } - } // namespace SimFIBODetails - - template - void - SimulatorFullyImplicitBlackoilPolymer:: - Impl::computeRESV(const std::size_t step, - const Wells* wells, - const BlackoilState& x, - WellStateFullyImplicitBlackoilPolymer& xw) - { - typedef SimFIBODetails::WellMap WellMap; - - const std::vector& w_ecl = eclipse_state_->getSchedule()->getWells(step); - const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); - - const std::vector& resv_wells = SimFIBODetails::resvWells(wells, step, wmap); - - if (! resv_wells.empty()) { - const PhaseUsage& pu = props_.phaseUsage(); - const std::vector::size_type np = props_.numPhases(); - - rateConverter_.defineState(x); - - std::vector distr (np); - std::vector hrates(np); - std::vector prates(np); - - for (std::vector::const_iterator - rp = resv_wells.begin(), e = resv_wells.end(); - rp != e; ++rp) - { - WellControls* ctrl = wells->ctrls[*rp]; - const bool is_producer = wells->type[*rp] == PRODUCER; - - // RESV control mode, all wells - { - const int rctrl = SimFIBODetails::resv_control(ctrl); - - if (0 <= rctrl) { - const std::vector::size_type off = (*rp) * np; - - if (is_producer) { - // Convert to positive rates to avoid issues - // in coefficient calculations. - std::transform(xw.wellRates().begin() + (off + 0*np), - xw.wellRates().begin() + (off + 1*np), - prates.begin(), std::negate()); - } else { - std::copy(xw.wellRates().begin() + (off + 0*np), - xw.wellRates().begin() + (off + 1*np), - prates.begin()); - } - - const int fipreg = 0; // Hack. Ignore FIP regions. - rateConverter_.calcCoeff(prates, fipreg, distr); - - well_controls_iset_distr(ctrl, rctrl, & distr[0]); - } - } - - // RESV control, WCONHIST wells. A bit of duplicate - // work, regrettably. - if (is_producer && wells->name[*rp] != 0) { - WellMap::const_iterator i = wmap.find(wells->name[*rp]); - - if (i != wmap.end()) { - WellConstPtr wp = i->second; - - const WellProductionProperties& p = - wp->getProductionProperties(step); - - if (! p.predictionMode) { - // History matching (WCONHIST/RESV) - SimFIBODetails::historyRates(pu, p, hrates); - - const int fipreg = 0; // Hack. Ignore FIP regions. - rateConverter_.calcCoeff(hrates, fipreg, distr); - - // WCONHIST/RESV target is sum of all - // observed phase rates translated to - // reservoir conditions. Recall sign - // convention: Negative for producers. - const double target = - - std::inner_product(distr.begin(), distr.end(), - hrates.begin(), 0.0); - - well_controls_clear(ctrl); - well_controls_assert_number_of_phases(ctrl, int(np)); - - const int ok_resv = - well_controls_add_new(RESERVOIR_RATE, target, - & distr[0], ctrl); - - // For WCONHIST/RESV the BHP limit is set to 1 atm. - // TODO: Make it possible to modify the BHP limit using - // the WELTARG keyword - const int ok_bhp = - well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm), - NULL, ctrl); - - if (ok_resv != 0 && ok_bhp != 0) { - xw.currentControls()[*rp] = 0; - well_controls_set_current(ctrl, 0); - } - } - } - } - } - } - } } // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp index bcdcee66d..e69de29bb 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp @@ -1,473 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2014 STATOIL ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - - -#if HAVE_CONFIG_H -#include "config.h" -#endif // HAVE_CONFIG_H - -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - - -namespace Opm -{ - - - - namespace - { - static void outputStateVtk(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir); - static void outputStateMatlab(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir); - static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir); - } // anonymous namespace - - class SimulatorFullyImplicitCompressiblePolymer::Impl - { - public: - Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const DerivedGeology& geo, - const BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - std::shared_ptr eclipse_state, - EclipseWriter& output_writer, - Opm::DeckConstPtr& deck, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity); - - SimulatorReport run(SimulatorTimer& timer, - PolymerBlackoilState& state); - - private: - // Data. - - // Parameters for output. - bool output_; - bool output_vtk_; - std::string output_dir_; - int output_interval_; - // Parameters for well control - bool check_well_controls_; - int max_well_control_iterations_; - // Observed objects. - const UnstructuredGrid& grid_; - const BlackoilPropsAdInterface& props_; - const PolymerPropsAd& polymer_props_; - const RockCompressibility* rock_comp_props_; - std::shared_ptr eclipse_state_; - EclipseWriter& output_writer_; - Opm::DeckConstPtr& deck_; - NewtonIterationBlackoilInterface& linsolver_; - const double* gravity_; - // Solvers - DerivedGeology geo_; - // Misc. data - std::vector allcells_; - }; - - - - - SimulatorFullyImplicitCompressiblePolymer:: - SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const DerivedGeology& geo, - const BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - std::shared_ptr eclipse_state, - EclipseWriter& output_writer, - Opm::DeckConstPtr& deck, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity) - - { - pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, eclipse_state, output_writer, deck, linsolver, gravity)); - } - - - - - - SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer, - PolymerBlackoilState& state) - { - return pimpl_->run(timer, state); - } - - - - - - // \TODO: Treat bcs. - SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param, - const UnstructuredGrid& grid, - const DerivedGeology& geo, - const BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, - const RockCompressibility* rock_comp_props, - std::shared_ptr eclipse_state, - EclipseWriter& output_writer, - Opm::DeckConstPtr& deck, - NewtonIterationBlackoilInterface& linsolver, - const double* gravity) - : grid_(grid), - props_(props), - polymer_props_(polymer_props), - rock_comp_props_(rock_comp_props), - eclipse_state_(eclipse_state), - output_writer_(output_writer), - deck_(deck), - linsolver_(linsolver), - gravity_(gravity), - geo_(geo) - { - // For output. - output_ = param.getDefault("output", true); - if (output_) { - output_vtk_ = param.getDefault("output_vtk", true); - output_dir_ = param.getDefault("output_dir", std::string("output")); - // Ensure that output dir exists - boost::filesystem::path fpath(output_dir_); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - output_interval_ = param.getDefault("output_interval", 1); - } - - // Well control related init. - check_well_controls_ = param.getDefault("check_well_controls", false); - max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10); - - // Misc init. - const int num_cells = grid.number_of_cells; - allcells_.resize(num_cells); - for (int cell = 0; cell < num_cells; ++cell) { - allcells_[cell] = cell; - } - } - - - - - SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer, - PolymerBlackoilState& state) - { - WellStateFullyImplicitBlackoil prev_well_state; - // Initialisation. - std::vector porevol; - if (rock_comp_props_ && rock_comp_props_->isActive()) { - computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); - } else { - computePorevolume(grid_, props_.porosity(), porevol); - } - std::vector initial_porevol = porevol; - - std::vector polymer_inflow_c(grid_.number_of_cells); - // Main simulation loop. - Opm::time::StopWatch solver_timer; - double stime = 0.0; - Opm::time::StopWatch step_timer; - Opm::time::StopWatch total_timer; - total_timer.start(); - std::string tstep_filename = output_dir_ + "/step_timing.txt"; - std::ofstream tstep_os(tstep_filename.c_str()); - - //Main simulation loop. - while (!timer.done()) { -#if 0 - double tot_injected[2] = { 0.0 }; - double tot_produced[2] = { 0.0 }; - Opm::Watercut watercut; - watercut.push(0.0, 0.0, 0.0); - std::vector fractional_flows; - std::vector well_resflows_phase; - if (wells_) { - well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0); - } - std::fstream tstep_os; - if (output_) { - std::string filename = output_dir_ + "/step_timing.param"; - tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app); - } -#endif - // Report timestep and (optionally) write state to disk. - - step_timer.start(); - timer.report(std::cout); - - WellsManager wells_manager(eclipse_state_, - timer.currentStepNum(), - Opm::UgGridHelpers::numCells(grid_), - Opm::UgGridHelpers::globalCell(grid_), - Opm::UgGridHelpers::cartDims(grid_), - Opm::UgGridHelpers::dimensions(grid_), - Opm::UgGridHelpers::cell2Faces(grid_), - Opm::UgGridHelpers::beginFaceCentroids(grid_), - props_.permeability()); - const Wells* wells = wells_manager.c_wells(); - WellStateFullyImplicitBlackoil well_state; - well_state.init(wells, state.blackoilState(), prev_well_state); - //Compute polymer inflow. - std::unique_ptr polymer_inflow_ptr; - if (deck_->hasKeyword("WPOLYMER")) { - if (wells_manager.c_wells() == 0) { - OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells."); - } - polymer_inflow_ptr.reset(new PolymerInflowFromDeck(deck_, eclipse_state_, *wells, Opm::UgGridHelpers::numCells(grid_), timer.currentStepNum())); - } else { - polymer_inflow_ptr.reset(new PolymerInflowBasic(0.0*Opm::unit::day, - 1.0*Opm::unit::day, - 0.0)); - } - std::vector polymer_inflow_c(Opm::UgGridHelpers::numCells(grid_)); - polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), - timer.simulationTimeElapsed() + timer.currentStepLength(), - polymer_inflow_c); - - if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - } - if (output_) { - if (timer.currentStepNum() == 0) { - output_writer_.writeInit(timer); - } - output_writer_.writeTimeStep(timer, state.blackoilState(), well_state); - } - // Run solver. - solver_timer.start(); - FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_); - solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); - // Stop timer and report. - solver_timer.stop(); - const double st = solver_timer.secsSinceStart(); - std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; - - stime += st; - // Update pore volumes if rock is compressible. - if (rock_comp_props_ && rock_comp_props_->isActive()) { - initial_porevol = porevol; - computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); - } -/* - double injected[2] = { 0.0 }; - double produced[2] = { 0.0 }; - double polyinj = 0; - double polyprod = 0; - Opm::computeInjectedProduced(props_, polymer_props_, - state, - transport_src, polymer_inflow_c, timer.currentStepLength(), - injected, produced, - polyinj, polyprod); - tot_injected[0] += injected[0]; - tot_injected[1] += injected[1]; - tot_produced[0] += produced[0]; - tot_produced[1] += produced[1]; - watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(), - produced[0]/(produced[0] + produced[1]), - tot_produced[0]/tot_porevol_init); - std::cout.precision(5); - const int width = 18; - std::cout << "\nMass balance report.\n"; - std::cout << " Injected reservoir volumes: " - << std::setw(width) << injected[0] - << std::setw(width) << injected[1] << std::endl; - std::cout << " Produced reservoir volumes: " - << std::setw(width) << produced[0] - << std::setw(width) << produced[1] << std::endl; - std::cout << " Total inj reservoir volumes: " - << std::setw(width) << tot_injected[0] - << std::setw(width) << tot_injected[1] << std::endl; - std::cout << " Total prod reservoir volumes: " - << std::setw(width) << tot_produced[0] - << std::setw(width) << tot_produced[1] << std::endl; -*/ - if (output_) { - SimulatorReport step_report; - step_report.pressure_time = st; - step_report.total_time = step_timer.secsSinceStart(); - step_report.reportParam(tstep_os); - } - ++timer; - prev_well_state = well_state; - } - // Write final simulation state. - if (output_) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state); - } - - total_timer.stop(); - SimulatorReport report; - report.pressure_time = stime; - report.transport_time = 0.0; - report.total_time = total_timer.secsSinceStart(); - return report; - } - - - - namespace - { - - static void outputStateVtk(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - // Write data in VTK format. - std::ostringstream vtkfilename; - vtkfilename << output_dir << "/vtk_files"; - boost::filesystem::path fpath(vtkfilename.str()); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; - std::ofstream vtkfile(vtkfilename.str().c_str()); - if (!vtkfile) { - OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); - } - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - dm["cmax"] = &state.maxconcentration(); - dm["concentration"] = &state.concentration(); - std::vector cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - Opm::writeVtkData(grid, dm, vtkfile); - } - - - static void outputStateMatlab(const UnstructuredGrid& grid, - const Opm::PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - dm["cmax"] = &state.maxconcentration(); - dm["concentration"] = &state.concentration(); - dm["surfvolume"] = &state.surfacevol(); - std::vector cell_velocity; - Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - -#if 0 - static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) - { - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - watercut.write(os); - } -#endif - } - -} // namespace Opm diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp index df868c8af..ecd33c71b 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -21,67 +21,70 @@ #ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED -#include -#include -#include +#include +#include -struct UnstructuredGrid; -struct Wells; +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include namespace Opm { - namespace parameter { class ParameterGroup; } - class BlackoilPropsAdInterface; - class RockCompressibility; - class DerivedGeology; - class WellStateFullyImplicitBlackoil; - class WellsManager; - class EclipseWriter; - class EclipseState; - class NewtonIterationBlackoilInterface; - class SimulatorTimer; - class PolymerBlackoilState; - class PolymerPropsAd; - class PolymerInflowInterface; - struct SimulatorReport; - /// Class collecting all necessary components for a two-phase simulation. + template class SimulatorFullyImplicitCompressiblePolymer + : public SimulatorBase > { + typedef SimulatorFullyImplicitCompressiblePolymer ThisType; + typedef SimulatorBase BaseType; + public: /// Initialise from parameters and objects to observe. - /// \param[in] param parameters, this class accepts the following: - /// parameter (default) effect - /// ----------------------------------------------------------- - /// output (true) write output to files? - /// output_dir ("output") output directoty - /// output_interval (1) output every nth step - /// nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal) - /// nl_pressure_change_tolerance (1.0) pressure solver change tolerance (in Pascal) - /// nl_pressure_maxiter (10) max nonlinear iterations in pressure - /// nl_maxiter (30) max nonlinear iterations in transport - /// nl_tolerance (1e-9) transport solver absolute residual tolerance - /// num_transport_substeps (1) number of transport steps per pressure step - /// use_segregation_split (false) solve for gravity segregation (if false, - /// segregation is ignored). - /// - /// \param[in] grid grid data structure - /// \param[in] props fluid and rock properties - /// \param[in] polymer_props polymer properties - /// \param[in] rock_comp_props if non-null, rock compressibility properties - /// \param[in] eclipse_state - /// \param[in] eclipse_writer - /// \param[in] deck - /// \param[in] linsolver linear solver - /// \param[in] gravity if non-null, gravity vector SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param, const UnstructuredGrid& grid, const DerivedGeology& geo, - const BlackoilPropsAdInterface& props, + BlackoilPropsAdInterface& props, const PolymerPropsAd& polymer_props, const RockCompressibility* rock_comp_props, std::shared_ptr eclipse_state, - EclipseWriter& eclipse_writer, + BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, NewtonIterationBlackoilInterface& linsolver, const double* gravity); @@ -95,12 +98,21 @@ namespace Opm SimulatorReport run(SimulatorTimer& timer, PolymerBlackoilState& state); - private: - class Impl; - // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl. - std::shared_ptr pimpl_; +private: + Opm::DeckConstPtr deck_; + const PolymerPropsAd& polymer_props_; + +#warning "To remove" + bool output_; + int output_interval_; + bool output_vtk_; + std::string output_dir_; + bool check_well_controls_; + int max_well_control_iterations_; }; } // namespace Opm +#include "SimulatorFullyImplicitCompressiblePolymer_impl.hpp" + #endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED