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.

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp

@@ -195,13 +195,14 @@ namespace {
 
 
 
-    void
+    int
     FullyImplicitCompressiblePolymerSolver::
     step(const double          dt,
          PolymerBlackoilState& x ,
-         WellStateFullyImplicitBlackoil& xw,
+         WellStateFullyImplicitBlackoilPolymer& xw)
-         const std::vector<double>& polymer_inflow)
     {
+        const std::vector<double>& polymer_inflow = xw.polymerInflow();
+
         // Initial max concentration of this time step from PolymerBlackoilState.
         cmax_ = Eigen::Map<V>(&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<Eigen::Infinity>();
-        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;
     }
 
 

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -27,6 +27,7 @@
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
 #include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
 #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
 
 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,
+             WellStateFullyImplicitBlackoilPolymer& wstate);
-             const std::vector<double>& polymer_inflow);
+
+        int newtonIterations() const
+        { return newtonIterations_; }
+        int linearIterations() const
+        { return linearIterations_; }
 
     private:
         typedef AutoDiffBlock<double> 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,

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp

@@ -88,6 +88,7 @@ namespace Opm
         typedef BlackoilOutputWriter OutputWriter;
         typedef UnstructuredGrid Grid;
         typedef BlackoilPolymerModel<Grid> Model;
+        typedef NewtonSolver<Model> Solver;
     };
 
     /// Class collecting all necessary components for a blackoil simulation with polymer
@@ -98,6 +99,8 @@ namespace Opm
         typedef SimulatorFullyImplicitBlackoilPolymer ThisType;
         typedef SimulatorBase<ThisType> BaseType;
 
+        typedef SimulatorTraits<ThisType> Traits;
+
     public:
         SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
                                               const typename BaseType::Grid& grid,
@@ -115,21 +118,33 @@ namespace Opm
                                               Opm::DeckConstPtr& deck,
                                               const std::vector<double>& threshold_pressures_by_face);
 
-        std::shared_ptr<Model> createModel(const typename Model::ModelParameters &modelParams,
+        std::shared_ptr<Solver> createSolver(const Wells* wells)
-                                           const Wells* wells)
         {
-            return std::make_shared<Model>(modelParams,
+            typedef typename Traits::Model Model;
-                                           BaseType::grid_,
+            typedef typename Model::ModelParameters ModelParams;
-                                           BaseType::props_,
+            ModelParams modelParams( param_ );
-                                           BaseType::geo_,
+            typedef NewtonSolver<Model> Solver;
-                                           BaseType::rock_comp_props_,
+
-                                           polymer_props_,
+            auto model = std::make_shared<Model>(modelParams,
-                                           wells,
+                                                 BaseType::grid_,
-                                           BaseType::solver_,
+                                                 BaseType::props_,
-                                           BaseType::has_disgas_,
+                                                 BaseType::geo_,
-                                           BaseType::has_vapoil_,
+                                                 BaseType::rock_comp_props_,
-                                           has_polymer_,
+                                                 polymer_props_,
-                                           BaseType::terminal_output_);
+                                                 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<Solver>(solverParams, model);
         }
 
         void handleAdditionalWellInflow(SimulatorTimer& timer,

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp

@@ -74,11 +74,10 @@ namespace Opm
     struct SimulatorTraits<SimulatorFullyImplicitCompressiblePolymer>
     {
         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 FullyImplicitCompressiblePolymerSolver Solver;
-        typedef BlackoilPolymerModel<Grid> Model;
     };
 
     /// Class collecting all necessary components for a two-phase simulation.
@@ -102,26 +101,45 @@ namespace Opm
                                        			  NewtonIterationBlackoilInterface& linsolver,
                                        			  const double* gravity);
 
-        /// Run the simulation.
+        std::shared_ptr<Solver> createSolver(const Wells* wells)
-        /// This will run succesive timesteps until timer.done() is true. It will
+        {
-        /// modify the reservoir and well states.
+            return std::make_shared<Solver>(BaseType::grid_,
-        /// \param[in,out] timer       governs the requested reporting timesteps
+                                            BaseType::props_,
-        /// \param[in,out] state       state of reservoir: pressure, fluxes
+                                            BaseType::geo_,
-        /// \return                    simulation report, with timing data
+                                            BaseType::rock_comp_props_,
-        SimulatorReport run(SimulatorTimer& timer,
+                                            polymer_props_,
-                            typename BaseType::ReservoirState& state);
+                                            *wells,
+                                            BaseType::solver_);
+        }
+
+        void handleAdditionalWellInflow(SimulatorTimer& timer,
+                                        WellsManager& wells_manager,
+                                        typename BaseType::WellState& well_state,
+                                        const Wells* wells)
+        {
+            // compute polymer inflow
+            std::unique_ptr<PolymerInflowInterface> 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<double> 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

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<double> 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<double> 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<double>& d = *(it->second);
-        std::copy(d.begin(), d.end(), std::ostream_iterator<double>(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<double> 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<double> initial_porevol = porevol;
-
-    std::vector<double> 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<double> fractional_flows;
-        std::vector<double> 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<PolymerInflowInterface> 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<double> 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