Merge pull request #73 from qilicun/bo_polymer

Fully implicit Two phase polymer solver and blackoil polymer solver

.gitignore

@@ -68,4 +68,10 @@ examples/upscale_*
 examples/grdecldips
 data
 steadystate_test_implicit
-steadystate_test_explicit
+steadystate_test_explicit
+=======
+*.iml
+.idea
+*/opm-simulation
+
+

examples/sim_poly_fi2p_comp_ad.cpp

@@ -0,0 +1,225 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 Statoil ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#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>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+
+#include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+#include <opm/polymer/PolymerState.hpp>
+
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/GeoProps.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <memory>
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for fully implicit three-phase black-oil flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    if (!use_deck) {
+        OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
+                  "Specify the deck with deck_filename=deckname.data (for example).");
+    }
+    std::shared_ptr<GridManager> grid;
+    std::shared_ptr<BlackoilPropertiesInterface> props;
+    std::shared_ptr<BlackoilPropsAdInterface> new_props;
+    std::shared_ptr<RockCompressibility> rock_comp;
+    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");
+
+    Opm::ParserPtr newParser(new Opm::Parser());
+    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));
+    auto &cGrid = *grid->c_grid();
+    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
+    Opm::EclipseWriter outputWriter(param,
+                                    eclipseState,
+                                    pu,
+                                    cGrid.number_of_cells,
+                                    cGrid.global_cell);
+    // Rock and fluid init
+    props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
+    new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
+    PolymerProperties polymer_props(deck, eclipseState);
+    PolymerPropsAd polymer_props_ad(polymer_props);
+    // Rock compressibility.
+    rock_comp.reset(new RockCompressibility(deck, eclipseState));
+    // Gravity.
+    gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+    // Init state variables (saturation and pressure).
+    if (param.has("init_saturation")) {
+        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        initBlackoilSurfvol(*grid->c_grid(), *props, state);
+    } else {
+        initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+    }
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+    // Solver for Newton iterations.
+    std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
+    if (param.getDefault("use_cpr", true)) {
+        fis_solver.reset(new NewtonIterationBlackoilCPR(param));
+    } else {
+        fis_solver.reset(new NewtonIterationBlackoilSimple(param));
+    }
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
+    SimulatorTimer simtimer;
+    simtimer.init(timeMap);
+
+
+    SimulatorReport rep;
+    // With a deck, we may have more epochs etc.
+    WellState well_state;
+    // Check for WPOLYMER presence in last epoch to decide
+    // polymer injection control type.
+    const bool use_wpolymer = deck->hasKeyword("WPOLYMER");
+    if (use_wpolymer) {
+        if (param.has("poly_start_days")) {
+            OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck. "
+                        "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
+        }
+    }
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n"
+              << std::flush;
+    SimulatorReport fullReport;
+    // Create and run simulator.
+    Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
+    SimulatorFullyImplicitCompressiblePolymer simulator(param,
+                                             *grid->c_grid(),
+                                             geology,
+                                             *new_props,
+                                             polymer_props_ad,
+                                             rock_comp->isActive() ? rock_comp.get() : 0,
+                                             eclipseState,
+                                             outputWriter,
+                                             deck,
+                                             *fis_solver,
+                                             grav);
+    fullReport= simulator.run(simtimer, state);
+
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    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);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
+    }
+
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}
+

examples/sim_poly_fibo_ad.cpp

@@ -0,0 +1,253 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#include "config.h"
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/initStateEquil.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/thresholdPressures.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>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilSimple.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilCPR.hpp>
+
+#include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+#include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/algorithm/string.hpp>
+
+#include <memory>
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for fully implicit three-phase black-oil flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    if (!use_deck) {
+        OPM_THROW(std::runtime_error, "This program must be run with an input deck. "
+                  "Specify the deck with deck_filename=deckname.data (for example).");
+    }
+    std::shared_ptr<GridManager> grid;
+    std::shared_ptr<BlackoilPropertiesInterface> props;
+    std::shared_ptr<BlackoilPropsAdInterface> new_props;
+    std::shared_ptr<RockCompressibility> rock_comp;
+    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");
+
+    Opm::ParserPtr newParser(new Opm::Parser() );
+    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));
+    auto &cGrid = *grid->c_grid();
+    const PhaseUsage pu = Opm::phaseUsageFromDeck(deck);
+    Opm::EclipseWriter outputWriter(param,
+                                    eclipseState,
+                                    pu,
+                                    cGrid.number_of_cells,
+                                    cGrid.global_cell);
+
+    // Rock and fluid init
+    props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param));
+    new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid()));
+    const bool polymer = deck->hasKeyword("POLYMER");
+    const bool use_wpolymer = deck->hasKeyword("WPOLYMER");
+    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.
+    gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity;
+
+    // Init state variables (saturation and pressure).
+    if (param.has("init_saturation")) {
+        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        initBlackoilSurfvol(*grid->c_grid(), *props, state);
+        enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
+        if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
+            const int np = props->numPhases();
+            const int nc = grid->c_grid()->number_of_cells;
+            for (int c = 0; c < nc; ++c) {
+                state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]]
+                    / state.surfacevol()[c*np + pu.phase_pos[Oil]];
+            }
+        }
+    } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) {
+        state.init(*grid->c_grid(), props->numPhases());
+        const double grav = param.getDefault("gravity", unit::gravity);
+        initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState());
+        state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0);
+    } else {
+        state.init(*grid->c_grid(), props->numPhases());
+        initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState());
+    }
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Solver for Newton iterations.
+    std::unique_ptr<NewtonIterationBlackoilInterface> fis_solver;
+    if (param.getDefault("use_cpr", true)) {
+        fis_solver.reset(new NewtonIterationBlackoilCPR(param));
+    } else {
+        fis_solver.reset(new NewtonIterationBlackoilSimple(param));
+    }
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::string output_dir;
+    if (output) {
+        // Create output directory if needed.
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        // Write simulation parameters.
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+    Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap());
+    SimulatorTimer simtimer;
+
+    // initialize variables
+    simtimer.init(timeMap);
+    if (polymer){
+        if (!use_wpolymer) {
+            OPM_MESSAGE("Warning: simulate polymer injection without WPOLYMER.");
+        } else {
+            if (param.has("polymer_start_days")) {
+                OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck."
+                            "You seem to be trying to control it via parameter poly_start_days (etc.) as well.");
+            }
+        }
+    } else {
+        if (use_wpolymer) {
+            OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario.");
+        }
+    }
+    std::cout << "\n\n================ Starting main simulation loop ===============\n"
+              << std::flush;
+    SimulatorReport fullReport;
+    Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav);
+
+    std::vector<double> threshold_pressures = thresholdPressures(deck, eclipseState, *grid->c_grid());
+    SimulatorFullyImplicitBlackoilPolymer<UnstructuredGrid> simulator(param,
+                                             *grid->c_grid(),
+                                             geology,
+                                             *new_props,
+                                             polymer_props_ad,
+                                             rock_comp->isActive() ? rock_comp.get() : 0,
+                                             *fis_solver,
+                                             grav,
+                                             deck->hasKeyword("DISGAS"),
+                                             deck->hasKeyword("VAPOIL"),
+                                             polymer,
+                                             eclipseState,
+                                             outputWriter,
+                                             deck,
+                                             threshold_pressures);
+
+
+    fullReport = simulator.run(simtimer, state);
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    fullReport.report(std::cout);
+
+    if (output) {
+        std::string filename = output_dir + "/walltime.txt";
+        std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+        fullReport.reportParam(tot_os);
+        warnIfUnusedParams(param);
+    }
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}
+

opm/polymer/PolymerBlackoilState.hpp

@@ -44,6 +44,10 @@ namespace Opm
             concentration_.resize(number_of_cells, 0.0);
             cmax_.resize(number_of_cells, 0.0);
         }
+        int numPhases() const
+        {
+            return state_blackoil_.numPhases();
+        }
 
         enum ExtremalSat { MinSat = BlackoilState::MinSat, MaxSat = BlackoilState::MaxSat };
 
@@ -61,6 +65,8 @@ namespace Opm
         std::vector<double>& facepressure()     { return state_blackoil_.facepressure(); }
         std::vector<double>& faceflux    ()     { return state_blackoil_.faceflux(); }
         std::vector<double>& saturation  ()     { return state_blackoil_.saturation(); }
+        std::vector<double>& gasoilratio ()     { return state_blackoil_.gasoilratio(); }
+        std::vector<double>& rv          ()     { return state_blackoil_.rv(); }
         std::vector<double>& concentration()    { return concentration_; }
         std::vector<double>& maxconcentration() { return cmax_; }
 
@@ -69,6 +75,8 @@ namespace Opm
         const std::vector<double>& facepressure() const     { return state_blackoil_.facepressure(); }
         const std::vector<double>& faceflux    () const     { return state_blackoil_.faceflux(); }
         const std::vector<double>& saturation  () const     { return state_blackoil_.saturation(); }
+        const std::vector<double>& gasoilratio() const     { return state_blackoil_.gasoilratio(); }
+        const std::vector<double>& rv          () const     { return state_blackoil_.rv(); }
         const std::vector<double>& concentration() const    { return concentration_; }
         const std::vector<double>& maxconcentration() const { return cmax_; }
 

opm/polymer/PolymerState.hpp

@@ -55,6 +55,10 @@ namespace Opm
             state2p_.setFirstSat(cells, props, static_cast<TwophaseState::ExtremalSat>(es));
         }
 
+        inline int numPhases() const
+        {
+            return state2p_.numPhases();
+        }
         std::vector<double>& pressure    ()     { return state2p_.pressure(); }
         std::vector<double>& facepressure()     { return state2p_.facepressure(); }
         std::vector<double>& faceflux    ()     { return state2p_.faceflux(); }
@@ -71,7 +75,7 @@ namespace Opm
 
         TwophaseState& twophaseState() { return state2p_; }
         const TwophaseState& twophaseState() const { return state2p_; }
-
+        
     private:
         TwophaseState state2p_;
         std::vector<double> concentration_;

opm/polymer/SimulatorCompressiblePolymer.cpp

@@ -292,207 +292,206 @@ namespace Opm
             wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
                             state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
         }
-        for (; !timer.done(); ++timer) {
-            // Report timestep and (optionally) write state to disk.
-            timer.report(std::cout);
-            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
-                if (output_vtk_) {
-                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+        // Report timestep and (optionally) write state to disk.
+        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_);
+        }
+
+        initial_pressure = state.pressure();
+
+        // Solve pressure equation.
+        if (check_well_controls_) {
+            computeFractionalFlow(props_, poly_props_, allcells_,
+                                  state.pressure(), state.surfacevol(), state.saturation(),
+                                  state.concentration(), state.maxconcentration(),
+                                  fractional_flows);
+            wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
+        }
+        bool well_control_passed = !check_well_controls_;
+        int well_control_iteration = 0;
+        do {
+            // Run solver
+            pressure_timer.start();
+            psolver_.solve(timer.currentStepLength(), state, well_state);
+
+            // Renormalize pressure if both fluids and rock are
+            // incompressible, and there are no pressure
+            // conditions (bcs or wells).  It is deemed sufficient
+            // for now to renormalize using geometric volume
+            // instead of pore volume.
+            if (psolver_.singularPressure()) {
+                // Compute average pressures of previous and last
+                // step, and total volume.
+                double av_prev_press = 0.0;
+                double av_press = 0.0;
+                double tot_vol = 0.0;
+                const int num_cells = grid_.number_of_cells;
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                    av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                    tot_vol       += grid_.cell_volumes[cell];
+                }
+                // Renormalization constant
+                const double ren_const = (av_prev_press - av_press)/tot_vol;
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    state.pressure()[cell] += ren_const;
+                }
+                const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                for (int well = 0; well < num_wells; ++well) {
+                    well_state.bhp()[well] += ren_const;
                 }
-                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
             }
 
-            initial_pressure = state.pressure();
+            // Stop timer and report
+            pressure_timer.stop();
+            double pt = pressure_timer.secsSinceStart();
+            std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
+            ptime += pt;
 
-            // Solve pressure equation.
+            // Optionally, check if well controls are satisfied.
             if (check_well_controls_) {
-                computeFractionalFlow(props_, poly_props_, allcells_,
-                                      state.pressure(), state.surfacevol(), state.saturation(),
-                                      state.concentration(), state.maxconcentration(),
-                                      fractional_flows);
-                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
-            }
-            bool well_control_passed = !check_well_controls_;
-            int well_control_iteration = 0;
-            do {
-                // Run solver
-                pressure_timer.start();
-                psolver_.solve(timer.currentStepLength(), state, well_state);
-
-                // Renormalize pressure if both fluids and rock are
-                // incompressible, and there are no pressure
-                // conditions (bcs or wells).  It is deemed sufficient
-                // for now to renormalize using geometric volume
-                // instead of pore volume.
-                if (psolver_.singularPressure()) {
-                    // Compute average pressures of previous and last
-                    // step, and total volume.
-                    double av_prev_press = 0.0;
-                    double av_press = 0.0;
-                    double tot_vol = 0.0;
-                    const int num_cells = grid_.number_of_cells;
-                    for (int cell = 0; cell < num_cells; ++cell) {
-                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
-                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
-                        tot_vol       += grid_.cell_volumes[cell];
-                    }
-                    // Renormalization constant
-                    const double ren_const = (av_prev_press - av_press)/tot_vol;
-                    for (int cell = 0; cell < num_cells; ++cell) {
-                        state.pressure()[cell] += ren_const;
-                    }
-                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
-                    for (int well = 0; well < num_wells; ++well) {
-                        well_state.bhp()[well] += ren_const;
-                    }
+                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.");
                 }
-
-                // Stop timer and report
-                pressure_timer.stop();
-                double pt = pressure_timer.secsSinceStart();
-                std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
-                ptime += pt;
-
-                // Optionally, check if well controls are satisfied.
-                if (check_well_controls_) {
-                    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.
-            if (rock_comp_props_ && rock_comp_props_->isActive()) {
-                initial_porevol = porevol;
-                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
-            }
-
-            // Process transport sources (to include bdy terms and well flows).
-            Opm::computeTransportSource(props_, wells_, well_state, transport_src);
-
-            // Find inflow rate.
-            const double current_time = timer.simulationTimeElapsed();
-            double stepsize = timer.currentStepLength();
-            polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
-
-            // Solve transport.
-            transport_timer.start();
-            if (num_transport_substeps_ != 1) {
-                stepsize /= double(num_transport_substeps_);
-                std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
-            }
-            double injected[2] = { 0.0 };
-            double produced[2] = { 0.0 };
-            double polyinj = 0.0;
-            double polyprod = 0.0;
-            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_.solve(&state.faceflux()[0], initial_pressure,
-                               state.pressure(), &initial_porevol[0], &porevol[0],
-                               &transport_src[0], &polymer_inflow_c[0], stepsize,
-                               state.saturation(), state.surfacevol(),
-                               state.concentration(), state.maxconcentration());
-                double substep_injected[2] = { 0.0 };
-                double substep_produced[2] = { 0.0 };
-                double substep_polyinj = 0.0;
-                double substep_polyprod = 0.0;
-                Opm::computeInjectedProduced(props_, poly_props_,
-                                             state,
-                                             transport_src, polymer_inflow_c, stepsize,
-                                             substep_injected, substep_produced,
-                                             substep_polyinj, substep_polyprod);
-                injected[0] += substep_injected[0];
-                injected[1] += substep_injected[1];
-                produced[0] += substep_produced[0];
-                produced[1] += substep_produced[1];
-                polyinj += substep_polyinj;
-                polyprod += substep_polyprod;
-                if (gravity_ != 0 && use_segregation_split_) {
-                    tsolver_.solveGravity(columns_, stepsize,
-                                          state.saturation(), state.surfacevol(),
-                                          state.concentration(), state.maxconcentration());
+                if (!well_control_passed) {
+                    std::cout << "Well controls not passed, solving again." << std::endl;
+                } else {
+                    std::cout << "Well conditions met." << std::endl;
                 }
             }
-            transport_timer.stop();
-            double tt = transport_timer.secsSinceStart();
-            std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
-            ttime += tt;
+        } while (!well_control_passed);
 
-            // Report volume balances.
-            Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
-            polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
-            polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
-                                                            state, rock_comp_props_);
-            tot_injected[0] += injected[0];
-            tot_injected[1] += injected[1];
-            tot_produced[0] += produced[0];
-            tot_produced[1] += produced[1];
-            tot_polyinj += polyinj;
-            tot_polyprod += polyprod;
-            std::cout.precision(5);
-            const int width = 18;
-            std::cout << "\nMass balance:        "
-                "                   water(surfvol)      oil(surfvol)       polymer(kg)\n";
-            std::cout << "    In-place:                       "
-                      << std::setw(width) << inplace_surfvol[0]
-                      << std::setw(width) << inplace_surfvol[1]
-                      << std::setw(width) << polymass << std::endl;
-            std::cout << "    Adsorbed:                       "
-                      << std::setw(width) << 0.0
-                      << std::setw(width) << 0.0
-                      << std::setw(width) << polymass_adsorbed << std::endl;
-            std::cout << "    Injected:                       "
-                      << std::setw(width) << injected[0]
-                      << std::setw(width) << injected[1]
-                      << std::setw(width) << polyinj << std::endl;
-            std::cout << "    Produced:                       "
-                      << std::setw(width) << produced[0]
-                      << std::setw(width) << produced[1]
-                      << std::setw(width) << polyprod << std::endl;
-            std::cout << "    Total inj:                      "
-                      << std::setw(width) << tot_injected[0]
-                      << std::setw(width) << tot_injected[1]
-                      << std::setw(width) << tot_polyinj << std::endl;
-            std::cout << "    Total prod:                     "
-                      << std::setw(width) << tot_produced[0]
-                      << std::setw(width) << tot_produced[1]
-                      << std::setw(width) << tot_polyprod << std::endl;
-            const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
-                                        init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1],
-                                        init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed };
-            std::cout << "    Initial - inplace + inj - prod: "
-                      << std::setw(width) << balance[0]
-                      << std::setw(width) << balance[1]
-                      << std::setw(width) << balance[2]
-                      << std::endl;
-            std::cout << "    Relative mass error:            "
-                      << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
-                      << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
-                      << std::setw(width) << balance[2]/(init_polymass + tot_polyinj)
-                      << std::endl;
-            std::cout.precision(8);
+        // Update pore volumes if rock is compressible.
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            initial_porevol = porevol;
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+        }
 
-            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
-                          produced[0]/(produced[0] + produced[1]),
-                          tot_produced[0]/tot_porevol_init);
-            if (wells_) {
-            wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
-                            state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(),
-                            well_state.bhp(), well_state.perfRates());
+        // Process transport sources (to include bdy terms and well flows).
+        Opm::computeTransportSource(props_, wells_, well_state, transport_src);
+
+        // Find inflow rate.
+        const double current_time = timer.simulationTimeElapsed();
+        double stepsize = timer.currentStepLength();
+        polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
+
+
+        // Solve transport.
+        transport_timer.start();
+        if (num_transport_substeps_ != 1) {
+            stepsize /= double(num_transport_substeps_);
+            std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
+        }
+        double injected[2] = { 0.0 };
+        double produced[2] = { 0.0 };
+        double polyinj = 0.0;
+        double polyprod = 0.0;
+        for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
+            tsolver_.solve(&state.faceflux()[0], initial_pressure,
+                           state.pressure(), &initial_porevol[0], &porevol[0],
+                           &transport_src[0], &polymer_inflow_c[0], stepsize,
+                           state.saturation(), state.surfacevol(),
+                           state.concentration(), state.maxconcentration());
+            double substep_injected[2] = { 0.0 };
+            double substep_produced[2] = { 0.0 };
+            double substep_polyinj = 0.0;
+            double substep_polyprod = 0.0;
+            Opm::computeInjectedProduced(props_, poly_props_,
+                                         state,
+                                         transport_src, polymer_inflow_c, stepsize,
+                                         substep_injected, substep_produced,
+                                         substep_polyinj, substep_polyprod);
+            injected[0] += substep_injected[0];
+            injected[1] += substep_injected[1];
+            produced[0] += substep_produced[0];
+            produced[1] += substep_produced[1];
+            polyinj += substep_polyinj;
+            polyprod += substep_polyprod;
+            if (gravity_ != 0 && use_segregation_split_) {
+                tsolver_.solveGravity(columns_, stepsize,
+                                      state.saturation(), state.surfacevol(),
+                                      state.concentration(), state.maxconcentration());
             }
         }
+        transport_timer.stop();
+        double tt = transport_timer.secsSinceStart();
+        std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+        ttime += tt;
+
+        // Report volume balances.
+        Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
+                                                        state, rock_comp_props_);
+        tot_injected[0] += injected[0];
+        tot_injected[1] += injected[1];
+        tot_produced[0] += produced[0];
+        tot_produced[1] += produced[1];
+        tot_polyinj += polyinj;
+        tot_polyprod += polyprod;
+        std::cout.precision(5);
+        const int width = 18;
+        std::cout << "\nMass balance:        "
+            "                   water(surfvol)      oil(surfvol)       polymer(kg)\n";
+        std::cout << "    In-place:                       "
+                  << std::setw(width) << inplace_surfvol[0]
+                  << std::setw(width) << inplace_surfvol[1]
+                  << std::setw(width) << polymass << std::endl;
+        std::cout << "    Adsorbed:                       "
+                  << std::setw(width) << 0.0
+                  << std::setw(width) << 0.0
+                  << std::setw(width) << polymass_adsorbed << std::endl;
+        std::cout << "    Injected:                       "
+                  << std::setw(width) << injected[0]
+                  << std::setw(width) << injected[1]
+                  << std::setw(width) << polyinj << std::endl;
+        std::cout << "    Produced:                       "
+                  << std::setw(width) << produced[0]
+                  << std::setw(width) << produced[1]
+                  << std::setw(width) << polyprod << std::endl;
+        std::cout << "    Total inj:                      "
+                  << std::setw(width) << tot_injected[0]
+                  << std::setw(width) << tot_injected[1]
+                  << std::setw(width) << tot_polyinj << std::endl;
+        std::cout << "    Total prod:                     "
+                  << std::setw(width) << tot_produced[0]
+                  << std::setw(width) << tot_produced[1]
+                  << std::setw(width) << tot_polyprod << std::endl;
+        const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
+                                    init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1],
+                                    init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed };
+        std::cout << "    Initial - inplace + inj - prod: "
+                  << std::setw(width) << balance[0]
+                  << std::setw(width) << balance[1]
+                  << std::setw(width) << balance[2]
+                  << std::endl;
+        std::cout << "    Relative mass error:            "
+                  << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
+                  << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
+                  << std::setw(width) << balance[2]/(init_polymass + tot_polyinj)
+                  << std::endl;
+        std::cout.precision(8);
+
+        watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
+                      produced[0]/(produced[0] + produced[1]),
+                      tot_produced[0]/tot_porevol_init);
+        if (wells_) {
+        wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
+                        state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(),
+                        well_state.bhp(), well_state.perfRates());
+        }
 
         if (output_) {
             if (output_vtk_) {

opm/polymer/SimulatorPolymer.cpp

@@ -315,196 +315,194 @@ namespace Opm
             well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
             wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
         }
-        for (; !timer.done(); ++timer) {
-            // Report timestep and (optionally) write state to disk.
-            timer.report(std::cout);
-            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
-                if (output_vtk_) {
-                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+        // Report timestep and (optionally) write state to disk.
+        timer.report(std::cout);
+        if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+            if (output_binary_) {
+                outputStateBinary(grid_, state, timer, output_dir_);
+            }
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+        }
+
+        // Solve pressure.
+        if (check_well_controls_) {
+            computeFractionalFlow(props_, poly_props_, allcells_,
+                                  state.saturation(), state.concentration(), state.maxconcentration(),
+                                  fractional_flows);
+            wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
+        }
+        bool well_control_passed = !check_well_controls_;
+        int well_control_iteration = 0;
+        do {
+            // Run solver.
+            pressure_timer.start();
+            std::vector<double> initial_pressure = state.pressure();
+            psolver_.solve(timer.currentStepLength(), state, well_state);
+
+            // Renormalize pressure if rock is incompressible, and
+            // there are no pressure conditions (bcs or wells).
+            // It is deemed sufficient for now to renormalize
+            // using geometric volume instead of pore volume.
+            if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
+                && allNeumannBCs(bcs_) && allRateWells(wells_)) {
+                // Compute average pressures of previous and last
+                // step, and total volume.
+                double av_prev_press = 0.0;
+                double av_press = 0.0;
+                double tot_vol = 0.0;
+                const int num_cells = grid_.number_of_cells;
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                    av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                    tot_vol       += grid_.cell_volumes[cell];
                 }
-                if (output_binary_) {
-                    outputStateBinary(grid_, state, timer, output_dir_);
+                // Renormalization constant
+                const double ren_const = (av_prev_press - av_press)/tot_vol;
+                for (int cell = 0; cell < num_cells; ++cell) {
+                    state.pressure()[cell] += ren_const;
+                }
+                const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                for (int well = 0; well < num_wells; ++well) {
+                    well_state.bhp()[well] += ren_const;
                 }
-                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
             }
 
-            // Solve pressure.
+            // Stop timer and report.
+            pressure_timer.stop();
+            double pt = pressure_timer.secsSinceStart();
+            std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
+            ptime += pt;
+
+            // Optionally, check if well controls are satisfied.
             if (check_well_controls_) {
-                computeFractionalFlow(props_, poly_props_, allcells_,
-                                      state.saturation(), state.concentration(), state.maxconcentration(),
-                                      fractional_flows);
-                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
-            }
-            bool well_control_passed = !check_well_controls_;
-            int well_control_iteration = 0;
-            do {
-                // Run solver.
-                pressure_timer.start();
-                std::vector<double> initial_pressure = state.pressure();
-                psolver_.solve(timer.currentStepLength(), state, well_state);
-
-                // Renormalize pressure if rock is incompressible, and
-                // there are no pressure conditions (bcs or wells).
-                // It is deemed sufficient for now to renormalize
-                // using geometric volume instead of pore volume.
-                if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
-                    && allNeumannBCs(bcs_) && allRateWells(wells_)) {
-                    // Compute average pressures of previous and last
-                    // step, and total volume.
-                    double av_prev_press = 0.0;
-                    double av_press = 0.0;
-                    double tot_vol = 0.0;
-                    const int num_cells = grid_.number_of_cells;
-                    for (int cell = 0; cell < num_cells; ++cell) {
-                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
-                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
-                        tot_vol       += grid_.cell_volumes[cell];
-                    }
-                    // Renormalization constant
-                    const double ren_const = (av_prev_press - av_press)/tot_vol;
-                    for (int cell = 0; cell < num_cells; ++cell) {
-                        state.pressure()[cell] += ren_const;
-                    }
-                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
-                    for (int well = 0; well < num_wells; ++well) {
-                        well_state.bhp()[well] += ren_const;
-                    }
+                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.");
                 }
-
-                // Stop timer and report.
-                pressure_timer.stop();
-                double pt = pressure_timer.secsSinceStart();
-                std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
-                ptime += pt;
-
-                // Optionally, check if well controls are satisfied.
-                if (check_well_controls_) {
-                    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.
-            if (rock_comp_props_ && rock_comp_props_->isActive()) {
-                initial_porevol = porevol;
-                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
-            }
-
-            // Process transport sources (to include bdy terms and well flows).
-            Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
-                                        wells_, well_state.perfRates(), transport_src);
-
-            // Find inflow rate.
-            const double current_time = timer.simulationTimeElapsed();
-            double stepsize = timer.currentStepLength();
-            polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
-
-            // Solve transport.
-            transport_timer.start();
-            if (num_transport_substeps_ != 1) {
-                stepsize /= double(num_transport_substeps_);
-                std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
-            }
-            double substep_injected[2] = { 0.0 };
-            double substep_produced[2] = { 0.0 };
-            double substep_polyinj = 0.0;
-            double substep_polyprod = 0.0;
-            injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
-            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
-                               state.saturation(), state.concentration(), state.maxconcentration());
-                Opm::computeInjectedProduced(props_, poly_props_,
-                                             state,
-                                             transport_src, polymer_inflow_c, stepsize,
-                                             substep_injected, substep_produced, substep_polyinj, substep_polyprod);
-                injected[0] += substep_injected[0];
-                injected[1] += substep_injected[1];
-                produced[0] += substep_produced[0];
-                produced[1] += substep_produced[1];
-                polyinj += substep_polyinj;
-                polyprod += substep_polyprod;
-                if (use_segregation_split_) {
-                    tsolver_.solveGravity(columns_, &porevol[0], stepsize,
-                                          state.saturation(), state.concentration(), state.maxconcentration());
+                if (!well_control_passed) {
+                    std::cout << "Well controls not passed, solving again." << std::endl;
+                } else {
+                    std::cout << "Well conditions met." << std::endl;
                 }
             }
-            transport_timer.stop();
-            double tt = transport_timer.secsSinceStart();
-            std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
-            ttime += tt;
+        } while (!well_control_passed);
 
-            // Report volume balances.
-            Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
-            polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
-            polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
-            tot_injected[0] += injected[0];
-            tot_injected[1] += injected[1];
-            tot_produced[0] += produced[0];
-            tot_produced[1] += produced[1];
-            tot_polyinj += polyinj;
-            tot_polyprod += polyprod;
-            std::cout.precision(5);
-            const int width = 18;
-            std::cout << "\nVolume and polymer mass balance: "
-                "   water(pv)           oil(pv)       polymer(kg)\n";
-            std::cout << "    Saturated volumes:     "
-                      << std::setw(width) << satvol[0]/tot_porevol_init
-                      << std::setw(width) << satvol[1]/tot_porevol_init
-                      << std::setw(width) << polymass << std::endl;
-            std::cout << "    Adsorbed volumes:      "
-                      << std::setw(width) << 0.0
-                      << std::setw(width) << 0.0
-                      << std::setw(width) << polymass_adsorbed << std::endl;
-            std::cout << "    Injected volumes:      "
-                      << std::setw(width) << injected[0]/tot_porevol_init
-                      << std::setw(width) << injected[1]/tot_porevol_init
-                      << std::setw(width) << polyinj << std::endl;
-            std::cout << "    Produced volumes:      "
-                      << std::setw(width) << produced[0]/tot_porevol_init
-                      << std::setw(width) << produced[1]/tot_porevol_init
-                      << std::setw(width) << polyprod << std::endl;
-            std::cout << "    Total inj volumes:     "
-                      << std::setw(width) << tot_injected[0]/tot_porevol_init
-                      << std::setw(width) << tot_injected[1]/tot_porevol_init
-                      << std::setw(width) << tot_polyinj << std::endl;
-            std::cout << "    Total prod volumes:    "
-                      << std::setw(width) << tot_produced[0]/tot_porevol_init
-                      << std::setw(width) << tot_produced[1]/tot_porevol_init
-                      << std::setw(width) << tot_polyprod << std::endl;
-            std::cout << "    In-place + prod - inj: "
-                      << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
-                      << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init
-                      << std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
-            std::cout << "    Init - now - pr + inj: "
-                      << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
-                      << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
-                      << std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
-                      << std::endl;
-            std::cout.precision(8);
+        // Update pore volumes if rock is compressible.
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            initial_porevol = porevol;
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+        }
 
-            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
-                          produced[0]/(produced[0] + produced[1]),
-                          tot_produced[0]/tot_porevol_init);
-            if (wells_) {
-                wellreport.push(props_, *wells_, state.saturation(),
-                                timer.simulationTimeElapsed() + timer.currentStepLength(),
-                                well_state.bhp(), well_state.perfRates());
+        // Process transport sources (to include bdy terms and well flows).
+        Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
+                                    wells_, well_state.perfRates(), transport_src);
+
+        // Find inflow rate.
+        const double current_time = timer.simulationTimeElapsed();
+        double stepsize = timer.currentStepLength();
+        polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
+
+        // Solve transport.
+        transport_timer.start();
+        if (num_transport_substeps_ != 1) {
+            stepsize /= double(num_transport_substeps_);
+            std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
+        }
+        double substep_injected[2] = { 0.0 };
+        double substep_produced[2] = { 0.0 };
+        double substep_polyinj = 0.0;
+        double substep_polyprod = 0.0;
+        injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
+        for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
+            tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
+                           state.saturation(), state.concentration(), state.maxconcentration());
+            Opm::computeInjectedProduced(props_, poly_props_,
+                                         state,
+                                         transport_src, polymer_inflow_c, stepsize,
+                                         substep_injected, substep_produced, substep_polyinj, substep_polyprod);
+            injected[0] += substep_injected[0];
+            injected[1] += substep_injected[1];
+            produced[0] += substep_produced[0];
+            produced[1] += substep_produced[1];
+            polyinj += substep_polyinj;
+            polyprod += substep_polyprod;
+            if (use_segregation_split_) {
+                tsolver_.solveGravity(columns_, &porevol[0], stepsize,
+                                      state.saturation(), state.concentration(), state.maxconcentration());
             }
         }
+        transport_timer.stop();
+        double tt = transport_timer.secsSinceStart();
+        std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+        ttime += tt;
+
+        // Report volume balances.
+        Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
+        tot_injected[0] += injected[0];
+        tot_injected[1] += injected[1];
+        tot_produced[0] += produced[0];
+        tot_produced[1] += produced[1];
+        tot_polyinj += polyinj;
+        tot_polyprod += polyprod;
+        std::cout.precision(5);
+        const int width = 18;
+        std::cout << "\nVolume and polymer mass balance: "
+            "   water(pv)           oil(pv)       polymer(kg)\n";
+        std::cout << "    Saturated volumes:     "
+                  << std::setw(width) << satvol[0]/tot_porevol_init
+                  << std::setw(width) << satvol[1]/tot_porevol_init
+                  << std::setw(width) << polymass << std::endl;
+        std::cout << "    Adsorbed volumes:      "
+                  << std::setw(width) << 0.0
+                  << std::setw(width) << 0.0
+                  << std::setw(width) << polymass_adsorbed << std::endl;
+        std::cout << "    Injected volumes:      "
+                  << std::setw(width) << injected[0]/tot_porevol_init
+                  << std::setw(width) << injected[1]/tot_porevol_init
+                  << std::setw(width) << polyinj << std::endl;
+        std::cout << "    Produced volumes:      "
+                  << std::setw(width) << produced[0]/tot_porevol_init
+                  << std::setw(width) << produced[1]/tot_porevol_init
+                  << std::setw(width) << polyprod << std::endl;
+        std::cout << "    Total inj volumes:     "
+                  << std::setw(width) << tot_injected[0]/tot_porevol_init
+                  << std::setw(width) << tot_injected[1]/tot_porevol_init
+                  << std::setw(width) << tot_polyinj << std::endl;
+        std::cout << "    Total prod volumes:    "
+                  << std::setw(width) << tot_produced[0]/tot_porevol_init
+                  << std::setw(width) << tot_produced[1]/tot_porevol_init
+                  << std::setw(width) << tot_polyprod << std::endl;
+        std::cout << "    In-place + prod - inj: "
+                  << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
+                  << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init
+                  << std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
+        std::cout << "    Init - now - pr + inj: "
+                  << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
+                  << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
+                  << std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
+                  << std::endl;
+        std::cout.precision(8);
+
+        watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
+                      produced[0]/(produced[0] + produced[1]),
+                      tot_produced[0]/tot_porevol_init);
+        if (wells_) {
+            wellreport.push(props_, *wells_, state.saturation(),
+                            timer.simulationTimeElapsed() + timer.currentStepLength(),
+                            well_state.bhp(), well_state.perfRates());
+        }
 
         if (output_) {
             if (output_vtk_) {

opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp

@@ -0,0 +1,372 @@
+/*
+  Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+  
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+    namespace parameter { class ParameterGroup; }
+    class DerivedGeology;
+    class RockCompressibility;
+    class NewtonIterationBlackoilInterface;
+    class PolymerBlackoilState;
+    class WellStateFullyImplicitBlackoil;
+
+
+    /// A fully implicit solver for the black-oil-polymer problem.
+    ///
+    /// The simulator is capable of handling three-phase problems
+    /// where gas can be dissolved in oil (but not vice versa). It
+    /// uses an industry-standard TPFA discretization with per-phase
+    /// upwind weighting of mobilities.
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    template<class T>
+    class FullyImplicitBlackoilPolymerSolver
+    {
+    public:
+        /// \brief The type of the grid that we use.
+        typedef T Grid;
+        /// Construct a solver. It will retain references to the
+        /// arguments of this functions, and they are expected to
+        /// remain in scope for the lifetime of the solver.
+        /// \param[in] param       parameters
+        /// \param[in] grid             grid data structure
+        /// \param[in] fluid            fluid properties
+        /// \param[in] geo              rock properties
+        /// \param[in] rock_comp_props  if non-null, rock compressibility properties
+        /// \param[in] wells            well structure
+        /// \param[in] linsolver        linear solver
+        FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param,
+                                           const Grid&                     grid ,
+                                           const BlackoilPropsAdInterface& fluid,
+                                           const DerivedGeology&           geo  ,
+                                           const RockCompressibility*      rock_comp_props,
+                                           const PolymerPropsAd&           polymer_props_ad,
+                                           const Wells&                    wells,
+                                           const NewtonIterationBlackoilInterface& linsolver,
+                                           const bool has_disgas,
+                                           const bool has_vapoil,
+                                           const bool has_polymer);
+
+        /// \brief Set threshold pressures that prevent or reduce flow.
+        /// This prevents flow across faces if the potential
+        /// difference is less than the threshold. If the potential
+        /// difference is greater, the threshold value is subtracted
+        /// before calculating flow. This is treated symmetrically, so
+        /// flow is prevented or reduced in both directions equally.
+        /// \param[in]  threshold_pressures_by_face   array of size equal to the number of faces
+        ///                                   of the grid passed in the constructor.
+        void setThresholdPressures(const std::vector<double>& threshold_pressures_by_face);
+
+        /// Take a single forward step, modifiying
+        ///   state.pressure()
+        ///   state.faceflux()
+        ///   state.saturation()
+        ///   state.gasoilratio()
+        ///   wstate.bhp()
+        /// \param[in] dt        time step size
+        /// \param[in] state     reservoir state
+        /// \param[in] wstate    well state
+        void
+        step(const double   dt    ,
+             PolymerBlackoilState& state ,
+             WellStateFullyImplicitBlackoil& wstate,
+             const std::vector<double>& polymer_inflow);
+
+    private:
+        // Types and enums
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+
+        struct ReservoirResidualQuant {
+            ReservoirResidualQuant();
+            std::vector<ADB> accum; // Accumulations
+            ADB              mflux; // Mass flux (surface conditions)
+            ADB              b;     // Reciprocal FVF
+            ADB              head;  // Pressure drop across int. interfaces
+            ADB              mob;   // Phase mobility (per cell)
+        };
+
+        struct SolutionState {
+            SolutionState(const int np);
+            ADB              pressure;
+            std::vector<ADB> saturation;
+            ADB              rs;
+            ADB              rv;
+            ADB              concentration;
+            ADB              qs;
+            ADB              bhp;       
+        };
+
+        struct WellOps {
+            WellOps(const Wells& wells);
+            M w2p;              // well -> perf (scatter)
+            M p2w;              // perf -> well (gather)
+        };
+
+        enum { Water = BlackoilPropsAdInterface::Water,
+               Oil   = BlackoilPropsAdInterface::Oil  ,
+               Gas   = BlackoilPropsAdInterface::Gas  };
+
+        // the Newton relaxation type
+        enum RelaxType { DAMPEN, SOR };
+        enum PrimalVariables { Sg = 0, RS = 1, RV = 2 };
+
+        // Member data
+        const Grid&         grid_;
+        const BlackoilPropsAdInterface& fluid_;
+        const DerivedGeology&           geo_;
+        const RockCompressibility*      rock_comp_props_;
+        const PolymerPropsAd&           polymer_props_ad_; 
+        const Wells&                    wells_;
+        const NewtonIterationBlackoilInterface&    linsolver_;
+        // For each canonical phase -> true if active
+        const std::vector<bool>         active_;
+        // Size = # active phases. Maps active -> canonical phase indices.
+        const std::vector<int>          canph_;
+        const std::vector<int>          cells_;  // All grid cells
+        HelperOps                       ops_;
+        const WellOps                   wops_;
+        V                               cmax_;
+        const bool has_disgas_;
+        const bool has_vapoil_;
+        const bool has_polymer_;
+        const int  poly_pos_;
+        double                          dp_max_rel_;
+        double                          ds_max_;
+        double                          drs_max_rel_;
+        enum RelaxType                  relax_type_;
+        double                          relax_max_;
+        double                          relax_increment_;
+        double                          relax_rel_tol_;
+        int                             max_iter_;
+        bool use_threshold_pressure_;
+        V threshold_pressures_by_interior_face_;
+
+        std::vector<ReservoirResidualQuant> rq_;
+        std::vector<PhasePresence> phaseCondition_;
+        V well_perforation_pressure_diffs_; // Diff to bhp for each well perforation.
+
+        LinearisedBlackoilResidual residual_;
+
+        std::vector<int>         primalVariable_;
+
+        // Private methods.
+        SolutionState
+        constantState(const PolymerBlackoilState& x,
+                      const WellStateFullyImplicitBlackoil& xw);
+
+        SolutionState
+        variableState(const PolymerBlackoilState& x,
+                      const WellStateFullyImplicitBlackoil& xw);
+
+        void
+        computeAccum(const SolutionState& state,
+                     const int            aix  );
+
+        void computeWellConnectionPressures(const SolutionState& state,
+                                            const WellStateFullyImplicitBlackoil& xw);
+
+        void
+        addWellControlEq(const SolutionState& state,
+                         const WellStateFullyImplicitBlackoil& xw,
+                         const V& aliveWells);
+
+        void
+        addWellEq(const SolutionState& state,
+                  WellStateFullyImplicitBlackoil& xw,
+                  V& aliveWells,
+                  const std::vector<double>& polymer_inflow);
+
+        void updateWellControls(ADB& bhp,
+                                ADB& well_phase_flow_rate,
+                                WellStateFullyImplicitBlackoil& xw) const;
+
+        void
+        assemble(const V&             dtpv,
+                 const PolymerBlackoilState& x,
+                 WellStateFullyImplicitBlackoil& xw,
+                 const std::vector<double>& polymer_inflow);
+
+        V solveJacobianSystem() const;
+
+        void updateState(const V& dx,
+                         PolymerBlackoilState& state,
+                         WellStateFullyImplicitBlackoil& well_state);
+
+        std::vector<ADB>
+        computePressures(const SolutionState& state) const;
+
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+
+        std::vector<ADB>
+        computeRelPermWells(const SolutionState& state,
+                            const DataBlock& well_s,
+                            const std::vector<int>& well_cells) const;
+
+        void
+        computeMassFlux(const int               actph ,
+                        const V&                transi,
+                        const ADB&              kr    ,
+                        const ADB&              p     ,
+                        const SolutionState&    state );
+
+        void
+        computeMassFlux(const V&                trans,
+                        const std::vector<ADB>& kr,
+                        const std::vector<ADB>& phasePressure,
+                        const SolutionState&    state);
+
+        void
+        computeCmax(PolymerBlackoilState& state,
+                    const ADB&            c);
+
+        ADB
+        computeMc(const SolutionState& state) const;
+
+        ADB
+        rockPorosity(const ADB& p) const;
+        
+        ADB
+        rockPermeability(const ADB& p) const;
+        void applyThresholdPressures(ADB& dp);
+
+        double
+        residualNorm() const;
+
+        std::vector<double> residuals() const;
+
+        ADB
+        fluidViscosity(const int               phase,
+                       const ADB&              p    ,
+                       const ADB&              rs   ,
+                       const ADB&              rv   ,
+                       const std::vector<PhasePresence>& cond,
+                       const std::vector<int>& cells) const;
+
+        ADB
+        fluidReciprocFVF(const int               phase,
+                         const ADB&              p    ,
+                         const ADB&              rs   ,
+                         const ADB&              rv   ,
+                         const std::vector<PhasePresence>& cond,
+                         const std::vector<int>& cells) const;
+
+        ADB
+        fluidDensity(const int               phase,
+                     const ADB&              p    ,
+                     const ADB&              rs   ,
+                     const ADB&              rv   ,
+                     const std::vector<PhasePresence>& cond,
+                     const std::vector<int>& cells) const;
+
+        V
+        fluidRsSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRsSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+        V
+        fluidRvSat(const V&                p,
+                   const V&                so,
+                   const std::vector<int>& cells) const;
+
+        ADB
+        fluidRvSat(const ADB&              p,
+                   const ADB&              so,
+                   const std::vector<int>& cells) const;
+
+
+        ADB
+        poroMult(const ADB& p) const;
+
+        ADB
+        transMult(const ADB& p) const;
+
+        void
+        classifyCondition(const SolutionState&        state,
+                          std::vector<PhasePresence>& cond ) const;
+
+        const std::vector<PhasePresence>
+        phaseCondition() const {return phaseCondition_;}
+
+        void
+        classifyCondition(const PolymerBlackoilState&        state);
+
+
+        /// update the primal variable for Sg, Rv or Rs. The Gas phase must
+        /// be active to call this method.
+        void
+        updatePrimalVariableFromState(const PolymerBlackoilState&        state);
+
+        /// Update the phaseCondition_ member based on the primalVariable_ member.
+        void
+        updatePhaseCondFromPrimalVariable();
+
+        /// Compute convergence based on total mass balance (tol_mb) and maximum
+        /// residual mass balance (tol_cnv).
+        bool getConvergence(const double dt);
+
+        void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
+                                      const int it, const double relaxRelTol,
+                                      bool& oscillate, bool& stagnate) const;
+
+        void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const;
+
+        double dpMaxRel() const { return dp_max_rel_; }
+        double dsMax() const { return ds_max_; }
+        double drsMaxRel() const { return drs_max_rel_; }
+        enum RelaxType relaxType() const { return relax_type_; }
+        double relaxMax() const { return relax_max_; };
+        double relaxIncrement() const { return relax_increment_; };
+        double relaxRelTol() const { return relax_rel_tol_; };
+        double maxIter() const { return max_iter_; }
+
+    };
+} // namespace Opm
+
+#include "FullyImplicitBlackoilPolymerSolver_impl.hpp"
+
+#endif // OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -0,0 +1,246 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
+#include <opm/autodiff/LinearisedBlackoilResidual.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+    class DerivedGeology;
+    class RockCompressibility;
+    class NewtonIterationBlackoilInterface;
+    class PolymerBlackoilState;
+    class WellStateFullyImplicitBlackoil;
+
+    /// A fully implicit solver for the oil-water with polymer problem.
+    ///
+    /// The simulator is capable of handling oil-water-polymer problems
+    /// It uses an industry-standard TPFA discretization with per-phase
+    /// upwind weighting of mobilities.
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    class FullyImplicitCompressiblePolymerSolver
+    {
+    public:
+        /// Construct a solver. It will retain references to the
+        /// arguments of this functions, and they are expected to
+        /// remain in scope for the lifetime of the solver.
+        /// \param[in] grid             grid data structure
+        /// \param[in] fluid            fluid properties
+        /// \param[in] geo              rock properties
+        /// \param[in] rock_comp_props  if non-null, rock compressibility properties
+        /// \param[in] polymer_props_ad polymer properties
+        /// \param[in] wells            well structure
+        /// \param[in] linsolver        linear solver
+        FullyImplicitCompressiblePolymerSolver(const UnstructuredGrid&         grid ,
+        		                               const BlackoilPropsAdInterface& fluid,
+                   			                   const DerivedGeology&           geo  ,
+                              			       const RockCompressibility*      rock_comp_props,
+                                    		   const PolymerPropsAd&           polymer_props_ad,
+                                    		   const Wells&                    wells,
+                                    		   const NewtonIterationBlackoilInterface&    linsolver);
+
+        /// Take a single forward step, modifiying
+        ///   state.pressure()
+        ///   state.faceflux()
+        ///   state.saturation()
+        ///   state.concentration()
+        ///   wstate.bhp()
+        /// \param[in] dt        time step size
+        /// \param[in] state     reservoir state
+        /// \param[in] wstate    well state
+        /// \param[in] polymer_inflow	polymer influx
+        void
+        step(const double   			dt,
+             PolymerBlackoilState& 		state ,
+             WellStateFullyImplicitBlackoil& wstate,
+             const std::vector<double>& polymer_inflow);
+
+    private:
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+
+        struct ReservoirResidualQuant {
+            ReservoirResidualQuant();
+            std::vector<ADB> accum; // Accumulations
+            ADB              mflux; // Mass flux (surface conditions)
+            ADB              b;     // Reciprocal FVF
+            ADB              head;  // Pressure drop across int. interfaces
+            ADB              mob;   // Phase mobility (per cell)
+			std::vector<ADB> ads;   // Adsorption term.
+        };
+
+        struct SolutionState {
+            SolutionState(const int np);
+            ADB              pressure;
+            std::vector<ADB> saturation;
+            ADB              concentration;
+            ADB              qs;
+            ADB              bhp;
+        };
+
+        struct WellOps {
+            WellOps(const Wells& wells);
+            M w2p;              // well -> perf (scatter)
+            M p2w;              // perf -> well (gather)
+        };
+
+        enum { Water = BlackoilPropsAdInterface::Water,
+               Oil   = BlackoilPropsAdInterface::Oil  };
+
+        // Member data
+        const UnstructuredGrid&         grid_;
+        const BlackoilPropsAdInterface& fluid_;
+        const DerivedGeology&           geo_;
+        const RockCompressibility*      rock_comp_props_;
+        const PolymerPropsAd&           polymer_props_ad_;
+        const Wells&                    wells_;
+        const NewtonIterationBlackoilInterface&    linsolver_;
+        const std::vector<int>          cells_;  // All grid cells
+        HelperOps                       ops_;
+        const WellOps                   wops_;
+        const M                         grav_;
+		V    			 				cmax_;
+        std::vector<PhasePresence> phaseCondition_;
+        std::vector<ReservoirResidualQuant> rq_;
+        // The mass_balance vector has one element for each active phase,
+        // each of which has size equal to the number of cells.
+        // The well_eq has size equal to the number of wells.
+        LinearisedBlackoilResidual  residual_;
+        // Private methods.
+        SolutionState
+        constantState(const PolymerBlackoilState& x,
+                      const WellStateFullyImplicitBlackoil&     xw);
+
+        SolutionState
+        variableState(const PolymerBlackoilState& x,
+                      const WellStateFullyImplicitBlackoil&     xw);
+
+        void
+        computeAccum(const SolutionState& state,
+                     const int            aix  );
+
+        void
+        assemble(const double             	 dt,
+                 const PolymerBlackoilState& x,
+                 const WellStateFullyImplicitBlackoil& xw,  
+                 const std::vector<double>& polymer_inflow);
+
+        V solveJacobianSystem() const;
+
+        void updateState(const V& dx,
+                         PolymerBlackoilState& state,
+                         WellStateFullyImplicitBlackoil& well_state) const;
+
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+
+        std::vector<ADB>
+        computeRelPermWells(const SolutionState& state,
+                            const DataBlock& well_s,
+                            const std::vector<int>& well_cells) const;
+		std::vector<ADB>
+		computePressures(const SolutionState& state) const;
+
+        void
+        computeMassFlux(const int               actph ,
+                        const V&                transi,
+                        const std::vector<ADB>& kr    ,
+                        const SolutionState&    state );
+
+        void
+        computeMassFlux(const V&                trans,
+                        const ADB&              mc,
+                        const ADB&              kro,
+                        const ADB&              krw_eff,
+                        const SolutionState&    state);
+
+        std::vector<ADB>
+        computeFracFlow(const ADB&              kro,
+                        const ADB&              krw_eff,
+                        const ADB&              c) const;
+
+        void
+        computeCmax(PolymerBlackoilState& state,
+					const ADB& c);
+
+        ADB
+        computeMc(const SolutionState&  state) const;
+
+        ADB
+        rockPorosity(const ADB& p) const;
+
+        ADB
+        rockPermeability(const ADB& p) const;
+
+        double
+        residualNorm() const;
+
+        ADB
+        fluidViscosity(const int                         phase,
+                       const ADB&                        p    ,
+                       const std::vector<PhasePresence>& cond,
+                       const std::vector<int>&           cells) const;
+
+        ADB
+        fluidReciprocFVF(const int                         phase,
+                         const ADB&                        p    ,
+                         const std::vector<PhasePresence>& cond,
+                         const std::vector<int>&           cells) const;
+
+        ADB
+        fluidDensity(const int                         phase,
+                     const ADB&                        p    ,
+                     const std::vector<PhasePresence>& cond,
+                     const std::vector<int>&           cells) const;
+
+        ADB
+        poroMult(const ADB& p) const;
+
+        ADB
+        transMult(const ADB& p) const;
+        
+        const std::vector<PhasePresence>
+        phaseCondition() const { return phaseCondition_; }
+        
+        void
+        classifyCondition(const PolymerBlackoilState& state);
+    };
+} // namespace Opm
+
+
+#endif // OPM_FULLYIMPLICITBLACKOILSOLVER_HEADER_INCLUDED

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.cpp

@@ -0,0 +1,896 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
+
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/IncompPropsAdInterface.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/PolymerState.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/well_controls.h>
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <iomanip>
+#include <Eigen/Eigen>
+#include <algorithm>
+
+namespace Opm {
+
+
+
+
+
+typedef AutoDiffBlock<double> ADB;
+typedef ADB::V V;
+typedef ADB::M M;
+typedef Eigen::Array<double,
+                     Eigen::Dynamic,
+                     Eigen::Dynamic,
+                     Eigen::RowMajor> DataBlock;
+
+
+
+
+
+namespace {
+    
+    std::vector<int>
+    buildAllCells(const int nc)
+    {
+        std::vector<int> all_cells(nc);
+        for (int c = 0; c < nc; ++c) { all_cells[c] = c; }
+
+        return all_cells;
+    }
+    struct Chop01 {
+        double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); }
+    };
+
+
+
+
+
+    V computePerfPress(const UnstructuredGrid& grid, const Wells& wells, const V& rho, const double grav)
+    {
+        const int nw = wells.number_of_wells;
+        const int nperf = wells.well_connpos[nw];
+        const int dim = grid.dimensions;
+        V wdp = V::Zero(nperf,1);
+        assert(wdp.size() == rho.size());
+
+        // Main loop, iterate over all perforations,
+        // using the following formula:
+        //    wdp(perf) = g*(perf_z - well_ref_z)*rho(perf)
+        // where the total density rho(perf) is taken to be
+        //    sum_p (rho_p*saturation_p) in the perforation cell.
+        // [although this is computed on the outside of this function].
+        for (int w = 0; w < nw; ++w) {
+            const double ref_depth = wells.depth_ref[w];
+            for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) {
+                const int cell = wells.well_cells[j];
+                const double cell_depth = grid.cell_centroids[dim * cell + dim - 1];
+                wdp[j] = rho[j]*grav*(cell_depth - ref_depth);
+            }
+        }
+        return wdp;
+    }
+
+} //anonymous namespace
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::
+    FullyImplicitTwophasePolymerSolver(const UnstructuredGrid&         grid,
+                                       const IncompPropsAdInterface&   fluid,
+                                       const PolymerPropsAd&           polymer_props_ad,
+                                       const LinearSolverInterface&    linsolver,
+                                       const Wells&                    wells,
+                                       const double*                   gravity)
+        : grid_ (grid)
+        , fluid_(fluid)
+        , polymer_props_ad_ (polymer_props_ad)
+        , linsolver_(linsolver)
+        , wells_(wells)
+        , gravity_(gravity)
+        , cells_ (buildAllCells(grid.number_of_cells))
+        , ops_(grid)
+        , wops_(wells)
+//        , mob_(std::vector<ADB>(fluid.numPhases() + 1, ADB::null()))
+		, cmax_(V::Zero(grid.number_of_cells))
+        , rq_(fluid.numPhases() + 1)
+        , residual_( { std::vector<ADB>(fluid.numPhases() + 1, ADB::null()), ADB::null(), ADB::null()})
+     {
+     }
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::
+    WellOps::WellOps(const Wells& wells)
+        : w2p(wells.well_connpos[ wells.number_of_wells ],
+              wells.number_of_wells)
+        , p2w(wells.number_of_wells,
+              wells.well_connpos[ wells.number_of_wells ])
+    {
+        const int        nw   = wells.number_of_wells;
+        const int* const wpos = wells.well_connpos;
+
+        typedef Eigen::Triplet<double> Tri;
+
+        std::vector<Tri> scatter, gather;
+        scatter.reserve(wpos[nw]);
+        gather .reserve(wpos[nw]);
+
+        for (int w = 0, i = 0; w < nw; ++w) {
+            for (; i < wpos[ w + 1 ]; ++i) {
+                scatter.push_back(Tri(i, w, 1.0));
+                gather .push_back(Tri(w, i, 1.0));
+            }
+        }
+
+        w2p.setFromTriplets(scatter.begin(), scatter.end());
+        p2w.setFromTriplets(gather .begin(), gather .end());
+    }
+
+
+
+
+
+    void
+    FullyImplicitTwophasePolymerSolver::
+    step(const double   dt,
+         PolymerState&  x,
+         WellState&     xw,
+         const std::vector<double>& polymer_inflow,
+		 std::vector<double>& src)
+    {
+        
+        V pvol(grid_.number_of_cells);
+        // Pore volume
+        const V::Index nc = grid_.number_of_cells;
+        V rho = V::Constant(pvol.size(), 1, *fluid_.porosity());
+        std::transform(grid_.cell_volumes, grid_.cell_volumes + nc,
+                       rho.data(), pvol.data(),
+                       std::multiplies<double>());
+
+        const V pvdt = pvol / dt;
+
+        const SolutionState old_state = constantState(x, xw);
+		computeCmax(x, old_state.concentration);
+		computeAccum(old_state, 0);
+        const double atol  = 1.0e-12;
+        const double rtol  = 5.0e-8;
+        const int    maxit = 40;
+        assemble(pvdt, x, xw, polymer_inflow, src);
+
+        const double r0  = residualNorm();
+        int          it  = 0;
+        std::cout << "\nIteration         Residual\n"
+                  << std::setw(9) << it << std::setprecision(9)
+                  << std::setw(18) << r0 << std::endl;
+        bool resTooLarge = r0 > atol;
+        while (resTooLarge && (it < maxit)) {
+            const V dx = solveJacobianSystem();
+            updateState(dx, x, xw);
+
+            assemble(pvdt, x, xw, polymer_inflow, src);
+
+            const double r = residualNorm();
+
+            resTooLarge = (r > atol) && (r > rtol*r0);
+
+            it += 1;
+            std::cout << std::setw(9) << it << std::setprecision(9)
+                      << std::setw(18) << r << std::endl;
+        }
+
+        if (resTooLarge) {
+            std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n";
+            // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations.");
+        }
+    }
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::ReservoirResidualQuant::ReservoirResidualQuant()
+        : accum(2, ADB::null())
+        , mflux(   ADB::null())
+        , b    (   ADB::null())
+        , head (   ADB::null())
+        , mob  (   ADB::null())
+    {
+    }
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::SolutionState::SolutionState(const int np)
+        : pressure   (    ADB::null())
+        , saturation (np, ADB::null())
+        , concentration ( ADB::null())
+        , qs            ( ADB::null())
+        , bhp           ( ADB::null())
+    {
+    }
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::SolutionState
+    FullyImplicitTwophasePolymerSolver::constantState(const PolymerState& x,
+                                                      const WellState&    xw)
+    {
+        const int nc = grid_.number_of_cells;
+        const int np = x.numPhases();
+        // The block pattern assumes the following primary variables:
+        //    pressure
+        //    water saturation 
+        //    polymer concentration
+        //    well surface rates
+        //    well bottom-hole pressure
+        // Note that oil is assumed to always be present, but is never
+        // a primary variable.
+        std::vector<int> bpat(np + 1, nc);
+        bpat.push_back(xw.bhp().size() * np);
+        bpat.push_back(xw.bhp().size());
+        
+        SolutionState state(np);
+
+        // Pressure.
+        assert (not x.pressure().empty());
+        const V p = Eigen::Map<const V>(& x.pressure()[0], nc);
+        state.pressure = ADB::constant(p);
+
+        // Saturation.
+        assert (not x.saturation().empty());
+        const DataBlock s_all = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
+        for (int phase = 0; phase < np; ++phase) {
+            state.saturation[phase] = ADB::constant(s_all.col(phase));
+        }
+
+        // Concentration
+        assert(not x.concentration().empty());
+        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
+        state.concentration = ADB::constant(c);
+
+        // Well rates.
+        assert (not xw.wellRates().empty());
+        // Need to reshuffle well rates, from ordered by wells, then phase,
+        // to ordered by phase, then wells.
+        const int nw = wells_.number_of_wells;
+        // The transpose() below switches the ordering.
+        const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
+        const V qs = Eigen::Map<const V>(wrates.data(), nw * np);
+        state.qs = ADB::constant(qs, bpat);
+
+        // Bottom hole pressure.
+        assert (not xw.bhp().empty());
+        const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
+        state.bhp = ADB::constant(bhp, bpat);
+        return state;
+    }
+
+
+
+
+
+    FullyImplicitTwophasePolymerSolver::SolutionState
+    FullyImplicitTwophasePolymerSolver::variableState(const PolymerState& x,
+                                                      const WellState&    xw)
+    {
+        const int nc = grid_.number_of_cells;
+        const int np = x.numPhases();
+
+        std::vector<V> vars0;
+        vars0.reserve(np + 3); 
+
+        // Initial pressure.
+        assert (not x.pressure().empty());
+        const V p = Eigen::Map<const V>(& x.pressure()[0], nc);
+        vars0.push_back(p);
+
+        // Initial saturation.
+        assert (not x.saturation().empty());
+        const DataBlock s_all = Eigen::Map<const DataBlock>(& x.saturation()[0], nc, np);
+        const V sw = s_all.col(0);
+        vars0.push_back(sw);
+
+        // Initial concentration.
+        assert (not x.concentration().empty());
+        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
+        vars0.push_back(c);
+
+        // Initial well rates.
+        assert (not xw.wellRates().empty());
+        // Need to reshuffle well rates, from ordered by wells, then phase,
+        // to ordered by phase, then wells.
+        const int nw = wells_.number_of_wells;
+        // The transpose() below switches the ordering.
+        const DataBlock wrates = Eigen::Map<const DataBlock>(& xw.wellRates()[0], nw, np).transpose();
+        const V qs = Eigen::Map<const V>(wrates.data(), nw * np);
+        vars0.push_back(qs);
+
+        // Initial well bottom hole pressure.
+        assert (not xw.bhp().empty());
+        const V bhp = Eigen::Map<const V>(& xw.bhp()[0], xw.bhp().size());
+        vars0.push_back(bhp);
+
+        std::vector<ADB> vars = ADB::variables(vars0);
+
+        SolutionState state(np);
+
+        // Pressure.
+        int nextvar = 0;
+        state.pressure = vars[ nextvar++ ];
+
+        // Saturation.
+        const std::vector<int>& bpat = vars[0].blockPattern();
+        {
+            ADB so  = ADB::constant(V::Ones(nc, 1), bpat);
+            ADB sw = vars[ nextvar++ ];
+            state.saturation[0] = sw;
+            so = so - sw;
+            state.saturation[1] = so;
+        }
+        
+        // Concentration.
+        state.concentration = vars[nextvar++];
+
+        // Qs.
+        state.qs = vars[ nextvar++ ];
+
+        // BHP.
+        state.bhp = vars[ nextvar++ ];
+        assert(nextvar == int(vars.size()));
+
+        return state;
+    }
+  
+
+
+
+
+    void
+    FullyImplicitTwophasePolymerSolver::
+    computeCmax(PolymerState& state,
+				const ADB& c)
+    {
+        const int nc = grid_.number_of_cells;
+        for (int i = 0; i < nc; ++i) {
+		    cmax_(i) = std::max(cmax_(i), c.value()(i));
+        }
+		std::copy(&cmax_[0], &cmax_[0] + nc, state.maxconcentration().begin());
+
+    }
+
+
+
+
+
+    void
+    FullyImplicitTwophasePolymerSolver::
+	computeAccum(const SolutionState& state,
+                 const int            aix  )
+    {
+
+        const std::vector<ADB>& sat   = state.saturation;
+        const ADB&              c     = state.concentration;
+
+        rq_[0].accum[aix] = sat[0];
+        rq_[1].accum[aix] = sat[1];
+
+		const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
+        const ADB ads = polymer_props_ad_.adsorption(state.concentration, cmax);
+        const double rho_rock = polymer_props_ad_.rockDensity();
+        const V phi = Eigen::Map<const V>(&fluid_.porosity()[0], grid_.number_of_cells, 1);
+        const double dead_pore_vol = polymer_props_ad_.deadPoreVol();
+
+        rq_[2].accum[aix] = sat[0] * c * (1. - dead_pore_vol) + rho_rock * (1. - phi) / phi * ads;
+    }
+
+
+
+
+
+    void
+    FullyImplicitTwophasePolymerSolver::
+    assemble(const V&             pvdt,
+             const PolymerState&  x,
+             const WellState&     xw,
+             const std::vector<double>& polymer_inflow,
+			 std::vector<double>& src)
+    {
+        // Create the primary variables.
+        const SolutionState state = variableState(x, xw);
+		computeAccum(state, 1);
+        // -------- Mass balance equations for water and oil --------
+        const V trans = subset(transmissibility(), ops_.internal_faces);
+        const std::vector<ADB> kr = computeRelPerm(state);
+		const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern());
+        const ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration, cmax, kr[0], state.saturation[0]);
+        const ADB mc = computeMc(state);
+        computeMassFlux(trans, mc, kr[1], krw_eff, state);
+        residual_.mass_balance[0] = pvdt*(rq_[0].accum[1] - rq_[0].accum[0])
+                                    + ops_.div*rq_[0].mflux;
+        residual_.mass_balance[1] = pvdt*(rq_[1].accum[1] - rq_[1].accum[0])
+                                    + ops_.div*rq_[1].mflux;
+        residual_.mass_balance[2] = pvdt*(rq_[2].accum[1] - rq_[2].accum[0])
+                                    + ops_.div*rq_[2].mflux;
+
+        // -------- Well equation, and well contributions to the mass balance equations --------
+
+        // Contribution to mass balance will have to wait.
+
+        const int nc = grid_.number_of_cells;
+        const int np = wells_.number_of_phases;
+        const int nw = wells_.number_of_wells;
+        const int nperf = wells_.well_connpos[nw];
+
+        const std::vector<int> well_cells(wells_.well_cells, wells_.well_cells + nperf);
+        const V transw = Eigen::Map<const V>(wells_.WI, nperf);
+
+        const ADB& bhp = state.bhp;
+
+        const DataBlock well_s = wops_.w2p * Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np).matrix();
+
+        // Extract variables for perforation cell pressures
+        // and corresponding perforation well pressures.
+        const ADB p_perfcell = subset(state.pressure, well_cells);
+        // Finally construct well perforation pressures and well flows.
+
+        // Compute well pressure differentials.
+        // Construct pressure difference vector for wells.
+        const int dim = grid_.dimensions;
+        if (gravity_) {
+            for (int dd = 0; dd < dim -1; ++dd) {
+                assert(g[dd] == 0.0);
+            }
+        }
+        ADB cell_rho_total = ADB::constant(V::Zero(nc), state.pressure.blockPattern());
+        for (int phase = 0; phase < 2; ++phase) {
+            // For incompressible flow cell rho is the same.
+            const ADB cell_rho = fluidDensity(phase, state.pressure);
+            cell_rho_total += state.saturation[phase] * cell_rho;
+        }
+        ADB inj_rho_total = ADB::constant(V::Zero(nperf), state.pressure.blockPattern());
+        assert(np == wells_.number_of_phases);
+        const DataBlock compi = Eigen::Map<const DataBlock>(wells_.comp_frac, nw, np);
+        for (int phase = 0; phase < 2; ++phase) {
+            const ADB cell_rho = fluidDensity(phase, state.pressure);
+            const V fraction = compi.col(phase);
+            inj_rho_total += (wops_.w2p * fraction.matrix()).array() * subset(cell_rho, well_cells);
+        }
+        const V rho_perf_cell = subset(cell_rho_total, well_cells).value();
+        const V rho_perf_well = inj_rho_total.value();
+        V prodperfs = V::Constant(nperf, -1.0);
+        for (int w = 0; w < nw; ++w) {
+            if (wells_.type[w] == PRODUCER) {
+                std::fill(prodperfs.data() + wells_.well_connpos[w],
+                          prodperfs.data() + wells_.well_connpos[w+1], 1.0);
+            }
+        }
+        const Selector<double> producer(prodperfs);
+        const V rho_perf = producer.select(rho_perf_cell, rho_perf_well);
+        const V well_perf_dp = computePerfPress(grid_, wells_, rho_perf, gravity_ ? gravity_[dim - 1] : 0.0);
+
+        const ADB p_perfwell = wops_.w2p * bhp + well_perf_dp;
+        const ADB nkgradp_well = transw * (p_perfcell - p_perfwell);
+        // DUMP(nkgradp_well);
+        const Selector<double> cell_to_well_selector(nkgradp_well.value());
+        ADB well_rates_all = ADB::constant(V::Zero(nw*np), state.bhp.blockPattern());
+
+        ADB perf_total_mob = subset(rq_[0].mob, well_cells) + subset(rq_[1].mob, well_cells);
+
+        std::vector<ADB> well_contribs(np, ADB::null());
+        std::vector<ADB> well_perf_rates(np, ADB::null());
+        for (int phase = 0; phase < np; ++phase) {
+//            const ADB& cell_b = rq_[phase].b;
+  //          const ADB perf_b = subset(cell_b, well_cells);
+            const ADB& cell_mob = rq_[phase].mob;
+            const V well_fraction = compi.col(phase);
+            // Using total mobilities for all phases for injection.
+            const ADB perf_mob_injector = (wops_.w2p * well_fraction.matrix()).array() * perf_total_mob;
+            const ADB perf_mob = producer.select(subset(cell_mob, well_cells),
+                                                 perf_mob_injector);
+            const ADB perf_flux = perf_mob * (nkgradp_well); // No gravity term for perforations.
+            well_perf_rates[phase] = perf_flux;
+            const ADB well_rates = wops_.p2w * well_perf_rates[phase];
+            well_rates_all += superset(well_rates, Span(nw, 1, phase*nw), nw*np);
+
+            // const ADB well_contrib = superset(perf_flux*perf_b, well_cells, nc);
+            well_contribs[phase] = superset(perf_flux, well_cells, nc);
+            // DUMP(well_contribs[phase]);
+            residual_.mass_balance[phase] += well_contribs[phase];
+			for (int i = 0; i < nc; ++i) {
+				src[i] += well_contribs[phase].value()[i];
+			}
+        }
+
+        // well rates contribs to polymer mass balance eqn.
+        // for injection wells.
+        const V polyin = Eigen::Map<const V>(& polymer_inflow[0], nc);
+        const V poly_in_perf = subset(polyin, well_cells);
+        const V poly_mc_cell = subset(mc, well_cells).value();
+        const V poly_c = producer.select(poly_mc_cell, poly_in_perf);
+        residual_.mass_balance[2] += superset(well_perf_rates[0] * poly_c, well_cells, nc);
+
+        // Set the well flux equation
+        residual_.well_flux_eq = state.qs + well_rates_all;
+        // DUMP(residual_.well_flux_eq);
+
+        // Handling BHP and SURFACE_RATE wells.
+        V bhp_targets(nw);
+        V rate_targets(nw);
+        M rate_distr(nw, np*nw);
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells_.ctrls[w];
+			if (well_controls_get_current_type(wc) == BHP) {
+				bhp_targets[w] = well_controls_get_current_target(wc);
+                rate_targets[w] = -1e100;
+            } else if (well_controls_get_current_type(wc) == SURFACE_RATE) {
+                bhp_targets[w] = -1e100;
+				rate_targets[w] = well_controls_get_current_target(wc);
+				{	
+					const double* distr = well_controls_get_current_distr(wc);
+            	    for (int phase = 0; phase < np; ++phase) {
+                	    rate_distr.insert(w, phase*nw + w) = distr[phase];
+                	}
+				}
+            } else {
+                OPM_THROW(std::runtime_error, "Can only handle BHP type controls.");
+            }
+        }
+        const ADB bhp_residual = bhp - bhp_targets;
+        const ADB rate_residual = rate_distr * state.qs - rate_targets;
+        // Choose bhp residual for positive bhp targets.
+        Selector<double> bhp_selector(bhp_targets);
+        residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual);
+ //       residual_.well_eq = bhp_residual;
+
+    }
+
+
+
+
+   
+    std::vector<ADB>
+    FullyImplicitTwophasePolymerSolver::
+	computePressures(const SolutionState& state) const
+    {
+        const ADB sw = state.saturation[0];
+        const ADB so = state.saturation[1];
+
+        // convert the pressure offsets to the capillary pressures
+        std::vector<ADB> pressure = fluid_.capPress(sw, so, cells_);
+		pressure[0] = pressure[0] - pressure[1];
+
+        // add the total pressure to the capillary pressures
+        for (int phaseIdx = 0; phaseIdx < 2; ++phaseIdx) {
+            pressure[phaseIdx] += state.pressure;
+        }
+
+        return pressure;
+    }
+
+
+
+
+
+    void
+    FullyImplicitTwophasePolymerSolver::computeMassFlux(const V&                trans,
+                                                        const ADB&              mc,
+                                                        const ADB&              kro,
+                                                        const ADB&              krw_eff,
+                                                        const SolutionState&    state )
+    {
+        const double* mus = fluid_.viscosity();
+        ADB inv_wat_eff_vis = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mus);
+        rq_[0].mob = krw_eff * inv_wat_eff_vis;
+        rq_[1].mob = kro / V::Constant(kro.size(), 1, mus[1]);
+        rq_[2].mob =  mc * krw_eff * inv_wat_eff_vis;
+
+        const int nc = grid_.number_of_cells; 
+        V z(nc);
+        // Compute z coordinates
+        for (int c = 0; c < nc; ++c){
+            z[c] = grid_.cell_centroids[c * 3 + 2];
+        }
+		std::vector<ADB> press = computePressures(state);
+        for (int phase = 0; phase < 2; ++phase) {
+            const ADB rho = fluidDensity(phase, state.pressure);
+			ADB& head = rq_[phase].head;
+            const ADB rhoavg = ops_.caver * rho;
+            const ADB dp = ops_.ngrad * press[phase]
+                           - gravity_[2] * (rhoavg * (ops_.ngrad * z.matrix()));
+            head = trans * dp;
+            UpwindSelector<double> upwind(grid_, ops_, head.value());
+			const ADB& mob = rq_[phase].mob;
+            rq_[phase].mflux = upwind.select(mob) * head;
+        }
+		rq_[2].head = rq_[0].head;
+		UpwindSelector<double> upwind(grid_, ops_, rq_[2].head.value());
+		rq_[2].mflux = upwind.select(rq_[2].mob) * rq_[2].head;
+    }
+
+
+
+
+
+    std::vector<ADB>
+    FullyImplicitTwophasePolymerSolver::accumSource(const ADB&                 kro,
+                                                    const ADB&                 krw_eff,
+                                                    const ADB&                 c,
+                                                    const std::vector<double>& src,
+                                                    const std::vector<double>& polymer_inflow_c) const
+    {
+        //extract the source to out and in source.
+        std::vector<double> outsrc;
+        std::vector<double> insrc;
+        std::vector<double>::const_iterator it;
+        for (it = src.begin(); it != src.end(); ++it) {
+            if (*it < 0) {
+                outsrc.push_back(*it);
+                insrc.push_back(0.0);
+            } else if (*it > 0) {
+                insrc.push_back(*it);
+                outsrc.push_back(0.0);
+            } else {
+                outsrc.push_back(0);
+                insrc.push_back(0);
+            }
+        }
+        const V outSrc = Eigen::Map<const V>(& outsrc[0], grid_.number_of_cells);
+        const V inSrc = Eigen::Map<const V>(& insrc[0], grid_.number_of_cells);
+        const V polyin = Eigen::Map<const V>(& polymer_inflow_c[0], grid_.number_of_cells);
+        // compute the out-fracflow.
+        const std::vector<ADB> f = computeFracFlow();
+        // compute the in-fracflow.
+        V   zero = V::Zero(grid_.number_of_cells);
+        V   one  = V::Ones(grid_.number_of_cells);
+
+        std::vector<ADB> source;
+        //water source
+        source.push_back(f[0] * outSrc + one * inSrc);
+        //oil source
+        source.push_back(f[1] * outSrc + zero * inSrc);
+        //polymer source
+        source.push_back(f[0] * outSrc * c + one * inSrc * polyin);
+
+        return source;
+     }
+
+
+
+
+
+    std::vector<ADB>
+    FullyImplicitTwophasePolymerSolver::computeFracFlow() const
+    {
+        ADB total_mob = rq_[0].mob + rq_[1].mob;
+
+        std::vector<ADB> fracflow;
+
+        fracflow.push_back(rq_[0].mob / total_mob);
+        fracflow.push_back(rq_[1].mob / total_mob);
+
+        return fracflow;
+    }
+
+
+
+
+
+    V 
+    FullyImplicitTwophasePolymerSolver::solveJacobianSystem() const
+    {
+        const int np = fluid_.numPhases();
+    	if (np != 2) {
+	        OPM_THROW(std::logic_error, "Only two-phase ok in FullyImplicitTwophasePolymerSolver.");
+	    }
+	    ADB mass_res = vertcat(residual_.mass_balance[0], residual_.mass_balance[1]);
+        mass_res = vertcat(mass_res, residual_.mass_balance[2]);
+        ADB well_res = vertcat(residual_.well_flux_eq, residual_.well_eq);
+	    ADB total_res = collapseJacs(vertcat(mass_res, well_res));
+
+        const Eigen::SparseMatrix<double, Eigen::RowMajor> matr = total_res.derivative()[0];
+        V dx(V::Zero(total_res.size()));
+        Opm::LinearSolverInterface::LinearSolverReport rep
+            = linsolver_.solve(matr.rows(), matr.nonZeros(),
+                               matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
+                               total_res.value().data(), dx.data());
+        if (!rep.converged) {
+            OPM_THROW(std::runtime_error,
+                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
+                      "Linear solver convergence failure.");
+        }
+        return dx;
+    }
+
+
+
+
+
+    void FullyImplicitTwophasePolymerSolver::updateState(const V& dx,
+                                                  PolymerState& state,
+                                                  WellState&    well_state) const
+    {
+        const int np = fluid_.numPhases();
+        const int nc = grid_.number_of_cells;
+        const int nw = wells_.number_of_wells;
+        const V one = V::Constant(nc, 1.0);
+        const V zero = V::Zero(nc);
+
+        // Extract parts of dx corresponding to each part.
+        const V dp = subset(dx, Span(nc));
+        int varstart = nc;
+        const V dsw = subset(dx, Span(nc, 1, varstart));
+        varstart += dsw.size();
+        const V dc = subset(dx, Span(nc, 1, varstart));
+        varstart += dc.size();
+        const V dqs = subset(dx, Span(np*nw, 1, varstart));
+        varstart += dqs.size();
+        const V dbhp = subset(dx, Span(nw, 1, varstart));
+        varstart += dbhp.size();
+
+        assert(varstart == dx.size());
+
+        // Pressure update.
+        const V p_old = Eigen::Map<const V>(&state.pressure()[0], nc);
+        const V p = p_old - dp;
+        std::copy(&p[0], &p[0] + nc, state.pressure().begin());
+
+        // Saturation updates.
+        const double dsmax = 0.3;
+        const DataBlock s_old = Eigen::Map<const DataBlock>(& state.saturation()[0], nc, np);
+        V so = one;
+        const V sw_old = s_old.col(0);
+        const V dsw_limited = sign(dsw) * dsw.abs().min(dsmax);
+        const V sw = (sw_old - dsw_limited).unaryExpr(Chop01());
+        so -= sw;
+        for (int c = 0; c < nc; ++c) {
+            state.saturation()[c*np] = sw[c];
+            state.saturation()[c*np + 1] = so[c];
+        }
+        
+        // Concentration updates.
+        const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc);
+        const V c = (c_old - dc).max(zero);
+        std::copy(&c[0], &c[0] + nc, state.concentration().begin());
+
+        // Qs update.
+        // Since we need to update the wellrates, that are ordered by wells,
+        // from dqs which are ordered by phase, the simplest is to compute
+        // dwr, which is the data from dqs but ordered by wells.
+        const DataBlock wwr = Eigen::Map<const DataBlock>(dqs.data(), np, nw).transpose();
+        const V dwr = Eigen::Map<const V>(wwr.data(), nw*np);
+        const V wr_old = Eigen::Map<const V>(&well_state.wellRates()[0], nw*np);
+        const V wr = wr_old - dwr;
+        std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin());
+
+        // Bhp update.
+        const V bhp_old = Eigen::Map<const V>(&well_state.bhp()[0], nw, 1);
+        const V bhp = bhp_old - dbhp;
+        std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin());
+    }
+   
+
+
+
+
+    std::vector<ADB>
+    FullyImplicitTwophasePolymerSolver::computeRelPerm(const SolutionState& state) const
+    {
+
+        const ADB sw = state.saturation[0];
+        const ADB so = state.saturation[1];
+
+        return fluid_.relperm(sw, so, cells_);
+    }
+    
+    
+    
+    
+    
+    double
+    FullyImplicitTwophasePolymerSolver::residualNorm() const
+    {
+        double r = 0;
+        for (std::vector<ADB>::const_iterator
+                 b = residual_.mass_balance.begin(),
+                 e = residual_.mass_balance.end();
+             b != e; ++b)
+        {
+            r = std::max(r, (*b).value().matrix().lpNorm<Eigen::Infinity>());
+        }
+
+        r = std::max(r, residual_.well_flux_eq.value().matrix().lpNorm<Eigen::Infinity>());
+        r = std::max(r, residual_.well_eq.value().matrix().lpNorm<Eigen::Infinity>());
+
+        return r;
+    }
+   
+   
+   
+
+
+    ADB
+    FullyImplicitTwophasePolymerSolver::fluidDensity(const int phase,
+		                                             const ADB p) const
+    {
+        const double* rhos = fluid_.surfaceDensity();
+        ADB rho = ADB::constant(V::Constant(grid_.number_of_cells, 1, rhos[phase]),
+                               p.blockPattern());
+        
+        return rho;
+    }
+  
+
+
+ 
+   
+    V
+    FullyImplicitTwophasePolymerSolver::transmissibility() const
+    {
+        const V::Index nc = grid_.number_of_cells;
+        V htrans(grid_.cell_facepos[nc]);
+        V trans(grid_.cell_facepos[nc]);
+        UnstructuredGrid* ug = const_cast<UnstructuredGrid*>(& grid_);
+        tpfa_htrans_compute(ug, fluid_.permeability(), htrans.data());
+        tpfa_trans_compute (ug, htrans.data(), trans.data());
+        
+        return trans;
+    }
+
+
+
+
+
+    // here mc means m(c) * c. 
+    ADB
+    FullyImplicitTwophasePolymerSolver::computeMc(const SolutionState& state) const
+    {
+        ADB c = state.concentration;
+        return polymer_props_ad_.polymerWaterVelocityRatio(c);
+    }
+
+
+
+
+
+} //namespace Opm

opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp

@@ -0,0 +1,217 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED
+#define OPM_FULLYIMPLICITTWOPHASEPOLYMERSOLVER_HEADER_INCLUDED
+
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/autodiff/IncompPropsAdInterface.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/core/pressure/tpfa/trans_tpfa.h>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm {
+
+    class LinearSolverInterface;
+    class PolymerState;
+    class WellState;
+    
+    /// A fully implicit solver for the incompressible oil-water flow wtih polymer problem.
+    ///
+    /// The simulator is capable of handling oil-water with polymer problems
+    /// It uses an industry-standard TPFA discretization with per-phase
+    /// upwind weighting of mobilities.
+    ///
+    /// It uses automatic differentiation via the class AutoDiffBlock
+    /// to simplify assembly of the jacobian matrix.
+    class FullyImplicitTwophasePolymerSolver
+    {
+    public:
+        /// Construct a solver. It will retain references to the
+        /// arguments of this functions, and they are expected to
+        /// remain in scope for the lifetime of the solver.
+        /// \param[in] grid             grid data structure
+        /// \param[in] fluid            fluid properties
+        /// \param[in] polymer_props_ad polymer properties
+        /// \param[in] wells            well structure
+        /// \param[in] linsolver        linear solver
+        /// \param[in] gravity          gravity
+        FullyImplicitTwophasePolymerSolver(const UnstructuredGrid&        grid,
+                                           const IncompPropsAdInterface&  fluid,
+                                           const PolymerPropsAd&          polymer_props_ad,
+                                           const LinearSolverInterface&   linsolver,
+                                           const Wells&                   wells,
+                                           const double*                  gravity);
+
+        /// Take a single forward step, modifiying
+        ///   state.pressure()
+        ///   state.faceflux()
+        ///   state.saturation()
+        ///   state.concentration()
+        ///   wstate.bhp()
+        /// \param[in] dt        		time step size
+        /// \param[in] state     		reservoir state with polymer
+        /// \param[in] wstate    		well state
+        /// \param[in] polymer_inflow	polymer influx
+        /// \param[in] src				to caculate wc
+        void step(const double   dt,
+                  PolymerState& state,
+                  WellState&    well_state,
+                  const std::vector<double>& polymer_inflow,
+				  std::vector<double>& src);
+    private:
+        typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+        typedef ADB::M M;
+        typedef Eigen::Array<double,
+                             Eigen::Dynamic,
+                             Eigen::Dynamic,
+                             Eigen::RowMajor> DataBlock;
+
+        struct ReservoirResidualQuant {
+            ReservoirResidualQuant();
+            std::vector<ADB> accum; // Accumulations
+            ADB              mflux; // Mass flux (surface conditions)
+            ADB              b;     // Reciprocal FVF
+            ADB              head;  // Pressure drop across int. interfaces
+            ADB              mob;   // Phase mobility (per cell)
+        };
+
+        struct SolutionState {
+            SolutionState(const int np);
+            ADB             pressure;
+            std::vector<ADB> saturation;
+            ADB             concentration;
+            ADB             qs;
+            ADB             bhp;
+        };
+
+        struct WellOps {
+            WellOps(const Wells& wells);
+            M w2p;              // well -> perf (scatter)
+            M p2w;              // perf -> well (gather)
+        };
+
+        const UnstructuredGrid&         grid_;
+        const IncompPropsAdInterface&   fluid_;
+        const PolymerPropsAd&           polymer_props_ad_;
+        const LinearSolverInterface&    linsolver_;
+        const Wells&                    wells_;
+        const double*                   gravity_;
+        const std::vector<int>          cells_;
+        HelperOps                       ops_;
+        const WellOps                   wops_;
+      	V								cmax_;
+        std::vector<ReservoirResidualQuant> rq_;
+
+        struct {
+            std::vector<ADB>     mass_balance;
+            ADB                  well_eq;
+            ADB                  well_flux_eq;
+        } residual_;
+
+        SolutionState
+        constantState(const PolymerState& x,
+                      const WellState&    xw);
+        SolutionState
+        variableState(const PolymerState& x,
+                      const WellState&    xw);
+        void
+        assemble(const V&               pvdt,
+                 const PolymerState&    x,
+                 const WellState&       xw,
+                 const std::vector<double>& polymer_inflow,
+				 std::vector<double>& src);
+
+        V solveJacobianSystem() const;
+
+        void updateState(const V&             dx,
+                         PolymerState& x,
+                         WellState&    xw) const;
+        std::vector<ADB>
+        computeRelPerm(const SolutionState& state) const;
+
+        V
+        transmissibility() const;
+        
+		std::vector<ADB>
+		computePressures(const SolutionState& state) const;
+
+		void
+        computeMassFlux(const V&                trans,
+                        const ADB&              mc,
+                        const ADB&              kro,
+                        const ADB&              krw_eff,
+                        const SolutionState&    state );
+    
+        std::vector<ADB>
+        accumSource(const ADB&                 kro,
+                    const ADB&                 krw_eff,
+                    const ADB&                 c,
+                    const std::vector<double>& src,
+                    const std::vector<double>& polymer_inflow_c) const;
+
+        
+        std::vector<ADB>
+        computeFracFlow() const;
+
+        double
+        residualNorm() const;
+
+        ADB
+        polymerSource(const std::vector<ADB>& kr,
+                      const std::vector<double>& src,
+                      const std::vector<double>& polymer_inflow_c,
+                      const SolutionState& state) const;
+
+        void
+        computeCmax(PolymerState& state,
+					const ADB& c);
+    	void
+		computeAccum(const SolutionState& state,
+                 	const int            aix  );
+        ADB 
+        computeMc(const SolutionState& state) const;
+
+        ADB
+        rockPorosity(const ADB& p) const;
+
+        ADB
+        rockPermeability(const ADB& p) const;
+
+        const double
+        fluidDensity(const int phase) const;
+
+        ADB
+        fluidDensity(const int phase,
+                     const ADB p) const;
+
+        ADB
+        transMult(const ADB& p) const;
+    };
+
+
+} // namespace Opm
+
+#endif// OPM_FULLYIMPLICITTWOPHASESOLVER_HEADER_INCLUDED

opm/polymer/fullyimplicit/PolymerPropsAd.cpp

@@ -0,0 +1,272 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cmath>
+#include <vector>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+
+namespace Opm {
+
+    typedef PolymerPropsAd::ADB ADB;
+    typedef PolymerPropsAd::V V;
+
+
+
+
+
+    double
+    PolymerPropsAd::rockDensity() const
+    {
+        return polymer_props_.rockDensity();
+    }
+
+
+
+
+
+    double
+    PolymerPropsAd::deadPoreVol() const
+    {
+        return polymer_props_.deadPoreVol();
+    }
+
+
+
+
+
+	double
+	PolymerPropsAd::cMax() const
+	{
+		return polymer_props_.cMax();
+	}
+
+
+
+	
+
+    PolymerPropsAd::PolymerPropsAd(const PolymerProperties& polymer_props)
+        : polymer_props_ (polymer_props)
+    {
+    }
+
+
+
+
+
+    PolymerPropsAd::~PolymerPropsAd()
+    {
+    }
+
+
+
+
+
+    V PolymerPropsAd::effectiveInvWaterVisc(const V& c,
+                                            const double* visc) const
+    {
+        const int nc = c.size();
+        V inv_mu_w_eff(nc);
+        for (int i = 0; i < nc; ++i) {
+            double im = 0;
+            polymer_props_.effectiveInvVisc(c(i), visc, im);
+            inv_mu_w_eff(i) = im;
+        }
+
+        return inv_mu_w_eff;
+    }
+
+
+
+
+
+    ADB PolymerPropsAd::effectiveInvWaterVisc(const ADB& c,
+	                    				      const double* visc) const
+    {
+	    const int nc = c.size();
+    	V inv_mu_w_eff(nc);
+    	V dinv_mu_w_eff(nc);
+    	for (int i = 0; i < nc; ++i) {
+    	    double im = 0, dim = 0;
+    	    polymer_props_.effectiveInvViscWithDer(c.value()(i), visc, im, dim);
+    	    inv_mu_w_eff(i) = im;
+    	    dinv_mu_w_eff(i) = dim;
+    	}
+        ADB::M dim_diag = spdiag(dinv_mu_w_eff);
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dim_diag * c.derivative()[block];
+        }
+        return ADB::function(inv_mu_w_eff, jacs);
+    }
+
+
+
+
+
+    V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
+    {
+        const int nc = c.size();
+        V mc(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            double m = 0;
+            polymer_props_.computeMc(c(i), m);
+            mc(i) = m;
+        }
+       
+       return mc;
+    }
+
+
+
+
+
+    ADB PolymerPropsAd::polymerWaterVelocityRatio(const ADB& c) const
+    {
+    
+        const int nc = c.size();
+        V mc(nc);
+        V dmc(nc);
+        
+        for (int i = 0; i < nc; ++i) {
+            double m = 0;
+            double dm = 0;
+            polymer_props_.computeMcWithDer(c.value()(i), m, dm);
+
+            mc(i) = m;
+            dmc(i) = dm;
+        }
+
+        ADB::M dmc_diag = spdiag(dmc);
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dmc_diag * c.derivative()[block];
+        }
+
+        return ADB::function(mc, jacs);
+    }
+
+
+
+
+
+    V PolymerPropsAd::adsorption(const V& c, const V& cmax_cells) const
+    {
+        const int nc = c.size();
+        V ads(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            double c_ads = 0;
+            polymer_props_.adsorption(c(i), cmax_cells(i), c_ads);
+
+            ads(i) = c_ads;
+        }
+
+        return ads;
+    }
+
+
+
+
+
+    ADB PolymerPropsAd::adsorption(const ADB& c, const ADB& cmax_cells) const
+    {
+        const int nc = c.value().size();
+
+        V ads(nc);
+        V dads(nc);
+
+        for (int i = 0; i < nc; ++i) {
+            double c_ads = 0;
+            double dc_ads = 0;
+            polymer_props_.adsorptionWithDer(c.value()(i), cmax_cells.value()(i), c_ads, dc_ads);
+            ads(i) = c_ads;
+            dads(i) = dc_ads;
+        }
+
+        ADB::M dads_diag = spdiag(dads);
+        int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dads_diag * c.derivative()[block];
+        }
+
+        return ADB::function(ads, jacs);
+    }
+
+
+
+
+
+    V
+    PolymerPropsAd::effectiveRelPerm(const V& c, 
+                                     const V& cmax_cells,
+                                     const V& krw) const
+    {
+        const int nc = c.size();
+
+        V one  = V::Ones(nc);
+        V ads = adsorption(c, cmax_cells);
+        double max_ads = polymer_props_.cMaxAds();
+        double res_factor = polymer_props_.resFactor();
+        double factor = (res_factor -1.) / max_ads;
+        V rk = one + factor * ads;
+        V krw_eff = krw / rk;
+
+        return krw_eff;
+    }
+
+
+
+
+
+    ADB
+    PolymerPropsAd::effectiveRelPerm(const ADB& c,
+                                     const ADB& cmax_cells,
+                                     const ADB& krw,
+                                     const ADB& sw) const
+    {
+        const int nc = c.value().size();
+        V one = V::Ones(nc);
+        ADB ads = adsorption(c, cmax_cells);
+        V krw_eff = effectiveRelPerm(c.value(), cmax_cells.value(), krw.value());
+
+        double max_ads = polymer_props_.cMaxAds();
+        double res_factor = polymer_props_.resFactor();
+        double factor = (res_factor - 1.) / max_ads;
+        ADB rk = one + ads * factor; 
+		ADB dkrw_ds = krw / rk;
+        ADB dkrw_dc = -factor * krw / (rk * rk) * ads ;
+
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dkrw_ds.derivative()[block] * sw.derivative()[block] 
+						+ dkrw_dc.derivative()[block] * c.derivative()[block];
+        }
+
+        return ADB::function(krw_eff, jacs);
+    }
+
+}// namespace Opm

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -0,0 +1,112 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_POLYMERPROPSAD_HEADED_INLCUDED
+#define OPM_POLYMERPROPSAD_HEADED_INLCUDED
+
+#include <cmath>
+#include <vector>
+#include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+#include <opm/polymer/PolymerProperties.hpp>
+
+namespace Opm {
+    
+    class PolymerPropsAd 
+    {
+    public:
+		/// \return		Reference rock density.
+        double rockDensity() const;
+		
+		/// \return 	The value of dead pore volume.
+        double deadPoreVol() const;
+
+		/// \return 	The max concentration injected.
+       	double cMax() const; 
+
+		typedef AutoDiffBlock<double> ADB;
+        typedef ADB::V V;
+
+		/// Constructor wrapping a polymer props.	
+        PolymerPropsAd(const PolymerProperties& polymer_props);
+
+		/// Destructor.
+        ~PolymerPropsAd();
+		
+		/// \param[in] c		Array of n polymer concentraion values.
+		/// \param[in] visc		Array of 2 viscosity value.
+		/// \return 			value of inverse effective water viscosity.
+        V 
+        effectiveInvWaterVisc(const V& c,const double* visc) const;
+
+		/// \param[in] c		Array of n polymer concentraion values.
+		/// \param[in] visc		Array of 2 viscosity value 
+		/// \return 			value of inverse effective water viscosity.
+        ADB 
+        effectiveInvWaterVisc(const ADB& c,const double* visc) const;
+		
+		/// \param[in] c		Array of n polymer concentraion values.
+		/// \return 			Array of n mc values, here mc means m(c) * c.
+        V 
+        polymerWaterVelocityRatio(const V& c) const;
+
+		/// \param[in] c		Array of n polymer concentraion values.
+		/// \return 			Array of n mc values, here mc means m(c) * c.
+        ADB
+        polymerWaterVelocityRatio(const ADB& c) const;
+
+		/// \param[in] c				Array of n polymer concentraion values.
+		/// \param[in] cmax_cells		Array of n polymer concentraion values
+		///								that the cell experienced.
+		/// \return						Array of n adsorption values.
+        V
+        adsorption(const V& c, const V& cmax_cells) const;
+
+		/// \param[in] c				Array of n polymer concentraion values.
+		/// \param[in] cmax_cells		Array of n polymer concentraion values
+		///								that the cell experienced.
+		/// \return						Array of n adsorption values.
+        ADB
+        adsorption(const ADB& c, const ADB& cmax_cells) const;
+
+		/// \param[in] c				Array of n polymer concentraion values.
+		/// \param[in] cmax_cells		Array of n polymer concentraion values
+		///								that the cell experienced.
+		/// \param[in] relperm			Array of n relative water relperm values.
+		/// \return						Array of n adsorption values.
+        V
+        effectiveRelPerm(const V& c, const V& cmax_cells, const V& relperm) const;
+
+
+		/// \param[in] c				Array of n polymer concentraion values.
+		/// \param[in] cmax_cells		Array of n polymer concentraion values
+		///								that the cell experienced.
+		/// \param[in] relperm			Array of n relative water relperm values.
+		/// \return						Array of n adsorption values.
+        ADB
+        effectiveRelPerm(const ADB& c, const ADB& cmax_cells, const ADB& krw, const ADB& sw) const;
+
+    private:
+        const PolymerProperties& polymer_props_;
+    };
+    
+} //namespace Opm
+
+#endif// OPM_POLYMERPROPSAD_HEADED_INLCUDED

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp

@@ -0,0 +1,117 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED
+
+#include <memory>
+#include <vector>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class BlackoilPropsAdInterface;
+    class RockCompressibility;
+    class DerivedGeology;
+    class NewtonIterationBlackoilInterface;
+    class SimulatorTimer;
+    class PolymerBlackoilState;
+    class WellStateFullyImplicitBlackoil;
+    class WellsManager;
+    class EclipseState;
+    class EclipseWriter;
+    class PolymerPropsAd;
+    class PolymerInflowInterface;
+    struct SimulatorReport;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    template<class T>
+    class SimulatorFullyImplicitBlackoilPolymer
+    {
+    public:
+        /// \brief The type of the grid that we use.
+        typedef T Grid;
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] geo           derived geological properties
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] polymer_props polymer properties
+        /// \param[in] rock_comp_props     if non-null, rock compressibility properties
+        /// \param[in] linsolver     linear solver
+        /// \param[in] gravity       if non-null, gravity vector
+        /// \param[in] disgas        true for dissolved gas option
+        /// \param[in] vapoil        true for vaporized oil option
+        /// \param[in] eclipse_state
+        /// \param[in] output_writer
+        /// \param[in] threshold_pressures_by_face   if nonempty, threshold pressures that inhibit flow
+        SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
+                                              const Grid& grid,
+                                              const DerivedGeology& geo,
+                                              BlackoilPropsAdInterface& props,
+                                              const PolymerPropsAd& polymer_props,
+                                              const RockCompressibility* rock_comp_props,
+                                              NewtonIterationBlackoilInterface& linsolver,
+                                              const double* gravity,
+                                              const bool disgas,
+                                              const bool vapoil,
+                                              const bool polymer,
+                                              std::shared_ptr<EclipseState> eclipse_state,
+                                              EclipseWriter& output_writer,
+                                              Opm::DeckConstPtr& deck,
+                                              const std::vector<double>& threshold_pressures_by_face);
+
+        /// Run the simulation.
+        /// This will run succesive timesteps until timer.done() is true. It will
+        /// modify the reservoir and well states.
+        /// \param[in,out] timer       governs the requested reporting timesteps
+        /// \param[in,out] state       state of reservoir: pressure, fluxes
+        /// \param[in,out] well_state  state of wells: bhp, perforation rates
+        /// \return                    simulation report, with timing data
+        SimulatorReport run(SimulatorTimer& timer,
+                            PolymerBlackoilState& state);
+
+    private:
+        class Impl;
+        // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
+        std::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#include "SimulatorFullyImplicitBlackoilPolymer_impl.hpp"
+#endif // OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMER_HEADER_INCLUDED

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp

@@ -0,0 +1,211 @@
+#include "config.h"
+
+#include "SimulatorFullyImplicitBlackoilPolymerOutput.hpp"
+
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <opm/core/utility/platform_dependent/disable_warnings.h>
+#include <dune/common/version.hh>
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+#include <opm/core/utility/platform_dependent/reenable_warnings.h>
+#endif
+namespace Opm
+{
+
+
+    void outputStateVtk(const UnstructuredGrid& grid,
+                        const 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["concentration"] = &state.concentration();
+        dm["maxconcentration"] = &state.maxconcentration();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  AutoDiffGrid::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(grid),
+                                  state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+
+
+    void outputStateMatlab(const UnstructuredGrid& grid,
+                           const PolymerBlackoilState& state,
+                           const int step,
+                           const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        dm["surfvolume"] = &state.surfacevol();
+        dm["rs"] = &state.gasoilratio();
+        dm["rv"] = &state.rv();
+        dm["concentration"] = &state.concentration();
+        dm["maxconcentration"] = &state.maxconcentration();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  UgGridHelpers::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(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"));
+        }
+    }
+    void outputWellStateMatlab(const Opm::WellState& well_state,
+                               const int step,
+                               const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
+
+        // 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
+    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);
+    }
+
+    void outputWellReport(const Opm::WellReport& wellreport,
+                          const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+#endif
+
+#ifdef HAVE_DUNE_CORNERPOINT
+    void outputStateVtk(const Dune::CpGrid& grid,
+                        const PolymerBlackoilState& state,
+                        const int step,
+                        const std::string& output_dir)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        std::ostringstream vtkpath;
+        vtkpath << output_dir << "/vtk_files";
+        vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step;
+        boost::filesystem::path fpath(vtkpath.str());
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step;
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3)
+        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafGridView(), Dune::VTK::nonconforming);
+#else
+        Dune::VTKWriter<Dune::CpGrid::LeafGridView> writer(grid.leafView(), Dune::VTK::nonconforming);
+#endif
+        writer.addCellData(state.saturation(), "saturation", state.numPhases());
+        writer.addCellData(state.pressure(), "pressure", 1);
+        writer.addCellData(state.concentration(), "concentration", 1);
+        writer.addCellData(state.maxconcentration(), "maxconcentration", 1);
+        
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  AutoDiffGrid::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(grid),
+                                  state.faceflux(), cell_velocity);
+        writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension);
+        writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii);
+    }
+#endif
+
+}

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp

@@ -0,0 +1,96 @@
+#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED
+#include <opm/core/grid.h>
+#include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/DataMap.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/autodiff/GridHelpers.hpp>
+
+#include <string>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+
+#include <boost/filesystem.hpp>
+
+#ifdef HAVE_DUNE_CORNERPOINT
+#include <dune/grid/CpGrid.hpp>
+#endif
+namespace Opm
+{
+
+    void outputStateVtk(const UnstructuredGrid& grid,
+                        const PolymerBlackoilState& state,
+                        const int step,
+                        const std::string& output_dir);
+
+
+    void outputStateMatlab(const UnstructuredGrid& grid,
+                           const PolymerBlackoilState& state,
+                           const int step,
+                           const std::string& output_dir);
+
+    void outputWellStateMatlab(const Opm::WellState& well_state,
+                               const int step,
+                               const std::string& output_dir);
+#ifdef HAVE_DUNE_CORNERPOINT
+    void outputStateVtk(const Dune::CpGrid& grid,
+                        const PolymerBlackoilState& state,
+                        const int step,
+                        const std::string& output_dir);
+#endif
+
+    template<class Grid>
+    void outputStateMatlab(const Grid& grid,
+                           const PolymerBlackoilState& state,
+                           const int step,
+                           const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        dm["surfvolume"] = &state.surfacevol();
+        dm["rs"] = &state.gasoilratio();
+        dm["rv"] = &state.rv();
+        dm["concentration"] = &state.concentration();
+        dm["maxconcentration"] = &state.maxconcentration();
+
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid),
+                                  AutoDiffGrid::numFaces(grid),
+                                  AutoDiffGrid::beginFaceCentroids(grid),
+                                  UgGridHelpers::faceCells(grid),
+                                  AutoDiffGrid::beginCellCentroids(grid),
+                                  AutoDiffGrid::beginCellVolumes(grid),
+                                  AutoDiffGrid::dimensions(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"));
+        }
+    }
+
+
+}
+#endif

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp

@@ -0,0 +1,621 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp>
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp>
+#include <opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp>
+#include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/RateConverter.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/well_controls.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
+
+#include <opm/core/props/rock/RockCompressibility.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/lexical_cast.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <cassert>
+#include <functional>
+#include <memory>
+#include <numeric>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace Opm
+{
+    template<class T>
+    class SimulatorFullyImplicitBlackoilPolymer<T>::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const Grid& grid,
+             const DerivedGeology& geo,
+             BlackoilPropsAdInterface& props,
+             const PolymerPropsAd&  polymer_props,
+             const RockCompressibility* rock_comp_props,
+             NewtonIterationBlackoilInterface& linsolver,
+             const double* gravity,
+             bool has_disgas,
+             bool has_vapoil,
+             bool has_polymer,
+             std::shared_ptr<EclipseState> eclipse_state,
+             EclipseWriter& output_writer,
+             Opm::DeckConstPtr& deck,
+             const std::vector<double>& threshold_pressures_by_face);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            PolymerBlackoilState& state);
+
+    private:
+        // Data.
+        typedef RateConverter::
+        SurfaceToReservoirVoidage< BlackoilPropsAdInterface,
+                                   std::vector<int> > RateConverterType;
+
+        const parameter::ParameterGroup param_;
+
+        // Parameters for output.
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Observed objects.
+        const Grid& grid_;
+        BlackoilPropsAdInterface& props_;
+        const PolymerPropsAd& polymer_props_;
+        const RockCompressibility* rock_comp_props_;
+        const double* gravity_;
+        // Solvers
+        const DerivedGeology& geo_;
+        NewtonIterationBlackoilInterface& solver_;
+        // Misc. data
+        std::vector<int> allcells_;
+        const bool has_disgas_;
+        const bool has_vapoil_;
+        const bool has_polymer_;
+        // eclipse_state
+        std::shared_ptr<EclipseState> eclipse_state_;
+        // output_writer
+        EclipseWriter& output_writer_;
+        Opm::DeckConstPtr& deck_;
+        RateConverterType rateConverter_;
+        // Threshold pressures.
+        std::vector<double> threshold_pressures_by_face_;
+
+        void
+        computeRESV(const std::size_t               step,
+                    const Wells*                    wells,
+                    const PolymerBlackoilState&     x,
+                    WellStateFullyImplicitBlackoil& xw);
+    };
+
+
+
+
+    template<class T>
+    SimulatorFullyImplicitBlackoilPolymer<T>::SimulatorFullyImplicitBlackoilPolymer(const parameter::ParameterGroup& param,
+                                                                                    const Grid& grid,
+                                                                                    const DerivedGeology& geo,
+                                                                                    BlackoilPropsAdInterface& props,
+                                                                                    const PolymerPropsAd& polymer_props,
+                                                                                    const RockCompressibility* rock_comp_props,
+                                                                                    NewtonIterationBlackoilInterface& linsolver,
+                                                                                    const double* gravity,
+                                                                                    const bool has_disgas,
+                                                                                    const bool has_vapoil,
+                                                                                    const bool has_polymer,
+                                                                                    std::shared_ptr<EclipseState> eclipse_state,
+                                                                                    EclipseWriter& output_writer,
+                                                                                    Opm::DeckConstPtr& deck,
+                                                                                    const std::vector<double>& threshold_pressures_by_face)
+
+    {
+        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, 
+                              has_vapoil, has_polymer, eclipse_state, output_writer, deck, threshold_pressures_by_face));
+    }
+
+
+
+
+
+    template<class T>
+    SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::run(SimulatorTimer& timer,
+                                                                  PolymerBlackoilState& state)
+    {
+        return pimpl_->run(timer, state);
+    }
+
+
+
+    static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state,
+                                      const int step,
+                                      const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
+
+        // 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);
+    }
+
+    static void outputWellReport(const Opm::WellReport& wellreport,
+                                 const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+#endif
+
+
+    // \TODO: Treat bcs.
+    template<class T>
+    SimulatorFullyImplicitBlackoilPolymer<T>::Impl::Impl(const parameter::ParameterGroup& param,
+                                                         const Grid& grid,
+                                                         const DerivedGeology& geo,
+                                                         BlackoilPropsAdInterface& props,
+                                                         const PolymerPropsAd& polymer_props,
+                                                         const RockCompressibility* rock_comp_props,
+                                                         NewtonIterationBlackoilInterface& linsolver,
+                                                         const double* gravity,
+                                                         const bool has_disgas,
+                                                         const bool has_vapoil,
+                                                         const bool has_polymer,
+                                                         std::shared_ptr<EclipseState> eclipse_state,
+                                                         EclipseWriter& output_writer,
+                                                         Opm::DeckConstPtr& deck,
+                                                         const std::vector<double>& threshold_pressures_by_face)
+        : param_(param),
+          grid_(grid),
+          props_(props),
+          polymer_props_(polymer_props),
+          rock_comp_props_(rock_comp_props),
+          gravity_(gravity),
+          geo_(geo),
+          solver_(linsolver),
+          has_disgas_(has_disgas),
+          has_vapoil_(has_vapoil),
+          has_polymer_(has_polymer),
+          eclipse_state_(eclipse_state),
+          output_writer_(output_writer),
+          deck_(deck),
+          rateConverter_(props_, std::vector<int>(AutoDiffGrid::numCells(grid_), 0)),
+          threshold_pressures_by_face_(threshold_pressures_by_face)
+    {
+        // 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);
+        }
+
+        // Misc init.
+        const int num_cells = AutoDiffGrid::numCells(grid);
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+
+
+
+    template<class T>
+    SimulatorReport SimulatorFullyImplicitBlackoilPolymer<T>::Impl::run(SimulatorTimer& timer,
+                                                                        PolymerBlackoilState& state)
+    {
+        WellStateFullyImplicitBlackoil prev_well_state;
+
+        // Create timers and file for writing timing info.
+        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()) {
+            // Report timestep.
+            step_timer.start();
+            timer.report(std::cout);
+
+            // Create wells and well state.
+            WellsManager wells_manager(eclipse_state_,
+                                       timer.currentStepNum(),
+                                       Opm::UgGridHelpers::numCells(grid_),
+                                       Opm::UgGridHelpers::globalCell(grid_),
+                                       Opm::UgGridHelpers::cartDims(grid_),
+                                       Opm::UgGridHelpers::dimensions(grid_),
+                                       Opm::UgGridHelpers::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(), 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_, *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);
+            // Output state at start of time step.
+            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+                if (output_vtk_) {
+                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+                }
+                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+            }
+            if (output_) {
+                if (timer.currentStepNum() == 0) {
+                    output_writer_.writeInit(timer);
+                }
+                output_writer_.writeTimeStep(timer, state.blackoilState(), well_state.basicWellState());
+            }
+
+            // Max oil saturation (for VPPARS), hysteresis update.
+            props_.updateSatOilMax(state.saturation());
+            props_.updateSatHyst(state.saturation(), allcells_);
+
+            // Compute reservoir volumes for RESV controls.
+            computeRESV(timer.currentStepNum(), wells, state, well_state);
+
+            // Run a single step of the solver.
+            solver_timer.start();
+            FullyImplicitBlackoilPolymerSolver<T> solver(param_, grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells, solver_, has_disgas_, has_vapoil_, has_polymer_);
+            if (!threshold_pressures_by_face_.empty()) {
+                solver.setThresholdPressures(threshold_pressures_by_face_);
+            }
+            solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
+            solver_timer.stop();
+
+            // Report timing.
+            const double st = solver_timer.secsSinceStart();
+            std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
+            stime += st;
+            if (output_) {
+                SimulatorReport step_report;
+                step_report.pressure_time = st;
+                step_report.total_time =  step_timer.secsSinceStart();
+                step_report.reportParam(tstep_os);
+            }
+
+            // Increment timer, remember well state.
+            ++timer;
+            prev_well_state = well_state;
+        }
+
+        // Write final simulation state.
+        if (output_) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_);
+            output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState());
+        } 
+
+        // Stop timer and create timing report
+        total_timer.stop();
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        return report;
+    }
+
+    namespace SimFIBODetails {
+        typedef std::unordered_map<std::string, WellConstPtr> WellMap;
+
+        inline WellMap
+        mapWells(const std::vector<WellConstPtr>& wells)
+        {
+            WellMap wmap;
+
+            for (std::vector<WellConstPtr>::const_iterator
+                     w = wells.begin(), e = wells.end();
+                 w != e; ++w)
+            {
+                wmap.insert(std::make_pair((*w)->name(), *w));
+            }
+
+            return wmap;
+        }
+
+        inline int
+        resv_control(const WellControls* ctrl)
+        {
+            int i, n = well_controls_get_num(ctrl);
+
+            bool match = false;
+            for (i = 0; (! match) && (i < n); ++i) {
+                match = well_controls_iget_type(ctrl, i) == RESERVOIR_RATE;
+            }
+
+            if (! match) { i = 0; }
+
+            return i - 1; // -1 if no match, undo final "++" otherwise
+        }
+
+        inline bool
+        is_resv_prod(const Wells& wells,
+                     const int    w)
+        {
+            return ((wells.type[w] == PRODUCER) &&
+                    (0 <= resv_control(wells.ctrls[w])));
+        }
+
+        inline bool
+        is_resv_prod(const WellMap&     wmap,
+                     const std::string& name,
+                     const std::size_t  step)
+        {
+            bool match = false;
+
+            WellMap::const_iterator i = wmap.find(name);
+
+            if (i != wmap.end()) {
+                WellConstPtr wp = i->second;
+
+                match = (wp->isProducer(step) &&
+                         wp->getProductionProperties(step)
+                         .hasProductionControl(WellProducer::RESV));
+            }
+
+            return match;
+        }
+
+        inline std::vector<int>
+        resvProducers(const Wells&      wells,
+                      const std::size_t step,
+                      const WellMap&    wmap)
+        {
+            std::vector<int> resv_prod;
+
+            for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) {
+                if (is_resv_prod(wells, w) ||
+                    ((wells.name[w] != 0) &&
+                     is_resv_prod(wmap, wells.name[w], step)))
+                {
+                    resv_prod.push_back(w);
+                }
+            }
+
+            return resv_prod;
+        }
+
+        inline void
+        historyRates(const PhaseUsage&               pu,
+                     const WellProductionProperties& p,
+                     std::vector<double>&            rates)
+        {
+            assert (! p.predictionMode);
+            assert (rates.size() ==
+                    std::vector<double>::size_type(pu.num_phases));
+
+            if (pu.phase_used[ BlackoilPhases::Aqua ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Aqua ];
+
+                rates[i] = p.WaterRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Liquid ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Liquid ];
+
+                rates[i] = p.OilRate;
+            }
+
+            if (pu.phase_used[ BlackoilPhases::Vapour ]) {
+                const std::vector<double>::size_type
+                    i = pu.phase_pos[ BlackoilPhases::Vapour ];
+
+                rates[i] = p.GasRate;
+            }
+        }
+    } // namespace SimFIBODetails
+
+    template <class T>
+    void
+    SimulatorFullyImplicitBlackoilPolymer<T>::
+    Impl::computeRESV(const std::size_t               step,
+                      const Wells*                    wells,
+                      const PolymerBlackoilState&     x,
+                      WellStateFullyImplicitBlackoil& xw)
+    {
+        typedef SimFIBODetails::WellMap WellMap;
+
+        const std::vector<WellConstPtr>& w_ecl = eclipse_state_->getSchedule()->getWells(step);
+        const WellMap& wmap = SimFIBODetails::mapWells(w_ecl);
+
+        const std::vector<int>& resv_prod =
+            SimFIBODetails::resvProducers(*wells, step, wmap);
+
+        if (! resv_prod.empty()) {
+            const PhaseUsage&                    pu = props_.phaseUsage();
+            const std::vector<double>::size_type np = props_.numPhases();
+
+            rateConverter_.defineState(x.blackoilState());
+
+            std::vector<double> distr (np);
+            std::vector<double> hrates(np);
+            std::vector<double> prates(np);
+
+            for (std::vector<int>::const_iterator
+                     rp = resv_prod.begin(), e = resv_prod.end();
+                 rp != e; ++rp)
+            {
+                WellControls* ctrl = wells->ctrls[*rp];
+
+                // RESV control mode, all wells
+                {
+                    const int rctrl = SimFIBODetails::resv_control(ctrl);
+
+                    if (0 <= rctrl) {
+                        const std::vector<double>::size_type off = (*rp) * np;
+
+                        // Convert to positive rates to avoid issues
+                        // in coefficient calculations.
+                        std::transform(xw.wellRates().begin() + (off + 0*np),
+                                       xw.wellRates().begin() + (off + 1*np),
+                                       prates.begin(), std::negate<double>());
+
+                        const int fipreg = 0; // Hack.  Ignore FIP regions.
+                        rateConverter_.calcCoeff(prates, fipreg, distr);
+
+                        well_controls_iset_distr(ctrl, rctrl, & distr[0]);
+                    }
+                }
+
+                // RESV control, WCONHIST wells.  A bit of duplicate
+                // work, regrettably.
+                if (wells->name[*rp] != 0) {
+                    WellMap::const_iterator i = wmap.find(wells->name[*rp]);
+
+                    if (i != wmap.end()) {
+                        WellConstPtr wp = i->second;
+
+                        const WellProductionProperties& p =
+                            wp->getProductionProperties(step);
+
+                        if (! p.predictionMode) {
+                            // History matching (WCONHIST/RESV)
+                            SimFIBODetails::historyRates(pu, p, hrates);
+
+                            const int fipreg = 0; // Hack.  Ignore FIP regions.
+                            rateConverter_.calcCoeff(hrates, fipreg, distr);
+
+                            // WCONHIST/RESV target is sum of all
+                            // observed phase rates translated to
+                            // reservoir conditions.  Recall sign
+                            // convention: Negative for producers.
+                            const double target =
+                                - std::inner_product(distr.begin(), distr.end(),
+                                                     hrates.begin(), 0.0);
+
+                            well_controls_clear(ctrl);
+                            well_controls_assert_number_of_phases(ctrl, int(np));
+
+                            const int ok =
+                                well_controls_add_new(RESERVOIR_RATE, target,
+                                                      & distr[0], ctrl);
+
+                            if (ok != 0) {
+                                xw.currentControls()[*rp] = 0;
+                                well_controls_set_current(ctrl, 0);
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+} // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.cpp

@@ -0,0 +1,474 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/autodiff/GeoProps.hpp>
+#include <opm/autodiff/BlackoilPropsAdInterface.hpp>
+#include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+
+#include <opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
+
+#include <opm/core/wells/WellsManager.hpp>
+
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/grid/ColumnExtract.hpp>
+#include <opm/polymer/PolymerBlackoilState.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+#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>
+
+#include <numeric>
+#include <fstream>
+#include <iostream>
+
+
+namespace Opm
+{
+
+
+
+    namespace 
+    {
+        static void outputStateVtk(const UnstructuredGrid& grid,
+                                   const Opm::PolymerBlackoilState& state,
+                                   const int step,
+                                   const std::string& output_dir);
+        static void outputStateMatlab(const UnstructuredGrid& grid,
+                                      const Opm::PolymerBlackoilState& state,
+                                      const int step,
+                                      const std::string& output_dir);
+        static  void outputWaterCut(const Opm::Watercut& watercut,
+                	            const std::string& output_dir);
+    } // anonymous namespace
+
+    class SimulatorFullyImplicitCompressiblePolymer::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const DerivedGeology& geo,
+             const BlackoilPropsAdInterface& props,
+             const PolymerPropsAd&          polymer_props,
+             const RockCompressibility* rock_comp_props,
+             std::shared_ptr<EclipseState> eclipse_state,
+             EclipseWriter& output_writer,
+             Opm::DeckConstPtr& deck,
+             NewtonIterationBlackoilInterface& linsolver,
+             const double* gravity);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            PolymerBlackoilState& state);
+
+    private:
+        // Data.
+
+        // Parameters for output.
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for well control
+        bool check_well_controls_;
+        int max_well_control_iterations_;
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const BlackoilPropsAdInterface& props_;
+        const PolymerPropsAd&       polymer_props_;
+        const RockCompressibility* rock_comp_props_;
+        std::shared_ptr<EclipseState> eclipse_state_;
+        EclipseWriter& output_writer_;
+        Opm::DeckConstPtr& deck_;
+        NewtonIterationBlackoilInterface& linsolver_;
+        const double* gravity_;
+        // Solvers
+        DerivedGeology geo_;
+        // Misc. data
+        std::vector<int> allcells_;
+    };
+
+
+
+
+    SimulatorFullyImplicitCompressiblePolymer::
+    SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
+                                              const UnstructuredGrid& grid,
+                                              const DerivedGeology& geo,
+                                              const BlackoilPropsAdInterface& props,
+                                              const PolymerPropsAd& polymer_props,
+                                              const RockCompressibility* rock_comp_props,
+                                              std::shared_ptr<EclipseState> eclipse_state,
+                                              EclipseWriter& output_writer,
+                                              Opm::DeckConstPtr& deck,
+                                              NewtonIterationBlackoilInterface& linsolver,
+                                              const double* gravity)
+
+    {
+        pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, eclipse_state, output_writer, deck, linsolver, gravity));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorFullyImplicitCompressiblePolymer::run(SimulatorTimer& timer,
+                                                        PolymerBlackoilState& state)
+    {
+        return pimpl_->run(timer, state);
+    }
+
+
+
+
+
+    // \TODO: Treat bcs.
+    SimulatorFullyImplicitCompressiblePolymer::Impl::Impl(const parameter::ParameterGroup& param,
+                   			                              const UnstructuredGrid& grid,
+                                                          const DerivedGeology& geo,
+                              			                  const BlackoilPropsAdInterface& props,
+                                         			      const PolymerPropsAd& polymer_props,
+                                  			              const RockCompressibility* rock_comp_props,
+                                                          std::shared_ptr<EclipseState> eclipse_state,
+                                                          EclipseWriter& output_writer,
+                                                          Opm::DeckConstPtr& deck,
+                         			                      NewtonIterationBlackoilInterface& linsolver,
+                                    		              const double* gravity)
+        : grid_(grid),
+          props_(props),
+          polymer_props_(polymer_props),
+          rock_comp_props_(rock_comp_props),
+          eclipse_state_(eclipse_state),
+          output_writer_(output_writer),
+          deck_(deck),
+          linsolver_(linsolver),
+          gravity_(gravity),
+          geo_(geo)
+    {
+        // For output.
+        output_ = param.getDefault("output", true);
+        if (output_) {
+            output_vtk_ = param.getDefault("output_vtk", true);
+            output_dir_ = param.getDefault("output_dir", std::string("output"));
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(output_dir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            output_interval_ = param.getDefault("output_interval", 1);
+        }
+
+        // Well control related init.
+        check_well_controls_ = param.getDefault("check_well_controls", false);
+        max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
+
+        // Misc init.
+        const int num_cells = grid.number_of_cells;
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+
+
+
+    SimulatorReport SimulatorFullyImplicitCompressiblePolymer::Impl::run(SimulatorTimer& timer,
+                                                                         PolymerBlackoilState& state)
+    {
+        WellStateFullyImplicitBlackoil prev_well_state;
+        // Initialisation.
+        std::vector<double> porevol;
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+        } else {
+            computePorevolume(grid_, props_.porosity(), porevol);
+        }
+        std::vector<double> initial_porevol = porevol;
+
+        std::vector<double> polymer_inflow_c(grid_.number_of_cells);
+        // Main simulation loop.
+        Opm::time::StopWatch solver_timer;
+        double stime = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        std::string tstep_filename = output_dir_ + "/step_timing.txt";
+        std::ofstream tstep_os(tstep_filename.c_str());
+
+        //Main simulation loop.
+        while (!timer.done()) {
+#if 0
+            double tot_injected[2] = { 0.0 };
+            double tot_produced[2] = { 0.0 };
+            Opm::Watercut watercut;
+            watercut.push(0.0, 0.0, 0.0);
+            std::vector<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(eclipse_state_,
+                                       timer.currentStepNum(),
+                                       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(), 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_, *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.
+            solver_timer.start();
+            FullyImplicitCompressiblePolymerSolver solver(grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells_manager.c_wells(), linsolver_);
+            solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c);
+            // Stop timer and report.
+            solver_timer.stop();
+            const double st = solver_timer.secsSinceStart();
+            std::cout << "Fully implicit solver took:  " << st << " seconds." << std::endl;
+
+            stime += st;
+            // Update pore volumes if rock is compressible.
+            if (rock_comp_props_ && rock_comp_props_->isActive()) {
+                initial_porevol = porevol;
+                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+            }
+/*
+            double injected[2] = { 0.0 };
+            double produced[2] = { 0.0 };
+    		double polyinj = 0;
+    		double polyprod = 0;
+            Opm::computeInjectedProduced(props_, polymer_props_,
+                                         state,
+                                         transport_src, polymer_inflow_c, timer.currentStepLength(),
+                                         injected, produced,
+                                         polyinj, polyprod);
+            tot_injected[0] += injected[0];
+            tot_injected[1] += injected[1];
+            tot_produced[0] += produced[0];
+            tot_produced[1] += produced[1];
+            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
+                          	  produced[0]/(produced[0] + produced[1]),
+                          	  tot_produced[0]/tot_porevol_init);
+            std::cout.precision(5);
+            const int width = 18;
+            std::cout << "\nMass balance report.\n";
+            std::cout << "    Injected reservoir volumes:      "
+                      << std::setw(width) << injected[0]
+                      << std::setw(width) << injected[1] << std::endl;
+            std::cout << "    Produced reservoir volumes:      "
+                      << std::setw(width) << produced[0]
+                      << std::setw(width) << produced[1] << std::endl;
+            std::cout << "    Total inj reservoir volumes:     "
+                      << std::setw(width) << tot_injected[0]
+                      << std::setw(width) << tot_injected[1] << std::endl;
+            std::cout << "    Total prod reservoir volumes:    "
+                      << std::setw(width) << tot_produced[0]
+                      << std::setw(width) << tot_produced[1] << std::endl;
+*/
+            if (output_) {
+                SimulatorReport step_report;
+                step_report.pressure_time = st;
+                step_report.total_time =  step_timer.secsSinceStart();
+                step_report.reportParam(tstep_os);
+            }
+            ++timer;
+            prev_well_state = well_state;
+        }
+        // Write final simulation state.
+        if (output_) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state.basicWellState());
+        }
+
+        total_timer.stop();
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        return report;
+    }
+
+
+
+    namespace
+    {
+
+        static void outputStateVtk(const UnstructuredGrid& grid,
+                                   const Opm::PolymerBlackoilState& state,
+                                   const int step,
+                                   const std::string& output_dir)
+        {
+            // Write data in VTK format.
+            std::ostringstream vtkfilename;
+            vtkfilename << output_dir << "/vtk_files";
+            boost::filesystem::path fpath(vtkfilename.str());
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+            std::ofstream vtkfile(vtkfilename.str().c_str());
+            if (!vtkfile) {
+                OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
+            }
+            Opm::DataMap dm;
+            dm["saturation"] = &state.saturation();
+            dm["pressure"] = &state.pressure();
+            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.concentration();
+            std::vector<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
+    }
+
+} // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp

@@ -0,0 +1,106 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED
+
+#include <memory>
+#include <vector>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+
+struct UnstructuredGrid;
+struct Wells;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class BlackoilPropsAdInterface;
+    class RockCompressibility;
+    class DerivedGeology;
+    class WellStateFullyImplicitBlackoil;
+    class WellsManager;
+    class EclipseWriter;
+    class EclipseState;
+    class NewtonIterationBlackoilInterface;
+    class SimulatorTimer;
+    class PolymerBlackoilState;
+    class PolymerPropsAd;
+    class PolymerInflowInterface;
+    struct SimulatorReport;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorFullyImplicitCompressiblePolymer
+    {
+    public:
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] polymer_props polymer properties
+        /// \param[in] rock_comp_props if non-null, rock compressibility properties
+        /// \param[in] eclipse_state  
+        /// \param[in] eclipse_writer
+        /// \param[in] deck
+        /// \param[in] linsolver     linear solver
+        /// \param[in] gravity       if non-null, gravity vector
+        SimulatorFullyImplicitCompressiblePolymer(const parameter::ParameterGroup& param,
+                        		                  const UnstructuredGrid& grid,
+                                                  const DerivedGeology& geo,
+                                   				  const BlackoilPropsAdInterface& props,
+                                       			  const PolymerPropsAd&    polymer_props,
+                                       			  const RockCompressibility* rock_comp_props,
+                                                  std::shared_ptr<EclipseState> eclipse_state,
+                                                  EclipseWriter& eclipse_writer,
+                                                  Opm::DeckConstPtr& deck,
+                                       			  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,
+                            PolymerBlackoilState& state);
+
+    private:
+        class Impl;
+        // Using shared_ptr instead of scoped_ptr since scoped_ptr requires complete type for Impl.
+        std::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORFULLYIMPLICITCOMPRESSIBLEPOLYMER_HEADER_INCLUDED

opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.cpp

@@ -0,0 +1,513 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/polymer/fullyimplicit/FullyImplicitTwophasePolymerSolver.hpp>
+#include <opm/autodiff/IncompPropsAdInterface.hpp>
+#include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
+#include <opm/polymer/fullyimplicit/utilities.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/io/vtk/writeVtkData.hpp>
+
+
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/core/grid/ColumnExtract.hpp>
+#include <opm/polymer/PolymerState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/polymer/PolymerInflow.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/lexical_cast.hpp>
+
+#include <numeric>
+#include <fstream>
+#include <iostream>
+#include <Eigen/Eigen>
+namespace Opm
+{
+
+
+
+    class SimulatorFullyImplicitTwophasePolymer::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const IncompPropsAdInterface& props,
+             const PolymerPropsAd&  polymer_props,
+             LinearSolverInterface& linsolver,
+             WellsManager&    wells_manager,
+             PolymerInflowInterface& polymer_inflow,
+             const double* gravity);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            PolymerState& state,
+                            WellState&    well_state);
+
+    private:
+
+        // Parameters for output.
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for well control
+        bool check_well_controls_;
+        int max_well_control_iterations_;
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const IncompPropsAdInterface& props_;
+        const PolymerPropsAd&  polymer_props_;
+        WellsManager&   wells_manager_;
+        const Wells*    wells_;
+        PolymerInflowInterface& polymer_inflow_;
+        // Solvers
+        FullyImplicitTwophasePolymerSolver solver_;
+        // Misc. data
+        std::vector<int> allcells_;
+    };
+
+
+
+
+    SimulatorFullyImplicitTwophasePolymer::
+    SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param,
+                                          const UnstructuredGrid& grid,
+                                          const IncompPropsAdInterface& props,
+                                          const PolymerPropsAd& polymer_props,
+                                          LinearSolverInterface& linsolver,
+                                          WellsManager& wells_manager,
+                                          PolymerInflowInterface&  polymer_inflow,
+                                          const double* gravity)
+    {
+        pimpl_.reset(new Impl(param, grid, props, polymer_props, linsolver, wells_manager, polymer_inflow, gravity));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorFullyImplicitTwophasePolymer::run(SimulatorTimer& timer,
+                                                        PolymerState& state,
+                                                        WellState&    well_state)
+    {
+        return pimpl_->run(timer, state, well_state);
+    }
+
+
+
+    static void outputStateVtk(const UnstructuredGrid& grid,
+                               const Opm::PolymerState& 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::PolymerState& 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();
+        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"));
+        }
+    }
+
+
+    
+    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);
+    }
+/*
+    static void outputWellReport(const Opm::WellReport& wellreport,
+                                 const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+  */  
+    
+    static void outputWellStateMatlab(WellState& well_state,
+                                      const int step,
+                                      const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["bhp"] = &well_state.bhp();
+        dm["wellrates"] = &well_state.wellRates();
+
+        // 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"));
+        }
+    }
+
+
+    
+    
+    SimulatorFullyImplicitTwophasePolymer::Impl::Impl(const parameter::ParameterGroup& param,
+                                               const UnstructuredGrid& grid,
+                                               const IncompPropsAdInterface& props,
+                                               const PolymerPropsAd&    polymer_props,
+                                               LinearSolverInterface& linsolver,
+                                               WellsManager&  wells_manager,
+                                               PolymerInflowInterface& polymer_inflow,
+                                               const double* gravity)
+        : grid_(grid),
+          props_(props),
+          polymer_props_(polymer_props),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
+          polymer_inflow_(polymer_inflow),
+          solver_(grid_, props_, polymer_props_, linsolver, *wells_manager.c_wells(), gravity)
+
+    {
+        // 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);
+        }
+        // Misc init.
+        const int num_cells = grid.number_of_cells;
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+    SimulatorReport SimulatorFullyImplicitTwophasePolymer::Impl::run(SimulatorTimer& timer,
+                                                              PolymerState& state,
+                                                              WellState&    well_state)
+    {
+
+        // Initialisation.
+        std::vector<double> porevol;
+        Opm::computePorevolume(grid_, props_.porosity(), porevol);
+
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        std::vector<double> polymer_inflow_c(grid_.number_of_cells);
+		std::vector<double> transport_src(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();
+        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.
+        double init_surfvol[2] = { 0.0 };
+        double inplace_surfvol[2] = { 0.0 };
+        double tot_injected[2] = { 0.0 };
+        double tot_produced[2] = { 0.0 };
+        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
+        Opm::Watercut watercut;
+        watercut.push(0.0, 0.0, 0.0);
+#endif
+        std::vector<double> fractional_flows;
+        std::vector<double> well_resflows_phase;
+        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);
+        }
+        while (!timer.done()) {
+            // 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_);
+        //        outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+
+            }
+
+            SimulatorReport sreport;
+
+            bool well_control_passed = !check_well_controls_;
+            int well_control_iteration = 0;
+            do {
+                // Run solver.
+                const double current_time = timer.currentTime();
+                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();
+                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;
+
+                // Optionally, check if well controls are satisfied.
+                if (check_well_controls_) {
+                    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);
+         	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.currentTime() + 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;
+
+
+            // Update pore volumes if rock is compressible.
+
+            // The reports below are geared towards two phases only.
+#if 0
+            // Report mass balances.
+            double injected[2] = { 0.0 };
+            double produced[2] = { 0.0 };
+            Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
+                                         injected, produced);
+            Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
+            tot_injected[0] += injected[0];
+            tot_injected[1] += injected[1];
+            tot_produced[0] += produced[0];
+            tot_produced[1] += produced[1];
+            std::cout.precision(5);
+            const int width = 18;
+            std::cout << "\nMass balance report.\n";
+            std::cout << "    Injected surface volumes:      "
+                      << std::setw(width) << injected[0]
+                      << std::setw(width) << injected[1] << std::endl;
+            std::cout << "    Produced surface volumes:      "
+                      << std::setw(width) << produced[0]
+                      << std::setw(width) << produced[1] << std::endl;
+            std::cout << "    Total inj surface volumes:     "
+                      << std::setw(width) << tot_injected[0]
+                      << std::setw(width) << tot_injected[1] << std::endl;
+            std::cout << "    Total prod surface volumes:    "
+                      << std::setw(width) << tot_produced[0]
+                      << std::setw(width) << tot_produced[1] << std::endl;
+            const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
+                                        init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] };
+            std::cout << "    Initial - inplace + inj - prod: "
+                      << std::setw(width) << balance[0]
+                      << std::setw(width) << balance[1]
+                      << std::endl;
+            std::cout << "    Relative mass error:            "
+                      << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
+                      << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
+                      << std::endl;
+            std::cout.precision(8);
+
+            // Make well reports.
+            watercut.push(timer.currentTime() + timer.currentStepLength(),
+                          produced[0]/(produced[0] + produced[1]),
+                          tot_produced[0]/tot_porevol_init);
+            if (wells_) {
+                wellreport.push(props_, *wells_,
+                                state.pressure(), state.surfacevol(), state.saturation(),
+                                timer.currentTime() + timer.currentStepLength(),
+                                well_state.bhp(), well_state.perfRates());
+            }
+#endif
+            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_);
+#if 0
+                outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
+                if (wells_) {
+                    outputWellReport(wellreport, output_dir_);
+                }
+#endif
+                tstep_os.close();
+            }
+
+            // advance to next timestep before reporting at this location
+            ++timer;
+
+            // write an output file for later inspection
+        }
+
+        total_timer.stop();
+
+        SimulatorReport report;
+        report.pressure_time = stime;
+        report.transport_time = 0.0;
+        report.total_time = total_timer.secsSinceStart();
+        return report;
+    }
+
+
+
+
+
+} // namespace Opm

opm/polymer/fullyimplicit/SimulatorFullyImplicitTwophasePolymer.hpp

@@ -0,0 +1,97 @@
+/*
+  Copyright 2014 SINTEF ICT, Applied Mathematics.
+  Copyright 2014 STATOIL ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED
+#define OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+
+struct UnstructuredGrid;
+struct Wells;
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class IncompPropsAdInterface;
+    class LinearSolverInterface;
+    class SimulatorTimer;
+    class PolymerState;
+    class PolymerPropsAd;
+    class PolymerInflowInterface;
+    class WellsManager;
+    class WellState;
+    struct SimulatorReport;
+    
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorFullyImplicitTwophasePolymer
+    {
+    public:
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] polymer_props polymer properties
+        /// \param[in] linsolver     linear solver
+        /// \param[in] well_manager  well manager, may manage no (null) wells
+        /// \param[in] polymer_inflow polymer influx.
+        /// \param[in] gravity       if non-null, gravity vector
+        SimulatorFullyImplicitTwophasePolymer(const parameter::ParameterGroup& param,
+        	                                  const UnstructuredGrid& grid,
+               		                          const IncompPropsAdInterface& props,
+                      		                  const PolymerPropsAd&    polymer_props,
+                             		          LinearSolverInterface& linsolver,
+                                    		  WellsManager&      wells_manager,
+	                                          PolymerInflowInterface& polymer_inflow,
+       		                                  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
+        /// \param[in,out] well_state  state of wells: bhp, perforation rates
+        /// \return                    simulation report, with timing data
+        SimulatorReport run(SimulatorTimer& timer,
+                            PolymerState& 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_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORFULLYIMPLICITTWOPHASEPOLYMER_HEADER_INCLUDED