Merge pull request #22 from atgeirr/master

Multiple changes to simulators and solvers

Makefile.am

@@ -92,6 +92,7 @@ opm/core/pressure/tpfa/compr_source.c			\
 opm/core/pressure/tpfa/ifs_tpfa.c			\
 opm/core/pressure/tpfa/trans_tpfa.c			\
 opm/core/pressure/well.c				\
+opm/core/simulator/SimulatorCompressibleTwophase.cpp	\
 opm/core/simulator/SimulatorIncompTwophase.cpp		\
 opm/core/simulator/SimulatorReport.cpp			\
 opm/core/simulator/SimulatorTimer.cpp			\
@@ -195,6 +196,7 @@ opm/core/pressure/tpfa/compr_source.h			\
 opm/core/pressure/tpfa/ifs_tpfa.h			\
 opm/core/pressure/tpfa/trans_tpfa.h			\
 opm/core/simulator/BlackoilState.hpp			\
+opm/core/simulator/SimulatorCompressibleTwophase.hpp	\
 opm/core/simulator/SimulatorReport.hpp			\
 opm/core/simulator/SimulatorIncompTwophase.hpp		\
 opm/core/simulator/SimulatorTimer.hpp			\

README

@@ -1,28 +1,42 @@
-These are the release notes for opm-core
-
 Open Porous Media Core Library
+==============================
+
+These are release notes for opm-core.
+
 
 CONTENT
+-------
+
 opm-core is the core library within OPM and contains the following 
 
-* Eclipse preprosessing
+* Eclipse deck input and preprosessing
 * Fluid properties (basic PVT models and rock properties)
-* Grid (generates different grids)
+* Grid handling (cornerpoint grids, unstructured grid interface)
-* Linear Algerbra (interface to different linear solvers)
+* Linear Algebra (interface to different linear solvers)
-* Pressure solvers (different discretization schemes, different flow models)
+* Pressure solvers (various discretization schemes, flow models)
-* Simulator (some basic examples of simulators based on sequential schemes)
+* Simulators (some basic examples of simulators based on sequential splitting schemes)
-* Transport solvers (different discretization schemes, different flow models)
+* Transport solvers (various discretization schemes, flow models)
-* Utility (input and output processing, unit conversion)
+* Utilities (input and output processing, unit conversion)
 * Wells (basic well handling)
 
-ON WHAT PLATFORMS DOES IT RUN?
-The opm-core module is designed to run on linux platforms. No efforts have
-been made to ensure that the code will compile and run on windows platforms.
 
-DOCUMENTATION
+LICENSE
-Efforts have been made to document the code with Doxygen.
+-------
+
+The library is distributed under the GNU General Public License,
+version 3 or later (GPLv3+).
+
+
+PLATFORMS
+---------
+
+The opm-core module is designed to run on Linux platforms. It is also
+regularly run on Mac OS X. No efforts have been made to ensure that
+the code will compile and run on windows platforms.
+
 
 DEPENDENCIES FOR DEBIAN BASED DISTRIBUTIONS (Debian Squeeze/Ubuntu Precise)
+---------------------------------------------------------------------------
 
 # packages necessary for building
 sudo apt-get install -y build-essential gfortran pkg-config libtool \
@@ -40,7 +54,9 @@ sudo apt-get install -y libboost-all-dev libsuperlu3-dev libsuitesparse-dev
 # libraries necessary for OPM
 sudo apt-get install -y libxml0-dev
 
+
 DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS
+-----------------------------------------
 
 # libraries
 sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack
@@ -48,7 +64,11 @@ sudo zypper install blas libblas3 lapack liblapack3 libboost libxml2 umfpack
 # tools
 sudo zypper install gcc automake autoconf git doxygen
 
+
 RETRIEVING AND BUILDING DUNE PREREQUISITES
+------------------------------------------
+
+(only necessary if you want to use opm-core as a dune module)
 
 # trust DUNE certificate (sic)
 echo p | svn list https://svn.dune-project.org/svn/dune-common
@@ -67,20 +87,39 @@ for module in common istl geometry grid localfunctions; do
 	--make-opts="-j -l 0.8" autogen : configure : make
 done
 
-DOWNLOAD opm-core
 
+DOWNLOADING
+-----------
+
+For a read-only download:
 git clone git://github.com/OPM/opm-core.git
 
-BUILDING opm-core
+If you want to contribute, fork OPM/opm-core on github.
 
-# note: this is done from the parent directory of opm-core
-env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \
-  --configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make
 
-or, without using dunecontrol:
+BUILDING
+--------
 
-pushd opm-core
+(standalone opm-core:)
-autoreconf -i
+
-./configure
+ cd ../opm-core
-make
+ autoreconf -i
-sudo make install
+ ./configure
+ make
+ sudo make install
+
+(using opm-core as a dune module:)
+
+ # note: this is done from the parent directory of opm-core
+ env CCACHE_DISABLE=1 dune-common/bin/dunecontrol --only=opm-core \
+   --configure-opts="" --make-opts="-j -l 0.8" autogen : configure : make
+
+
+
+DOCUMENTATION
+-------------
+
+Efforts have been made to document the code with Doxygen.
+In order to build the documentation, enter the command
+$ doxygen
+in the topmost directory.

examples/Makefile.am

@@ -17,6 +17,7 @@ $(BOOST_SYSTEM_LIB)
 noinst_PROGRAMS =				\
 refine_wells					\
 scaneclipsedeck					\
+sim_2p_comp_reorder				\
 sim_2p_incomp_reorder				\
 sim_wateroil					\
 wells_example
@@ -29,6 +30,7 @@ wells_example
 # Please maintain sort order from "noinst_PROGRAMS".
 
 refine_wells_SOURCES          = refine_wells.cpp
+sim_2p_comp_reorder_SOURCES   = sim_2p_comp_reorder.cpp
 sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp
 sim_wateroil_SOURCES          = sim_wateroil.cpp
 wells_example_SOURCES         = wells_example.cpp

examples/sim_2p_comp_reorder.cpp

@@ -0,0 +1,285 @@
+/*
+  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/GridManager.hpp>
+#include <opm/core/newwells.h>
+#include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/utility/initState.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>
+#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/fluid/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 <boost/scoped_ptr.hpp>
+#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)
+{
+    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");
+    boost::scoped_ptr<EclipseGridParser> deck;
+    boost::scoped_ptr<GridManager> grid;
+    boost::scoped_ptr<BlackoilPropertiesInterface> props;
+    boost::scoped_ptr<RockCompressibility> rock_comp;
+    BlackoilState state;
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    if (use_deck) {
+        std::string deck_filename = param.get<std::string>("deck_filename");
+        deck.reset(new EclipseGridParser(deck_filename));
+        // Grid init
+        grid.reset(new GridManager(*deck));
+        // Rock and fluid init
+        props.reset(new BlackoilPropertiesFromDeck(*deck, *grid->c_grid()));
+        // check_well_controls = param.getDefault("check_well_controls", false);
+        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+        // Rock compressibility.
+        rock_comp.reset(new RockCompressibility(*deck));
+        // Gravity.
+        gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
+        // Init state variables (saturation and pressure).
+        if (param.has("init_saturation")) {
+            initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        } else {
+            initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
+        }
+        initBlackoilSurfvol(*grid->c_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));
+        // Rock and fluid init.
+        props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        // Rock compressibility.
+        rock_comp.reset(new RockCompressibility(param));
+        // Gravity.
+        gravity[2] = param.getDefault("gravity", 0.0);
+        // Init state variables (saturation and pressure).
+        initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
+        initBlackoilSurfvol(*grid->c_grid(), *props, state);
+    }
+
+    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 (...) {
+            THROW("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"
+              << "                        (number of epochs: "
+              << (use_deck ? deck->numberOfEpochs() : 1) << ")\n\n" << std::flush;
+
+    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.
+        deck->setCurrentEpoch(deck->numberOfEpochs() - 1);
+        simtimer.init(*deck);
+        const double total_time = simtimer.totalTime();
+        for (int epoch = 0; epoch < deck->numberOfEpochs(); ++epoch) {
+            // Set epoch index.
+            deck->setCurrentEpoch(epoch);
+
+            // Update the timer.
+            if (deck->hasField("TSTEP")) {
+                simtimer.init(*deck);
+            } else {
+                if (epoch != 0) {
+                    THROW("No TSTEP in deck for epoch " << epoch);
+                }
+                simtimer.init(param);
+            }
+            simtimer.setCurrentStepNum(step);
+            simtimer.setTotalTime(total_time);
+
+            // Report on start of epoch.
+            std::cout << "\n\n--------------    Starting epoch " << epoch << "    --------------"
+                      << "\n                  (number of steps: "
+                      << simtimer.numSteps() - step << ")\n\n" << std::flush;
+
+            // Create new wells, well_state
+            WellsManager wells(*deck, *grid->c_grid(), props->permeability());
+            // @@@ HACK: we should really make a new well state and
+            // properly transfer old well state to it every epoch,
+            // since number of wells may change etc.
+            if (epoch == 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 (epoch == 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);
+    }
+
+}

examples/sim_wateroil.cpp

@@ -141,7 +141,6 @@ main(int argc, char** argv)
     std::cout << "---------------    Reading parameters     ---------------" << std::endl;
 
     // Reading various control parameters.
-    const bool use_reorder = param.getDefault("use_reorder", true);
     const bool output = param.getDefault("output", true);
     std::string output_dir;
     int output_interval = 1;
@@ -229,28 +228,8 @@ main(int argc, char** argv)
         }
     }
     bool use_segregation_split = false;
-    bool use_column_solver = false;
+    if (use_gravity) {
-    bool use_gauss_seidel_gravity = false;
-    if (use_gravity && use_reorder) {
         use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
-        if (use_segregation_split) {
-            use_column_solver = param.getDefault("use_column_solver", use_column_solver);
-            if (use_column_solver) {
-                use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity);
-            }
-        }
-    }
-
-    // Check that rock compressibility is not used with solvers that do not handle it.
-    // int nl_pressure_maxiter = 0;
-    // double nl_pressure_tolerance = 0.0;
-    if (rock_comp->isActive()) {
-        THROW("No rock compressibility in comp. pressure solver yet.");
-        if (!use_reorder) {
-            THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
-        }
-        // nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
-        // nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
     }
 
     // Source-related variables init.
@@ -262,11 +241,11 @@ main(int argc, char** argv)
     // Extra rock init.
     std::vector<double> porevol;
     if (rock_comp->isActive()) {
-        THROW("CompressibleTpfa solver does not handle this.");
         computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
     } else {
         computePorevolume(*grid->c_grid(), props->porosity(), porevol);
     }
+    std::vector<double> initial_porevol = porevol;
     double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
 
 
@@ -293,20 +272,20 @@ main(int argc, char** argv)
     const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0);
     const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20);
     const double *grav = use_gravity ? &gravity[0] : 0;
-    Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, linsolver,
+    Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
                                   nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
                                   grav, wells->c_wells());
     // Reordering solver.
     const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
     const int nl_maxiter = param.getDefault("nl_maxiter", 30);
     Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
-    if (use_gauss_seidel_gravity) {
+    if (use_segregation_split) {
-        reorder_model.initGravity();
+        reorder_model.initGravity(grav);
     }
 
     // Column-based gravity segregation solver.
     std::vector<std::vector<int> > columns;
-    if (use_column_solver) {
+    if (use_segregation_split) {
         Opm::extractColumn(*grid->c_grid(), columns);
     }
 
@@ -409,6 +388,11 @@ main(int argc, char** argv)
         Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
                                     wells->c_wells(), well_state.perfRates(), reorder_src);
 
+        // Compute new porevolumes after pressure solve, if necessary.
+        if (rock_comp->isActive()) {
+            initial_porevol = porevol;
+            computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
+        }
         // Solve transport.
         transport_timer.start();
         double stepsize = simtimer.currentStepLength();
@@ -419,11 +403,13 @@ main(int argc, char** argv)
         for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
             // Note that for now we do not handle rock compressibility,
             // although the transport solver should be able to.
-            reorder_model.solve(&state.faceflux()[0], &state.pressure()[0], &state.surfacevol()[0],
+            reorder_model.solve(&state.faceflux()[0], &state.pressure()[0],
-                                &porevol[0],  &porevol[0], &reorder_src[0], stepsize, state.saturation());
+                                &porevol[0],  &initial_porevol[0], &reorder_src[0], stepsize, 
+                                state.saturation(), state.surfacevol());
             // Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
             if (use_segregation_split) {
-                reorder_model.solveGravity(columns, &state.pressure()[0], &porevol[0], stepsize, grav, state.saturation());
+                reorder_model.solveGravity(columns, &state.pressure()[0], &initial_porevol[0],
+                                           stepsize, state.saturation(), state.surfacevol());
             }
         }
         transport_timer.stop();

generate_doc_figures.py

@@ -1,6 +1,24 @@
+
+# This script generate the illustration pictures for the documentation.
+#
+# To run this script, you have to install paraview, see:
+#
+# http://www.paraview.org/paraview/resources/software.php
+#
+# Eventually, set up the paths (figure_path, tutorial_data_path, tutorial_path) according to your own installation.
+# (The default values should be ok.)
+#
+# Make sure that pvpython is in your path of executables.
+#
+# Run the following command at the root of the directory where
+# opm-core is installed (which also corresponds to the directory where
+# generate_doc_figures is located):
+#
+#   pvpython generate_doc_figures.py
+#
+
 from subprocess import call
 from paraview.simple import *
-# from paraview import servermanager
 from os import remove, mkdir, curdir
 from os.path import join, isdir
 
@@ -13,7 +31,6 @@ collected_garbage_file = []
 if not isdir(figure_path):
     mkdir(figure_path)
     
-# connection = servermanager.Connect()
 
 # tutorial 1
 call(join(tutorial_path, "tutorial1"))

opm/core/GridManager.cpp

@@ -166,6 +166,16 @@ namespace Opm
     // Construct tensor grid from deck.
     void GridManager::initFromDeckTensorgrid(const Opm::EclipseGridParser& deck)
     {
+        // Extract logical cartesian size.
+        std::vector<int> dims;
+        if (deck.hasField("DIMENS")) {
+            dims = deck.getIntegerValue("DIMENS");
+        } else if (deck.hasField("SPECGRID")) {
+            dims = deck.getSPECGRID().dimensions;
+        } else {
+            THROW("Deck must have either DIMENS or SPECGRID.");
+        }
+
         // Extract coordinates (or offsets from top, in case of z).
         const std::vector<double>& dxv = deck.getFloatingPointValue("DXV");
         const std::vector<double>& dyv = deck.getFloatingPointValue("DYV");
@@ -174,6 +184,17 @@ namespace Opm
         std::vector<double> y = coordsFromDeltas(dyv);
         std::vector<double> z = coordsFromDeltas(dzv);
 
+        // Check that number of cells given are consistent with DIMENS/SPECGRID.
+        if (dims[0] != int(dxv.size())) {
+            THROW("Number of DXV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[1] != int(dyv.size())) {
+            THROW("Number of DYV data points do not match DIMENS or SPECGRID.");
+        }
+        if (dims[2] != int(dzv.size())) {
+            THROW("Number of DZV data points do not match DIMENS or SPECGRID.");
+        }
+
         // Extract top corner depths, if available.
         const double* top_depths = 0;
         std::vector<double> top_depths_vec;

opm/core/fluid/RockFromDeck.cpp

@@ -36,8 +36,6 @@ namespace Opm
         PermeabilityKind fillTensor(const EclipseGridParser&                 parser,
                                     std::vector<const std::vector<double>*>& tensor,
                                     std::tr1::array<int,9>&                     kmap);
-
-        int numGlobalCells(const EclipseGridParser& parser);
     } // anonymous namespace
 
 
@@ -87,7 +85,7 @@ namespace Opm
                                           double perm_threshold)
     {
         const int dim              = 3;
-        const int num_global_cells = numGlobalCells(parser);
+        const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2];
         const int nc = grid.number_of_cells;
 
         ASSERT (num_global_cells > 0);
@@ -334,26 +332,6 @@ namespace Opm
             return kind;
         }
 
-        int numGlobalCells(const EclipseGridParser& parser)
-        {
-            int ngc = -1;
-
-            if (parser.hasField("DIMENS")) {
-                const std::vector<int>&
-                    dims = parser.getIntegerValue("DIMENS");
-
-                ngc = dims[0] * dims[1] * dims[2];
-            }
-            else if (parser.hasField("SPECGRID")) {
-                const SPECGRID& sgr = parser.getSPECGRID();
-
-                ngc  = sgr.dimensions[ 0 ];
-                ngc *= sgr.dimensions[ 1 ];
-                ngc *= sgr.dimensions[ 2 ];
-            }
-
-            return ngc;
-        }
     } // anonymous namespace
 
 } // namespace Opm

opm/core/pressure/CompressibleTpfa.cpp

@@ -30,6 +30,7 @@
 #include <opm/core/newwells.h>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/fluid/RockCompressibility.hpp>
 
 #include <algorithm>
 #include <cmath>
@@ -58,6 +59,7 @@ namespace Opm
     ///                                to change.
     CompressibleTpfa::CompressibleTpfa(const UnstructuredGrid& grid,
                                        const BlackoilPropertiesInterface& props,
+                                       const RockCompressibility* rock_comp_props,
                                        const LinearSolverInterface& linsolver,
                                        const double residual_tol,
                                        const double change_tol,
@@ -66,6 +68,7 @@ namespace Opm
                                        const struct Wells* wells)
         : grid_(grid),
           props_(props),
+          rock_comp_props_(rock_comp_props),
           linsolver_(linsolver),
           residual_tol_(residual_tol),
           change_tol_(change_tol),
@@ -74,8 +77,8 @@ namespace Opm
           wells_(wells),
           htrans_(grid.cell_facepos[ grid.number_of_cells ]),
           trans_ (grid.number_of_faces),
-          porevol_(grid.number_of_cells),
+          allcells_(grid.number_of_cells),
-          allcells_(grid.number_of_cells)
+          singular_(false)
     {
         if (wells_ && (wells_->number_of_phases != props.numPhases())) {
             THROW("Inconsistent number of phases specified (wells vs. props): "
@@ -86,7 +89,12 @@ namespace Opm
         UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
         tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
         tpfa_trans_compute(gg, &htrans_[0], &trans_[0]);
-        computePorevolume(grid_, props.porosity(), porevol_);
+        // If we have rock compressibility, pore volumes are updated
+        // in the compute*() methods, otherwise they are constant and
+        // hence may be computed here.
+        if (rock_comp_props_ == NULL || !rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props.porosity(), porevol_);
+        }
         for (int c = 0; c < grid.number_of_cells; ++c) {
             allcells_[c] = c;
         }
@@ -182,6 +190,21 @@ namespace Opm
 
 
 
+    /// @brief After solve(), was the resulting pressure singular.
+    /// Returns true if the pressure is singular in the following
+    /// sense: if everything is incompressible and there are no
+    /// pressure conditions, the absolute values of the pressure
+    /// solution are arbitrary. (But the differences in pressure
+    /// are significant.)
+    bool CompressibleTpfa::singularPressure() const
+    {
+        return singular_;
+    }
+
+
+
+
+
     /// Compute well potentials.
     void CompressibleTpfa::computeWellPotentials(const BlackoilState& state)
     {
@@ -230,6 +253,9 @@ namespace Opm
                                                       const WellState& /*well_state*/)
     {
         computeWellPotentials(state);
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), initial_porevol_);
+        }
     }
 
 
@@ -252,6 +278,8 @@ namespace Opm
         // std::vector<double> face_gravcap_;
         // std::vector<double> wellperf_A_;
         // std::vector<double> wellperf_phasemob_;
+        // std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        // std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
         computeCellDynamicData(dt, state, well_state);
         computeFaceDynamicData(dt, state, well_state);
         computeWellDynamicData(dt, state, well_state);
@@ -273,6 +301,8 @@ namespace Opm
         // std::vector<double> cell_viscosity_;
         // std::vector<double> cell_phasemob_;
         // std::vector<double> cell_voldisc_;
+        // std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        // std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
         const int nc = grid_.number_of_cells;
         const int np = props_.numPhases();
         const double* cell_p = &state.pressure()[0];
@@ -296,6 +326,14 @@ namespace Opm
         // TODO: Check this!
         cell_voldisc_.clear();
         cell_voldisc_.resize(nc, 0.0);
+
+        if (rock_comp_props_ && rock_comp_props_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol_);
+            rock_comp_.resize(nc);
+            for (int cell = 0; cell < nc; ++cell) {
+                rock_comp_[cell] = rock_comp_props_->rockComp(state.pressure()[cell]);
+            }
+        }
     }
 
 
@@ -465,9 +503,20 @@ namespace Opm
         cq.Af = &face_A_[0];
         cq.phasemobf = &face_phasemob_[0];
         cq.voldiscr = &cell_voldisc_[0];
-        cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
+        int was_adjusted = 0;
-                              &face_gravcap_[0], cell_press, well_bhp,
+        if (! (rock_comp_props_ && rock_comp_props_->isActive())) {
-                              &porevol_[0], h_);
+            was_adjusted =
+                cfs_tpfa_res_assemble(gg, dt, &forces, z, &cq, &trans_[0],
+                                      &face_gravcap_[0], cell_press, well_bhp,
+                                      &porevol_[0], h_);
+        } else {
+            was_adjusted =
+                cfs_tpfa_res_comprock_assemble(gg, dt, &forces, z, &cq, &trans_[0],
+                                               &face_gravcap_[0], cell_press, well_bhp,
+                                               &porevol_[0], &initial_porevol_[0],
+                                               &rock_comp_[0], h_);
+        }
+        singular_ = (was_adjusted == 1);
     }
 
 

opm/core/pressure/CompressibleTpfa.hpp

@@ -33,6 +33,7 @@ namespace Opm
 
     class BlackoilState;
     class BlackoilPropertiesInterface;
+    class RockCompressibility;
     class LinearSolverInterface;
     class WellState;
 
@@ -44,23 +45,25 @@ namespace Opm
     {
     public:
         /// Construct solver.
-        /// \param[in] grid          A 2d or 3d grid.
+        /// \param[in] grid             A 2d or 3d grid.
-        /// \param[in] props         Rock and fluid properties.
+        /// \param[in] props            Rock and fluid properties.
-        /// \param[in] linsolver     Linear solver to use.
+        /// \param[in] rock_comp_props  Rock compressibility properties. May be null.
-        /// \param[in] residual_tol  Solution accepted if inf-norm of residual is smaller.
+        /// \param[in] linsolver        Linear solver to use.
-        /// \param[in] change_tol    Solution accepted if inf-norm of change in pressure is smaller.
+        /// \param[in] residual_tol     Solution accepted if inf-norm of residual is smaller.
-        /// \param[in] maxiter       Maximum acceptable number of iterations.
+        /// \param[in] change_tol       Solution accepted if inf-norm of change in pressure is smaller.
-        /// \param[in] gravity       Gravity vector. If non-null, the array should
+        /// \param[in] maxiter          Maximum acceptable number of iterations.
-        ///                          have D elements.
+        /// \param[in] gravity          Gravity vector. If non-null, the array should
-        /// \param[in] wells         The wells argument. Will be used in solution,
+        ///                             have D elements.
-        ///                          is ignored if NULL.
+        /// \param[in] wells            The wells argument. Will be used in solution,
-        ///                          Note: this class observes the well object, and
+        ///                             is ignored if NULL.
-        ///                                makes the assumption that the well topology
+        ///                             Note: this class observes the well object, and
-        ///                                and completions does not change during the
+        ///                                   makes the assumption that the well topology
-        ///                                run. However, controls (only) are allowed
+        ///                                   and completions does not change during the
-        ///                                to change.
+        ///                                   run. However, controls (only) are allowed
-	CompressibleTpfa(const UnstructuredGrid& grid,
+        ///                                   to change.
+        CompressibleTpfa(const UnstructuredGrid& grid,
                          const BlackoilPropertiesInterface& props,
+                         const RockCompressibility* rock_comp_props,
                          const LinearSolverInterface& linsolver,
                          const double residual_tol,
                          const double change_tol,
@@ -68,8 +71,8 @@ namespace Opm
                          const double* gravity,
                          const Wells* wells);
 
-	/// Destructor.
+        /// Destructor.
-	~CompressibleTpfa();
+        ~CompressibleTpfa();
 
         /// Solve the pressure equation by Newton-Raphson scheme.
         /// May throw an exception if the number of iterations
@@ -78,6 +81,14 @@ namespace Opm
                    BlackoilState& state,
                    WellState& well_state);
 
+        /// @brief After solve(), was the resulting pressure singular.
+        /// Returns true if the pressure is singular in the following
+        /// sense: if everything is incompressible and there are no
+        /// pressure conditions, the absolute values of the pressure
+        /// solution are arbitrary. (But the differences in pressure
+        /// are significant.)
+        bool singularPressure() const;
+
     private:
         void computePerSolveDynamicData(const double dt,
                                         const BlackoilState& state,
@@ -101,28 +112,29 @@ namespace Opm
         void solveIncrement();
         double residualNorm() const;
         double incrementNorm() const;
-	void computeResults(BlackoilState& state,
+        void computeResults(BlackoilState& state,
                             WellState& well_state) const;
 
         // ------ Data that will remain unmodified after construction. ------
-	const UnstructuredGrid& grid_;
+        const UnstructuredGrid& grid_;
         const BlackoilPropertiesInterface& props_;
+        const RockCompressibility* rock_comp_props_;
         const LinearSolverInterface& linsolver_;
         const double residual_tol_;
         const double change_tol_;
         const int maxiter_;
         const double* gravity_; // May be NULL
         const Wells* wells_;    // May be NULL, outside may modify controls (only) between calls to solve().
-	std::vector<double> htrans_;
+        std::vector<double> htrans_;
-	std::vector<double> trans_ ;
+        std::vector<double> trans_ ;
-        std::vector<double> porevol_;
         std::vector<int> allcells_;
 
         // ------ Internal data for the cfs_tpfa_res solver. ------
-	struct cfs_tpfa_res_data* h_;
+        struct cfs_tpfa_res_data* h_;
 
         // ------ Data that will be modified for every solve. ------
         std::vector<double> wellperf_gpot_;
+        std::vector<double> initial_porevol_;
 
         // ------ Data that will be modified for every solver iteration. ------
         std::vector<double> cell_A_;
@@ -135,13 +147,15 @@ namespace Opm
         std::vector<double> face_gravcap_;
         std::vector<double> wellperf_A_;
         std::vector<double> wellperf_phasemob_;
+        std::vector<double> porevol_;   // Only modified if rock_comp_props_ is non-null.
+        std::vector<double> rock_comp_; // Empty unless rock_comp_props_ is non-null.
         // The update to be applied to the pressures (cell and bhp).
         std::vector<double> pressure_increment_;
-
+        // True if the matrix assembled would be singular but for the
-
+        // adjustment made in the cfs_*_assemble() calls. This happens
-
+        // if everything is incompressible and there are no pressure
-
+        // conditions.
-
+        bool singular_;
     };
 
 } // namespace Opm

opm/core/pressure/tpfa/cfs_tpfa_residual.c

@@ -1156,7 +1156,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid   *G      ,
 
 
 /* ---------------------------------------------------------------------- */
-void
+int
 cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
                       double                       dt       ,
                       struct cfs_tpfa_res_forces  *forces   ,
@@ -1170,7 +1170,7 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
                       struct cfs_tpfa_res_data    *h        )
 /* ---------------------------------------------------------------------- */
 {
-    int res_is_neumann, well_is_neumann, c, np2;
+    int res_is_neumann, well_is_neumann, c, np2, singular;
 
     csrmatrix_zero(         h->J);
     vector_zero   (h->J->m, h->F);
@@ -1207,9 +1207,76 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
         assemble_sources(dt, forces->src, h);
     }
 
-    if (res_is_neumann && well_is_neumann && h->pimpl->is_incomp) {
+    singular = res_is_neumann && well_is_neumann && h->pimpl->is_incomp;
-        h->J->sa[0] *= 2;
+    if (singular) {
+        h->J->sa[0] *= 2.0;
     }
+
+    return singular;
+}
+
+
+/* ---------------------------------------------------------------------- */
+int
+cfs_tpfa_res_comprock_assemble(
+                      struct UnstructuredGrid     *G        ,
+                      double                       dt       ,
+                      struct cfs_tpfa_res_forces  *forces   ,
+                      const double                *zc       ,
+                      struct compr_quantities_gen *cq       ,
+                      const double                *trans    ,
+                      const double                *gravcap_f,
+                      const double                *cpress   ,
+                      const double                *wpress   ,
+                      const double                *porevol  ,
+                      const double                *porevol0 ,
+                      const double                *rock_comp,
+                      struct cfs_tpfa_res_data    *h        )
+/* ---------------------------------------------------------------------- */
+{
+    /* We want to add this term to the usual residual:
+     *
+     * (porevol(pressure)-porevol(initial_pressure))/dt.
+     *
+     * Its derivative (for the diagonal term of the Jacobian) is:
+     *
+     * porevol(pressure)*rock_comp(pressure)/dt
+     */
+
+    int     c, rock_is_incomp, singular;
+    size_t  j;
+    double  dpv;
+
+    /* Assemble usual system (without rock compressibility). */
+    singular = cfs_tpfa_res_assemble(G, dt, forces, zc, cq, trans, gravcap_f,
+                                     cpress, wpress, porevol0, h);
+
+    /* If we made a singularity-removing adjustment in the
+       regular assembly, we undo it here. */
+    if (singular) {
+        h->J->sa[0] /= 2.0;
+    }
+
+    /* Add new terms to residual and Jacobian. */
+    rock_is_incomp = 1;
+    for (c = 0; c < G->number_of_cells; c++) {
+        j = csrmatrix_elm_index(c, c, h->J);
+
+        dpv = (porevol[c] - porevol0[c]);
+        if (dpv != 0.0 || rock_comp[c] != 0.0) {
+            rock_is_incomp = 0;
+        }
+
+        h->J->sa[j] += porevol[c] * rock_comp[c];
+        h->F[c]     += dpv;
+    }
+
+    /* Re-do the singularity-removing adjustment if necessary */
+    if (rock_is_incomp && singular) {
+        h->J->sa[0] *= 2.0;
+    }
+
+    return rock_is_incomp && singular;
 }
 
 

opm/core/pressure/tpfa/cfs_tpfa_residual.h

@@ -59,7 +59,11 @@ cfs_tpfa_res_construct(struct UnstructuredGrid   *G      ,
 void
 cfs_tpfa_res_destroy(struct cfs_tpfa_res_data *h);
 
-void
+/* Return value is 1 if the assembled matrix was adjusted to remove a
+   singularity. This happens if all fluids are incompressible and
+   there are no pressure conditions on wells or boundaries.
+   Otherwise return 0.  */
+int
 cfs_tpfa_res_assemble(struct UnstructuredGrid     *G,
                       double                       dt,
                       struct cfs_tpfa_res_forces  *forces,
@@ -72,6 +76,27 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G,
                       const double                *porevol,
                       struct cfs_tpfa_res_data    *h);
 
+/* Return value is 1 if the assembled matrix was adjusted to remove a
+   singularity. This happens if all fluids are incompressible, the
+   rock is incompressible, and there are no pressure conditions on
+   wells or boundaries.
+   Otherwise return 0.  */
+int
+cfs_tpfa_res_comprock_assemble(
+                      struct UnstructuredGrid     *G,
+                      double                       dt,
+                      struct cfs_tpfa_res_forces  *forces,
+                      const double                *zc,
+                      struct compr_quantities_gen *cq,
+                      const double                *trans,
+                      const double                *gravcap_f,
+                      const double                *cpress,
+                      const double                *wpress,
+                      const double                *porevol,
+                      const double                *porevol0,
+                      const double                *rock_comp,
+                      struct cfs_tpfa_res_data    *h);
+
 void
 cfs_tpfa_res_flux(struct UnstructuredGrid    *G        ,
                   struct cfs_tpfa_res_forces *forces   ,

opm/core/pressure/tpfa/ifs_tpfa.c

@@ -771,7 +771,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid      *G        ,
     assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok);
 
     /* We want to solve a Newton step for the residual
-     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible
+     * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_incompressible
      *
      */
 

opm/core/simulator/SimulatorCompressibleTwophase.cpp

@@ -0,0 +1,533 @@
+/*
+  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/core/utility/ErrorMacros.hpp>
+
+#include <opm/core/pressure/CompressibleTpfa.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/newwells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/simulator/SimulatorReport.hpp>
+#include <opm/core/simulator/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/writeVtkData.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
+
+#include <opm/core/wells/WellsManager.hpp>
+
+#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
+#include <opm/core/fluid/RockCompressibility.hpp>
+
+#include <opm/core/utility/ColumnExtract.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/lexical_cast.hpp>
+
+#include <numeric>
+#include <fstream>
+
+
+namespace Opm
+{
+
+    class SimulatorCompressibleTwophase::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const BlackoilPropertiesInterface& props,
+             const RockCompressibility* rock_comp,
+             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_;
+        WellsManager& wells_manager_;
+        const Wells* wells_;
+        const std::vector<double>& src_;
+        const FlowBoundaryConditions* bcs_;
+        const double* gravity_;
+        // Solvers
+        CompressibleTpfa psolver_;
+        TransportModelCompressibleTwophase 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,
+                                                                 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, wells_manager, src, bcs, linsolver, gravity));
+    }
+
+
+
+
+
+    SimulatorReport SimulatorCompressibleTwophase::run(SimulatorTimer& timer,
+                                                       BlackoilState& state,
+                                                       WellState& well_state)
+    {
+        return pimpl_->run(timer, state, well_state);
+    }
+
+
+
+    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 (...) {
+          THROW("Creating directories failed: " << fpath);
+        }
+        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
+        std::ofstream vtkfile(vtkfilename.str().c_str());
+        if (!vtkfile) {
+            THROW("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);
+    }
+
+
+    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();
+        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 (...) {
+              THROW("Creating directories failed: " << fpath);
+            }
+            fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
+            std::ofstream file(fname.str().c_str());
+            if (!file) {
+                THROW("Failed to open " << fname.str());
+            }
+            const std::vector<double>& d = *(it->second);
+            std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
+        }
+    }
+
+
+    static void outputWaterCut(const Opm::Watercut& watercut,
+                               const std::string& output_dir)
+    {
+        // Write water cut curve.
+        std::string fname = output_dir  + "/watercut.txt";
+        std::ofstream os(fname.c_str());
+        if (!os) {
+            THROW("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) {
+            THROW("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,
+                                              WellsManager& wells_manager,
+                                              const std::vector<double>& src,
+                                              const FlowBoundaryConditions* bcs,
+                                              LinearSolverInterface& linsolver,
+                                              const double* gravity)
+        : grid_(grid),
+          props_(props),
+          rock_comp_(rock_comp),
+          wells_manager_(wells_manager),
+          wells_(wells_manager.c_wells()),
+          src_(src),
+          bcs_(bcs),
+          psolver_(grid, props, rock_comp, 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 (...) {
+                THROW("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_ && rock_comp_->isActive()) {
+            computePorevolume(grid_, props_.porosity(), *rock_comp_, 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_satvol[2] = { 0.0 };
+        double satvol[2] = { 0.0 };
+        double injected[2] = { 0.0 };
+        double produced[2] = { 0.0 };
+        double tot_injected[2] = { 0.0 };
+        double tot_produced[2] = { 0.0 };
+        Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
+        std::cout << "\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.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_) {
+                    outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+                }
+                outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+
+            SimulatorReport sreport;
+
+            // Solve pressure equation.
+            if (check_well_controls_) {
+                computeFractionalFlow(props_, allcells_,
+                                      state.pressure(), 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_,
+                                                      fractional_flows,
+                                                      well_state.perfRates(),
+                                                      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_) {
+                        THROW("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_ && rock_comp_->isActive()) {
+                initial_porevol = porevol;
+                computePorevolume(grid_, props_.porosity(), *rock_comp_, 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_);
+                std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
+            }
+            for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
+                tsolver_.solve(&state.faceflux()[0], &state.pressure()[0],
+                               &porevol[0],  &initial_porevol[0], &transport_src[0], stepsize,
+                               state.saturation(), state.surfacevol());
+                Opm::computeInjectedProduced(props_,
+                                             state.pressure(), state.surfacevol(), state.saturation(),
+                                             transport_src, stepsize, injected, produced);
+                if (gravity_ != 0 && use_segregation_split_) {
+                    tsolver_.solveGravity(columns_, &state.pressure()[0], &initial_porevol[0],
+                                          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.saturation(), satvol);
+            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 << "\nVolume balance report (all numbers relative to total pore volume).\n";
+            std::cout << "    Saturated volumes:     "
+                      << std::setw(width) << satvol[0]/tot_porevol_init
+                      << std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
+            std::cout << "    Injected volumes:      "
+                      << std::setw(width) << injected[0]/tot_porevol_init
+                      << std::setw(width) << injected[1]/tot_porevol_init << std::endl;
+            std::cout << "    Produced volumes:      "
+                      << std::setw(width) << produced[0]/tot_porevol_init
+                      << std::setw(width) << produced[1]/tot_porevol_init << std::endl;
+            std::cout << "    Total inj volumes:     "
+                      << std::setw(width) << tot_injected[0]/tot_porevol_init
+                      << std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
+            std::cout << "    Total prod volumes:    "
+                      << std::setw(width) << tot_produced[0]/tot_porevol_init
+                      << std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
+            std::cout << "    In-place + prod - inj: "
+                      << std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
+                      << std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
+            std::cout << "    Init - now - pr + inj: "
+                      << std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
+                      << std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
+                      << std::endl;
+            std::cout.precision(8);
+
+            watercut.push(timer.currentTime() + timer.currentStepLength(),
+                          produced[0]/(produced[0] + produced[1]),
+                          tot_produced[0]/tot_porevol_init);
+            if (wells_) {
+                wellreport.push(props_, *wells_,
+                                state.pressure(), state.surfacevol(), state.saturation(),
+                                timer.currentTime() + timer.currentStepLength(),
+                                well_state.bhp(), well_state.perfRates());
+            }
+            sreport.total_time =  step_timer.secsSinceStart();
+            if (output_) {
+                sreport.reportParam(tstep_os);
+            }
+        }
+
+        if (output_) {
+            if (output_vtk_) {
+                outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+            outputStateMatlab(grid_, state, 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();
+        return report;
+    }
+
+
+} // namespace Opm

opm/core/simulator/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 <boost/shared_ptr.hpp>
+#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     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,
+                                     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 scoped_ptr since scoped_ptr requires complete type for Impl.
+        boost::shared_ptr<Impl> pimpl_;
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORCOMPRESSIBLETWOPHASE_HEADER_INCLUDED

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -132,7 +132,7 @@ namespace Opm
         return pimpl_->run(timer, state, well_state);
     }
 
-    static void reportVolumes(std::ostream &os,double satvol[2],double tot_porevol_init,
+    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])
@@ -164,7 +164,7 @@ namespace Opm
            << std::endl;
         os.precision(8);
     }
-    
+
     static void outputStateVtk(const UnstructuredGrid& grid,
                                const Opm::TwophaseState& state,
                                const int step,
@@ -180,7 +180,7 @@ namespace Opm
         catch (...) {
             THROW("Creating directories failed: " << fpath);
         }
-        vtkfilename << "/" << std::setw(3) << std::setfill('0') << step << ".vtu";
+        vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
         std::ofstream vtkfile(vtkfilename.str().c_str());
         if (!vtkfile) {
             THROW("Failed to open " << vtkfilename.str());
@@ -193,7 +193,8 @@ namespace Opm
         dm["velocity"] = &cell_velocity;
         Opm::writeVtkData(grid, dm, vtkfile);
     }
-    static void outputVectorMatlab(const std::string name,
+
+    static void outputVectorMatlab(const std::string& name,
                                    const std::vector<int>& vec,
                                    const int step,
                                    const std::string& output_dir)
@@ -214,29 +215,6 @@ namespace Opm
         }
         std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(file, "\n"));
     }
-    
-    static void outputVectorMatlab(const std::string name,
-                                   const std::vector<double> 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 (...) {
-            THROW("Creating directories failed: " << fpath);
-        }
-        fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt";
-        std::ofstream file(fname.str().c_str());
-        if (!file) {
-            THROW("Failed to open " << fname.str());
-        }
-        std::copy(vec.begin(), vec.end(), std::ostream_iterator<double>(file, "\n"));
-    }
-
 
     static void outputStateMatlab(const UnstructuredGrid& grid,
                                   const Opm::TwophaseState& state,
@@ -408,7 +386,7 @@ namespace Opm
             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;
@@ -528,6 +506,7 @@ namespace Opm
 
             // Update pore volumes if rock is compressible.
             if (rock_comp_ && rock_comp_->isActive()) {
+                initial_porevol = porevol;
                 computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol);
             }
 
@@ -543,11 +522,11 @@ namespace Opm
                 std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
             }
             for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_.solve(&state.faceflux()[0], &porevol[0], &transport_src[0],
+                tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0],
-                               stepsize, state.saturation());
+                              stepsize, state.saturation());
                 Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced);
                 if (use_segregation_split_) {
-                    tsolver_.solveGravity(columns_, &porevol[0], stepsize, state.saturation());
+                    tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
                 }
                 watercut.push(timer.currentTime() + timer.currentStepLength(),
                               produced[0]/(produced[0] + produced[1]),
@@ -556,7 +535,7 @@ namespace Opm
                     wellreport.push(props_, *wells_, state.saturation(),
                                     timer.currentTime() + timer.currentStepLength(),
                                     well_state.bhp(), well_state.perfRates());
-                }                            
+                }
             }
             transport_timer.stop();
             double tt = transport_timer.secsSinceStart();
@@ -569,11 +548,10 @@ namespace Opm
             tot_injected[1] += injected[1];
             tot_produced[0] += produced[0];
             tot_produced[1] += produced[1];
-            //reportVolumes(std::cout,satvol[0],&tot_porevol_init[0],&tot_injected[0]);
+            reportVolumes(std::cout,satvol, tot_porevol_init,
-            Opm::reportVolumes(std::cout,satvol, tot_porevol_init,
+                          tot_injected, tot_produced,
-                               tot_injected, tot_produced,
+                          injected, produced,
-                               injected, produced,
+                          init_satvol);
-                               init_satvol);
             sreport.total_time =  step_timer.secsSinceStart();
             if (output_) {
                 sreport.reportParam(tstep_os);

opm/core/simulator/SimulatorIncompTwophase.hpp

@@ -96,4 +96,4 @@ namespace Opm
 
 } // namespace Opm
 
-#endif // OPM_SIMULATORTWOPHASE_HEADER_INCLUDED
+#endif // OPM_SIMULATORINCOMPTWOPHASE_HEADER_INCLUDED

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -24,6 +24,7 @@
 #include <opm/core/transport/reorder/reordersequence.h>
 #include <opm/core/utility/RootFinders.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 
 #include <fstream>
@@ -38,7 +39,8 @@ namespace Opm
     typedef RegulaFalsi<WarnAndContinueOnError> RootFinder;
 
 
-    TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(const UnstructuredGrid& grid,
+    TransportModelCompressibleTwophase::TransportModelCompressibleTwophase(
+                                                   const UnstructuredGrid& grid,
                                                    const Opm::BlackoilPropertiesInterface& props,
                                                    const double tol,
                                                    const int maxit)
@@ -51,6 +53,7 @@ namespace Opm
           dt_(0.0),
           saturation_(grid.number_of_cells, -1.0),
           fractionalflow_(grid.number_of_cells, -1.0),
+          gravity_(0),
           mob_(2*grid.number_of_cells, -1.0),
           ia_upw_(grid.number_of_cells + 1, -1),
           ja_upw_(grid.number_of_faces, -1),
@@ -75,15 +78,15 @@ namespace Opm
 
     void TransportModelCompressibleTwophase::solve(const double* darcyflux,
                                                    const double* pressure,
-                                                   const double* surfacevol0,
                                                    const double* porevolume0,
                                                    const double* porevolume,
                                                    const double* source,
                                                    const double dt,
-                                                   std::vector<double>& saturation)
+                                                   std::vector<double>& saturation,
+                                                   std::vector<double>& surfacevol)
     {
         darcyflux_ = darcyflux;
-        surfacevol0_ = surfacevol0;
+        surfacevol0_ = &surfacevol[0];
         porevolume0_ = porevolume0;
         porevolume_ = porevolume;
         source_ = source;
@@ -93,6 +96,11 @@ namespace Opm
         props_.viscosity(props_.numCells(), pressure, NULL, &allcells_[0], &visc_[0], NULL);
         props_.matrix(props_.numCells(), pressure, NULL, &allcells_[0], &A_[0], NULL);
 
+        // Check immiscibility requirement (only done for first cell).
+        if (A_[1] != 0.0 || A_[2] != 0.0) {
+            THROW("TransportModelCompressibleTwophase requires a property object without miscibility.");
+        }
+
         std::vector<int> seq(grid_.number_of_cells);
         std::vector<int> comp(grid_.number_of_cells + 1);
         int ncomp;
@@ -107,6 +115,9 @@ namespace Opm
                                &ia_downw_[0], &ja_downw_[0]);
         reorderAndTransport(grid_, darcyflux);
         toBothSat(saturation_, saturation);
+
+        // Compute surface volume as a postprocessing step from saturation and A_
+        computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]);
     }
 
     // Residual function r(s) for a single-cell implicit Euler transport
@@ -375,20 +386,24 @@ namespace Opm
 
 
 
-    void TransportModelCompressibleTwophase::initGravity()
+    void TransportModelCompressibleTwophase::initGravity(const double* grav)
     {
         // Set up transmissibilities.
         std::vector<double> htrans(grid_.cell_facepos[grid_.number_of_cells]);
         const int nf = grid_.number_of_faces;
         trans_.resize(nf);
+        gravflux_.resize(nf);
         tpfa_htrans_compute(const_cast<UnstructuredGrid*>(&grid_), props_.permeability(), &htrans[0]);
         tpfa_trans_compute(const_cast<UnstructuredGrid*>(&grid_), &htrans[0], &trans_[0]);
+
+        // Remember gravity vector.
+        gravity_ = grav;
     }
 
 
 
 
-    void TransportModelCompressibleTwophase::initGravityDynamic(const double* grav)
+    void TransportModelCompressibleTwophase::initGravityDynamic()
     {
         // Set up gravflux_ = T_ij g [   (b_w,i rho_w,S - b_o,i rho_o,S) (z_i - z_f)
         //                             + (b_w,j rho_w,S - b_o,j rho_o,S) (z_f - z_j) ]
@@ -410,7 +425,7 @@ namespace Opm
                 for (int ci = 0; ci < 2; ++ci) {
                     double gdz = 0.0;
                     for (int d = 0; d < dim; ++d) {
-                        gdz += grav[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
+                        gdz += gravity_[d]*(grid_.cell_centroids[dim*c[ci] + d] - grid_.face_centroids[dim*f + d]);
                     }
                     gravflux_[f] += signs[ci]*trans_[f]*gdz*(density_[2*c[ci]] - density_[2*c[ci] + 1]);
                 }
@@ -494,11 +509,11 @@ namespace Opm
                                                           const double* pressure,
                                                           const double* porevolume0,
                                                           const double dt,
-                                                          const double* grav,
+                                                          std::vector<double>& saturation,
-                                                          std::vector<double>& saturation)
+                                                          std::vector<double>& surfacevol)
     {
         // Assume that solve() has already been called, so that A_ is current.
-        initGravityDynamic(grav);
+        initGravityDynamic();
 
         // Initialize mobilities.
         const int nc = grid_.number_of_cells;
@@ -531,6 +546,9 @@ namespace Opm
         std::cout << "Gauss-Seidel column solver average iterations: "
                   << double(num_iters)/double(columns.size()) << std::endl;
         toBothSat(saturation_, saturation);
+
+        // Compute surface volume as a postprocessing step from saturation and A_
+        computeSurfacevol(grid_.number_of_cells, props_.numPhases(), &A_[0], &saturation[0], &surfacevol[0]);
     }
 
 } // namespace Opm

opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp

@@ -30,6 +30,8 @@ namespace Opm
 
     class BlackoilPropertiesInterface;
 
+    /// Implements a reordering transport solver for compressible,
+    /// non-miscible two-phase flow.
     class TransportModelCompressibleTwophase : public TransportModelInterface
     {
     public:
@@ -52,17 +54,19 @@ namespace Opm
         /// \param[in] source            Transport source term.
         /// \param[in] dt                Time step.
         /// \param[in, out] saturation   Phase saturations.
+        /// \param[in, out] surfacevol   Surface volume densities for each phase.
         void solve(const double* darcyflux,
                    const double* pressure,
-                   const double* surfacevol0,
                    const double* porevolume0,
                    const double* porevolume,
                    const double* source,
                    const double dt,
-                   std::vector<double>& saturation);
+                   std::vector<double>& saturation,
+                   std::vector<double>& surfacevol);
 
         /// Initialise quantities needed by gravity solver.
-        void initGravity();
+        /// \param[in] grav    Gravity vector
+        void initGravity(const double* grav);
 
         /// Solve for gravity segregation.
         /// This uses a column-wise nonlinear Gauss-Seidel approach.
@@ -72,14 +76,14 @@ namespace Opm
         /// \param[in] columns           Vector of cell-columns.
         /// \param[in] porevolume0       Array of pore volumes at start of timestep.
         /// \param[in] dt                Time step.
-        /// \param[in] grav              Gravity vector.
         /// \param[in, out] saturation   Phase saturations.
+        /// \param[in, out] surfacevol   Surface volume densities for each phase.
         void solveGravity(const std::vector<std::vector<int> >& columns,
                           const double* pressure,
                           const double* porevolume0,
                           const double dt,
-                          const double* grav,
+                          std::vector<double>& saturation,
-                          std::vector<double>& saturation);
+                          std::vector<double>& surfacevol);
 
     private:
         virtual void solveSingleCell(const int cell);
@@ -88,7 +92,7 @@ namespace Opm
                                     const int pos,
                                     const double* gravflux);
         int solveGravityColumn(const std::vector<int>& cells);
-        void initGravityDynamic(const double* grav);
+        void initGravityDynamic();
 
     private:
         const UnstructuredGrid& grid_;
@@ -110,6 +114,7 @@ namespace Opm
         std::vector<double> saturation_;        // P (= num. phases) per cell
         std::vector<double> fractionalflow_;  // = m[0]/(m[0] + m[1]) per cell
         // For gravity segregation.
+        const double* gravity_;
         std::vector<double> trans_;
         std::vector<double> density_;
         std::vector<double> gravflux_;

opm/core/transport/reorder/TransportModelTwophase.cpp

@@ -52,8 +52,8 @@ namespace Opm
           source_(0),
           dt_(0.0),
           saturation_(grid.number_of_cells, -1.0),
-          reorder_iterations_(grid.number_of_cells, 0),
           fractionalflow_(grid.number_of_cells, -1.0),
+          reorder_iterations_(grid.number_of_cells, 0),
           mob_(2*grid.number_of_cells, -1.0)
 #ifdef EXPERIMENT_GAUSS_SEIDEL
         , ia_upw_(grid.number_of_cells + 1, -1),
@@ -107,6 +107,13 @@ namespace Opm
         toBothSat(saturation_, saturation);
     }
 
+
+    const std::vector<int>& TransportModelTwophase::getReorderIterations() const
+    {
+        return reorder_iterations_;
+    }
+
+
     // Residual function r(s) for a single-cell implicit Euler transport
     //
     //     r(s) = s - s0 + dt/pv*( influx + outflux*f(s) )
@@ -645,10 +652,7 @@ namespace Opm
 
         toBothSat(saturation_, saturation);
     }
-    void TransportModelTwophase::getReorderIterations()
+
-    {
-        return reorder_iterations_;
-    }; 
 } // namespace Opm
 
 

opm/core/transport/reorder/TransportModelTwophase.hpp

@@ -74,10 +74,11 @@ namespace Opm
                           const double* porevolume,
                           const double dt,
                           std::vector<double>& saturation);
-        //// Return reorder iterations
+
-        ////
+        //// Return the number of iterations used by the reordering solver.
-        //// \param[out] vector of iteration per cell
+        //// \return vector of iteration per cell
-        const std::vector<int>& getReorderIterations(){return reorder_iterations_;};      
+        const std::vector<int>& getReorderIterations() const;
+
     private:
         virtual void solveSingleCell(const int cell);
         virtual void solveMultiCell(const int num_cells, const int* cells);
@@ -86,7 +87,7 @@ namespace Opm
                                     const int pos,
                                     const double* gravflux);
         int solveGravityColumn(const std::vector<int>& cells);
-    private:         
+    private:
         const UnstructuredGrid& grid_;
         const IncompPropertiesInterface& props_;
         const double* visc_;

opm/core/utility/miscUtilitiesBlackoil.cpp

@@ -48,39 +48,39 @@ namespace Opm
     void computeInjectedProduced(const BlackoilPropertiesInterface& props,
                                  const std::vector<double>& press,
                                  const std::vector<double>& z,
-				 const std::vector<double>& s,
+                                 const std::vector<double>& s,
-				 const std::vector<double>& src,
+                                 const std::vector<double>& src,
-				 const double dt,
+                                 const double dt,
-				 double* injected,
+                                 double* injected,
-				 double* produced)
+                                 double* produced)
     {
-	const int num_cells = src.size();
+        const int num_cells = src.size();
-	const int np = s.size()/src.size();
+        const int np = s.size()/src.size();
-	if (int(s.size()) != num_cells*np) {
+        if (int(s.size()) != num_cells*np) {
-	    THROW("Sizes of s and src vectors do not match.");
+            THROW("Sizes of s and src vectors do not match.");
-	}
+        }
-	std::fill(injected, injected + np, 0.0);
+        std::fill(injected, injected + np, 0.0);
-	std::fill(produced, produced + np, 0.0);
+        std::fill(produced, produced + np, 0.0);
         std::vector<double> visc(np);
-	std::vector<double> mob(np);
+        std::vector<double> mob(np);
-	for (int c = 0; c < num_cells; ++c) {
+        for (int c = 0; c < num_cells; ++c) {
-	    if (src[c] > 0.0) {
+            if (src[c] > 0.0) {
-		injected[0] += src[c]*dt;
+                injected[0] += src[c]*dt;
-	    } else if (src[c] < 0.0) {
+            } else if (src[c] < 0.0) {
-		const double flux = -src[c]*dt;
+                const double flux = -src[c]*dt;
-		const double* sat = &s[np*c];
+                const double* sat = &s[np*c];
-		props.relperm(1, sat, &c, &mob[0], 0);
+                props.relperm(1, sat, &c, &mob[0], 0);
                 props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
-		double totmob = 0.0;
+                double totmob = 0.0;
-		for (int p = 0; p < np; ++p) {
+                for (int p = 0; p < np; ++p) {
-		    mob[p] /= visc[p];
+                    mob[p] /= visc[p];
-		    totmob += mob[p];
+                    totmob += mob[p];
-		}
+                }
-		for (int p = 0; p < np; ++p) {
+                for (int p = 0; p < np; ++p) {
-		    produced[p] += (mob[p]/totmob)*flux;
+                    produced[p] += (mob[p]/totmob)*flux;
-		}
+                }
-	    }
+            }
-	}
+        }
     }
 
 
@@ -93,11 +93,11 @@ namespace Opm
     /// @param[in]  s         saturation values (for all phases)
     /// @param[out] totmob    total mobilities.
     void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
-			      const std::vector<int>& cells,
+                              const std::vector<int>& cells,
                               const std::vector<double>& press,
                               const std::vector<double>& z,
-			      const std::vector<double>& s,
+                              const std::vector<double>& s,
-			      std::vector<double>& totmob)
+                              std::vector<double>& totmob)
     {
         std::vector<double> pmobc;
 
@@ -126,12 +126,12 @@ namespace Opm
     /// @param[out] totmob    total mobility
     /// @param[out] omega     fractional-flow weighted fluid densities.
     void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props,
-				   const std::vector<int>& cells,
+                                   const std::vector<int>& cells,
                                    const std::vector<double>& p,
                                    const std::vector<double>& z,
-				   const std::vector<double>& s,
+                                   const std::vector<double>& s,
-				   std::vector<double>& totmob,
+                                   std::vector<double>& totmob,
-				   std::vector<double>& omega)
+                                   std::vector<double>& omega)
     {
         std::vector<double> pmobc;
 
@@ -185,10 +185,10 @@ namespace Opm
         props.relperm(nc, &s[0], &cells[0],
                       &pmobc[0], dpmobc);
 
-	std::transform(pmobc.begin(), pmobc.end(),
+        std::transform(pmobc.begin(), pmobc.end(),
-		       mu.begin(),
+                       mu.begin(),
-		       pmobc.begin(),
+                       pmobc.begin(),
-		       std::divides<double>());
+                       std::divides<double>());
     }
 
     /// Computes the fractional flow for each cell in the cells argument
@@ -220,5 +220,35 @@ namespace Opm
         }
     }
 
+    /// Computes the surface volume densities from saturations by the formula
+    ///     z = A s
+    /// for a number of data points, where z is the surface volume density,
+    /// s is the saturation (both as column vectors) and A is the
+    /// phase-to-component relation matrix.
+    /// @param[in]  n            number of data points
+    /// @param[in]  np           number of phases, must be 2 or 3
+    /// @param[in]  A            array containing n square matrices of size num_phases^2,
+    ///                          in Fortran ordering, typically the output of a call
+    ///                          to the matrix() method of a BlackoilProperties* class.
+    /// @param[in]  saturation   concatenated saturation values (for all P phases)
+    /// @param[out] surfacevol   concatenated surface-volume values (for all P phases)
+    void computeSurfacevol(const int n,
+                           const int np,
+                           const double* A,
+                           const double* saturation,
+                           double* surfacevol)
+    {
+        // Note: since this is a simple matrix-vector product, it can
+        // be done by a BLAS call, but then we have to reorder the A
+        // matrix data.
+        std::fill(surfacevol, surfacevol + n*np, 0.0);
+        for (int i = 0; i < n; ++i) {
+            for (int col = 0; col < np; ++col) {
+                for (int row = 0; row < np; ++row) {
+                    surfacevol[i*np + row] += A[i*np*np + row + col*np] * saturation[i*np + col];
+                }
+            }
+        }
+    }
 
 } // namespace Opm

opm/core/utility/miscUtilitiesBlackoil.hpp

@@ -46,11 +46,11 @@ namespace Opm
     void computeInjectedProduced(const BlackoilPropertiesInterface& props,
                                  const std::vector<double>& p,
                                  const std::vector<double>& z,
-				 const std::vector<double>& s,
+                                 const std::vector<double>& s,
-				 const std::vector<double>& src,
+                                 const std::vector<double>& src,
-				 const double dt,
+                                 const double dt,
-				 double* injected,
+                                 double* injected,
-				 double* produced);
+                                 double* produced);
 
     /// @brief Computes total mobility for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
@@ -60,11 +60,11 @@ namespace Opm
     /// @param[in]  s         saturation values (for all phases)
     /// @param[out] totmob    total mobilities.
     void computeTotalMobility(const Opm::BlackoilPropertiesInterface& props,
-			      const std::vector<int>& cells,
+                              const std::vector<int>& cells,
                               const std::vector<double>& p,
                               const std::vector<double>& z,
-			      const std::vector<double>& s,
+                              const std::vector<double>& s,
-			      std::vector<double>& totmob);
+                              std::vector<double>& totmob);
 
     /// @brief Computes total mobility and omega for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
@@ -75,12 +75,12 @@ namespace Opm
     /// @param[out] totmob    total mobility
     /// @param[out] omega     fractional-flow weighted fluid densities.
     void computeTotalMobilityOmega(const Opm::BlackoilPropertiesInterface& props,
-				   const std::vector<int>& cells,
+                                   const std::vector<int>& cells,
                                    const std::vector<double>& p,
                                    const std::vector<double>& z,
-				   const std::vector<double>& s,
+                                   const std::vector<double>& s,
-				   std::vector<double>& totmob,
+                                   std::vector<double>& totmob,
-				   std::vector<double>& omega);
+                                   std::vector<double>& omega);
 
 
     /// @brief Computes phase mobilities for a set of saturation values.
@@ -96,7 +96,7 @@ namespace Opm
                                 const std::vector<double>&              z,
                                 const std::vector<double>&              s,
                                 std::vector<double>&                    pmobc);
-    
+
 
     /// Computes the fractional flow for each cell in the cells argument
     /// @param[in]  props            rock and fluid properties
@@ -112,6 +112,25 @@ namespace Opm
                                const std::vector<double>& s,
                                std::vector<double>& fractional_flows);
 
+
+    /// Computes the surface volume densities from saturations by the formula
+    ///     z = A s
+    /// for a number of data points, where z is the surface volume density,
+    /// s is the saturation (both as column vectors) and A is the
+    /// phase-to-component relation matrix.
+    /// @param[in]  n            number of data points
+    /// @param[in]  np           number of phases, must be 2 or 3
+    /// @param[in]  A            array containing n square matrices of size num_phases,
+    ///                          in Fortran ordering, typically the output of a call
+    ///                          to the matrix() method of a BlackoilProperties* class.
+    /// @param[in]  saturation   concatenated saturation values (for all P phases)
+    /// @param[out] surfacevol   concatenated surface-volume values (for all P phases)
+    void computeSurfacevol(const int n,
+                           const int np,
+                           const double* A,
+                           const double* saturation,
+                           double* surfacevol);
+
 } // namespace Opm
 
 #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED