merge.

Doxyfile

@@ -610,8 +610,8 @@ WARN_LOGFILE           =
 # directories like "/usr/src/myproject". Separate the files or directories
 # with spaces.
 
-#INPUT                  = opm/core tutorials examples
+INPUT                  = opm/core tutorials/tutorial1.cpp tutorials/tutorial2.cpp tutorials/tutorial3.cpp examples
-INPUT                  = tutorials
+#INPUT                  = tutorials
 
 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is

Makefile.am

@@ -51,6 +51,8 @@ opm/core/fluid/RockFromDeck.cpp                         \
 opm/core/fluid/SaturationPropsBasic.cpp                 \
 opm/core/fluid/SaturationPropsFromDeck.cpp              \
 opm/core/grid.c                                         \
+opm/core/newwells.c                                     \
+opm/core/GridManager.cpp                                \
 opm/core/grid/cart_grid.c                               \
 opm/core/grid/cornerpoint_grid.c                        \
 opm/core/grid/cpgpreprocess/geometry.c                  \
@@ -72,14 +74,12 @@ opm/core/utility/miscUtilities.cpp                      \
 opm/core/utility/miscUtilitiesBlackoil.cpp              \
 opm/core/utility/SimulatorTimer.cpp                     \
 opm/core/utility/StopWatch.cpp                          \
-opm/core/utility/newwells.c                             \
 opm/core/utility/writeVtkData.cpp			\
-opm/core/GridManager.cpp                                \
+opm/core/wells/WellsManager.cpp                         \
-opm/core/WellsManager.cpp                               \
+opm/core/wells/WellsGroup.cpp				\
-opm/core/WellsGroup.cpp					\
+opm/core/wells/WellCollection.cpp			\
-opm/core/WellCollection.cpp				\
+opm/core/wells/InjectionSpecification.cpp               \
-opm/core/InjectionSpecification.cpp                     \
+opm/core/wells/ProductionSpecification.cpp		\
-opm/core/ProductionSpecification.cpp			\
 opm/core/linalg/sparse_sys.c                            \
 opm/core/linalg/LinearSolverInterface.cpp               \
 opm/core/linalg/LinearSolverFactory.cpp                 \
@@ -109,6 +109,7 @@ opm/core/pressure/ifsh.c                                \
 opm/core/pressure/mimetic/mimetic.c                     \
 opm/core/pressure/mimetic/hybsys_global.c               \
 opm/core/pressure/mimetic/hybsys.c                      \
+opm/core/simulator/SimulatorTwophase.cpp                \
 opm/core/transport/spu_explicit.c                       \
 opm/core/transport/spu_implicit.c			\
 opm/core/transport/transport_source.c			\
@@ -158,6 +159,7 @@ opm/core/grid/cpgpreprocess/geometry.h                  \
 opm/core/grid/cpgpreprocess/facetopology.h              \
 opm/core/grid/cpgpreprocess/grdecl.h                    \
 opm/core/utility/Average.hpp                            \
+opm/core/utility/ColumnExtract.hpp			\
 opm/core/utility/ErrorMacros.hpp                        \
 opm/core/utility/Factory.hpp                            \
 opm/core/utility/MonotCubicInterpolator.hpp             \
@@ -193,17 +195,13 @@ opm/core/linalg/LinearSolverFactory.hpp                 \
 opm/core/grid.h                                         \
 opm/core/well.h                                         \
 opm/core/newwells.h                                     \
-opm/core/WellsGroup.hpp					\
+opm/core/wells/WellsManager.hpp                         \
-opm/core/WellCollection.hpp				\
+opm/core/wells/WellsGroup.hpp				\
-opm/core/InjectionSpecification.hpp			\
+opm/core/wells/WellCollection.hpp			\
-opm/core/ProductionSpecification.hpp			\
+opm/core/wells/InjectionSpecification.hpp		\
-opm/core/BlackoilState.hpp                              \
+opm/core/wells/ProductionSpecification.hpp		\
-opm/core/ColumnExtract.hpp				\
 opm/core/GridAdapter.hpp                                \
 opm/core/GridManager.hpp                                \
-opm/core/TwophaseState.hpp                              \
-opm/core/WellsManager.hpp                               \
-opm/core/WellState.hpp                                  \
 opm/core/pressure/fsh.h                                 \
 opm/core/pressure/HybridPressureSolver.hpp              \
 opm/core/pressure/IncompTpfa.hpp                        \
@@ -230,6 +228,10 @@ opm/core/pressure/fsh_common_impl.h                     \
 opm/core/pressure/mimetic/hybsys_global.h               \
 opm/core/pressure/mimetic/hybsys.h                      \
 opm/core/pressure/mimetic/mimetic.h                     \
+opm/core/simulator/BlackoilState.hpp                    \
+opm/core/simulator/SimulatorTwophase.hpp                \
+opm/core/simulator/TwophaseState.hpp                    \
+opm/core/simulator/WellState.hpp                        \
 opm/core/transport/ImplicitTransport.hpp                \
 opm/core/transport/CSRMatrixUmfpackSolver.hpp           \
 opm/core/transport/CSRMatrixBlockAssembler.hpp          \
@@ -271,3 +273,17 @@ opm/core/linalg/LinearSolverIstl.cpp
 nobase_include_HEADERS +=			\
 opm/core/linalg/LinearSolverIstl.hpp
 endif
+
+
+if BUILD_AGMG
+libopmcore_la_SOURCES +=			\
+$(AGMG_SRCDIR)/dagmg.f90			\
+$(AGMG_SRCDIR)/dagmg_mumps.f90			\
+opm/core/linalg/LinearSolverAGMG.cpp
+
+nobase_include_HEADERS +=			\
+opm/core/linalg/LinearSolverAGMG.hpp
+
+libopmcore_la_LDFLAGS +=			\
+$(FCLIBS)
+endif

configure.ac

@@ -8,32 +8,33 @@ AM_INIT_AUTOMAKE([-Wall -Werror foreign subdir-objects])
 
 m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])])
 
-m4_ifdef([LT_INIT], [LT_INIT])
-
 AC_CONFIG_MACRO_DIR([m4])
 AC_CONFIG_SRCDIR([opm/core/grid.h])
 AC_CONFIG_HEADERS([config.h])
 
 # Checks for programs.
-AC_PROG_CXX
 AC_PROG_CC
 AM_PROG_CC_C_O
-m4_ifdef([LT_INIT], [], [AC_PROG_LIBTOOL])
-AC_PROG_RANLIB
 
-# AX_LAPACK requires a working F77 compiler or, rather, its runtime
+m4_ifdef([LT_INIT],
-# support libraries.
+         [LT_INIT[]dnl
-AC_PROG_F77
+          LT_LANG([C++])dnl
-
+          LT_LANG([Fortran 77])dnl
-# F77 name mangling
+          LT_LANG([Fortran])dnl
-AC_F77_WRAPPERS
+         ],dnl
+         [AC_PROG_LIBTOOL[]dnl
+          AC_PROG_CXX[]dnl
+          AC_PROG_F77[]dnl
+          AC_PROG_FC[]dnl
+         ])[]dnl
 
 # Checks for libraries.
 
 # Bring in numerics support (standard library component)
 AC_SEARCH_LIBS([sqrt], [m])
 
-AX_LAPACK
+OPM_LAPACK
+
 AX_BOOST_BASE([1.37])
 AX_BOOST_SYSTEM
 AX_BOOST_DATE_TIME
@@ -41,6 +42,7 @@ AX_BOOST_FILESYSTEM
 AX_BOOST_UNIT_TEST_FRAMEWORK
 
 AX_DUNE_ISTL
+OPM_AGMG
 
 # Checks for header files.
 AC_CHECK_HEADERS([float.h limits.h stddef.h stdlib.h string.h])

examples/Makefile.am

@@ -1,31 +1,40 @@
+# Build-time flags needed to form example programs
 AM_CPPFLAGS =					\
 -I$(top_srcdir)					\
 $(BOOST_CPPFLAGS)
 
+# All targets link to the library
 LDADD = $(top_builddir)/libopmcore.la
 
+# ----------------------------------------------------------------------
+# Declare products (i.e., the example programs).
+#
+# Please keep the list sorted.
+
 noinst_PROGRAMS =				\
 refine_wells					\
 scaneclipsedeck					\
+sim_2p_incomp_reorder				\
 sim_wateroil					\
 wells_example
 
+# ----------------------------------------------------------------------
+# Product constituents.  Must be specified for every product that's
+# built from more than a single ".c" file and/or that link to anything
+# more than the OPM-Core library.
+#
+# Please maintain sort order from "noinst_PROGRAMS".
+
 refine_wells_SOURCES          = refine_wells.cpp
-
+sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp
 sim_wateroil_SOURCES          = sim_wateroil.cpp
-sim_wateroil_LDADD    =				\
-$(LDADD) $(LIBS)				\
-$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB)	\
-$(LAPACK_LIBS) $(LIBS) $(LIBS)
-
 wells_example_SOURCES         = wells_example.cpp
 
+# ----------------------------------------------------------------------
+# Optional examples, or examples that use optional add-on components.
+
 if UMFPACK
 noinst_PROGRAMS += spu_2p
 
 spu_2p_SOURCES   = spu_2p.cpp
-spu_2p_LDADD   =				\
-$(LDADD) $(LIBS)				\
-$(BOOST_SYSTEM_LIB) $(BOOST_FILESYSTEM_LIB)	\
-$(LAPACK_LIBS) $(LIBS) $(LIBS)
 endif

examples/sim_2p_incomp_reorder.cpp

@@ -0,0 +1,220 @@
+/*
+  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/IncompTpfa.hpp>
+#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/utility/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/utility/writeVtkData.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+
+#include <opm/core/fluid/IncompPropertiesBasic.hpp>
+#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
+#include <opm/core/fluid/RockCompressibility.hpp>
+
+#include <opm/core/linalg/LinearSolverFactory.hpp>
+
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/simulator/SimulatorTwophase.hpp>
+
+#include <boost/filesystem/convenience.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/lexical_cast.hpp>
+
+#include <cassert>
+#include <cstddef>
+
+#include <algorithm>
+#include <tr1/array>
+#include <functional>
+#include <iostream>
+#include <iomanip>
+#include <fstream>
+#include <iterator>
+#include <vector>
+#include <numeric>
+
+
+
+
+// ----------------- Main program -----------------
+int
+main(int argc, char** argv)
+{
+    using namespace Opm;
+
+    std::cout << "\n================    Test program for incompressible two-phase flow     ===============\n\n";
+    Opm::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<Opm::GridManager> grid;
+    boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
+    boost::scoped_ptr<Opm::WellsManager> wells;
+    boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
+    Opm::SimulatorTimer simtimer;
+    Opm::TwophaseState state;
+    // bool check_well_controls = false;
+    // int max_well_control_iterations = 0;
+    double gravity[3] = { 0.0 };
+    if (use_deck) {
+        std::string deck_filename = param.get<std::string>("deck_filename");
+        Opm::EclipseGridParser deck(deck_filename);
+        // Grid init
+        grid.reset(new Opm::GridManager(deck));
+        // Rock and fluid init
+        const int* gc = grid->c_grid()->global_cell;
+        std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
+        props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell));
+        // Wells init.
+        wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
+        // check_well_controls = param.getDefault("check_well_controls", false);
+        // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
+        // Timer init.
+        if (deck.hasField("TSTEP")) {
+            simtimer.init(deck);
+        } else {
+            simtimer.init(param);
+        }
+        // Rock compressibility.
+        rock_comp.reset(new Opm::RockCompressibility(deck));
+        // Gravity.
+        gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::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);
+        }
+    } 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 Opm::GridManager(nx, ny, nz, dx, dy, dz));
+        // Rock and fluid init.
+        props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
+        // Wells init.
+        wells.reset(new Opm::WellsManager());
+        // Timer init.
+        simtimer.init(param);
+        // Rock compressibility.
+        rock_comp.reset(new Opm::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);
+    }
+
+    // Warn if gravity but no density difference.
+    bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
+    if (use_gravity) {
+        if (props->density()[0] == props->density()[1]) {
+            std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
+        }
+    }
+
+    // Source-related variables init.
+    int num_cells = grid->c_grid()->number_of_cells;
+    std::vector<double> totmob;
+    std::vector<double> omega; // Will remain empty if no gravity.
+    std::vector<double> rc; // Will remain empty if no rock compressibility.
+
+    // Extra rock init.
+    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);
+    }
+    double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
+
+    // Initialising src
+    std::vector<double> src(num_cells, 0.0);
+    if (wells->c_wells()) {
+        // Do nothing, wells will be the driving force, not source terms.
+        // Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
+    } else {
+        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/Opm::unit::day;
+        src[0] = flow_per_sec;
+        src[num_cells - 1] = -flow_per_sec;
+    }
+
+    // Boundary conditions.
+    Opm::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(), Opm::FlowBCManager::Side(pside), pside_pressure);
+    }
+
+    // Linear solver.
+    Opm::LinearSolverFactory linsolver(param);
+
+    const double *grav = use_gravity ? &gravity[0] : 0;
+
+
+    Opm::SimulatorTwophase simulator(param,
+                                     *grid->c_grid(),
+                                     *props,
+                                     rock_comp->isActive() ? rock_comp.get() : 0,
+                                     wells->c_wells(),
+                                     src,
+                                     bcs.c_bcs(),
+                                     linsolver,
+                                     grav);
+
+    // Warn if any parameters are unused.
+    if (param.anyUnused()) {
+        std::cout << "--------------------   Unused parameters:   --------------------\n";
+        param.displayUsage();
+        std::cout << "----------------------------------------------------------------" << std::endl;
+    }
+
+    // Write parameters used for later reference.
+    // if (output) {
+    //     param.writeParam(output_dir + "/spu_2p.param");
+    // }
+
+    WellState well_state;
+    well_state.init(wells->c_wells(), state);
+
+    simulator.run(simtimer, state, well_state);
+}

examples/sim_wateroil.cpp

@@ -27,7 +27,7 @@
 #include <opm/core/grid.h>
 #include <opm/core/GridManager.hpp>
 #include <opm/core/newwells.h>
-#include <opm/core/WellsManager.hpp>
+#include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/initState.hpp>
 #include <opm/core/utility/SimulatorTimer.hpp>
@@ -44,9 +44,9 @@
 
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 
-#include <opm/core/ColumnExtract.hpp>
+#include <opm/core/utility/ColumnExtract.hpp>
-#include <opm/core/BlackoilState.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/WellState.hpp>
+#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/transport/GravityColumnSolver.hpp>
 
 #include <opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp>

examples/spu_2p.cpp

@@ -43,7 +43,7 @@
 #include <opm/core/grid.h>
 #include <opm/core/GridManager.hpp>
 #include <opm/core/newwells.h>
-#include <opm/core/WellsManager.hpp>
+#include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/initState.hpp>
 #include <opm/core/utility/SimulatorTimer.hpp>
@@ -69,8 +69,8 @@
 #include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
 #include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
 
-#include <opm/core/ColumnExtract.hpp>
+#include <opm/core/utility/ColumnExtract.hpp>
-#include <opm/core/TwophaseState.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/transport/GravityColumnSolver.hpp>
 
 #include <opm/core/transport/reorder/TransportModelTwophase.hpp>

examples/wells_example.cpp

@@ -5,14 +5,14 @@
 
 #include "opm/core/utility/initState.hpp"
 #include "opm/core/utility/SimulatorTimer.hpp"
-#include <opm/core/WellsManager.hpp>
+#include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/GridManager.hpp>
 #include <opm/core/pressure/IncompTpfa.hpp>
 #include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
 #include <opm/core/newwells.h>
 #include <opm/core/grid.h>
 #include <opm/core/utility/miscUtilities.hpp>
-#include <opm/core/TwophaseState.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/pressure/FlowBCManager.hpp>
 #include <opm/core/linalg/LinearSolverFactory.hpp>
 #include <opm/core/fluid/RockCompressibility.hpp>

m4/agmg.m4

@@ -0,0 +1,26 @@
+AC_DEFUN([OPM_AGMG],dnl
+[AC_ARG_WITH([agmg],dnl
+ [AS_HELP_STRING([--with-agmg=SRCDIR],dnl
+  [Include DOUBLE PRECISION version Notay's of AGMG Algebraic
+   Multigrid solver from specified source location SRCDIR. Note: this
+   option requires a complete, working Fortran 90 environment.])],
+ [],dnl
+ [with_agmg=no])[]dnl
+
+ AS_IF([test x"$with_agmg" != x"no" -a -d "$with_agmg"],dnl
+ [AS_IF([test -f "$with_agmg/dagmg.f90"],dnl
+  [AC_SUBST([AGMG_SRCDIR], [$with_agmg])[]dnl
+   AC_DEFINE([HAVE_AGMG], [1],dnl
+             [Define to `1' if Notay's AGMG solver is included])[]dnl
+   build_agmg="yes"],dnl
+  [build_agmg="no"])],dnl
+ [build_agmg="no"])[]dnl
+
+ AS_IF([test x"$build_agmg" = x"yes"],dnl
+   [AC_PROG_FC_C_O[]dnl
+    AC_FC_WRAPPERS[]dnl
+    AC_FC_LIBRARY_LDFLAGS[]dnl
+   ], [:])[]dnl
+
+ AM_CONDITIONAL([BUILD_AGMG], [test x"$build_agmg" = x"yes"])
+])[]dnl

m4/opm_lapack.m4

@@ -0,0 +1,4 @@
+AC_DEFUN([OPM_LAPACK],
+[AC_REQUIRE([AC_F77_WRAPPERS])dnl
+ AC_REQUIRE([AX_LAPACK])dnl
+])[]dnl

opm/core/doxygen_main.hpp

@@ -0,0 +1,35 @@
+/*
+  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_DOXYGEN_MAIN_HEADER_INCLUDED
+#define OPM_DOXYGEN_MAIN_HEADER_INCLUDED
+
+
+/** \mainpage Documentation for the opm-core library.
+
+The following are the main library features:
+<h3>Grid handling</h3>
+<h3>Well handling</h3>
+<h3>Pressure solvers</h3>
+<h3>Transport solvers</h3>
+<h3>Various utilities</h3>
+
+*/
+
+#endif // OPM_DOXYGEN_MAIN_HEADER_INCLUDED

opm/core/eclipse/EclipseGridParser.cpp

@@ -41,6 +41,7 @@
 #include <exception>
 #include <algorithm>
 #include <limits>
+#include <numeric>
 #include <cfloat>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 #include <opm/core/eclipse/EclipseGridParserHelpers.hpp>
@@ -97,6 +98,7 @@ namespace EclipseKeywords
           string("SGOF"),     string("SWOF"),   string("ROCK"),
           string("ROCKTAB"),  string("WELSPECS"), string("COMPDAT"),
           string("WCONINJE"), string("WCONPROD"), string("WELTARG"),
+          string("WELOPEN"),
           string("EQUIL"),    string("PVCDO"),    string("TSTEP"),
           string("PLYVISC"),  string("PLYROCK"),  string("PLYADS"),
           string("PLYMAX"),   string("TLMIXPAR"), string("WPOLYMER"),
@@ -147,6 +149,7 @@ namespace {
     enum FieldType {
         Integer,
         FloatingPoint,
+        Timestepping,
         SpecialField,
         IgnoreWithData,
         IgnoreNoData,
@@ -161,6 +164,8 @@ namespace {
             return Integer;
         } else if (count(floating_fields, floating_fields + num_floating_fields, keyword)) {
             return FloatingPoint;
+        } else if (keyword == "TSTEP" || keyword == "DATES") {
+            return Timestepping;
         } else if (count(special_fields, special_fields + num_special_fields, keyword)) {
             return SpecialField;
         } else if (count(ignore_with_data, ignore_with_data + num_ignore_with_data, keyword)) {
@@ -276,6 +281,10 @@ namespace {
 //---------------------------------------------------------------------------
 EclipseGridParser::EclipseGridParser()
 //---------------------------------------------------------------------------
+    : current_reading_mode_(Regular),
+      start_date_(boost::date_time::not_a_date_time),
+      current_time_days_(0.0),
+      current_epoch_(0)
 {
 }
 
@@ -284,6 +293,10 @@ EclipseGridParser::EclipseGridParser()
 //---------------------------------------------------------------------------
 EclipseGridParser::EclipseGridParser(const string& filename, bool convert_to_SI)
 //---------------------------------------------------------------------------
+    : current_reading_mode_(Regular),
+      start_date_(boost::date_time::not_a_date_time),
+      current_time_days_(0.0),
+      current_epoch_(0)
 {
     // Store directory of filename
     boost::filesystem::path p(filename);
@@ -304,13 +317,39 @@ void EclipseGridParser::read(istream& is, bool convert_to_SI)
 {
     integer_field_map_.clear();
     floating_field_map_.clear();
-    special_field_map_.clear();
+    special_field_by_epoch_.clear();
+    special_field_by_epoch_.push_back(SpecialMap());
 
     readImpl(is);
+
+    current_epoch_ = 0;
+
     computeUnits();
     if (convert_to_SI) {
         convertToSI();
     }
+
+#define VERBOSE_LIST_FIELDS 1
+#if VERBOSE_LIST_FIELDS
+    std::cout << "\nInteger fields:\n";
+    for (std::map<string, std::vector<int> >::iterator
+             i = integer_field_map_.begin(); i != integer_field_map_.end(); ++i)
+        std::cout << '\t' << i->first << '\n';
+
+    std::cout << "\nFloat fields:\n";
+    for (std::map<string, std::vector<double> >::iterator
+             i = floating_field_map_.begin(); i != floating_field_map_.end(); ++i)
+        std::cout << '\t' << i->first << '\n';
+
+    std::cout << "\nSpecial fields:\n";
+    for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) {
+        std::cout << "Epoch " << epoch << '\n';
+        const SpecialMap& sm = special_field_by_epoch_[epoch];
+        for (SpecialMap::const_iterator i = sm.begin(); i != sm.end(); ++i) {
+            std::cout << '\t' << i->first << '\n';
+        }
+    }
+#endif
 }
 
 //---------------------------------------------------------------------------
@@ -329,7 +368,6 @@ void EclipseGridParser::readImpl(istream& is)
     //       though (of course retaining the basic guarantee).
     map<string, vector<int> >& intmap = integer_field_map_;
     map<string, vector<double> >& floatmap = floating_field_map_;
-    map<string, tr1::shared_ptr<SpecialBase> >& specialmap = special_field_map_;
 
     // Actually read the data
     std::string keyword;
@@ -341,17 +379,77 @@ void EclipseGridParser::readImpl(istream& is)
             cout << "Keyword found: " << keyword << endl;
             //#endif
             FieldType type = classifyKeyword(keyword);
+            // std::cout << "Classification: " << type << std::endl;
             switch (type) {
-            case Integer:
+            case Integer: {
                 readVectorData(is, intmap[keyword]);
                 break;
-            case FloatingPoint:
+            }
+            case FloatingPoint: {
                 readVectorData(is, floatmap[keyword]);
                 break;
+            }
+            case Timestepping: {
+                SpecialMap& sm = special_field_by_epoch_[current_epoch_];
+                if (start_date_.is_not_a_date()) {
+                    // Set it to START date, or default if no START.
+                    // This will only ever happen in the first epoch,
+                    // upon first encountering a timestepping keyword.
+                    SpecialMap::const_iterator it = sm.find("START");
+                    if (hasField("START")) {
+                        start_date_ = getSTART().date;
+                    } else {
+                        start_date_ = boost::gregorian::date(1983, 1, 1);
+                    }
+                }
+                if (current_reading_mode_ == Regular) {
+                    current_reading_mode_ = Timesteps;
+                }
+                // Get current epoch's TSTEP, if it exists, create new if not.
+                SpecialMap::iterator it = sm.find("TSTEP");
+                TSTEP* tstep = 0;
+                if (it != sm.end()) {
+                    tstep = dynamic_cast<TSTEP*>(it->second.get());
+                } else {
+                    SpecialFieldPtr sb_ptr(new TSTEP());
+                    tstep = dynamic_cast<TSTEP*>(sb_ptr.get());
+                    sm["TSTEP"] = sb_ptr;
+                }
+                ASSERT(tstep != 0);
+                // Append new steps to current TSTEP object
+                if (keyword == "TSTEP") {
+                    const int num_steps_old = tstep->tstep_.size();
+                    tstep->read(is); // This will append to the TSTEP object.
+                    const double added_days
+                        = std::accumulate(tstep->tstep_.begin() + num_steps_old, tstep->tstep_.end(), 0.0);
+                    current_time_days_ += added_days;
+                } else if (keyword == "DATES") {
+                    DATES dates;
+                    dates.read(is);
+                    for (std::size_t dix = 0; dix < dates.dates.size(); ++dix) {
+                        boost::gregorian::date_duration since_start = dates.dates[dix] - start_date_;
+                        double step = double(since_start.days()) - current_time_days_;
+                        tstep->tstep_.push_back(step);
+                        current_time_days_ = double(since_start.days());
+                    }
+                } else {
+                    THROW("Keyword " << keyword << " cannot be handled here.");
+                }
+                break;
+            }
             case SpecialField: {
-                std::tr1::shared_ptr<SpecialBase> sb_ptr = createSpecialField(is, keyword);
+                if (current_reading_mode_ == Timesteps) {
+                    // We have been reading timesteps, but have
+                    // now encountered something else.
+                    // That means we are in a new epoch.
+                    current_reading_mode_ = Regular;
+                    special_field_by_epoch_.push_back(SpecialMap());
+                    ++current_epoch_;
+                    ASSERT(int(special_field_by_epoch_.size()) == current_epoch_ + 1);
+                }
+                SpecialFieldPtr sb_ptr = createSpecialField(is, keyword);
                 if (sb_ptr) {
-                    specialmap[keyword] = sb_ptr;
+                    special_field_by_epoch_[current_epoch_][keyword] = sb_ptr;
                 } else {
                     THROW("Could not create field " << keyword);
                 }
@@ -398,24 +496,6 @@ void EclipseGridParser::readImpl(istream& is)
             is >> ignoreLine;
         }
     }
-
-#define VERBOSE_LIST_FIELDS 0
-#if VERBOSE_LIST_FIELDS
-    std::cout << "\nInteger fields:\n";
-    for (std::map<string, std::vector<int> >::iterator
-             i = intmap.begin(); i != intmap.end(); ++i)
-        std::cout << '\t' << i->first << '\n';
-
-    std::cout << "\nFloat fields:\n";
-    for (std::map<string, std::vector<double> >::iterator
-             i = floatmap.begin(); i != floatmap.end(); ++i)
-        std::cout << '\t' << i->first << '\n';
-
-    std::cout << "\nSpecial fields:\n";
-    for (std::map<string, std::tr1::shared_ptr<SpecialBase> >::iterator
-              i = specialmap.begin(); i != specialmap.end(); ++i)
-        std::cout << '\t' << i->first << '\n';
-#endif
 }
 
 
@@ -425,10 +505,13 @@ void EclipseGridParser::convertToSI()
 //---------------------------------------------------------------------------
 {
     // Convert all special fields.
-    typedef std::map<string, std::tr1::shared_ptr<SpecialBase> >::iterator SpecialIt;
+    typedef SpecialMap::iterator SpecialIt;
-    for (SpecialIt i = special_field_map_.begin(); i != special_field_map_.end(); ++i) {
+    for (int epoch = 0; epoch < numberOfEpochs(); ++epoch) {
+        SpecialMap& sm = special_field_by_epoch_[epoch];
+        for (SpecialIt i = sm.begin(); i != sm.end(); ++i) {
             i->second->convertToSI(units_);
         }
+    }
 
     // Convert all floating point fields.
     typedef std::map<string, std::vector<double> >::iterator FloatIt;
@@ -479,7 +562,7 @@ bool EclipseGridParser::hasField(const string& keyword) const
 {
     string ukey = upcase(keyword);
     return integer_field_map_.count(ukey) || floating_field_map_.count(ukey) ||
-        special_field_map_.count(ukey) || ignored_fields_.count(ukey);
+        special_field_by_epoch_[current_epoch_].count(ukey) || ignored_fields_.count(ukey);
 }
 
 
@@ -505,7 +588,7 @@ vector<string> EclipseGridParser::fieldNames() const
     vector<string> names;
     names.reserve(integer_field_map_.size() +
                   floating_field_map_.size() +
-                  special_field_map_.size() +
+                  special_field_by_epoch_[current_epoch_].size() +
                   ignored_fields_.size());
     {
         map<string, vector<int> >::const_iterator it = integer_field_map_.begin();
@@ -520,8 +603,8 @@ vector<string> EclipseGridParser::fieldNames() const
         }
     }
     {
-        map<string, std::tr1::shared_ptr<SpecialBase> >::const_iterator it = special_field_map_.begin();
+        SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].begin();
-        for (; it != special_field_map_.end(); ++it) {
+        for (; it != special_field_by_epoch_[current_epoch_].end(); ++it) {
             names.push_back(it->first);
         }
     }
@@ -534,6 +617,24 @@ vector<string> EclipseGridParser::fieldNames() const
     return names;
 }
 
+
+//---------------------------------------------------------------------------
+int EclipseGridParser::numberOfEpochs() const
+//---------------------------------------------------------------------------
+{
+    return special_field_by_epoch_.size();
+}
+
+
+//---------------------------------------------------------------------------
+void EclipseGridParser::setCurrentEpoch(int epoch)
+//---------------------------------------------------------------------------
+{
+    ASSERT(epoch >= 0 && epoch < numberOfEpochs());
+    current_epoch_ = epoch;
+}
+
+
 //---------------------------------------------------------------------------
 const std::vector<int>& EclipseGridParser::getIntegerValue(const std::string& keyword) const
 //---------------------------------------------------------------------------
@@ -574,8 +675,8 @@ const std::vector<double>& EclipseGridParser::getFloatingPointValue(const std::s
 const std::tr1::shared_ptr<SpecialBase> EclipseGridParser::getSpecialValue(const std::string& keyword) const
 //---------------------------------------------------------------------------
 {
-    map<string, std::tr1::shared_ptr<SpecialBase> >::const_iterator it = special_field_map_.find(keyword);
+    SpecialMap::const_iterator it = special_field_by_epoch_[current_epoch_].find(keyword);
-    if (it == special_field_map_.end()) {
+    if (it == special_field_by_epoch_[current_epoch_].end()) {
         THROW("No such field: " << keyword);
     } else {
         return it->second;
@@ -617,7 +718,7 @@ void EclipseGridParser::setSpecialField(const std::string& keyword,
                                         std::tr1::shared_ptr<SpecialBase> field)
 //---------------------------------------------------------------------------
 {
-    special_field_map_[keyword] = field;
+    special_field_by_epoch_[current_epoch_][keyword] = field;
 }
 
 //---------------------------------------------------------------------------

opm/core/eclipse/EclipseGridParser.hpp

@@ -94,6 +94,14 @@ public:
     /// The keywords/fields found in the file.
     std::vector<std::string> fieldNames() const;
 
+    /// Returns the number of distinct epochs found in the deck SCHEDULE.
+    int numberOfEpochs() const;
+
+    /// Sets the current epoch.
+    /// Valid arguments are in [0, ..., numberOfEpochs() - 1].
+    /// After reading, current epoch always starts at 0.
+    void setCurrentEpoch(int epoch);
+
     /// Returns a reference to a vector containing the values
     /// corresponding to the given integer keyword.
     const std::vector<int>& getIntegerValue(const std::string& keyword) const;
@@ -102,10 +110,12 @@ public:
     /// corresponding to the given floating-point keyword.
     const std::vector<double>& getFloatingPointValue(const std::string& keyword) const;
 
+    typedef std::tr1::shared_ptr<SpecialBase> SpecialFieldPtr;
+
     /// Returns a reference to a vector containing pointers to the values 
     /// corresponding to the given keyword when the values are not only integers
     /// or floats.
-    const std::tr1::shared_ptr<SpecialBase> getSpecialValue(const std::string& keyword) const;
+    const SpecialFieldPtr getSpecialValue(const std::string& keyword) const;
 
     // This macro implements support for a special field keyword. It requires that a subclass
     // of SpecialBase exists, that has the same name as the keyword.
@@ -140,6 +150,7 @@ public:
     SPECIAL_FIELD(WCONINJE);
     SPECIAL_FIELD(WCONPROD);
     SPECIAL_FIELD(WELTARG);
+    SPECIAL_FIELD(WELOPEN);
     SPECIAL_FIELD(EQUIL);
     SPECIAL_FIELD(PVCDO);
     SPECIAL_FIELD(TSTEP);
@@ -170,7 +181,7 @@ public:
     void setFloatingPointField(const std::string& keyword, const std::vector<double>& field);
 
     /// Sets a special field to have a particular value.
-    void setSpecialField(const std::string& keyword, std::tr1::shared_ptr<SpecialBase> field);
+    void setSpecialField(const std::string& keyword, SpecialFieldPtr field);
 
     /// Compute the units used by the deck, depending on the presence
     /// of keywords such as METRIC, FIELD etc.  It is an error to call
@@ -181,17 +192,25 @@ public:
     const EclipseUnits& units() const;
 
 private:
-    std::tr1::shared_ptr<SpecialBase> createSpecialField(std::istream& is, const std::string& fieldname);
+    SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname);
     void readImpl(std::istream& is);
 
 
     std::string directory_;
     std::map<std::string, std::vector<int> > integer_field_map_;
     std::map<std::string, std::vector<double> > floating_field_map_;
-    std::map<std::string, std::tr1::shared_ptr<SpecialBase> > special_field_map_;
+    // std::map<std::string, SpecialFieldPtr> special_field_map_;
     std::set<std::string> ignored_fields_;
-    std::tr1::shared_ptr<SpecialBase> empty_special_field_;
     EclipseUnits units_;
+
+    // For SCHEDULE handling.
+    enum ReadingMode { Regular, Timesteps };
+    ReadingMode current_reading_mode_;
+    boost::gregorian::date start_date_;
+    double current_time_days_;
+    int current_epoch_;
+    typedef std::map<std::string, SpecialFieldPtr> SpecialMap;
+    std::vector<SpecialMap> special_field_by_epoch_;
 };
 
 

opm/core/eclipse/SpecialEclipseFields.hpp

@@ -1658,6 +1658,67 @@ struct WELTARG : public SpecialBase
 };
 
 
+
+struct WelopenLine
+{
+    std::string well_;             // Well name or well name root
+    std::string openshutflag_;     // What to do with the well.
+};
+
+/// Class for keyword WELOPEN
+struct WELOPEN : public SpecialBase
+{
+    std::vector<WelopenLine> welopen;
+
+    WELOPEN()
+    {
+    }
+
+    virtual ~WELOPEN()
+    {}
+
+    virtual std::string name() const {return std::string("WELOPEN");}
+
+    virtual void read(std::istream& is)
+    {
+	while(is) {
+	    std::string wellname = readString(is); 
+	    if (wellname[0] == '/') {
+		is >> ignoreLine;
+		break;
+	    }
+	    while (wellname.find("--") == 0) {
+		// This line is a comment
+		is >> ignoreLine;
+		wellname = readString(is);
+	    }
+	    WelopenLine welopen_line;
+	    welopen_line.well_ = wellname;
+	    welopen_line.openshutflag_ = readString(is);
+	    ignoreSlashLine(is);
+	    welopen.push_back(welopen_line);
+	}
+    }
+
+    virtual void write(std::ostream& os) const
+    {
+	os << name() << std::endl;
+	for (int i=0; i<(int)welopen.size(); ++i) {
+	    os << welopen[i].well_ << "  " 
+	       << welopen[i].openshutflag_
+	       << std::endl;
+	}
+	os << std::endl;
+    }
+
+    virtual void convertToSI(const EclipseUnits& /*units*/)
+    {
+    }
+};
+
+
+
+
 /// Class holding a data line of keyword EQUIL
 struct EquilLine
 {

opm/core/linalg/LinearSolverAGMG.cpp

@@ -0,0 +1,124 @@
+/*
+  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
+
+#include <opm/core/linalg/LinearSolverAGMG.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <iterator>
+#include <stdexcept>
+#include <vector>
+
+#if HAVE_AGMG
+// Manual prototype of main gateway routine to DOUBLE PRECISION,
+// serial version of the AGMG software package.
+//
+// Note that both the matrix entries and column indices are writable.
+// The solver may permute the matrix entries within each row during
+// the setup phase.
+#define DAGMG_ FC_FUNC(dagmg, DAGMG)
+
+extern "C"
+void
+DAGMG_(const int*    N     ,    // System size
+       double*       sa    ,    // Non-zero elements
+       int*          ja    ,    // Column indices
+       const int*    ia    ,    // Row pointers
+       const double* f     ,    // Right-hand side
+       double*       x     ,    // Solution
+       int*          ijob  ,    // Main control parameter
+       int*          iprint,    // Message output unit
+       int*          nrest ,    // Number of GCR restarts
+       int*          iter  ,    // Maximum (and actual) number of iterations
+       const double* tol   );   // Residual reduction tolerance
+
+#endif // HAVE_AGMG
+
+namespace Opm
+{
+    LinearSolverAGMG::LinearSolverAGMG(const int    max_it,
+                                       const double rtol  ,
+                                       const bool   is_spd)
+        : max_it_(max_it),
+          rtol_  (rtol)  ,
+          is_spd_(is_spd)
+    {
+#if !HAVE_AGMG
+        THROW("AGMG support is not enabled in this library");
+#endif  // HAVE_AGMG
+    }
+
+    LinearSolverAGMG::~LinearSolverAGMG() {}
+    
+    LinearSolverInterface::LinearSolverReport
+    LinearSolverAGMG::solve(const int     size    ,
+                            const int     nonzeros,
+                            const int*    ia      ,
+                            const int*    ja      ,
+                            const double* sa      ,
+                            const double* rhs     ,
+                            double*       solution) const
+    {
+        const std::vector<double>::size_type nnz = ia[size];
+
+        assert (nnz == std::vector<double>::size_type(nonzeros));
+
+#if defined(NDEBUG)
+        // Suppress warning about unused parameter.
+        static_cast<void>(nonzeros);
+#endif
+
+        std::vector<double> a(sa, sa + nnz);
+
+        // Account for 1-based indexing.
+        std::vector<int>    i(ia, ia + std::vector<int>::size_type(size + 1));
+        std::transform(i.begin(), i.end(), i.begin(),
+                       std::bind2nd(std::plus<int>(), 1));
+
+        std::vector<int>    j(ja, ja + nnz);
+        std::transform(j.begin(), j.end(), j.begin(),
+                       std::bind2nd(std::plus<int>(), 1));
+
+        LinearSolverInterface::LinearSolverReport rpt = {};
+        rpt.iterations = max_it_;
+
+        int ijob = 0;           // Setup + solution + cleanup, x0==0.
+
+        int nrest;
+        if (is_spd_) {
+            nrest = 1;          // Use CG algorithm
+        }
+        else {
+            nrest = 10;         // Suggested default number of GCR restarts.
+        }
+
+        int iprint = 0;         // Suppress most output
+        DAGMG_(& size, & a[0], & j[0], & i[0], rhs, solution,
+               & ijob, & iprint, & nrest, & rpt.iterations, & rtol_);
+
+        rpt.converged = rpt.iterations <= max_it_;
+        return rpt;
+    }
+} // namespace Opm

opm/core/linalg/LinearSolverAGMG.hpp

@@ -0,0 +1,96 @@
+/*
+  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_LINEARSOLVERAGMG_HEADER_INCLUDED
+#define OPM_LINEARSOLVERAGMG_HEADER_INCLUDED
+
+/**
+ * \file
+ * Gateway to Notay's AGMG package implementing an aggregation-based
+ * algebraic multigrid method.
+ *
+ * References:
+ *   * Y. Notay,
+ *     An aggregation-based algebraic multigrid method,
+ *     Electronic Transactions on Numerical Analysis, vol 37,
+ *     pp. 123-146, 2010.
+ *
+ *   * A. Napov and Y. Notay,
+ *     An algebraic multigrid method with guaranteed convergence rate,
+ *     Report GANMN 10-03, Université Libre de Bruxelles, Brussels,
+ *     Belgium, 2010 (Revised 2011).
+ *
+ *   * Y. Notay,
+ *     Aggregation-based algebraic multigrid for convection-diffusion
+ *     equations, Report GANMN 11-01, Université Libre de Bruxelles,
+ *     Brussels, Belgium, 2011.
+ */
+
+#include <opm/core/linalg/LinearSolverInterface.hpp>
+
+namespace Opm
+{
+    /// Concrete class encapsulating Notay's AGMG package
+    class LinearSolverAGMG : public LinearSolverInterface
+    {
+    public:
+        /**
+         * Constructor.
+         * \param[in] max_it  Maximum number of solver iterations.
+         * \param[in] rtol    Residual reduction tolerance.
+         * \param[in] is_spd  Whether or not the matrix is SPD.  SPD
+         *                    systems are solved using CG while other
+         *                    systems are solved using GCR.
+         */
+        LinearSolverAGMG(const int    max_it = 100   ,
+                         const double rtol   = 1.0e-6,
+                         const bool   is_spd = false);
+
+        /// Destructor.
+        virtual ~LinearSolverAGMG();
+
+        using LinearSolverInterface::solve;
+
+        /// Solve a linear system, with a matrix given in compressed
+        /// sparse row format.
+        /// \param[in] size        Number of rows (and colums).
+        /// \param[in] nonzeros    Number of (structural) non-zeros.
+        /// \param[in] ia          Row pointers.
+        /// \param[in] ja          Column indices.
+        /// \param[in] sa          (structurally) non-zero elements.
+        /// \param[in] rhs         System right-hand side.
+        /// \param[inout] solution System solution.
+        /// \return Solver meta-data concerning most recent system solve.
+        virtual LinearSolverInterface::LinearSolverReport
+        solve(const int size, const int nonzeros,
+              const int* ia, const int* ja, const double* sa,
+              const double* rhs, double* solution) const;
+
+    private:
+        int    max_it_;
+        double rtol_  ;
+        bool   is_spd_;
+    };
+
+
+} // namespace Opm
+
+
+
+#endif // OPM_LINEARSOLVERAGMG_HEADER_INCLUDED

opm/core/pressure/CompressibleTpfa.cpp

@@ -28,8 +28,8 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/newwells.h>
-#include <opm/core/BlackoilState.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
-#include <opm/core/WellState.hpp>
+#include <opm/core/simulator/WellState.hpp>
 
 #include <algorithm>
 #include <cmath>

opm/core/simulator/SimulatorTwophase.cpp

@@ -0,0 +1,461 @@
+/*
+  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/SimulatorTwophase.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+
+#include <opm/core/pressure/IncompTpfa.hpp>
+
+#include <opm/core/grid.h>
+#include <opm/core/newwells.h>
+#include <opm/core/pressure/flow_bc.h>
+
+#include <opm/core/utility/SimulatorTimer.hpp>
+#include <opm/core/utility/StopWatch.hpp>
+#include <opm/core/utility/writeVtkData.hpp>
+#include <opm/core/utility/miscUtilities.hpp>
+
+#include <opm/core/fluid/IncompPropertiesInterface.hpp>
+#include <opm/core/fluid/RockCompressibility.hpp>
+
+#include <opm/core/utility/ColumnExtract.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
+#include <opm/core/simulator/WellState.hpp>
+#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
+
+#include <boost/filesystem/convenience.hpp>
+#include <boost/scoped_ptr.hpp>
+#include <boost/lexical_cast.hpp>
+
+#include <numeric>
+#include <fstream>
+
+
+namespace Opm
+{
+
+    class SimulatorTwophase::Impl
+    {
+    public:
+        Impl(const parameter::ParameterGroup& param,
+             const UnstructuredGrid& grid,
+             const IncompPropertiesInterface& props,
+             const RockCompressibility* rock_comp,
+             const Wells* wells,
+             const std::vector<double>& src,
+             const FlowBoundaryConditions* bcs,
+             const LinearSolverInterface& linsolver,
+             const double* gravity);
+
+        void run(SimulatorTimer& timer,
+                 TwophaseState& state,
+                 WellState& well_state);
+
+    private:
+        // Data.
+
+        // Parameters for output.
+        bool output_;
+        std::string output_dir_;
+        int output_interval_;
+        // Parameters for pressure solver.
+        int nl_pressure_maxiter_;
+        double nl_pressure_tolerance_;
+        // Parameters for transport solver.
+        int nl_maxiter_;
+        double nl_tolerance_;
+        int num_transport_substeps_;
+        bool use_segregation_split_;
+        // Observed objects.
+        const UnstructuredGrid& grid_;
+        const IncompPropertiesInterface& props_;
+        const RockCompressibility* rock_comp_;
+        const Wells* wells_;
+        const std::vector<double>& src_;
+        const FlowBoundaryConditions* bcs_;
+        const LinearSolverInterface& linsolver_;
+        const double* gravity_;
+        // Solvers
+        IncompTpfa psolver_;
+        TransportModelTwophase tsolver_;
+        // Needed by column-based gravity segregation solver.
+        std::vector< std::vector<int> > columns_;
+        // Misc. data
+        std::vector<int> allcells_;
+    };
+
+
+
+
+    SimulatorTwophase::SimulatorTwophase(const parameter::ParameterGroup& param,
+                                         const UnstructuredGrid& grid,
+                                         const IncompPropertiesInterface& props,
+                                         const RockCompressibility* rock_comp,
+                                         const Wells* wells,
+                                         const std::vector<double>& src,
+                                         const FlowBoundaryConditions* bcs,
+                                         const LinearSolverInterface& linsolver,
+                                         const double* gravity)
+    {
+        pimpl_.reset(new Impl(param, grid, props, rock_comp, wells, src, bcs, linsolver, gravity));
+    }
+
+
+
+
+
+    void SimulatorTwophase::run(SimulatorTimer& timer,
+                                TwophaseState& state,
+                                WellState& well_state)
+    {
+        pimpl_->run(timer, state, well_state);
+    }
+
+
+
+    static void outputState(const UnstructuredGrid& grid,
+                            const Opm::TwophaseState& state,
+                            const int step,
+                            const std::string& output_dir)
+    {
+        // Write data in VTK format.
+        std::ostringstream vtkfilename;
+        vtkfilename << output_dir << "/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);
+
+        // 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 << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
+            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);
+    }
+
+
+
+
+
+    SimulatorTwophase::Impl::Impl(const parameter::ParameterGroup& param,
+                                  const UnstructuredGrid& grid,
+                                  const IncompPropertiesInterface& props,
+                                  const RockCompressibility* rock_comp,
+                                  const Wells* wells,
+                                  const std::vector<double>& src,
+                                  const FlowBoundaryConditions* bcs,
+                                  const LinearSolverInterface& linsolver,
+                                  const double* gravity)
+        : grid_(grid),
+          props_(props),
+          rock_comp_(rock_comp),
+          wells_(wells),
+          src_(src),
+          bcs_(bcs),
+          linsolver_(linsolver),
+          gravity_(gravity),
+          psolver_(grid, props.permeability(), gravity, linsolver),
+          tsolver_(grid, props, 1e-9, 30)
+    {
+        // For output.
+        output_ = param.getDefault("output", true);
+        if (output_) {
+            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);
+        }
+
+        // For pressure solver
+        nl_pressure_maxiter_ = param.getDefault("nl_pressure_maxiter", 10);
+        nl_pressure_tolerance_ = param.getDefault("nl_pressure_tolerance", 1.0); // Pascal
+
+        // For transport solver.
+        nl_maxiter_ = param.getDefault("nl_maxiter", 30);
+        nl_tolerance_ = param.getDefault("nl_tolerance", 1e-9);
+        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;
+        }
+    }
+
+
+
+
+    void SimulatorTwophase::Impl::run(SimulatorTimer& timer,
+                                      TwophaseState& state,
+                                      WellState& well_state)
+    {
+        std::vector<double> totmob;
+        std::vector<double> omega; // Will remain empty if no gravity.
+        std::vector<double> rc; // Will remain empty if no rock compressibility.
+        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);
+
+
+        // Main simulation loop.
+        Opm::time::StopWatch pressure_timer;
+        double ptime = 0.0;
+        Opm::time::StopWatch transport_timer;
+        double ttime = 0.0;
+        Opm::time::StopWatch total_timer;
+        total_timer.start();
+        std::cout << "\n\n================    Starting main simulation loop     ===============" << std::endl;
+        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;
+        int num_wells = 0;
+        if (wells_) {
+            num_wells = wells_->number_of_wells;
+            well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
+            wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
+        }
+        const int num_cells = grid_.number_of_cells;
+        for (; !timer.done(); ++timer) {
+            // Report timestep and (optionally) write state to disk.
+            timer.report(std::cout);
+            if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
+                outputState(grid_, state, timer.currentStepNum(), output_dir_);
+            }
+
+            // Solve pressure.
+            if (gravity_) {
+                computeTotalMobilityOmega(props_, allcells_, state.saturation(), totmob, omega);
+            } else {
+                computeTotalMobility(props_, allcells_, state.saturation(), totmob);
+            }
+            std::vector<double> wdp;
+            if (wells_) {
+                Opm::computeWDP(*wells_, grid_, state.saturation(), props_.density(),
+                                gravity_ ? gravity_[2] : 0.0, true, wdp);
+            }
+            do {
+                pressure_timer.start();
+                if (rock_comp_ && rock_comp_->isActive()) {
+                    rc.resize(num_cells);
+                    std::vector<double> initial_pressure = state.pressure();
+                    std::vector<double> initial_porevolume(num_cells);
+                    computePorevolume(grid_, props_.porosity(), *rock_comp_, initial_pressure, initial_porevolume);
+                    std::vector<double> pressure_increment(num_cells + num_wells);
+                    std::vector<double> prev_pressure(num_cells + num_wells);
+                    for (int iter = 0; iter < nl_pressure_maxiter_; ++iter) {
+                        for (int cell = 0; cell < num_cells; ++cell) {
+                            rc[cell] = rock_comp_->rockComp(state.pressure()[cell]);
+                        }
+                        computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol);
+                        std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
+                        std::copy(well_state.bhp().begin(), well_state.bhp().end(), prev_pressure.begin() + num_cells);
+                        // prev_pressure = state.pressure();
+
+                        // compute pressure increment
+                        psolver_.solveIncrement(totmob, omega, src_, wdp, bcs_, porevol, rc,
+                                                prev_pressure, initial_porevolume, timer.currentStepLength(),
+                                                pressure_increment);
+
+                        double max_change = 0.0;
+                        for (int cell = 0; cell < num_cells; ++cell) {
+                            state.pressure()[cell] += pressure_increment[cell];
+                            max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
+                        }
+                        for (int well = 0; well < num_wells; ++well) {
+                            well_state.bhp()[well] += pressure_increment[num_cells + well];
+                            max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
+                        }
+
+                        std::cout << "Pressure iter " << iter << "   max change = " << max_change << std::endl;
+                        if (max_change < nl_pressure_tolerance_) {
+                            break;
+                        }
+                    }
+                    psolver_.computeFaceFlux(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(),
+                                             well_state.bhp(), well_state.perfRates());
+                } else {
+                    psolver_.solve(totmob, omega, src_, wdp, bcs_, state.pressure(), state.faceflux(),
+                                   well_state.bhp(), well_state.perfRates());
+                }
+                pressure_timer.stop();
+                double pt = pressure_timer.secsSinceStart();
+                std::cout << "Pressure solver took:  " << pt << " seconds." << std::endl;
+                ptime += pt;
+            } while (false);
+
+            // 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], &porevol[0], &transport_src[0],
+                              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());
+                }
+            }
+            transport_timer.stop();
+            double tt = transport_timer.secsSinceStart();
+            std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
+            ttime += tt;
+
+            // Report volume balances.
+            Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
+            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.saturation(),
+                                timer.currentTime() + timer.currentStepLength(),
+                                well_state.bhp(), well_state.perfRates());
+            }
+        }
+        total_timer.stop();
+
+        std::cout << "\n\n================    End of simulation     ===============\n"
+                  << "Total time taken: " << total_timer.secsSinceStart()
+                  << "\n  Pressure time:  " << ptime
+                  << "\n  Transport time: " << ttime << std::endl;
+
+        if (output_) {
+            outputState(grid_, state, timer.currentStepNum(), output_dir_);
+            outputWaterCut(watercut, output_dir_);
+            if (wells_) {
+                outputWellReport(wellreport, output_dir_);
+            }
+        }
+
+
+
+    }
+
+} // namespace Opm

opm/core/simulator/SimulatorTwophase.hpp

@@ -0,0 +1,95 @@
+/*
+  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_SIMULATORTWOPHASE_HEADER_INCLUDED
+#define OPM_SIMULATORTWOPHASE_HEADER_INCLUDED
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+
+struct UnstructuredGrid;
+struct Wells;
+struct FlowBoundaryConditions;
+
+namespace Opm
+{
+    namespace parameter { class ParameterGroup; }
+    class IncompPropertiesInterface;
+    class RockCompressibility;
+    class LinearSolverInterface;
+    class SimulatorTimer;
+    class TwophaseState;
+    class WellState;
+
+    /// Class collecting all necessary components for a two-phase simulation.
+    class SimulatorTwophase
+    {
+    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_maxiter (10)       max nonlinear iterations in pressure
+        ///     nl_pressure_tolerance (1.0)    pressure solver nonlinear tolerance (in Pascal)
+        ///     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] wells       if non-null, wells data structure
+        /// \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
+       SimulatorTwophase(const parameter::ParameterGroup& param,
+                         const UnstructuredGrid& grid,
+                         const IncompPropertiesInterface& props,
+                         const RockCompressibility* rock_comp,
+                         const Wells* wells,
+                         const std::vector<double>& src,
+                         const FlowBoundaryConditions* bcs,
+                         const 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
+        void run(SimulatorTimer& timer,
+                 TwophaseState& 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_SIMULATORTWOPHASE_HEADER_INCLUDED

opm/core/simulator/WellState.hpp

@@ -26,9 +26,11 @@
 namespace Opm
 {
 
+    /// The state of a set of wells.
     class WellState
     {
     public:
+        /// Allocate and initialize if wells is non-null.
         template <class State>
         void init(const Wells* wells, const State& state)
         {
@@ -44,9 +46,11 @@ namespace Opm
             }
         }
 
+        /// One bhp pressure per well.
         std::vector<double>& bhp() { return bhp_; }
         const std::vector<double>& bhp() const { return bhp_; }
 
+        /// One rate per well connection.
         std::vector<double>& perfRates() { return perfrates_; }
         const std::vector<double>& perfRates() const { return perfrates_; }
 

opm/core/wells/InjectionSpecification.cpp

@@ -1,4 +1,5 @@
-#include <opm/core/InjectionSpecification.hpp>
+#include <opm/core/wells/InjectionSpecification.hpp>
+
 namespace Opm
 {
 

opm/core/wells/ProductionSpecification.cpp

@@ -1,4 +1,4 @@
-#include <opm/core/ProductionSpecification.hpp>
+#include <opm/core/wells/ProductionSpecification.hpp>
 
 namespace Opm
 {

opm/core/wells/WellCollection.cpp

@@ -18,7 +18,7 @@ You should have received a copy of the GNU General Public License
 along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include "WellCollection.hpp"
+#include <opm/core/wells/WellCollection.hpp>
 
 namespace Opm
 {

opm/core/wells/WellCollection.hpp

@@ -23,7 +23,7 @@ along with OpenRS.  If not, see <http://www.gnu.org/licenses/>.
 #define	OPM_WELLCOLLECTION_HPP
 
 #include <vector>
-#include <opm/core/WellsGroup.hpp>
+#include <opm/core/wells/WellsGroup.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 

opm/core/wells/WellsGroup.cpp

@@ -17,7 +17,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <opm/core/WellsGroup.hpp>
+#include <opm/core/wells/WellsGroup.hpp>
 #include <cmath>
 #include <opm/core/newwells.h>
 #include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>

opm/core/wells/WellsGroup.hpp

@@ -20,8 +20,8 @@
 #ifndef OPM_WELLSGROUP_HPP
 #define	OPM_WELLSGROUP_HPP
 
-#include <opm/core/InjectionSpecification.hpp>
+#include <opm/core/wells/InjectionSpecification.hpp>
-#include <opm/core/ProductionSpecification.hpp>
+#include <opm/core/wells/ProductionSpecification.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>

opm/core/wells/WellsManager.cpp

@@ -18,13 +18,13 @@
 */
 
 
-#include <opm/core/WellsManager.hpp>
+#include <opm/core/wells/WellsManager.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/newwells.h>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/Units.hpp>
-#include <opm/core/WellCollection.hpp>
+#include <opm/core/wells/WellCollection.hpp>
 #include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
 
 #include <tr1/array>
@@ -435,13 +435,13 @@ namespace Opm
                         // Add all controls that are present in well.
                         int ok = 1;
                         int control_pos[5] = { -1, -1, -1, -1, -1 };
-                        if (ok && wci_line.surface_flow_max_rate_ > 0.0) {
+                        if (ok && wci_line.surface_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num;
                             const double distr[3] = { 1.0, 1.0, 1.0 };
                             ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_,
                                                       distr, wix, w_);
                         }
-                        if (ok && wci_line.reservoir_flow_max_rate_ > 0.0) {
+                        if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num;
                             const double distr[3] = { 1.0, 1.0, 1.0 };
                             ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_,
@@ -515,7 +515,7 @@ namespace Opm
                         // Add all controls that are present in well.
                         int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 };
                         int ok = 1;
-                        if (ok && wcp_line.oil_max_rate_ > 0.0) {
+                        if (ok && wcp_line.oil_max_rate_ >= 0.0) {
                             if (!pu.phase_used[BlackoilPhases::Liquid]) {
                                 THROW("Oil phase not active and ORAT control specified.");
                             }
@@ -525,7 +525,7 @@ namespace Opm
                             ok = append_well_controls(SURFACE_RATE, -wcp_line.oil_max_rate_,
                                                       distr, wix, w_);
                         }
-                        if (ok && wcp_line.water_max_rate_ > 0.0) {
+                        if (ok && wcp_line.water_max_rate_ >= 0.0) {
                             if (!pu.phase_used[BlackoilPhases::Aqua]) {
                                 THROW("Water phase not active and WRAT control specified.");
                             }
@@ -535,7 +535,7 @@ namespace Opm
                             ok = append_well_controls(SURFACE_RATE, -wcp_line.water_max_rate_,
                                                       distr, wix, w_);
                         }
-                        if (ok && wcp_line.gas_max_rate_ > 0.0) {
+                        if (ok && wcp_line.gas_max_rate_ >= 0.0) {
                             if (!pu.phase_used[BlackoilPhases::Vapour]) {
                                 THROW("Gas phase not active and GRAT control specified.");
                             }
@@ -545,7 +545,7 @@ namespace Opm
                             ok = append_well_controls(SURFACE_RATE, -wcp_line.gas_max_rate_,
                                                       distr, wix, w_);
                         }
-                        if (ok && wcp_line.liquid_max_rate_ > 0.0) {
+                        if (ok && wcp_line.liquid_max_rate_ >= 0.0) {
                             if (!pu.phase_used[BlackoilPhases::Aqua]) {
                                 THROW("Water phase not active and LRAT control specified.");
                             }
@@ -559,7 +559,7 @@ namespace Opm
                             ok = append_well_controls(SURFACE_RATE, -wcp_line.liquid_max_rate_,
                                                       distr, wix, w_);
                         }
-                        if (ok && wcp_line.reservoir_flow_max_rate_ > 0.0) {
+                        if (ok && wcp_line.reservoir_flow_max_rate_ >= 0.0) {
                             control_pos[ProductionControl::RESV] = w_->ctrls[wix]->num;
                             double distr[3] = { 1.0, 1.0, 1.0 };
                             ok = append_well_controls(RESERVOIR_RATE, -wcp_line.reservoir_flow_max_rate_,

opm/core/wells/WellsManager.hpp

@@ -21,8 +21,8 @@
 #define OPM_WELLSMANAGER_HEADER_INCLUDED
 
 
-#include <opm/core/WellCollection.hpp>
+#include <opm/core/wells/WellCollection.hpp>
-#include <opm/core/WellsGroup.hpp>
+#include <opm/core/wells/WellsGroup.hpp>
 
 struct Wells;
 struct UnstructuredGrid;

tests/Makefile.am

@@ -47,3 +47,9 @@ noinst_PROGRAMS    += test_cfs_tpfa
 noinst_PROGRAMS    += test_jacsys
 test_jacsys_SOURCES = test_jacsys.cpp
 endif
+
+
+if BUILD_AGMG
+noinst_PROGRAMS   += test_agmg
+test_agmg_SOURCES  = test_agmg.cpp
+endif

tests/test_agmg.cpp

@@ -0,0 +1,126 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+  Copyright 2012 Statoil ASA.
+
+  This file is part of the Open Porous Media Project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <iomanip>
+#include <iostream>
+#include <vector>
+
+#include <boost/shared_ptr.hpp>
+
+#include <opm/core/linalg/sparse_sys.h>
+#include <opm/core/linalg/LinearSolverAGMG.hpp>
+
+namespace {
+    std::size_t
+    compute_nnz(const std::size_t m)
+    {
+        assert (m > 0);
+
+        std::size_t nnz = m;              // A(i,i)
+
+        nnz += (m > 1) ? 2           : 0; // A(0,1), A(m-1,m-2)
+        nnz += (m > 2) ? 2 * (m - 2) : 0; // A(i,i-1), A(i,i+1)
+
+        return nnz;
+    }
+
+
+    boost::shared_ptr<CSRMatrix>
+    build_laplace_1d(const std::size_t m)
+    {
+        assert (m >= 2);
+
+        const std::size_t nnz = compute_nnz(m);
+
+        boost::shared_ptr<CSRMatrix>
+            A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete);
+
+        A->ia[ 0 ] = 0;
+
+        // First row
+        A->ia[       0 + 1     ] = A->ia[ 0 ];
+        A->ja[ A->ia[0 + 1]    ] = 0 + 0;
+        A->sa[ A->ia[0 + 1] ++ ] =   2.0;
+        A->ja[ A->ia[0 + 1]    ] = 0 + 1;
+        A->sa[ A->ia[0 + 1] ++ ] = - 1.0;
+
+        // General rows
+        for (std::size_t i = 1; i < m - 1; ++i) {
+            A->ia[i + 1] = A->ia[i];
+
+            A->ja[ A->ia[i + 1]    ] = int(i) - 1;
+            A->sa[ A->ia[i + 1] ++ ] = - 1.0;
+
+            A->ja[ A->ia[i + 1]    ] = int(i)    ;
+            A->sa[ A->ia[i + 1] ++ ] =   2.0;
+
+            A->ja[ A->ia[i + 1]    ] = int(i) + 1;
+            A->sa[ A->ia[i + 1] ++ ] = - 1.0;
+        }
+
+        // Last row
+        A->ia[ (m - 1) + 1 ] = A->ia[ m - 1 ];
+
+        A->ja[ A->ia[ (m - 1) + 1 ]    ] = int(m - 1) - 1;
+        A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = - 1.0;
+
+        A->ja[ A->ia[ (m - 1) + 1 ]    ] = int(m - 1)    ;
+        A->sa[ A->ia[ (m - 1) + 1 ] ++ ] =   2.0;
+
+        return A;
+    }
+}
+
+
+int main()
+{
+    const std::size_t m = 10;
+
+    boost::shared_ptr<CSRMatrix> A = build_laplace_1d(m);
+
+    // Form right-hand side [1, 0, 0, ...., 0, 1]
+    std::vector<double> b(m, 0.0);
+    b[0] = 1.0; b.back() = 1.0;
+
+    // Allocate solution vector
+    std::vector<double> x(m);
+
+    // Create solver for SPD system.
+    Opm::LinearSolverAGMG linsolve(100, 1e-9, true);
+
+    Opm::LinearSolverInterface::LinearSolverReport
+        rpt = linsolve.solve(A.get(), & b[0], & x[0]);
+
+    double e = 0.0;
+    for (std::size_t i = 0; i < m; ++i) {
+        const double d = x[i] - 1.0;
+        e += d * d;
+    }
+
+    std::cerr << "|| e ||_2 = "
+              << std::scientific
+              << std::setprecision(5)
+              << std::sqrt(e) / double(m) << '\n';
+
+    return 0;
+}

tests/test_column_extract.cpp

@@ -6,7 +6,7 @@
 
 #define BOOST_TEST_MODULE ColumnExtractTest
 #include <boost/test/unit_test.hpp>
-#include <opm/core/ColumnExtract.hpp>
+#include <opm/core/utility/ColumnExtract.hpp>
 #include <opm/core/GridManager.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 

tutorials/tutorial3.cpp

@@ -37,7 +37,7 @@
 
 #include <opm/core/transport/reorder/TransportModelTwophase.hpp>
 
-#include <opm/core/TwophaseState.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
 
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/Units.hpp>

tutorials/tutorial4.cpp

@@ -44,12 +44,12 @@
 #include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
 #include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
 
-#include <opm/core/TwophaseState.hpp>
+#include <opm/core/simulator/TwophaseState.hpp>
 
 #include <opm/core/utility/miscUtilities.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
-#include <opm/core/WellCollection.hpp>
+#include <opm/core/wells/WellCollection.hpp>
 
 /// \page tutorial4 Well controls
 ///