let sim_poly_fi2p_comp_ad works as the same way as sim_poly_fibo_ad.

- move time iteration to simulator class
 - add EclipseWriter and output Eclipse binaries.
 - add DeckConstPtr for polymer inflow.

examples/sim_poly_fi2p_comp_ad.cpp

@@ -36,6 +36,7 @@
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
+#include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/props/BlackoilPropertiesBasic.hpp>
 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
@@ -58,7 +59,7 @@
 #include <opm/autodiff/GeoProps.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 
 #include <boost/filesystem.hpp>
 #include <boost/algorithm/string.hpp>
@@ -105,45 +106,34 @@ try
     std::shared_ptr<BlackoilPropertiesInterface> props;
     std::shared_ptr<BlackoilPropsAdInterface> new_props;
     std::shared_ptr<RockCompressibility> rock_comp;
-    Opm::DeckConstPtr deck;
-    EclipseStateConstPtr eclipseState;
     PolymerBlackoilState state;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     std::string deck_filename = param.get<std::string>("deck_filename");
-    ParserPtr parser(new Opm::Parser());
+
-    deck = parser->parseFile(deck_filename);
+    Opm::ParserPtr newParser(new Opm::Parser());
-    eclipseState.reset(new Opm::EclipseState(deck));
+    bool strict_parsing = param.getDefault("strict_parsing", true);
+    Opm::DeckConstPtr deck = newParser->parseFile(deck_filename, strict_parsing);
+    std::shared_ptr<EclipseState> eclipseState(new EclipseState(deck));
     // Grid init
     std::vector<double> porv;
     if (eclipseState->hasDoubleGridProperty("PORV")) {
         porv = eclipseState->getDoubleGridProperty("PORV")->getData();
     }
     grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv));
-//    grid.reset(new GridManager(deck));
+    auto &cGrid = *grid->c_grid();
-
+    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
-    // use the capitalized part of the deck's filename between the
+    Opm::EclipseWriter outputWriter(param,
-    // last '/' and the last '.' character as base name.
+                                    eclipseState,
-/*
+                                    pu,
-    std::string baseName = deck_filename;
+                                    cGrid.number_of_cells,
-    auto charPos = baseName.rfind('/');
+                                    cGrid.global_cell);
-    if (charPos != std::string::npos)
-        baseName = baseName.substr(charPos + 1);
-    charPos = baseName.rfind('.');
-    if (charPos != std::string::npos)
-        baseName = baseName.substr(0, charPos);
-    baseName = boost::to_upper_copy(baseName);
-
-    Opm::EclipseWriter outputWriter(param, share_obj(*deck), share_obj(*grid->c_grid()));
-*/
     // Rock and fluid init
-    props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
+    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);
-    // check_well_controls = param.getDefault("check_well_controls", false);
-    // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
     // Rock compressibility.
     rock_comp.reset(new RockCompressibility(deck, eclipseState));
     // Gravity.
@@ -182,18 +172,14 @@ try
         param.writeParam(output_dir + "/simulation.param");
     }
 
+    Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
+    SimulatorTimer simtimer;
+    simtimer.init(timeMap);
 
-    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
-              << std::flush;
 
     SimulatorReport rep;
     // With a deck, we may have more epochs etc.
     WellState well_state;
-    int step = 0;
-    Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
-    SimulatorTimer simtimer;
-    simtimer.init(timeMap);
-    const double total_time = simtimer.totalTime();
     // Check for WPOLYMER presence in last epoch to decide
     // polymer injection control type.
     const bool use_wpolymer = deck->hasKeyword("WPOLYMER");
@@ -203,62 +189,32 @@ try
                         "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
         }
     }
-    for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
-        simtimer.setCurrentStepNum(reportStepIdx);
+              << std::flush;
-
+    SimulatorReport fullReport;
-        // Report on start of step.
+    // Create and run simulator.
-        std::cout << "\n\n--------------    Starting report step " << reportStepIdx << "    --------------"
+    Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
-                  << "\n                  (number of remaining steps: "
+    SimulatorFullyImplicitCompressiblePolymer simulator(param,
-                  << simtimer.numSteps() - step << ")\n\n" << std::flush;
+                                             *grid->c_grid(),
-
+                                             geology,
-        // Create new wells, polymer inflow controls.
+                                             *new_props,
-        WellsManager wells(eclipseState, reportStepIdx, *grid->c_grid(), props->permeability());
+                                             polymer_props_ad,
-        std::unique_ptr<PolymerInflowInterface> polymer_inflow;
+                                             rock_comp->isActive() ? rock_comp.get() : 0,
-        if (use_wpolymer) {
+                                             eclipseState,
-            if (wells.c_wells() == 0) {
+                                             outputWriter,
-                OPM_THROW(std::runtime_error, "Cannot control polymer injection via WPOLYMER without wells.");
+                                             deck,
-            }
+                                             *fis_solver,
-            polymer_inflow.reset(new PolymerInflowFromDeck(deck, *wells.c_wells(), props->numCells()));
+                                             grav);
-        } else {
+    fullReport= simulator.run(simtimer, state);
-            polymer_inflow.reset(new PolymerInflowBasic(param.getDefault("poly_start_days", 300.0)*Opm::unit::day,
-                                                        param.getDefault("poly_end_days", 800.0)*Opm::unit::day,
-                                                        param.getDefault("poly_amount", polymer_props.cMax())));
-        }
-        // @@@ HACK: we should really make a new well state and
-        // properly transfer old well state to it every epoch,
-        // since number of wells may change etc.
-        if (reportStepIdx == 0) {
-            well_state.init(wells.c_wells(), state);
-        }
-
-        // 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,
-                                                 wells,
-                                                 *polymer_inflow,
-                                                 *fis_solver,
-                                                 grav);
-        if (reportStepIdx == 0) {
-            warnIfUnusedParams(param);
-        }
-        SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
-        // Update total timing report and remember step number.
-        rep += epoch_rep;
-        step = simtimer.currentStepNum();
-    }
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
-    rep.report(std::cout);
+    fullReport.report(std::cout);
 
     if (output) {
         std::string filename = output_dir + "/walltime.param";
         std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
-        rep.reportParam(tot_os);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
     }
 
 }

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp

@@ -29,6 +29,8 @@
 
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
 #include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
 #include <opm/polymer/fullyimplicit/utilities.hpp>
 #include <opm/core/grid.h>
@@ -53,6 +55,11 @@
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
 
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleEnums.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/WellProductionProperties.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/scoped_ptr.hpp>
 #include <boost/lexical_cast.hpp>
@@ -90,14 +97,14 @@ namespace Opm
              const BlackoilPropsAdInterface& props,
              const PolymerPropsAd&          polymer_props,
              const RockCompressibility* rock_comp_props,
-             WellsManager& wells_manager,
+             std::shared_ptr<EclipseState> eclipse_state,
-             PolymerInflowInterface& polymer_inflow,
+             EclipseWriter& output_writer,
+             Opm::DeckConstPtr& deck,
              NewtonIterationBlackoilInterface& linsolver,
              const double* gravity);
 
         SimulatorReport run(SimulatorTimer& timer,
-                            PolymerBlackoilState& state,
+                            PolymerBlackoilState& state);
-                            WellState& well_state);
 
     private:
         // Data.
@@ -115,13 +122,13 @@ namespace Opm
         const BlackoilPropsAdInterface& props_;
         const PolymerPropsAd&       polymer_props_;
         const RockCompressibility* rock_comp_props_;
-        WellsManager& wells_manager_;
+        std::shared_ptr<EclipseState> eclipse_state_;
-        const Wells* wells_;
+        EclipseWriter& output_writer_;
-        PolymerInflowInterface& polymer_inflow_;
+        Opm::DeckConstPtr& deck_;
+        NewtonIterationBlackoilInterface& linsolver_;
         const double* gravity_;
         // Solvers
         DerivedGeology geo_;
-        FullyImplicitCompressiblePolymerSolver solver_;
         // Misc. data
         std::vector<int> allcells_;
     };
@@ -131,18 +138,19 @@ namespace Opm
 
     SimulatorFullyImplicitCompressiblePolymer::
     SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
-                                                                   const UnstructuredGrid& grid,
+                                              const UnstructuredGrid& grid,
-                                                          const DerivedGeology& geo,
+                                              const DerivedGeology& geo,
-                                                                   const BlackoilPropsAdInterface& props,
+                                              const BlackoilPropsAdInterface& props,
-                                                                   const PolymerPropsAd&    polymer_props,
+                                              const PolymerPropsAd& polymer_props,
-                                                                   const RockCompressibility* rock_comp_props,
+                                              const RockCompressibility* rock_comp_props,
-                                                                   WellsManager& wells_manager,
+                                              std::shared_ptr<EclipseState> eclipse_state,
-                                                                   PolymerInflowInterface&  polymer_inflow,
+                                              EclipseWriter& output_writer,
-                                                                   NewtonIterationBlackoilInterface& linsolver,
+                                              Opm::DeckConstPtr& deck,
-                                                                   const double* gravity)
+                                              NewtonIterationBlackoilInterface& linsolver,
+                                              const double* gravity)
 
     {
-        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, wells_manager, polymer_inflow, linsolver, gravity));
+        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, eclipse_state, output_writer, deck, linsolver, gravity));
     }
 
 
@@ -150,10 +158,9 @@ namespace Opm
 
 
     SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer,
-                                                        PolymerBlackoilState& state,
+                                                        PolymerBlackoilState& state)
-                                                        WellState& well_state)
     {
-        return pimpl_->run(timer, state, well_state);
+        return pimpl_->run(timer, state);
     }
 
 
@@ -165,27 +172,23 @@ namespace Opm
                    			                              const UnstructuredGrid& grid,
                                                           const DerivedGeology& geo,
                               			                  const BlackoilPropsAdInterface& props,
-                                         			      const PolymerPropsAd&    polymer_props,
+                                         			      const PolymerPropsAd& polymer_props,
                                   			              const RockCompressibility* rock_comp_props,
-                                  		                  WellsManager& wells_manager,
+                                                          std::shared_ptr<EclipseState> eclipse_state,
-                                 			              PolymerInflowInterface&  polymer_inflow,
+                                                          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),
-          wells_manager_(wells_manager),
+          eclipse_state_(eclipse_state),
-          wells_(wells_manager.c_wells()),
+          output_writer_(output_writer),
-          polymer_inflow_(polymer_inflow),
+          deck_(deck),
+          linsolver_(linsolver),
           gravity_(gravity),
-          geo_(geo),
+          geo_(geo)
-          solver_(grid_, props_, geo_, rock_comp_props, polymer_props, *wells_manager.c_wells(), linsolver)
-
-          /*                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
-                               param.getDefault("nl_pressure_change_tolerance", 1.0),
-                               param.getDefault("nl_pressure_maxiter", 10),
-                               gravity,  */
     {
         // For output.
         output_ = param.getDefault("output", true);
@@ -219,10 +222,9 @@ namespace Opm
 
 
     SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer,
-                                                              PolymerBlackoilState& state,
+                                                                         PolymerBlackoilState& state)
-                                                              WellState& well_state)
     {
-
+        WellStateFullyImplicitBlackoil prev_well_state;
         // Initialisation.
         std::vector<double> porevol;
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -241,134 +243,146 @@ namespace Opm
         Opm::time::StopWatch step_timer;
         Opm::time::StopWatch total_timer;
         total_timer.start();
-        double tot_injected[2] = { 0.0 };
+        std::string tstep_filename = output_dir_ + "/step_timing.txt";
-        double tot_produced[2] = { 0.0 };
+        std::ofstream tstep_os(tstep_filename.c_str());
-        Opm::Watercut watercut;
+
-        watercut.push(0.0, 0.0, 0.0);
+        //Main simulation loop.
+        while (!timer.done()) {
+            double tot_injected[2] = { 0.0 };
+            double tot_produced[2] = { 0.0 };
+            Opm::Watercut watercut;
+            watercut.push(0.0, 0.0, 0.0);
 #if 0
-        // These must be changed for three-phase.
+            std::vector<double> fractional_flows;
-        double init_surfvol[2] = { 0.0 };
+            std::vector<double> well_resflows_phase;
-        double inplace_surfvol[2] = { 0.0 };
+            if (wells_) {
-        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
+                well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
-        Opm::WellReport wellreport;
-#endif
-        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);
-#if 0
-            wellreport.push(props_, *wells_,
-                            state.pressure(), state.surfacevol(), state.saturation(),
-                            0.0, well_state.bhp(), well_state.perfRates());
-#endif
-        }
-        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);
-        }
-            // Report timestep and (optionally) write state to disk.
-        step_timer.start();
-        timer.report(std::cout);
-        if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
-            if (output_vtk_) {
-                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
-            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            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.
 
-        SimulatorReport sreport;
+            step_timer.start();
+            timer.report(std::cout);
 
-        bool well_control_passed = !check_well_controls_;
+            WellsManager wells_manager(eclipse_state_,
-        int well_control_iteration = 0;
+                                       timer.currentStepNum(),
-        do {
+                                       Opm::UgGridHelpers::numCells(grid_),
+                                       Opm::UgGridHelpers::globalCell(grid_),
+                                       Opm::UgGridHelpers::cartDims(grid_),
+                                       Opm::UgGridHelpers::dimensions(grid_),
+                                       Opm::UgGridHelpers::beginCellCentroids(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());
+            if (timer.currentStepNum() != 0) {
+                // Transfer previous well state to current.
+                well_state.partialCopy(prev_well_state, *wells, prev_well_state.numWells());
+            }
+            //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_, *wells, Opm::UgGridHelpers::numCells(grid_)));
+            } 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(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.basicWellState());
+            }
             // Run solver.
-            const double current_time = timer.simulationTimeElapsed();
-            double stepsize = timer.currentStepLength();
-            polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
             solver_timer.start();
-            std::vector<double> initial_pressure = state.pressure();
+            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, transport_src);
+            solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
             // 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;
-            sreport.pressure_time = st;
+            // Update pore volumes if rock is compressible.
-
+            if (rock_comp_props_ && rock_comp_props_->isActive()) {
-            // Optionally, check if well controls are satisfied.
+                initial_porevol = porevol;
-            if (check_well_controls_) {
+                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
-                Opm::computePhaseFlowRatesPerWell(*wells_,
-                                                  well_state.perfRates(),
-                                                  fractional_flows,
-                                                  well_resflows_phase);
-                std::cout << "Checking well conditions." << std::endl;
-                // For testing we set surface := reservoir
-                well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
-                ++well_control_iteration;
-                if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
-                    OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
-                }
-                if (!well_control_passed) {
-                    std::cout << "Well controls not passed, solving again." << std::endl;
-                } else {
-                    std::cout << "Well conditions met." << std::endl;
-                }
             }
-        } while (!well_control_passed);
+
-        // Update pore volumes if rock is compressible.
+            double injected[2] = { 0.0 };
-        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            double produced[2] = { 0.0 };
-            initial_porevol = porevol;
+    		double polyinj = 0;
-            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+    		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);
+                outputWaterCut(watercut, output_dir_);
+            }
+            ++timer;
+            prev_well_state = well_state;
         }
-
+        // Write final simulation state.
-        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;
-        sreport.total_time =  step_timer.secsSinceStart();
         if (output_) {
-            sreport.reportParam(tstep_os);
-
             if (output_vtk_) {
                 outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
             outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
-            outputWaterCut(watercut, output_dir_);
+            output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState());
-            tstep_os.close();
         }
 
         total_timer.stop();
-
         SimulatorReport report;
         report.pressure_time = stime;
         report.transport_time = 0.0;

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp

@@ -18,11 +18,12 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
-#define OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
 
-#include <boost/shared_ptr.hpp>
+#include <memory>
 #include <vector>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
 
 struct UnstructuredGrid;
 struct Wells;
@@ -33,11 +34,13 @@ namespace Opm
     class BlackoilPropsAdInterface;
     class RockCompressibility;
     class DerivedGeology;
+    class WellStateFullyImplicitBlackoil;
     class WellsManager;
+    class EclipseWriter;
+    class EclipseState;
     class NewtonIterationBlackoilInterface;
     class SimulatorTimer;
     class PolymerBlackoilState;
-    class WellState;
     class PolymerPropsAd;
     class PolymerInflowInterface;
     struct SimulatorReport;
@@ -66,8 +69,9 @@ namespace Opm
         /// \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] well_manager  well manager, may manage no (null) wells
+        /// \param[in] eclipse_state  
-        /// \param[in] polymer_inflow polymer influx.
+        /// \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,
@@ -76,8 +80,9 @@ namespace Opm
                                    				  const BlackoilPropsAdInterface& props,
                                        			  const PolymerPropsAd&    polymer_props,
                                        			  const RockCompressibility* rock_comp_props,
-                                       			  WellsManager& wells_manager,
+                                                  std::shared_ptr<EclipseState> eclipse_state,
-                                       			  PolymerInflowInterface& polymer_inflow,
+                                                  EclipseWriter& eclipse_writer,
+                                                  Opm::DeckConstPtr& deck,
                                        			  NewtonIterationBlackoilInterface& linsolver,
                                        			  const double* gravity);
 
@@ -86,18 +91,16 @@ namespace Opm
         /// 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,
+                            PolymerBlackoilState& state);
-                            WellState& well_state);
 
     private:
         class Impl;
         // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
-        boost::shared_ptr<Impl> pimpl_;
+        std::shared_ptr<Impl> pimpl_;
     };
 
 } // namespace Opm
 
-#endif // OPM_SIMULATORFULLYIMPLICITBLACKOIL_HEADER_INCLUDED
+#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED