From 01555da823b0f8a88ed9c53fd176ac102611977c Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 28 May 2015 12:20:59 +0200 Subject: [PATCH] remove the run() method from SimulatorFullyImplicitCompressiblePolymer With this the simulator is basically done, but since FullyImplicitCompressiblePolymerSolver has not yet been converted to the NewtonSolver plus Model approach, the solver cannot be removed and thus still contains quite a bit of copy-and-pasted code. --- ...FullyImplicitCompressiblePolymerSolver.cpp | 16 +- ...FullyImplicitCompressiblePolymerSolver.hpp | 15 +- .../SimulatorFullyImplicitBlackoilPolymer.hpp | 43 ++- ...ulatorFullyImplicitCompressiblePolymer.hpp | 56 ++-- ...rFullyImplicitCompressiblePolymer_impl.hpp | 279 ------------------ 5 files changed, 90 insertions(+), 319 deletions(-) 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