Merge pull request #642 from akva2/import_files_from_opm_core

Import some files from opm core (with history)

CMakeLists_files.cmake

@@ -63,6 +63,8 @@ list (APPEND MAIN_SOURCE_FILES
   opm/polymer/TransportSolverTwophasePolymer.cpp
   opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
   opm/polymer/fullyimplicit/PolymerPropsAd.cpp
+  opm/simulators/SimulatorCompressibleTwophase.cpp
+  opm/simulators/SimulatorIncompTwophase.cpp
   )
 
 
@@ -82,6 +84,7 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_vfpproperties.cpp
   tests/test_singlecellsolves.cpp
   tests/test_solventprops_ad.cpp
+	# tests/test_thresholdpressure.cpp
   )
 
 list (APPEND TEST_DATA_FILES
@@ -99,6 +102,7 @@ list (APPEND EXAMPLE_SOURCE_FILES
   examples/flow_multisegment.cpp
   examples/flow_solvent.cpp
   examples/sim_2p_incomp_ad.cpp
+  examples/sim_2p_comp_reorder.cpp
   examples/sim_simple.cpp
   examples/opm_init_check.cpp
   examples/sim_poly2p_comp_reorder.cpp

examples/sim_2p_comp_reorder.cpp

@@ -0,0 +1,300 @@
+/*
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/pressure/FlowBCManager.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/core/grid/GridHelpers.hpp>
+#include <opm/core/wells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/core/simulator/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/props/BlackoilPropertiesBasic.hpp>
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/simulator/SimulatorCompressibleTwophase.hpp>
+
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
+
+#include <memory>
+#include <boost/filesystem.hpp>
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <numeric>
+
+
+namespace
+{
+    void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
+    {
+        if (param.anyUnused()) {
+            std::cout << "--------------------   Unused parameters:   --------------------\n";
+            param.displayUsage();
+            std::cout << "----------------------------------------------------------------" << std::endl;
+        }
+    }
+} // anon namespace
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+try
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for weakly compressible two-phase flow     ===============\n\n";
+    parameter::ParameterGroup param(argc, argv, false);
+    std::cout << "---------------    Reading parameters     ---------------" << std::endl;
+
+    // If we have a "deck_filename", grid and props will be read from that.
+    bool use_deck = param.has("deck_filename");
+    EclipseStateConstPtr eclipseState;
+    std::unique_ptr<GridManager> grid;
+    std::unique_ptr<BlackoilPropertiesInterface> props;
+    std::unique_ptr<RockCompressibility> rock_comp;
+    std::unique_ptr<BlackoilState> state;
+
+
+    ParserPtr parser(new Opm::Parser());
+    Opm::DeckConstPtr deck;
+
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    if (use_deck) {
+        ParseContext parseContext;
+        std::string deck_filename = param.get<std::string>("deck_filename");
+        deck = parser->parseFile(deck_filename , parseContext);
+        eclipseState.reset(new EclipseState(deck, parseContext));
+
+        // Grid init
+        grid.reset(new GridManager(deck));
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+            state.reset( new BlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) ,2));
+            // Rock and fluid init
+            props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, ug_grid, param));
+            // 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(ug_grid, *props, param, gravity[2], *state);
+            } else {
+                initStateFromDeck(ug_grid, *props, deck, gravity[2], *state);
+            }
+            initBlackoilSurfvol(ug_grid, *props, *state);
+        }
+    } else {
+        // Grid init.
+        const int nx = param.getDefault("nx", 100);
+        const int ny = param.getDefault("ny", 100);
+        const int nz = param.getDefault("nz", 1);
+        const double dx = param.getDefault("dx", 1.0);
+        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));
+        {
+            const UnstructuredGrid& ug_grid = *(grid->c_grid());
+            // Rock and fluid init.
+            props.reset(new BlackoilPropertiesBasic(param, ug_grid.dimensions,  UgGridHelpers::numCells( ug_grid )));
+
+            // State init
+            state.reset( new BlackoilState(  UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 3));
+            // 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);
+        }
+    }
+
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+    // Initialising src
+    int num_cells = grid->c_grid()->number_of_cells;
+    std::vector<double> src(num_cells, 0.0);
+    if (use_deck) {
+        // Do nothing, wells will be the driving force, not source terms.
+    } else {
+        // Compute pore volumes, in order to enable specifying injection rate
+        // terms of total pore volume.
+        std::vector<double> porevol;
+        if (rock_comp->isActive()) {
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol);
+        } else {
+            computePorevolume(*grid->c_grid(), props->porosity(), porevol);
+        }
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        const double default_injection = use_gravity ? 0.0 : 0.1;
+        const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
+            *tot_porevol_init/unit::day;
+        src[0] = flow_per_sec;
+        src[num_cells - 1] = -flow_per_sec;
+    }
+
+    // Boundary conditions.
+    FlowBCManager bcs;
+    if (param.getDefault("use_pside", false)) {
+        int pside = param.get<int>("pside");
+        double pside_pressure = param.get<double>("pside_pressure");
+        bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
+    }
+
+    // Linear solver.
+    LinearSolverFactory linsolver(param);
+
+    // Write parameters used for later reference.
+    bool output = param.getDefault("output", true);
+    std::ofstream epoch_os;
+    std::string output_dir;
+    if (output) {
+        output_dir =
+            param.getDefault("output_dir", std::string("output"));
+        boost::filesystem::path fpath(output_dir);
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        std::string filename = output_dir + "/epoch_timing.param";
+        epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        // open file to clean it. The file is appended to in SimulatorTwophase
+        filename = output_dir + "/step_timing.param";
+        std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
+        step_os.close();
+        param.writeParam(output_dir + "/simulation.param");
+    }
+
+
+    std::cout << "\n\n================    Starting main simulation loop     ===============\n";
+
+    SimulatorReport rep;
+    if (!use_deck) {
+        // Simple simulation without a deck.
+        WellsManager wells; // no wells.
+        SimulatorCompressibleTwophase simulator(param,
+                                                *grid->c_grid(),
+                                                *props,
+                                                rock_comp->isActive() ? rock_comp.get() : 0,
+                                                wells,
+                                                src,
+                                                bcs.c_bcs(),
+                                                linsolver,
+                                                grav);
+        SimulatorTimer simtimer;
+        simtimer.init(param);
+        warnIfUnusedParams(param);
+        WellState 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;
+        int step = 0;
+        SimulatorTimer simtimer;
+        // Use timer for last epoch to obtain total time.
+        Opm::TimeMapPtr timeMap(new Opm::TimeMap(deck));
+        simtimer.init(timeMap);
+        const double total_time = simtimer.totalTime();
+        for (size_t reportStepIdx = 0; reportStepIdx < timeMap->numTimesteps(); ++reportStepIdx) {
+            simtimer.setCurrentStepNum(step);
+            simtimer.setTotalTime(total_time);
+
+            // Report on start of report step.
+            std::cout << "\n\n--------------    Starting report step " << reportStepIdx << "    --------------"
+                      << "\n                  (number of steps: "
+                      << simtimer.numSteps() - step << ")\n\n" << std::flush;
+
+            // Create new wells, well_state
+            WellsManager wells(eclipseState , reportStepIdx , *grid->c_grid(), props->permeability());
+            // @@@ HACK: we should really make a new well state and
+            // 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);
+            }
+
+            // Create and run simulator.
+            SimulatorCompressibleTwophase simulator(param,
+                                                    *grid->c_grid(),
+                                                    *props,
+                                                    rock_comp->isActive() ? rock_comp.get() : 0,
+                                                    wells,
+                                                    src,
+                                                    bcs.c_bcs(),
+                                                    linsolver,
+                                                    grav);
+            if (reportStepIdx == 0) {
+                warnIfUnusedParams(param);
+            }
+            SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state);
+            if (output) {
+                epoch_rep.reportParam(epoch_os);
+            }
+            // Update total timing report and remember step number.
+            rep += epoch_rep;
+            step = simtimer.currentStepNum();
+        }
+    }
+
+    std::cout << "\n\n================    End of simulation     ===============\n\n";
+    rep.report(std::cout);
+
+    if (output) {
+        std::string filename = output_dir + "/walltime.param";
+        std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out);
+        rep.reportParam(tot_os);
+    }
+
+}
+catch (const std::exception &e) {
+    std::cerr << "Program threw an exception: " << e.what() << "\n";
+    throw;
+}

opm/simulators/SimulatorCompressibleTwophase.cpp

@@ -0,0 +1,555 @@
+/*
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/simulator/SimulatorCompressibleTwophase.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/core/pressure/CompressibleTpfa.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#endif
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
+
+#include <opm/core/wells/WellsManager.hpp>
+
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/grid/ColumnExtract.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
+
+#include <boost/filesystem.hpp>
+#include <memory>
+#include <boost/lexical_cast.hpp>
+
+#include <numeric>
+#include <fstream>
+#include <iostream>
+
+
+namespace Opm
+{
+
+    class SimulatorCompressibleTwophase::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const BlackoilPropertiesInterface& props,
+             const RockCompressibility* rock_comp_props,
+             WellsManager& wells_manager,
+             const std::vector<double>& src,
+             const FlowBoundaryConditions* bcs,
+             LinearSolverInterface& linsolver,
+             const double* gravity);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            BlackoilState& state,
+                            WellState& well_state);
+
+    private:
+        // Data.
+
+        // Parameters for output.
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for well control
+        bool check_well_controls_;
+        int max_well_control_iterations_;
+        // Parameters for transport solver.
+        int num_transport_substeps_;
+        bool use_segregation_split_;
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const BlackoilPropertiesInterface& props_;
+        const RockCompressibility* rock_comp_props_;
+        WellsManager& wells_manager_;
+        const Wells* wells_;
+        //const std::vector<double>& src_;
+        //const FlowBoundaryConditions* bcs_;
+        const double* gravity_;
+        // Solvers
+        CompressibleTpfa psolver_;
+        TransportSolverCompressibleTwophaseReorder tsolver_;
+        // Needed by column-based gravity segregation solver.
+        std::vector< std::vector<int> > columns_;
+        // Misc. data
+        std::vector<int> allcells_;
+    };
+
+
+
+
+    SimulatorCompressibleTwophase::SimulatorCompressibleTwophase(const parameter::ParameterGroup& param,
+                                                                 const UnstructuredGrid& grid,
+                                                                 const BlackoilPropertiesInterface& props,
+                                                                 const RockCompressibility* rock_comp_props,
+                                                                 WellsManager& wells_manager,
+                                                                 const std::vector<double>& src,
+                                                                 const FlowBoundaryConditions* bcs,
+                                                                 LinearSolverInterface& linsolver,
+                                                                 const double* gravity)
+    {
+        pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorCompressibleTwophase::run(SimulatorTimer& timer,
+                                                       BlackoilState& state,
+                                                       WellState& well_state)
+    {
+        return pimpl_->run(timer, state, well_state);
+    }
+
+
+
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+    static void outputStateVtk(const UnstructuredGrid& grid,
+                               const Opm::BlackoilState& state,
+                               const int step,
+                               const std::string& output_dir)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        vtkfilename << output_dir << "/vtk_files";
+        boost::filesystem::path fpath(vtkfilename.str());
+        try {
+          create_directories(fpath);
+        }
+        catch (...) {
+          OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        if (!vtkfile) {
+            OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
+        }
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+#endif
+
+
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+    static void outputStateMatlab(const UnstructuredGrid& grid,
+                                  const Opm::BlackoilState& state,
+                                  const int step,
+                                  const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        dm["surfvolume"] = &state.surfacevol();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+              create_directories(fpath);
+            }
+            catch (...) {
+              OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+#endif
+
+
+    static void outputWaterCut(const Opm::Watercut& watercut,
+                               const std::string& output_dir)
+    {
+        // Write water cut curve.
+        std::string fname = output_dir  + "/watercut.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        watercut.write(os);
+    }
+
+
+    static void outputWellReport(const Opm::WellReport& wellreport,
+                                 const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+
+
+
+    // \TODO: make CompressibleTpfa take src and bcs.
+    SimulatorCompressibleTwophase::Impl::Impl(const parameter::ParameterGroup& param,
+                                              const UnstructuredGrid& grid,
+                                              const BlackoilPropertiesInterface& props,
+                                              const RockCompressibility* rock_comp_props,
+                                              WellsManager& wells_manager,
+                                              const std::vector<double>&,
+                                              const FlowBoundaryConditions*,
+                                              LinearSolverInterface& linsolver,
+                                              const double* gravity)
+        : grid_(grid),
+          props_(props),
+          rock_comp_props_(rock_comp_props),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
+          //src_(src),
+          //bcs_(bcs),
+          gravity_(gravity),
+          psolver_(grid, props, rock_comp_props, linsolver,
+                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
+                   param.getDefault("nl_pressure_change_tolerance", 1.0),
+                   param.getDefault("nl_pressure_maxiter", 10),
+                   gravity, wells_manager.c_wells() /*, src, bcs*/),
+          tsolver_(grid, props,
+                   param.getDefault("nl_tolerance", 1e-9),
+                   param.getDefault("nl_maxiter", 30))
+    {
+        // For output.
+        output_ = param.getDefault("output", true);
+        if (output_) {
+            output_vtk_ = param.getDefault("output_vtk", true);
+            output_dir_ = param.getDefault("output_dir", std::string("output"));
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(output_dir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            output_interval_ = param.getDefault("output_interval", 1);
+        }
+
+        // Well control related init.
+        check_well_controls_ = param.getDefault("check_well_controls", false);
+        max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
+
+        // Transport related init.
+        num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
+        use_segregation_split_ = param.getDefault("use_segregation_split", false);
+        if (gravity != 0 && use_segregation_split_){
+            tsolver_.initGravity(gravity);
+            extractColumn(grid_, columns_);
+        }
+
+        // Misc init.
+        const int num_cells = grid.number_of_cells;
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+
+
+
+    SimulatorReport SimulatorCompressibleTwophase::Impl::run(SimulatorTimer& timer,
+                                                             BlackoilState& state,
+                                                             WellState& well_state)
+    {
+        std::vector<double> transport_src;
+
+        // Initialisation.
+        std::vector<double> porevol;
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+        } else {
+            computePorevolume(grid_, props_.porosity(), porevol);
+        }
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        std::vector<double> initial_porevol = porevol;
+
+        // Main simulation loop.
+        Opm::time::StopWatch pressure_timer;
+        double ptime = 0.0;
+        Opm::time::StopWatch transport_timer;
+        double ttime = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        double init_surfvol[2] = { 0.0 };
+        double inplace_surfvol[2] = { 0.0 };
+        double tot_injected[2] = { 0.0 };
+        double tot_produced[2] = { 0.0 };
+        Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
+        Opm::Watercut watercut;
+        watercut.push(0.0, 0.0, 0.0);
+        Opm::WellReport wellreport;
+        std::vector<double> fractional_flows;
+        std::vector<double> well_resflows_phase;
+        if (wells_) {
+            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
+            wellreport.push(props_, *wells_,
+                            state.pressure(), state.surfacevol(), state.saturation(),
+                            0.0, well_state.bhp(), well_state.perfRates());
+        }
+        std::fstream tstep_os;
+        if (output_) {
+            std::string filename = output_dir_ + "/step_timing.param";
+            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
+        }
+        for (; !timer.done(); ++timer) {
+            // Report timestep and (optionally) write state to disk.
+            step_timer.start();
+            timer.report(std::cout);
+            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+                if (output_vtk_) {
+#ifdef DISABLE_OUTPUT
+                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+                }
+#ifdef DISABLE_OUTPUT
+                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+            }
+
+            SimulatorReport sreport;
+
+            // Solve pressure equation.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, allcells_,
+                                      state.pressure(), state.temperature(), state.surfacevol(), state.saturation(),
+                                      fractional_flows);
+                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
+            }
+            bool well_control_passed = !check_well_controls_;
+            int well_control_iteration = 0;
+            do {
+                // Run solver.
+                pressure_timer.start();
+                std::vector<double> initial_pressure = state.pressure();
+                psolver_.solve(timer.currentStepLength(), state, well_state);
+
+                // Renormalize pressure if both fluids and rock are
+                // incompressible, and there are no pressure
+                // conditions (bcs or wells).  It is deemed sufficient
+                // for now to renormalize using geometric volume
+                // instead of pore volume.
+                if (psolver_.singularPressure()) {
+                    // Compute average pressures of previous and last
+                    // step, and total volume.
+                    double av_prev_press = 0.0;
+                    double av_press = 0.0;
+                    double tot_vol = 0.0;
+                    const int num_cells = grid_.number_of_cells;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                        tot_vol       += grid_.cell_volumes[cell];
+                    }
+                    // Renormalization constant
+                    const double ren_const = (av_prev_press - av_press)/tot_vol;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        state.pressure()[cell] += ren_const;
+                    }
+                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                    for (int well = 0; well < num_wells; ++well) {
+                        well_state.bhp()[well] += ren_const;
+                    }
+                }
+
+                // Stop timer and report.
+                pressure_timer.stop();
+                double pt = pressure_timer.secsSinceStart();
+                std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
+                ptime += pt;
+                sreport.pressure_time = pt;
+
+                // Optionally, check if well controls are satisfied.
+                if (check_well_controls_) {
+                    Opm::computePhaseFlowRatesPerWell(*wells_,
+                                                      well_state.perfRates(),
+                                                      fractional_flows,
+                                                      well_resflows_phase);
+                    std::cout << "Checking well conditions." << std::endl;
+                    // For testing we set surface := reservoir
+                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
+                    ++well_control_iteration;
+                    if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
+                        OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
+                    }
+                    if (!well_control_passed) {
+                        std::cout << "Well controls not passed, solving again." << std::endl;
+                    } else {
+                        std::cout << "Well conditions met." << std::endl;
+                    }
+                }
+            } while (!well_control_passed);
+
+            // Update pore volumes if rock is compressible.
+            if (rock_comp_props_ && rock_comp_props_->isActive()) {
+                initial_porevol = porevol;
+                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+            }
+
+            // Process transport sources from well flows.
+            Opm::computeTransportSource(props_, wells_, well_state, transport_src);
+
+            // Solve transport.
+            transport_timer.start();
+            double stepsize = timer.currentStepLength();
+            if (num_transport_substeps_ != 1) {
+                stepsize /= double(num_transport_substeps_);
+                std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
+            }
+            double injected[2] = { 0.0 };
+            double produced[2] = { 0.0 };
+            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
+                tsolver_.solve(&state.faceflux()[0], &state.pressure()[0], &state.temperature()[0],
+                               &initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
+                               state.saturation(), state.surfacevol());
+                double substep_injected[2] = { 0.0 };
+                double substep_produced[2] = { 0.0 };
+                Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
+                                             substep_injected, substep_produced);
+                injected[0] += substep_injected[0];
+                injected[1] += substep_injected[1];
+                produced[0] += substep_produced[0];
+                produced[1] += substep_produced[1];
+                if (gravity_ != 0 && use_segregation_split_) {
+                    tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol());
+                }
+            }
+            transport_timer.stop();
+            double tt = transport_timer.secsSinceStart();
+            sreport.transport_time = tt;
+            std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+            ttime += tt;
+            // Report volume balances.
+            Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
+            tot_injected[0] += injected[0];
+            tot_injected[1] += injected[1];
+            tot_produced[0] += produced[0];
+            tot_produced[1] += produced[1];
+            std::cout.precision(5);
+            const int width = 18;
+            std::cout << "\nMass balance report.\n";
+            std::cout << "    Injected surface volumes:      "
+                      << std::setw(width) << injected[0]
+                      << std::setw(width) << injected[1] << std::endl;
+            std::cout << "    Produced surface volumes:      "
+                      << std::setw(width) << produced[0]
+                      << std::setw(width) << produced[1] << std::endl;
+            std::cout << "    Total inj surface volumes:     "
+                      << std::setw(width) << tot_injected[0]
+                      << std::setw(width) << tot_injected[1] << std::endl;
+            std::cout << "    Total prod surface volumes:    "
+                      << std::setw(width) << tot_produced[0]
+                      << std::setw(width) << tot_produced[1] << std::endl;
+            const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
+                                        init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] };
+            std::cout << "    Initial - inplace + inj - prod: "
+                      << std::setw(width) << balance[0]
+                      << std::setw(width) << balance[1]
+                      << std::endl;
+            std::cout << "    Relative mass error:            "
+                      << std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
+                      << std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
+                      << std::endl;
+            std::cout.precision(8);
+
+            watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
+                          produced[0]/(produced[0] + produced[1]),
+                          tot_produced[0]/tot_porevol_init);
+            if (wells_) {
+                wellreport.push(props_, *wells_,
+                                state.pressure(), state.surfacevol(), state.saturation(),
+                                timer.simulationTimeElapsed() + timer.currentStepLength(),
+                                well_state.bhp(), well_state.perfRates());
+            }
+            sreport.total_time =  step_timer.secsSinceStart();
+            if (output_) {
+                sreport.reportParam(tstep_os);
+            }
+        }
+
+        if (output_) {
+            if (output_vtk_) {
+#ifdef DISABLE_OUTPUT
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+            }
+#ifdef DISABLE_OUTPUT
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+            outputWaterCut(watercut, output_dir_);
+            if (wells_) {
+                outputWellReport(wellreport, output_dir_);
+            }
+            tstep_os.close();
+        }
+
+        total_timer.stop();
+
+        SimulatorReport report;
+        report.pressure_time = ptime;
+        report.transport_time = ttime;
+        report.total_time = total_timer.secsSinceStart();
+        return report;
+    }
+
+
+} // namespace Opm

opm/simulators/SimulatorCompressibleTwophase.hpp

@@ -0,0 +1,99 @@
+/*
+  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/>.
+*/
+
+#ifndef OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED
+#define OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED
+
+#include <memory>
+#include <vector>
+
+struct UnstructuredGrid;
+struct Wells;
+struct FlowBoundaryConditions;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class BlackoilPropertiesInterface;
+    class RockCompressibility;
+    class WellsManager;
+    class LinearSolverInterface;
+    class SimulatorTimer;
+    class BlackoilState;
+    class WellState;
+    struct SimulatorReport;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorCompressibleTwophase
+    {
+    public:
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] rock_comp_props     if non-null, rock compressibility properties
+        /// \param[in] well_manager  well manager, may manage no (null) wells
+        /// \param[in] src           source terms
+        /// \param[in] bcs           boundary conditions, treat as all noflow if null
+        /// \param[in] linsolver     linear solver
+        /// \param[in] gravity       if non-null, gravity vector
+       SimulatorCompressibleTwophase(const parameter::ParameterGroup& param,
+                                     const UnstructuredGrid& grid,
+                                     const BlackoilPropertiesInterface& props,
+                                     const RockCompressibility* rock_comp_props,
+                                     WellsManager& wells_manager,
+                                     const std::vector<double>& src,
+                                     const FlowBoundaryConditions* bcs,
+                                     LinearSolverInterface& linsolver,
+                                     const double* gravity);
+
+        /// Run the simulation.
+        /// This will run succesive timesteps until timer.done() is true. It will
+        /// modify the reservoir and well states.
+        /// \param[in,out] timer       governs the requested reporting timesteps
+        /// \param[in,out] state       state of reservoir: pressure, fluxes
+        /// \param[in,out] well_state  state of wells: bhp, perforation rates
+        /// \return                    simulation report, with timing data
+        SimulatorReport run(SimulatorTimer& timer,
+                            BlackoilState& state,
+                            WellState& well_state);
+
+    private:
+        class Impl;
+        // Using shared_ptr instead of unique_ptr since unique_ptr requires complete type for Impl.
+        std::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED

opm/simulators/SimulatorIncompTwophase.cpp

@@ -0,0 +1,666 @@
+/*
+  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/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/simulator/SimulatorIncompTwophase.hpp>
+#include <opm/core/utility/NullStream.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/common/ErrorMacros.hpp>
+
+#include <opm/core/pressure/IncompTpfa.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/wells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+#include <opm/core/io/vtk/writeVtkData.hpp>
+#endif
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/Event.hpp>
+
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/well_controls.h>
+#include <opm/core/wells.h>
+
+#include <opm/core/props/IncompPropertiesInterface.hpp>
+#include <opm/core/props/rock/RockCompressibility.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
+#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
+#include <boost/filesystem.hpp>
+#include <memory>
+
+#include <iostream>
+#include <numeric>
+#include <fstream>
+
+
+namespace Opm
+{
+
+    struct SimulatorIncompTwophase::Impl
+    {
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const IncompPropertiesInterface& props,
+             const RockCompressibility* rock_comp_props,
+             WellsManager& wells_manager,
+             const std::vector<double>& src,
+             const FlowBoundaryConditions* bcs,
+             LinearSolverInterface& linsolver,
+             const double* gravity);
+
+        SimulatorReport run(SimulatorTimer& timer,
+                            TwophaseState& state,
+                            WellState& well_state);
+
+        // Data.
+        // Parameters for output.
+        std::ostream* log_;
+        bool output_;
+        bool output_vtk_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for well control
+        bool check_well_controls_;
+        int max_well_control_iterations_;
+        // Parameters for transport solver.
+        int num_transport_substeps_;
+        bool use_reorder_;
+        bool use_segregation_split_;
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const IncompPropertiesInterface& props_;
+        const RockCompressibility* rock_comp_props_;
+        WellsManager& wells_manager_;
+        const Wells* wells_;
+        const std::vector<double>& src_;
+        const FlowBoundaryConditions* bcs_;
+        // Solvers
+        IncompTpfa psolver_;
+        std::unique_ptr<TransportSolverTwophaseInterface> tsolver_;
+        // Misc. data
+        std::vector<int> allcells_;
+
+        // list of hooks that are notified when a timestep completes
+        EventSource timestep_completed_;
+    };
+
+
+
+
+    SimulatorIncompTwophase::SimulatorIncompTwophase(const parameter::ParameterGroup& param,
+                                                     const UnstructuredGrid& grid,
+                                                     const IncompPropertiesInterface& props,
+                                                     const RockCompressibility* rock_comp_props,
+                                                     WellsManager& wells_manager,
+                                                     const std::vector<double>& src,
+                                                     const FlowBoundaryConditions* bcs,
+                                                     LinearSolverInterface& linsolver,
+                                                     const double* gravity)
+    {
+        pimpl_.reset(new Impl(param, grid, props, rock_comp_props, wells_manager, src, bcs, linsolver, gravity));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorIncompTwophase::run(SimulatorTimer& timer,
+                                                 TwophaseState& state,
+                                                 WellState& well_state)
+    {
+        return pimpl_->run(timer, state, well_state);
+    }
+
+    // connect the hook to the signal in the implementation class
+    Event& SimulatorIncompTwophase::timestep_completed () {
+        return pimpl_->timestep_completed_;
+    }
+
+    // empty default implementation, but provided in module; it is dangerous ta have
+    // this inlined in clients from the header, because then it can't be updated!
+    void SimulatorIncompTwophase::sync () {
+    }
+
+    static void reportVolumes(std::ostream &os, double satvol[2], double tot_porevol_init,
+                              double tot_injected[2], double tot_produced[2],
+                              double injected[2], double produced[2],
+                              double init_satvol[2])
+    {
+        os.precision(5);
+        const int width = 18;
+        os << "\nVolume balance report (all numbers relative to total pore volume).\n";
+        os << "    Saturated volumes:     "
+           << std::setw(width) << satvol[0]/tot_porevol_init
+           << std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
+        os << "    Injected volumes:      "
+           << std::setw(width) << injected[0]/tot_porevol_init
+           << std::setw(width) << injected[1]/tot_porevol_init << std::endl;
+        os << "    Produced volumes:      "
+           << std::setw(width) << produced[0]/tot_porevol_init
+           << std::setw(width) << produced[1]/tot_porevol_init << std::endl;
+        os << "    Total inj volumes:     "
+           << std::setw(width) << tot_injected[0]/tot_porevol_init
+           << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
+        os << "    Total prod volumes:    "
+           << std::setw(width) << tot_produced[0]/tot_porevol_init
+           << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
+        os << "    In-place + prod - inj: "
+           << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
+           << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
+        os << "    Init - now - pr + inj: "
+           << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
+           << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
+           << std::endl;
+        os.precision(8);
+    }
+
+
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+    static void outputStateVtk(const UnstructuredGrid& grid,
+                               const Opm::TwophaseState& state,
+                               const int step,
+                               const std::string& output_dir)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        vtkfilename << output_dir << "/vtk_files";
+        boost::filesystem::path fpath(vtkfilename.str());
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        if (!vtkfile) {
+            OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str());
+        }
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+        Opm::writeVtkData(grid, dm, vtkfile);
+    }
+#endif
+    static void outputVectorMatlab(const std::string& name,
+                                   const std::vector<int>& vec,
+                                   const int step,
+                                   const std::string& output_dir)
+    {
+        std::ostringstream fname;
+        fname << output_dir << "/" << name;
+        boost::filesystem::path fpath = fname.str();
+        try {
+            create_directories(fpath);
+        }
+        catch (...) {
+            OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+        }
+        fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+        std::ofstream file(fname.str().c_str());
+        if (!file) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+        }
+        std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(file, "\n"));
+    }
+
+// 17.03.2016 Temporarily removed while moving functionality to opm-output
+#ifdef DISABLE_OUTPUT
+    static void outputStateMatlab(const UnstructuredGrid& grid,
+                                  const Opm::TwophaseState& state,
+                                  const int step,
+                                  const std::string& output_dir)
+    {
+        Opm::DataMap dm;
+        dm["saturation"] = &state.saturation();
+        dm["pressure"] = &state.pressure();
+        std::vector<double> cell_velocity;
+        Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
+        dm["velocity"] = &cell_velocity;
+
+        // Write data (not grid) in Matlab format
+        for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
+            std::ostringstream fname;
+            fname << output_dir << "/" << it->first;
+            boost::filesystem::path fpath = fname.str();
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
+            }
+            file.precision(15);
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+#endif
+
+
+    static void outputWaterCut(const Opm::Watercut& watercut,
+                               const std::string& output_dir)
+    {
+        // Write water cut curve.
+        std::string fname = output_dir  + "/watercut.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        watercut.write(os);
+    }
+
+
+    static void outputWellReport(const Opm::WellReport& wellreport,
+                                 const std::string& output_dir)
+    {
+        // Write well report.
+        std::string fname = output_dir  + "/wellreport.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            OPM_THROW(std::runtime_error, "Failed to open " << fname);
+        }
+        wellreport.write(os);
+    }
+
+
+    static bool allNeumannBCs(const FlowBoundaryConditions* bcs)
+    {
+        if (bcs == NULL) {
+            return true;
+        } else {
+            return std::find(bcs->type, bcs->type + bcs->nbc, BC_PRESSURE)
+                == bcs->type + bcs->nbc;
+        }
+    }
+
+
+    static bool allRateWells(const Wells* wells)
+    {
+        if (wells == NULL) {
+            return true;
+        }
+        const int nw = wells->number_of_wells;
+        for (int w = 0; w < nw; ++w) {
+            const WellControls* wc = wells->ctrls[w];
+            if (well_controls_well_is_open( wc )) {
+                if (well_controls_get_current_type(wc) == BHP) {
+                    return false;
+                }
+            }
+        }
+        return true;
+    }
+
+
+
+
+
+    SimulatorIncompTwophase::Impl::Impl(const parameter::ParameterGroup& param,
+                                        const UnstructuredGrid& grid,
+                                        const IncompPropertiesInterface& props,
+                                        const RockCompressibility* rock_comp_props,
+                                        WellsManager& wells_manager,
+                                        const std::vector<double>& src,
+                                        const FlowBoundaryConditions* bcs,
+                                        LinearSolverInterface& linsolver,
+                                        const double* gravity)
+        : use_reorder_(param.getDefault("use_reorder", true)),
+          use_segregation_split_(param.getDefault("use_segregation_split", false)),
+          grid_(grid),
+          props_(props),
+          rock_comp_props_(rock_comp_props),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
+          src_(src),
+          bcs_(bcs),
+          psolver_(grid, props, rock_comp_props, linsolver,
+                   param.getDefault("nl_pressure_residual_tolerance", 0.0),
+                   param.getDefault("nl_pressure_change_tolerance", 1.0),
+                   param.getDefault("nl_pressure_maxiter", 10),
+                   gravity, wells_manager.c_wells(), src, bcs)
+    {
+        // Initialize transport solver.
+        if (use_reorder_) {
+            tsolver_.reset(new Opm::TransportSolverTwophaseReorder(grid,
+                                                                   props,
+                                                                   use_segregation_split_ ? gravity : NULL,
+                                                                   param.getDefault("nl_tolerance", 1e-9),
+                                                                   param.getDefault("nl_maxiter", 30)));
+
+        } else {
+            if (rock_comp_props && rock_comp_props->isActive()) {
+                OPM_THROW(std::runtime_error, "The implicit pressure solver cannot handle rock compressibility.");
+            }
+            if (use_segregation_split_) {
+                OPM_THROW(std::runtime_error, "The implicit pressure solver is not set up to use segregation splitting.");
+            }
+            std::vector<double> porevol;
+            computePorevolume(grid, props.porosity(), porevol);
+            tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(grid,
+                                                                    props,
+                                                                    porevol,
+                                                                    gravity,
+                                                                    psolver_.getHalfTrans(),
+                                                                    param));
+        }
+
+        // For output.
+        log_ = param.getDefault("quiet", false) ? &Opm::null_stream : &std::cout;
+        output_ = param.getDefault("output", true);
+        if (output_) {
+            output_vtk_ = param.getDefault("output_vtk", true);
+            output_dir_ = param.getDefault("output_dir", std::string("output"));
+            // Ensure that output dir exists
+            boost::filesystem::path fpath(output_dir_);
+            try {
+                create_directories(fpath);
+            }
+            catch (...) {
+                OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
+            }
+            output_interval_ = param.getDefault("output_interval", 1);
+        }
+
+        // Well control related init.
+        check_well_controls_ = param.getDefault("check_well_controls", false);
+        max_well_control_iterations_ = param.getDefault("max_well_control_iterations", 10);
+
+        // Transport related init.
+        num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
+
+        // Misc init.
+        const int num_cells = grid.number_of_cells;
+        allcells_.resize(num_cells);
+        for (int cell = 0; cell < num_cells; ++cell) {
+            allcells_[cell] = cell;
+        }
+    }
+
+
+
+
+    SimulatorReport SimulatorIncompTwophase::Impl::run(SimulatorTimer& timer,
+                                                       TwophaseState& state,
+                                                       WellState& well_state)
+    {
+        std::vector<double> transport_src;
+
+        // Initialisation.
+        std::vector<double> porevol;
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+        } else {
+            computePorevolume(grid_, props_.porosity(), porevol);
+        }
+        const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+        std::vector<double> initial_porevol = porevol;
+
+        // Main simulation loop.
+        Opm::time::StopWatch pressure_timer;
+        double ptime = 0.0;
+        Opm::time::StopWatch transport_timer;
+        double ttime = 0.0;
+        Opm::time::StopWatch callback_timer;
+        double time_in_callbacks = 0.0;
+        Opm::time::StopWatch step_timer;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        double init_satvol[2] = { 0.0 };
+        double satvol[2] = { 0.0 };
+        double tot_injected[2] = { 0.0 };
+        double tot_produced[2] = { 0.0 };
+        Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
+        *log_ << "\nInitial saturations are    " << init_satvol[0]/tot_porevol_init
+              << "    " << init_satvol[1]/tot_porevol_init << std::endl;
+        Opm::Watercut watercut;
+        watercut.push(0.0, 0.0, 0.0);
+        Opm::WellReport wellreport;
+        std::vector<double> fractional_flows;
+        std::vector<double> well_resflows_phase;
+        if (wells_) {
+            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
+            wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
+        }
+        std::fstream tstep_os;
+        if (output_) {
+            std::string filename = output_dir_ + "/step_timing.param";
+            tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
+        }
+        while (!timer.done()) {
+            // Report timestep and (optionally) write state to disk.
+            step_timer.start();
+            timer.report(*log_);
+            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+                if (output_vtk_) {
+#ifdef DISABLE_OUTPUT
+                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+                }
+#ifdef DISABLE_OUTPUT
+                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+                if (use_reorder_) {
+                    // This use of dynamic_cast is not ideal, but should be safe.
+                    outputVectorMatlab(std::string("reorder_it"),
+                                       dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
+                                       timer.currentStepNum(), output_dir_);
+                }
+            }
+
+            SimulatorReport sreport;
+
+            // Solve pressure equation.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, allcells_, state.saturation(), fractional_flows);
+                wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
+            }
+            bool well_control_passed = !check_well_controls_;
+            int well_control_iteration = 0;
+            do {
+                // Run solver.
+                pressure_timer.start();
+                std::vector<double> initial_pressure = state.pressure();
+                psolver_.solve(timer.currentStepLength(), state, well_state);
+
+                // Renormalize pressure if rock is incompressible, and
+                // there are no pressure conditions (bcs or wells).
+                // It is deemed sufficient for now to renormalize
+                // using geometric volume instead of pore volume.
+                if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
+                    && allNeumannBCs(bcs_) && allRateWells(wells_)) {
+                    // Compute average pressures of previous and last
+                    // step, and total volume.
+                    double av_prev_press = 0.0;
+                    double av_press = 0.0;
+                    double tot_vol = 0.0;
+                    const int num_cells = grid_.number_of_cells;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
+                        av_press      += state.pressure()[cell]*grid_.cell_volumes[cell];
+                        tot_vol       += grid_.cell_volumes[cell];
+                    }
+                    // Renormalization constant
+                    const double ren_const = (av_prev_press - av_press)/tot_vol;
+                    for (int cell = 0; cell < num_cells; ++cell) {
+                        state.pressure()[cell] += ren_const;
+                    }
+                    const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
+                    for (int well = 0; well < num_wells; ++well) {
+                        well_state.bhp()[well] += ren_const;
+                    }
+                }
+
+                // Stop timer and report.
+                pressure_timer.stop();
+                double pt = pressure_timer.secsSinceStart();
+                *log_ << "Pressure solver took:  " << pt << " seconds." << std::endl;
+                ptime += pt;
+                sreport.pressure_time = pt;
+
+                // Optionally, check if well controls are satisfied.
+                if (check_well_controls_) {
+                    Opm::computePhaseFlowRatesPerWell(*wells_,
+                                                      well_state.perfRates(),
+                                                      fractional_flows,
+                                                      well_resflows_phase);
+                    *log_ << "Checking well conditions." << std::endl;
+                    // For testing we set surface := reservoir
+                    well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
+                    ++well_control_iteration;
+                    if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
+                        OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
+                    }
+                    if (!well_control_passed) {
+                        *log_ << "Well controls not passed, solving again." << std::endl;
+                    } else {
+                        *log_ << "Well conditions met." << std::endl;
+                    }
+                }
+            } while (!well_control_passed);
+
+            // Update pore volumes if rock is compressible.
+            if (rock_comp_props_ && rock_comp_props_->isActive()) {
+                initial_porevol = porevol;
+                computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
+            }
+
+            // Process transport sources (to include bdy terms and well flows).
+            Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
+                                        wells_, well_state.perfRates(), transport_src);
+
+            // Solve transport.
+            transport_timer.start();
+            double stepsize = timer.currentStepLength();
+            if (num_transport_substeps_ != 1) {
+                stepsize /= double(num_transport_substeps_);
+                *log_ << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
+            }
+            double injected[2] = { 0.0 };
+            double produced[2] = { 0.0 };
+            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
+                tsolver_->solve(&initial_porevol[0], &transport_src[0], stepsize, state);
+
+                double substep_injected[2] = { 0.0 };
+                double substep_produced[2] = { 0.0 };
+                Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize,
+                                             substep_injected, substep_produced);
+                injected[0] += substep_injected[0];
+                injected[1] += substep_injected[1];
+                produced[0] += substep_produced[0];
+                produced[1] += substep_produced[1];
+                if (use_reorder_ && use_segregation_split_) {
+                    // Again, unfortunate but safe use of dynamic_cast.
+                    // Possible solution: refactor gravity solver to its own class.
+                    dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
+                        .solveGravity(&initial_porevol[0], stepsize, state);
+                }
+                watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
+                              produced[0]/(produced[0] + produced[1]),
+                              tot_produced[0]/tot_porevol_init);
+                if (wells_) {
+                    wellreport.push(props_, *wells_, state.saturation(),
+                                    timer.simulationTimeElapsed() + timer.currentStepLength(),
+                                    well_state.bhp(), well_state.perfRates());
+                }
+            }
+            transport_timer.stop();
+            double tt = transport_timer.secsSinceStart();
+            sreport.transport_time = tt;
+            *log_ << "Transport solver took: " << tt << " seconds." << std::endl;
+            ttime += tt;
+            // Report volume balances.
+            Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
+            tot_injected[0] += injected[0];
+            tot_injected[1] += injected[1];
+            tot_produced[0] += produced[0];
+            tot_produced[1] += produced[1];
+            reportVolumes(*log_, satvol, tot_porevol_init,
+                          tot_injected, tot_produced,
+                          injected, produced,
+                          init_satvol);
+            sreport.total_time =  step_timer.secsSinceStart();
+            if (output_) {
+                sreport.reportParam(tstep_os);
+            }
+
+            // advance the timer to the end of the timestep *before* notifying
+            // the client that the timestep is done
+            ++timer;
+
+            // notify all clients that we are done with the timestep
+            callback_timer.start ();
+            timestep_completed_.signal ();
+            callback_timer.stop ();
+            time_in_callbacks += callback_timer.secsSinceStart ();
+        }
+
+        if (output_) {
+            if (output_vtk_) {
+#ifdef DISABLE_OUTPUT
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+            }
+#ifdef DISABLE_OUTPUT
+            outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+#endif
+            if (use_reorder_) {
+                // This use of dynamic_cast is not ideal, but should be safe.
+                outputVectorMatlab(std::string("reorder_it"),
+                                   dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
+                                   timer.currentStepNum(), output_dir_);
+                }
+            outputWaterCut(watercut, output_dir_);
+            if (wells_) {
+                outputWellReport(wellreport, output_dir_);
+            }
+            tstep_os.close();
+        }
+
+        total_timer.stop();
+
+        SimulatorReport report;
+        report.pressure_time = ptime;
+        report.transport_time = ttime;
+        report.total_time = total_timer.secsSinceStart() - time_in_callbacks;
+        return report;
+    }
+
+
+} // namespace Opm

opm/simulators/SimulatorIncompTwophase.hpp

@@ -0,0 +1,144 @@
+/*
+  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/>.
+*/
+
+#ifndef OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED
+#define OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED
+
+#include <memory>
+#include <vector>
+#include <iostream>
+
+struct UnstructuredGrid;
+struct Wells;
+struct FlowBoundaryConditions;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class IncompPropertiesInterface;
+    class RockCompressibility;
+    class WellsManager;
+    class LinearSolverInterface;
+    class SimulatorTimer;
+    class TwophaseState;
+    class WellState;
+    struct SimulatorReport;
+    struct Event;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorIncompTwophase
+    {
+    public:
+        /// Initialise from parameters and objects to observe.
+        /// \param[in] param       parameters, this class accepts the following:
+        ///     parameter (default)            effect
+        ///     -----------------------------------------------------------
+        ///     output (true)                  write output to files?
+        ///     output_dir ("output")          output directoty
+        ///     output_interval (1)            output every nth step
+        ///     nl_pressure_residual_tolerance (0.0) pressure solver residual tolerance (in Pascal)
+        ///     nl_pressure_change_tolerance (1.0)   pressure solver change tolerance (in Pascal)
+        ///     nl_pressure_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_maxiter (30)                max nonlinear iterations in transport
+        ///     nl_tolerance (1e-9)            transport solver absolute residual tolerance
+        ///     num_transport_substeps (1)     number of transport steps per pressure step
+        ///     use_segregation_split (false)  solve for gravity segregation (if false,
+        ///                                    segregation is ignored).
+        ///
+        /// \param[in] grid          grid data structure
+        /// \param[in] props         fluid and rock properties
+        /// \param[in] rock_comp_props     if non-null, rock compressibility properties
+        /// \param[in] well_manager  well manager, may manage no (null) wells
+        /// \param[in] src           source terms
+        /// \param[in] bcs           boundary conditions, treat as all noflow if null
+        /// \param[in] linsolver     linear solver
+        /// \param[in] gravity       if non-null, gravity vector
+       SimulatorIncompTwophase(const parameter::ParameterGroup& param,
+                               const UnstructuredGrid& grid,
+                               const IncompPropertiesInterface& props,
+                               const RockCompressibility* rock_comp_props,
+                               WellsManager& wells_manager,
+                               const std::vector<double>& src,
+                               const FlowBoundaryConditions* bcs,
+                               LinearSolverInterface& linsolver,
+                               const double* gravity);
+
+        /// Run the simulation.
+        /// This will run succesive timesteps until timer.done() is true. It will
+        /// modify the reservoir and well states.
+        /// \param[in,out] timer       governs the requested reporting timesteps
+        /// \param[in,out] state       state of reservoir: pressure, fluxes
+        /// \param[in,out] well_state  state of wells: bhp, perforation rates
+        /// \return                    simulation report, with timing data
+        SimulatorReport run(SimulatorTimer& timer,
+                            TwophaseState& state,
+                            WellState& well_state);
+
+        /// Event that is signaled every time the simulator has completed a
+        /// a timestep.
+        ///
+        /// Register a callback with this event to do processing at the end
+        /// of every timestep, for instance to do reporting.
+        ///
+        /// \note
+        /// If you want to know the current timestep, the callback must
+        /// also monitor the timer object which was passed to run().
+        ///
+        /// \example
+        /// \code{.cpp}
+        /// struct Foo {
+        ///   void bar () { cout << "Called!" << endl; }
+        /// };
+        ///
+        /// SimulatorIncompTwophase sim (...);
+        /// Foo f;
+        /// sim.timestep_completed ().add <Foo, &Foo::bar> (f);
+        /// sim.run (...);
+        /// \endcode
+        ///
+        /// \note
+        /// Registered callbacks should call the sync() method before
+        /// accessing the state that was passed into the run() method.
+        ///
+        /// \see Opm::SimulatorIncompTwophase::sync
+        Event& timestep_completed ();
+
+        /// Notify the simulator that a callback has an interest in reading
+        /// for reporting purposes the contents of the state argument that
+        /// was passed to the run() method. The simulator will then flush
+        /// any internal state which is currently not reflected in it.
+        ///
+        /// \note
+        /// This should only be called from within a notification which has
+        /// been setup with timestep_completed(). Avoid calling this method
+        /// outside of run().
+        ///
+        /// \see Opm::SimulatorIncompTwophase::run,
+        ///      Opm::SimulatorIncompTwophase::timestep_completed
+        void sync ();
+
+    private:
+        struct Impl;
+        // Using shared_ptr instead of unique_ptr since unique_ptr requires complete type for Impl.
+        std::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED

opm/simulators/thresholdPressures.hpp

@@ -0,0 +1,413 @@
+/*
+  Copyright 2014 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/>.
+*/
+
+#ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
+#define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
+
+#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <vector>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
+#include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
+
+
+namespace Opm
+{
+/// \brief Compute the maximum gravity corrected pressure difference of all
+///        equilibration regions given a reservoir state.
+/// \tparam    Grid           Type of grid object (UnstructuredGrid or CpGrid).
+/// \param[out] maxDp         The resulting pressure difference between equilibration regions
+/// \param[in] deck           Input deck, EQLOPTS and THPRES are accessed from it.
+/// \param[in] eclipseState   Processed eclipse state, EQLNUM is accessed from it.
+/// \param[in] grid           The grid to which the thresholds apply.
+/// \param[in] initialState   The state of the reservoir
+/// \param[in] props          The object which calculates fluid properties
+/// \param[in] gravity        The gravity constant
+template <class Grid>
+void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
+                  const DeckConstPtr& deck,
+                  EclipseStateConstPtr eclipseState,
+                  const Grid& grid,
+                  const BlackoilState& initialState,
+                  const BlackoilPropertiesFromDeck& props,
+                  const double gravity)
+{
+
+    const PhaseUsage& pu = props.phaseUsage();
+
+    std::shared_ptr<const GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
+    const auto& eqlnumData = eqlnum->getData();
+
+    const int numPhases = initialState.numPhases();
+    const int numCells = UgGridHelpers::numCells(grid);
+    const int numPvtRegions = deck->getKeyword("TABDIMS").getRecord(0).getItem("NTPVT").get< int >(0);
+
+    // retrieve the minimum (residual!?) and the maximum saturations for all cells
+    std::vector<double> minSat(numPhases*numCells);
+    std::vector<double> maxSat(numPhases*numCells);
+    std::vector<int> allCells(numCells);
+    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+        allCells[cellIdx] = cellIdx;
+    }
+    props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());
+
+    // retrieve the surface densities
+    std::vector<std::vector<double> > surfaceDensity(numPvtRegions);
+    const auto& densityKw = deck->getKeyword("DENSITY");
+    for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) {
+        surfaceDensity[regionIdx].resize(numPhases);
+
+        if (pu.phase_used[BlackoilPhases::Aqua]) {
+            const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
+            surfaceDensity[regionIdx][wpos] =
+                densityKw.getRecord(regionIdx).getItem("WATER").getSIDouble(0);
+        }
+
+        if (pu.phase_used[BlackoilPhases::Liquid]) {
+            const int opos = pu.phase_pos[BlackoilPhases::Liquid];
+            surfaceDensity[regionIdx][opos] =
+                densityKw.getRecord(regionIdx).getItem("OIL").getSIDouble(0);
+        }
+
+        if (pu.phase_used[BlackoilPhases::Vapour]) {
+            const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+            surfaceDensity[regionIdx][gpos] =
+                densityKw.getRecord(regionIdx).getItem("GAS").getSIDouble(0);
+        }
+    }
+
+    // retrieve the PVT region of each cell. note that we need c++ instead of
+    // Fortran indices.
+    const int* gc = UgGridHelpers::globalCell(grid);
+    std::vector<int> pvtRegion(numCells);
+    const auto& cartPvtRegion = eclipseState->getIntGridProperty("PVTNUM")->getData();
+    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+        const int cartCellIdx = gc ? gc[cellIdx] : cellIdx;
+        pvtRegion[cellIdx] = std::max(0, cartPvtRegion[cartCellIdx] - 1);
+    }
+
+    // compute the initial "phase presence" of each cell (required to calculate
+    // the inverse formation volume factors
+    std::vector<PhasePresence> cond(numCells);
+    for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
+        const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
+        const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
+        const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
+
+        if (pu.phase_used[BlackoilPhases::Aqua] && sw > 0.0) {
+            cond[cellIdx].setFreeWater();
+        }
+
+        if (pu.phase_used[BlackoilPhases::Liquid] && so > 0.0) {
+            cond[cellIdx].setFreeOil();
+        }
+
+        if (pu.phase_used[BlackoilPhases::Vapour] && sg > 0.0) {
+            cond[cellIdx].setFreeGas();
+        }
+    }
+
+    // calculate the initial fluid densities for the gravity correction.
+    std::vector<std::vector<double>> rho(numPhases);
+    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+        rho[phaseIdx].resize(numCells);
+    }
+
+    // compute the capillary pressures of the active phases
+    std::vector<double> capPress(numCells*numPhases);
+    std::vector<int> cellIdxArray(numCells);
+    for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
+        cellIdxArray[cellIdx] = cellIdx;
+    }
+    props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL);
+
+    // compute the absolute pressure of each active phase: for some reason, E100
+    // defines the capillary pressure for the water phase as p_o - p_w while it
+    // uses p_g - p_o for the gas phase. (it would be more consistent to use the
+    // oil pressure as reference for both the other phases.) probably this is
+    // done to always have a positive number for the capillary pressure (as long
+    // as the medium is hydrophilic)
+    std::vector<std::vector<double> > phasePressure(numPhases);
+    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+        phasePressure[phaseIdx].resize(numCells);
+    }
+
+    for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
+        // we currently hard-code the oil phase as the reference phase!
+        assert(pu.phase_used[BlackoilPhases::Liquid]);
+
+        const int opos = pu.phase_pos[BlackoilPhases::Liquid];
+        phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];
+
+        if (pu.phase_used[BlackoilPhases::Aqua]) {
+            const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
+            phasePressure[wpos][cellIdx] =
+                initialState.pressure()[cellIdx]
+                + (capPress[cellIdx*numPhases + opos] - capPress[cellIdx*numPhases + wpos]);
+        }
+
+        if (pu.phase_used[BlackoilPhases::Vapour]) {
+            const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+            phasePressure[gpos][cellIdx] =
+                initialState.pressure()[cellIdx]
+                + (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]);
+        }
+    }
+
+    // calculate the densities of the active phases for each cell
+    if (pu.phase_used[BlackoilPhases::Aqua]) {
+        const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
+        const auto& pvtw = props.waterPvt();
+        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
+            int pvtRegionIdx = pvtRegion[cellIdx];
+
+            double T = initialState.temperature()[cellIdx];
+            double p = initialState.pressure()[cellIdx];
+            double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+
+            rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
+        }
+    }
+
+    if (pu.phase_used[BlackoilPhases::Liquid]) {
+        const int opos = pu.phase_pos[BlackoilPhases::Liquid];
+        const auto& pvto = props.oilPvt();
+        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
+            int pvtRegionIdx = pvtRegion[cellIdx];
+
+            double T = initialState.temperature()[cellIdx];
+            double p = initialState.pressure()[cellIdx];
+            double Rs = initialState.gasoilratio()[cellIdx];
+            double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p);
+
+            if (Rs >= RsSat) {
+                double b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
+                rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
+            }
+            else {
+                double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
+                rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
+                if (pu.phase_used[BlackoilPhases::Vapour]) {
+                    int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+                    rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
+                }
+            }
+        }
+    }
+
+    if (pu.phase_used[BlackoilPhases::Vapour]) {
+        const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
+        const auto& pvtg = props.gasPvt();
+        for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
+            int pvtRegionIdx = pvtRegion[cellIdx];
+
+            double T = initialState.temperature()[cellIdx];
+            double p = initialState.pressure()[cellIdx];
+            double Rv = initialState.rv()[cellIdx];
+            double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p);
+
+            if (Rv >= RvSat) {
+                double b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
+                rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
+            }
+            else {
+                double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
+                rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
+
+                if (pu.phase_used[BlackoilPhases::Liquid]) {
+                    int opos = pu.phase_pos[BlackoilPhases::Liquid];
+                    rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
+                }
+            }
+        }
+    }
+
+    // Calculate the maximum pressure potential difference between all PVT region
+    // transitions of the initial solution.
+    const int num_faces = UgGridHelpers::numFaces(grid);
+    const auto& fc = UgGridHelpers::faceCells(grid);
+    for (int face = 0; face < num_faces; ++face) {
+        const int c1 = fc(face, 0);
+        const int c2 = fc(face, 1);
+        if (c1 < 0 || c2 < 0) {
+            // Boundary face, skip this.
+            continue;
+        }
+        const int gc1 = (gc == 0) ? c1 : gc[c1];
+        const int gc2 = (gc == 0) ? c2 : gc[c2];
+        const int eq1 = eqlnumData[gc1];
+        const int eq2 = eqlnumData[gc2];
+
+        if (eq1 == eq2) {
+            // not an equilibration region boundary. skip this.
+            continue;
+        }
+
+        // update the maximum pressure potential difference between the two
+        // regions
+        const auto barrierId = std::make_pair(eq1, eq2);
+        if (maxDp.count(barrierId) == 0) {
+            maxDp[barrierId] = 0.0;
+        }
+
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
+            const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
+            const double zAvg = (z1 + z2)/2; // average depth
+
+            const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
+
+            const double s1 = initialState.saturation()[numPhases*c1 + phaseIdx];
+            const double s2 = initialState.saturation()[numPhases*c2 + phaseIdx];
+
+            const double sResid1 = minSat[numPhases*c1 + phaseIdx];
+            const double sResid2 = minSat[numPhases*c2 + phaseIdx];
+
+            // compute gravity corrected pressure potentials at the average depth
+            const double p1 = phasePressure[phaseIdx][c1] + rhoAvg*gravity*(zAvg - z1);
+            const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(zAvg - z2);
+
+            if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
+                maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
+        }
+    }
+}
+
+    /// \brief Get a vector of pressure thresholds from EclipseState.
+    /// This function looks at EQLOPTS, THPRES and EQLNUM to determine
+    /// pressure thresholds.  It does not consider the case where the
+    /// threshold values are defaulted, in which case they should be
+    /// determined from the initial, equilibrated simulation state.
+    /// \tparam    Grid           Type of grid object (UnstructuredGrid or CpGrid).
+    /// \param[in] deck           Input deck, EQLOPTS and THPRES are accessed from it.
+    /// \param[in] eclipseState   Processed eclipse state, EQLNUM is accessed from it.
+    /// \param[in] maxDp          The maximum gravity corrected pressure differences between
+    ///                           the equilibration regions.
+    /// \param[in] grid           The grid to which the thresholds apply.
+    /// \return                   A vector of pressure thresholds, one
+    ///                           for each face in the grid. A value
+    ///                           of zero means no threshold for that
+    ///                           particular face. An empty vector is
+    ///                           returned if there is no THPRES
+    ///                           feature used in the deck.
+
+
+
+    template <class Grid>
+    std::vector<double> thresholdPressures(const DeckConstPtr& /* deck */,
+                                           EclipseStateConstPtr eclipseState,
+                                           const Grid& grid,
+                                           const std::map<std::pair<int, int>, double>& maxDp)
+    {
+        SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig();
+        std::vector<double> thpres_vals;
+        if (simulationConfig->hasThresholdPressure()) {
+            std::shared_ptr<const ThresholdPressure> thresholdPressure = simulationConfig->getThresholdPressure();
+            std::shared_ptr<const GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
+            const auto& eqlnumData = eqlnum->getData();
+
+            // Set threshold pressure values for each cell face.
+            const int num_faces = UgGridHelpers::numFaces(grid);
+            const auto& fc = UgGridHelpers::faceCells(grid);
+            const int* gc = UgGridHelpers::globalCell(grid);
+            thpres_vals.resize(num_faces, 0.0);
+            for (int face = 0; face < num_faces; ++face) {
+                const int c1 = fc(face, 0);
+                const int c2 = fc(face, 1);
+                if (c1 < 0 || c2 < 0) {
+                    // Boundary face, skip it.
+                    continue;
+                }
+                const int gc1 = (gc == 0) ? c1 : gc[c1];
+                const int gc2 = (gc == 0) ? c2 : gc[c2];
+                const int eq1 = eqlnumData[gc1];
+                const int eq2 = eqlnumData[gc2];
+
+                if (thresholdPressure->hasRegionBarrier(eq1,eq2)) {
+                    if (thresholdPressure->hasThresholdPressure(eq1,eq2)) {
+                        thpres_vals[face] = thresholdPressure->getThresholdPressure(eq1,eq2);
+                    }
+                    else {
+                        // set the threshold pressure for faces of PVT regions where the third item
+                        // has been defaulted to the maximum pressure potential difference between
+                        // these regions
+                        const auto barrierId = std::make_pair(eq1, eq2);
+                        thpres_vals[face] = maxDp.at(barrierId);
+                    }
+                }
+
+            }
+        }
+        return thpres_vals;
+    }
+
+    /// \brief Get a vector of pressure thresholds from either EclipseState
+    /// or maxDp (for defaulted values) for all Non-neighbour connections (NNCs).
+    /// \param[in] nnc            The NNCs,
+    /// \param[in] eclipseState   Processed eclipse state, EQLNUM is accessed from it.
+    /// \param[in] maxDp          The maximum gravity corrected pressure differences between
+    ///                           the equilibration regions.
+    /// \return                   A vector of pressure thresholds, one
+    ///                           for each NNC in the grid. A value
+    ///                           of zero means no threshold for that
+    ///                           particular connection. An empty vector is
+    ///                           returned if there is no THPRES
+    ///                           feature used in the deck.
+     std::vector<double> thresholdPressuresNNC(EclipseStateConstPtr eclipseState,
+                                               const NNC& nnc,
+                                               const std::map<std::pair<int, int>, double>& maxDp)
+    {
+        SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig();
+        std::vector<double> thpres_vals;
+        if (simulationConfig->hasThresholdPressure()) {
+            std::shared_ptr<const ThresholdPressure> thresholdPressure = simulationConfig->getThresholdPressure();
+            std::shared_ptr<const GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
+            auto eqlnumData = eqlnum->getData();
+
+            // Set values for each NNC
+
+            thpres_vals.resize(nnc.numNNC(), 0.0);
+            for (size_t i = 0 ; i < nnc.numNNC(); ++i) {
+                const int gc1 = nnc.nncdata()[i].cell1;
+                const int gc2 = nnc.nncdata()[i].cell2;
+                const int eq1 = eqlnumData[gc1];
+                const int eq2 = eqlnumData[gc2];
+
+                if (thresholdPressure->hasRegionBarrier(eq1,eq2)) {
+                    if (thresholdPressure->hasThresholdPressure(eq1,eq2)) {
+                        thpres_vals[i] = thresholdPressure->getThresholdPressure(eq1,eq2);
+                    } else {
+                        // set the threshold pressure for NNC of PVT regions where the third item
+                        // has been defaulted to the maximum pressure potential difference between
+                        // these regions
+                        const auto barrierId = std::make_pair(eq1, eq2);
+                        thpres_vals[i] = maxDp.at(barrierId);
+                    }
+                }
+            }
+        }
+        return thpres_vals;
+    }
+}
+
+#endif // OPM_THRESHOLDPRESSURES_HEADER_INCLUDED

tests/test_thresholdpressure.cpp

@@ -0,0 +1,56 @@
+#include <opm/parser/eclipse/Parser/Parser.hpp>
+#include <opm/parser/eclipse/Parser/ParseContext.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/grid/cornerpoint_grid.h>
+#include <opm/core/grid/GridManager.hpp>
+#include <opm/autodiff/GridHelpers.hpp>
+#include <opm/simulators/thresholdPressures.hpp> // Note: the GridHelpers must be included before this (to make overloads available)
+
+
+static DeckPtr createDeckSimConfig() {
+    const std::string& inputStr = "RUNSPEC\n"
+                                  "EQLOPTS\n"
+                                  "THPRES /\n "
+                                  "DIMENS\n"
+                                  "10 3 4 /\n"
+                                  "\n"
+                                  "GRID\n"
+                                  "REGIONS\n"
+                                  "EQLNUM\n"
+                                  "10*1 10*2 100*3 /\n "
+                                  "\n"
+
+                                  "SOLUTION\n"
+                                  "THPRES\n"
+                                  "1 2 12.0/\n"
+                                  "1 3 5.0/\n"
+                                  "2 3 7.0/\n"
+                                  "/\n"
+                                  "\n";
+
+
+    ParserPtr parser(new Parser());
+    return parser->parseString(inputStr, ParseContext()) ;
+}
+
+/*
+  The main purpose of this is to get the code in
+  thresholdPressures.hpp to compile; it currently has no users in the
+  opm-core codebase.
+*/
+
+BOOST_AUTO_TEST_CASE(CreateSimulationConfig) {
+    ParseContext parseContext;
+    typedef UnstructuredGrid  Grid;
+    DeckPtr deck = createDeckSimConfig();
+    EclipseState state(deck, parseContext);
+    EclipseGridConstPtr eclipseGrid = state.getEclipseGrid();
+    std::vector<double> porv = eclipseState->getDoubleGridProperty("PORV")->getData();
+    GridManager gridManager( eclipseState->getEclipseGrid(), porv );
+    const Grid& grid = *(gridManager.c_grid());
+
+    std::vector<double> threshold_pressures = thresholdPressures(parseContext, eclipseState, grid);
+    BOOST_CHECK( threshold_pressures.size() > 0 );
+}