diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 7ceec001e..c9728d8c3 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -195,13 +195,14 @@ namespace { - void + int FullyImplicitCompressiblePolymerSolver:: step(const double dt, PolymerBlackoilState& x , - WellStateFullyImplicitBlackoil& xw, - const std::vector& polymer_inflow) + WellStateFullyImplicitBlackoilPolymer& xw) { + const std::vector& polymer_inflow = xw.polymerInflow(); + // Initial max concentration of this time step from PolymerBlackoilState. cmax_ = Eigen::Map(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_)); @@ -215,7 +216,7 @@ namespace { const double r0 = residualNorm(); const double r_polymer = residual_.material_balance_eq[2].value().matrix().lpNorm(); - int it = 0; + int it = 0; std::cout << "\nIteration Residual Polymer Res\n" << std::setw(9) << it << std::setprecision(9) << std::setw(18) << r0 << std::setprecision(9) @@ -224,6 +225,9 @@ namespace { while (resTooLarge && (it < maxit)) { const V dx = solveJacobianSystem(); + // update the number of linear iterations used. + linearIterations_ += linsolver_.iterations(); + updateState(dx, x, xw); assemble(dt, x, xw, polymer_inflow); @@ -233,6 +237,7 @@ namespace { resTooLarge = (r > atol) && (r > rtol*r0); it += 1; + newtonIterations_ += 1; std::cout << std::setw(9) << it << std::setprecision(9) << std::setw(18) << r << std::setprecision(9) << std::setw(18) << rr_polymer << std::endl; @@ -240,11 +245,14 @@ namespace { if (resTooLarge) { std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; + return -1; // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); } // Update max concentration. computeCmax(x); + + return it; } diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index d2db12108..43e5f27a5 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include struct UnstructuredGrid; @@ -79,11 +80,15 @@ namespace Opm { /// \param[in] state reservoir state /// \param[in] wstate well state /// \param[in] polymer_inflow polymer influx - void + int step(const double dt, PolymerBlackoilState& state , - WellStateFullyImplicitBlackoil& wstate, - const std::vector& polymer_inflow); + WellStateFullyImplicitBlackoilPolymer& wstate); + + int newtonIterations() const + { return newtonIterations_; } + int linearIterations() const + { return linearIterations_; } private: typedef AutoDiffBlock ADB; @@ -142,6 +147,10 @@ namespace Opm { // each of which has size equal to the number of cells. // The well_eq has size equal to the number of wells. LinearisedBlackoilResidual residual_; + + unsigned int newtonIterations_; + unsigned int linearIterations_; + // Private methods. SolutionState constantState(const PolymerBlackoilState& x, diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 48c3a6e1a..ab78bcbb2 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -88,6 +88,7 @@ namespace Opm typedef BlackoilOutputWriter OutputWriter; typedef UnstructuredGrid Grid; typedef BlackoilPolymerModel Model; + typedef NewtonSolver Solver; }; /// Class collecting all necessary components for a blackoil simulation with polymer @@ -98,6 +99,8 @@ namespace Opm typedef SimulatorFullyImplicitBlackoilPolymer ThisType; typedef SimulatorBase BaseType; + typedef SimulatorTraits Traits; + public: SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param, const typename BaseType::Grid& grid, @@ -115,21 +118,33 @@ namespace Opm Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face); - std::shared_ptr createModel(const typename Model::ModelParameters &modelParams, - const Wells* wells) + std::shared_ptr createSolver(const Wells* wells) { - return std::make_shared(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_); + typedef typename Traits::Model Model; + typedef typename Model::ModelParameters ModelParams; + ModelParams modelParams( param_ ); + typedef NewtonSolver Solver; + + auto model = std::make_shared(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 (!threshold_pressures_by_face_.empty()) { + model->setThresholdPressures(threshold_pressures_by_face_); + } + + typedef typename Solver::SolverParameters SolverParams; + SolverParams solverParams( param_ ); + return std::make_shared(solverParams, model); } void handleAdditionalWellInflow(SimulatorTimer& timer, diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp index f3c7e38dd..94a8db36c 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -74,11 +74,10 @@ namespace Opm struct SimulatorTraits { typedef PolymerBlackoilState ReservoirState; - typedef WellStateFullyImplicitBlackoil WellState; + typedef WellStateFullyImplicitBlackoilPolymer WellState; typedef BlackoilOutputWriter OutputWriter; typedef UnstructuredGrid Grid; -#warning TODO: the 2p polymer solver does not yet adhere to The New Order! - typedef BlackoilPolymerModel Model; + typedef FullyImplicitCompressiblePolymerSolver Solver; }; /// Class collecting all necessary components for a two-phase simulation. @@ -102,26 +101,45 @@ namespace Opm NewtonIterationBlackoilInterface& linsolver, const double* gravity); - /// 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 - /// \return simulation report, with timing data - SimulatorReport run(SimulatorTimer& timer, - typename BaseType::ReservoirState& state); + std::shared_ptr createSolver(const Wells* wells) + { + return std::make_shared(BaseType::grid_, + BaseType::props_, + BaseType::geo_, + BaseType::rock_comp_props_, + polymer_props_, + *wells, + BaseType::solver_); + } + + void handleAdditionalWellInflow(SimulatorTimer& timer, + WellsManager& wells_manager, + typename BaseType::WellState& well_state, + const Wells* wells) + { + // 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_, 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(BaseType::grid_)); + polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), + timer.simulationTimeElapsed() + timer.currentStepLength(), + polymer_inflow_c); + well_state.polymerInflow() = polymer_inflow_c; + } + 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 diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp index ab973b5fa..e8af9fce5 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer_impl.hpp @@ -22,92 +22,6 @@ namespace Opm { -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 -} /// Class collecting all necessary components for a two-phase simulation. SimulatorFullyImplicitCompressiblePolymer:: @@ -138,199 +52,6 @@ SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param , polymer_props_(polymer_props) { - // 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; - BaseType::allcells_.resize(num_cells); - for (int cell = 0; cell < num_cells; ++cell) { - BaseType::allcells_[cell] = cell; - } -} - -SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer, - typename BaseType::ReservoirState& state) -{ - WellStateFullyImplicitBlackoil prev_well_state; - // Initialisation. - std::vector porevol; - if (BaseType::rock_comp_props_ && BaseType::rock_comp_props_->isActive()) { - computePorevolume(BaseType::grid_, - BaseType::props_.porosity(), - *BaseType::rock_comp_props_, - state.pressure(), - porevol); - } else { - computePorevolume(BaseType::grid_, - BaseType::props_.porosity(), - porevol); - } - std::vector initial_porevol = porevol; - - std::vector polymer_inflow_c(BaseType::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(BaseType::eclipse_state_, - timer.currentStepNum(), - 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(); - WellStateFullyImplicitBlackoil well_state; - well_state.init(wells, state, 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_, 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(BaseType::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(BaseType::grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(BaseType::grid_, state, timer.currentStepNum(), output_dir_); - } - if (output_) { - if (timer.currentStepNum() == 0) { - BaseType::output_writer_.writeInit(timer); - } - BaseType::output_writer_.writeTimeStep(timer, state, well_state); - } - // Run solver. - solver_timer.start(); - FullyImplicitCompressiblePolymerSolver solver(BaseType::grid_, BaseType::props_, BaseType::geo_, BaseType::rock_comp_props_, polymer_props_, *wells_manager.c_wells(), BaseType::solver_); - 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 (BaseType::rock_comp_props_ && BaseType::rock_comp_props_->isActive()) { - initial_porevol = porevol; - computePorevolume(BaseType::grid_, BaseType::props_.porosity(), *BaseType::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(BaseType::grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(BaseType::grid_, state, timer.currentStepNum(), output_dir_); - BaseType::output_writer_.writeTimeStep(timer, state, 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 Opm