Merge remote-tracking branch 'remotes/opm/master' into change-PRT-folder

CMakeLists_files.cmake

@@ -47,6 +47,9 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/autodiff/VFPProdProperties.cpp
 	opm/autodiff/VFPInjProperties.cpp
 	opm/autodiff/WellMultiSegment.cpp
+        opm/autodiff/BlackoilSolventState.cpp
+        opm/polymer/PolymerState.cpp
+        opm/polymer/PolymerBlackoilState.cpp
 	opm/polymer/CompressibleTpfaPolymer.cpp
 	opm/polymer/IncompTpfaPolymer.cpp
 	opm/polymer/PolymerInflow.cpp

examples/opm_init_check.cpp

@@ -31,7 +31,7 @@
 #include <ert/ecl/ecl_nnc_export.h>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/TransMult.hpp>
@@ -345,10 +345,10 @@ int main(int argc, char** argv) {
     
     ParserPtr parser(new Parser());
 
-    ParseMode parseMode;
+    ParseContext parseContext;
     std::cout << "Parsing input file ............: " << input_file << std::endl;
-    DeckConstPtr deck = parser->parseFile(input_file, parseMode);
-    std::shared_ptr<EclipseState> state = std::make_shared<EclipseState>( deck , parseMode );
+    DeckConstPtr deck = parser->parseFile(input_file, parseContext);
+    std::shared_ptr<EclipseState> state = std::make_shared<EclipseState>( deck , parseContext );
     
     std::cout << "Loading eclipse INIT file .....: " << init_file << std::endl;
     ecl_file_type * ecl_init = ecl_file_open( init_file.c_str() , 0 );

examples/sim_2p_incomp_ad.cpp

@@ -47,7 +47,7 @@
 #include <opm/autodiff/SimulatorIncompTwophaseAd.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
 #include <boost/filesystem.hpp>
@@ -105,32 +105,39 @@ try
     std::unique_ptr<GridManager> grid;
     std::unique_ptr<IncompPropertiesInterface> props;
     std::unique_ptr<RockCompressibility> rock_comp;
+    std::unique_ptr<TwophaseState> state;
     EclipseStateConstPtr eclipseState;
-    TwophaseState state;
+
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     if (use_deck) {
         std::string deck_filename = param.get<std::string>("deck_filename");
-        Opm::ParseMode parseMode;
-        deck = parser->parseFile(deck_filename, parseMode);
-        eclipseState.reset(new EclipseState(deck , parseMode));
+        Opm::ParseContext parseContext;
+        deck = parser->parseFile(deck_filename, parseContext);
+        eclipseState.reset(new EclipseState(deck , parseContext));
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // 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);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init
+            props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+
+            state.reset( new TwophaseState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
+            // 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(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
         }
     } else {
         // Grid init.
@@ -141,14 +148,20 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+
+            state.reset( new TwophaseState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid )));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+        }
     }
 
     // Warn if gravity but no density difference.
@@ -170,7 +183,7 @@ try
         // terms of total pore volume.
         std::vector<double> porevol;
         if (rock_comp->isActive()) {
-            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
         } else {
             computePorevolume(*grid->c_grid(), props->porosity(), porevol);
         }
@@ -236,8 +249,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more report steps etc.
         WellState well_state;
@@ -255,7 +268,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             simtimer.setCurrentStepNum(reportStepIdx);
@@ -273,7 +286,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
             if (output) {
                 epoch_rep.reportParam(epoch_os);
             }

examples/sim_poly2p_comp_reorder.cpp

@@ -48,7 +48,7 @@
 #include <opm/polymer/PolymerProperties.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
 #include <boost/scoped_ptr.hpp>
@@ -92,7 +92,7 @@ try
     boost::scoped_ptr<RockCompressibility> rock_comp;
     Opm::DeckConstPtr deck;
     EclipseStateConstPtr eclipseState;
-    PolymerBlackoilState state;
+    std::unique_ptr<PolymerBlackoilState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -100,29 +100,35 @@ try
     if (use_deck) {
         std::string deck_filename = param.get<std::string>("deck_filename");
         ParserPtr parser(new Opm::Parser());
-        Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
-        deck = parser->parseFile(deck_filename , parseMode);
-        eclipseState.reset(new Opm::EclipseState(deck , parseMode));
+        Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+        deck = parser->parseFile(deck_filename , parseContext);
+        eclipseState.reset(new Opm::EclipseState(deck , parseContext));
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // 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);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init
+            props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, ug_grid));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+
+            state.reset( new PolymerBlackoilState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+            // 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(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
+            initBlackoilSurfvol(ug_grid, *props, *state);
+            // Init polymer properties.
+            poly_props.readFromDeck(deck, eclipseState);
         }
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
-        // Init polymer properties.
-        poly_props.readFromDeck(deck, eclipseState);
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
@@ -132,31 +138,43 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        initBlackoilSurfvol(*grid->c_grid(), *props, state);
-        // Init Polymer state
-        if (param.has("poly_init")) {
-            double poly_init = param.getDefault("poly_init", 0.0);
-            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-                double smin[2], smax[2];
-                props->satRange(1, &cell, smin, smax);
-                if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                    state.concentration()[cell] = poly_init;
-                    state.maxconcentration()[cell] = poly_init;
-                } else {
-                    state.saturation()[2*cell + 0] = 0.;
-                    state.saturation()[2*cell + 1] = 1.;
-                    state.concentration()[cell] = 0.;
-                    state.maxconcentration()[cell] = 0.;
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new BlackoilPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+            state.reset( new PolymerBlackoilState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            initBlackoilSurfvol(ug_grid, *props, *state);
+            // Init Polymer state
+
+            if (param.has("poly_init")) {
+                double poly_init = param.getDefault("poly_init", 0.0);
+                for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
+                    double smin[2], smax[2];
+
+                    auto& saturation = state->saturation();
+                    auto& concentration = state->getCellData( state->CONCENTRATION );
+                    auto& max_concentration = state->getCellData( state->CMAX );
+
+                    props->satRange(1, &cell, smin, smax);
+                    if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                        concentration[cell] = poly_init;
+                        max_concentration[cell] = poly_init;
+                    } else {
+                        saturation[2*cell + 0] = 0.;
+                        saturation[2*cell + 1] = 1.;
+                        concentration[cell] = 0.;
+                        max_concentration[cell] = 0.;
+                    }
                 }
             }
+
         }
         // Init polymer properties.
         // Setting defaults to provide a simple example case.
@@ -240,8 +258,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more epochs etc.
         WellState well_state;
@@ -284,7 +302,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             // Create and run simulator.
@@ -300,7 +318,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
 
             // Update total timing report and remember step number.
             rep += epoch_rep;

examples/sim_poly2p_incomp_reorder.cpp

@@ -48,7 +48,7 @@
 #include <opm/polymer/PolymerProperties.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
 #include <boost/scoped_ptr.hpp>
@@ -92,36 +92,41 @@ try
     boost::scoped_ptr<IncompPropertiesInterface> props;
     boost::scoped_ptr<RockCompressibility> rock_comp;
     EclipseStateConstPtr eclipseState;
-    PolymerState state;
+    std::unique_ptr<PolymerState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
     if (use_deck) {
         std::string deck_filename = param.get<std::string>("deck_filename");
-        Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+        Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
         ParserPtr parser(new Opm::Parser());
-        deck = parser->parseFile(deck_filename , parseMode);
-        eclipseState.reset(new Opm::EclipseState(deck , parseMode));
+        deck = parser->parseFile(deck_filename , parseContext);
+        eclipseState.reset(new Opm::EclipseState(deck , parseContext));
 
         // Grid init
         grid.reset(new GridManager(deck));
-        // Rock and fluid init
-        props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid()));
-        // 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);
-        } else {
-            initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+            // Rock and fluid init
+            props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid ));
+            // check_well_controls = param.getDefault("check_well_controls", false);
+            // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+            state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+
+            // 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(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
+            // Init polymer properties.
+            poly_props.readFromDeck(deck, eclipseState);
         }
-        // Init polymer properties.
-        poly_props.readFromDeck(deck, eclipseState);
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
@@ -131,28 +136,38 @@ try
         const double dy = param.getDefault("dy", 1.0);
         const double dz = param.getDefault("dz", 1.0);
         grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
-        // Rock and fluid init.
-        props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-        // Rock compressibility.
-        rock_comp.reset(new RockCompressibility(param));
-        // Gravity.
-        gravity[2] = param.getDefault("gravity", 0.0);
-        // Init state variables (saturation and pressure).
-        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
-        // Init Polymer state
-        if (param.has("poly_init")) {
-            double poly_init = param.getDefault("poly_init", 0.0);
-            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-                double smin[2], smax[2];
-                props->satRange(1, &cell, smin, smax);
-                if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                    state.concentration()[cell] = poly_init;
-                    state.maxconcentration()[cell] = poly_init;
-                } else {
-                    state.saturation()[2*cell + 0] = 0.;
-                    state.saturation()[2*cell + 1] = 1.;
-                    state.concentration()[cell] = 0.;
-                    state.maxconcentration()[cell] = 0.;
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+
+            // Rock and fluid init.
+            props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));;
+            state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2));
+            // Rock compressibility.
+            rock_comp.reset(new RockCompressibility(param));
+            // Gravity.
+            gravity[2] = param.getDefault("gravity", 0.0);
+            // Init state variables (saturation and pressure).
+            initStateBasic(ug_grid, *props, param, gravity[2], *state);
+            // Init Polymer state
+            if (param.has("poly_init")) {
+                double poly_init = param.getDefault("poly_init", 0.0);
+                for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) {
+                    double smin[2], smax[2];
+
+                    auto& saturation = state->saturation();
+                    auto& concentration = state->getCellData( state->CONCENTRATION );
+                    auto& max_concentration = state->getCellData( state->CMAX );
+
+                    props->satRange(1, &cell, smin, smax);
+                    if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                        concentration[cell] = poly_init;
+                        max_concentration[cell] = poly_init;
+                    } else {
+                        saturation[2*cell + 0] = 0.;
+                        saturation[2*cell + 1] = 1.;
+                        concentration[cell] = 0.;
+                        max_concentration[cell] = 0.;
+                    }
                 }
             }
         }
@@ -212,7 +227,7 @@ try
         // terms of total pore volume.
         std::vector<double> porevol;
         if (rock_comp->isActive()) {
-            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
         } else {
             computePorevolume(*grid->c_grid(), props->porosity(), porevol);
         }
@@ -276,8 +291,8 @@ try
         simtimer.init(param);
         warnIfUnusedParams(param);
         WellState well_state;
-        well_state.init(0, state);
-        rep = simulator.run(simtimer, state, well_state);
+        well_state.init(0, *state);
+        rep = simulator.run(simtimer, *state, well_state);
     } else {
         // With a deck, we may have more epochs etc.
 
@@ -321,7 +336,7 @@ try
             // properly transfer old well state to it every report step,
             // since number of wells may change etc.
             if (reportStepIdx == 0) {
-                well_state.init(wells.c_wells(), state);
+                well_state.init(wells.c_wells(), *state);
             }
 
             // Create and run simulator.
@@ -339,7 +354,7 @@ try
             if (reportStepIdx == 0) {
                 warnIfUnusedParams(param);
             }
-            SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state);
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
 
             // Update total timing report and remember step number.
             rep += epoch_rep;

examples/sim_poly_fi2p_comp_ad.cpp

@@ -37,7 +37,7 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
-#include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/output/eclipse/EclipseWriter.hpp>
 #include <opm/core/props/BlackoilPropertiesBasic.hpp>
 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
@@ -61,12 +61,12 @@
 #include <opm/autodiff/GeoProps.hpp>
 #include <opm/autodiff/GridHelpers.hpp>
 
-#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
-#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
-#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/OpmLog/StreamLog.hpp>
+#include <opm/common/OpmLog/CounterLog.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
@@ -115,7 +115,7 @@ try
     std::shared_ptr<BlackoilPropertiesInterface> props;
     std::shared_ptr<BlackoilPropsAdInterface> new_props;
     std::shared_ptr<RockCompressibility> rock_comp;
-    PolymerBlackoilState state;
+    std::unique_ptr<PolymerBlackoilState> state;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
     double gravity[3] = { 0.0 };
@@ -141,7 +141,7 @@ try
     }
 
     std::string logFile = output_dir + "/LOGFILE.txt";
-    Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+    Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
     Opm::ParserPtr parser(new Opm::Parser());
     {
         std::shared_ptr<Opm::StreamLog> streamLog = std::make_shared<Opm::StreamLog>(logFile , Opm::Log::DefaultMessageTypes);
@@ -154,9 +154,9 @@ try
     Opm::DeckConstPtr deck;
     std::shared_ptr<EclipseState> eclipseState;
     try {
-        deck = parser->parseFile(deck_filename , parseMode);
+        deck = parser->parseFile(deck_filename , parseContext);
         Opm::checkDeck(deck, parser);
-        eclipseState.reset(new Opm::EclipseState(deck , parseMode));
+        eclipseState.reset(new Opm::EclipseState(deck , parseContext));
     }
     catch (const std::invalid_argument& e) {
         std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl;
@@ -187,6 +187,8 @@ try
                                                 Opm::UgGridHelpers::globalCell(cGrid),
                                                 Opm::UgGridHelpers::cartDims(cGrid),
                                                 param));
+
+    state.reset( new PolymerBlackoilState( Opm::UgGridHelpers::numCells(cGrid), Opm::UgGridHelpers::numFaces(cGrid), 2));
     new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, materialLawManager, cGrid));
     PolymerProperties polymer_props(deck, eclipseState);
     PolymerPropsAd polymer_props_ad(polymer_props);
@@ -199,10 +201,10 @@ try
 
     // 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);
+        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);
+        initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], *state);
     }
 
     bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
@@ -254,7 +256,7 @@ try
                   deck,
                   *fis_solver,
                   grav);
-    fullReport= simulator.run(simtimer, state);
+    fullReport= simulator.run(simtimer, *state);
 
     std::cout << "\n\n================    End of simulation     ===============\n\n";
     fullReport.report(std::cout);

opm/autodiff/AutoDiffBlock.hpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -623,25 +624,20 @@ namespace Opm
         return rhs * lhs; // Commutative operation.
     }
 
+
     /**
-     * @brief Computes the value of base raised to the power of exp elementwise
+     * @brief Computes the value of base raised to the power of exponent
      *
      * @param base The AD forward block
-     * @param exp  array of exponents
-     * @return The value of base raised to the power of exp elementwise
+     * @param exponent  double
+     * @return The value of base raised to the power of exponent
      */
     template <typename Scalar>
     AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
-                              const typename AutoDiffBlock<Scalar>::V& exp)
+                              const double exponent)
     {
-        const int num_elem = base.value().size();
-        typename AutoDiffBlock<Scalar>::V val (num_elem);
-        typename AutoDiffBlock<Scalar>::V derivative = exp;
-        assert(exp.size() == num_elem);
-        for (int i = 0; i < num_elem; ++i) {
-            val[i] = std::pow(base.value()[i], exp[i]);
-            derivative[i] *= std::pow(base.value()[i], exp[i] - 1.0);
-        }
+        const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exponent);
+        const typename AutoDiffBlock<Scalar>::V derivative = exponent * base.value().pow(exponent - 1.0);
         const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
 
         std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
@@ -652,28 +648,95 @@ namespace Opm
         return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
     }
 
-
     /**
-     * @brief Computes the value of base raised to the power of exp
+     * @brief Computes the value of base raised to the power of exponent
      *
      * @param base The AD forward block
-     * @param exp  exponent
+     * @param exponent  Array of exponents
+     * @return The value of base raised to the power of exponent elementwise
+     */
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
+                              const typename AutoDiffBlock<Scalar>::V& exponent)
+    {
+        // Add trivial derivatives and use the AD pow function
+        return  pow(base,AutoDiffBlock<Scalar>::constant(exponent));
+    }
+
+    /**
+     * @brief Computes the value of base raised to the power of exponent
+     *
+     * @param base      Array of base values
+     * @param exponent  The AD forward block
+     * @return The value of base raised to the power of exponent elementwise
+     */
+    template <typename Scalar>
+    AutoDiffBlock<Scalar> pow(const typename AutoDiffBlock<Scalar>::V& base,
+                              const AutoDiffBlock<Scalar>& exponent)
+    {
+        // Add trivial derivatives and use the AD pow function
+        return pow(AutoDiffBlock<Scalar>::constant(base),exponent);
+    }
+
+    /**
+     * @brief Computes the value of base raised to the power of exponent
+     *
+     * @param base The base AD forward block
+     * @param exponent  The exponent AD forward block
      * @return The value of base raised to the power of exp
      */
     template <typename Scalar>
     AutoDiffBlock<Scalar> pow(const AutoDiffBlock<Scalar>& base,
-                              const double exp)
-    {
-        const typename AutoDiffBlock<Scalar>::V val = base.value().pow(exp);
-        const typename AutoDiffBlock<Scalar>::V derivative = exp * base.value().pow(exp - 1.0);
-        const typename AutoDiffBlock<Scalar>::M derivative_diag(derivative.matrix().asDiagonal());
-
-        std::vector< typename AutoDiffBlock<Scalar>::M > jac (base.numBlocks());
-        for (int block = 0; block < base.numBlocks(); block++) {
-             fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]);
+                              const AutoDiffBlock<Scalar>& exponent)
+    {                       
+        const int num_elem = base.value().size();
+        assert(exponent.size() == num_elem);
+        typename AutoDiffBlock<Scalar>::V val (num_elem);
+        for (int i = 0; i < num_elem; ++i) {
+            val[i] = std::pow(base.value()[i], exponent.value()[i]);
         }
 
-        return AutoDiffBlock<Scalar>::function( std::move(val), std::move(jac) );
+        // (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' = der1 + der2
+        // if f' is empty only der1 is calculated
+        // if g' is empty only der2 is calculated
+        // if f' and g' are non empty they should have the same size
+        int num_blocks = std::max (base.numBlocks(), exponent.numBlocks());
+        if (!base.derivative().empty() && !exponent.derivative().empty()) {
+            assert(exponent.numBlocks() == base.numBlocks());
+        }
+        std::vector< typename AutoDiffBlock<Scalar>::M > jac (num_blocks);
+
+        if ( !exponent.derivative().empty() ) {
+            typename AutoDiffBlock<Scalar>::V der1 = val;
+            for (int i = 0; i < num_elem; ++i) {
+                der1[i] *= std::log(base.value()[i]);
+            }
+            std::vector< typename AutoDiffBlock<Scalar>::M > jac1 (exponent.numBlocks());
+            const typename AutoDiffBlock<Scalar>::M der1_diag(der1.matrix().asDiagonal());
+            for (int block = 0; block < exponent.numBlocks(); block++) {
+                fastSparseProduct(der1_diag, exponent.derivative()[block], jac1[block]);
+                jac[block] = jac1[block];
+            }
+        }
+
+        if ( !base.derivative().empty() ) {
+            typename AutoDiffBlock<Scalar>::V der2 = exponent.value();
+            for (int i = 0; i < num_elem; ++i) {
+                der2[i] *= std::pow(base.value()[i], exponent.value()[i] - 1.0);
+            }
+            std::vector< typename AutoDiffBlock<Scalar>::M > jac2 (base.numBlocks());
+            const typename AutoDiffBlock<Scalar>::M der2_diag(der2.matrix().asDiagonal());
+            for (int block = 0; block < base.numBlocks(); block++) {
+                fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]);
+                if (!exponent.derivative().empty()) {
+                    jac[block] += jac2[block];
+                } else {
+                    jac[block] = jac2[block];
+                }
+            }
+        }
+
+        return AutoDiffBlock<Scalar>::function(std::move(val), std::move(jac));
     }
 
 

opm/autodiff/BackupRestore.hpp

@@ -21,8 +21,9 @@
 #define OPM_BACKUPRESTORE_HEADER_INCLUDED
 
 #include <iostream>
+#include <opm/common/data/SimulationDataContainer.hpp>
+
 
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
 
@@ -159,7 +160,7 @@ namespace Opm {
                BlackoilStateId = 2,
                WellStateFullyImplicitBackoilId = 3 };
 
-        inline int objectId( const SimulatorState& /* state */) {
+        inline int objectId( const SimulationDataContainer& /* state */) {
             return SimulatorStateId;
         }
         inline int objectId( const WellState& /* state */) {
@@ -184,9 +185,9 @@ namespace Opm {
         }
     }
 
-    // SimulatorState
+    // SimulationDataContainer
     inline
-    std::ostream& operator << (std::ostream& out, const SimulatorState& state )
+    std::ostream& operator << (std::ostream& out, const SimulationDataContainer& state )
     {
         // write id of object to stream
         writeValue( out, objectId( state ) );
@@ -205,7 +206,7 @@ namespace Opm {
     }
 
     inline
-    std::istream& operator >> (std::istream& in, SimulatorState& state )
+    std::istream& operator >> (std::istream& in, SimulationDataContainer& state )
     {
         // check id of stored object
         checkObjectId( in, state );
@@ -213,7 +214,7 @@ namespace Opm {
         int numPhases = 0;
         readValue( in, numPhases );
 
-        if( numPhases != state.numPhases() )
+        if( numPhases != (int) state.numPhases() )
             OPM_THROW(std::logic_error,"num phases wrong");
 
         // read variables
@@ -234,7 +235,7 @@ namespace Opm {
         writeValue( out, objectId( state ) );
 
         // backup simulator state
-        const SimulatorState& simstate = static_cast< const SimulatorState& > (state);
+        const SimulationDataContainer& simstate = static_cast< const SimulationDataContainer& > (state);
         out << simstate;
 
         // backup additional blackoil state variables
@@ -252,7 +253,7 @@ namespace Opm {
         checkObjectId( in, state );
 
         // restore simulator state
-        SimulatorState& simstate = static_cast< SimulatorState& > (state);
+        SimulationDataContainer& simstate = static_cast< SimulationDataContainer& > (state);
         in >> simstate;
 
         // restore additional blackoil state variables

opm/autodiff/BlackoilModelBase.hpp

@@ -45,7 +45,7 @@ namespace Opm {
     class RockCompressibility;
     class NewtonIterationBlackoilInterface;
     class VFPProperties;
-
+    class SimulationDataContainer;
 
     /// Struct for containing iteration variables.
     struct DefaultBlackoilSolutionState
@@ -207,7 +207,7 @@ namespace Opm {
 
         /// \brief compute the relative change between to simulation states
         //  \return || u^n+1 - u^n || / || u^n+1 ||
-        double relativeChange( const SimulatorState& previous, const SimulatorState& current ) const;
+        double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const;
 
         /// The size (number of unknowns) of the nonlinear system of equations.
         int sizeNonLinear() const;

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -48,6 +48,7 @@
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
 
+#include <opm/common/data/SimulationDataContainer.hpp>
 #include <cassert>
 #include <cmath>
 #include <iostream>
@@ -1610,11 +1611,37 @@ namespace detail {
                 break;
 
             case SURFACE_RATE:
-                for (int phase = 0; phase < np; ++phase) {
-                    if (distr[phase] > 0.0) {
-                        xw.wellRates()[np*w + phase] = target * distr[phase];
+                // assign target value as initial guess for injectors and
+                // single phase producers (orat, grat, wrat)
+                const WellType& well_type = wells().type[w];
+                if (well_type == INJECTOR) {
+                    for (int phase = 0; phase < np; ++phase) {
+                        const double& compi = wells().comp_frac[np * w + phase];
+                        if (compi > 0.0) {
+                            xw.wellRates()[np*w + phase] = target * compi;
+                        }
                     }
+                } else if (well_type == PRODUCER) {
+
+                    // only set target as initial rates for single phase
+                    // producers. (orat, grat and wrat, and not lrat)
+                    // lrat will result in numPhasesWithTargetsUnderThisControl == 2
+                    int numPhasesWithTargetsUnderThisControl = 0;
+                    for (int phase = 0; phase < np; ++phase) {
+                        if (distr[phase] > 0.0) {
+                            numPhasesWithTargetsUnderThisControl += 1;
+                        }
+                    }
+                    for (int phase = 0; phase < np; ++phase) {
+                        if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
+                            xw.wellRates()[np*w + phase] = target * distr[phase];
+                        }
+                    }
+                } else {
+                    OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
                 }
+
+
                 break;
             }
 
@@ -2520,8 +2547,8 @@ namespace detail {
     template <class Grid, class Implementation>
     double
     BlackoilModelBase<Grid, Implementation>::
-    relativeChange(const SimulatorState& previous,
-                   const SimulatorState& current ) const
+    relativeChange(const SimulationDataContainer& previous,
+                   const SimulationDataContainer& current ) const
     {
         std::vector< double > p0  ( previous.pressure() );
         std::vector< double > sat0( previous.saturation() );

opm/autodiff/BlackoilMultiSegmentModel.hpp

@@ -205,6 +205,7 @@ namespace Opm {
             AutoDiffMatrix eliminate_topseg;              // change the top segment related to be zero
             std::vector<int> well_cells;                  // the set of perforated cells
             V conn_trans_factors;                         // connection transmissibility factors
+            bool has_multisegment_wells;                  // flag indicating whether there is any muli-segment well
         };
 
         MultiSegmentWellOps wops_ms_;

opm/autodiff/BlackoilMultiSegmentModel_impl.hpp

@@ -109,6 +109,9 @@ namespace Opm {
     BlackoilMultiSegmentModel<Grid>::
     MultiSegmentWellOps::MultiSegmentWellOps(const std::vector<WellMultiSegmentConstPtr>& wells_ms)
     {
+        // no multi-segment wells are involved by default.
+        has_multisegment_wells = false;
+
         if (wells_ms.empty()) {
             return;
         }
@@ -120,6 +123,9 @@ namespace Opm {
         for (int w = 0; w < nw; ++w) {
             total_nperf += wells_ms[w]->numberOfPerforations();
             total_nseg += wells_ms[w]->numberOfSegments();
+            if (wells_ms[w]->isMultiSegmented()) {
+                has_multisegment_wells = true;
+            }
         }
 
         // Create well_cells and conn_trans_factors.
@@ -230,6 +236,12 @@ namespace Opm {
         top_well_segments_ = well_state.topSegmentLoc();
 
         const int nw = wellsMultiSegment().size();
+
+        if ( !wops_ms_.has_multisegment_wells ) {
+            segvdt_ = V::Zero(nw);
+            return;
+        }
+
         const int nseg_total = well_state.numSegments();
         std::vector<double> segment_volume;
         segment_volume.reserve(nseg_total);
@@ -462,6 +474,12 @@ namespace Opm {
         well_perforation_densities_ = Eigen::Map<const V>(cd.data(), nperf_total); // This one is not useful for segmented wells at all
         well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf_total);
 
+        if ( !wops_ms_.has_multisegment_wells ) {
+            well_perforation_cell_densities_ = V::Zero(nperf_total);
+            well_perforation_cell_pressure_diffs_ = V::Zero(nperf_total);
+            return;
+        }
+
         // compute the average of the fluid densites in the well blocks.
         // the average is weighted according to the fluid relative permeabilities.
         const std::vector<ADB> kr_adb = Base::computeRelPerm(state);
@@ -680,7 +698,7 @@ namespace Opm {
 
             // Compute drawdown.
             ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_,
-                                                    ADB::constant(well_perforation_pressure_diffs_));
+                                              ADB::constant(well_perforation_pressure_diffs_));
             const V h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, V::Zero(nperf));
 
             // Special handling for when we are called from solveWellEq().
@@ -890,23 +908,27 @@ namespace Opm {
 
         std::vector<ADB> segment_volume_change_dt(np, ADB::null());
         for (int phase = 0; phase < np; ++phase) {
-            // Gain of the surface volume of each component in the segment by dt
-            segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
-                                              segment_comp_surf_volume_initial_[phase];
+            if ( wops_ms_.has_multisegment_wells ) {
+                // Gain of the surface volume of each component in the segment by dt
+                segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] -
+                                                  segment_comp_surf_volume_initial_[phase];
 
-            // Special handling for when we are called from solveWellEq().
-            // TODO: restructure to eliminate need for special treatmemt.
-            if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
-                assert(segment_volume_change_dt[phase].numBlocks() > 2);
-                assert(segqs.numBlocks() == 2);
-                segment_volume_change_dt[phase] = detail::onlyWellDerivs(segment_volume_change_dt[phase]);
-                assert(segment_volume_change_dt[phase].numBlocks() == 2);
+                // Special handling for when we are called from solveWellEq().
+                // TODO: restructure to eliminate need for special treatmemt.
+                if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) {
+                    assert(segment_volume_change_dt[phase].numBlocks() > 2);
+                    assert(segqs.numBlocks() == 2);
+                    segment_volume_change_dt[phase] = detail::onlyWellDerivs(segment_volume_change_dt[phase]);
+                    assert(segment_volume_change_dt[phase].numBlocks() == 2);
+                }
+
+                const ADB cq_s_seg = wops_ms_.p2s * cq_s[phase];
+                const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total));
+                segqs -= superset(cq_s_seg + wops_ms_.s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
+                                  Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
+            } else {
+                segqs -= superset(wops_ms_.p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
             }
-
-            const ADB cq_s_seg = wops_ms_.p2s * cq_s[phase];
-            const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total));
-            segqs -= superset(cq_s_seg + wops_ms_.s2s_inlets * segqs_phase + segment_volume_change_dt[phase],
-                              Span(nseg_total, 1, phase * nseg_total), np * nseg_total);
         }
 
         residual_.well_flux_eq = segqs;
@@ -1169,13 +1191,17 @@ namespace Opm {
 
         ADB others_residual = ADB::constant(V::Zero(nseg_total));
 
-        // Special handling for when we are called from solveWellEq().
-        // TODO: restructure to eliminate need for special treatmemt.
-        ADB wspd = (state.segp.numBlocks() == 2)
-            ? detail::onlyWellDerivs(well_segment_pressures_delta_)
-            : well_segment_pressures_delta_;
+        if ( wops_ms_.has_multisegment_wells ) {
+            // Special handling for when we are called from solveWellEq().
+            // TODO: restructure to eliminate need for special treatmemt.
+            ADB wspd = (state.segp.numBlocks() == 2)
+                ? detail::onlyWellDerivs(well_segment_pressures_delta_)
+                : well_segment_pressures_delta_;
 
-        others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp + wspd);
+            others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp + wspd);
+        } else {
+            others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp);
+        }
 
         //       all the control equations
         // TODO: can be optimized better
@@ -1284,10 +1310,25 @@ namespace Opm {
     void
     BlackoilMultiSegmentModel<Grid>::computeSegmentFluidProperties(const SolutionState& state)
     {
+
         const int nw = wellsMultiSegment().size();
         const int nseg_total = state.segp.size();
         const int np = numPhases();
 
+        if ( !wops_ms_.has_multisegment_wells ){
+            // not sure if this is needed actually
+            // TODO: to check later if this is really necessary.
+            segment_mass_flow_rates_ = ADB::constant(V::Zero(nseg_total));
+            well_segment_densities_ = ADB::constant(V::Zero(nseg_total));
+            segment_mass_flow_rates_ = ADB::constant(V::Zero(nseg_total));
+            segment_viscosities_ = ADB::constant(V::Zero(nseg_total));
+            for (int phase = 0; phase < np; ++phase) {
+                segment_comp_surf_volume_current_[phase] = ADB::constant(V::Zero(nseg_total));
+                segment_comp_surf_volume_initial_[phase] = V::Zero(nseg_total);
+            }
+            return;
+        }
+
         // although we will calculate segment density for non-segmented wells at the same time,
         // while under most of the cases, they will not be used,
         // since for most of the cases, the density calculation for non-segment wells are
@@ -1475,6 +1516,12 @@ namespace Opm {
         const int nw = wellsMultiSegment().size();
         const int nseg_total = state.segp.size();
 
+        if ( !wops_ms_.has_multisegment_wells ) {
+            well_segment_pressures_delta_ = ADB::constant(V::Zero(nseg_total));
+            well_segment_perforation_pressure_diffs_ = wops_ms_.s2p * well_segment_pressures_delta_;
+            return;
+        }
+
         // calculate the depth difference of the segments
         // TODO: we need to store the following values somewhere to avoid recomputation.
         V segment_depth_delta = V::Zero(nseg_total);

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -22,27 +22,20 @@
 #include <config.h>
 
 #include <opm/autodiff/BlackoilPropsAdFromDeck.hpp>
+#include <opm/autodiff/AutoDiffHelpers.hpp>
+
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/extractPvtTableIndex.hpp>
 
 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
 
-#include <opm/autodiff/AutoDiffHelpers.hpp>
-#include <opm/core/props/BlackoilPropertiesInterface.hpp>
-#include <opm/core/props/BlackoilPhases.hpp>
-#include <opm/core/props/pvt/PvtInterface.hpp>
-#include <opm/core/props/pvt/PvtConstCompr.hpp>
-#include <opm/core/props/pvt/PvtDead.hpp>
-#include <opm/core/props/pvt/PvtDeadSpline.hpp>
-#include <opm/core/props/pvt/PvtLiveOil.hpp>
-#include <opm/core/props/pvt/PvtLiveGas.hpp>
-#include <opm/core/props/pvt/ThermalWaterPvtWrapper.hpp>
-#include <opm/core/props/pvt/ThermalOilPvtWrapper.hpp>
-#include <opm/core/props/pvt/ThermalGasPvtWrapper.hpp>
-#include <opm/common/ErrorMacros.hpp>
-#include <opm/core/utility/Units.hpp>
-
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
+#include <opm/common/ErrorMacros.hpp>
+
 namespace Opm
 {
     // Making these typedef to make the code more readable.

opm/autodiff/BlackoilSolventModel.hpp

@@ -144,7 +144,6 @@ namespace Opm {
         using Base::wellsActive;
         using Base::wells;
         using Base::variableState;
-        using Base::computePressures;
         using Base::computeGasPressure;
         using Base::applyThresholdPressures;
         using Base::fluidRsSat;
@@ -157,7 +156,6 @@ namespace Opm {
         using Base::dsMax;
         using Base::drMaxRel;
         using Base::maxResidualAllowed;
-
         using Base::updateWellControls;
         using Base::computeWellConnectionPressures;
         using Base::addWellControlEq;
@@ -237,6 +235,14 @@ namespace Opm {
         // compute density and viscosity using the ToddLongstaff mixing model
         void computeToddLongstaffMixing(std::vector<ADB>& viscosity, std::vector<ADB>& density, const std::vector<ADB>& saturations, const Opm::PhaseUsage pu);
 
+        // compute phase pressures.
+        std::vector<ADB>
+        computePressures(const ADB& po,
+                         const ADB& sw,
+                         const ADB& so,
+                         const ADB& sg,
+                         const ADB& ss) const;
+
     };
 
 

opm/autodiff/BlackoilSolventModel_impl.hpp

@@ -137,11 +137,14 @@ namespace Opm {
 
         // Initial polymer concentration.
         if (has_solvent_) {
-            assert (not x.solvent_saturation().empty());
-            const int nc = x.solvent_saturation().size();
-            const V ss = Eigen::Map<const V>(&x.solvent_saturation()[0], nc);
+            const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL );
+            const int nc = solvent_saturation.size();
+            const V ss = Eigen::Map<const V>(solvent_saturation.data() , nc);
+
+            // This is some insanely detailed flickety flackety code;
             // Solvent belongs after other reservoir vars but before well vars.
             auto solvent_pos = vars0.begin() + fluid_.numPhases();
+            assert (not solvent_saturation.empty());
             assert(solvent_pos == vars0.end() - 2);
             vars0.insert(solvent_pos, ss);
         }
@@ -177,17 +180,72 @@ namespace Opm {
     BlackoilSolventModel<Grid>::variableStateExtractVars(const ReservoirState& x,
                                                          const std::vector<int>& indices,
                                                          std::vector<ADB>& vars) const
-    {
-        SolutionState state = Base::variableStateExtractVars(x, indices, vars);
-        if (has_solvent_) {
-            state.solvent_saturation = std::move(vars[indices[Solvent]]);
+    {         
+        // This is more or less a copy of the base class. Refactoring is needed in the base class
+        // to avoid unnecessary copying.
+
+        //using namespace Opm::AutoDiffGrid;
+        const int nc = Opm::AutoDiffGrid::numCells(grid_);
+        const Opm::PhaseUsage pu = fluid_.phaseUsage();
+
+        SolutionState state(fluid_.numPhases());
+
+        // Pressure.
+        state.pressure = std::move(vars[indices[Pressure]]);
+
+        // Temperature cannot be a variable at this time (only constant).
+        const V temp = Eigen::Map<const V>(& x.temperature()[0], x.temperature().size());
+        state.temperature = ADB::constant(temp);
+
+        // Saturations
+        {
+            ADB so = ADB::constant(V::Ones(nc, 1));
+
+            if (active_[ Water ]) {
+                state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]);
+                const ADB& sw = state.saturation[pu.phase_pos[ Water ]];
+                so -= sw;
+            }
+            if (has_solvent_) {
+                state.solvent_saturation = std::move(vars[indices[Solvent]]);
+                so -= state.solvent_saturation;
+            }
+
+            if (active_[ Gas ]) {
+                // Define Sg Rs and Rv in terms of xvar.
+                // Xvar is only defined if gas phase is active
+                const ADB& xvar = vars[indices[Xvar]];
+                ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ];
+                sg = Base::isSg_*xvar + Base::isRv_*so;
+                so -= sg;
+
+                if (active_[ Oil ]) {
+                    // RS and RV is only defined if both oil and gas phase are active.
+                    state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation);
+                    const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_);
+                    if (has_disgas_) {
+                        state.rs = (1-Base::isRs_)*rsSat + Base::isRs_*xvar;
+                    } else {
+                        state.rs = rsSat;
+                    }
+                    const ADB rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_);
+                    if (has_vapoil_) {
+                        state.rv = (1-Base::isRv_)*rvSat + Base::isRv_*xvar;
+                    } else {
+                        state.rv = rvSat;
+                    }
+                }
+            }
+
             if (active_[ Oil ]) {
                 // Note that so is never a primary variable.
-                const Opm::PhaseUsage pu = fluid_.phaseUsage();
-                state.saturation[pu.phase_pos[ Oil ]] -= state.solvent_saturation;
+                state.saturation[pu.phase_pos[ Oil ]] = std::move(so);
             }
         }
+        // wells
+        Base::variableStateExtractWellsVars(indices, vars, state);
         return state;
+
     }
 
 
@@ -493,9 +551,10 @@ namespace Opm {
             Base::updateState(modified_dx, reservoir_state, well_state);
 
             // Update solvent.
-            const V ss_old = Eigen::Map<const V>(&reservoir_state.solvent_saturation()[0], nc, 1);
+            auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL );
+            const V ss_old = Eigen::Map<const V>(&solvent_saturation[0], nc, 1);
             const V ss = (ss_old - dss).max(zero);
-            std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin());
+            std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin());
 
             // adjust oil saturation
             const Opm::PhaseUsage& pu = fluid_.phaseUsage();
@@ -665,7 +724,9 @@ namespace Opm {
 
                 Selector<double> zero_selector(ss.value() + sg.value(), Selector<double>::Zero);
                 ADB F_solvent = zero_selector.select(ss, ss / (ss + sg));
-                const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_);
+                const ADB& po = state.canonical_phase_pressures[ Oil ];
+                const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_)
+                        * solvent_props_.pressureMiscibilityFunction(po, cells_);
 
                 assert(active_[ Oil ]);
                 assert(active_[ Gas ]);
@@ -769,16 +830,27 @@ namespace Opm {
         // Compute effective viscosities and densities
         computeToddLongstaffMixing(viscosity, density, effective_saturations, pu);
 
-        // Store the computed volume factors and viscosities
-        b_eff_[pu.phase_pos[ Water ]] = bw;
-        b_eff_[pu.phase_pos[ Oil ]] = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ],  cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs);
-        b_eff_[pu.phase_pos[ Gas ]] = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ],  cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv);
-        b_eff_[solvent_pos_] = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_);
+        // compute the volume factors from the densities
+        const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ],  cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs);
+        const ADB b_eff_g = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ],  cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv);
+        const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_);
 
+        // account for pressure effects and store the computed volume factors and viscosities
+        const V ones = V::Constant(nc, 1.0);
+        const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
+
+        b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo;
+        b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg;
+        b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs;
+
+        // keep the mu*b interpolation
+        mu_eff_[pu.phase_pos[ Oil ]] = b_eff_[pu.phase_pos[ Oil ]] / (pmisc * b_eff_o / viscosity[pu.phase_pos[ Oil ]] + (ones - pmisc) * bo / mu_o);
+        mu_eff_[pu.phase_pos[ Gas ]] = b_eff_[pu.phase_pos[ Gas ]] / (pmisc * b_eff_g / viscosity[pu.phase_pos[ Gas ]] + (ones - pmisc) * bg / mu_g);
+        mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s);
+
+        // for water the pure values are used
         mu_eff_[pu.phase_pos[ Water ]] = mu_w;
-        mu_eff_[pu.phase_pos[ Oil ]] = viscosity[pu.phase_pos[ Oil ]];
-        mu_eff_[pu.phase_pos[ Gas ]] = viscosity[pu.phase_pos[ Gas ]];
-        mu_eff_[solvent_pos_] = viscosity[solvent_pos_];
+        b_eff_[pu.phase_pos[ Water ]] = bw;
     }
 
     template <class Grid>
@@ -881,6 +953,30 @@ namespace Opm {
     }
 
 
+    template <class Grid>
+    std::vector<ADB>
+    BlackoilSolventModel<Grid>::
+    computePressures(const ADB& po,
+                     const ADB& sw,
+                     const ADB& so,
+                     const ADB& sg,
+                     const ADB& ss) const
+    {
+        std::vector<ADB> pressures = Base::computePressures(po, sw, so, sg);
+        // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss)
+        std::vector<ADB> pressures_imisc = Base::computePressures(po, sw, so, sg + ss);
+
+        // Pressure effects on capillary pressure miscibility
+        const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_);
+        // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent
+        // to changing the gas pressure directly.
+        const int  nc = cells_.size();
+        const V ones = V::Constant(nc, 1.0);
+        pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas]));
+        return pressures;
+    }
+
+
     template <class Grid>
     void
     BlackoilSolventModel<Grid>::assemble(const ReservoirState& reservoir_state,

opm/autodiff/BlackoilSolventState.cpp

@@ -0,0 +1,34 @@
+/*
+  Copyright 2015 IRIS AS, Applied Mathematics.
+
+  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/autodiff/BlackoilSolventState.hpp>
+
+namespace Opm
+{
+    const std::string BlackoilSolventState::SSOL = "SSOL";
+
+    BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases)
+        : BlackoilState( number_of_cells , number_of_faces , number_of_phases)
+    {
+        registerCellData( SSOL , 1 );
+    };
+
+} // namespace Opm
+

opm/autodiff/BlackoilSolventState.hpp

@@ -20,9 +20,9 @@
 #ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
 #define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED
 
-#include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/grid.h>
-#include <vector>
+#include <string>
+
+ #include <opm/core/simulator/BlackoilState.hpp>
 
 
 namespace Opm
@@ -33,22 +33,9 @@ namespace Opm
     class BlackoilSolventState : public BlackoilState
     {
     public:
-        void init(const UnstructuredGrid& g, int num_phases)
-        {
-            this->init(g.number_of_cells, g.number_of_faces, num_phases);
-        }
+        static const std::string SSOL;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            BlackoilState::init(number_of_cells, number_of_faces, num_phases);
-            solventId_  = SimulatorState::registerCellData( "SSOL", 1 );
-        }
-
-        std::vector<double>& solvent_saturation()    { return cellData()[ solventId_ ]; }
-        const std::vector<double>& solvent_saturation() const    { return cellData()[ solventId_ ]; }
-
-    private:
-        int solventId_;
+        BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases);
     };
 
 } // namespace Opm

opm/autodiff/FlowMain.hpp

@@ -70,12 +70,11 @@
 
 #include <opm/core/utility/share_obj.hpp>
 
-#include <opm/parser/eclipse/OpmLog/OpmLog.hpp>
-#include <opm/parser/eclipse/OpmLog/StreamLog.hpp>
-#include <opm/parser/eclipse/OpmLog/CounterLog.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/OpmLog/EclipsePRTLog.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
@@ -95,6 +94,7 @@
 #include <vector>
 #include <numeric>
 #include <cstdlib>
+#include <stdexcept>
 
 
 
@@ -102,6 +102,33 @@
 namespace Opm
 {
 
+    boost::filesystem::path simulationCaseName( const std::string& casename ) {
+        namespace fs = boost::filesystem;
+
+        const auto exists = []( const fs::path& f ) -> bool {
+            if( !fs::exists( f ) ) return false;
+
+            if( fs::is_regular_file( f ) ) return true;
+
+            return fs::is_symlink( f )
+                && fs::is_regular_file( fs::read_symlink( f ) );
+        };
+
+        auto simcase = fs::path( casename );
+
+        if( exists( simcase ) ) {
+            return simcase;
+        }
+
+        for( const auto& ext : { std::string("data"), std::string("DATA") } ) {
+            if( exists( simcase.replace_extension( ext ) ) ) {
+                return simcase;
+            }
+        }
+
+        throw std::invalid_argument( "Cannot find input case " + casename );
+    }
+
     /// This class encapsulates the setup and running of
     /// a simulator based on an input deck.
     template <class Implementation, class Grid, class Simulator>
@@ -184,7 +211,8 @@ namespace Opm
         bool use_local_perm_ = true;
         std::unique_ptr<DerivedGeology> geoprops_;
         // setupState()
-        ReservoirState state_;
+        std::unique_ptr<ReservoirState> state_;
+
         std::vector<double> threshold_pressures_;
         // distributeData()
         boost::any parallel_information_;
@@ -259,10 +287,6 @@ namespace Opm
             }
         }
 
-
-
-
-
         // Read parameters, see if a deck was specified on the command line, and if
         // it was, insert it into parameters.
         // Writes to:
@@ -282,7 +306,8 @@ namespace Opm
                     std::cerr << "You can only specify a single input deck on the command line.\n";
                     return false;
                 } else {
-                    param_.insertParameter("deck_filename", param_.unhandledArguments()[0]);
+                    const auto casename = simulationCaseName( param_.unhandledArguments()[ 0 ] );
+                    param_.insertParameter("deck_filename", casename.string() );
                 }
             }
 
@@ -352,19 +377,16 @@ namespace Opm
             // Create Parser
             ParserPtr parser(new Parser());
             {
-                std::shared_ptr<StreamLog> streamLog = std::make_shared<StreamLog>(logFile_ , Log::DefaultMessageTypes);
-                std::shared_ptr<CounterLog> counterLog = std::make_shared<CounterLog>(Log::DefaultMessageTypes);
-
-                OpmLog::addBackend( "STREAM" , streamLog );
-                OpmLog::addBackend( "COUNTER" , counterLog );
+                std::shared_ptr<EclipsePRTLog> prtLog = std::make_shared<EclipsePRTLog>(logFile_ , Log::DefaultMessageTypes);
+                OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog );
             }
 
             // Create Deck and EclipseState.
             try {
-                ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
-                deck_ = parser->parseFile(deck_filename, parseMode);
+                ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }});
+                deck_ = parser->parseFile(deck_filename, parseContext);
                 checkDeck(deck_, parser);
-                eclipse_state_.reset(new EclipseState(deck_, parseMode));
+                eclipse_state_.reset(new EclipseState(deck_, parseContext));
             }
             catch (const std::invalid_argument& e) {
                 std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile_ << std::endl;
@@ -445,8 +467,13 @@ namespace Opm
                                               Opm::UgGridHelpers::cartDims(grid),
                                               param_);
 
+
             // Init state variables (saturation and pressure).
             if (param_.has("init_saturation")) {
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases() ));
+
                 initStateBasic(Opm::UgGridHelpers::numCells(grid),
                                Opm::UgGridHelpers::globalCell(grid),
                                Opm::UgGridHelpers::cartDims(grid),
@@ -455,26 +482,35 @@ namespace Opm
                                Opm::UgGridHelpers::beginFaceCentroids(grid),
                                Opm::UgGridHelpers::beginCellCentroids(grid),
                                Opm::UgGridHelpers::dimensions(grid),
-                               props, param_, gravity_[2], state_);
+                               props, param_, gravity_[2], *state_);
 
-                initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state_);
+                initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_);
 
                 enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour };
                 if (pu.phase_used[Oil] && pu.phase_used[Gas]) {
                     const int numPhases = props.numPhases();
                     const int numCells  = Opm::UgGridHelpers::numCells(grid);
+
+                    // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState.
+                    auto& gor = state_->getCellData( BlackoilState::GASOILRATIO );
+                    const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL );
                     for (int c = 0; c < numCells; ++c) {
-                        state_.gasoilratio()[c] = state_.surfacevol()[c*numPhases + pu.phase_pos[Gas]]
-                            / state_.surfacevol()[c*numPhases + pu.phase_pos[Oil]];
+                        // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field.
+                        gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]];
                     }
                 }
             } else if (deck_->hasKeyword("EQUIL") && props.numPhases() == 3) {
-                state_.init(Opm::UgGridHelpers::numCells(grid),
-                           Opm::UgGridHelpers::numFaces(grid),
-                           props.numPhases());
-                initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], state_);
-                state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
+                // Which state class are we really using - what a f... mess?
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases()));
+
+                initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], *state_);
+                //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0);
             } else {
+                state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid),
+                                                  Opm::UgGridHelpers::numFaces(grid),
+                                                  props.numPhases()));
                 initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid),
                                           Opm::UgGridHelpers::globalCell(grid),
                                           Opm::UgGridHelpers::numFaces(grid),
@@ -482,12 +518,12 @@ namespace Opm
                                           Opm::UgGridHelpers::beginFaceCentroids(grid),
                                           Opm::UgGridHelpers::beginCellCentroids(grid),
                                           Opm::UgGridHelpers::dimensions(grid),
-                                          props, deck_, gravity_[2], state_);
+                                          props, deck_, gravity_[2], *state_);
             }
 
             // Threshold pressures.
             std::map<std::pair<int, int>, double> maxDp;
-            computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]);
+            computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]);
             threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp);
             std::vector<double> threshold_pressures_nnc = thresholdPressuresNNC(eclipse_state_, geoprops_->nnc(), maxDp);
             threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end());
@@ -497,9 +533,9 @@ namespace Opm
                 const int numCells = Opm::UgGridHelpers::numCells(grid);
                 std::vector<int> cells(numCells);
                 for (int c = 0; c < numCells; ++c) { cells[c] = c; }
-                std::vector<double> pc = state_.saturation();
-                props.capPress(numCells, state_.saturation().data(), cells.data(), pc.data(), nullptr);
-                fluidprops_->setSwatInitScaling(state_.saturation(), pc);
+                std::vector<double> pc = state_->saturation();
+                props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr);
+                fluidprops_->setSwatInitScaling(state_->saturation(), pc);
             }
         }
 
@@ -522,7 +558,7 @@ namespace Opm
             // and initilialize new properties and states for it.
             if (must_distribute_) {
                 distributeGridAndData(grid_init_->grid(), deck_, eclipse_state_,
-                                      state_, *fluidprops_, *geoprops_,
+                                      *state_, *fluidprops_, *geoprops_,
                                       material_law_manager_, threshold_pressures_,
                                       parallel_information_, use_local_perm_);
             }
@@ -621,7 +657,7 @@ namespace Opm
                               << std::flush;
                 }
 
-                SimulatorReport fullReport = simulator_->run(simtimer, state_);
+                SimulatorReport fullReport = simulator_->run(simtimer, *state_);
 
                 if (output_cout_) {
                     std::cout << "\n\n================    End of simulation     ===============\n\n";

opm/autodiff/FlowMainPolymer.hpp

@@ -38,6 +38,10 @@ namespace Opm
     {
     protected:
         using Base = FlowMainBase<FlowMainPolymer<Grid, Simulator>, Grid, Simulator>;
+        using Base::eclipse_state_;
+        using Base::param_;
+        using Base::fis_solver_;
+        using Base::parallel_information_;
         friend Base;
 
         // Set in setupGridAndProps()
@@ -90,10 +94,34 @@ namespace Opm
         // Setup linear solver.
         // Writes to:
         //   fis_solver_
+        // Currently, the CPR solver is not ready for polymer solver yet
         void setupLinearSolver()
         {
-            OPM_MESSAGE("Caution: polymer solver currently only works with direct linear solver.");
-            Base::fis_solver_.reset(new NewtonIterationBlackoilSimple(Base::param_, Base::parallel_information_));
+            const std::string cprSolver = "cpr";
+            const std::string interleavedSolver = "interleaved";
+            const std::string directSolver = "direct";
+            const std::string flowDefaultSolver = interleavedSolver;
+
+            std::shared_ptr<const Opm::SimulationConfig> simCfg = eclipse_state_->getSimulationConfig();
+            std::string solver_approach = flowDefaultSolver;
+
+            if (param_.has("solver_approach")) {
+                solver_approach = param_.template get<std::string>("solver_approach");
+            }  else {
+                if (simCfg->useCPR()) {
+                    solver_approach = cprSolver;
+                }
+            }
+
+            if (solver_approach == cprSolver) {
+                OPM_THROW( std::runtime_error , "CPR solver is not ready for use with polymer solver yet.");
+            } else if (solver_approach == interleavedSolver) {
+                fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_));
+            } else if (solver_approach == directSolver) {
+                fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_));
+            } else {
+                OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
+            }
         }
 
 

opm/autodiff/FlowMainSolvent.hpp

@@ -107,7 +107,6 @@ namespace Opm
                                                  Base::deck_->hasKeyword("SOLVENT")));
         }
 
-
     };
 
 

opm/autodiff/ParallelDebugOutput.hpp

@@ -19,8 +19,10 @@
 #ifndef OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
 #define OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+
 #include <opm/core/grid.h>
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/wells/WellsManager.hpp>
 
@@ -35,6 +37,7 @@ namespace Opm
 {
 
     class ParallelDebugOutputInterface
+
     {
     protected:
         ParallelDebugOutputInterface () {}
@@ -42,11 +45,11 @@ namespace Opm
         virtual ~ParallelDebugOutputInterface() {}
 
         // gather solution to rank 0 for EclipseWriter
-        virtual bool collectToIORank( const SimulatorState& localReservoirState,
+        virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
                                       const int reportStep ) = 0;
 
-        virtual const SimulatorState& globalReservoirState() const = 0 ;
+        virtual const SimulationDataContainer& globalReservoirState() const = 0 ;
         virtual const WellState& globalWellState() const = 0 ;
         virtual bool isIORank() const = 0;
         virtual bool isParallel() const = 0;
@@ -60,7 +63,7 @@ namespace Opm
     protected:
         const GridImpl& grid_;
 
-        const SimulatorState* globalState_;
+        const SimulationDataContainer* globalState_;
         const WellState*      wellState_;
 
     public:
@@ -71,7 +74,7 @@ namespace Opm
             : grid_( grid ) {}
 
         // gather solution to rank 0 for EclipseWriter
-        virtual bool collectToIORank( const SimulatorState& localReservoirState,
+        virtual bool collectToIORank( const SimulationDataContainer& localReservoirState,
                                       const WellState& localWellState,
                                       const int /* reportStep */)
         {
@@ -80,7 +83,7 @@ namespace Opm
             return true ;
         }
 
-        virtual const SimulatorState& globalReservoirState() const { return *globalState_; }
+        virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; }
         virtual const WellState& globalWellState() const { return *wellState_; }
         virtual bool isIORank () const { return true; }
         virtual bool isParallel () const { return false; }
@@ -243,7 +246,7 @@ namespace Opm
                     Dune::CpGrid& globalGrid = *grid_;
 
                     // initialize global state with correct sizes
-                    globalReservoirState_.init( globalGrid.numCells(), globalGrid.numFaces(), numPhases );
+                    globalReservoirState_.reset( new SimulationDataContainer( globalGrid.numCells(), globalGrid.numFaces(), numPhases ));
 
                     // copy global cartesian index
                     globalIndex_ = globalGrid.globalCell();
@@ -306,18 +309,18 @@ namespace Opm
             }
         }
 
-        class PackUnPackSimulatorState : public P2PCommunicatorType::DataHandleInterface
+        class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface
         {
-            const SimulatorState& localState_;
-            SimulatorState& globalState_;
+            const SimulationDataContainer& localState_;
+            SimulationDataContainer& globalState_;
             const WellState& localWellState_;
             WellState& globalWellState_;
             const IndexMapType& localIndexMap_;
             const IndexMapStorageType& indexMaps_;
 
         public:
-            PackUnPackSimulatorState( const SimulatorState& localState,
-                                      SimulatorState& globalState,
+            PackUnPackSimulationDataContainer( const SimulationDataContainer& localState,
+                                      SimulationDataContainer& globalState,
                                       const WellState& localWellState,
                                       WellState& globalWellState,
                                       const IndexMapType& localIndexMap,
@@ -333,12 +336,11 @@ namespace Opm
                 if( isIORank )
                 {
                     // add missing data to global state
-                    for( size_t i=globalState_.cellData().size();
-                         i<localState.cellData().size(); ++i )
-                    {
-                        const size_t components = localState.cellData()[ i ].size() / localState.numCells();
-                        assert( components * localState.numCells() == localState.cellData()[ i ].size() );
-                        globalState_.registerCellData( localState.cellDataNames()[ i ], components );
+   		    for (const auto& pair : localState.cellData()) {
+		        const std::string& key = pair.first;
+                        if (!globalState_.hasCellData( key )) {
+			    globalState_.registerCellData( key , localState.numCellDataComponents( key ));
+                        }
                     }
 
                     MessageBufferType buffer;
@@ -357,13 +359,10 @@ namespace Opm
                 }
 
                 // write all cell data registered in local state
-                const size_t numCells = localState_.numCells();
-                const size_t numCellData = localState_.cellData().size();
-                for( size_t d=0; d<numCellData; ++d )
-                {
-                    const std::vector< double >& data = localState_.cellData()[ d ];
-                    const size_t stride = data.size() / numCells ;
-                    assert( numCells * stride == data.size() );
+		for (const auto& pair : localState_.cellData()) {
+		    const std::string& key = pair.first;
+                    const auto& data = pair.second;
+                    const size_t stride = localState_.numCellDataComponents( key );
 
                     for( size_t i=0; i<stride; ++i )
                     {
@@ -378,22 +377,20 @@ namespace Opm
 
             void doUnpack( const IndexMapType& indexMap, MessageBufferType& buffer )
             {
-                // read all cell data registered in local state
-                const size_t numCells = globalState_.numCells();
-                const size_t numCellData = globalState_.cellData().size();
-                for( size_t d=0; d<numCellData; ++d )
-                {
-                    std::vector< double >& data = globalState_.cellData()[ d ];
-                    const size_t stride = data.size() / numCells ;
-                    assert( numCells * stride == data.size() );
+                // write all cell data registered in local state
+	        for (auto& pair : globalState_.cellData()) {
+                    const std::string& key = pair.first;
+		    auto& data = pair.second;
+                    const size_t stride = globalState_.numCellDataComponents( key );
 
                     for( size_t i=0; i<stride; ++i )
                     {
-                        // write all data from local state to buffer
-                        read( buffer, indexMap, data, i, stride );
+                        //write all data from local state to buffer
+		        read( buffer, indexMap, data, i, stride );
                     }
                 }
 
+
                 // read well data from buffer
                 readWells( buffer );
             }
@@ -516,7 +513,7 @@ namespace Opm
         };
 
         // gather solution to rank 0 for EclipseWriter
-        bool collectToIORank( const SimulatorState& localReservoirState,
+        bool collectToIORank( const SimulationDataContainer& localReservoirState,
                               const WellState& localWellState,
                               const int reportStep )
         {
@@ -536,10 +533,10 @@ namespace Opm
                                            false);
 
                 const Wells* wells = wells_manager.c_wells();
-                globalWellState_.init(wells, globalReservoirState_, globalWellState_ );
+                globalWellState_.init(wells, *globalReservoirState_, globalWellState_ );
             }
 
-            PackUnPackSimulatorState packUnpack( localReservoirState, globalReservoirState_,
+            PackUnPackSimulationDataContainer packUnpack( localReservoirState, *globalReservoirState_,
                                                  localWellState, globalWellState_,
                                                  localIndexMap_, indexMaps_, isIORank() );
 
@@ -552,7 +549,7 @@ namespace Opm
             return isIORank();
         }
 
-        const SimulatorState& globalReservoirState() const { return globalReservoirState_; }
+        const SimulationDataContainer& globalReservoirState() const { return *globalReservoirState_; }
         const WellState& globalWellState() const { return globalWellState_; }
 
         bool isIORank() const
@@ -573,18 +570,18 @@ namespace Opm
         }
 
     protected:
-        std::unique_ptr< Dune::CpGrid > grid_;
-        Opm::EclipseStateConstPtr       eclipseState_;
-        const double*                   permeability_;
-        P2PCommunicatorType             toIORankComm_;
-        IndexMapType                    globalIndex_;
-        IndexMapType                    localIndexMap_;
-        IndexMapStorageType             indexMaps_;
-        SimulatorState                  globalReservoirState_;
+        std::unique_ptr< Dune::CpGrid >           grid_;
+        Opm::EclipseStateConstPtr                 eclipseState_;
+        const double*                             permeability_;
+        P2PCommunicatorType                       toIORankComm_;
+        IndexMapType                              globalIndex_;
+        IndexMapType                              localIndexMap_;
+        IndexMapStorageType                       indexMaps_;
+        std::unique_ptr<SimulationDataContainer>  globalReservoirState_;
         // this needs to be revised
-        WellStateFullyImplicitBlackoil  globalWellState_;
+        WellStateFullyImplicitBlackoil            globalWellState_;
         // true if we are on I/O rank
-        const bool                      isIORank_;
+        const bool                                isIORank_;
     };
 #endif // #if HAVE_DUNE_CORNERPOINT
 

opm/autodiff/RedistributeDataHandles.hpp

@@ -232,7 +232,7 @@ public:
     {
         assert( T::codimension == 0);
 
-        for ( int i=0; i<sendState_.numPhases(); ++i )
+        for ( size_t i=0; i<sendState_.numPhases(); ++i )
         {
             buffer.write(sendState_.surfacevol()[e.index()*sendState_.numPhases()+i]);
         }
@@ -240,7 +240,7 @@ public:
         buffer.write(sendState_.rv()[e.index()]);
         buffer.write(sendState_.pressure()[e.index()]);
         buffer.write(sendState_.temperature()[e.index()]);
-        for ( int i=0; i<sendState_.numPhases(); ++i )
+        for ( size_t i=0; i<sendState_.numPhases(); ++i )
         {
             buffer.write(sendState_.saturation()[e.index()*sendState_.numPhases()+i]);
         }
@@ -257,11 +257,11 @@ public:
     void scatter(B& buffer, const T& e, std::size_t size_arg)
     {
         assert( T::codimension == 0);
-        assert( int(size_arg) == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
+        assert( size_arg == 2 * recvState_.numPhases() +4+2*recvGrid_.numCellFaces(e.index()));
         static_cast<void>(size_arg);
 
         double val;
-        for ( int i=0; i<recvState_.numPhases(); ++i )
+        for ( size_t i=0; i<recvState_.numPhases(); ++i )
         {
             buffer.read(val);
             recvState_.surfacevol()[e.index()*sendState_.numPhases()+i]=val;
@@ -274,7 +274,7 @@ public:
         recvState_.pressure()[e.index()]=val;
         buffer.read(val);
         recvState_.temperature()[e.index()]=val;
-        for ( int i=0; i<recvState_.numPhases(); ++i )
+        for ( size_t i=0; i<recvState_.numPhases(); ++i )
         {
             buffer.read(val);
             recvState_.saturation()[e.index()*sendState_.numPhases()+i]=val;
@@ -411,8 +411,6 @@ void distributeGridAndData( Dune::CpGrid& grid,
                             boost::any& parallelInformation,
                             const bool useLocalPerm)
 {
-    BlackoilState distributed_state;
-
     Dune::CpGrid global_grid ( grid );
     global_grid.switchToGlobalView();
 
@@ -441,8 +439,8 @@ void distributeGridAndData( Dune::CpGrid& grid,
     BlackoilPropsAdFromDeck distributed_props(properties,
                                               distributed_material_law_manager,
                                               grid.numCells());
-    distributed_state.init(grid.numCells(), grid.numFaces(), state.numPhases());
-    // init does not resize surfacevol. Do it manually.
+    BlackoilState distributed_state(grid.numCells(), grid.numFaces(), state.numPhases());
+    // construction does not resize surfacevol. Do it manually.
     distributed_state.surfacevol().resize(grid.numCells()*state.numPhases(),
                                           std::numeric_limits<double>::max());
     BlackoilStateDataHandle state_handle(global_grid, grid,

opm/autodiff/SimulatorBase.hpp

@@ -40,7 +40,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 #include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp

@@ -21,8 +21,13 @@
 
 #include "SimulatorFullyImplicitBlackoilOutput.hpp"
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
+#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/utility/DataMap.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/Units.hpp>
@@ -30,9 +35,6 @@
 #include <opm/autodiff/GridHelpers.hpp>
 #include <opm/autodiff/BackupRestore.hpp>
 
-#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
-
-
 #include <sstream>
 #include <iomanip>
 #include <fstream>
@@ -50,7 +52,7 @@ namespace Opm
 
 
     void outputStateVtk(const UnstructuredGrid& grid,
-                        const SimulatorState& state,
+                        const SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir)
     {
@@ -87,16 +89,15 @@ namespace Opm
 
 
     void outputStateMatlab(const UnstructuredGrid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir)
     {
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
-        for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
-        {
-            const std::string& name = state.cellDataNames()[ i ];
+        for (const auto& pair : state.cellData()) {
+            const std::string& name = pair.first;
             std::string key;
             if( name == "SURFACEVOL" ) {
                 key = "surfvolume";
@@ -111,7 +112,7 @@ namespace Opm
                 continue;
             }
             // set data to datmap
-            dm[ key ] = &state.cellData()[ i ];
+            dm[ key ] = &pair.second;
         }
 
         std::vector<double> cell_velocity;
@@ -204,7 +205,7 @@ namespace Opm
 
 #ifdef HAVE_DUNE_CORNERPOINT
     void outputStateVtk(const Dune::CpGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir)
     {
@@ -255,7 +256,7 @@ namespace Opm
     void
     BlackoilOutputWriter::
     writeTimeStep(const SimulatorTimerInterface& timer,
-                  const SimulatorState& localState,
+                  const SimulationDataContainer& localState,
                   const WellState& localWellState,
                   bool substep)
     {
@@ -271,7 +272,7 @@ namespace Opm
             isIORank = parallelOutput_->collectToIORank( localState, localWellState, timer.reportStepNum() );
         }
 
-        const SimulatorState& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
+        const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState;
         const WellState& wellState  = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState;
 
         // output is only done on I/O rank
@@ -306,6 +307,7 @@ namespace Opm
                     // write resport step number
                     backupfile_.write( (const char *) &reportStep, sizeof(int) );
 
+                    /*
                     try {
                         const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state);
                         backupfile_ << boState;
@@ -317,6 +319,7 @@ namespace Opm
                     {
 
                     }
+                    */
 
                     /*
                     const WellStateFullyImplicitBlackoil* boWellState =

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -20,16 +20,15 @@
 #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
 #define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED
 #include <opm/core/grid.h>
-#include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/DataMap.hpp>
 #include <opm/common/ErrorMacros.hpp>
-#include <opm/core/io/eclipse/EclipseReader.hpp>
+#include <opm/output/eclipse/EclipseReader.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
-#include <opm/core/io/OutputWriter.hpp>
-#include <opm/core/io/eclipse/EclipseWriter.hpp>
+#include <opm/output/OutputWriter.hpp>
+#include <opm/output/eclipse/EclipseWriter.hpp>
 
 #include <opm/autodiff/GridHelpers.hpp>
 #include <opm/autodiff/ParallelDebugOutput.hpp>
@@ -53,14 +52,17 @@
 namespace Opm
 {
 
+    class SimulationDataContainer;
+    class BlackoilState;
+
     void outputStateVtk(const UnstructuredGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir);
 
 
     void outputStateMatlab(const UnstructuredGrid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir);
 
@@ -69,23 +71,23 @@ namespace Opm
                                const std::string& output_dir);
 #ifdef HAVE_DUNE_CORNERPOINT
     void outputStateVtk(const Dune::CpGrid& grid,
-                        const Opm::SimulatorState& state,
+                        const Opm::SimulationDataContainer& state,
                         const int step,
                         const std::string& output_dir);
 #endif
 
     template<class Grid>
     void outputStateMatlab(const Grid& grid,
-                           const Opm::SimulatorState& state,
+                           const Opm::SimulationDataContainer& state,
                            const int step,
                            const std::string& output_dir)
     {
         Opm::DataMap dm;
         dm["saturation"] = &state.saturation();
         dm["pressure"] = &state.pressure();
-        for( unsigned int i=0; i<state.cellDataNames().size(); ++ i )
+        for (const auto& pair : state.cellData()) 
         {
-            const std::string& name = state.cellDataNames()[ i ];
+            const std::string& name = pair.first;
             std::string key;
             if( name == "SURFACEVOL" ) {
                 key = "surfvolume";
@@ -100,7 +102,7 @@ namespace Opm
                 continue;
             }
             // set data to datmap
-            dm[ key ] = &state.cellData()[ i ];
+            dm[ key ] = &pair.second;
         }
 
         std::vector<double> cell_velocity;
@@ -154,7 +156,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& state,
+                           const SimulationDataContainer& state,
                            const WellState&,
                            bool /*substep*/ = false)
         {
@@ -184,7 +186,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& reservoirState,
+                           const SimulationDataContainer& reservoirState,
                            const WellState& wellState,
                            bool /*substep*/ = false)
         {
@@ -214,7 +216,7 @@ namespace Opm
 
         /** \copydoc Opm::OutputWriter::writeTimeStep */
         void writeTimeStep(const SimulatorTimerInterface& timer,
-                           const SimulatorState& reservoirState,
+                           const SimulationDataContainer& reservoirState,
                            const Opm::WellState& wellState,
                            bool substep = false);
 
@@ -235,7 +237,7 @@ namespace Opm
         void initFromRestartFile(const PhaseUsage& phaseusage,
                                  const double* permeability,
                                  const Grid& grid,
-                                 SimulatorState& simulatorstate,
+                                 SimulationDataContainer& simulatorstate,
                                  WellStateFullyImplicitBlackoil& wellstate);
 
         bool isRestart() const;
@@ -315,7 +317,7 @@ namespace Opm
     initFromRestartFile( const PhaseUsage& phaseusage,
                          const double* permeability,
                          const Grid& grid,
-                         SimulatorState& simulatorstate,
+                         SimulationDataContainer& simulatorstate,
                          WellStateFullyImplicitBlackoil& wellstate)
     {
         WellsManager wellsmanager(eclipseState_,

opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp

@@ -44,7 +44,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 //#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 

opm/autodiff/SimulatorIncompTwophaseAd.cpp

@@ -37,7 +37,7 @@
 #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/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 
 #include <opm/core/wells/WellsManager.hpp>

opm/autodiff/SolventPropsAdFromDeck.cpp

@@ -21,12 +21,13 @@
 #include <opm/autodiff/SolventPropsAdFromDeck.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
 
+#include <opm/core/utility/extractPvtTableIndex.hpp>
+
 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/PvdsTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
 
-
 namespace Opm
 {
 // Making these typedef to make the code more readable.
@@ -172,6 +173,25 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck,
                 OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n");
             }
 
+            const TableContainer& pmiscTables = tables->getPmiscTables();
+            if (!pmiscTables.empty()) {
+
+                int numRegions = pmiscTables.size();
+
+                // resize the attributes of the object
+                pmisc_.resize(numRegions);
+                for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
+                    const Opm::PmiscTable& pmiscTable = pmiscTables.getTable<PmiscTable>(regionIdx);
+
+                    // Copy data
+                    const auto& po = pmiscTable.getOilPhasePressureColumn();
+                    const auto& pmisc = pmiscTable.getMiscibilityColumn();
+
+                    pmisc_[regionIdx] = NonuniformTableLinear<double>(po, pmisc);
+
+                }
+            }
+
             // miscible relative permeability multipleiers
             const TableContainer& msfnTables = tables->getMsfnTables();
             if (!msfnTables.empty()) {
@@ -339,6 +359,16 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction,
     return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_);
 }
 
+ADB SolventPropsAdFromDeck::pressureMiscibilityFunction(const ADB& po,
+                                                        const Cells& cells) const
+{
+    if (pmisc_.size() > 0) {
+        return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, pmisc_);
+    }
+    // return ones if not specified i.e. no effect.
+    return ADB::constant(V::Constant(po.size(), 1.0));
+}
+
 
 ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw,
                                                                    const Cells& cells) const {
@@ -383,7 +413,6 @@ ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD,
 }
 
 
-
 V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const {
     const int n = cells.size();
     V density(n);

opm/autodiff/SolventPropsAdFromDeck.hpp

@@ -23,8 +23,7 @@
 #include <opm/autodiff/BlackoilPropsAdInterface.hpp>
 #include <opm/autodiff/AutoDiffBlock.hpp>
 
-#include <opm/core/props/pvt/PvtDead.hpp>
-#include <opm/core/props/pvt/PvtInterface.hpp>
+#include <opm/core/utility/NonuniformTableLinear.hpp>
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -108,6 +107,13 @@ public:
     ADB miscibilityFunction(const ADB& solventFraction,
                             const Cells& cells) const;
 
+    /// Pressure dependent miscibility function
+    /// \param[in]  solventFraction Array of n oil phase pressure .
+    /// \param[in]  cells           Array of n cell indices to be associated with the pressure values.
+    /// \return                     Array of n miscibility values
+    ADB pressureMiscibilityFunction(const ADB& po,
+                                    const Cells& cells) const;
+
     /// Miscible critical gas saturation function
     /// \param[in]  Sw              Array of n water saturation values.
     /// \param[in]  cells           Array of n cell indices to be associated with the saturation values.
@@ -139,7 +145,6 @@ public:
 
 
 private:
-
     /// Makes ADB from table values
     /// \param[in]  X               Array of n table lookup values.
     /// \param[in]  cells           Array of n cell indices to be associated with the fraction values.
@@ -177,6 +182,7 @@ private:
     std::vector<NonuniformTableLinear<double> > mkro_;
     std::vector<NonuniformTableLinear<double> > mkrsg_;
     std::vector<NonuniformTableLinear<double> > misc_;
+    std::vector<NonuniformTableLinear<double> > pmisc_;
     std::vector<NonuniformTableLinear<double> > sorwmis_;
     std::vector<NonuniformTableLinear<double> > sgcwmis_;
     std::vector<double> mix_param_viscosity_;

opm/autodiff/TransportSolverTwophaseAd.cpp

@@ -22,6 +22,7 @@
 #include <opm/autodiff/TransportSolverTwophaseAd.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/linalg/LinearSolverInterface.hpp>
+#include <opm/core/props/IncompPropertiesInterface.hpp>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/common/ErrorMacros.hpp>

opm/polymer/CompressibleTpfaPolymer.cpp

@@ -88,8 +88,8 @@ namespace Opm
                                   PolymerBlackoilState& state,
                                   WellState& well_state)
     {
-        c_ = &state.concentration();
-        cmax_ = &state.maxconcentration();
+        c_ = &state.getCellData( state.CONCENTRATION );
+        cmax_ = &state.getCellData( state.CMAX );
         CompressibleTpfa::solve(dt, state, well_state);
     }
 

opm/polymer/IncompTpfaPolymer.cpp

@@ -20,6 +20,8 @@
 
 #include <config.h>
 
+#include <opm/common/data/SimulationDataContainer.hpp>
+
 #include <opm/polymer/IncompTpfaPolymer.hpp>
 
 #include <opm/core/props/IncompPropertiesInterface.hpp>
@@ -99,8 +101,8 @@ namespace Opm
                                   PolymerState& state,
                                   WellState& well_state)
     {
-        c_ = &state.concentration();
-        cmax_ = &state.maxconcentration();
+        c_ = &state.getCellData( state.CONCENTRATION );
+        cmax_ = &state.getCellData( state.CMAX) ;
         if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) {
             solveRockComp(dt, state, well_state);
         } else {
@@ -116,7 +118,7 @@ namespace Opm
 
     /// Compute per-solve dynamic properties.
     void IncompTpfaPolymer::computePerSolveDynamicData(const double /*dt*/,
-                                                       const SimulatorState& state,
+                                                       const SimulationDataContainer& state,
                                                        const WellState& /*well_state*/)
     {
         // Computed here:

opm/polymer/IncompTpfaPolymer.hpp

@@ -37,6 +37,7 @@ namespace Opm
     class LinearSolverInterface;
     class PolymerState;
     class WellState;
+    class SimulationDataContainer;
 
     /// Encapsulating a tpfa pressure solver for the incompressible-fluid case with polymer.
     /// Supports gravity, wells controlled by bhp or reservoir rates,
@@ -96,7 +97,7 @@ namespace Opm
 
     private:
         virtual void computePerSolveDynamicData(const double dt,
-                                                const SimulatorState& state,
+                                                const SimulationDataContainer& state,
                                                 const WellState& well_state);
     private:
         // ------ Data that will remain unmodified after construction. ------

opm/polymer/PolymerBlackoilState.cpp

@@ -0,0 +1,42 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  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/common/data/SimulationDataContainer.hpp>
+
+#include <opm/polymer/PolymerBlackoilState.hpp>
+
+namespace Opm
+{
+    const std::string PolymerBlackoilState::CONCENTRATION = "CONCENTRATION";
+    const std::string PolymerBlackoilState::CMAX = "CMAX";
+
+    PolymerBlackoilState::PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases) :
+        BlackoilState( number_of_cells , number_of_faces, num_phases)
+    {
+        registerCellData(CONCENTRATION , 1 , 0 );
+        registerCellData(CMAX , 1 , 0 );
+    }
+
+
+} // namespace Opm
+
+
+
+

opm/polymer/PolymerBlackoilState.hpp

@@ -33,27 +33,10 @@ namespace Opm
     class PolymerBlackoilState : public BlackoilState
     {
     public:
-        void init(const UnstructuredGrid& g, int num_phases)
-        {
-            this->init(g.number_of_cells, g.number_of_faces, num_phases);
-        }
+        static const std::string CONCENTRATION;
+        static const std::string CMAX;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            BlackoilState::init(number_of_cells, number_of_faces, num_phases);
-            concentration_.resize(number_of_cells, 0.0);
-            cmax_.resize(number_of_cells, 0.0);
-        }
-
-        std::vector<double>& concentration()    { return concentration_; }
-        std::vector<double>& maxconcentration() { return cmax_; }
-
-        const std::vector<double>& concentration() const    { return concentration_; }
-        const std::vector<double>& maxconcentration() const { return cmax_; }
-
-    private:
-        std::vector<double> concentration_;
-        std::vector<double> cmax_;
+        PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases);
     };
 
 } // namespace Opm

opm/polymer/PolymerState.cpp

@@ -0,0 +1,42 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  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/common/data/SimulationDataContainer.hpp>
+
+#include <opm/polymer/PolymerState.hpp>
+
+namespace Opm
+{
+    const std::string PolymerState::CONCENTRATION = "CONCENTRATION";
+    const std::string PolymerState::CMAX = "CMAX";
+
+    PolymerState::PolymerState(int number_of_cells, int number_of_faces, int num_phases) :
+        SimulationDataContainer( number_of_cells , number_of_faces , num_phases )
+    {
+        registerCellData(CONCENTRATION , 1 , 0 );
+        registerCellData(CMAX , 1 , 0 );
+    }
+
+
+} // namespace Opm
+
+
+
+

opm/polymer/PolymerState.hpp

@@ -20,8 +20,8 @@
 #ifndef OPM_POLYMERSTATE_HEADER_INCLUDED
 #define OPM_POLYMERSTATE_HEADER_INCLUDED
 
+#include <opm/common/data/SimulationDataContainer.hpp>
 
-#include <opm/core/simulator/SimulatorState.hpp>
 #include <opm/core/grid.h>
 #include <vector>
 
@@ -29,34 +29,13 @@ namespace Opm
 {
 
     /// Simulator state for a two-phase simulator with polymer.
-    class PolymerState : public SimulatorState
+    class PolymerState : public SimulationDataContainer
     {
     public:
+        static const std::string CONCENTRATION;
+        static const std::string CMAX;
 
-        void init(int number_of_cells, int number_of_faces, int num_phases)
-        {
-            SimulatorState::init( number_of_cells , number_of_faces , num_phases );
-            registerCellData("CONCENTRATION" , 1 , 0 );
-            registerCellData("CMAX" , 1 , 0 );
-        }
-
-        const std::vector<double>& concentration() const {
-            return getCellData("CONCENTRATION");
-        }
-
-        const std::vector<double>& maxconcentration() const {
-            return getCellData("CMAX");
-        }
-
-
-        std::vector<double>& concentration()  {
-            return getCellData("CONCENTRATION");
-        }
-
-        std::vector<double>& maxconcentration()  {
-            return getCellData("CMAX");
-        }
-
+        PolymerState(int number_of_cells, int number_of_faces, int num_phases);
     };
 
 } // namespace Opm

opm/polymer/SimulatorCompressiblePolymer.cpp

@@ -34,7 +34,7 @@
 #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/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
@@ -268,7 +268,7 @@ namespace Opm
         total_timer.start();
         double init_surfvol[2] = { 0.0 };
         double inplace_surfvol[2] = { 0.0 };
-        double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
         double polymass_adsorbed = computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_);
         double init_polymass = polymass + polymass_adsorbed;
         double tot_injected[2] = { 0.0 };
@@ -301,7 +301,7 @@ namespace Opm
         if (check_well_controls_) {
             computeFractionalFlow(props_, poly_props_, allcells_,
                                   state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
-                                  state.concentration(), state.maxconcentration(),
+                                  state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ) ,
                                   fractional_flows);
             wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
         }
@@ -397,7 +397,7 @@ namespace Opm
                            state.pressure(), state.temperature(), &initial_porevol[0], &porevol[0],
                            &transport_src[0], &polymer_inflow_c[0], stepsize,
                            state.saturation(), state.surfacevol(),
-                           state.concentration(), state.maxconcentration());
+                           state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             double substep_injected[2] = { 0.0 };
             double substep_produced[2] = { 0.0 };
             double substep_polyinj = 0.0;
@@ -416,7 +416,7 @@ namespace Opm
             if (gravity_ != 0 && use_segregation_split_) {
                 tsolver_.solveGravity(columns_, stepsize,
                                       state.saturation(), state.surfacevol(),
-                                      state.concentration(), state.maxconcentration());
+                                      state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             }
         }
         transport_timer.stop();
@@ -426,7 +426,7 @@ namespace Opm
 
         // Report volume balances.
         Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
-        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
         polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
                                                         state, rock_comp_props_);
         tot_injected[0] += injected[0];
@@ -539,8 +539,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             dm["surfvol"] = &state.surfacevol();
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
@@ -556,8 +556,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             dm["surfvol"] = &state.surfacevol();
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);

opm/polymer/SimulatorPolymer.cpp

@@ -36,7 +36,7 @@
 #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/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 
 #include <opm/core/wells/WellsManager.hpp>
@@ -62,7 +62,7 @@
 #include <iostream>
 
 #ifdef HAVE_ERT
-#include <opm/core/io/eclipse/writeECLData.hpp>
+#include <opm/output/eclipse/writeECLData.hpp>
 #endif
 
 
@@ -292,8 +292,8 @@ namespace Opm
         total_timer.start();
         double init_satvol[2] = { 0.0 };
         double satvol[2] = { 0.0 };
-        double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
-        double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
+        double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
+        double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
         double init_polymass = polymass + polymass_adsorbed;
         double injected[2] = { 0.0 };
         double produced[2] = { 0.0 };
@@ -330,7 +330,7 @@ namespace Opm
         // Solve pressure.
         if (check_well_controls_) {
             computeFractionalFlow(props_, poly_props_, allcells_,
-                                  state.saturation(), state.concentration(), state.maxconcentration(),
+                                  state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ),
                                   fractional_flows);
             wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
         }
@@ -425,7 +425,7 @@ namespace Opm
         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());
+                           state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             Opm::computeInjectedProduced(props_, poly_props_,
                                          state,
                                          transport_src, polymer_inflow_c, stepsize,
@@ -438,7 +438,7 @@ namespace Opm
             polyprod += substep_polyprod;
             if (use_segregation_split_) {
                 tsolver_.solveGravity(columns_, &porevol[0], stepsize,
-                                      state.saturation(), state.concentration(), state.maxconcentration());
+                                      state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ));
             }
         }
         transport_timer.stop();
@@ -448,8 +448,8 @@ namespace Opm
 
         // 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());
+        polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol());
+        polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX ));
         tot_injected[0] += injected[0];
         tot_injected[1] += injected[1];
         tot_produced[0] += produced[0];
@@ -559,8 +559,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;
@@ -575,8 +575,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;
@@ -611,8 +611,8 @@ namespace Opm
             Opm::DataMap dm;
             dm["saturation"] = &state.saturation();
             dm["pressure"] = &state.pressure();
-            dm["concentration"] = &state.concentration();
-            dm["cmax"] = &state.maxconcentration();
+            dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ;
+            dm["cmax"] = &state.getCellData( state.CMAX ) ;
             std::vector<double> cell_velocity;
             Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
             dm["velocity"] = &cell_velocity;

opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp

@@ -234,6 +234,8 @@ namespace Opm {
                                            const SolutionState& state,
                                            WellState& xw);
 
+        void updateEquationsScaling();
+
         void
         computeMassFlux(const int               actph ,
                         const V&                transi,

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -107,6 +107,7 @@ namespace Opm {
             if (!active_[Water]) {
                 OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
             }
+            residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase
             // If deck has polymer, residual_ should contain polymer equation.
             rq_.resize(fluid_.numPhases() + 1);
             residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
@@ -125,8 +126,10 @@ namespace Opm {
                 WellState& well_state)
     {
         Base::prepareStep(dt, reservoir_state, well_state);
+        auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX );
         // Initial max concentration of this time step from PolymerBlackoilState.
-        cmax_ = Eigen::Map<const V>(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_));
+
+        cmax_ = Eigen::Map<const V>(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_));
     }
 
 
@@ -168,9 +171,10 @@ namespace Opm {
 
         // Initial polymer concentration.
         if (has_polymer_) {
-            assert (not x.concentration().empty());
-            const int nc = x.concentration().size();
-            const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
+            const auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert (not concentration.empty());
+            const int nc = concentration.size();
+            const V c = Eigen::Map<const V>(concentration.data() , nc);
             // Concentration belongs after other reservoir vars but before well vars.
             auto concentration_pos = vars0.begin() + fluid_.numPhases();
             assert(concentration_pos == vars0.end() - 2);
@@ -255,12 +259,14 @@ namespace Opm {
     template <class Grid>
     void BlackoilPolymerModel<Grid>::computeCmax(ReservoirState& state)
     {
-        const int nc = AutoDiffGrid::numCells(grid_);
-        V tmp = V::Zero(nc);
-        for (int i = 0; i < nc; ++i) {
-            tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
-        }
-        std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
+        auto& max_concentration = state.getCellData( state.CMAX );
+        const auto& concentration = state.getCellData( state.CONCENTRATION );
+        std::transform( max_concentration.begin() ,
+                        max_concentration.end() ,
+                        concentration.begin() ,
+                        max_concentration.begin() ,
+                        [](double c_max , double c) { return std::max( c_max , c ); });
+
     }
 
 
@@ -339,6 +345,26 @@ namespace Opm {
             residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
                                                + ops_.div*rq_[poly_pos_].mflux;
         }
+
+
+        if (param_.update_equations_scaling_) {
+            updateEquationsScaling();
+        }
+
+    }
+
+
+
+
+
+    template <class Grid>
+    void BlackoilPolymerModel<Grid>::updateEquationsScaling()
+    {
+        Base::updateEquationsScaling();
+        if (has_polymer_) {
+            const int water_pos = fluid_.phaseUsage().phase_pos[Water];
+            residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
+        }
     }
 
 
@@ -397,10 +423,13 @@ namespace Opm {
             // Call base version.
             Base::updateState(modified_dx, reservoir_state, well_state);
 
-            // Update concentration.
-            const V c_old = Eigen::Map<const V>(&reservoir_state.concentration()[0], nc, 1);
-            const V c = (c_old - dc).max(zero);
-            std::copy(&c[0], &c[0] + nc, reservoir_state.concentration().begin());
+            {
+                auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION );
+                // Update concentration.
+                const V c_old = Eigen::Map<const V>(concentration.data(), nc, 1);
+                const V c = (c_old - dc).max(zero);
+                std::copy(&c[0], &c[0] + nc, concentration.begin());
+            }
         } else {
             // Just forward call to base version.
             Base::updateState(dx, reservoir_state, well_state);

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp

@@ -206,7 +206,7 @@ namespace {
         const std::vector<double>& polymer_inflow = xw.polymerInflow();
 
         // Initial max concentration of this time step from PolymerBlackoilState.
-        cmax_ = Eigen::Map<V>(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_));
+        cmax_ = Eigen::Map<V>(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_));
 
         const SolutionState state = constantState(x, xw);
         computeAccum(state, 0);
@@ -366,9 +366,18 @@ namespace {
         state.saturation[1] = ADB::constant(so, bpat);
 
         // Concentration
-        assert(not x.concentration().empty());
-        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
-        state.concentration = ADB::constant(c);
+        {
+            auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert(concentration.empty());
+            const V c = Eigen::Map<const V>(concentration.data(), nc);
+            // Do not understand:
+            //concentration = ADB::constant(c);
+            // Old code based on concentraton() method had the statement:
+            //
+            //   state.concentration = ADB::constant(c)
+            //
+            // This looks like it was a method assignment - how did it even compile?
+        }
 
         // Well rates.
         assert (not xw.wellRates().empty());
@@ -413,9 +422,12 @@ namespace {
         vars0.push_back(sw0);
 
         // Initial concentration.
-        assert (not x.concentration().empty());
-        const V c = Eigen::Map<const V>(&x.concentration()[0], nc);
-        vars0.push_back(c);
+        {
+            auto& concentration = x.getCellData( x.CONCENTRATION );
+            assert (not concentration.empty());
+            const V c = Eigen::Map<const V>(concentration.data() , nc);
+            vars0.push_back(c);
+        }
 
         // Initial well rates.
         assert (not xw.wellRates().empty());
@@ -511,11 +523,14 @@ namespace {
     {
         const int nc = grid_.number_of_cells;
         V tmp = V::Zero(nc);
+        const auto& concentration = state.getCellData( state.CONCENTRATION );
+        auto& cmax = state.getCellData( state.CMAX );
+
         for (int i = 0; i < nc; ++i) {
-            tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]);
+            tmp[i] = std::max(cmax[i], concentration[i]);
         }
 
-        std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin());
+        std::copy(&tmp[0], &tmp[0] + nc, cmax.begin());
     }
 
 
@@ -764,11 +779,14 @@ namespace {
         // Concentration updates.
 //        const double dcmax = 0.3 * polymer_props_ad_.cMax();
 //		std::cout << "\n the max concentration: " << dcmax / 0.3 << std::endl;
-        const V c_old = Eigen::Map<const V>(&state.concentration()[0], nc, 1);
-//        const V dc_limited = sign(dc) * dc.abs().min(dcmax);
-//        const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
-        const V c = (c_old - dc).max(zero);
-        std::copy(&c[0], &c[0] + nc, state.concentration().begin());
+        {
+            auto& concentration = state.getCellData( state.CONCENTRATION );
+            const V c_old = Eigen::Map<const V>(concentration.data() , nc, 1);
+            //        const V dc_limited = sign(dc) * dc.abs().min(dcmax);
+            //        const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02());
+            const V c = (c_old - dc).max(zero);
+            std::copy(&c[0], &c[0] + nc, concentration.begin());
+        }
 
         // Qs update.
         // Since we need to update the wellrates, that are ordered by wells,

opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp

@@ -45,7 +45,7 @@
 #include <opm/core/simulator/SimulatorTimer.hpp>
 //#include <opm/core/simulator/AdaptiveSimulatorTimer.hpp>
 #include <opm/core/utility/StopWatch.hpp>
-#include <opm/core/io/vtk/writeVtkData.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 

opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp

@@ -38,8 +38,8 @@
 #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/output/eclipse/EclipseWriter.hpp>
+#include <opm/output/vtk/writeVtkData.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 

opm/polymer/polymerUtilities.cpp

@@ -211,8 +211,8 @@ namespace Opm
             OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells.");
         }
         const std::vector<double>& s = state.saturation();
-        const std::vector<double>& c = state.concentration();
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
         std::fill(injected, injected + np, 0.0);
         std::fill(produced, produced + np, 0.0);
         polyinj = 0.0;
@@ -282,8 +282,8 @@ namespace Opm
         const std::vector<double>& temp = state.temperature();
         const std::vector<double>& s = state.saturation();
         const std::vector<double>& z = state.surfacevol();
-        const std::vector<double>& c = state.concentration();
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& c = state.getCellData( state.CONCENTRATION );
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
         std::fill(injected, injected + np, 0.0);
         std::fill(produced, produced + np, 0.0);
         polyinj = 0.0;
@@ -406,7 +406,7 @@ namespace Opm
             porosity.assign(props.porosity(), props.porosity() + num_cells);
         }
         double abs_mass = 0.0;
-        const std::vector<double>& cmax = state.maxconcentration();
+        const std::vector<double>& cmax = state.getCellData( state.CMAX );
 	for (int cell = 0; cell < num_cells; ++cell) {
             double c_ads;
             polyprops.simpleAdsorption(cmax[cell], c_ads);

tests/test_block.cpp

@@ -1,5 +1,6 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
+  Copyright 2016 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -293,7 +294,7 @@ BOOST_AUTO_TEST_CASE(Pow)
     vx << 0.2, 1.2, 13.4;
 
     ADB::V vy(3);
-    vy << 1.0, 2.2, 3.4;
+    vy << 2.0, 3.0, 0.5;
 
     std::vector<ADB::V> vals{ vx, vy };
     std::vector<ADB> vars = ADB::variables(vals);
@@ -303,6 +304,7 @@ BOOST_AUTO_TEST_CASE(Pow)
 
     const double tolerance = 1e-14;
 
+    // test exp = double
     const ADB xx = x * x;
     ADB xxpow2 = Opm::pow(x,2.0);
     checkClose(xxpow2, xx, tolerance);
@@ -322,6 +324,48 @@ BOOST_AUTO_TEST_CASE(Pow)
     for (int i = 0 ; i < 3; ++i){
         BOOST_CHECK_CLOSE(xpowhalf.value()[i], x_sqrt[i], 1e-4);
     }
+
+    // test exp = ADB::V
+    ADB xpowyval = Opm::pow(x,y.value());
+
+    // each of the component of y is tested in the test above
+    // we compare with the results from the above tests.
+    ADB::V pick1(3);
+    pick1 << 1,0,0;
+    ADB::V pick2(3);
+    pick2 << 0,1,0;
+    ADB::V pick3(3);
+    pick3 << 0,0,1;
+
+    ADB compare = pick1 * xx + pick2 * xxx + pick3 * xpowhalf;
+    checkClose(xpowyval, compare, tolerance);
+
+    // test exponent = ADB::V and base = ADB
+    ADB xvalpowy = Opm::pow(x.value(),y);
+
+    // the value should be equal to xpowyval
+    // the first jacobian should be trivial
+    // the second jacobian is hand calculated
+    // log(0.2)*0.2^2.0, log(1.2) * 1.2^3.0, log(13.4) * 13.4^0.5
+    ADB::V jac2(3);
+    jac2 << -0.0643775165 , 0.315051650 , 9.50019208855;
+    for (int i = 0 ; i < 3; ++i){
+        BOOST_CHECK_CLOSE(xvalpowy.value()[i], xpowyval.value()[i], tolerance);
+        BOOST_CHECK_CLOSE(xvalpowy.derivative()[0].coeff(i,i), 0.0, tolerance);
+        BOOST_CHECK_CLOSE(xvalpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4);
+    }
+
+    // test exp = ADB
+    ADB xpowy = Opm::pow(x,y);
+
+    // the first jacobian should be equal to the xpowyval
+    // the second jacobian should be equal to the xvalpowy
+    for (int i = 0 ; i < 3; ++i){
+        BOOST_CHECK_CLOSE(xpowy.value()[i], xpowyval.value()[i], tolerance);
+        BOOST_CHECK_CLOSE(xpowy.derivative()[0].coeff(i,i), xpowyval.derivative()[0].coeff(i,i), tolerance);
+        BOOST_CHECK_CLOSE(xpowy.derivative()[1].coeff(i,i), xvalpowy.derivative()[1].coeff(i,i), tolerance);
+    }
+
 }
 
 

tests/test_boprops_ad.cpp

@@ -37,7 +37,7 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 
@@ -48,10 +48,10 @@
 struct SetupSimple {
     SetupSimple()
     {
-        Opm::ParseMode parseMode;
+        Opm::ParseContext parseContext;
         Opm::ParserPtr parser(new Opm::Parser());
-        deck = parser->parseFile("fluid.data", parseMode);
-        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        deck = parser->parseFile("fluid.data", parseContext);
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
 
         param.disableOutput();
         param.insertParameter("init_rock"       , "false" );

tests/test_rateconverter.cpp

@@ -32,6 +32,7 @@
 
 #include <boost/test/unit_test.hpp>
 
+#include <opm/core/grid/GridHelpers.hpp>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
 #include <opm/core/utility/Units.hpp>
@@ -39,7 +40,7 @@
 #include <opm/core/simulator/BlackoilState.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 
@@ -48,9 +49,9 @@ struct SetupSimple {
     SetupSimple()
     {
         Opm::ParserPtr parser(new Opm::Parser());
-        Opm::ParseMode parseMode;
-        deck = parser->parseFile("fluid.data" , parseMode);
-        eclState.reset(new Opm::EclipseState(deck , parseMode));
+        Opm::ParseContext parseContext;
+        deck = parser->parseFile("fluid.data" , parseContext);
+        eclState.reset(new Opm::EclipseState(deck , parseContext));
 
         param.disableOutput();
         param.insertParameter("init_rock"       , "false" );
@@ -107,8 +108,7 @@ BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture<SetupSimple>)
     Region reg{ 0 };
     RCvrt  cvrt(ad_props, reg);
 
-    Opm::BlackoilState x;
-    x.init(*grid.c_grid(), 3);
+    Opm::BlackoilState x( Opm::UgGridHelpers::numCells( *grid.c_grid()) , Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3);
 
     cvrt.defineState(x);
 

tests/test_singlecellsolves.cpp

@@ -71,7 +71,7 @@ try
 
     boost::scoped_ptr<GridManager> grid;
     boost::scoped_ptr<IncompPropertiesInterface> props;
-    PolymerState state;
+    std::unique_ptr<PolymerState> state;
     Opm::PolymerProperties poly_props;
     // bool check_well_controls = false;
     // int max_well_control_iterations = 0;
@@ -81,23 +81,35 @@ try
     // Grid init.
     grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0));
     // Rock and fluid init.
-    props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
-    // Init state variables (saturation and pressure).
-    initStateBasic(*grid->c_grid(), *props, param, 0.0, state);
-    // Init Polymer state
-    if (param.has("poly_init")) {
-        double poly_init = param.getDefault("poly_init", 0.0);
-        for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
-            double smin[2], smax[2];
-            props->satRange(1, &cell, smin, smax);
-            if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) {
-                state.concentration()[cell] = poly_init;
-                state.maxconcentration()[cell] = poly_init;
-            } else {
-                state.saturation()[2*cell + 0] = 0.;
-                state.saturation()[2*cell + 1] = 1.;
-                state.concentration()[cell] = 0.;
-                state.maxconcentration()[cell] = 0.;
+    {
+        const UnstructuredGrid& ug_grid = *(grid->c_grid());
+        state.reset( new PolymerState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2));
+
+        props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));
+        // Init state variables (saturation and pressure).
+        initStateBasic(*grid->c_grid(), *props, param, 0.0, *state);
+        // Init Polymer state
+        if (param.has("poly_init")) {
+            double poly_init = param.getDefault("poly_init", 0.0);
+            for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) {
+                double smin[2], smax[2];
+                props->satRange(1, &cell, smin, smax);
+
+                auto& saturation = state->saturation();
+                auto& concentration = state->getCellData( state->CONCENTRATION );
+                auto& max_concentration = state->getCellData( state->CMAX );
+
+                props->satRange(1, &cell, smin, smax);
+                if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) {
+                    concentration[cell] = poly_init;
+                    max_concentration[cell] = poly_init;
+                } else {
+                    saturation[2*cell + 0] = 0.;
+                    saturation[2*cell + 1] = 1.;
+                    concentration[cell] = 0.;
+                    max_concentration[cell] = 0.;
+                }
+
             }
         }
     }
@@ -210,7 +222,7 @@ try
     if (face01 == -1) {
         OPM_THROW(std::runtime_error, "Could not find face adjacent to cells [0 1]");
     }
-    state.faceflux()[face01] = src[0];
+    state->faceflux()[face01] = src[0];
     for (int sats = 0; sats < num_sats; ++sats) {
         const double s = double(sats)/double(num_sats - 1);
         const double ff = s; // Simplified a lot...
@@ -220,20 +232,26 @@ try
             // std::cout << "(s, c) = (" << s << ", " << c << ")\n";
             transport_src[0] = src[0]*ff;
             // Resetting the state for next run.
-            state.saturation()[0] = 0.0;
-            state.saturation()[1] = 0.0;
-            state.concentration()[0] = 0.0;
-            state.concentration()[1] = 0.0;
-            state.maxconcentration()[0] = 0.0;
-            state.maxconcentration()[1] = 0.0;
-            reorder_model.solve(&state.faceflux()[0],
+            auto& saturation = state->saturation();
+            auto& concentration = state->getCellData( state->CONCENTRATION );
+            auto& max_concentration = state->getCellData( state->CMAX );
+
+
+            saturation[0] = 0.0;
+            saturation[1] = 0.0;
+            concentration[0] = 0.0;
+            concentration[1] = 0.0;
+            max_concentration[0] = 0.0;
+            max_concentration[1] = 0.0;
+
+            reorder_model.solve(&state->faceflux()[0],
                                 &porevol[0],
                                 &transport_src[0],
                                 &polymer_inflow_c[0],
                                 dt,
-                                state.saturation(),
-                                state.concentration(),
-                                state.maxconcentration());
+                                saturation,
+                                concentration,
+                                max_concentration);
 
 #ifdef PROFILING
             // Extract residual counts.

tests/test_transmissibilitymultipliers.cpp

@@ -31,7 +31,7 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 #include <opm/parser/eclipse/Parser/Parser.hpp>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
@@ -158,12 +158,12 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
 {
     Opm::parameter::ParameterGroup param;
     Opm::ParserPtr parser(new Opm::Parser() );
-    Opm::ParseMode parseMode;
+    Opm::ParseContext parseContext;
 
     /////
     // create a DerivedGeology object without any multipliers involved
-    Opm::DeckConstPtr origDeck = parser->parseString(origDeckString, parseMode);
-    Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseMode));
+    Opm::DeckConstPtr origDeck = parser->parseString(origDeckString, parseContext);
+    Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseContext));
 
     auto origGridManager = std::make_shared<Opm::GridManager>(origEclipseState->getEclipseGrid());
     auto origProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(origDeck, origEclipseState, *(origGridManager->c_grid()));
@@ -173,8 +173,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
 
     /////
     // create a DerivedGeology object _with_ transmissibility multipliers involved
-    Opm::DeckConstPtr multDeck = parser->parseString(multDeckString, parseMode);
-    Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseMode));
+    Opm::DeckConstPtr multDeck = parser->parseString(multDeckString, parseContext);
+    Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseContext));
 
     auto multGridManager = std::make_shared<Opm::GridManager>(multEclipseState->getEclipseGrid());
     auto multProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multDeck, multEclipseState, *(multGridManager->c_grid()));
@@ -185,8 +185,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
     /////
     // create a DerivedGeology object _with_ transmissibility multipliers involved for
     // the negative faces
-    Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString, parseMode);
-    Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck , parseMode));
+    Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString, parseContext);
+    Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck , parseContext));
 
     auto multMinusGridManager = std::make_shared<Opm::GridManager>(multMinusEclipseState->getEclipseGrid());
     auto multMinusProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(multMinusDeck, multMinusEclipseState, *(multMinusGridManager->c_grid()));
@@ -196,8 +196,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface)
 
     /////
     // create a DerivedGeology object with the NTG keyword involved
-    Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseMode);
-    Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseMode));
+    Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseContext);
+    Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseContext));
 
     auto ntgGridManager = std::make_shared<Opm::GridManager>(ntgEclipseState->getEclipseGrid());
     auto ntgProps = std::make_shared<Opm::BlackoilPropsAdFromDeck>(ntgDeck, ntgEclipseState, *(ntgGridManager->c_grid()));
@@ -273,12 +273,12 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
 
     Opm::parameter::ParameterGroup param;
     Opm::ParserPtr parser(new Opm::Parser() );
-    Opm::ParseMode parseMode;
+    Opm::ParseContext parseContext;
 
     /////
     // create a DerivedGeology object without any multipliers involved
-    Opm::DeckConstPtr origDeck = parser->parseString(origDeckString , parseMode);
-    Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseMode));
+    Opm::DeckConstPtr origDeck = parser->parseString(origDeckString , parseContext);
+    Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseContext));
 
     auto origGrid = std::make_shared<Dune::CpGrid>();
     origGrid->processEclipseFormat(origEclipseState->getEclipseGrid(), 0.0, false);
@@ -292,8 +292,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
 
     /////
     // create a DerivedGeology object _with_ transmissibility multipliers involved
-    Opm::DeckConstPtr multDeck = parser->parseString(multDeckString,parseMode);
-    Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseMode));
+    Opm::DeckConstPtr multDeck = parser->parseString(multDeckString,parseContext);
+    Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseContext));
 
     auto multGrid = std::make_shared<Dune::CpGrid>();
     multGrid->processEclipseFormat(multEclipseState->getEclipseGrid(), 0.0, false);
@@ -306,8 +306,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
     /////
     // create a DerivedGeology object _with_ transmissibility multipliers involved for
     // the negative faces
-    Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString , parseMode);
-    Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck, parseMode));
+    Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString , parseContext);
+    Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck, parseContext));
 
     auto multMinusGrid = std::make_shared<Dune::CpGrid>();
     multMinusGrid->processEclipseFormat(multMinusEclipseState->getEclipseGrid(), 0.0, false);
@@ -320,8 +320,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid)
 
     /////
     // create a DerivedGeology object with the NTG keyword involved
-    Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseMode);
-    Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseMode));
+    Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseContext);
+    Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseContext));
 
     auto ntgGrid = std::make_shared<Dune::CpGrid>();
     ntgGrid->processEclipseFormat(ntgEclipseState->getEclipseGrid(), 0.0, false);

tests/test_vfpproperties.cpp

@@ -39,7 +39,7 @@
 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 #include <opm/core/wells.h>
-#include <opm/parser/eclipse/Parser/ParseMode.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
 #include <opm/parser/eclipse/Parser/Parser.hpp>
 #include <opm/parser/eclipse/EclipseState/checkDeck.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -1047,7 +1047,7 @@ VFPPROD \n\
     std::shared_ptr<Opm::UnitSystem> units(Opm::UnitSystem::newFIELD());
 
     Opm::ParserPtr parser(new Opm::Parser());
-    Opm::ParseMode parse_mode;
+    Opm::ParseContext parse_mode;
     deck = parser->parseString(table_str, parse_mode);
 
     BOOST_REQUIRE(deck->hasKeyword("VFPPROD"));
@@ -1108,7 +1108,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD)
     std::shared_ptr<Opm::UnitSystem> units(Opm::UnitSystem::newMETRIC());
 
     Opm::ParserPtr parser(new Opm::Parser());
-    Opm::ParseMode parse_mode;
+    Opm::ParseContext parse_mode;
     boost::filesystem::path file("VFPPROD2");
 
     deck = parser->parseFile(file.string(), parse_mode);